[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.12, Sat Mar 6 23:54:11 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 18  Line 20 
20  package ProtSims;  package ProtSims;
21    
22  use strict;  use strict;
23    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            if (defined $FIG_Config::ext_bin) {
36                $blat_cmd = "$FIG_Config::ext_bin/blat";
37                $blastall_cmd = "$FIG_Config::ext_bin/blastall";
38                $formatdb_cmd = "$FIG_Config::ext_bin/formatdb";
39            }
40        };
41    }
42    
43  sub blastP {  sub blastP {
44      my($q,$db,$min_hits) = @_;      my($q,$db,$min_hits,$use_blast) = @_;
45    
46      my(@q,@db);      my(@q,@db);
47    
# Line 39  Line 56 
56          $qF = $q;          $qF = $q;
57          if (-e $qF)          if (-e $qF)
58          {          {
             if (! ((-e "$qF.psq") && ((-M "$qF.psq") < (-M $qF))))  
             {  
                 system "$FIG_Config::ext_bin/formatdb -i $qF";  
             }  
59              @q = &gjoseqlib::read_fasta($qF);              @q = &gjoseqlib::read_fasta($qF);
60              $q = \@q;              $q = \@q;
61          }          }
# Line 57  Line 70 
70      {      {
71          $dbF = "/tmp/db.$$";          $dbF = "/tmp/db.$$";
72          &gjoseqlib::print_alignment_as_fasta($dbF,$db);          &gjoseqlib::print_alignment_as_fasta($dbF,$db);
73          system "$FIG_Config::ext_bin/formatdb -i $dbF";          system "$formatdb_cmd -i $dbF";
74      }      }
75      else      else
76      {      {
# Line 66  Line 79 
79          {          {
80              if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))              if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))
81              {              {
82                  system "$FIG_Config::ext_bin/formatdb -i $dbF";                  system "formatdb_cmd -i $dbF";
83              }              }
84              @db = &gjoseqlib::read_fasta($dbF);              @db = &gjoseqlib::read_fasta($dbF);
85              $db = \@db;              $db = \@db;
# Line 78  Line 91 
91      }      }
92    
93      my $tmpF = "/tmp/sim.out.$$";      my $tmpF = "/tmp/sim.out.$$";
94      system "$FIG_Config::ext_bin/blat $dbF $qF -prot -out=blast8 $tmpF > /dev/null";      my $cmd = $use_blast ? "touch $tmpF" : "$blat_cmd $dbF $qF -prot -out=blast8 $tmpF > /dev/null";
95        #print STDERR "$cmd\n";
96        my $rc = system $cmd;
97        if ($rc != 0)
98        {
99            die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";
100        }
101    
102      my @sims1 = map { chomp; [split(/\s+/,$_)] } `cat $tmpF`;      my @sims1 = ();
103        open(BLAT,"<$tmpF") || die "could not open $tmpF";
104        my $sim = <BLAT>;
105        while ($sim && ($sim=~ /^(\S+\t\S+)/))
106        {
107            my $pegs = $1;
108            my @pieces = ();
109            while ($sim && ($sim=~ /^((\S+\t\S+).*\S)/) && ($2 eq $pegs))
110            {
111                push(@pieces,[split(/\t/,$1)]);
112                $sim = <BLAT>;
113            }
114            push(@sims1,&condense(\@pieces));
115        }
116        close(BLAT);
117      unlink $tmpF;      unlink $tmpF;
118  #    print STDERR &Dumper(\@sims1);      #print STDERR &Dumper(sims1 => \@sims1);
119    
120      my @rerun = ();      my @rerun = ();
121      my @sims  = ();      my @sims  = ();
# Line 105  Line 138 
138  #       print STDERR "hits=$hits min_hits=$min_hits\n";  #       print STDERR "hits=$hits min_hits=$min_hits\n";
139          if ($hits >= $min_hits)          if ($hits >= $min_hits)
140          {          {
141              push(@sims,@sims1[$simI..$simI1]);              push(@sims,@sims1[$simI..$simI1-1]);
142          }          }
143          else          else
144          {          {
# Line 114  Line 147 
147          $simI = $simI1;          $simI = $simI1;
148          $qI++;          $qI++;
149      }      }
150  #    print STDERR &Dumper(\@rerun);      #print STDERR &Dumper(rerun => \@rerun);
151    
152      if (@rerun > 0)      if (@rerun > 0)
153      {      {
154          my $tmpQ = "/tmp/tmpQ.$$";          my $tmpQ = "/tmp/tmpQ.$$";
155          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
156          my @blastout = map { chomp; [split(/\s+/,$_)] }          my $cmd = "$blastall_cmd -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";
157                         `$FIG_Config::ext_bin/blastall -m 8 -i $tmpQ -d $dbF -FF -p blastp`;          #print STDERR "$cmd\n";
158            open(BL, "$cmd|") or die "ProtSims::blastP: pipe to blast failed with $!: $cmd\n";
159            my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
160            close(BL);
161          push(@sims,@blastout);          push(@sims,@blastout);
162          unlink $tmpQ;          unlink $tmpQ;
163      }      }
164    
165      my %lnQ   = map { $_->[0] => length($_->[2]) } @q;      my %lnQ   = map { $_->[0] => length($_->[2]) } @$q;
166      my %lnDB  = map { $_->[0] => length($_->[2]) } @db;      my %lnDB  = map { $_->[0] => length($_->[2]) } @$db;
167      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;
168    
169      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }
170      if ($dbF     eq "/tmp/db.$$")      { unlink $dbF }      if ($dbF     eq "/tmp/db.$$")      { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }
171      return @sims;      return sort { ($a->id1 cmp $b->id1) or ($a->psc <=> $b->psc) or ($a->id2 cmp $b->id2) } @sims;
172    }
173    
174    sub condense {
175        my($pieces) = @_;
176    
177        my $best_sc = $pieces->[0]->[10];
178        my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;
179        while (@sims > 1)
180        {
181            my $gap1 = $sims[1]->[6] - $sims[0]->[7];
182            my $gap2 = $sims[1]->[8] - $sims[0]->[9];
183            my $diff = abs($gap1 - $gap2);
184            if (($gap1 <= 300) && ($gap2 <= 300) && ($diff <= (0.1 * &min($gap1,$gap2))))
185            {
186                $sims[0]->[7] = $sims[1]->[7];
187                $sims[0]->[9] = $sims[1]->[9];
188            }
189            splice(@sims,1,1);
190        }
191        $sims[0]->[10] = $best_sc;
192        return $sims[0];
193    }
194    
195    sub min {
196        my($x,$y) = @_;
197        return ($x <= $y) ? $x : $y;
198  }  }
199    
200  1;  1;

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3