DNA Methylation中的Detect P value解析

继续研究IDAT文件的读取问题，目前遇到的问题的是如何提取detectP value. Detect P value是测序过后，每一个芯片上的数据的可信度，也就是说，芯片上有那么多CpG，有时候可能会因为人为的失误，导致有些数据出现错误，在那种时候，我们就需要先评测一下每一个CpG的可信度有多高，如果可信度太低，也就是说，这个位点很有可能出错了，我们就把这个位点直接删除，以免它影响了我们的分析。

具体来说，应该如何判断一个CpG是不是“测序质量差”呢？在illumina的DNA Methylation芯片中，有一系列的Control Probe，就是说对照位点，那些对照位点的作用就是让你用它们来比较，测出来的CpG位点是不是正确的。

比如说，以450K的Annotation为例，在Annotation中，包含了850个Probe，这些Probe的作用各异，不过主要的作用，就是让人们与之做出比较。

比如上边的截图，第一列是各种不同功能的对照，有的是用来看Extension的，有的是用来看Bisulfile Conversion的，很多我也不知道是做什么的，太复杂了，改天我认真看看。

Anyway，有关于CpG Probe的Detect Pvalue用到的是第一列中名称是NEGATIVE的位点，这些位点一共有613个

> sum(Anno$ControlProbe[,1]=="NEGATIVE")[1] 613

下面是很重要的个部分：

在芯片中，交叉着Type-II和Type-I两种芯片，不同的芯片对应的红绿数据的值是不一样的。每一个CpG，都会有一个测序的intensity（其实就是芯片上一共测到了CpG），这个数值就是Meth和UnMeth的加和。所以每一个CpG的Detect P value，其实就是这一个CpG所测的颜色通道Intensity，与这个颜色的NEGATIVE Probe的Intensity分布的比较值。

距离来说，假设有一个CpG Probe是TypeII的。对于所有Type-II的CpG的，它们的Meth值是红色通道上的TypeII$AddressA值，他们的UnMeth是绿色通道上的TypeII$AddressA。所以这个位点的Intensity就是：

intensity <- r[TypeII$AddressA, i] + g[TypeII$AddressA, i]

假设i是那一个CpG。

然后我们去找NEGATIVE Probe，一共有613个，直接从Annotation里可以找出来。针对Type-II的特性（它的数值用到了两个颜色通道），我们计算出绿色通道里的NEGATIVE Control Probe的平均值和标准差，再计算出红色通道里的NEGATIVE Control Probe的平均值和标准差。之所以要计算这两个平均值和标准差，是因为我们需要它来生成对应的分布：

detP[TypeII$Name, i] <- 1-pnorm(intensity, mean=rMu[i]+gMu[i], sd=rSd[i]+gSd[i])

在上边的公式中，pnorm生成了给予rMu和gMu的均值和SD的分布（并且Assume其是一个正态分布），然后将Intensity的数值与这个分布做比较，看这个数值出现的比例是多少。

讲完上边的流程，这个问题就很直观了：NEGATIVE Control Probe是一系列数值会比较低的位点，通过获得它们的均值和标准差，我们就可以知道一个测序质量不及格的分布，那是一个正态分布，然后我们用我们需要的CpG的测序质量（Intensity）去与这个分布作比较，如果它处在这个分布的中间位置，或者说不够显著偏高的位置，我们就可以假设，这个位点是不及格的。

上边写的是Type-II Probe的问题，如果是Type-I，还分为两种，一种是Type-I红色通道，另一种是Type-II绿色通道。针对红色通道，intensity是：

 intensity <- r[TypeI.Red$AddressA, i] + r[TypeI.Red$AddressB, i]

而针对绿色通道，intensity是：

intensity <- g[TypeI.Green$AddressA, i] + g[TypeI.Green$AddressB, i]

可见，他们都只用到了单一通道，所以在计算它们的显著性（Detect P value）的时候，我们只需要单一的通道的分布就行了：

detP[TypeI.Red$Name, i] <- 1-pnorm(intensity, mean=rMu[i]*2, sd=rSd[i]*2)detP[TypeI.Green$Name, i] <- 1-pnorm(intensity, mean=gMu[i]*2, sd=gSd[i]*2)

这样我们就可以获得所有的Probe的Detect P value。这个数值绝对不是什么可有可无的数值，再做任何的分析之前，必须要认真检查，每一个Probe的这个数值，确保无误才能继续，因为质量不好的Probe，可能会带来错误的信号。