bioconductor的minfi处理甲基化数据

学习最快乐

下载工具包：

[AppleScript] 纯文本查看 复制代码

source(“[url=http://www.bioconductor.org/biocLite.R]http://www.bioconductor.org/biocLite.R[/url]“)
biocLite(c(“minfi”, “GEOquery”))#如果想同时下载多个包

载入工具包:

[AppleScript] 纯文本查看 复制代码

library(minfi)
library(GEOquery)

Overview

Illumina 450k DNA甲基化微阵列是用于测定DNA甲基化的廉价且相对全面的阵列。 它是DNA甲基化大样本分析的首选平台，特别是所谓的EWAS（表观基因组范围的关联研究）。这些研究像全基因组相关研究而不是具有基因型（通常在SNP阵列或外显子/全基因组测序上测量）的与表观型相关联表型研究。

DNA methylation

含有很多CpG 结构，2CpG 和2GPC 中两个胞嘧啶的5 位碳原子通常被甲基化，且两个甲基集团在DNA 双链大沟中呈特定三维结构。基因组中60%～ 90% 的CpG 都被甲基化，未甲基化的CpG 成簇地组成CpG 岛，位于结构基因启动子的核心序列和转录起始点。有实验证明超甲基化阻遏转录的进行。DNA 甲基化可引起基因组中相应区域染色质结构变化，使DNA 失去核酶ö限制性内切酶的切割位点，以及DNA 酶的敏感位点，使染色质高度螺旋化，凝缩成团，失去转录活性。
5 位C 甲基化的胞嘧啶脱氨基生成胸腺嘧啶（C-T转换），由此可能导致基因置换突变，发生碱基错配，如果在细胞分裂过程中不被纠正，就会诱发遗传病或癌症。

• DNA 甲基化芯片中 的甲基化水平值则代表了甲基化等位基因密度(M) 基化等位基因密度(U)的比值, 常称作beta 可见,beta 之间的连续值.
  

DNA甲基化的大多数分析的目标是将表型变化与一个或多个基因座中甲基化的变化相关联。(有问题)

Array Design

450k阵列具有非常不寻常的设计，在一定程度上影响分析。 它真的是一个双色阵列和两个单色阵列的混合。 探头有两种主要类型（I型和II型），探针设计会影响探头的信号分布。

450k阵列的原始数据格式称为IDAT。 因为数组是用两种不同的颜色来测量的，所以每个样本都有两个文件，通常是扩展名_Grn.idat和_Red.idat。 Illumina用于分析该阵列的软件套件称为GenomeStudio。 从业者只能访问GenomeStudio的处理数据而不是原始的IDAT文件，但是经验已经表明，在IDAT文件中有信息有利于分析。

Data

我们将访问一个旨在研究急性躁狂症的数据集。 从患有急性躁狂症住院的个体以及未受影响的对照组获得血清样品。

我们要获取可用的IDAT文件。 首先下载数据

[AppleScript] 纯文本查看 复制代码

library(GEOquery)
getGEOSuppFiles(“GSE68777”)
untar(“GSE68777/GSE68777_RAW.tar”, exdir = “GSE68777/idat”)#解压
head(list.files(“GSE68777/idat”, pattern = “idat”))

目前的R不支持压缩的IDAT files。所以用如下方法解决

[AppleScript] 纯文本查看 复制代码

idatFiles <- list.files(“GSE68777/idat”, pattern = “idat.gz$”, full = TRUE)
sapply(idatFiles, gunzip, overwrite = TRUE)

现在我们使用read.450k.exp（）/现在是read.metharray.exp读取IDAT文件，这个文件读取目录中的所有IDAT文件。

[AppleScript] 纯文本查看 复制代码

rgSet <- read.metharray.exp(“GSE68777/idat”)
rgSet
pData(rgSet)
head(sampleNames(rgSet))

现在我们有数据，但请注意，我们没有pheno数据。 文件名在这里非常无益。 这些名称由GEO标识符（“GSM”部分）组成，其后跟标准IDAT命名约定，其中10位数字是数组标识符，后跟“R01C01”形式的标识符。 这是因为每个阵列实际上允许以6x2排列的12个样品的杂交。 “5958091020_R01C0”表示芯片“5958091020”上的第1行和第1列。 这一切都很好，但不能帮助我们了解哪些样本是案例，哪些是控件。

我们现在获得标准的GEO表示，以获得存储在GEO中的表型数据。 该表型数据中的大多列是无关紧要的（包含数据，例如提交数据的人的地址）; 在此只保留的信息。

[AppleScript] 纯文本查看 复制代码

geoMat <- getGEO(“GSE68777”,getGPL=FALSE)#这里我多设置了一个参数，getGPL=FALSE，就不会出错了。
pD.all <- pData(geoMat[[1]])
pD <- pD.all[, c(“title”, “geo_accession”, “characteristics_ch1.1”, “characteristics_ch1.2”)]
head(pD)
names(pD)[c(3,4)] <- c(“group”, “sex”)
pD[img]https://chart.googleapis.com/chart?cht=tx&chl=group%20%3C-%20sub(%22%5Ediagnosis%3A%20%22%2C%20%22%22%2C%20pD[/img]group)
pD[img]https://chart.googleapis.com/chart?cht=tx&chl=sex%20%3C-%20sub(%22%5ESex%3A%20%22%2C%20%22%22%2C%20pD[/img]sex)

我们现在需要将该pheno数据合并到甲基化数据中。 为此，我们需要一个常见的样本标识符，我们确保我们以与甲基化数据相同的顺序重新排序表型数据。 最后我们将表型数据放入甲基化数据。

[AppleScript] 纯文本查看 复制代码

sampleNames(rgSet) <- sub(“.*_5”, “5”, sampleNames(rgSet))
rownames(pD) <- pD$title
pD <- pD[sampleNames(rgSet),]
pData(rgSet) <- pD
rgSet

前处理

rgSet对象是一个叫做RGChannelSet的类，它表示两个颜色数据，一个绿色和一个红色通道，非常类似于’ExpressionSet’。

第一步通常是使用多个函数来预处理数据

• preprocessRaw（）：什么都不做。
• preprocessIllumina（）：使用Illumina的标准处理选择。
• “preprocessQuantile（）”：使用适合甲基化数组的分位数归一化版本。
• preprocessNoob（）：使用NOOB背景校正方法。
• preprocessSWAN（）：使用SWAN方法。
• preprocessFunnorm（）：使用函数规范化。
  

这些函数输出不同类型的对象。

minfi中的类层次结构如下：数据可以存储在甲基化和非甲基化通道中或甲基化百分比（称为Beta）通道中。对于第一种情况，有’MethylSet`类，对于第二种情况，有’RatioSet’类。当（有甲基化/非甲基化值）时，仍然可以即时计算Beta值。可以将“MethylSet”转换为“RatioSet”，并使用“compareConvert（）”。

除了这两个类之外，还有“GenomicMethylSet”和“GenomicRatioSet”。 “Genomic”表示数据已经使用mapToGenome（）函数与基因组坐标相关联。

大多数分析的起点应该是一个“GenomicRatioSet”类。如果您选择的预处理方法没有得到您的使用，请使用ratioConvert（）和mapToGenome（）去到最后一步。

运行preprocessQuantile（），它到达GenomicRatioSet：

[AppleScript] 纯文本查看 复制代码

grSet <- preprocessQuantile(rgSet)#前处理
grSet

这个类似 SummarizedExperiment; 我们可以通过下面的函数得到 CpGs 的位置

[AppleScript] 纯文本查看 复制代码

granges(grSet)

The usual methylation measure is called “Beta” values; equal to percent methylation and defined as Meth divided by Meth + Unmeth.
通常的甲基化衡量方法被称为“Beta”值; 等于甲基化百分比，并定义为“Meth”除以“Meth + Unmeth”。

[AppleScript] 纯文本查看 复制代码

getBeta(grSet)[1:3,1:3]

CpG形成的簇称为“CpG岛”的簇。 CpG群岛附近的地区被称为CpGShores，其次是CpG Shelfs，最后是CpG Open Sea探测器。 一个简单的方法如下

[AppleScript] 纯文本查看 复制代码

head(getIslandStatus(grSet))

差异甲基化

一旦数据被标准化，一种可能的方法是基于Beta值上使用“r Biocpkg（”limma“）来识别差异甲基化的CpG。
另一种可能的方法是找到所有改变方向相同的CpG簇。 通过r Biocpkg("bumphunter")包的bumphunter()函数可以实现

[size=0em]​

大表哥

写的很好，受教了