[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.1, Tue Aug 18 15:11:51 2009 UTC revision 1.8, Sat Nov 7 18:37:44 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 26  Line 27 
27  sub blastP {  sub blastP {
28      my($q,$db,$min_hits) = @_;      my($q,$db,$min_hits) = @_;
29    
30        my(@q,@db);
31    
32      my $qF;      my $qF;
33      if (ref $q)      if (ref $q)
34      {      {
# Line 37  Line 40 
40          $qF = $q;          $qF = $q;
41          if (-e $qF)          if (-e $qF)
42          {          {
43              if (! ((-e "$qF.psq") && ((-M "$qF.psq") < (-M $qF))))              @q = &gjoseqlib::read_fasta($qF);
             {  
                 system "$FIG_Config::ext_bin/formatdb -i $qF";  
             }  
             my @q = &gjoseqlib::read_fasta($qF);  
44              $q = \@q;              $q = \@q;
45          }          }
46          else          else
# Line 66  Line 65 
65              {              {
66                  system "$FIG_Config::ext_bin/formatdb -i $dbF";                  system "$FIG_Config::ext_bin/formatdb -i $dbF";
67              }              }
68              my @db = &gjoseqlib::read_fasta($dbF);              @db = &gjoseqlib::read_fasta($dbF);
69              $db = \@db;              $db = \@db;
70          }          }
71          else          else
# Line 76  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 = ();
87        open(BLAT,"<$tmpF") || die "could not open $tmpF";
88        my $sim = <BLAT>;
89        while ($sim && ($sim=~ /^(\S+\t\S+)/))
90        {
91            my $pegs = $1;
92            my @pieces = ();
93            while ($sim && ($sim=~ /^((\S+\t\S+).*\S)/) && ($2 eq $pegs))
94            {
95                push(@pieces,[split(/\t/,$1)]);
96                $sim = <BLAT>;
97            }
98            push(@sims1,&condense(\@pieces));
99        }
100        close(BLAT);
101      unlink $tmpF;      unlink $tmpF;
102  #    print STDERR &Dumper(\@sims1);      #print STDERR &Dumper(sims1 => \@sims1);
103    
104      my @rerun = ();      my @rerun = ();
105      my @sims  = ();      my @sims  = ();
# Line 103  Line 122 
122  #       print STDERR "hits=$hits min_hits=$min_hits\n";  #       print STDERR "hits=$hits min_hits=$min_hits\n";
123          if ($hits >= $min_hits)          if ($hits >= $min_hits)
124          {          {
125              push(@sims,@sims1[$simI..$simI1]);              push(@sims,@sims1[$simI..$simI1-1]);
126          }          }
127          else          else
128          {          {
# Line 112  Line 131 
131          $simI = $simI1;          $simI = $simI1;
132          $qI++;          $qI++;
133      }      }
134  #    print STDERR &Dumper(\@rerun);      #print STDERR &Dumper(rerun => \@rerun);
135    
136      if (@rerun > 0)      if (@rerun > 0)
137      {      {
138          my $tmpQ = "/tmp/tmpQ.$$";          my $tmpQ = "/tmp/tmpQ.$$";
139          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
140          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";
141                         `$FIG_Config::ext_bin/blastall -m 8 -i $tmpQ -d $dbF -FF -p blastp`;          #print STDERR "$cmd\n";
142            open(BL, "$cmd|") or die "ProtSims::blastP: pipe to blast failed with $!: $cmd\n";
143            my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
144            close(BL);
145          push(@sims,@blastout);          push(@sims,@blastout);
146          unlink $tmpQ;          unlink $tmpQ;
147      }      }
148    
149      my %lnQ   = map { $_->[0] => length($_->[2]) } @q;      my %lnQ   = map { $_->[0] => length($_->[2]) } @$q;
150      my %lnDB  = map { $_->[0] => length($_->[2]) } @db;      my %lnDB  = map { $_->[0] => length($_->[2]) } @$db;
151      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;
152    
153      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }
154      if ($dbF     eq "/tmp/db.$$")      { unlink $dbF }      if ($dbF     eq "/tmp/db.$$")      { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }
155      return @sims;      return sort { ($a->id1 cmp $b->id1) or ($a->id2 cmp $b->id2) or ($b->bsc <=> $a->bsc) } @sims;
156    }
157    
158    sub condense {
159        my($pieces) = @_;
160    
161        my $best_sc = $pieces->[0]->[10];
162        my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;
163        while (@sims > 1)
164        {
165            my $gap1 = $sims[1]->[6] - $sims[0]->[7];
166            my $gap2 = $sims[1]->[8] - $sims[0]->[9];
167            my $diff = abs($gap1 - $gap2);
168            if (($gap1 <= 300) && ($gap2 <= 300) && ($diff <= (0.1 * &min($gap1,$gap2))))
169            {
170                $sims[0]->[7] = $sims[1]->[7];
171                $sims[0]->[9] = $sims[1]->[9];
172            }
173            splice(@sims,1,1);
174        }
175        $sims[0]->[10] = $best_sc;
176        return $sims[0];
177    }
178    
179    sub min {
180        my($x,$y) = @_;
181        return ($x <= $y) ? $x : $y;
182  }  }
183    
184  1;  1;

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.8

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3