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序列”有了一定的了解，如果想了解更多相关知识，欢迎关注亿速云行业资讯频道，感谢各位的阅读！