植物甲基化位点状态的判定（bismark之后，call DMR之

狂躁的两天，什么都不想干。写点东西吧。

甲基化原始数据在质控和mapping之后，需要进行状态的判定，也就是说extract出的这个位点是不是真的发生了甲基化，网上几乎没有教程，我也摸索了好久，请教了一些大神，分享出来给大家共勉。

这一步需要用二项检验进行统计学分析，然后用FDR去校正。先讲一下怎么做，后面会附上我基于R语言呕心沥血写出来的代码（其实也就是十几行，毕竟是个生信小白，讲究能用就行）。

无论你用bismark，BSMAP还是BS-seeker，在mapping之后，都会得到一个包含有染色体位置，以及该位置支持甲基化的reads数和不支持甲基化的reads数，如下（bismark后得到的cov文件）：

前三列是位置，第四列是该位点甲基化率，第五列是支持甲基化的reads数，第六列是不支持甲基化的reads数。这说的啥玩意？别着急，给你拿4023位置举个例子，五六两列为3和1，意思就是这个位点在测序时（深度30X，最多可以被检测到30次）被检测到了4次，其中3次发现该位点有甲基化，所以是75%的甲基化率。懂了吧，那就继续。

由于WGBS并不是百分百的转化率，因此我们需要在计算测序数据转化率的基础上（如何计算转化率？我在文章后面会安利一篇文章），N个reads中被检测到X个甲基化reads是否可靠。需要用二项分布来证明，FDR来校正。如转化率为99.98%，检测到了20个reads，18个reads是甲基化的，怎么做呢？如下：

为什么FDR校正后的数值和二项分布的p值一样？因为只有一个数据。数据多了，自然就不一样了。一般来说，FDR校正后，只保留<0.01的甲基化位点就好了，当然，阙值也可以自定，cell文章1K拟南芥的甲基化位点就选用了FDR<0.001。

附R代码：

M <- read.table("C:/Users/SSD/Desktop/222.txt",header=F,sep="")

N1 <- M$V4 #甲基化列

N2 <- M$V3 #测序深度列

len <- length(N1)

pv <- numeric(len) #设置相同长度的变量

P <- 0.032 #测序深度，需要自己计算

for(i in 1:len){pv[i]=binom.test(N1[i], N2[i], P)$p.value} #计算P值

qv <- p.adjust(pv,"fdr") #对p值进行FDR检验

total =cbind(M,qv) #将FDR(qv)列加入到原文件后面

write.table(total, file = "C:/Users/SSD/Desktop/111.txt", row.names=FALSE, quote=FALSE) #输出合并后的文件

这一步是必须的，因为比对出来的文件中，有一些位点的甲基化reads为0，或者甲基化率非常低，这些都会被二项检验过滤掉。保留了FDR<0.01的位点之后，再过滤掉reads数比较低的位点，一般是<5。可以想象一下，30X数据，某位点总共才检测到一两个reads，可信度也太低了，必须过滤掉。接下来，就快乐的去call DMR吧。对了，植物中，CG，CHG和CHH要分别过滤。

甲基化分析还处于摸索阶段，欢迎各路大神提问题并交流。

甲基化率计算参考简书：https://www.jianshu.com/p/7e4903232b3f