文章來源：http://www.cnblogs.com/chenwenyan/
1 文件準備
首先準備兩個文件，一個是基因的cds序列，一個是蛋白質(zhì)序列。
cds序列和蛋白質(zhì)可以在ensembl網(wǎng)站找到：http://ftp.ensembl.org/pub/current_fasta/
這兩個文件的示例如下：
cds序列文件cds.fa
>gene1
ATGGAGGTTGCAATGGTGAGTGCGGAGAGCTCAGGATGCAACAGTCACATGCCTTACGGT
TATGCTGCCCAGGCCCGGGCCCGGGAGCGGGAGAGGCTTGCTCACTCCAGGGCAGCTGCG
GCAGCTGCCGTTGCAGCGGCCACAGCTGCCGTCGAAGGAAGTGGGGGTTCTGGTGGGGGG
>gene2
ATGGAGGTTGCAATGGTGAGTGCGGAGAGCTCAGGGTGCAACAGTCACATGCCTTATGGT
TATGCTGCCCAGGCCCGGGCCCGGGAGCGGGAGAGGCTTGCTCACTCCAGGGCAGCTGCA
GCAGCTGCTGTTGCAGCGGCCACAGCTGCTGTCGAAGGTAGCGGGGGTTCTGGTGGGGGC
TCCCAC
>gene3
ATGGAGGTGGCGATGGTGAGTGCGGAGAGCTCAGGGTGCAACAGTCACATGCCTTACGGG
TACGCGGCCCAGGCCCGGGCCCGGGAGCGGGAGAGGCTGGCTCACTCCCGGGCGGCGGCG
GCCGCCGCCGTCGCGGCTGCCACGGCTGCCGTGGAAGGCAGTGGGGGGCCTGG
蛋白質(zhì)序列pep.fa
>gene1
MEVAMVSAESSGCNSHMPYGYAAQARARERERLAHSRAAAAAAVAAATAAVEGSGGSGGG
>gene2
MEVAMVSAESSRCNSHMPYGYAAQARARERERLAHSRAAAAAAVAAATAAVEGSGSSGGGSH
>gene3
MEVAMVSAESSGCNSHMPYGYAAQARARERERLAHSRAAAAAAVAAAKAAVEGSGGP
注意：cds.fa和pep.fa文件的序列ID號（gene1,2,3）要一致。
2 對蛋白質(zhì)序列pep.fa進行比對
進行蛋白質(zhì)序列比對前，需要先安裝mafft軟件。
下載mafft軟件：
wget https://mafft.cbrc.jp/alignment/software/mafft-7.453-with-extensions-src.tgz
tar -zxvf mafft-7.453-with-extensions-src.tgz
cd mafft-7.453-with-extensions/core
安裝：
1）有root權(quán)限用戶安裝法：
make clean
make
su
make install
2）無root權(quán)限用戶安裝法：
vi Makefile
進入makefile文件后，修改第一行和第三行，如下所示兩張圖，分別為修改前和修改后（請務(wù)必修改你有權(quán)限的路徑）：
安裝成功后，輸入命令mafft --auto pep.fa > alig_pep.fa
生成的alig_pep.fa文件如下：
3 將比對好的蛋白質(zhì)序列與cds序列比對
這一步需要下載pal2nal.pl文件
wget http://www.bork.embl.de/pal2nal/distribution/pal2nal.v14.tar.gz
tar -zxvf pal2nal.v14.tar.gz
cd pal2nal.v14/
下載后就能看見pal2nal.pl這個文件
隨后將蛋白質(zhì)序列與cds序列比對
pal2nal.pl alig_pep.fa cds.fa -output fasta > cds_pep_aln.fa
比對好的cds_pep_aln.fa文件如下所示：
4 生成計算kaks值時的基因?qū)?/p>
計算kaks值前需要先將cds_pep_aln.fa文件拆分：
csplit cds_pep_aln.fa /\>/ -n2 -s {*} -f gene -b "%1d.fa" ; rm gene0.fa
拆分后，會生成gene1.fa?、gene2.fa、?gene3.fa三個文件。
接下來，將生成的gene1.fa、?gene2.fa、?gene3.fa組成新的基因?qū)ene1-gene2、gene1-gene3、gene2-gene3，見如下命令：
fori in $(seq13)docatgene1.fa gene${i}.fa > gene1_${i}.fadonecatgene2.fa gene3.fa > gene2_3.fa
生成如下幾個文件：
gene1_1.fa
gene1_2.fa
gene1_3.fa
gene2_3.fa
其中，gene1_2.fa、gene1_3.fa、gene2_3.fa為我們所需的基因?qū)Α?/p>
這里將他們提取成基因?qū)?，而不是多條序列進行計算的原因是，KaKs_Calculator這個軟件在處理多序列的輸入文件時，會報錯：Error. The size of two sequences in 'gene1&gene2' is not equal。我嘗試了很多遍，發(fā)現(xiàn)只能提取成基因?qū)Σ挪粫筮@種錯誤。當然，偶爾運氣好的時候，KaKs_Calculator是能處理多序列的kaks值的，為了防止出錯，我們還是將他們拆開計算好一點。
5 將gene1_2.fa、gene1_3.fa、gene2_3.fa文件轉(zhuǎn)化為axt文件
轉(zhuǎn)化為axt文件需要下載parseFastaIntoAXT.pl文件，文件地址：https://gitee.com/liaochenlanruo/kaks_pupline/blob/master/parseFastaIntoAXT.pl
下載后，輸入如下命令:
foriin*.fadoecho$iperl parseFastaIntoAXT.pl$idone
生成三個文件：
gene1_2.fa.axt
gene1_3.fa.axt
gene2_3.fa.axt
6 計算kaks值
下載安裝kakscalculator
下載鏈接：https://sourceforge.net/projects/kakscalculator2/
輸入以下命令：
foriin*.fa.axtdoecho$iKaKs_Calculator -i$i-o${i%%.*}.kaksdone
生成三個文件：gene1_2.kaks、gene1_3.kaks、gene2_3.kaks
到這一步，批量計算kaks值的工作就已經(jīng)搞定。
這里附上安裝kaks_calculator軟件可能會遇到報錯：
g++ -O -o AXTConvertor AXTConvertor.cpp -lstdc++ -lm
AXTConvertor.cpp: In function ?.ool readPhylip()?.
AXTConvertor.cpp:210:22: error: ?.toi?.was not declared in this scope
if (atoi(num.c_str())!=sequence.size()){
AXTConvertor.cpp: In function ?.ool readNexus()?.
AXTConvertor.cpp:338:39: error: ?.toi?.was not declared in this scope
if (sequence.size()!=atoi(num.c_str())) >{
make: *** [AXTConvertor] Error 1
解決方法在這里：
編輯KaKs.cpp文件，加上include "string.h"
編輯AXTConvertor.cpp文件，加上include "stdlib.h"
編輯GY94.cpp文件，加上?include "string.h"
如無報錯請忽略上述內(nèi)容。
