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