[Bio] / FigKernelPackages / ProtSims.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/ProtSims.pm

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

revision 1.2, Tue Aug 18 15:16:21 2009 UTC revision 1.4, Fri Aug 21 21:31:01 2009 UTC
# Line 18  Line 18 
18  package ProtSims;  package ProtSims;
19    
20  use strict;  use strict;
21    use Sim;
22  use Data::Dumper;  use Data::Dumper;
23  use Carp;  use Carp;
24  use FIG_Config;  use FIG_Config;
# Line 39  Line 40 
40          $qF = $q;          $qF = $q;
41          if (-e $qF)          if (-e $qF)
42          {          {
             if (! ((-e "$qF.psq") && ((-M "$qF.psq") < (-M $qF))))  
             {  
                 system "$FIG_Config::ext_bin/formatdb -i $qF";  
             }  
43              @q = &gjoseqlib::read_fasta($qF);              @q = &gjoseqlib::read_fasta($qF);
44              $q = \@q;              $q = \@q;
45          }          }
# Line 78  Line 75 
75      }      }
76    
77      my $tmpF = "/tmp/sim.out.$$";      my $tmpF = "/tmp/sim.out.$$";
78      system "$FIG_Config::ext_bin/blat $dbF $qF -prot -out=blast8 $tmpF > /dev/null";      my $cmd = "$FIG_Config::ext_bin/blat $dbF $qF -prot -out=blast8 $tmpF > /dev/null";
79        #print STDERR "$cmd\n";
80        my $rc = system $cmd;
81        if ($rc != 0)
82        {
83            die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";
84        }
85    
86      my @sims1 = map { chomp; [split(/\s+/,$_)] } `cat $tmpF`;      my @sims1 = map { chomp; [split(/\s+/,$_)] } `cat $tmpF`;
87      unlink $tmpF;      unlink $tmpF;
88  #    print STDERR &Dumper(\@sims1);      #print STDERR &Dumper(sims1 => \@sims1);
89    
90      my @rerun = ();      my @rerun = ();
91      my @sims  = ();      my @sims  = ();
# Line 105  Line 108 
108  #       print STDERR "hits=$hits min_hits=$min_hits\n";  #       print STDERR "hits=$hits min_hits=$min_hits\n";
109          if ($hits >= $min_hits)          if ($hits >= $min_hits)
110          {          {
111              push(@sims,@sims1[$simI..$simI1]);              push(@sims,@sims1[$simI..$simI1-1]);
112          }          }
113          else          else
114          {          {
# Line 114  Line 117 
117          $simI = $simI1;          $simI = $simI1;
118          $qI++;          $qI++;
119      }      }
120  #    print STDERR &Dumper(\@rerun);      #print STDERR &Dumper(rerun => \@rerun);
121    
122      if (@rerun > 0)      if (@rerun > 0)
123      {      {
124          my $tmpQ = "/tmp/tmpQ.$$";          my $tmpQ = "/tmp/tmpQ.$$";
125          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
126          my @blastout = map { chomp; [split(/\s+/,$_)] }          my $cmd = "$FIG_Config::ext_bin/blastall -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";
127                         `$FIG_Config::ext_bin/blastall -m 8 -i $tmpQ -d $dbF -FF -p blastp`;          #print STDERR "$cmd\n";
128            open(BL, "$cmd|") or die "ProtSims::blastP: pipe to blast failed with $!: $cmd\n";
129            my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
130            close(BL);
131          push(@sims,@blastout);          push(@sims,@blastout);
132          unlink $tmpQ;          unlink $tmpQ;
133      }      }
134    
135      my %lnQ   = map { $_->[0] => length($_->[2]) } @q;      my %lnQ   = map { $_->[0] => length($_->[2]) } @$q;
136      my %lnDB  = map { $_->[0] => length($_->[2]) } @db;      my %lnDB  = map { $_->[0] => length($_->[2]) } @$db;
137      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;
138    
139      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }
140      if ($dbF     eq "/tmp/db.$$")      { unlink $dbF }      if ($dbF     eq "/tmp/db.$$")      { unlink $dbF }
141      return @sims;      return sort { ($a->id1 cmp $b->id1) or ($b->bsc <=> $a->bsc) or ($a->id2 cmp $b->id2) } @sims;
142  }  }
143    
144  1;  1;

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.4

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3