KaKs计算之BUG— 序列长度不相等！

在遗传学中，Ka/Ks表示的是两个蛋白编码基因的非同义替换率（Ka）和同义替换率（Ks）之间的比例，这个比例可以判断是否有选择压力作用于这个蛋白质编码基因。

要计算两个基因之间的Ka/Ks值，一般步骤：1.利用clustalw进行多序列比对（结果为aln格式），2.将比对结果转换成Ka/Ks_Calculator2.0需要的格式axt文件（转换用到的软件为Ka/Ks_Calculator2.0提供的AXTConvertor工具），之后就可以计算Ka/Ks了。具体方法可参考如何使用KaKs_Calculator计算Ka/Ks。

格式转换

在使用axt文件计算ka/ks时，有时候计算会报错：如下图所示：

1.jpg

软件提示两个输入fasta序列不相等，但是，我通过clustalw多序列比对的之后结果应该是相等的，怎么转换过后就不相等呢？

其实AXTConvertor工具在转换过程中会出现错误，也就是软件bug，导致转换为axt格式后序列长度不相等。那该怎么办呢？这里我们提供一个Perl脚本，可以将aln文件转换为axt文件。

使用方法：
perl   aln2axt.pl   a.aln   b.axt
脚本代码：
#email: huangls@biomics.com.cn
die "perl $0 <fa> <OUT>" unless ( @ARGV == 2 );
use Bio::AlignIO;
use Bio::SimpleAlign;
$in = Bio::AlignIO->new(
-file   => "$ARGV[0]",
-format => 'clustalw'
);
open OUT, ">$ARGV[1]" or die "$!";
while(my $aln = $in->next_aln() ){
$seq1 = $aln->get_seq_by_pos(1);
($id1, $sequence1) = ( $seq1->id, $seq1->seq);
$seq2 = $aln->get_seq_by_pos(2);
($id2, $sequence2) = ( $seq2->id, $seq2->seq);
print OUT "$id1&$id2\n$sequence1\n$sequence2\n\n";
}
$in->close();
close(OUT);
这样最后得到的序列长度就是相等的了，如果提示依然不相等，则可能是：
1. 序列比对总长度是不是3的整倍数，如果不是3的倍数可能计算不出结果；
2. 看看两条序列是否都含有起始密码子ATG，如果没有可能也会报错；
遇到上述两条原因解析：可能基因组注释的不好，有些基因不是很完整，建议筛选条件严格一些，比如相似性和比对率都提高一些。

序列提取

计算ka/ks时，如果有多组串联重复基因，则需要分别提取每组串联重复基因序列并进行序列比对。如果串联重复基因数量太多，手动提取就会比较麻烦，这里我们有一个快速提取的脚本。
使用方法：
1.准备串联重复基因对文件。每一行是一对串联重复基因，中间用Tab键分割。如下图：

2.jpg

2.准备所有基因cds序列。

3.jpg

3.运行命令
perl get_fa_by_id_kaks.pl WRKY.tandem Sspon.v20180123.cds.fa ./
WRKY.tandem：串联重复基因对文件；
Sspon.v20180123.cds.fa：所有基因cds序列文件；
脚本代码：
#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();
好啦，有需要的小伙伴拿去用吧！
分享