1简介

的missMethyl包包含分析Illumina公司HumanMethylation450和MethylationEPIC串珠芯片甲基化数据的功能。这些阵列是全基因组亚硫酸氢盐测序的成本效益替代方案，因此被广泛用于DNA甲基化分析。具体地说,missMethyl包含执行SWAN归一化的函数(Maksimovic, Gordon, and Oshlack 2012)，进行差异甲基化分析RUVm(Maksimovic et al. 2015)，差异变率分析(Phipson and Oshlack 2014)基因集分析(Phipson, Maksimovic, and Oshlack 2016)．随着我们实验室对这些阵列的专门分析的研究继续进行，我们预计包将不断更新新的功能。

原始数据文件是IDAT格式的，可以使用minfi包(Aryee et al. 2014)．统计分析通常对m值进行β\ (\ \)值用于可视化，两者都可以从对象中提取，这是由创建的一类对象minfi．为了检测差异变量cpg，我们建议对m值进行分析。这里描述的所有分析都是在CpG站点级别上执行的。

2将数据读入R

我们将使用的数据minfiData包中演示的函数missMethyl．示例数据集在两张幻灯片上有6个样本。示例信息在targets文件中。目标文件中的一个基本列是Basename列，该列表示要读入的idat文件的位置。用于读入数据的R命令取自minfi用户指南。有关如何将IDAT文件读入R的其他详细信息，以及有关质量控制的信息，请参阅minfi用户指南。

library(missMethyl) library(limma) library(minfi)

baseDir <- system. library(minfiData)file("extdata"， package = "minfiData") target <- read.metharray.sheet(baseDir)

## [1] "/home/biocbuild/bb -3.6-bio /R/library/minfiData/extdata/SampleSheet.csv"

目标[1:9]

# # Sample_Name Sample_Well Sample_Plate Sample_Group Pool_ID人年龄性别# # 1 GroupA_3 H5 < NA > GroupA < NA > id3 83 # # 2 GroupA_2 D5 < NA > GroupA F < NA > id2 58 # # 3 GroupB_3 C6 < NA > GroupB < NA > id3 83 # # 4 GroupB_1 F7 < NA > GroupB < NA > id1 75 F # # 5 GroupA_1七国集团(G7) < NA > GroupA < NA > id1 75 F # # 6 GroupB_2 H7 < NA > GroupB F < NA > id2 58 # # # # 1 # # 2 #正常状态# 3癌症# # 4 # 5 # #正常6癌症

目标(,十12)

##阵列幻灯片## 1 R02C02 5723646052 ## 2 R04C01 5723646052 ## 3 R05C02 5723646053 ## 5 R05C02 5723646053 ## 6 R06C02 5723646053 ## baseename ## 1 /home/biocbuild/ bbbs -3.6-bioc/R/library/minfiData/extdata/5723646052/5723646052 _r02c02 ## 2/ home/biocbuild/ bbbs -3.6-bioc/R/library/minfiData/extdata/5723646052/5723646052 _r04c01 ## 3 /home/biocbuild/ bbbs -3.6-bioc/R/library/minfiData/extdata/5723646052/5723646052 _r05c02 ## 45 /home/biocbuild/论坛-3.6-bioc/R/library/minfiData/extdata/5723646053/5723646053_R04C02 ## 6 /home/biocbuild/论坛-3.6-bioc/R/library/minfiData/extdata/5723646053/5723646053_R05C02 ##

rgSet <- read.metharray。Exp(目标=目标)

数据现在是RGChannelSet对象，需要规范化并转换为MethylSet对象。

3.数组归一化中的子分位数(SWAN)

