[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.18, Mon Oct 4 21:31:34 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    use DB_File;
28    
29    use SeedAware;
30    
31    my $blat_cmd = SeedAware::executable_for("blat");
32    my $blastall_cmd = SeedAware::executable_for("blastall");
33    my $formatdb_cmd = SeedAware::executable_for("formatdb");
34    
35    #
36    # Need to call temps different names on each invocation -
37    # if a subprocess hung on windows, we get problems with
38    # "file in use" errors.
39    #
40    my $tmp_serial = 0;
41    
42  sub blastP {  sub blastP {
43      my($q,$db,$min_hits) = @_;      my($q,$db,$min_hits,$use_blast) = @_;
44    
45      my(@q,@db);      my(@q,@db);
46    
47        my $tmp_dir = SeedAware::location_of_tmp();
48    
49        my $tmp_suffix = $$ . "." . $tmp_serial++ . "." . time;
50    
51      my $qF;      my $qF;
52      if (ref $q)      if (ref $q)
53      {      {
54          $qF = "/tmp/query.$$";          $qF = "$tmp_dir/query.$tmp_suffix";
55          &gjoseqlib::print_alignment_as_fasta($qF,$q);          &gjoseqlib::print_alignment_as_fasta($qF,$q);
56      }      }
57      else      else
# Line 49  Line 68 
68          }          }
69      }      }
70    
71        my $db_lengths = {};
72    
73      my $dbF;      my $dbF;
74      if (ref $db)      if (ref $db)
75      {      {
76          $dbF = "/tmp/db.$$";          $dbF = "$tmp_dir/db.$tmp_suffix";
77          &gjoseqlib::print_alignment_as_fasta($dbF,$db);          &gjoseqlib::print_alignment_as_fasta($dbF,$db);
78          system "$FIG_Config::ext_bin/formatdb -i $dbF";          system($formatdb_cmd, '-i', $dbF);
79            $db_lengths->{$_->[0]} = length($_->[2]) for @$db;
80      }      }
81      else      else
82      {      {
# Line 63  Line 85 
85          {          {
86              if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))              if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))
87              {              {
88                  system "$FIG_Config::ext_bin/formatdb -i $dbF";                  system($formatdb_cmd, '-i', $dbF);;
89                }
90    
91                my $db_tie;
92                if (-f "$dbF.lengths")
93                {
94                    $db_tie = tie %$db_lengths, 'DB_File', "$dbF.lengths", O_RDONLY, 0, $DB_BTREE;
95                    $db_tie or warn "Cannot tie $dbF.lengths: $!";
96                }
97                if (!defined($db_tie))
98                {
99                    if (open(my $dbfh, "<", $dbF))
100                    {
101                        while (my($id, $def, $seq) = gjoseqlib::read_next_fasta_seq($dbfh))
102                        {
103                            $db_lengths->{$id} = length($seq);
104                        }
105                        close($dbfh);
106                    }
107                    else
108                    {
109                        warn "Cannot open $dbF: $!";
110                    }
111              }              }
             @db = &gjoseqlib::read_fasta($dbF);  
             $db = \@db;  
