曾经QQ上有人叫我帮忙，跑某个R包做个分析，我一看那包，一堆bug，显然不能用，那就不是简单的事了，我可没空帮他写代码，于是他就经常恶心我，有事没事就来说我写ggtree没意义，不是实质科研，不如跟他做点牛逼的，不如再写个R包干他那事。昨天用马甲在进化群里问画树，我还热心贴个代码给他，说看不懂，我在群里说了，我写个文详细介绍一下。然而他的马甲身份暴露啦，所以今日跳票，我写好的文也不发了，你牛逼就不要用马甲来骗代码呀，ggtree没用你不要用呀，叔就是这么任性😎

之前讲过某个R包我一看一堆bug，直接放弃，今天倒是拿出来晒一下，不为别的，我只想说一句，一知半解是很危险的，盲目相信软件也是很危险的。

比对？不比对？这是个严肃的问题

话说此包的输入是alignment，我刚开始也没看代码，照着示例跑一下，但读序列的时候，我就惊呆了，怎么读个fasta这么慢，于是我开始看代码，不看不知道，一看吓一跳：

seq_f01=as.list(sapply(allSeq,function(x)c2s(s2c(x)[ref01!="-"])))

这个s2c把string变成character vector，c2s则反过来，这一句就干一件简单的事情，就是把所有string（fasta读进来的sequences）中的-全部删除。就这么一切一贴，非常非常慢，简单地用gsub(“-“, “”, allSeq)就可以非常快速做同样的事情。 然后这时候我们似乎觉得那里不妥！你把-全删了，等于是把比对序列变成了非比对序列。那么你要求输入是比对序列是为什么？！ 确实序列必须得是比对好的，因为这包就是要检测correlated mutations，它比的是codon, 这必须是比对好的序列才行，但在读序列的时候，就把-给删了，序列不再是比对好的，本来对应的codon现在必然很多都移位了，这在一开始就挖了巨大的坑。后面所有的计算必须都是辣鸡了。

codon比较，数数字你会吗？

https://github.com/Bioconductor-mirror/CorMut/blob/master/R/CorMut-internal.R#L125-L129:

for(j in 1:length(seqf[[i]])){
	if(j%%3==1){ tris=seqf[[i]][j:(j+2)]
	}else if(j%%3==2){ tris=seqf[[i]][(j-1):(j+1)]
	}else if(j%%3==0){ tris=seqf[[i]][(j-2):j]
}

这循环竟然不是一个codon一个codon来，而是一个base一个base地迭代，这里移动三个位置，拿到的都是同一个codon，然后它就去比较是否不同，数数目了，一个不同，数它3遍，我也是醉了。

p值你会算吗？

https://github.com/Bioconductor-mirror/CorMut/blob/release-3.4/R/kaksAA.R#L67-L77:

LOD=0
Nnum=NS+NY
for(ii in NYs:Nnum){
	per_item=choose(Nnum,ii)*(q^ii)*((1-q)^(Nnum-ii))
	#avoid the appear of NaN
	if(per_item!="NaN"){
		LOD=per_item+LOD
	}
	#print(choose(dim(mat)[2],i)*(q^i)*((1-q)^(dim(mat)[2]-i)))
}
LOD=-log10(LOD)

我们知道算p值是算积分面积，它倒好，离散性分布这一个个数字来计算，然后求合。 这里用的是二项式分布：

per_item=choose(Nnum,ii)(q^ii)((1-q)^(Nnum-ii))

这个基本上你算出来会有很多0和NaN，精度问题嘛，这样算必须是不准的。然后他必须是跑了自己代码出不来p值，发现有很多NaN，也不去想为什么，直接把NaN给扔了。这根本就是错的，也不知道用pbinom去算，我也是服了。

总结

我只是粗略地看了一下代码，就发现很多错误，而且几乎每个重要的部分都出错。问题肯定比我粗略看一下所能发现地多得多。

  • 读序列，把比对的序列搞成了非比对。
  • 比较codon，序列非比对这比毛线啊，然后连数数目都数错。
  • p值，简单二项式，你都不会算啊，简单粗暴地对NaN视而不见。

这货中国CDC出品，竟然还发了Bioinformatics，reviewer简直瞎眼了！

​