SWAN(数组内子集-分位数归一化)是Illumina 450k和EPIC串珠芯片的数组内归一化方法。已证明Infinium I和Infinium II在单个Illumina人类甲基化阵列上的分析方法存在技术差异(Bibikova et al. 2011，Dedeurwaerder, Defrance, and Calonne (2011)）．使用SWAN方法在保持重要生物学差异的同时，大大减少了分析设计之间的技术变异性。SWAN方法假设50bp探针序列中cpg的数量反映了所调查区域的潜在生物学。因此，无论检测类型如何，探针体中具有相同数量cpg的探针的总体强度分布应该是相同的。然后，该方法使用子集分位数归一化方法来调整每个数组的强度(Maksimovic, Gordon, and Oshlack 2012)．

天鹅可以MethylSet，RGChannelSet或MethyLumiSet作为输入。需要注意的是，为了创建规范化子集，天鹅随机选择有1、2和3个cpg的Infinium I和II探针;因此，我们建议使用set.seed之前，以确保归一化强度将是相同的，如果是重复归一化。

Infinium I和II分析设计之间的技术差异可能导致beta值分布异常(图1，面板“原始”)。使用SWAN纠正了Infinium I和II分析设计之间的技术差异，并产生了更平滑的整体β\ (\ \)值分布(图1，面板“SWAN”)。

mSet <- preprocessRaw(rgSet)

mSetSw <- SWAN(mSet,verbose=TRUE)

# # # # 450 k(天鹅)准备正常化子集# #(天鹅)正常化甲基化频道# #(天鹅)正常化数组1 6 # #(天鹅)正常化数组2的6 # #(天鹅)正常化数组3的6 # #(天鹅)正常化数组4 6 # #(天鹅)正常化阵列5 6 # #(天鹅)正常化数组6 6 # #(天鹅)正常化unmethylated频道# #(天鹅)正常化数组1 6 # #(天鹅)正常化数组2的6 # #(天鹅)正常化数组3的6 # #(天鹅)正常化数组4 6 # #(天鹅)正规范化数组6 ## [SWAN]正规范化数组6 (6)

par(mfrow=c(1,2)， cex=1.25) densityByProbeType(mSet[，1]， main = "Raw")

图1:密度分布\ (埃塔\)使用SWAN前后的值

4过滤掉质量差的探针

根据检测p值可以过滤出质量差的探针。对于本例，为了保留CpG用于进一步分析，我们要求所有样品的检测p值小于0.01。

detP <- detectionP(rgSet) keep <- rowsum (detP < 0.01) == ncol(rgSet) mSetSw <- mSetSw[keep，]

5提取Beta和m值

现在数据已经天鹅正常化，我们可以提取β\ (\ \)和来自对象的m值。在计算m值时，我们倾向于在甲基化和非甲基化强度上增加一个偏移量，因此我们分别提取甲基化和非甲基化通道并执行我们自己的计算。对于所有后续分析，我们使用20000 cpg的随机选择来减少计算时间。

mset_reduced <- mSetSw[sample(1:nrow(mSetSw)， 20000)，] meth <- getMeth(mset_reduced) unmeth <- getUnmeth(mset_reduced) Mval <- log2((meth + 100)/(unmeth + 100)) beta <- getBeta(mset_reduced) dim(Mval)

## [1] 20000

par(mfrow=c(1,1)) plotMDS(Mval, labels=目标$Sample_Name, col=as.integer(factor(目标$status))) legend("topleft"，legend=c("癌症"，"正常")，pch=16,cex=1.2,col=1:2)

图2:癌症和正常样本的多维标度(MDS)图

MDS图(图2)是一个很好的完整性检查，以确保样本根据主要感兴趣的因素聚集在一起，在本例中是癌症和正常。

6测试差异甲基化使用

为了测试差异甲基化，我们使用limma包(G. K.史密斯2005)，它采用了基于高斯模型理论的经验贝叶斯框架。首先，我们需要建立设计矩阵。有许多方法可以做到这一点，最直接的是直接从目标文件。有很多变量状态列显示癌症/正常样本。从人列的目标文件，我们看到癌症/正常样本是匹配的，3个人每个贡献两个a癌症而且正常的样本。自limma模型框架可以处理任何可以通过设计矩阵总结的实验设计，我们可以在分析中考虑数据的配对性质。如需更复杂的实验设计，请参阅limma用户指南。

组<- factor(目标$status，水平=c("正常"，"癌症"))id <- factor(目标$person)设计<-模型。矩阵(~id +组)设计

##(拦截)idid2 idid3 groupcancer ## 1 1 0 1 0 0 # 2 1 1 0 1 1 # 4 1 0 0 1 # 5 1 0 0 0 0 0 # 6 1 1 0 1 ## attr(，“分配”)## [1]0 1 1 2 ## attr(，“对比”)## attr(，“对比”)$id ##[1]“对照。治疗”## ## attr(，“对照”)$group ##[1]“对照。治疗”

现在我们可以测试差异甲基化使用lmFit而且易趣函数从limma．我们使用m值矩阵作为输入数据。

健康。减小<- lmFit(Mval，设计)适合度。reduced <- eBayes(fit.reduced)

的值可以显示超甲基化(1)和低甲基化(-1)的数目decideTests函数limma和前10个差异甲基化CpGs癌症与正常的提取使用topTable．

总结(decideTests (fit.reduced))

##(拦截)idid2 idid3 groupcancer ## 1 7064 0 97 665 ## 0 3359 20000 19898 18809 ## 1 9577 0 5 526

前< -topTable (fit.reduced系数= 4)

## logFC AveExpr t P.Value adj.P.Val B ## cg21938148 4.482621 -0.4749969 15.89282 1.15771e -05 0.03141079 3.836719 ## cg13272280 4.263743 -2.1495705 15.21106 1.453251e-05 0.03141079 3.435331 ## cg10471437 4.872744 - 1.1454806 14.26313 2.027309e-05 0.03141079 3.424614 ## cg23664459 3.539586 -2.0681570 14.22077 2.058703e-05 0.03141079 3.423856 ##Cg26532358 3.587568 -0.8751333 13.68149 3.272654 ## cg00995327 5.068487 -0.7270585 13.47634 2.716959e-05 0.03141079 3.212541 ## cg01134185 3.692594 -0.7147741 13.34658 2.855829e-05 0.03141079 3.173752 ## cg00262031 3.773122 -2.0652330 13.15865 3.072153e-05 0.03141079 3.116489

请注意，由于我们对m值进行了分析，因此logFC而且AveExpr列按m值尺度计算。对于可解释性和可视化，我们可以看看β\ (\ \)值。图3显示了前4种不同甲基化cpg的beta值。

cpgs <- rownames(top) par(mfrow=c(2,2)) for(i in 1:4){stripchart(beta[rownames(beta)==cpgs[i]，]~design[，4]，method="jitter"， group.names=c("Normal"，"Cancer")，pch=16,cex=1.5,col=c(4,2)，ylab=" beta values"， vertical=TRUE,cex.axis=1.5,cex.lab=1.5) title(cpgs[i]，cex.main=1.5)}

图3:的\ (埃塔\)前4个差异甲基化CpGs的值

7在测试差异甲基化时去除不需要的变异

与其他平台一样，450k阵列研究也会受到不必要的技术变化的影响，例如批量效应和其他通常未知的变化来源。在基因表达阵列研究中，不需要的变异的不利影响已被广泛记录，并已被证明既能降低检测真实差异的能力，又能增加错误发现的数量。因此，当数据明显受到不需要的变化的显著影响时，建议进行调整以减轻其影响。

missMethyl提供了一个limma类似接口的函数从CRAN包ruv这样就可以在执行微分分析时去除不必要的变化(Maksimovic et al. 2015)．所有的方法都依赖于负控制特征来准确估计不需要的变异成分。负控制特征是探针/基因等。都是已知的先天的不是真正与感兴趣的生物因素相关，而是受不需要的变异影响。例如，在微阵列基因表达研究中，这些可能是管家基因或一组峰值控制。Gagnon-Bartsch和Speed广泛讨论了负控制特性(2012)和Gagnon-Bartsch等人。(2013)．一旦从数据中准确估计出不需要的因素，就会在描述差分分析的线性模型中对它们进行调整。

如果不知道负控制特性先天的，它们可以根据经验加以识别。这可以通过两阶段的方法来实现，RUVm，基于RUV-inverse．阶段1涉及使用RUV-inverse613个Illumina阴性对照(INCs)为阴性对照特征。这将产生一个cpg列表，根据它们与感兴趣因素的关联程度，根据p值进行排名。然后，该列表可用于识别一组经验控制探针(ecp)，这将捕获比单独使用INCs更多的不需要的变化。通过指定与感兴趣因素关联最少的cpg的比例作为负控制特征来选择ecp;这可以基于FDR截止或从排名列表底部取固定百分比的探测来完成。阶段2涉及对原始数据执行第二差分甲基化分析RUV-inverse和ecp。的配对性质为简单起见，我们忽略了癌症而且正常的本例中的示例。

# get所有探针的M值meth <- getMeth(mSet) unmeth <- getUnmeth(mSet) M <- log2((meth + 100)/(unmeth + 100)) grp <- factor(targets$status,levels=c("normal"，"cancer")) des <- model.matrix(~grp) des

##(拦截)grpcancer ## 1 1 0 ## 2 1 0 ## 3 1 1 ## 4 1 1 ## 5 1 0 ## 6 1 1 ## attr(，“分配”)## [1]0 1 ## attr(，“对比”)## attr(，“对比”)$grp ##[1]“对照。治疗”

INCs <- getINCs(rgSet) head(INCs)

## 13792480 -0.3299654 - 1.0955488 -1.4943396 -1.0067050 ## 34772371 -1.1286422 -0.2995603 -0.8192636 ## 28715352 -0.5553373 -0.7599489 -0.7186973 ## # 33730459 -0.7714684 -0.5622424 -0.7724825 ## # 5723646053 -1.586679 - 0.9217393 ## #-1.161155 -0.6186795 -1.7903619 -1.348105 -1.0067259 ## 74737439 -0.8872082 -1.064986 -0.9841833 ## 33730459 -1.5623138 -2.079184 -1.0445246

Mc <- rbind(M,INCs) ctl <- rownames(Mc) %in% rownames(INCs) table(ctl)

## ctl ## FALSE TRUE ## 485512 613

rfit1 <- RUVfit(data=Mc, design=des, coef=2, ctl=ctl) #第一阶段分析rfit2 <- RUVadj(rfit1)

现在我们已经进行了初始差异甲基化分析，对cpg与感兴趣因子的相关性进行了排序，我们可以根据fdr调整的p值将与感兴趣因子相关性最小的cpg指定为ecp。

top1 <- topRUV(rfit2, num=Inf)

##参数t p. bh .ebayes ## cg04743961 4.838190 26.74467 3.812882e-05 0.140103 1.608653e-04 0.1401969 3.583107e-07 ## cg20925841 4.790211 4.69524 3.837354e-05 0.1401969 4.721205e-07 ## cg10607359 4.394397 34.74068 4.721205e-07 ## cg07655636 4.571758 22.99708 6.424216e-05 0.1401969 6.080091e-07 ## p.ebayes. bh ## cg04743961 0.01017357 ## cg071553360.01017357 ## cg20925841 0.01017357 ## cg03607359 0.01017357 ## cg07655636 0.01017357

ctl <- rownames(M) %in% rownames(top1[top1$p.ebayes. txt]BH > 0.5，])表(ctl)

