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

Annotation of /FigWebServices/diffsF.cgi

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (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.6 push @both, [ $changed, $fid1, $link1, $len1, $fid2, $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.6 push @just1, [ $fid1, $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.6 push @just2, [ $fid2, $link2, $org2T->{$key}->[LEN] ];
202 : overbeek 1.1 }
203 :     }
204 : overbeek 1.3
205 : overbeek 1.4 push @$html, map {
206 : overbeek 1.6 $align = qq(<A HREF="diffsF.cgi?align1=$_->[1]&align2=$_->[4]&org1=$org1&org2=$org2">align</A>);
207 :     $cgi->h2("$align $_->[0]: $_->[2] ($_->[3] bp) => $_->[5] ($_->[6] bp), diff=$_->[7]")
208 : overbeek 1.4 }
209 :     sort {
210 : overbeek 1.6 abs($b->[7]) <=> abs($a->[7])
211 : overbeek 1.4 } @both;
212 :    
213 :    
214 :     push @$html, map {
215 : overbeek 1.6 $cgi->h2("just2: $_->[1] ($_->[2] bp)")
216 : overbeek 1.4 }
217 :     sort {
218 : overbeek 1.6 $b->[2] <=> $a->[2]
219 : overbeek 1.4 } @just2;
220 :    
221 :    
222 :     push @$html, map {
223 : overbeek 1.6 $cgi->h2("just1: $_->[1] ($_->[2] bp)")
224 : overbeek 1.4 }
225 :     sort {
226 : overbeek 1.6 $b->[2] <=> $a->[2]
227 : overbeek 1.4 } @just1;
228 :     }
229 :    
230 :     sub show_alignment {
231 :     my ($fig, $cgi, $html, $gene1, $gene2, $org1, $org2) = @_;
232 :    
233 : overbeek 1.7 my @neighbors;
234 :     my %neighbors;
235 :     my $neighbors;
236 :     if (-s "$org2/neighbors") {
237 :     @neighbors = sort { $a <=> $b} map { m/^(\d+\.\d+)/o ? ($1) : () } `cat $org2/neighbors`;
238 :     %neighbors = map { $_ => 1 } @neighbors;
239 :     $neighbors = join(" ", @neighbors);
240 :     @neighbors ? push(@$html, $cgi->h2("Close neighbors are: " . join(" ", @neighbors)))
241 :     : push(@$html, $cgi->h2("No close nieghbors found")) && return
242 :     ;
243 :     }
244 :     else {
245 :     push(@$html, $cgi->h2("Close neighbors file $org2/neighbors does not exist"))
246 :     && return;
247 :     }
248 :    
249 :    
250 :     my $md5exec = (($_ = `which md5sum`) !~ m/^no/) ? "md5sum"
251 :     : (($_ = `which md5`) !~ m/^no/) ? "md5"
252 :     : push(@$html, $cgi->h2("Could not find an MD5 checksum tool")) && return
253 :     ;
254 :     my $chksum = `echo $neighbors | $md5exec`; chomp $chksum;
255 :     my $tmp_close = "$FIG_Config::temp/tmp_close_$chksum";
256 :    
257 :     my ($size, $mod);
258 :     if ((not ($size = (-s $tmp_close))) || ($mod = (-M $tmp_close > 0.1))) {
259 :     push @$html, $cgi->h2("Rebuilding DB $tmp_close, size=$size, mod=$mod");
260 :    
261 :     system("rm -f $tmp_close");
262 :     foreach my $org (@neighbors) {
263 :     push @$html, $cgi->h2("Appending fasta for $org to $tmp_close");
264 :     system("cat $FIG_Config::organisms/$org/Features/peg/fasta >> $tmp_close")
265 :     && ((push @$html
266 :     , $cgi->h2("Could not append $FIG_Config::organisms/$org/Features/peg/fasta to $tmp_close"))
267 :     && return);
268 :     }
269 :     system("$FIG_Config::ext_bin/formatdb -i $tmp_close")
270 :     && ((push @$html, $cgi->h2("Could not run $FIG_Config::ext_bin/formatdb on $tmp_close, err=$?"))
271 :     && return);
272 :     }
273 :    
274 :     my $tmp_fid = "$FIG_Config::temp/tmp_fid_$$";
275 :    
276 :     open(TMP_FID, ">$tmp_fid")
277 :     || ((push @$html, $cgi->h2("Could not write-open $tmp_fid")) && &cleanup($tmp_fid, $tmp_seqs) && return);
278 :     print TMP_FID "$gene1\n";
279 :     $seq_gene1 = `$FIG_Config::bin/pull_fasta_entries $org1/Features/peg/fasta < $tmp_fid`;
280 :     $seq_gene1 =~ s/^>fig/>old/;
281 :    
282 :     open(TMP_FID, ">$tmp_fid")
283 :     || ((push @$html, $cgi->h2("Could not write-open $tmp_fid")) && &cleanup($tmp_fid, $tmp_seqs) && return);
284 :     print TMP_FID "$gene2\n";
285 :     $seq_gene2 = `$FIG_Config::bin/pull_fasta_entries $org2/Features/peg/fasta < $tmp_fid`;
286 :     $seq_gene2 =~ s/^>fig/>new/;
287 : overbeek 1.4
288 : overbeek 1.7 my ($figfam, $fam_id, @members);
289 :     my @tmp = grep { m/^(\S+)/o && ($1 eq $gene2) } `cat $org2/found`;
290 :     if (not @tmp) {
291 :     push @$html, $cgi->h2("No family found for gene2=$gene2")
292 :     }
293 :     else {
294 : overbeek 1.4 if ($tmp[0] =~ m/^\S+\t(FIG\d{6})/o) {
295 :     $figfam = new FigFam($fig, $1, $default_dir);
296 : overbeek 1.7 $fam_id = $figfam->family_id();
297 : overbeek 1.6 push @$html, $cgi->h2("$gene2 is in Family $fam_id");
298 : overbeek 1.7 @members = $figfam->list_members();
299 :     }
300 :     else {
301 :     $tmp[0] =~ s/\t/ /go;
302 :     push(@$html, $cgi->h2("Could not parse family from entry: $tmp[0]"))
303 :     && &cleanup($tmp_fid, $tmp_seqs) && return;
304 :     }
305 :     }
306 :    
307 :     if (@members) {
308 :     push @$html, $cgi->h2(qq($fam_id has ) . (scalar @members) . qq( members));
309 :     if (@members > 15) {
310 :     @members = grep { $neighbors{$fig->genome_of($_)} } @members;
311 :     my $s = (@members == 1) ? qq() : qq(s);
312 :     my $were = (@members == 1) ? qq(was from a) : qq(were from);
313 :     push @$html, $cgi->h2((scalar @members) . qq( member$s of $fam_id $were close neighbor$s));
314 :     }
315 :     }
316 :     else {
317 :     push @$html, $cgi->h2("No family members found --- searching with BLASTP");
318 :     if (@members = &find_top_blast_hits($fig, $cgi, $html, $seq_gene2, $tmp_close)) {
319 :     if (@members > 15) {
320 :     @members = grep { $neighbors{$fig->genome_of($_)} } @members;
321 :     my $s = (@members == 1) ? qq() : qq(s);
322 :     my $were = (@members == 1) ? qq(was from a) : qq(were from);
323 :     push @$html, $cgi->h2((scalar @members) . qq( BLASTP hit$s $were close neighbor$s));
324 : overbeek 1.4 }
325 :     }
326 :     else {
327 : overbeek 1.7 push( @$html, $cgi->h2("No BLAST sims for $gene2 against neighbors") )
328 :     && &cleanup($tmp_fid, $tmp_seqs) && return;
329 : overbeek 1.4 }
330 :     }
331 : overbeek 1.7
332 :     my $tmp_seqs = "$FIG_Config::temp/tmp_seqs_$$";
333 :     open(TMP_SEQS, ">$tmp_seqs")
334 :     || ((push @$html, $cgi->h2("Could not write-open $tmp_seqs")) && &cleanup($tmp_fid, $tmp_seqs) && return);
335 :     foreach my $id (@members) {
336 :     my $seq = $fig->get_translation($id);
337 :     if ($seq) {
338 :     # push @$html, $cgi->h2("Got seq for $id");
339 :     &FIG::display_id_and_seq($id, \$seq, \*TMP_SEQS);
340 :     }
341 :     else {
342 :     push @$html, $cgi->h2("No seq for $id");
343 :     }
344 : overbeek 1.4 }
345 : overbeek 1.7
346 :     print TMP_SEQS $seq_gene1;
347 :     print TMP_SEQS $seq_gene2;
348 :    
349 :     close(TMP_SEQS)
350 :     || ((push @$html, $cgi->h2("Could not close $tmp_seqs"))
351 :     && &cleanup($tmp_fid, $tmp_seqs) && return);
352 :    
353 :     system("$FIG_Config::ext_bin/clustalw -infile=$tmp_seqs -align -outorder=aligned > /dev/null")
354 :     && ((push @$html, $cgi->h2("clustalw failed on file $tmp_seqs"))
355 :     && &cleanup($tmp_fid, $tmp_seqs) && return);
356 :    
357 :     push @$html, qq(<DIV><FONT FACE="Courier">\n);
358 :     @tmp = `cat $tmp_seqs.aln`;
359 :     @tmp = map { &fid_link($_) } @tmp;
360 :     push @$html, (@tmp ? (@tmp) : ("No lines read from $tmp_seqs.aln"));
361 :     push @$html, qq(</FONT></DIV>\n);
362 :    
363 :     &cleanup($tmp_fid, $tmp_seqs);
364 :     }
365 :    
366 :     sub cleanup {
367 :     my ($tmp_seqs, $tmp_fid) =@_;
368 :     return if $ENV{DEBUG};
369 :    
370 :     system("rm -f $tmp_fid $tmp_seqs*")
371 :     && ((push @$html, $cgi->h2("Could not remove temporary files"))
372 :     && return);
373 :     }
374 :    
375 :     sub find_top_blast_hits {
376 :     my ($fig, $cgi, $html, $query_seq, $tmp_close) = @_;
377 :    
378 :     my $tmp_query = "$FIG_Config::temp/tmp_query_$$";
379 :     open( TMP_QUERY, ">$tmp_query" )
380 :     || ((push @$html, $cgi->h2("Could not write-open $tmp_query")) && return);
381 :     print TMP_QUERY $query_seq;
382 :     close(TMP_QUERY)
383 :     || ((push @$html, $cgi->h2("Could not close $tmp_query")) && return);
384 :    
385 :     my @hits = `$FIG_Config::ext_bin/blastall -i $tmp_query -d $tmp_close -p blastp -m8 -FF -e1.0e-10 | cut -f2`;
386 :     chomp @hits;
387 :    
388 :     push @$html, $cgi->h2(qq(BLAST found ) . (scalar @hits) . qq( hits\n));
389 :     return @hits;
390 : overbeek 1.1 }
391 :    
392 :     sub load_tbl1
393 :     {
394 :     my ($dir) = @_;
395 :     my ($entry, $id, $locus, $contig, $beg, $end, $len, $strand, $taxid, $type);
396 :    
397 :     open(TBL, "cat $dir/Features/*/tbl |") || die "Could not open $dir";
398 :    
399 :     my $tbl = {};
400 :     while (defined($entry = <TBL>))
401 :     {
402 :     chomp $entry;
403 :    
404 :     ($id, $locus) = split /\t/, $entry;
405 :     $id =~ m/^[^\|]+\|(\d+\.\d+)\.([^\.]+)/;
406 :     ($taxid, $type) = ($1, $2);
407 :    
408 :     if ((($contig, $beg, $end, $len, $strand) = &from_locus($locus))
409 :     && defined($contig) && $contig
410 :     && defined($beg) && $beg
411 :     && defined($end) && $end
412 :     && defined($len) && $len
413 :     && defined($strand) && $strand
414 :     )
415 :     {
416 :     $tbl->{"$contig\t$strand$end"} = [ $id, $locus, $contig, $beg, $end, $len, $strand, $type, $taxid ];
417 :     }
418 :     else
419 :     {
420 :     warn "INVALID ENTRY:\t$entry\n";
421 :     }
422 :     }
423 :    
424 :     return $tbl;
425 :     }
426 :    
427 :     sub load_tbl2
428 :     {
429 :     my ($dir,$gene) = @_;
430 :    
431 : overbeek 1.3 my ($entry, $id, $locus, $contig, $beg, $end, $len, $strand, $taxid, $type, $keep);
432 : overbeek 1.1
433 :     open(TBL, "cat $dir/Features/*/tbl |") || die "Could not open $dir";
434 :    
435 :     my $tbl = {};
436 :     while (defined($entry = <TBL>))
437 :     {
438 :     chomp $entry;
439 :    
440 :     ($id, $locus) = split /\t/, $entry;
441 :     $id =~ m/^[^\|]+\|(\d+\.\d+)\.([^\.]+)/;
442 :     ($taxid, $type) = ($1, $2);
443 :    
444 :     if ((($contig, $beg, $end, $len, $strand) = &from_locus($locus))
445 :     && defined($contig) && $contig
446 :     && defined($beg) && $beg
447 :     && defined($end) && $end
448 :     && defined($len) && $len
449 :     && defined($strand) && $strand
450 :     )
451 :     {
452 : overbeek 1.3 push(@{$tbl->{$contig}}, [ $id, $locus, $contig, $beg, $end, $len, $strand, $type, $taxid ]);
453 : overbeek 1.1 if ($id eq $gene) { $keep = $contig }
454 :     }
455 :     else
456 :     {
457 :     warn "INVALID ENTRY:\t$entry\n";
458 :     }
459 :     }
460 : overbeek 1.5 $keep || die "no contig for gene=$gene";
461 : overbeek 1.1 return sort { by_locus($a,$b) } @{$tbl->{$keep}};
462 :     }
463 :    
464 :     sub from_locus
465 :     {
466 :     my ($locus) = @_;
467 :     my ($contig, $beg, $end);
468 :    
469 :     if ( ($locus =~ m/^([^,]+)_(\d+)_\d+/) && ($contig = $1) && ($beg = $2)
470 :     && ($locus =~ m/[^,]+_\d+_(\d+)$/) && ($end = $1)
471 :     )
472 :     {
473 :     return ($contig, $beg, $end, (1+abs($end-$beg)), (($beg < $end) ? '+' : '-'));
474 :     }
475 :     else
476 :     {
477 :     die "Invalid locus $locus";
478 :     }
479 :    
480 :     return ();
481 :     }
482 :    
483 :     sub by_locus
484 :     {
485 :     my ($a, $b) = @_;
486 :    
487 :     my (undef, undef, $A_contig, $A_beg, $A_end, $A_len, $A_strand) = @$a;
488 :     my (undef, undef, $B_contig, $B_beg, $B_end, $B_len, $B_strand) = @$b;
489 :    
490 :     return ( ($A_contig cmp $B_contig)
491 :     || (&FIG::min($A_beg, $A_end) <=> &FIG::min($B_beg, $B_end))
492 :     || ($B_len <=> $A_len)
493 :     || ($A_strand cmp $B_strand)
494 :     );
495 :     }
496 :    
497 :     sub link {
498 :     my($org,$fid) = @_;
499 :    
500 :     return "<a href=diffsF.cgi?org1=$org&gene=$fid>$fid</a>";
501 :     }
502 : overbeek 1.7
503 :     sub fid_link {
504 :     my($line) = @_;
505 :    
506 :     chomp $line;
507 :     if ($line =~ m/^(fig\S+)\b(.*)$/o) {
508 :     my ($id, $rest) = ($1, $2);
509 :     $rest =~ s/\s/\&nbsp\;/go;
510 :     $line = qq(<A HREF="http://anno-3.nmpdr.org/anno/FIG/protein.cgi?prot=$id">$id</A>$rest);
511 :     }
512 :     else {
513 :     $line =~ s/\s/\&nbsp\;/go;
514 :     }
515 :    
516 :     return $line . qq(<BR>\n);
517 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3