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

Diff of /FigWebServices/diffsF.cgi

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.5, Wed Oct 18 07:38:44 2006 UTC revision 1.7, Wed Oct 18 23:12:25 2006 UTC
# Line 178  Line 178 
178                      $changed = qq(shortened);                      $changed = qq(shortened);
179                  }                  }
180    
181                  push @both, [ $changed, $link1, $len1, $link2, $len2, ($len2 - $len1) ];                  push @both, [ $changed, $fid1, $link1, $len1, $fid2, $link2, $len2, ($len2 - $len1) ];
182              }              }
183          }          }
184          else          else
185          {          {
186              $fid1 = $org1T->{$key}->[FID];              $fid1 = $org1T->{$key}->[FID];
187              $link1 = &link($org1,$fid1);              $link1 = &link($org1,$fid1);
188              push @just1, [ $link1, $org1T->{$key}->[LEN] ];              push @just1, [ $fid1, $link1, $org1T->{$key}->[LEN] ];
189          }          }
190      }      }
191    
# Line 198  Line 198 
198              $locus2 = $org2T->{$key}->[LOCUS];              $locus2 = $org2T->{$key}->[LOCUS];
199              $fid2 = $org2T->{$key}->[FID];              $fid2 = $org2T->{$key}->[FID];
200              $link2 = &link($org2,$fid2);              $link2 = &link($org2,$fid2);
201              push @just2, [ $link2, $org2T->{$key}->[LEN] ];              push @just2, [ $fid2, $link2, $org2T->{$key}->[LEN] ];
202          }          }
203      }      }
204    
205      push @$html, map {      push @$html, map {
206          $align = qq(<A HREF="diffsF.cgi?align1=$_->[1]&align2=$_->[3]&org1=$org1&org2=$org2">align</A>);          $align = qq(<A HREF="diffsF.cgi?align1=$_->[1]&align2=$_->[4]&org1=$org1&org2=$org2">align</A>);
207          $cgi->h2("$align $_->[0]: $_->[1] ($_->[2] bp) => $_->[3] ($_->[4] bp), diff=$_->[5]")          $cgi->h2("$align $_->[0]: $_->[2] ($_->[3] bp) => $_->[5] ($_->[6] bp), diff=$_->[7]")
208          }          }
209      sort {      sort {
210          abs($b->[5]) <=> abs($a->[5])          abs($b->[7]) <=> abs($a->[7])
211          } @both;          } @both;
212    
213    
214      push @$html, map {      push @$html, map {
215          $cgi->h2("just2: $_->[0] ($_->[1] bp)")          $cgi->h2("just2: $_->[1] ($_->[2] bp)")
216          }          }
217      sort {      sort {
218          $b->[1] <=> $a->[1]          $b->[2] <=> $a->[2]
219          } @just2;          } @just2;
220    
221    
222      push @$html, map {      push @$html, map {
223          $cgi->h2("just1: $_->[0] ($_->[1] bp)")          $cgi->h2("just1: $_->[1] ($_->[2] bp)")
224          }          }
225      sort {      sort {
226          $b->[1] <=> $a->[1]          $b->[2] <=> $a->[2]
227          } @just1;          } @just1;
228  }  }
229    
230  sub show_alignment {  sub show_alignment {
231      my ($fig, $cgi, $html, $gene1, $gene2, $org1, $org2) = @_;      my ($fig, $cgi, $html, $gene1, $gene2, $org1, $org2) = @_;
232    
233      my %neighbors = map { m/^(\d+\.\d+)/o ? ($1 => 1) : () } `cat $org2/neighbors`;      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      my (@tmp, $tmp);  
250      if (@tmp = grep { m/^(\S+)/o && ($1 eq $gene2) } `cat $org2/found`) {      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    
288        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          if ($tmp[0] =~ m/^\S+\t(FIG\d{6})/o) {          if ($tmp[0] =~ m/^\S+\t(FIG\d{6})/o) {
295              $figfam = new FigFam($fig, $1, $default_dir);              $figfam = new FigFam($fig, $1, $default_dir);
296              my @members = grep { $neighbors{$fig->genome_of($_)} } $figfam->list_members();              $fam_id = $figfam->family_id();
297                push @$html, $cgi->h2("$gene2 is in Family $fam_id");
298                @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              my $tmp_seqs = "$FIG_Config::temp/tmp.$$";      if (@members) {
308              open(TMP, ">$tmp_seqs")          push @$html, $cgi->h2(qq($fam_id has ) . (scalar @members) . qq( members));
309                  || ((push @$html, $cgi->h2("Could not write-open $tmp_seqs")) && return);          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                }
325            }
326            else {
327                push( @$html, $cgi->h2("No BLAST sims for $gene2 against neighbors") )
328                    && &cleanup($tmp_fid, $tmp_seqs) && return;
329            }
330        }
331    
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) {              foreach my $id (@members) {
336                  my $seq = $fig->get_translation($id);                  my $seq = $fig->get_translation($id);
337                  &FIG::display_id_and_seq($id, \$seq, \*TMP);          if ($seq) {
338              }  #           push @$html, $cgi->h2("Got seq for $id");
339              close(TMP)              &FIG::display_id_and_seq($id, \$seq, \*TMP_SEQS);
                 || ((push @$html, $cgi->h2("Could not close $tmp_seqs")) && return);  
   
             $tmp = $gene1;  
             $tmp =~ s{\|}{\\\|}o;  
             system("echo $tmp1 | pull_fasta_entries $org1/Features/peg/fasta >> $tmp_seqs")  
                 && ((push @$html, $cgi->h2("Could not append seq of $gene1 to $tmp_seqs")) && return);  
   
             $tmp = $gene2;  
             $tmp =~ s{\|}{\\\|}o;  
             system("echo $tmp1 | pull_fasta_entries $org2/Features/peg/fasta >> $tmp_seqs")  
                 && ((push @$html, $cgi->h2("Could not append seq of $gene2 to $tmp_seqs")) && return);  
   
             system("clustalw -infile=$tmp_seqs -align -outorder=aligned > /dev/null")  
                 && ((push @$html, $cgi->h2("clustalw failed on file $tmp_seqs")) && return);  
   
             push @$html, qq(<pre>\n);  
             push @$html, `cat $tmp_seqs.aln`;  
             push @$html, qq(</pre>\n);  
340          }          }
341          else {          else {
342              $tmp[0] =~ s/\t/ /go;              push @$html, $cgi->h2("No seq for $id");
             push @$html, $cgi->h2("Could not parse family from entry: $tmp[0]");  
343          }          }
344      }      }
345      else {  
346          push @$html, $cgi->h2("No family found for gene2=$gene2");      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  }  }
391    
392  sub load_tbl1  sub load_tbl1
# Line 386  Line 499 
499    
500      return "<a href=diffsF.cgi?org1=$org&gene=$fid>$fid</a>";      return "<a href=diffsF.cgi?org1=$org&gene=$fid>$fid</a>";
501  }  }
502    
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    }

Legend:
Removed from v.1.5  
changed lines
  Added in v.1.7

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3