小编给大家分享一下perl如何提取指定基因的fasta序列，希望大家阅读完这篇文章之后都有所收获，下面让我们一起去探讨吧！
简便好用的序列提取的perl脚本
这里，介绍一个非常简便好用的序列提取的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遍历序列文件，判断每条序列是否是要提取的序列，是的话就输出。
看完了这篇文章，相信你对“perl如何提取指定基因的fasta序列”有了一定的了解，如果想了解更多相关知识，欢迎关注云搜网行业资讯频道，感谢各位的阅读！
