PERL ONLY. Use the BioPerl modules to write a script that takes a protein sequen

Question

PERL ONLY.

Use the BioPerl modules to write a script that takes a protein sequence accession number as input (on the command line), then fetches this sequence from GenBank, BLASTs it against the non-redundant database, and then retrieves all the hits that are above a certain e-value cutoff (also given as input on the command line). The retrieved sequences should be written to files, one sequence per file, in the current directory. If no hits satisfy the e-value cutoff, the script should note that instead of writing any files.

Explanation / Answer

use warnings;
use strict;
use Bio::Perl,
use Bio::Tools::Run::RemoteBlast;

my ($accession, $evalue) = @ARGV;

my $db = 'genepept';
my $sequence ($db , $accession);
my $factory = Bio::Tolls:Run::RemoteBlast->new ( -prog=> 'blastp' , -data => 'nr' , -expect => $evalue , -readmethod = > 'SearchIO' );

$factory->submit_blast ($sequence);
my $final = "";
while (my $rids == $facotry-> rids) {
   foreach my $rid (@rids) {
       $final = $facotry->get_blast($rid);
       if (ref( $final)) {
           my $output = $final->next();
           $facctory->remove( $rid);
           foreach my $return ($output->hits) {
               my $acc = $hit->accession;
               my $seq = get_sequence($db, $acc);
               my $file = $acc.".out";
               open (mmy $OUT, '>' , $file) or die ("Cannot open this $file");
               print $OUT $seq->sequence;
           }
       }
       elsif ( $final < 0 ) {
       $factory->remove( $rid);
       }
       else {
           sleep 10;
       }
   }
}
if ($final eq "") {
   print" No hits return with the e-value specified";
}

Related Questions

Navigate

• Browse (All)
• Browse P
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.