1 介绍

参考文献：Zhang Z . KaKs_calculator 3.0: Calculating selective pressure on coding and non-coding sequences. 2022.
KaKs_Calculator 3.0能够计算编码和非编码序列的选择压力。与编码序列的非同义/同义替换率比类似，对非编码序列的选择可以量化为非编码核苷酸替换率与相邻编码序列的同义替换率的比值。

2 安装

Linux版本：https://ngdc.cncb.ac.cn/biocode/tools/1/releases/3.0/file/download?filename=KaKs_Calculator3.0.zip 通过wget命令下载
使用以下命令解压 KaKs_CalculatorXXX.tar.gz 的包。
unzip KaKs_CalculatorXXX.zip
使用 g++/gcc 编译器编译源代码文件夹中的程序。
cd KaKs_CalculatorXXX/src
make
对编码序列运行“KaKs”，对非编码序列运行“KnKs”。

windows版本：https://ngdc.cncb.ac.cn/biocode/tools/1/releases/3.0/file/download?filename=KaKs_Calculator3.0-Windows-setup.zip
对于 Windows，名为“bin”的文件夹中有两个已编译的可执行文件。

• 进入名为“GUI”的文件夹。
• 单击设置应用程序并安装 KaKs_Calculator 3.0。
• 请在“开始菜单”中找到 KaKs_Calculator。
  ! 值得注意的是，在 Windows 10 上测试的两个示例文件（coding.axt 和 noncoding.axt）位于以下文件夹中。
  C:\Program Files (x86)\CNCB-NGDC\KaKs_Calculator

3 使用

使用手册：https://ngdc.cncb.ac.cn/biocode/tools/BT000001/manual
以下使用Linux版本的 KaKs_Calculator

计算选择压力有两种方法：1近似法（更省时）；2最大似然法（准确度更高）

