双序列比对 | 尝试数次，调试千遍
写在前面
这两天，“Molecular Plant - 华南农业大学作物科学前沿暨中衡山论坛”在华南农业大学举办，地点是学校红满堂，亦即华南农大最高规格的报告厅。做报告的均是植物学领域的大牛级人物，报告各个精彩。在此期间，我有幸与樊龙江老师攀谈，同时也要来了签名。
其中《生物信息学》对应的是樊老师 2001 年编写的《生物信息学札记》一书。该书为电子书。我相信可能国内绝大多数接触生物信息的朋友都可能翻阅过。事实上，小时候（应该也是研二研三，接触生信时），我也翻了。后来，我运营《生信札记》微信公众号。尽管我现在只记得是临时想的名字，但多少应是受樊老师该本电子书书名影响。《生物信息学札记》多次更新迭代，我印象我看到时似乎是第4版。其整理成纸质版，即《生物信息学》一书。
昨天半夜我在公众号发了图片后，不少朋友问起哪里可以购置该书。这里给出答复：
《植物基因组学》一书是2020年1月份出版，内容还是紧靠领域发展前沿，可以到各大图书网站购买；
《生物信息学》一书，樊老师与我提及，第二版将会在下个月（2020年05月）正式出版，我个人建议大伙 5月份直接购买新版。书出来了，我会第一时间在公众号上推送。
回到主题
与樊老师一起聊到了生物信息学这二十多年来的发展变化，也请教了生物信息学教学，更或者学习上的一些问题。非常碰巧，我们回到的最经典的一个生信编程入门题：Needlman-Wunsch算法，即 双序列全局比对。对于这个算法（以及SW），我个人还是有比较深的感触：
从算法原理或理论来说，讲到底就是动态对话，打分矩阵；
从程序实现来说，就是开数组 或 做递归。
可以说，看起来是一个“说起来简单，学起来简单，做起来也简单”的事情。但这并非事实。
为什么一定要自己写双序列全局比对
很久以前，我写过一个推文《高速！大规模 (10w+) 基因对 Ka/Ks 计算工具 - TBtools》，推文中，我提到了一个场景，即需要计算数十万个基因对的Ka/Ks值，这意味着需要做大量双序列比对。经过我的 Emboss软件套件的Needle Man是表现最佳的双序列比对工具，但他并不支持规模化的指定双序列比对，只能逐对调用。所以我需要另找路子：
看了下Biopython等，本质是调用 Emboss Needleman，效率必然是低下，甚至是最低下的。
BioJava，似乎找不到支持 Gap Extension 的实现（似乎很多 BioXXX 库并没有做足够的优化，只能确保有这个功能，却并不保证效率和准确度）
自己写
这里继续用下之前推文的截图。
从打分上来说（这里只测序了十来对序列），之前我写的版本效果仅次于Emboss Needle，但明显优于MUSCLE 和 Biojava，尤其是后者。Biojava的实现似乎有点问题（Manual Edit意指手动给他矫正了一部分明显比对错误的部分），具体不清楚。
直接代码原生实现，我们可以达到很快的速度
biopython的emboss调用实现，一秒钟跑一对（注意到，调用一般开多线程也快不起来）而我们自己写的，放在TBtools里了单线程 1 min 10000个基因对，换句话说，一秒钟大概跑 167 对....四线程 ~20 s 10000个基因对，多线程更快....权衡之下，只能选择自己写。期间我也看了不少资料，可以说，几乎所有资料，都只会提及Gap Open（开Gap），根本不会提及Gap Extension（Gap延伸）。这事完全可以理解。作为经典的动态规划在生信上的应用，只讲Gap Open，完全可以把最直接的原理呈现出来，比如 wiki 上的。
在没有 Gap Extension 的情况下，延伸发生=开Gap，
按照上图，我们可以很快写出程序，因为逻辑简单，只需要开一个数组来打分就可以了。前一步的最优状态直接确定，从而每一步都可以考虑。而一旦考虑 Gap Extension，此时上一步如果是Gap，那么这一步可以是Gap Extension，两步累计分数存在一些情况会正好高于上一步非Gap而这一步Open Gap的状态，亦即，上一步的几个状态均需要保存。当然，想明白这一点就简单了，只需要开多两个数组，分别存储横向和纵向的Gap即可。
是的，以前我没想明白，昨天听完报告，想想就再写写。发现还有一些问题。调试了一会，终于还是解决了，说到底，具体算法实现时，最麻烦的或许边界处理。正如我以前写 NG86 KaKs计算逻辑一样。文献一读就懂，但具体代码实现时还是会发现不少细节处，边界处需要自己思考和处理。
尝试数次，调试千遍
说实话，我从来没想过这个哪哪都能翻阅到学习资料的基本算法，真正动手去写有这么麻烦。
可以看到的记录是四次，前前后后跨度不少于三年。事实上，我有时想起来就会重写一遍，换个实现逻辑，再写一遍，远不止四次。不过以后不会再写了，这一次搞定了。
测试了一下，拿 1w 对同源蛋白序列，长度在450~500aa，运行时间是：2min。比之前慢了一倍，可能是我一次开了四个数组。Anyway，不关心这个。
整体上，我还是比较满意。毕竟现在已经完全达到和Emboss Needle一样的准确度，但以后相关功能实现上，TBtools会有绝对的速度优势。
写在最后
Emmm，我觉得似乎可以进行一波 Ka/Ks 批量计算大规模宣传，毕竟之前 TBtools 上的多线程运算使用的双序列比对效果准确度尽管高于Muscle，但不及Emboss，现在可以了。接下来，TBtools是最准的，也是最快的。
至于，这个双序列比对界面化接口的开放，那就回头再说吧。
最后，这是樊老师主编的第一版的《生物信息学》（意味着不会再印刷，感兴趣的还是可以购买收藏的，推荐当当....不要问为啥），而下个月，会有第二版出来。