112          }          }
113          else          else
114          {          {
# Line 74  Line 116 
116          }          }
117      }      }
118    
119      my $tmpF = "/tmp/sim.out.$$";      my $tmpF = "$tmp_dir/sim.out.$tmp_suffix";
120      my $cmd = "$FIG_Config::ext_bin/blat $dbF $qF -prot -out=blast8 $tmpF > /dev/null";  
121      #print STDERR "$cmd\n";      if ($use_blast)
122      my $rc = system $cmd;      {
123            open(my $fh, ">", $tmpF);
124            close($fh);
125        }
126        else
127        {
128            my @cmd = ($blat_cmd, $dbF, $qF, "-prot", "-out=blast8", $tmpF);
129    
130            #
131            # When running under FCGI, the system_with_redirect fails due to
132            # the FCGI library redefining open.
133            #
134    
135            my $rc;
136            if (defined($ENV{FCGI_ROLE}))
137            {
138                $rc = system(@cmd);
139            }
140            else
141            {
142                $rc = system_with_redirect(\@cmd, { stdout => '/dev/null' } );
143            }
144    
145            print STDERR "Blat returns $rc: @cmd\n";
146      if ($rc != 0)      if ($rc != 0)
147      {      {
148          die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";              warn "Blat run failed with rc=$rc\n";
149                open(my $fh, ">", $tmpF);
150                close($fh);
151      }      }
152        }
153    
154    #     my $cmd = $use_blast ? "touch $tmpF" : "$blat_cmd $dbF $qF -prot -out=blast8 $tmpF > /dev/null";
155    #     #print STDERR "$cmd\n";
156    #     my $rc = system $cmd;
157    #     if ($rc != 0)
158    #     {
159    #       die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";
160    #     }
161    
162      my @sims1 = ();      my @sims1 = ();
163      open(BLAT,"<$tmpF") || die "could not open $tmpF";      open(BLAT,"<$tmpF") || die "could not open $tmpF";
# Line 135  Line 211 
211    
212      if (@rerun > 0)      if (@rerun > 0)
213      {      {
214          my $tmpQ = "/tmp/tmpQ.$$";          my $tmpQ = "$tmp_dir/tmpQ.$tmp_suffix";
215          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
216          my $cmd = "$FIG_Config::ext_bin/blastall -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";  
217    
218            #
219            # If we're under FCGI (and thus under the servers), and the loadavg is low, do a small parallel run.
220            #
221    
222            my @par = ();
223            if (defined($ENV{FCGI_ROLE}) && open(my $la, "/proc/loadavg"))
224            {
225                my $l = <$la>;
226                chomp $l;
227                my @vals = split(/\s+/, $l);
228                my @procs = split(/\//, $vals[3]);
229                if ($vals[0] < 4 && $procs[0] < 8)
230                {
231                    @par = ("-a", 4);
232                }
233            }
234    
235    
236            # my $cmd = "$blastall_cmd -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";
237            my @cmd = ($blastall_cmd,  '-m', 8, '-i', $tmpQ, '-d', $dbF, '-FF', '-p', 'blastp', '-e', '1e-5', @par);
238          #print STDERR "$cmd\n";          #print STDERR "$cmd\n";
239          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";
240          my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;  
241            #
242            # It'd be nice to do this but windows doesn't support it.
243            #
244            #open(BL, "-|", @cmd) or die "ProtSims::blastP: pipe to blast failed with $!: @cmd\n";
245            my $out_tmp = "$tmp_dir/blast_out.$tmp_suffix";
246            push(@cmd, "-o", $out_tmp);
247            print STDERR "@cmd\n";
248            my $rc = system(@cmd);
249            if ($rc != 0)
250            {
251                warn "Error rc=$rc running blast: @cmd\n";
252            }
253            else
254            {
255                if (open(BL, "<", $out_tmp))
256                {
257                    while (<BL>)
258                    {
259                        chomp;
260                        push @sims, [ split(/\s+/, $_) ];
261                    }
262                    #my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
263                    # push(@sims,@blastout);
264          close(BL);          close(BL);
265          push(@sims,@blastout);              }
266            }
267            unlink $out_tmp;
268          unlink $tmpQ;          unlink $tmpQ;
269      }      }
270    
271      my %lnQ   = map { $_->[0] => length($_->[2]) } @$q;      my %lnQ   = map { $_->[0] => length($_->[2]) } @$q;
     my %lnDB  = map { $_->[0] => length($_->[2]) } @$db;  
     @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;  
272    
273      if ($qF      eq "/tmp/query.$$")   { unlink $qF  }      @sims = map { push(@$_, $lnQ{$_->[0]}, $db_lengths->{$_->[1]}); bless($_,'Sim') } @sims;
274      if ($dbF     eq "/tmp/db.$$")      { unlink $dbF }  
275      return sort { ($a->id1 cmp $b->id1) or ($a->id2 cmp $b->id2) or ($b->bsc <=> $a->bsc) } @sims;      if ($qF      eq "$tmp_dir/query.$tmp_suffix")   { unlink $qF  }
276        if ($dbF     eq "$tmp_dir/db.$tmp_suffix")      { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }
277        return sort { ($a->id1 cmp $b->id1) or ($a->psc <=> $b->psc) or ($a->id2 cmp $b->id2) } @sims;
278  }  }
279    
280  sub condense {  sub condense {
281      my($pieces) = @_;      my($pieces) = @_;
282    
283      my $best_sc = $pieces->[0]->[9];      my $best_sc = $pieces->[0]->[10];
284      my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;      my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;
285      while ((@sims > 1) &&      while (@sims > 1)
286             ($sims[1]->[6] <= ($sims[0]->[7] + 100)) &&      {
287             ($sims[1]->[8] <= ($sims[0]->[9] + 100)) &&          my $gap1 = $sims[1]->[6] - $sims[0]->[7];
288             (abs(($sims[1]->[6] - $sims[0]->[7]) - ($sims[1]->[8] - $sims[0]->[9])) < 20))          my $gap2 = $sims[1]->[8] - $sims[0]->[9];
289            my $diff = abs($gap1 - $gap2);
290            if (($gap1 <= 300) && ($gap2 <= 300) && ($diff <= (0.1 * &min($gap1,$gap2))))
291      {      {
292          $sims[0]->[7] = $sims[1]->[7];          $sims[0]->[7] = $sims[1]->[7];
293          $sims[0]->[9] = $sims[1]->[9];          $sims[0]->[9] = $sims[1]->[9];
294            }
295          splice(@sims,1,1);          splice(@sims,1,1);
296      }      }
297        $sims[0]->[10] = $best_sc;
298      return $sims[0];      return $sims[0];
299  }  }
300    
301    sub min {
302        my($x,$y) = @_;
303        return ($x <= $y) ? $x : $y;
304    }
305    
306  1;  1;

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3