Perl脚本练习
bioperl读入写出fasta
要求
根据序列ID,从fasta文件中提取目标序列并输出
数据
序列ID
fasta文件
思路
- 以序列ID为键,构建哈希
- 用bioperl读入fasta,获得序列id
- 如果id存在于哈希中,输出序列
代码
die "perl $0 <id> <fa> <OUT>" unless(@ARGV==3);
#$0程序名
use Bio::SeqIO;
use Bio::Seq;
$in = Bio::SeqIO->new(-file=>"$ARGV[1]", -format =>'Fasta');
$out = Bio::SeqIO->new(-file=>"$ARGV[2]", -format=>'Fasta');
my %keep=();
open (IN, "$ARGV[0]") or die "$!";
while(<IN>){
chomp;
next if ($_=~/^#/);
$keep{$_}=1;
}
close IN;
while(my$seqobj = $in->next_seq()){
my($id, $desc) = ($seqobj->id, $seqobj->desc);
if(exists $keep{$id}){
$out->wirte_seq($seqobj);
}
}
$in->close();
$out->close();