超级详细的ka、ks分析
1 文件准备
首先准备两个文件,一个是基因的cds序列,一个是蛋白质序列。
cds序列和蛋白质可以在ensembl网站找到:http://ftp.ensembl.org/pub/current_fasta/
这两个文件的示例如下:
cds序列文件cds.fa
>gene1
ATGGAGGTTGCAATGGTGAGTGCGGAGAGCTCAGGATGCAACAGTCACATGCCTTACGGT
TATGCTGCCCAGGCCCGGGCCCGGGAGCGGGAGAGGCTTGCTCACTCCAGGGCAGCTGCG
GCAGCTGCCGTTGCAGCGGCCACAGCTGCCGTCGAAGGAAGTGGGGGTTCTGGTGGGGGG
>gene2
ATGGAGGTTGCAATGGTGAGTGCGGAGAGCTCAGGGTGCAACAGTCACATGCCTTATGGT
TATGCTGCCCAGGCCCGGGCCCGGGAGCGGGAGAGGCTTGCTCACTCCAGGGCAGCTGCA
GCAGCTGCTGTTGCAGCGGCCACAGCTGCTGTCGAAGGTAGCGGGGGTTCTGGTGGGGGC
TCCCAC
>gene3
ATGGAGGTGGCGATGGTGAGTGCGGAGAGCTCAGGGTGCAACAGTCACATGCCTTACGGG
TACGCGGCCCAGGCCCGGGCCCGGGAGCGGGAGAGGCTGGCTCACTCCCGGGCGGCGGCG
GCCGCCGCCGTCGCGGCTGCCACGGCTGCCGTGGAAGGCAGTGGGGGGCCTGG
蛋白质序列pep.fa
>gene1
MEVAMVSAESSGCNSHMPYGYAAQARARERERLAHSRAAAAAAVAAATAAVEGSGGSGGG
>gene2
MEVAMVSAESSRCNSHMPYGYAAQARARERERLAHSRAAAAAAVAAATAAVEGSGSSGGGSH
>gene3
MEVAMVSAESSGCNSHMPYGYAAQARARERERLAHSRAAAAAAVAAAKAAVEGSGGP
注意:cds.fa
和pep.fa
文件的序列ID号(gene1,2,3)要一致。
2 对蛋白质序列pep.fa进行比对
进行蛋白质序列比对前,需要先安装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权限用户安装法:
make clean
make
su
make install
2)无root权限用户安装法:
vi Makefile
进入makefile文件后,修改第一行和第三行,如下所示两张图,分别为修改前和修改后(请务必修改你有权限的路径):
图片.png图片.png
安装成功后,输入命令
mafft --auto pep.fa > alig_pep.fa
生成的alig_pep.fa文件如下
图片.png
3 将比对好的蛋白质序列与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
这个文件
随后将蛋白质序列与cds序列比对
pal2nal.pl alig_pep.fa cds.fa -output fasta > cds_pep_aln.fa
比对好的cds_pep_aln.fa
文件如下所示:
4 生成计算kaks值时的基因对
计算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
组成新的基因对gene1-gene2
、gene1-gene3
、gene2-gene3
,见如下命令:
for i in $(seq 1 3)
do
cat gene1.fa gene${i}.fa > gene1_${i}.fa
done
cat gene2.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
为我们所需的基因对。
这里将他们提取成基因对,而不是多条序列进行计算的原因是,
KaKs_Calculator
这个软件在处理多序列的输入文件时,会报错:Error. The size of two sequences in 'gene1&gene2' is not equal。我尝试了很多遍,发现只能提取成基因对才不会报这种错误。当然,偶尔运气好的时候,KaKs_Calculator
是能处理多序列的kaks值的,为了防止出错,我们还是将他们拆开计算好一点。
5 将gene1_2.fa
、gene1_3.fa
、gene2_3.fa
文件转化为axt文件
转化为axt
文件需要下载parseFastaIntoAXT.pl文件,文件地址:https://gitee.com/liaochenlanruo/kaks_pupline/blob/master/parseFastaIntoAXT.pl
下载后,输入如下命令:
for i in *.fa
do
echo $i
perl parseFastaIntoAXT.pl $i
done
生成三个文件:
gene1_2.fa.axt
gene1_3.fa.axt
gene2_3.fa.axt
6 计算ka/ks值
下载安装kakscalculator
下载链接:https://sourceforge.net/projects/kakscalculator2/
输入以下命令:
for i in *.fa.axt
do
echo $i
KaKs_Calculator -i $i -o ${i%%.*}.kaks
done
生成三个文件:gene1_2.kaks
、gene1_3.kaks
、gene2_3.kaks
到这一步,批量计算kaks值的工作就已经搞定。
这里附上安装
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"
如无报错请忽略上述内容。