perl如何提取合并基因家族的domain序列

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

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

使用方法：

命令如下：

perl   domain_hebing.fa.pl   hmmsearch.out.txt   Arabidopsis_thaliana.TAIR10.31.pep.all.fa   domain.fa   0.001

参数介绍：

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

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

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

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

脚本代码如下：

die "perl $0    " unless ( @ARGV == 4 );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 %fasta;while ( my $seq = $in->next_seq() ) {my ( $id, $sequence, $desc ) = ( $seq->id, $seq->seq, $seq->desc );$fasta{$id} = $seq;}$in->close();my %keep = ();open IN, "$ARGV[0]" or die "$!";while () {chomp;next if /^#/;my @a = split /\s+/;next if $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序列"这篇文章就分享到这里了，希望以上内容可以对大家有一定的帮助，让大家可以学到更多知识，如果觉得文章不错，可以把它分享出去让更多的人看到吧！