[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.5, Sun Aug 23 01:20:23 2009 UTC revision 1.13, Mon Aug 30 16:09:00 2010 UTC
# Line 1  Line 1 
1    
2    # This is a SAS component
3    
4  #  #
5  # Copyright (c) 2003-2006 University of Chicago and Fellowship  # Copyright (c) 2003-2006 University of Chicago and Fellowship
6  # for Interpretations of Genomes. All Rights Reserved.  # for Interpretations of Genomes. All Rights Reserved.
# Line 21  Line 23 
23  use Sim;  use Sim;
24  use Data::Dumper;  use Data::Dumper;
25  use Carp;  use Carp;
 use FIG_Config;  
26  use gjoseqlib;  use gjoseqlib;
27    
28    use SeedAware;
29    
30    my $blat_cmd = SeedAware::executable_for("blat");
31    my $blastall_cmd = SeedAware::executable_for("blastall");
32    my $formatdb_cmd = SeedAware::executable_for("formatdb");
33    
34  sub blastP {  sub blastP {
35      my($q,$db,$min_hits) = @_;      my($q,$db,$min_hits,$use_blast) = @_;
36    
37      my(@q,@db);      my(@q,@db);
38    
39        my $tmp_dir = SeedAware::location_of_tmp();
40    
41      my $qF;      my $qF;
42      if (ref $q)      if (ref $q)
43      {      {
44          $qF = "/tmp/query.$$";          $qF = "$tmp_dir/query.$$";
45          &gjoseqlib::print_alignment_as_fasta($qF,$q);          &gjoseqlib::print_alignment_as_fasta($qF,$q);
46      }      }
47      else      else
# Line 52  Line 61 
61      my $dbF;      my $dbF;
62      if (ref $db)      if (ref $db)
63      {      {
64          $dbF = "/tmp/db.$$";          $dbF = "$tmp_dir/db.$$";
65          &gjoseqlib::print_alignment_as_fasta($dbF,$db);          &gjoseqlib::print_alignment_as_fasta($dbF,$db);
66          system "$FIG_Config::ext_bin/formatdb -i $dbF";          system($formatdb_cmd, '-i', $dbF);
67      }      }
68      else      else
69      {      {
# Line 63  Line 72 
72          {          {
73              if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))              if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))
74              {              {
75                  system "$FIG_Config::ext_bin/formatdb -i $dbF";                  system($formatdb_cmd, '-i', $dbF);;
76              }              }
77              @db = &gjoseqlib::read_fasta($dbF);              @db = &gjoseqlib::read_fasta($dbF);
78              $db = \@db;              $db = \@db;
# Line 74  Line 83 
83          }          }
84      }      }
85    
86      my $tmpF = "/tmp/sim.out.$$";      my $tmpF = "$tmp_dir/sim.out.$$";
87      my $cmd = "$FIG_Config::ext_bin/blat $dbF $qF -prot -out=blast8 $tmpF > /dev/null";  
88      #print STDERR "$cmd\n";      if ($use_blast)
89      my $rc = system $cmd;      {
90            open(my $fh, ">", $tmpF);
91            close($fh);
92        }
93        else
94        {
95            my @cmd = ($blat_cmd, $dbF, $qF, "-prot", "-out=blast8", $tmpF);
96    
97            my $rc = system_with_redirect(\@cmd, { stdout => '/dev/null' } );
98            print STDERR "Blat returns $rc: @cmd\n";
99      if ($rc != 0)      if ($rc != 0)
100      {      {
101          die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";              warn "Blat run failed with rc=$rc\n";
102                open(my $fh, ">", $tmpF);
103                close($fh);
104            }
105      }      }
106    
107    #     my $cmd = $use_blast ? "touch $tmpF" : "$blat_cmd $dbF $qF -prot -out=blast8 $tmpF > /dev/null";
108    #     #print STDERR "$cmd\n";
109    #     my $rc = system $cmd;
110    #     if ($rc != 0)
111    #     {
112    #       die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";
113    #     }
114    
115      my @sims1 = ();      my @sims1 = ();
116      open(BLAT,"<$tmpF") || die "could not open $tmpF";      open(BLAT,"<$tmpF") || die "could not open $tmpF";
117      my $sim = <BLAT>;      my $sim = <BLAT>;
# Line 135  Line 164 
164    
165      if (@rerun > 0)      if (@rerun > 0)
166      {      {
167          my $tmpQ = "/tmp/tmpQ.$$";          my $tmpQ = "$tmp_dir/tmpQ.$$";
168          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
169          my $cmd = "$FIG_Config::ext_bin/blastall -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";          # my $cmd = "$blastall_cmd -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";
170            my @cmd = ($blastall_cmd,  '-m', 8, '-i', $tmpQ, '-d', $dbF, '-FF', '-p', 'blastp', '-e', '1e-5');
171            print STDERR "@cmd\n";
172          #print STDERR "$cmd\n";          #print STDERR "$cmd\n";
173          open(BL, "$cmd|") or die "ProtSims::blastP: pipe to blast failed with $!: $cmd\n";          #open(BL, "$cmd|") or die "ProtSims::blastP: pipe to blast failed with $!: $cmd\n";
174    
175            open(BL, "-|", @cmd) or die "ProtSims::blastP: pipe to blast failed with $!: @cmd\n";
176    
177          my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;          my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
178          close(BL);          close(BL);
179          push(@sims,@blastout);          push(@sims,@blastout);
# Line 150  Line 184 
184      my %lnDB  = map { $_->[0] => length($_->[2]) } @$db;      my %lnDB  = map { $_->[0] => length($_->[2]) } @$db;
185      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;
186    
187      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }      if ($qF      eq "$tmp_dir/query.$$")   { unlink $qF  }
188      if ($dbF     eq "/tmp/db.$$")      { unlink $dbF }      if ($dbF     eq "$tmp_dir/db.$$")      { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }
189      return sort { ($a->id1 cmp $b->id1) or ($a->id2 cmp $b->id2) or ($b->bsc <=> $a->bsc) } @sims;      return sort { ($a->id1 cmp $b->id1) or ($a->psc <=> $b->psc) or ($a->id2 cmp $b->id2) } @sims;
190  }  }
191    
192  sub condense {  sub condense {
193      my($pieces) = @_;      my($pieces) = @_;
194    
195      my $best_sc = $pieces->[0]->[9];      my $best_sc = $pieces->[0]->[10];
196      my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;      my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;
197      while ((@sims > 1) &&      while (@sims > 1)
198             ($sims[1]->[6] <= ($sims[0]->[7] + 100)) &&      {
199             ($sims[1]->[8] <= ($sims[0]->[9] + 100)) &&          my $gap1 = $sims[1]->[6] - $sims[0]->[7];
200             (abs(($sims[1]->[6] - $sims[0]->[7]) - ($sims[1]->[8] - $sims[0]->[9])) < 20))          my $gap2 = $sims[1]->[8] - $sims[0]->[9];
201            my $diff = abs($gap1 - $gap2);
202            if (($gap1 <= 300) && ($gap2 <= 300) && ($diff <= (0.1 * &min($gap1,$gap2))))
203      {      {
204          $sims[0]->[7] = $sims[1]->[7];          $sims[0]->[7] = $sims[1]->[7];
205          $sims[0]->[9] = $sims[1]->[9];          $sims[0]->[9] = $sims[1]->[9];
206            }
207          splice(@sims,1,1);          splice(@sims,1,1);
208      }      }
209        $sims[0]->[10] = $best_sc;
210      return $sims[0];      return $sims[0];
211  }  }
212    
213    sub min {
214        my($x,$y) = @_;
215        return ($x <= $y) ? $x : $y;
216    }
217    
218  1;  1;

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3