[Bio] / FigWebServices / diffsF.cgi Repository:
ViewVC logotype

Annotation of /FigWebServices/diffsF.cgi

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (view) (download)

1 : overbeek 1.2 # -*- perl -*-
2 :     #
3 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     #
18 :    
19 : overbeek 1.1 use GenoGraphics;
20 :    
21 :     use FIG;
22 :     my $fig = new FIG;
23 :    
24 : overbeek 1.4 use FigFam;
25 :     my $default_dir = "$FIG_Config::data/FigfamsData";
26 :    
27 : overbeek 1.1 use HTML;
28 :    
29 :     use CGI;
30 :     my $cgi = new CGI;
31 :    
32 :     if (0)
33 :     {
34 :     my $VAR1;
35 :     eval(join("",`cat /tmp/diffsF_parms`));
36 :     $cgi = $VAR1;
37 :     # print STDERR &Dumper($cgi);
38 :     }
39 :    
40 :     if (0)
41 :     {
42 :     print $cgi->header;
43 :     my @params = $cgi->param;
44 :     print "<pre>\n";
45 :     foreach $_ (@params)
46 :     {
47 :     print "$_\t:",join(",",$cgi->param($_)),":\n";
48 :     }
49 :    
50 :     if (0)
51 :     {
52 :     if (open(TMP,">/tmp/diffsF_parms"))
53 :     {
54 :     print TMP &Dumper($cgi);
55 :     close(TMP);
56 :     }
57 :     }
58 :     exit;
59 :     }
60 :    
61 :     use constant FID => 0;
62 :     use constant LOCUS => 1;
63 :     use constant CONTIG => 2;
64 :     use constant START => 3;
65 :     use constant STOP => 4;
66 :     use constant LEN => 5;
67 :     use constant STRAND => 6;
68 :     use constant TYPE => 7;
69 :     use constant TAXID => 8;
70 :    
71 :     my $html = [];
72 :     unshift @$html, "<TITLE>Feature Differences</TITLE>\n";
73 : overbeek 1.4 $org1 = $cgi->param('org1');
74 :     $org2 = $cgi->param('org2');
75 :     $gene = $cgi->param('gene');
76 :     $align1 = $cgi->param('align1');
77 :     $align2 = $cgi->param('align2');
78 : overbeek 1.1
79 :     if (!((-d $org1) || (-d $org2)))
80 :     {
81 : overbeek 1.4 push(@$html,$cgi->h1("bad parms: gene=$gene align1=$align1 align2=$align2 org1=$org1 org2=$org2"));
82 : overbeek 1.1 }
83 :     elsif ($gene = $cgi->param('gene'))
84 :     {
85 :     &show_region($fig,$cgi,$html,$org1,$gene);
86 :     }
87 : overbeek 1.4 elsif (($gene1 = $cgi->param('align1')) && ($gene2 = $cgi->param('align2')))
88 :     {
89 :     &show_alignment($fig, $cgi, $html, $gene1, $gene2, $org1, $org2);
90 :     }
91 : overbeek 1.1 else
92 :     {
93 : overbeek 1.4 &list_diffs($fig, $cgi, $html, $org1, $org2);
94 : overbeek 1.1 }
95 :     &HTML::show_page($cgi,$html);
96 :    
97 :     sub show_region {
98 :     my($fig,$cgi,$html,$org,$gene) = @_;
99 :     my @feat = &load_tbl2($org,$gene);
100 :     for ($i=0; ($i < @feat) && ($feat[$i]->[FID] ne $gene); $i++) {}
101 :     if ($i < @feat)
102 :     {
103 :     $first = ($i > 4) ? $i-5 : 0;
104 :     $last = ($i < (@feat - 5)) ? $i+5 : @feat - 1;
105 :     my $gg = [];
106 :     my $genes = [];
107 :     $min = 10000000;
108 :     $max = 0;
109 :     while ($first <= $last)
110 :     {
111 :     $f = $feat[$first];
112 :     $beg = $f->[START];
113 :     $end = $f->[STOP];
114 :     $min = &FIG::min($min,$beg);
115 :     $min = &FIG::min($min,$end);
116 :     $max = &FIG::max($max,$beg);
117 :     $max = &FIG::max($max,$end);
118 :     $fid = $f->[FID];
119 : overbeek 1.4 # $fid =~ /\.([a-z]+\.\d+)$/;
120 :     # $info = $1;
121 : overbeek 1.1 push(@$genes,[ &FIG::min($beg,$end),
122 :     &FIG::max($beg,$end),
123 :     ($beg < $end) ? "rightArrow" : "leftArrow",
124 :     ($fid eq $gene) ? "red" : "blue",
125 :     "",
126 :     "",
127 :     "$fid: $beg to $end",
128 :     ""
129 :     ] );
130 :     $first++;
131 :     }
132 :     my $map = ['Region',0,$max+1-$min,&decr_coords($genes,$min)];
133 :     push(@$gg,$map);
134 :     push( @$html, @{ &GenoGraphics::render( $gg, 1000, 4, 1 ) } );
135 :     }
136 :     }
137 :    
138 :     sub decr_coords {
139 :     my($genes,$min) = @_;
140 :     my($gene);
141 :    
142 :     foreach $gene (@$genes)
143 :     {
144 :     $gene->[0] -= $min;
145 :     $gene->[1] -= $min;
146 :     }
147 :     return $genes;
148 :     }
149 :    
150 :     sub list_diffs {
151 : overbeek 1.3 my ($fig, $cgi, $html, $org1, $org2) = @_;
152 :    
153 : overbeek 1.1 $org1T = &load_tbl1($org1);
154 :     $org2T = &load_tbl1($org2);
155 : overbeek 1.3
156 :     my @both;
157 :     my @just1;
158 :     foreach $key (sort { by_locus($org1T->{$a}, $org1T->{$b}) } keys %$org1T)
159 : overbeek 1.1 {
160 :     $type1 = $org1T->{$key}->[TYPE];
161 :     $locus1 = $org1T->{$key}->[LOCUS];
162 :    
163 :     if ($org2T->{$key})
164 :     {
165 :     $type2 = $org2T->{$key}->[TYPE];
166 :     $locus2 = $org2T->{$key}->[LOCUS];
167 :     if (($type1 ne $type2) || ($locus1 ne $locus2))
168 :     {
169 :     $fid1 = $org1T->{$key}->[FID];
170 :     $fid2 = $org2T->{$key}->[FID];
171 :     $link1 = &link($org1,$fid1);
172 :     $link2 = &link($org2,$fid2);
173 : overbeek 1.2
174 :     if (($len1 = $org1T->{$key}->[LEN]) < ($len2 = $org2T->{$key}->[LEN])) {
175 :     $changed = qq(lengthened);
176 :     }
177 :     else {
178 :     $changed = qq(shortened);
179 :     }
180 : overbeek 1.3
181 : overbeek 1.4 push @both, [ $changed, $link1, $len1, $link2, $len2, ($len2 - $len1) ];
182 : overbeek 1.1 }
183 :     }
184 :     else
185 :     {
186 :     $fid1 = $org1T->{$key}->[FID];
187 :     $link1 = &link($org1,$fid1);
188 : overbeek 1.4 push @just1, [ $link1, $org1T->{$key}->[LEN] ];
189 : overbeek 1.1 }
190 :     }
191 : overbeek 1.3
192 :     my @just2;
193 : overbeek 1.1 foreach $key (sort { by_locus($org2T->{$a},$org2T->{$b}) } keys %$org2T)
194 :     {
195 :     if (! $org1T->{$key})
196 :     {
197 :     $type2 = $org2T->{$key}->[TYPE];
198 :     $locus2 = $org2T->{$key}->[LOCUS];
199 :     $fid2 = $org2T->{$key}->[FID];
200 :     $link2 = &link($org2,$fid2);
201 : overbeek 1.3 push @just2, [ $link2, $org2T->{$key}->[LEN] ];
202 : overbeek 1.1 }
203 :     }
204 : overbeek 1.3
205 : overbeek 1.4 push @$html, map {
206 :     $align = qq(<A HREF="diffsF.cgi?align1=$_->[1]&align2=$_->[3]&org1=$org1&org2=$org2">align</A>);
207 :     $cgi->h2("$align $_->[0]: $_->[1] ($_->[2] bp) => $_->[3] ($_->[4] bp), diff=$_->[5]")
208 :     }
209 :     sort {
210 :     abs($b->[5]) <=> abs($a->[5])
211 :     } @both;
212 :    
213 :    
214 :     push @$html, map {
215 :     $cgi->h2("just2: $_->[0] ($_->[1] bp)")
216 :     }
217 :     sort {
218 :     $b->[1] <=> $a->[1]
219 :     } @just2;
220 :    
221 :    
222 :     push @$html, map {
223 :     $cgi->h2("just1: $_->[0] ($_->[1] bp)")
224 :     }
225 :     sort {
226 :     $b->[1] <=> $a->[1]
227 :     } @just1;
228 :     }
229 :    
230 :     sub show_alignment {
231 :     my ($fig, $cgi, $html, $gene1, $gene2, $org1, $org2) = @_;
232 :    
233 :     my %neighbors = map { m/^(\d+\.\d+)/o ? ($1 => 1) : () } `cat $org2/neighbors`;
234 :    
235 :     my (@tmp, $tmp);
236 :     if (@tmp = grep { m/^(\S+)/o && ($1 eq $gene2) } `cat $org2/found`) {
237 :     if ($tmp[0] =~ m/^\S+\t(FIG\d{6})/o) {
238 :     $figfam = new FigFam($fig, $1, $default_dir);
239 :     my @members = grep { $neighbors{$fig->genome_of($_)} } $figfam->list_members();
240 :    
241 :     my $tmp_seqs = "$FIG_Config::temp/tmp.$$";
242 :     open(TMP, ">$tmp_seqs")
243 :     || ((push @$html, $cgi->h2("Could not write-open $tmp_seqs")) && return);
244 :    
245 :     foreach my $id (@members) {
246 :     my $seq = $fig->get_translation($id);
247 :     &FIG::display_id_and_seq($id, \$seq, \*TMP);
248 :     }
249 :     close(TMP)
250 :     || ((push @$html, $cgi->h2("Could not close $tmp_seqs")) && return);
251 :    
252 :     $tmp = $gene1;
253 :     $tmp =~ s{\|}{\\\|}o;
254 :     system("echo $tmp1 | pull_fasta_entries $org1/Features/peg/fasta >> $tmp_seqs")
255 :     && ((push @$html, $cgi->h2("Could not append seq of $gene1 to $tmp_seqs")) && return);
256 :    
257 :     $tmp = $gene2;
258 :     $tmp =~ s{\|}{\\\|}o;
259 :     system("echo $tmp1 | pull_fasta_entries $org2/Features/peg/fasta >> $tmp_seqs")
260 :     && ((push @$html, $cgi->h2("Could not append seq of $gene2 to $tmp_seqs")) && return);
261 :    
262 :     system("clustalw -infile=$tmp_seqs -align -outorder=aligned > /dev/null")
263 :     && ((push @$html, $cgi->h2("clustalw failed on file $tmp_seqs")) && return);
264 :    
265 :     push @$html, qq(<pre>\n);
266 :     push @$html, `cat $tmp_seqs.aln`;
267 :     push @$html, qq(</pre>\n);
268 :     }
269 :     else {
270 :     $tmp[0] =~ s/\t/ /go;
271 :     push @$html, $cgi->h2("Could not parse family from entry: $tmp[0]");
272 :     }
273 :     }
274 :     else {
275 :     push @$html, $cgi->h2("No family found for gene2=$gene2");
276 :     }
277 : overbeek 1.1 }
278 :    
279 :     sub load_tbl1
280 :     {
281 :     my ($dir) = @_;
282 :     my ($entry, $id, $locus, $contig, $beg, $end, $len, $strand, $taxid, $type);
283 :    
284 :     open(TBL, "cat $dir/Features/*/tbl |") || die "Could not open $dir";
285 :    
286 :     my $tbl = {};
287 :     while (defined($entry = <TBL>))
288 :     {
289 :     chomp $entry;
290 :    
291 :     ($id, $locus) = split /\t/, $entry;
292 :     $id =~ m/^[^\|]+\|(\d+\.\d+)\.([^\.]+)/;
293 :     ($taxid, $type) = ($1, $2);
294 :    
295 :     if ((($contig, $beg, $end, $len, $strand) = &from_locus($locus))
296 :     && defined($contig) && $contig
297 :     && defined($beg) && $beg
298 :     && defined($end) && $end
299 :     && defined($len) && $len
300 :     && defined($strand) && $strand
301 :     )
302 :     {
303 :     $tbl->{"$contig\t$strand$end"} = [ $id, $locus, $contig, $beg, $end, $len, $strand, $type, $taxid ];
304 :     }
305 :     else
306 :     {
307 :     warn "INVALID ENTRY:\t$entry\n";
308 :     }
309 :     }
310 :    
311 :     return $tbl;
312 :     }
313 :    
314 :     sub load_tbl2
315 :     {
316 :     my ($dir,$gene) = @_;
317 :    
318 : overbeek 1.3 my ($entry, $id, $locus, $contig, $beg, $end, $len, $strand, $taxid, $type, $keep);
319 : overbeek 1.1
320 :     open(TBL, "cat $dir/Features/*/tbl |") || die "Could not open $dir";
321 :    
322 :     my $tbl = {};
323 :     while (defined($entry = <TBL>))
324 :     {
325 :     chomp $entry;
326 :    
327 :     ($id, $locus) = split /\t/, $entry;
328 :     $id =~ m/^[^\|]+\|(\d+\.\d+)\.([^\.]+)/;
329 :     ($taxid, $type) = ($1, $2);
330 :    
331 :     if ((($contig, $beg, $end, $len, $strand) = &from_locus($locus))
332 :     && defined($contig) && $contig
333 :     && defined($beg) && $beg
334 :     && defined($end) && $end
335 :     && defined($len) && $len
336 :     && defined($strand) && $strand
337 :     )
338 :     {
339 : overbeek 1.3 push(@{$tbl->{$contig}}, [ $id, $locus, $contig, $beg, $end, $len, $strand, $type, $taxid ]);
340 : overbeek 1.1 if ($id eq $gene) { $keep = $contig }
341 :     }
342 :     else
343 :     {
344 :     warn "INVALID ENTRY:\t$entry\n";
345 :     }
346 :     }
347 : overbeek 1.5 $keep || die "no contig for gene=$gene";
348 : overbeek 1.1 return sort { by_locus($a,$b) } @{$tbl->{$keep}};
349 :     }
350 :    
351 :     sub from_locus
352 :     {
353 :     my ($locus) = @_;
354 :     my ($contig, $beg, $end);
355 :    
356 :     if ( ($locus =~ m/^([^,]+)_(\d+)_\d+/) && ($contig = $1) && ($beg = $2)
357 :     && ($locus =~ m/[^,]+_\d+_(\d+)$/) && ($end = $1)
358 :     )
359 :     {
360 :     return ($contig, $beg, $end, (1+abs($end-$beg)), (($beg < $end) ? '+' : '-'));
361 :     }
362 :     else
363 :     {
364 :     die "Invalid locus $locus";
365 :     }
366 :    
367 :     return ();
368 :     }
369 :    
370 :     sub by_locus
371 :     {
372 :     my ($a, $b) = @_;
373 :    
374 :     my (undef, undef, $A_contig, $A_beg, $A_end, $A_len, $A_strand) = @$a;
375 :     my (undef, undef, $B_contig, $B_beg, $B_end, $B_len, $B_strand) = @$b;
376 :    
377 :     return ( ($A_contig cmp $B_contig)
378 :     || (&FIG::min($A_beg, $A_end) <=> &FIG::min($B_beg, $B_end))
379 :     || ($B_len <=> $A_len)
380 :     || ($A_strand cmp $B_strand)
381 :     );
382 :     }
383 :    
384 :     sub link {
385 :     my($org,$fid) = @_;
386 :    
387 :     return "<a href=diffsF.cgi?org1=$org&gene=$fid>$fid</a>";
388 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3