## ctl ## FALSE TRUE ## 172540 312972

然后我们可以使用ECPs来执行第二次差异甲基化RUV-inverse，根据从数据中估计出的不需要的变化进行调整。

#执行RUV调整和拟合rfit1 <- RUVfit(data=M, design=des, coef=2, ctl=ctl) #第二阶段分析rfit2 <- RUVadj(rfit1) #查看topRUV(rfit2)的顶级结果表

##系数t p. p. bh .ebayes ## cg07155336 5.769286 15.345069 0.002005546 0.3431163 1.43483434e -55 ## cg06463958 5.733093 15.434797 0.001978272 0.3431163 6.749298e- 53 ## cg02040433 5.651399 10.054445 0.005389436 0.3431163 2.146210e-53 ## cg02040433 5.651399 10.054445 0.005389436 0.3431163 2.146210e-53 ##3.710480e- 45 ## cg16306898 5.466780 5.573935 0.0203431163 7.700032e-49 ## p.ebayes.BH ## cg07155336 6.966293e-50 ## cg06463958 1.638433e-49 ## cg02040433 2.135266e-48 ## cg02040433 2.022589e-47 ## cg00306898 2.399554e- 46 ## cg03735888 3.738458e-44

注意，目前不支持对比，因此使用带有截距项的设计矩阵一次只能查询一个感兴趣的因素。

