KaKs_Calculator使用中需要用到axt格式的比對(duì)文件,在KaKs_Calculator的文件夾中自帶一個(gè)生成axt格式文件的腳本,?AXTConvertor ,可以批量把比對(duì)結(jié)果文件轉(zhuǎn)化為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. 然而在計(jì)算ks/ks值的時(shí)候發(fā)現(xiàn)了一些問(wèn)題。
首先在計(jì)算時(shí)出現(xiàn)了極高的ka,ks值,而且p值忽高忽低,有很多的s-sites和a-sites。但這是不應(yīng)該的,按照比對(duì)時(shí)的經(jīng)驗(yàn)看大部分堿基其實(shí)都一樣,align超過(guò)90%。
觀察aln文件也沒(méi)有什么問(wèn)題,文件都是對(duì)齊的。
打開(kāi)axt文件一看,果然有很大的問(wèn)題

序列后半段均有不同程度的錯(cuò)位,而非之前比對(duì)的結(jié)果,似乎部分位置的gaps數(shù)量減少
2. 進(jìn)一步的改正
回到程序發(fā)布頁(yè),發(fā)現(xiàn)了paraAT這個(gè)軟件,可以把序列轉(zhuǎn)換成KaKs_Calculator使用的格式,于是就研究了一下
運(yùn)行環(huán)境:ubuntu/deepin
需要安裝的軟件:clustalw2 | t_coffee | mafft | muscle 任選其一,我安裝了mafft作為比對(duì)軟件。
生成axt需要的文件:
1? ? 多序列比對(duì)的fasta序列,兩端對(duì)齊
2? ? 多序列的蛋白質(zhì)序列
3? ? 序列名兩兩比對(duì)的名錄
(這里我刪去了我在研究的物種)

都放在paraAT文件夾下
運(yùn)行前的準(zhǔn)備 把paraAT放進(jìn)環(huán)境變量
export PATH=$PATH:/home/username/Desktop/ParaAT2.0
此時(shí),在終端的paraAT路徑下輸入如下命令(文件名是我自己的文件名)
perl ParaAT.pl -h hom.TXT -n mafft_kaks.fasta -a ksks_translation.fasta -p proc -o output -f axt
在output文件夾下即為axt兩兩比對(duì)的所有文件n*(n-1)/2個(gè)比對(duì)文件
cat *.axt > kaks.all.axt
kaks.all.axt即可作為KaKs_Calculator的輸入文件,打開(kāi)看一下,應(yīng)該沒(méi)有錯(cuò)位。

計(jì)算一下,數(shù)值沒(méi)有太大異常,p值也都小于0.05,應(yīng)該是真實(shí)值。
在excel中以熱圖的形式展示出來(lái)
當(dāng)然用這種方法忽視了親緣遠(yuǎn)近對(duì)ka/ks值的影響,更加合理的方式應(yīng)是使用paml中dnds計(jì)算的包,加入系統(tǒng)發(fā)育樹(shù)來(lái)獲得考慮了多次顛換倒換后的值。

具體來(lái)說(shuō),使用KaKs_Calculator2.0 做ka/ks的大致過(guò)程分以下幾步
1選擇每個(gè)物種相同數(shù)目的cds編碼區(qū),每個(gè)物種的cds應(yīng)首先去除終止子,然后按照相同的基因順序聯(lián)合拼接 拼接過(guò)程要注意數(shù)目始終為3的倍數(shù)。
2 導(dǎo)出每個(gè)物種cds拼接序列翻譯的對(duì)應(yīng)氨基酸序列
3 制作一一對(duì)應(yīng)的兩兩比較序列名稱(chēng)表
4 使用KaKs_Calculator 制作者的另一款軟件 paraAT軟件,通過(guò)堿基序列、氨基酸序列、序列名稱(chēng)表獲得對(duì)應(yīng)的axt文件用于KaKs_Calculator對(duì)kaks的計(jì)算。或者使用phylosuite軟件導(dǎo)出對(duì)應(yīng)的cds和氨基酸序列。
生成axt文件后要手動(dòng)查看序列是否兩兩對(duì)應(yīng),或者總序列長(zhǎng)度是否是3的整數(shù)。不可以直接使用多序列比對(duì)文件和conaxt腳本生成axt文件,這會(huì)使得多序列比對(duì)結(jié)果中添加的多余的空格影響到密碼子順序。如果paraAT結(jié)果不好,可以手動(dòng)兩兩序列比對(duì)并手動(dòng)排列為axt文件。
5 KaKs_Calculator 計(jì)算ka/ks值,同義突變和非同義突變的位點(diǎn)應(yīng)該較少,如果較多則說(shuō)明axt文件有錯(cuò)誤 需要返回查看。同時(shí),計(jì)算得到ka/ks的值時(shí),相應(yīng)的p值應(yīng)該遠(yuǎn)小于0.05,才具有統(tǒng)計(jì)學(xué)的意義。只有在序列非常短的情況下才會(huì)出現(xiàn)較高的ka和ka/ks值。