高速！大规模 (10w+) 基因对 Ka/Ks 计算工具 - TBtools

写在前面

目前还是在继续整一些擦边“比较基因组”的工具。前段时间停下来，原因只有一个：我发现 TBtools 的Simple Ka/Ks Calculator 速度确实不够快。在分析的过程中，我们常常会对两类基因对集合进行 Ks 计算：

共线性分析结果：比如 MCScan 和 I-ADHoRe 的分析结果，这系列软件已经挖掘了共线性，于是得到的基因对本身也不是很多，往往是 w+ 级别
blastP 结果：即相对初步很原始的，具有一定相似性的基因对，未经过共线性分析，往往是 10w+ 级别

无论哪一种，要对全部基因对进行 Ka/Ks 计算，那么对于普通 PC 来说，计算量也不是太小。主要计算分两步：

双序列的codon alignment，一般批量分析时，先做 protein alignment ，随后基于 CDS 转换为 codon alignment
进行 ka/ks 计算，其中在两基因间最快的最方便的最常用的是 NG86

很明显，看起来都很简单。但具体实践是另一回事。

旧版本的 Simple Ka/Ks Calculator

早前，我在推文已经提及，写这个功能最原始的目的，还是更好的理解大体计算逻辑。当然写完之后，我还是写了推文，具体见《批量计算dn/ds，粗略筛选正选择？所有人-技能Get !》，自认为还可以。后面课题组做相关分析其实也用上了。
当初的实现比较简单，双序列比对上，直接调用 muscle （大家都清楚，TBtools 打包了muscle）。整体效果其实还可以。就是：

1w+ 基因对，需要 ~20 min
10w+ 基因对，需要 ~4 hour

从某个角度来说，静置过夜似乎是不错的办法。尽管往往只有我们做实验的时候才会这么干。

增强版本

事实上，频繁地通过命令行调用外部程序，对于程序的开销极其大。于是，最好的解决办法是直接用原生的双序列比对逻辑。为了保证自己写的逻辑正确，我反反复复其实间断鼓捣果几次，整体周期应有一年以上，其实这不能怪我，具体如下：

EMBOSS Needleman-Wunsch，表现最佳，调用needle，主要作为准确参考
Muscle，多序列比对软件
mafft，多序列比对软件
clustal w2
biojava
我自己写的

具体得分如下

事实上，我是有点意外的。当然，我后续只做了 100 来个基因对。needle 表现最优，这个跟它本身软件开发目的有关。其次是我写的（其实我写的也是 needle 算法，但是我也发现了有一点瑕疵），随后才是 biojava 和 muscle。

各类 bioXXX 库当然很优秀，我也相信他的实现其实很好。具体没去看源码。但这一年来我并不是没有尝试在全网检索，看看是否有表现与 emboss needle 一样好的。结果是，没有。这个完全可以理解。
双序列比对是很多人的入门算法，尤其是动态规划。大家最多考虑是 gap+extension 的罚分。但是序列足够多，足够长的时候，常常会遇到连续的双元选择，这一步最优的结果，可能会导致下一步的局部最优。当然，这个跟具体打分实现有关（后续我会继续优化）。Emmm，简单的实现，到处都是，毕竟大家目的是理解动态规划的概念，而不是写一个着实可用的Code。

既然测试结果，我自己写的并不亚于 biojava，更不亚于原始的调用 muscle 版本，那么就可以做增强。
基于测试：

TBtools 原始版本单线程计算 1.6w 个基因对，20min
改进版本单线程计算 1.6w 个基因对，3min

多线程

市面上应是几乎找不到 PC 或者 laptop 还是只有单线程的CPU的... 多线程从某种角度可以加速一些分析任务。
在使用我自己写的双序列比对逻辑之前，我想到的加速方案是多线程。
但测试结果如下：

单线程 10000个基因对，12 min
四线程 10000个基因对，15 min

开了多线程，反而更慢，甚至比 15 min 还久。这个其实很好理解，因为开辟线程是需要开销的，另外本身就是一个重IO的过程，调用muscle中设计系列IO操作，尤其是程序之间的通讯。Sad。

不过，现在既然用自己写的双序列比对逻辑，那么就会得到

单线程 10000个基因对，1 min
四线程 10000个基因对， ~20 s

开心！完美，不仅单线程加速了，还可以使用多个线程，从某个角度来说，并行再加速...当然，线程不是越多越多....

增强版本的 Simple Ka/Ks Calculator

相比于旧版本，增加了一个 CPU 参数：1）保持为 0 ，那么仍然会自动使用旧版本的计算逻辑，调用muscle，速度慢；2）调整为 1 或者更高，比如 4，那么会使用我自己写的逻辑，速度快

使用实例

既然现在，我们可以进行 10w+ 基因对级别的 Ks 快速计算（也就10min），那么可以做的事情就很多，比如，看看香蕉基因组的 Ks Plot，（注：直接走一个蛋白全集BlastP，计算Ks 即可，无需经过共线性分析，如MCScan等）

看看多年前的香蕉基因组 Paper

感觉还不错，于是我们还可以看看拟南芥的

其实也很好。

还有更多

比如，在拟南芥上，你可以画一个原始的Blast结果，同时还画一个共线性分析之后的结果
大体参数

注：前面我们看到了，Ks分布大体就是两个Peaks，0.6 和 1.8。

总的来说，共线性分析之后确实看起来清晰了很多，尽管仁者见仁。

榕树基因组的情况

正好前几天测试 “Find Best Homology” 的时候，下载了榕树基因组文章的数据。那就跑跑看。
看看榕树基因组文章原本的点图

然后再看看 TBtools 目前版本的点图（整体跑下来，从输入数据到输出，大概 20 min）

比较一下，就会发现，基因组文章原图中，可以非常明显的看出两个榕树物种的共线性情况（应是进行了过滤，只相似度最高的基因对进行可视化，但具体就不清楚了）。而 TBtools 输出的这个图，其实也还可以。
拿到这个图，我一眼下去是觉得有多倍化事件的痕迹。于是找了课题组的沉睡的小五郎，聊了下这个图，也聊了下一些基因组分析的杂事。
针对这个图，我们一致的观点是：图中显示的应该似乎双子叶共有的古三倍化事件。可以看到，两个榕树物种基因组共线性很好（蓝色点线），同时存在一个基因又另外对应两个位置（红色点先）。
通过 Ks 着色，可以大体判断基因之间的分化时间。事实上，我是挺想画出一个还不错的图的。区分不同时期的倍化事件。不过香蕉基因组似乎不给我机会。

勉强调了下出图参数

多少可以看得出来，香蕉基因组有两次独立相对近期的全基因组复制，其中一个比较老（灰色），一个比较近（蓝色），当然，还有一次更老的，可以看到都是红色的点线。

写在后面

天下武功，唯快不破！
香蕉基因组注释，确实比较差。开展相关研究，或许得先从基因结构注释矫正做起。
写到这里，已经忘了为啥要写了。
不过，话说回来，如果没有完善，简便，高效的流程，我相信我几乎不可能在一两个小时内做完以上的分析，而且还可视化。同时，写了这个推文。