8差异变异性测试(DiffVar)

8.1甲基化数据

而不是测试平均甲基化的差异，我们可能对测试组方差之间的差异感兴趣。例如，有人假设癌症中高度可变的cpg对肿瘤进展很重要(K. D. Hansen et al. 2011)．因此，我们可能对在正常样本中一致甲基化，但在癌症样本中甲基化变化的CpG位点感兴趣。

一般来说，我们建议每组至少有10个样本来进行准确的方差估计，但为了本插图的目的，我们对3 vs 3进行分析。在这个例子中，我们感兴趣的是测试癌症组与正常组的差异变异性。注意，当我们指定系数参数，它对应于用于测试差异可变性的设计矩阵的列，我们需要指定截距和第四列。ID变量是一个讨厌的参数，在获得绝对偏差时不使用，但是它可以包括在线性建模步骤中。对于甲基化数据，该函数将取m值矩阵，β\ (\ \)值或对象作为输入。如果β\ (\ \)提供值后，执行logit转换。注意，默认情况下，varFit方法中的健壮设置limma框架，这需要使用statmod包中。

fitvar <- varFit(Mval, design = design, coef = c(1,4))

高变基因(1)和低变基因(-1)的数量癌症vs正常的可使用decideTests．

总结(decideTests (fitvar))

##(拦截)idid2 idid3 groupcancer ## 1 0 1 1 0 ## 0 19746 19998 19982 19999 ## 1 254 1 17 1

