[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.13, Mon Aug 30 16:09:00 2010 UTC revision 1.16, Tue Sep 14 21:34:35 2010 UTC
# Line 24  Line 24 
24  use Data::Dumper;  use Data::Dumper;
25  use Carp;  use Carp;
26  use gjoseqlib;  use gjoseqlib;
27    use DB_File;
28    
29  use SeedAware;  use SeedAware;
30    
# Line 58  Line 59 
59          }          }
60      }      }
61    
62        my $db_lengths = {};
63    
64      my $dbF;      my $dbF;
65      if (ref $db)      if (ref $db)
66      {      {
67          $dbF = "$tmp_dir/db.$$";          $dbF = "$tmp_dir/db.$$";
68          &gjoseqlib::print_alignment_as_fasta($dbF,$db);          &gjoseqlib::print_alignment_as_fasta($dbF,$db);
69          system($formatdb_cmd, '-i', $dbF);          system($formatdb_cmd, '-i', $dbF);
70            $db_lengths->{$_->[0]} = length($_->[2]) for @$db;
71      }      }
72      else      else
73      {      {
# Line 74  Line 78 
78              {              {
79                  system($formatdb_cmd, '-i', $dbF);;                  system($formatdb_cmd, '-i', $dbF);;
80              }              }
81              @db = &gjoseqlib::read_fasta($dbF);  
82              $db = \@db;              my $db_tie;
83                if (-f "$dbF.lengths")
84                {
85                    $db_tie = tie %$db_lengths, 'DB_File', "$dbF.lengths", O_RDONLY, 0, $DB_BTREE;
86                    $db_tie or warn "Cannot tie $dbF.lengths: $!";
87                }
88                if (!defined($db_tie))
89                {
90                    if (open(my $dbfh, "<", $dbF))
91                    {
92                        while (my($id, $def, $seq) = gjoseqlib::read_next_fasta_seq($dbfh))
93                        {
94                            $db_lengths->{$id} = length($seq);
95                        }
96                        close($dbfh);
97                    }
98                    else
99                    {
100                        warn "Cannot open $dbF: $!";
101                    }
102                }
103          }          }
104          else          else
105          {          {
# Line 94  Line 118 
118      {      {
119          my @cmd = ($blat_cmd, $dbF, $qF, "-prot", "-out=blast8", $tmpF);          my @cmd = ($blat_cmd, $dbF, $qF, "-prot", "-out=blast8", $tmpF);
120    
121          my $rc = system_with_redirect(\@cmd, { stdout => '/dev/null' } );          #
122            # When running under FCGI, the system_with_redirect fails due to
123            # the FCGI library redefining open.
124            #
125    
126            my $rc;
127            if (defined($ENV{FCGI_ROLE}))
128            {
129                $rc = system(@cmd);
130            }
131            else
132            {
133                $rc = system_with_redirect(\@cmd, { stdout => '/dev/null' } );
134            }
135    
136          print STDERR "Blat returns $rc: @cmd\n";          print STDERR "Blat returns $rc: @cmd\n";
137          if ($rc != 0)          if ($rc != 0)
138          {          {
# Line 166  Line 204 
204      {      {
205          my $tmpQ = "$tmp_dir/tmpQ.$$";          my $tmpQ = "$tmp_dir/tmpQ.$$";
206          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);          &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
207    
208    
209            #
210            # If we're under FCGI (and thus under the servers), and the loadavg is low, do a small parallel run.
211            #
212    
213            my @par = ();
214            if (defined($ENV{FCGI_ROLE}) && open(my $la, "/proc/loadavg"))
215            {
216                my $l = <$la>;
217                chomp $l;
218                my @vals = split(/\s+/, $l);
219                my @procs = split(/\//, $vals[3]);
220                if ($vals[0] < 4 && $procs[0] < 8)
221                {
222                    @par = ("-a", 4);
223                }
224            }
225    
226    
227          # my $cmd = "$blastall_cmd -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";
228          my @cmd = ($blastall_cmd,  '-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', @par);
         print STDERR "@cmd\n";  
229          #print STDERR "$cmd\n";          #print STDERR "$cmd\n";
230          #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";
231    
232          open(BL, "-|", @cmd) or die "ProtSims::blastP: pipe to blast failed with $!: @cmd\n";          #
233            # It'd be nice to do this but windows doesn't support it.
234          my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;          #
235            #open(BL, "-|", @cmd) or die "ProtSims::blastP: pipe to blast failed with $!: @cmd\n";
236            my $out_tmp = "$tmp_dir/blast_out.$$";
237            push(@cmd, "-o", $out_tmp);
238            print STDERR "@cmd\n";
239            my $rc = system(@cmd);
240            if ($rc != 0)
241            {
242                warn "Error rc=$rc running blast: @cmd\n";
243            }
244            else
245            {
246                if (open(BL, "<", $out_tmp))
247                {
248                    while (<BL>)
249                    {
250                        chomp;
251                        push @sims, [ split(/\s+/, $_) ];
252                    }
253                    #my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
254                    # push(@sims,@blastout);
255          close(BL);          close(BL);
256          push(@sims,@blastout);              }
257            }
258            unlink $out_tmp;
259          unlink $tmpQ;          unlink $tmpQ;
260      }      }
261    
262      my %lnQ   = map { $_->[0] => length($_->[2]) } @$q;      my %lnQ   = map { $_->[0] => length($_->[2]) } @$q;
263      my %lnDB  = map { $_->[0] => length($_->[2]) } @$db;  
264      @sims = map { push(@$_,$lnQ{$_->[0]},$lnDB{$_->[1]}); bless($_,'Sim') } @sims;      @sims = map { push(@$_, $lnQ{$_->[0]}, $db_lengths->{$_->[1]}); bless($_,'Sim') } @sims;
265    
266      if ($qF      eq "$tmp_dir/query.$$")   { unlink $qF  }      if ($qF      eq "$tmp_dir/query.$$")   { unlink $qF  }
267      if ($dbF     eq "$tmp_dir/db.$$")      { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }      if ($dbF     eq "$tmp_dir/db.$$")      { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }

Legend:
Removed from v.1.13  
changed lines
  Added in v.1.16

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3