[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.6, Wed Oct 18 09:12:24 2006 UTC revision 1.7, Wed Oct 18 23:12:25 2006 UTC
# Line 230  Line 230 
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    
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    
288      my (@tmp, $tmp);      my ($figfam, $fam_id, @members);
289      if (@tmp = grep { m/^(\S+)/o && ($1 eq $gene2) } `cat $org2/found`) {      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 $fam_id = $figfam->family_id();              $fam_id = $figfam->family_id();
297              push @$html, $cgi->h2("$gene2 is in Family $fam_id");              push @$html, $cgi->h2("$gene2 is in Family $fam_id");
298  #           my @members = grep { $neighbors{$fig->genome_of($_)} } $figfam->list_members();              @members = $figfam->list_members();
299              my @members = $figfam->list_members();          }
300              @members || push @$html, $cgi->h2("No members in $fam_id");          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                }
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_$$";              my $tmp_seqs = "$FIG_Config::temp/tmp_seqs_$$";
             my $tmp_fid  = "$FIG_Config::temp/tmp_fid_$$";  
333              open(TMP_SEQS, ">$tmp_seqs")              open(TMP_SEQS, ">$tmp_seqs")
334                  || ((push @$html, $cgi->h2("Could not write-open $tmp_seqs")) && return);          || ((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                  if ($seq) {                  if ($seq) {
# Line 258  Line 343 
343                  }                  }
344              }              }
345    
346              open(TMP_FID, ">$tmp_fid")      print TMP_SEQS $seq_gene1;
347                  || ((push @$html, $cgi->h2("Could not write-open $tmp_fid")) && return);      print TMP_SEQS $seq_gene2;
             print TMP_FID "$gene1\n";  
             $tmp =  `$FIG_Config::bin/pull_fasta_entries $org1/Features/peg/fasta < $tmp_fid`;  
             $tmp =~ s/^>fig/>old/;  
             print TMP_SEQS $tmp;  
   
             open(TMP_FID, ">$tmp_fid")  
                 || ((push @$html, $cgi->h2("Could not write-open $tmp_fid")) && return);  
             print TMP_FID "$gene2\n";  
             $tmp = `$FIG_Config::bin/pull_fasta_entries $org2/Features/peg/fasta < $tmp_fid`;  
             $tmp =~ s/^>fig/>new/;  
             print TMP_SEQS $tmp;  
348    
349              close(TMP_SEQS)              close(TMP_SEQS)
350                  || ((push @$html, $cgi->h2("Could not close $tmp_seqs")) && return);          || ((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")              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")) && return);          && ((push @$html, $cgi->h2("clustalw failed on file $tmp_seqs"))
355                && &cleanup($tmp_fid, $tmp_seqs) && return);
356    
357              push @$html, qq(<pre>\n);      push @$html, qq(<DIV><FONT FACE="Courier">\n);
358              @tmp = `cat $tmp_seqs.aln`;              @tmp = `cat $tmp_seqs.aln`;
359        @tmp = map { &fid_link($_) } @tmp;
360              push @$html, (@tmp ? (@tmp) : ("No lines read from $tmp_seqs.aln"));              push @$html, (@tmp ? (@tmp) : ("No lines read from $tmp_seqs.aln"));
361              push @$html, qq(</pre>\n);      push @$html, qq(</FONT></DIV>\n);
362    
363              system("rm -f $tmp_seqs* $tmp_fid")      &cleanup($tmp_fid, $tmp_seqs);
                 && ((push @$html, $cgi->h2("Could not remove temporary files")) && return);  
364          }          }
365          else {  
366              $tmp[0] =~ s/\t/ /go;  sub cleanup {
367              push @$html, $cgi->h2("Could not parse family from entry: $tmp[0]");      my ($tmp_seqs, $tmp_fid) =@_;
368          }      return if $ENV{DEBUG};
369      }  
370      else {      system("rm -f $tmp_fid $tmp_seqs*")
371          push @$html, $cgi->h2("No family found for gene2=$gene2");          && ((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 406  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.6  
changed lines
  Added in v.1.7

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3