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序列”這篇文章就分享到這里了，希望以上內容可以對大家有一定的幫助，讓大家可以學到更多知識，如果覺得文章不錯，可以把它分享出去讓更多的人看到吧！