这篇文章给大家分享的是有关perl如何提取合并基因家族的domain序列的内容。小编觉得挺实用的，因此分享给大家做个参考，一起跟随小编过来看看吧。

在做基因家族分析中，使用hmmsearch搜索到结构域后，如果基因有多个domain，并且需要将所有domain提取连接在一起，可以使用下面的脚本完成。

使用方法：

命令如下：

perldomain_hebing.fa.plhmmsearch.out.txtArabidopsis_thaliana.TAIR10.31.pep.all.fadomain.fa0.001

参数介绍：

hmmsearch.out.txt ：hmmsearch搜索的结果文件。

Arabidopsis_thaliana.TAIR10.31.pep.all.fa ：所有基因蛋白质序列。

domain.fa ：输出合并后的domain序列文件。

0.001 ：设置过滤的e-value值。

脚本代码如下：

die"perl$0<hmmoutfile><fa><OUT><E-value>"unless(@ARGV==4);useMath::BigFloat;useBio::SeqIO;useBio::Seq;$in=Bio::SeqIO->new(-file=>"$ARGV[1]",-format=>'Fasta');$out=Bio::SeqIO->new(-file=>">$ARGV[2]",-format=>'Fasta');my%fasta;while(my$seq=$in->next_seq()){my($id,$sequence,$desc)=($seq->id,$seq->seq,$seq->desc);$fasta{$id}=$seq;}$in->close();my%keep=();openIN,"$ARGV[0]"ordie"$!";while(<IN>){chomp;nextif/^#/;my@a=split/\s+/;nextif$a[6]>$ARGV[3];my$subseq=$fasta{$a[0]}->subseq($a[17],$a[18]);if(exists$keep{$a[0]}){$keep{$a[0]}.=$subseq;}else{$keep{$a[0]}=$subseq;}}close(IN);while(my($key,$value)=each%keep){my$newseqobj=Bio::Seq->new(-seq=>$value,-id=>$key,);$out->write_seq($newseqobj);}$out->close();

感谢各位的阅读！关于“perl如何提取合并基因家族的domain序列”这篇文章就分享到这里了，希望以上内容可以对大家有一定的帮助，让大家可以学到更多知识，如果觉得文章不错，可以把它分享出去让更多的人看到吧！