计算 Ka 和 Ks 通常包括三个步骤。让我们假设比较的两个 DNA 序列之间的长度数为 n，它们之间的替换数为 m。要计算 Ka 和 Ks，我们需要计算同义 (S) 和非同义 (N) 位点的数量 (S + N = n) 以及同义 (Sd) 和非同义 (Nd) 替换的数量 (Sd + Nd = m ）。然后在校正多个替换之后，(Nd/N) 和 (Sd/S) 可以分别代表 Ka 和 Ks，因为观察到的替换数低估了实际的替换数，因为序列随着时间的推移而发散。因此，我们可以从上面得出结论，这些方法通常包括三个步骤来估计 Ka 和 Ks：

• 计算 S 和 N
• 计算 Sd 和 Nd
• 多次替换的校正

3.1 蛋白质序列的获得

我们使用日本晴（后用Nip简称）、梗稻（后用Jap简称）及普通野生稻（后用Ruf简称）的蛋白质数据。

3.2 蛋白质比对并反翻译

使用章张老师团队开发的软件——ParaAT。该软件可以平行比对蛋白质序列并将比对的序列反向翻译成相应的 DNA 比对。
下载地址：https://ngdc.cncb.ac.cn/tools/paraat

示意图
ParaAT 整合了计算ka/ks所需的一整套分析，包括:

1. 蛋白序列比对（可选 clustalw2 | t_coffee | mafft | muscle）
2. 根据蛋白比对结果回译成codon对应的核酸比对结果
3. 计算kaks值（KaKs_Calculator实现）
   常用的kaks计算软件如paml和kaks_calculator需要准备的文件和指令都比较复杂，上手困难且容易出现问题，使用ParaAT能自动批量准备所需要的计算文件，十分便捷。

Q: 如果想在Windows下运行ParaAT 怎么办？
A： 有以下四个步骤：

• 安装 Cygwin，以运行 Linux/Unix 等程序；
• 在 Cygwin 中，安装多序列比对器（例如 ClustalW）；
• 使用“chmod”设置多序列比对器、Epal2nal.pl和ParaAT.pl的执行权限；
• 设置全局目录并确保 ParaAT.pl、Epal2nal.pl 和多序列比对器可从任何工作目录访问。

3.2.1 运行前准备工作

• 确保ParaAT.pl、Epal2nal.pl及多序列比对软件(clustalw2 | t_coffee | mafft | muscle)在环境变量中，以便使用

vim ~/.bashrc  #或者~/.bash_profile
按i进入编辑，在末尾追加写
export PATH=$PATH:$HOME/Path/To/File/ParaAT.pl
export PATH=$PATH:$HOME/Path/To/File/Epal2nal.pl
按Esc退出，输入:wq保存并退出

• 蛋白质序列文件、DNA序列文件和同源基因对文件的序列ID要保持一致！！！
• 需要的文件：同源基因列表；fasta格式的氨基酸序列文件和核苷酸序列文件；多线程运行文件（输入一个线程数就可以了）

3.2.2 OrthoFinder寻找同源基因

安装 conda install orthofinder
常规运用

1. 将目的物种的蛋白质序列文件都放在一个文件夹

orthofinder -f <dir>

得到的文件是orthogroups文件夹里的txt文件和csv文件。用txt比较好，因为paraAT的案例文件homolog需要用空格隔开。
因为得到的结果有很多都是自己和自己的同源基因，所以我做了2个版本——筛选出 只有3个同源基因的结果 和 未筛选结果。未筛选结果是得不到axt文件的，估计是因为同源基因对太多了（？）
如果数据量很大的话，建议用脚本来做，excel不合适。

f = open('Orthogroups.txt','r')  # 文件需要和脚本文件同一目录下，要不然需要写路径
fi = open('Orthogroups_revise_all.txt','w')
file = f.readlines()

for line in file:
    ls = line.split(': ') 
    if len(ls) == 3:
         fi.writelines(s) # 未筛选序列，只是把第一列的内容删除，要不然paraAT识别不了
f.close()
fi.close()

3.2.3 运行ParaAT

3.0版本中，把KaKs_Calculator加入环境变量时，是把/src/目录放进去，而不是/bin/（/bin/的是Windows在用）。但是在后续的ParaAT.pl的使用时，一直说KaKs_Calculator没有在环境变量，这是因为KaKs_Calculator 3.0的命令是“KaKs“，而ParaAT的“-k”命令使用的是“KaKs_Calculator”（2.0版本的命令是“KaKs_Calculator”）。

3.0是2022年发表的，而ParaAT是2012年发表的软件不适配吗55。所以我没有用ParaAT -k 来一体化流程。

命令

ParaAT.pl -h test.homologs -n test.cds -a test.pep -p proc -m muscle -f axt -g -k -o result_dir #proc文件必须与输出位置在同一个目录下，不然会报错

参数：
-h, 指定同源基因列表文件
-n, 指定核酸序列文件
-a, 指定蛋白序列文件
-p, 指定多线程文件
-m, 指定比对工具
-g, 去除比对有gap的密码子
-k, 用KaKs_Calculator 计算kaks值
-o, 输出结果的目录
-f, 输出比对文件的格式

但是得到的axt文件全都是野生稻之间的
原因是：三个文件的序列名格式要一样，不要有注释信息。蛋白质文件也不要有多余的符号。

Nip的序列名
Ruf的序列名
把注释信息去掉再做一次就可以了。

运行KaKs

常规运用

KaKs -i test.axt -o test.axt.kaks

参数
-i 输入axt文件
-o 输出文件
-m 选择的模型（默认MA）
-c 遗产密码表（默认1-标准代码）

因为有很多axt文件，所以用个脚本跑吧~

报错
报错：“[Error. The sequences are not equal in length.]”
原因：（参考https://www.omicsclass.com/article/700）

1.序列比对总长度是不是3的整倍数，如果不是3的倍数可能计算不出结果；
2.看看两条序列是否都含有起始密码子ATG，如果没有可能也会报错。
3.该软件只能算2条序列之间的kaks，但是它不会检查输入是几条序列。因为axt文件这个奇怪的输入格式，也没法检查。除了第一行，输入文件里面所有的序列是一起被读成一个字符串。然而，它会检查“两”条序列是不是一样长。方法是这样：

检查序列方式 所以，如果你输入的3条序列都是偶数长度，就pass了；如果是奇数长度，就报错......但是对于pass的那些，问题就大了。你得到了一个结果，但是这个结果比较的实际上是：
sequence_1 + sequence_2的前半段，和
sequence_2的后半段 + sequence_3

4 KaKs_Calculator 2.0和paraAT一起使用

下载：https://sourceforge.net/projects/kakscalculator2/
使用：

1. 安装：打开src/目录，输入make命令。
   有可能会报错：KaKs.cpp:756:6: error: ‘strlen’ was not declared in this scope
   这时，打开出错的那个文件，在开头加一行：#include<cstring>，保存并退出，重新make即可。
   报错
2. 把2.0加入环境变量。
3. 使用paraAT

ParaAT.pl -h Orthogroups.txt -n rice.cds -a rice.pep -p proc -m muscle -f axt -g -k -o ./axt_final/

5 总结

总体来说，KaKs_Calculator 2.0与paraAT结合使用更方便，能一体化流程。但3.0比2.0来说，功能、准确性等都有优化，只是要分2步走，麻烦了些，得到的结果是一样的。希望章老师也能尽快将paraAT适配3.0吧~
但是！KaKs_Calculator只能计算2条序列之间的kaks，3条以上的可以尝试用codeml或者hyphy等软件。