这篇文章给大家分享的是有关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 <hmmoutfile> <fa> <OUT> <E-value>" 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 (<IN>) {
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序列”这篇文章就分享到这里了，希望以上内容可以对大家有一定的帮助，让大家可以学到更多知识，如果觉得文章不错，可以把它分享出去让更多的人看到吧！