topDV <- topVar(fitvar, coef=4

## SampleVar LogVarRatio DiffLevene P.Value ## cg17900854 5.930237 4.881015 3.006691 5.174173 2.476912e -06 ## cg25587181 4.374441 5.370599 2.812563 4.390434 1.134141e-05 ## cg20214319 4.546082 6.251434 2.842861 4.353243 1.344707e-05 ## cg00231519 8.016708 -1.198117 -1.601571 -4.307668 1.653774e-05 ## cg17631972 5.893932 3.902126 2.599803 4.449899e -05 #### cg07035503 ## cg07035506 ## cg26029345 0.066150946 ## cg20214319 0.066150946 ## # cg17434540 0.070082638 ## cg17631972 0.070082638 ## cg12434587 0.072323608 ## cg07035503 0.111782841 ## Adj.P.Value ##

也可以使用不包括截距项的设计矩阵的另一种参数化，并与测试的具体对比contrasts.varFit．这里我们指定设计矩阵，使前两列对应于正常的而且癌症组,分别。

Design2 <- model.matrix(~0+group+id)contr <- varFit(Mval, design=design2, coef=c(1,2)) contr <- makeContrasts(groupcancer-groupnormal,levels=colnames(design2))contr <- constants . varfit (fitvar.contr，反差=contr)

结果与之前相同。

总结(decideTests (fitvar.contr))

## groupcancer - groupnormal ## -1 0 ## 0 19999 ## 1

托普瓦(fitvar.contr系数= 1)

## SampleVar LogVarRatio DiffLevene P.Value ## cg17900854 5.930237 4.881015 3.006691 5.174173 2.476912e -06 ## cg25587181 4.374441 5.370599 2.812563 4.390434 1.134141e-05 ## cg20214319 4.546082 6.251434 2.842861 4.353243 1.344707e-05 ## cg00231519 8.016708 -1.198117 -1.601571 -4.307668 1.653774e-05 ## cg17631972 5.893932 3.902126 2.599803 4.449899e -05 #### cg07035503 ## cg07035506 ## cg26029345 0.066150946 ## cg20214319 0.066150946 ## # cg17434540 0.070082638 ## cg17631972 0.070082638 ## cg12434587 0.072323608 ## cg07035503 0.111782841 ## Adj.P.Value ##

的β\ (\ \)前4个差分变量cpg的值可以在图4中看到。

cpgsDV <- rownames(topDV) par(mfrow=c(2,2)) for(i in 1:4){stripchart(beta[rownames(beta)==cpgsDV[i]，]~design[，4]，method="jitter"， group.names=c("Normal"，"Cancer")，pch=16,cex=1.5,col=c(4,2)，ylab=" beta values"， vertical=TRUE,cex.axis=1.5,cex.lab=1.5) title(cpgsDV[i]，cex.main=1.5)}

图4:的\ (埃塔\)前四个差分变量cpg的值

8.2RNA-Seq表达数据

如果技术是基因表达微阵列，那么对表达数据差异可变性的测试就很简单。表达式值的矩阵可以直接提供给varFit函数。对于RNA-Seq数据，需要考虑计数数据中出现的均值-方差关系。为了处理这个问题，我们应用了轰转换(Law et al. 2014)获取观测权重，然后在线性建模步骤中使用。对于RNA-Seq数据，使用varFit函数的值为DGElist对象作为输入。

为了证明这一点，我们使用tweeDEseqCountData包中。这些数据是国际HapMap项目的一部分，包括来自69个不相关的尼日利亚人的RNA-Seq配置文件(Pickrell et al. 2010)．唯一的协变量是性别，所以我们可以看看男性和女性之间的差异变量表达。我们遵循的代码limma在测试差异可变性之前读入和处理数据。

首先，我们加载数据并提取相关信息。

library(tweeDEseqCountData) data(pickrell1) counts<-exprs(pickrell1.eset) dim(counts)

## [1] 38415 69

性别<- pickrell1.eset$性别表(性别)

##性别##女男性## 40 29

rm(pickrell1.eset) data(genderGenes) data(annotEnsembl63) annot <- annotEnsembl63[，c("Symbol"，"Chr")] rm(annotEnsembl63)

我们现在有了每个Ensemble基因的计数、性别和注释(基因符号和染色体)。我们可以形成一个DGElist对象使用刨边机包中。

library(edgeR) y <- DGEList(counts=counts, genes=annot[rownames(counts)，])

我们通过在至少20个样本中保持每百万次读取至少1个计数的基因以及已经定义注释的基因来过滤低表达基因。最后，我们进行缩放归一化。

isexpr <- rowsum (cpm(y)>1) >= 20 hasannot <- rowsum (is.na(y$genes))==0 y <- y[isexpr & hasannot，，keep.lib. txt]。大小= FALSE)暗(y)

## [1] 17310 69

y <- calcNormFactors(y)

我们建立了设计矩阵并测试了差异变异性。在这种情况下，没有妨害参数，所以系数不需要显式指定。

设计。Hapmap <- model.matrix(~gender)hapmap <- varFit(y, design = design.hapmap)

##使用voom将计数转换为每百万日志计数。

fitvar。Hapmap $genes <- y$genes

我们可以显示测试的结果:

总结(decideTests (fitvar.hapmap))

##(拦截)性别男性## -1 0 2 ## 0 0 17308 ## 1 17310 0

topDV。hapmap <- topVar(fitvar.hapmap,coef=ncol(design.hapmap)) topDV.hapmap

##符号Chr SampleVar LogVarRatio DiffLevene t ## ENSG00000213318 RP11-331F4.1 16 5.69839463 -2.562939 -0.9859943 -8.031243 ## ENSG00000129824 RPS4Y1 Y 2.32497726 -2.087025 - 0.4585233 -4.957005 ## ENSG00000233864 TTTY15 Y 6.790044140 -2.245369 -0.6085233 -4.612934 ## ENSG00000176171 BNIP3 10 0.41317384 1.199292 0.3632133 4.219404 ## ENSG00000197358 BNIP3P1 14 0.39969125 1.149754 0.3353288 4.058147 ## ENSG00000025039 RRAGD 6 0.91837213 1.091229 0.4926839 3.977022 ## ENSG00000103671 TRIP4 150.07456448 -1.457139 -0.1520583 - 3.89131300 ## ENSG00000171100 MTM1 X 0.44049558 -1.133295 -0.3334619 -3.896490 ## ENSG00000149476 DAK 11 0.13289523 -1.470460 -0.1919880 -3.785893 ## ENSG00000064886 CHI3L2 1 2.70234584 1.468059 0.8449434 3.782010 ## P.Value Adj.P.Value ## ENSG00000213318 7.238039e-12 1.252905e-07 ## ensg00000233824 3.960855e-06 3.428120e-02 ## ENSG00000176171 6.441668e-05 3.97398e -01 #### eng00000103671 1.921104e-04 4.375736e-01 ## eng00000171100 2.22293e -04 4.375736e-01 ## eng00000149476 2.956364e-04 4.425050e-01 ## eng00000064886 4.425050e- 04

