[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.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 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            $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 39  Line 54 
54          $qF = $q;          $qF = $q;
55          if (-e $qF)          if (-e $qF)
56          {          {
             if (! ((-e "$qF.psq") && ((-M "$qF.psq") < (-M $qF))))  
             {  
                 system "$FIG_Config::ext_bin/formatdb -i $qF";  
             }  
57              @q = &gjoseqlib::read_fasta($qF);              @q = &gjoseqlib::read_fasta($qF);
58              $q = \@q;              $q = \@q;
59          }          }
# Line 57  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 66  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 78  Line 89 
89      }      }
90    
91      my $tmpF = "/tmp/sim.out.$$";      my $tmpF = "/tmp/sim.out.$$";
92      system "$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";
94        my $rc = system $cmd;
95        if ($rc != 0)
96        {
97            die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";
98        }
99    
100      my @sims1 = map { chomp; [split(/\s+/,$_)] } `cat $tmpF`;      my @sims1 = ();
101        open(BLAT,"<$tmpF") || die "could not open $tmpF";
102        my $sim = <BLAT>;
103        while ($sim && ($sim=~ /^(\S+\t\S+)/))
104        {
105            my $pegs = $1;
106            my @pieces = ();
107            while ($sim && ($sim=~ /^((\S+\t\S+).*\S)/) && ($2 eq $pegs))
108            {
109                push(@pieces,[split(/\t/,$1)]);
110                $sim = <BLAT>;
111            }
112            push(@sims1,&condense(\@pieces));
113        }
114        close(BLAT);
115      unlink $tmpF;      unlink $tmpF;
116  #    print STDERR &Dumper(\@sims1);      #print STDERR &Dumper(sims1 => \@sims1);
117    
118      my @rerun = ();      my @rerun = ();
119      my @sims  = ();      my @sims  = ();
# Line 105  Line 136 
136  #       print STDERR "hits=$hits min_hits=$min_hits\n";  #       print STDERR "hits=$hits min_hits=$min_hits\n";
137          if ($hits >= $min_hits)          if ($hits >= $min_hits)
138          {          {
139              push(@sims,@sims1[$simI..$simI1]);              push(@sims,@sims1[$simI..$simI1-1]);
140          }          }
141          else          else
142          {          {
# Line 114  Line 145 
145          $simI = $simI1;          $simI = $simI1;
146          $qI++;          $qI++;
147      }      }
148  #    print STDERR &Dumper(\@rerun);      #print STDERR &Dumper(rerun => \@rerun);
149    
150      if (@rerun > 0)      if (@rerun > 0)
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 @blastout = map { chomp; [split(/\s+/,$_)] }          my $cmd = "$blastall_cmd -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";
155                         `$FIG_Config::ext_bin/blastall -m 8 -i $tmpQ -d $dbF -FF -p blastp`;          #print STDERR "$cmd\n";
156            open(BL, "$cmd|") or die "ProtSims::blastP: pipe to blast failed with $!: $cmd\n";
157            my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
158            close(BL);
159          push(@sims,@blastout);          push(@sims,@blastout);
160          unlink $tmpQ;          unlink $tmpQ;
161      }      }
162    
163      my %lnQ   = map { $_->[0] => length($_->[2]) } @q;      my %lnQ   = map { $_->[0] => length($_->[2]) } @$q;
164      my %lnDB  = map { $_->[0] => length($_->[2]) } @db;      my %lnDB  = map { $_->[0] => length($_->[2]) } @$db;
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 @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 {
173        my($pieces) = @_;
174    
175        my $best_sc = $pieces->[0]->[10];
176        my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;
177        while (@sims > 1)
178        {
179            my $gap1 = $sims[1]->[6] - $sims[0]->[7];
180            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];
185                $sims[0]->[9] = $sims[1]->[9];
186            }
187            splice(@sims,1,1);
188        }
189        $sims[0]->[10] = $best_sc;
190        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.2  
changed lines
  Added in v.1.10

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3