KaKs_Calculator2.0 做ka/ks过程的一点记录
KaKs_Calculator使用中需要用到axt格式的比对文件,在KaKs_Calculator的文件夹中自带一个生成axt格式文件的脚本, AXTConvertor ,可以批量把比对结果文件转化为axt文件。
x@x-PC:~/Desktop/KaKs_Calculator2.0/bin/Linux/home/x/Desktop/KaKs_Calculator2.0/src/AXTConvertor /home/x/Desktop/KaKs_Calculator2.0/bin/Linux/20_spe_cds.aln 20_spe_cds.axt
1. 然而在计算ks/ks值的时候发现了一些问题。
首先在计算时出现了极高的ka,ks值,而且p值忽高忽低,有很多的s-sites和a-sites。但这是不应该的,按照比对时的经验看大部分碱基其实都一样,align超过90%。
观察aln文件也没有什么问题,文件都是对齐的。
打开axt文件一看,果然有很大的问题
序列后半段均有不同程度的错位,而非之前比对的结果,似乎部分位置的gaps数量减少
2. 进一步的改正
回到程序发布页,发现了paraAT这个软件,可以把序列转换成KaKs_Calculator使用的格式,于是就研究了一下
运行环境:ubuntu/deepin
需要安装的软件:clustalw2 | t_coffee | mafft | muscle 任选其一,我安装了mafft作为比对软件。
生成axt需要的文件:
1 多序列比对的fasta序列,两端对齐
2 多序列的蛋白质序列
3 序列名两两比对的名录
(这里我删去了我在研究的物种)
都放在paraAT文件夹下
运行前的准备 把paraAT放进环境变量
export PATH=$PATH:/home/username/Desktop/ParaAT2.0
此时,在终端的paraAT路径下输入如下命令(文件名是我自己的文件名)
perl ParaAT.pl -h hom.TXT -n mafft_kaks.fasta -a ksks_translation.fasta -p proc -o output -f axt
在output文件夹下即为axt两两比对的所有文件n*(n-1)/2个比对文件
cat *.axt > kaks.all.axt
kaks.all.axt即可作为KaKs_Calculator的输入文件,打开看一下,应该没有错位。
计算一下,数值没有太大异常,p值也都小于0.05,应该是真实值。
在excel中以热图的形式展示出来
当然用这种方法忽视了亲缘远近对ka/ks值的影响,更加合理的方式应是使用paml中dnds计算的包,加入系统发育树来获得考虑了多次颠换倒换后的值。
具体来说,使用KaKs_Calculator2.0 做ka/ks的大致过程分以下几步
1选择每个物种相同数目的cds编码区,每个物种的cds应首先去除终止子,然后按照相同的基因顺序联合拼接 拼接过程要注意数目始终为3的倍数。
2 导出每个物种cds拼接序列翻译的对应氨基酸序列
3 制作一一对应的两两比较序列名称表
4 使用KaKs_Calculator 制作者的另一款软件 paraAT软件,通过碱基序列、氨基酸序列、序列名称表获得对应的axt文件用于KaKs_Calculator对kaks的计算。或者使用phylosuite软件导出对应的cds和氨基酸序列。
生成axt文件后要手动查看序列是否两两对应,或者总序列长度是否是3的整数。不可以直接使用多序列比对文件和conaxt脚本生成axt文件,这会使得多序列比对结果中添加的多余的空格影响到密码子顺序。如果paraAT结果不好,可以手动两两序列比对并手动排列为axt文件。
5 KaKs_Calculator 计算ka/ks值,同义突变和非同义突变的位点应该较少,如果较多则说明axt文件有错误 需要返回查看。同时,计算得到ka/ks的值时,相应的p值应该远小于0.05,才具有统计学的意义。只有在序列非常短的情况下才会出现较高的ka和ka/ks值。