遗传学选择压力分析湿实验

KaKs_Calculator2.0 做ka/ks过程的一点记录

2019-06-26  本文已影响0人  徒唤奈何_c5f0

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值。

上一篇下一篇

猜你喜欢

热点阅读