根据id提取fasta序列

Perl脚本练习

bioperl读入写出fasta

要求

根据序列ID，从fasta文件中提取目标序列并输出

数据

序列ID
fasta文件

思路

1. 以序列ID为键，构建哈希
2. 用bioperl读入fasta，获得序列id
3. 如果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();