前4个差异可变基因的每百万日志计数见图5。

genesDV <- rownames(topDV.hapmap) par(mfrow=c(2,2)) for(i in 1:4){stripchart(cpm(y,log=TRUE)[rownames(y)==genesDV[i]，]~design.hapmap[，ncol(design.hapmap)]，method="jitter"， group.names=c("Female"，"Male")，pch=16,cex=1.5,col=c(4,2)，ylab=" log counts per million"， vertical=TRUE,cex.axis=1.5,cex.lab=1.5) title(genesDV[i]，cex.main=1.5)}

图5:前四个差异可变表达基因的对数计数为百万分之一

9基因本体分析

一旦进行了差异甲基化或差异变异性分析，了解哪些基因通路被重要的CpG位点靶向可能是有兴趣的。从文献中还不完全清楚如何最好地进行这样的分析，但是Geeleher等人。(Geeleher et al. 2013)表明用甲基化数据进行基因本体分析时存在严重偏倚。这是因为在几种不同的阵列技术上，每个基因都有不同数量的探针。对于Illumina Infinium HumanMethylation450阵列，每个基因的探针数量从1到1299，每个基因的中位数为15个探针。对于EPIC阵列，范围为1到1487，每个基因的中位数为20个探针。这意味着当将CpG位点映射到基因时，一个基因更有可能被定位选择是否有许多与该基因相关的CpG位点。

考虑这种选择偏差的一种方法是对每个基因的探针数量与被选择的概率之间的关系进行建模。可以通过调整goseq杨等人的方法。(Young et al. 2010)．然后，每个基因都有一个与之相关的先验概率，并且可以执行修改版的超几何测试，测试每个基因集中所选基因的过度代表性。

的gometh函数对GO类或KEGG通路进行基因集检测(Phipson, Maksimovic, and Oshlack 2016)．的gsameth函数是一个更通用的基因集测试函数，它可以将用户指定的基因集列表作为输入。注意，对于gsameth，则基因集中每个基因的基因id的格式必须为Entrez基因id．例如，布罗德的分子特征数据库中的整个策划基因集列表(C2)可以指定为输入。这些列表的R版本可以从http://bioinf.wehi.edu.au/software/MSigDB/index.html．这两个函数都以重要CpG探测名称的向量作为输入。

来说明如何使用gometh，考虑差异甲基化分析的结果RUVm．

topRUV (rfit2)

##系数t p. p. bh .ebayes ## cg07155336 5.769286 15.345069 0.002005546 0.3431163 1.43483434e -55 ## cg06463958 5.733093 15.434797 0.001978272 0.3431163 6.749298e- 53 ## cg02040433 5.651399 10.054445 0.005389436 0.3431163 2.146210e-53 ## cg02040433 5.651399 10.054445 0.005389436 0.3431163 2.146210e-53 ##3.710480e- 45 ## cg16306898 5.466780 5.573935 0.0203431163 7.700032e-49 ## p.ebayes.BH ## cg07155336 6.966293e-50 ## cg06463958 1.638433e-49 ## cg02040433 2.135266e-48 ## cg02040433 2.022589e-47 ## cg00306898 2.399554e- 46 ## cg03735888 3.738458e-44

表(rfit2 p.ebayes美元。Bh < 0.01)

## ##假真## 424168 61344

在1%的错误发现率截点下，仍有数万个CpG位点发生差异甲基化。毫无疑问，这些将映射到基因组中的几乎所有基因，使基因本体分析无关紧要。在这种情况下，选择cpg的一个选项是，不仅应用错误发现率下限，而且应用\(δβ\ \ \)截止。但是，对于这个数据集，取一个比较大的\(δβ\ \ \)0.25的临界值仍然有超过30000 CpGs的差异甲基化。

beta <- getBeta(mSet) beta_norm <- rowMeans(beta[，des[，2]==0]) beta_can <- rowMeans(beta[，des[，2]==1]) Delta_beta <- beta_can - beta_norm sigDM <- rfit2$p.ebayes。BH < 0.01 & abs(Delta_beta) > 0.25表(sigDM)

## sigDM ##错误TRUE ## 451760 33748

相反，我们将前10000个CpG站点作为输入gometh．

topCpGs<- topruv (rfit2,number=10000) sigCpGs <- rownames(topCpGs) sigCpGs[1:10]

