针对blast特定案例下批量处理程序

问题为：首先数据源是已经处理掉引物并且拼接好的单条碱基。需要通过BLAST比对得到其Features是否属于特定蛋白质？

BlAST比对采用perl连接远程BLAST从而进行批量比对。

首先对数据源处理将所有数据放到一个FASTA文件中。

#!/usr/bin/perl
use strict;
use Bio::SeqIO;
use Getopt::Long;
my $dirname="D:/sutest";#输入目录
my $targetFile="E:/sutest/";#输出目录

my $result='';
opendir ( DIR,$dirname ) || die "Error in opening dir $dirnamen";
open(OUTFILE,">>".$targetFile."result.fasta");
while( (my $filename = readdir(DIR))){

if(!($filename=~/^./)&&$filename=~/^[a-zA-Z]/){

            print OUTFILE (">$filename"."n");
open(IN,$dirname."/".$filename)||die "Can not open!";
while(my $seq= <IN>){
    print OUTFILE ($seq."n");


}



close(IN);
}
}
close(OUTFILE);

然后针对这些数据采用远程blast比对，其代码参考网络上的bioperl的方法

use strict;
use warnings;
use Bio::Tools::Run::RemoteBlast;#调用bioperl的RemoteBlast模块
my $prog='blastn';#选择blast的程序，核酸用blastn，还有另外几个blastp，blastx等，根据比对序列类型选择
my $db='nt';#选择比对的数据库，蛋白质的可以选swissport，除了nt，还可以用nr
my $e_value='1e-10';#设定e值
my $factory = Bio::Tools::Run::RemoteBlast->new('-prog' => $prog,'-data' => $db,'-expect' => $e_value,'-readmethod' => 'SearchIO' );#这是perl面向对象编程的用法，不做过多解释，一般写法是固定的
my $v = 1;
my $str = Bio::SeqIO->new(-file=>'E:/result.fasta',-format => 'fasta' );#读入序列，多条序列放在Query.fasta文件
#下面代码看起来有点复杂，不过基本原理是先将Query.fasta文件中的序列提交到NCBI（多条序列通过循环方式一次提交），再接收比对的结果（多个结果接收也是通过循环实现的），最后按序列的名字把结果写入到用户当前的文件夹下，比对结果文件均以.out结尾，可以用记事本打开。
while (my $input = $str->next_seq())
{
my $r = $factory->submit_blast($input);
print STDERR "waiting..." if( $v > 0 );
    while ( my @rids = $factory->each_rid )
{
    foreach my $rid ( @rids )
  {
    my $rc = $factory->retrieve_blast($rid);
     if( !ref($rc) )
   {
              if( $rc < 0 )
      {
                $factory->remove_rid($rid);
       }
              print STDERR "." if ( $v > 0 );
              sleep 5;
            } else
   {
              my $result = $rc->next_result();
              #save the output
              my $filename = $result->query_name().".out";
              $factory->save_output($filename);
              $factory->remove_rid($rid);
              print "nQuery Name: ",$result->query_name(),"n";
              while ( my $hit = $result->next_hit )
   {
                next unless ( $v > 0);
                print "thit name is ",$hit->name,"n";
                while( my $hsp = $hit->next_hsp )
    {
                  print "ttscore is ",$hsp->score,"n";
                }
              }
            }
          }
        }
}

#!/usr/bin/perl -w
use strict;
my $dirname="E:/sutest/";#输入目录
my $targetFile="D:/sutest/";#输出目录
my $flag=0;
opendir ( DIR,">>".$targetFile."result.fasta");
while(my $filename=readdir(DIR)){
 if(!($filename=~/^./) && $filename=~/^[a-zA-Z]/){

open(IN,$dirname.$filename) || die "Can not open!";

while(my $seq= <IN>){
    if($seq=~/.CheA/g){#CHEA即为目标蛋白质

      print OUTFILE (">$filename chemotaxis protein CheA"."n");
      $flag=1;
      last;
    }
}
if($flag==0){
  print OUTFILE (">$filename none"."n");

}
 $flag=0;
 close(IN);
}
}

closedir(DIR);
close(OUTFILE);