如何計算Kaks時批量提取多對基因的序列

小編給大家分享一下如何計算Kaks時批量提取多對基因的序列，希望大家閱讀完這篇文章之后都有所收獲，下面讓我們一起去探討吧！

計算Kaks時批量提取多對基因的序列

基因家族分析中，在計算kaks時需要將每組串聯(lián)重復(fù)基因序列提取出，然后分別進(jìn)行多序列比對。

我們給大家提供的腳本只能一對一對的提取序列，那數(shù)量少還可以，可如果串聯(lián)重復(fù)基因太多就有些麻煩了。這里提供一個批量提取串聯(lián)重復(fù)基因序列的腳本，可以將每組串聯(lián)重復(fù)基因序列提取到單獨的文件中。

首先我們需要準(zhǔn)備3個文件：

1.串聯(lián)重復(fù)基因?qū)ξ募?。每一行是一對串?lián)重復(fù)基因，中間用Tab鍵分割。

2.所有基因cds序列。

3.一個空的目錄。

腳本運行命令如下：

perl  get_fa_by_id_kaks.pl   WRKY.tandem   ../data/Sspon.v20180123.cds.fa   test/

WRKY.tandem：串聯(lián)重復(fù)基因?qū)ξ募?/p>

../data/Sspon.v20180123.cds.fa：所有基因cds序列文件；

test/：空目錄。

get_fa_by_id_kaks.pl就是腳本了，其代碼如下：

#email: huangls@biomics.com.cn
die "perl $0 <idlist> <fa> <OUT>" unless ( @ARGV == 3 );
use Math::BigFloat;
use Bio::SeqIO;
use Bio::Seq;
$in = Bio::SeqIO->new(
        -file   => "$ARGV[1]",
        -format => 'Fasta'
);
my %gene;
while ( my $seq = $in->next_seq() ) {
        my ( $id, $sequence, $desc ) = ( $seq->id, $seq->seq, $seq->desc );
        $gene{$id} = $seq;
}
open IN, "$ARGV[0]" or die "$!";
my $n = 1;
while (<IN>) {
        chomp;
        next if /^#/;
        my @a = split /\s+/;
        my $out = Bio::SeqIO->new(
                        -file   => ">$ARGV[2]/dup_gene_paired$n.fa",
                        -format => 'Fasta'
        );
        if ( exists $gene{$a[0]} ) {
                my ( $id, $sequence, $desc ) = ( $gene{$a[0]}->id, $gene{$a[0]}->seq, $gene{$a[0]}->desc );
                my $newSeqobj = Bio::Seq->new(-seq => $sequence,
                -id =>$id,
                );
                $out->write_seq($newSeqobj);
        }
        if ( exists $gene{$a[1]} ) {
                my ( $id, $sequence, $desc ) = ( $gene{$a[1]}->id, $gene{$a[1]}->seq, $gene{$a[1]}->desc );
                my $newSeqobj = Bio::Seq->new(-seq => $sequence,
                -desc => $desc,
                -id =>$id,
                );
                $out->write_seq($newSeqobj);
        }
        $n++;
        $out->close();
}
close(IN);
$in->close();

看完了這篇文章，相信你對“如何計算Kaks時批量提取多對基因的序列”有了一定的了解，如果想了解更多相關(guān)知識，歡迎關(guān)注億速云行業(yè)資訊頻道，感謝各位的閱讀！