##[1]“cg07155336”“cg06463958”“cg00024472”“cg02040433”“cg13355248”##[6]“cg02467990”“cg00817367”“cg11396157”“cg16306898”“cg03735888”

它以CpG名称的字符向量作为输入，还可以选择测试的所有CpG站点的字符向量。如果all.cpg参数省略，则数组上的所有cpg都用作背景。若要更改数组类型，请使用array.type参数可以指定为“450K”或“EPIC”。默认值为“450K”。

如果plot.bias参数是真正的，将显示每个基因被选中的概率与探针数量之间的关系图。

对于GO术语的测试，使用集合参数的值为“GO”，这是默认设置。对于KEGG通路分析，设置集合“KEGG”。这个函数topGO显示了最丰富的GO类别。对于KEGG测试，使用topKEGG函数。的函数goana而且kegga分别称为GO和KEGG通路分析。

在我们的示例数据集上进行GO测试:

库(illuminahumanmethylation450kano .ilmn12.hg19) gst <- gometh(sig。cpg=sigCpGs, all.cpg=rownames(rfit2)， collection="GO")

alias2SymbolTable(flat$symbol)中的##警告:一个##或多个alias2SymbolTable中忽略多个符号

topGO(销售税)

## Term Ont N DE P.DE ## GO:0007399神经系统发育BP 2178 582 9.314801e-92 ## GO:0048731系统发育BP 4377 916 1.132861e-84 ## GO:0048856解剖结构发育BP 5376 1051 4.352461e-82 ## GO:0044707单细胞多细胞生物发育BP 5989 1133 3.253996e-81 ## GO:0007275多细胞生物发育BP 4934 986 6.566254e-81 ## GO:0044767单细胞发育过程BP 5690 1082 1.864022e-77 ## GO:0032502发育过程BP 5777 10922.407186e-76 ## GO:0032501多细胞生物过程BP 6938 1218 1.891492e-65 ## GO:0022008神经发生BP 1430 394 6.475886e-65 ## GO:0048699神经元分化BP 1333 375 2.1621422e -64 ## GO:0030182神经元分化BP 1213 349 1.100265e-62 ## GO:0009653解剖结构形态发生BP 2368 547 1.418750e-59 ## GO:0030154细胞分化BP 3739 748 2.127294e-57 ## GO:0048869细胞发育过程BP 3909 770 3.774791e-56 ## GO:0048468细胞发育BP1875 448 1.303373e-54 ## GO:0007417中枢神经系统发育BP 893 273 2.303354e-54 ## GO: 00048666细胞-细胞信号BP 1541 387 5.728620e-53 ## GO:0048666神经元发育BP 963 275 3.813120e-50 ## GO: 005886质膜CC 4970 902 8.116557e-50 ## GO:0071944细胞外围CC 5067 913 5.574217e-49 ## FDR ## GO:0007399 2.031931e-87 ## GO:0048731 1.235612e-80 ## GO:0048856 3.164820e-78 ## GO:0044707 1.774567e-77 ## GO:0007275 2.864725e-77 ## GO:0044767 6.776964e-74 ##GO:0032502 7.501478e-73 ## GO:0032501 5.157627e-62 ## GO:0022008 1.569611e-61 ## GO:0048699 4.716497e- 59 ## GO:0030182 2.181925e-59 ## GO:0030154 3.569599e-54 ## GO:0048869 5.881664e-53 ## GO:0048869 1.895452e-51 ## GO:0007417 3.140335e-51 ## GO:0007267 7.350830e-50 ## GO:0048666 4.621078e-47 ## GO:0071944 6.079799e-46

对于甲基化数据中基因集测试的更通用版本，用户可以指定要测试的基因集，函数gsameth可以使用。要显示前20条路径，topGSA可称为。gsameth可以是单个基因集，也可以是一系列基因集。基因集中的基因标识符必须是Entrez基因id．为了演示gsameth，下面展示了一个简单的例子，其中的基因集由随机选择的基因组成org.Hs.eg.db包中。

库(org.Hs.eg.db)

##加载所需的包:AnnotationDbi

genes <- toTable(org.Hs.egSYMBOL2EG) set1 <- sample(genes$gene_id,size=80) set2 <- sample(genes$gene_id,size=100) set3 <- sample(genes$gene_id,size=30) genesets <- list(set1,set2,set3) gsa <- gsameth(sig. hs . egsymbol2eg)cpg=sigCpGs, all.cpg=rownames(rfit2)， collection=genesets)

alias2SymbolTable(flat$symbol)中的##警告:一个##或多个alias2SymbolTable中忽略多个符号

topGSA (gsa)

## n de p.de FDR ## [1，] 804 0.1945991 0.3028804 ## [2，] 100 4 0.2019202 0.3028804 ## [3，] 30 0 0.7460663 0.7460663

注意，如果是感兴趣的获取Entrez基因id即重要cpg被映射到的getMappedEntrezIDs可称为。

10会话信息

sessionInfo ()

