R包辣鸡之CorMut
曾经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比较，数数字你会吗？
:
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值你会算吗？
:
LOD=0Nnum=NS+NYfor(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简直瞎眼了！