如何计算Kaks时批量提取多对基因的序列

小编给大家分享一下如何计算Kaks时批量提取多对基因的序列，希望大家阅读完这篇文章之后都有所收获，下面让我们一起去探讨吧！

计算Kaks时批量提取多对基因的序列

基因家族分析中，在计算kaks时需要将每组串联重复基因序列提取出，然后分别进行多序列比对。

我们给大家提供的脚本只能一对一对的提取序列，那数量少还可以，可如果串联重复基因太多就有些麻烦了。这里提供一个批量提取串联重复基因序列的脚本，可以将每组串联重复基因序列提取到单独的文件中。

首先我们需要准备3个文件：

1.串联重复基因对文件。每一行是一对串联重复基因，中间用Tab键分割。

2.所有基因cds序列。

3.一个空的目录。

脚本运行命令如下：

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

WRKY.tandem：串联重复基因对文件；

../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时批量提取多对基因的序列”有了一定的了解，如果想了解更多相关知识，欢迎关注亿速云行业资讯频道，感谢各位的阅读！