R版本3.4.2(2017-09-28)平台:x86_64-pc-linux-gnu(64位)运行在Ubuntu 16.04.3 LTS下

矩阵产品:默认BLAS: /home/biocbuild/bbs-3.6-bioc/R/lib/libRblas。so LAPACK: /home/biocbuild/bbs-3.6-bioc/R/lib/libRlapack.so

locale: [1] LC_CTYPE=en_US。utf - 8 LC_NUMERIC = C
[3]而= en_US。utf - 8 LC_COLLATE = C
[5] LC_MONETARY = en_US。utf - 8LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER = en_US。utf - 8 LC_NAME = C
[9] lc_address = c lc_phone = c
[11] LC_MEASUREMENT = en_US。utf - 8 LC_IDENTIFICATION = C

附加的基本包:[1]stats4并行统计图形grDevices utils数据集[8]方法基础

其他附件:[1]org.Hs.eg.db_3.4.2
[2] AnnotationDbi_1.40.0
[3] edgeR_3.20.0
[4] tweeDEseqCountData_1.15.0
[5] minfiData_0.23.0
[6] illuminahumanmethylation450kano .ilmn12.hg19_0.6.0
[8] minfi_1.24.0
[9] bumphunter_1.20.0
[10] locfit_1.5 - 9.1
[11] iterators_1.0.8
[12] foreach_1.4.3
[13] Biostrings_2.46.0
[14] XVector_0.18.0
[15] SummarizedExperiment_1.8.0
[16] DelayedArray_0.4.0
[17] matrixStats_0.52.2
[18] Biobase_2.38.0
[19] GenomicRanges_1.30.0
[20] GenomeInfoDb_1.14.0
[21] IRanges_2.12.0
[22] S4Vectors_0.16.0
[23] BiocGenerics_0.24.0
[24] limma_3.34.0
[25] missMethyl_1.12.0
[26] BiocStyle_2.6.0

通过命名空间加载(且未附加):[1]nlme_1 .1-131
[2] bitops_1.0-6
[3] bit64_0.9-7
[4] httr_1.3.1
[5] RColorBrewer_1.1-2
[6] progress_1.1.2
[7] rprojroot_1.2
[8] tools_3.4.2
[9] backports_1.1.1
[10] doRNG_1.6.6
[11] nor1mix_1.2-3
[12] R6_2.2.2
[13] colorspace_1.3-2
[14] DBI_0.7
[15] methylumi_2.24.0
[16] prettyunits_1.0.2
[17] RMySQL_0.10.13
[18] base64_2.0
[19] bit_1.1-12
[20] compiler_3.4.2
[21] preprocessCore_1.40.0
[22] BiasedUrn_1.07
[23] xml2_1.1.1
[24] pkgmaker_0.22
[25] rtracklayer_1.38.0
[26] bookdown_0.5
[27] readr_1.1.1
[28] genefilter_1.60.0
[29] quadprog_1.5-5
[30] stringr_1.2.0
[31] digest_0.6.12
[32] Rsamtools_1.30.0
[33] illuminaio_0.20.0
[34] rmarkdown_1.6
[35] siggenes_1.52.0
[36] GEOquery_2.46.0
[37] pkgconfig_2.0.1
[38] htmltools_0.3.6
[39] highr_0.6
[40] illuminahumanmethylationepicano .ilm10b2.hg19_0.6.0 [41] ruv_0.9.6
[42] rlang_0.1.2
[43] RSQLite_2.0
[44] bindr_0.1
[45] mclust_5.3
[46] BiocParallel_1.12.0
[47] dplyr_0.7.4
[48] rcurl_1.95 - 4.8
[49] magrittr_1.5
[50] GO.db_3.4.2
[51] GenomeInfoDbData_0.99.1
[52] Matrix_1.2-11
[53] Rcpp_0.12.13
[54] stringi_1.1.5
[55] yaml_2.1.14
[56] MASS_7.3-47
[57] zlibbioc_1.24.0
[58] plyr_1.8.4
[59] grid_3.4.2
[60] blob_1.1.0
[61] lattice_0.20-35
[62] splines_3.4.2
[63] multtest_2.34.0
[64] GenomicFeatures_1.30.0
[65] annotate_1.56.0
[66] hms_0.3
[67] knitr_1.17
[68] beanplot_1.2
[69] rngtools_1.2.4
[70] codetools_0.2-15
[71] biomaRt_2.34.0
[72] xml_3.98 - 1.9
[73] glue_1.2.0
[74] evaluate_0.10.1
[75] data.table_1.10.4-3
[76] openssl_0.9.7
[77] purrr_0.2.4
[78] tidyr_0.7.2
[79] reshape_0.8.7
[80] assertthat_0.2.0
[81] xtable_1.8-2
[82] survival_2.41-3
[83] tibble_1.3.4
[84] GenomicAlignments_1.14.0
[85] registry_0.3
[86] memoise_1.1.0
[87] bindrcpp_0.2
[88] statmod_1.4.30