Solved! | 为什么我 Ka/Ks 计算出来的总是 NaN?
写在前面
很久以前,参考NG86计算逻辑提出的文章,我在TBtools中写了Simple Ka/Ks Calculator,旨在理解算法,具体可见 《批量计算dn/ds,粗略筛选正选择?所有人-技能Get !》。近期,不时有人在 TBtools 群提出,TBtools 的 Simple Ka/Ks Calculator 计算出来的有一堆 NaN。Emmm,这是一个比较有趣的事情。有趣之处有二:
TBtools中几乎每一个功能,我都会在从头实现算法之后,拿一些软件的运行结果进行比较。针对 Ka/Ks 的计算,我拿了 DNAsp 和 KaKs_Calculator2 的结果与TBtools的结果进行比较,无出入。
曾经有人提出 TBtools 的计算存在问题,经过一堆比较,TBtools的结果反而是对的。具体可见推文《Ks值用TBtools计算不出来?讲一个故事》。
但是今天的情况有所不同。一般十来个基因对,那么有两三个NaN可以理解。如果全部都是 NaN,那就要引起重视。如
不是用户使用问题
我要来数据,随后自己使用 TBtools 最新版本进行计算,结果与该用户的结果一致。
说明不是用户使用问题。
不是 TBtools 实现逻辑问题
既然如何,那么就有可能是 TBtools 的算法实现逻辑问题。为此,我直接让 TBtools 输出 AXT 格式的比对结果。随后使用 章张老师 的 KaKs_Calculator 计算,结果如下。
换句话说,至少在 Ka/Ks 的计算上,TBtools 完全没问题。
那么还存在一个可能,是不是 Codon Alignment 出现了问题?于是我拿了一个基因对,使用 Mega 进行 Codon Alignment,随后导出结果。再使用 KaKs_Calculator 进行结算,
结果与前述一致。
序列距离(差异)太大,导致计算出错
回到 TBtools 计算源码,在 Java 的数值计算中,除非是NaN参与运算,才会得到 NaN。那是咋回事?
重新过了一遍,发现问题应该是在 JC69 模型应用上。
当 pS 大于等于 3/4 的时候,计算结果就会得到负值。换成生物学角度,那就是...绝大部分可发生同义突变位点都发生了同义突变,亦即,序列分歧度太大,进化距离很远。
为此.... 我继续对 TBtools 做个改进,当计算出来的 pS>=3/4,补充一个信息,说明原因...那就不会有用户问我了。
写在后面
Emmm,有时候我觉得,之所以总是有用户提出问题,是因为大家能直接接触到我。毕竟为什么同样的问题,似乎没有人会去问其他软件开发的人呢?还是因为我不知道?
或者说,我太贴心了?哈哈哈哈哈哈哈啊哈哈哈啊哈。
尽管人生如戏,还是要诚以待人,与人为善,投之以桃,而报之以李。