perl如何提取指定基因的fasta序列

小編給大家分享一下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序列”有了一定的了解，如果想了解更多相關知識，歡迎關注億速云行業資訊頻道，感謝各位的閱讀！