[Bioperl-l] issue with closing tempfile and bl2seq
BHurwitz at twt.com
BHurwitz at twt.com
Wed Mar 5 09:23:38 EST 2003
Hi,
I am having some trouble using bl2seq within Bioperl. I have a fairly
large file of sequences that need to be pairwise compared (40K).
Everything runs fine for the first ~500 sequences, but then my program
errors out with the following message:
------------- EXCEPTION -------------
MSG: Could not open tempfile /tmp/bhurwitz.16882.3051: Too many open files
STACK Bio::Root::IO::tempfile
/usr/lib/perl5/site_perl/5.6.0/Bio/Root/IO.pm:500
STACK Bio::Tools::Run::StandAloneBlast::new
/usr/lib/perl5/site_perl/5.6.0/Bio/Tools/Run/StandAloneBlast.pm:352
STACK main::pairwise ./footprint_checker.pl:117
STACK toplevel ./footprint_checker.pl:82
--------------------------------------
I haven't had a chance to update this machine yet to Bioperl 1.2 yet ( it
is still using Bioperl 1.0 ) so maybe this is why? Has anyone ever seen
this error? Is there a way for me to push all of my sequences through and
close these tmp files as I go? I have tried to pass the sequences into my
subroutine as both a file of 40K and as files of only 1 sequence, but it
seems to error out either way. I am sure there must be an easy fix??
Thanks,
Bonnie
Here is my subroutine that calls bl2seq
sub pairwise {
my $infile1 = shift @_;
my $infile2 = shift @_;
my $set = shift @_;
my $query1_in = Bio::SeqIO->new ( -file => $infile1,
-format => 'fasta' );
while (my $q1 = $query1_in->next_seq() ) {
push(@query1_seqs, $q1);
}
my $query2_in = Bio::SeqIO->new ( -file => $infile2,
-format => 'fasta' );
while (my $q2 = $query2_in->next_seq() ) {
push(@query2_seqs, $q2);
}
my $seq_count = @query1_seqs;
print "seqcount: $seq_count\n";
my $i;
for ($i = 0 ; $i < $seq_count; $i++) {
my $query1 = shift @query1_seqs;
my $query2 = shift @query2_seqs;
my $factory = Bio::Tools::Run::StandAloneBlast->new( );
$factory->e($expectvalue);
$factory->W($wordsize);
$factory->F($repeat_mask);
$factory->r($match_reward);
$factory->q($mismatch_penalty);
$factory->p($blprogram);
$factory->g($gapped);
my $report = $factory->bl2seq($query1, $query2);
my $sbjct_name = $report->sbjctName;
$sbjct_name =~ s/>\s*//;
while(my $hsp = $report->next_feature) {
my $pvalue = $hsp->P;
my $match = $hsp->match;
my $length = $hsp->length;
my $q_len;
if ($set =~ /syn1/) {
$q_len = $syn1_len{$sbjct_name};
}
else {
$q_len = $syn2_len{$sbjct_name};
}
my $frac_iden = $match / $q_len;
my $s_start = $hsp->hit->start( ); # hit start
my $s_end = $hsp->hit->end( ); # hit end
print OUT "$sbjct_name\t$pvalue\t$s_start\t$s_end\t$frac_iden\n";
}
}
} #end of sub pairwise
More information about the Bioperl-l
mailing list