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 <script type="math/tex" id="MathJax-Element-1">i</script>是那一个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,可能会带来错误的信号。

已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页