[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.10, Thu Dec 17 21:35:24 2009 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    my $blat_cmd = "blat";
29    my $blastall_cmd = "blastall";
30    my $formatdb_cmd = "formatdb";
31    
32    BEGIN {
33        eval {
34            require FIG_Config;
35            $blat_cmd = "$FIG_Config::ext_bin/blat";
36            $blastall_cmd = "$FIG_Config::ext_bin/blastall";
37            $formatdb_cmd = "$FIG_Config::ext_bin/formatdb";
38        };
39    }
40    
41  sub blastP {  sub blastP {
42      my($q,$db,$min_hits) = @_;      my($q,$db,$min_hits) = @_;
43    
# Line 54  Line 68 
68      {      {
69          $dbF = "/tmp/db.$$";          $dbF = "/tmp/db.$$";
70          &gjoseqlib::print_alignment_as_fasta($dbF,$db);          &gjoseqlib::print_alignment_as_fasta($dbF,$db);
71          system "$FIG_Config::ext_bin/formatdb -i $dbF";          system "$formatdb_cmd -i $dbF";
72      }      }
73      else      else
74      {      {
# Line 63  Line 77 
77          {          {
78              if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))              if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))
79              {              {
80                  system "$FIG_Config::ext_bin/formatdb -i $dbF";                  system "formatdb_cmd -i $dbF";
81              }              }
82              @db = &gjoseqlib::read_fasta($dbF);              @db = &gjoseqlib::read_fasta($dbF);
83              $db = \@db;              $db = \@db;
# Line 75  Line 89 
89      }      }
90    
91      my $tmpF = "/tmp/sim.out.$$";      my $tmpF = "/tmp/sim.out.$$";
92      my $cmd = "$FIG_Config::ext_bin/blat $dbF $qF -prot -out=blast8 $tmpF > /dev/null";      my $cmd = "$blat_cmd $dbF $qF -prot -out=blast8 $tmpF > /dev/null";
93      #print STDERR "$cmd\n";      #print STDERR "$cmd\n";
94      my $rc = system $cmd;      my $rc = system $cmd;
95      if ($rc != 0)      if ($rc != 0)
# Line 137  Line 151 
151      {      {
152          my $tmpQ = "/tmp/tmpQ.$$";          my $tmpQ = "/tmp/tmpQ.$$";
153          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
154          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";
155          #print STDERR "$cmd\n";          #print STDERR "$cmd\n";
156          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";
157          my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;          my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
# Line 151  Line 165 
165      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;
166    
167      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }
168      if ($dbF     eq "/tmp/db.$$")      { unlink $dbF }      if ($dbF     eq "/tmp/db.$$")      { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }
169      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;
170  }  }
171    
172  sub condense {  sub condense {
173      my($pieces) = @_;      my($pieces) = @_;
174    
175      my $best_sc = $pieces->[0]->[9];      my $best_sc = $pieces->[0]->[10];
176      my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;      my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;
177      while ((@sims > 1) &&      while (@sims > 1)
178             ($sims[1]->[6] <= ($sims[0]->[7] + 100)) &&      {
179             ($sims[1]->[8] <= ($sims[0]->[9] + 100)) &&          my $gap1 = $sims[1]->[6] - $sims[0]->[7];
180             (abs(($sims[1]->[6] - $sims[0]->[7]) - ($sims[1]->[8] - $sims[0]->[9])) < 20))          my $gap2 = $sims[1]->[8] - $sims[0]->[9];
181            my $diff = abs($gap1 - $gap2);
182            if (($gap1 <= 300) && ($gap2 <= 300) && ($diff <= (0.1 * &min($gap1,$gap2))))
183      {      {
184          $sims[0]->[7] = $sims[1]->[7];          $sims[0]->[7] = $sims[1]->[7];
185          $sims[0]->[9] = $sims[1]->[9];          $sims[0]->[9] = $sims[1]->[9];
186            }
187          splice(@sims,1,1);          splice(@sims,1,1);
188      }      }
189        $sims[0]->[10] = $best_sc;
190      return $sims[0];      return $sims[0];
191  }  }
192    
193    sub min {
194        my($x,$y) = @_;
195        return ($x <= $y) ? $x : $y;
196    }
197    
198  1;  1;

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3