这里，介绍一个非常简便好用的序列提取的perl脚本，用法非常简单。
用法如下：
perl /share/work/huangls/piplines/01.script/get_fa_by_id.pl <id><fa><OUT>
例如：
perl /share/work/huangls/piplines/01.script/get_fa_by_id.pl id.txt input.fasta out.fa
其中 id.txt 为要提取的序列ID，input.fasta 为输入序列文件，out.fa 是输出提取的序列文件。
id.txt 格式如下：
TRINITY_DN116733_c6_g37
TRINITY_DN116733_c6_g70
TRINITY_DN95808_c0_g7
TRINITY_DN104586_c1_g2
TRINITY_DN108413_c2_g23
TRINITY_DN37223_c0_g1
TRINITY_DN107955_c0_g8
TRINITY_DN117047_c0_g2
TRINITY_DN78058_c0_g1
这里是脚本代码：
die "perl $0 <id><fa><OUT>" unless(@ARGV==3);
use Math::BigFloat;
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 /^#/;
#next unless />>/;
my@tmp=split(/\s+/);
$keep{$tmp[0]}=1;
}
close(IN);
my$i=0;
while ( my $seq = $in->next_seq() ) {
my($id,$sequence,$desc)=($seq->id,$seq->seq,$seq->desc);
if(exists $keep{$id}){
$out->write_seq($seq);
}
}
$in->close();
$out->close();
脚本使用了Bio::SeqIO模块来处理序列文件，简洁而高效，先使用哈希来存储要提取的序列ID，然后利用Bio::SeqIO遍历序列文件，判断每条序列是否是要提取的序列，是的话就输出。
希望对大家有帮助！
