设备

个人电脑与Windows，MacOS的，或基于Unix的操作系统

软件

R软件环境4.0.2。（https://www.r-project.org/）
RStudio集成开发环境1.3.1056（不是必需的，但强烈建议使用https://rstudio.com/）
GEOquery R包2.56.0 （Davis和Meltzer，2007），https: //www.bioconductor.org/packages/release/bioc/html/GEOquery.html
Minfi R package 1.34.0 （Aryee et al。，2014），https: //www.bioconductor.org/packages/release/bioc/html/minfi.html
M3C R包1.10.0 （John等，2020），https://www.bioconductor.org/packages/release/bioc/html/M3C.html
ggfortify R包0.4.11 （Yuan et al。，2016），https: //cran.r-project.org/web/packages/ggfortify/index.html
DMRcate R包2.2.3 （Peters等人，2015），https: //www.bioconductor.org/packages/release/bioc/html/DMRcate.html
wateRmelon R软件包1.32.0 （Pidsley等，2013），https: //www.bioconductor.org/packages/release/bioc/html/wateRmelon.html
FlowSorted.Blood.450k R软件包1.26.0，http： //bioconductor.org/packages/release/data/experiment/html/FlowSorted.Blood.450k.html
mCSEA R软件包1.8.0 （Martorell-Marugánet al。，2019），http: //bioconductor.org/packages/release/bioc/html/mCSEA.html

数据集

先前生成的甲基化数据（Strong等人，2015）。GEO标识符：GSE66552（https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66552）

程序

安装必要的R软件包
可以在R终端中使用以下命令从CRAN和Bioconductor存储库安装运行此管道所需的软件包：

install.packages （“ ggfortify ”）

if （！requireNamespace （“ BiocManager ”，静默= TRUE））

install.packages （“ BiocManager ”）

BiocManager :: install（c（“ GEOquery ”，“ minfi ”，“ M3C ”，“ DMRcate ”，“西瓜”，

“ FlowSorted.Blood.450k”，“ mCSEA ”））

下载原始数据
图1示出的图的数据分析的工作流程。Illumina的BeadChip阵列原始数据是idat文件，其中包含每个微阵列点的荧光强度。这些IDAT文件可以从公众基因表达进行下载[R epository GEO （埃德加等人，2002年）使用该GEOquery包（Davis和梅尔策，2007年）。在这个协议中，我们将到ANALY ž由强报告Ë数据等。（2015），可以使用GEO标识符GSE66552检索。在R会话中，执行以下代码：

图书馆（GEOquery ）

getGEOSuppFiles （“ GSE66552”）

下载完成后，当前目录中将存在一个包含压缩数据的文件夹。使用以下代码将数据解压缩到IDAT文件夹中：

解压缩（“ GSE66552 / GSE66552_RAW.tar”，exdir =“ IDAT”）

最后，将实验的元数据保存到名为pheno的变量中：

GEOData <-getGEO（“ GSE66552”）

现象<-pData（phenoData（GEOData [[1]]））

图1.数据分析工作流。黑盒子包含输入，输出，和中间文件。蓝色框包含此协议中描述的步骤。

读取idat文件
Minfi软件包（Aryee等，2014）可用于读取idat文件并处理甲基化数据。当t他的目的，目标数据框应该从实验的元数据来创建。表1包含此数据集的前5行，以显示此对象应具有的结构。

图书馆（minfi ）

basenames_split < -strsplit （pheno $ supplementary_file ，“ suppl /”）

basenames < -sapply （basenames_split ，function（x）{return（x [2]）}）

目标< -data.frame （“ Sample_Name ” = pheno $ title ，“ Basename ” = basenames ）

表1.第5行的所述目标数据帧

读取数据并将其存储在RGChannelSet对象中，该对象将用于以后的分析。RGChannelSet是minfi用于存储原始甲基化数据的类。此对象中的样本名称也可以在此步骤中指定。

rgchannelset < -read.metharray.exp （目标=目标，基准=“ IDAT”）

sampleNames （rgchannelset ）<- targets $ Sample_Name

执行质量控制（QC）
为确保没有影响任何样品的技术问题，应执行质量控制步骤。Minfi程序包包含一项功能，可探索每个样品的未甲基化通道和甲基化通道的中值强度。低强度可以指示问题小号与样品。图2 A显示了此QC的输出图。在这种情况下，没有任何样品未通过此质量控制。如果任何样品未通过此质量控制，则应将其从分析中丢弃。

msetraw < - preprocessRaw （rgchannelset ）

qc < -getQC （msetraw ）

plotQC （qc）

主成分小号分析（PCA）也是推荐的QC步骤，检测潜在的异常值。为此，可以使用M3C （John等人，2020）和ggfortify （Yuan等人，2016）R包。图2B显示了我们示例的PCA图。如可以观察到，存在有没有异常值样品小号和基团被部分地分离。组之间可观察到的混合可能是由于协变量，如性别，年龄，或患者的血统。

图书馆（M3C）

图书馆（ggfortify ）

betas.raw < -getBeta （msetraw ）

phenoPCA < -data.frame （“ Group” = pheno $`group：ch1`，stringsAsFactors = F ）

pcaResult < -prcomp （t（na.omit （betas.raw ）），scale。= TRUE）

自动绘图（pcaResult ，data = phenoPCA ，color =“ Group”）

图2.甲基化数据的质量控制。（A）Minfi QC ，每个样品的中位未甲基化和甲基化信号强度。（B）表示前两个主要成分的PCA图。样品按实验组着色。

预处理数据
牛逼Ø正确的技术偏见和提高数据质量，必须过滤掉一些探头和应用规范化方法。Minfi preprocessNoob函数应用Noob（正态指数带外），这是一种通过染料偏差归一化进行的背景校正（Triche等人，2013）。标准化输出是一个MethylSet对象，其中包含每个CpG探针的甲基化和未甲基化信号s 。可以将该对象转换为RatioSet类，在其中存储β值而不是信号。最后，应创建一个包含映射到基因组的β值的GenomicRatioSet对象。该GenomicRatioSet是用于进一步分析必要的类。

MSET < - preprocessNoob （rgchannelset ）

rset < -ratioConvert （mset ）

grset < -mapToGenome （rset ）

随后，可以获得每个探针和样品的检测P值。通过将探针信号与使用阴性对照探针测得的背景信号进行比较，可以计算出此P值。甲d etection P -值> 0.01指示š的不可信值得信号。过滤步骤可以应用于与以除去那些探针一个检测P样品中的百分比-值> 0.01（在这种情况下，样品的至少10％）。

detectionPvalues < - detectionP （rgchannelset ）

detectionPvalues < - detectionPvalues [ rownames （grset ），]

propGood < -apply （detectionPvalues ，1，function（x）{

长度（x [which （x <0.01）]）/ ncol （detectionPvals ）

}）

grset < -grset [which（propGood > = 0.9），]

我t被建议以除去那些探针位于性染色体（仅当被包括男性和女性的样本）（Maksimovic等人。，2017），含有单核苷酸多态性（SNP）的探针（纳伊姆等人，2014），和交叉反应性探针（陈等人。，2013）。DMRCate包（彼得斯。等人，2015）中包含的便利功能rmSNPandCH ，其可用于在一个步骤中所有这些探针筛选：

图书馆（DMRcate ）

beta.noob < -getBeta （grset ）

beta.noob.filtered < -rmSNPandCH （beta.noob ，rmXY = TRUE）

最后，可以使用Beta-Mixture Quantile（BMIQ）归一化来减轻450K和EPIC阵列中存在的I型和II型探针设计之间的系统差异（Teschendorff等，2013）。对于第是目的，所述西瓜包（Pidsley 。等人，2013），可以使用：

图书馆（wateRmelon ）

annot < -getAnnotation （rgchannelset ）

probeType < -as.data.frame （annot [rownames（beta.noob.filtered），c（“ Name”，“ Type”）]）

probeType < -ifelse （probeType $ Type ％in％“ I”，1,2）

beta.normalized < -apply （beta.noob.filtered ，2，function（x）{

wateRmelon :: BMIQ（x，design.v = probeType ）$ nbeta

}）

估计细胞类型异质性
血液甲基化数据变异性的主要来源之一是样本中不同的细胞类型比例（McGregor等，2016）。在这种情况下，可以使用FlowSorted.Blood.450k程序包中的参考数据来估计每个样本中的这些比例。在所述与EPIC平台产生的数据的情况下，FlowSorted.Blood.EPIC包应该被代替使用。

库（FlowSorted.Blood.450k）

cellCounts < - estimateCellCounts （rgchannelset ）

所述cellCounts对象是一个可与实验合并的矩阵人的元数据以考虑这些估计作为在进一步的分析协变量：

pheno = merge（pheno ，cellCounts ，by.x =“ title”，by.y =“ row.names ”）

行名（pheno ）<- pheno $ title

识别差异甲基化的启动子
最常见的甲基化分析之一是检测跨表型显示不同甲基化水平的区域，称为差异甲基化区域（DMR）。为此，不同的统计人程序已经开发出来，大部分的这些基于成熟的方法，如线性模型（里奇等人，2015）。在这个协议中，我们使用的mCSEA包（马托雷尔-马鲁甘等人，2019）来进行的DMRS分析集中于启动子区。第一步是使用线性模型根据不同条件下CpG探针的甲基化差异对其进行排名。在这些模型中，可以包括可用的协变量（例如，估计的细胞类型比例）以调整其对数据的可能影响。在这里，我们是第一ANALY ž荷兰国际集团的WS VS 。TD控件比较：

图书馆（mCSEA ）

phenoComp < -pheno [，c（“ group：ch1”，“ CD8T”，“ CD4T”，“ NK”，“ Bcell ”，“ Mono”，“ Gran”）]]

rankCpGs1 < -rankProbes （beta.normalized ，phenoComp ，

refGroup =“ TD控件”，caseGroup =“ WS”，

协变量= c（“ CD8T”，“ CD4T”，“ NK”，“ Bcell ”，“ Mono”，“ Gran”），

连续= c（“ CD8T”，“ CD4T”，“ NK”，“ Bcell ”，“ Mono”，“ Gran”））

下面步骤是使用于计算差异甲基化的启动子。输出对象是一个数据帧与基因符号行的名字，在估计的P值的PVAL柱，和P-值通过调整所述的Benjamini -Hochberg方法（和的Benjamini霍赫贝格，1995）中所述padj柱。此外，鉴于THA吨的mCSEA算法是基于基因集富集分析（GSEA）（萨勃拉曼尼亚等人，2005），富集得分（ES）被提供给每个基因。这个分数的符号指示是否启动子超甲基化的（正号）或低甲基化（负号）的情况下组中作为比较用的控制。最后，将归一化的ES（NES）和量的CpG中的每个启动子中报告的NES和大小的列，分别。表2显示了results_promoters1对象的前20行。在该EPIC数据的情况下，参数平台=“EPIC”应该被添加到mCSEATest功能。此外，该regionsTypes参数可以被修改以“基因”或“CGI”到ANALY Ž e本基因体和CGI ，分别。

种子（123）

results_mCSEA1 < -mCSEATest （rankCpGs1，beta.normalized ，

phenoComp ，regionsTypes =“发起人”）

results_promoters1 <-results_mCSEA1 $ promoters [order（results_mCSEA1 $ promoters $ padj），

- C（3,7）]

表2. WS vs. WS中排名前20位的差异甲基化启动子。TD控制比较。启动子按其调整后的P值排序。

对Dup7与进行相同的分析。TD控件比较：

rankCpGs2 < -rankProbes （beta.normalized ，phenoComp ，

refGroup =“ TD控件”，caseGroup =“ Dup7”，

协变量= c（“ CD8T”，“ CD4T”，“ NK”，“ Bcell ”，“ Mono”，“ Gran”），

连续= c（“ CD8T”，“ CD4T”，“ NK”，“ Bcell ”，“ Mono”，“ Gran”））

种子（123）

results_mCSEA2 < -mCSEATest （rankCpGs2，beta.normalized ，

phenoComp ，regionsTypes =“发起人”）

results_promoters2 <-results_mCSEA2 $ promoters [order（results_mCSEA2 $ promoters $ padj），

- C（3,7）]

此外，可以在其基因组背景下代表任何启动子的甲基化值。图3A显示了该CRISP2基因图的结果，该基因的启动子在WS患者中高甲基化，在Dup7患者中低甲基化。要构建这些图，请使用mCSEAPlot函数：

mCSEAPlot （results_mCSEA1， “启动子”， “CRISP2”，makePDF = F，前缘= F）

最后，还可以为每个启动子生成一个GSEA图。在该图中，所有的ANALY Ž ED的CpG被排序它们差异性甲基化的x轴，和垂直线标记的CpG属于表示启动子。图3B包含CRISP2 GSEA图，可以使用mCSEAPlotGSEA函数生成该图：

mCSEAPlotGSEA （排名为CpGs1，results_mCSEA1，“发起人”，“ CRISP2”）

图3. mCSEA地块进行了CRISP2基因。（A）基因组的上下文中的CRISP2启动子。点代表每个CpG和样本的β值。线显示每个实验组的平均甲基化。（B）GSEA为情节的在WS CRISP2启动子VS 。TD控制比较。垂直线标记CRISP2启动子的位置的CpG沿ANALY列表ž ED的CpG （水平黑线）。

结论：

在该协议中，我们展示了如何使用R软件包对Illumina BeadChip甲基化数据进行整体分析，以鉴定疾病样品中差异甲基化的启动子。该协议可以直接应用于任何病例对照研究设计。此外，它可以很容易地用于研究其他区域（如基因体和CGI）的差异甲基化。我们相信，这种工作流程将对研究与不同疾病相关的甲基化改变的研究人员有用。

致谢

尔迪Martorell的-马鲁甘被部分资金部：德Economía ，INDUSTRIA ÿ Competitividad 。这项工作是由部分支持FEDER /军政府的Andalucía- Consejería德Economía Ÿ Conocimiento （批准号：CV20-36723 ）和Consejería德Salud的，安达卢西亚委员会（格兰特PI - 0173 - 2017年）。我们承认原有的研究论文中，该协议是派生（Martorell的-马鲁甘等人。，2019）。

利益争夺

作者宣称他们没有竞争利益。

参考

1. Aryee，MJ，Jaffe，AE，Corrada-Bravo，H.，Ladd-Acosta，C.，Feinberg，AP，Hansen，KD，and Irizarry，RA（2014）。Minfi：一种灵活而全面的Bioconductor软件包，用于分析Infinium DNA甲基化微阵列。生物信息学30 ：1363 - 1369年。

2.的Benjamini ，ÿ 。和霍赫贝格，ÿ 。（1995）。控制˚F ALSE d iscovery ř吃：一个p ractical和p owerful一个接近角到米ultiple吨esting。JR统计志57 ：289 - 300。

3. Boyes ，J 。和Bird ，A 。（1992）。DNA甲基化对基因的抑制取决于CpG密度和启动子强度：甲基CpG结合蛋白参与的证据。EMBO J.11 （1）：327-333。

4. Cedar，H.和Bergman，Y.（2012）。DNA甲基化模式的编程。Annu启生物化学81 ：97 - 117。

5. Chen，Y.，Lemire，M.，Choufani，S.，Butcher，DT，Grafodatskaya，D.，Zanke，BW，Gallinger，S.，Hudson，TJ和Weksberg，R.（2013）。在Illumina Infinium HumanMethylation450微阵列中发现交叉反应探针和多态CpG。表观遗传学8 ：203 - 209。

6. Davis，S.和Meltzer，PS（2007）。GEOquery：基因表达综合（GEO）与BioConductor之间的桥梁。生物信息学23 ：1846年- 1847年。

7. Dozmorov，MG，Wren，JD和Alarcón-Riquelme，ME（2014年）。表观基因组元素富含自身免疫敏感性基因的启动子。表观遗传学9 ，276-285。

8. R. Edgar，M. Domrachev和AE，Lash（2002）。基因表达综合：NCBI基因表达和杂交阵列数据存储库。核酸研究。30 ：207 - 210。

9. Ehrlich，M.和Lacey，M.（2013）。癌症中的DNA次甲基化和半甲基化。生物医学进展754：31-56。

10.弗拉纳根，JM（2015）。表观基因组范围的关联研究（EWAS）：过去，现在和将来。方法分子生物学克利夫顿NJ 1238 ：51 - 63。

11.约翰·华润，华盛顿·沃森，罗斯，华盛顿，戈德曼，爱伦斯坦，皮茨卡里斯，路易斯和米纳斯·巴恩斯（2020）。M3C：基于蒙特卡罗参考的共识聚类。科学代表10 ：1816。

12. Maksimovic，J.，Phipson，B.和Oshlack，A.（2016）。用于分析甲基化阵列数据的交叉包装Bioconductor工作流程。F1000Res 5：1281。

13.马尔托雷尔-马鲁甘，J.，冈萨雷斯- Rumayor，V.，和卡莫纳-Sáez研究，P。（2019）。mCSEA：检测微妙的差异甲基化区域。Bioinforma OXF英格兰35 ：3257 - 3262。

14. K. McGregor，S。Bernatsky，I。Colmegna，M。Hudson，T。Pastinen，A。Labbe和CMT Greenwood（2016）。对DNA甲基化研究中纠正细胞类型异质性的方法的评估。基因组生物学杂志17 ：84。

15. Naeem，H.，Wong，NC，Chatterton，Z.，Hong，MKH，Pedersen，JS，Corcoran，NM，Hovens，CM和Macintyre，G.（2014年）。减少错误发现的风险，从而可以使用HumanMethylation450阵列鉴定生物学上重要的全基因组甲基化状态。BMC基因组学15 ：51。

16.彼得斯，新泽西州，巴克利，密西根州，斯坦森，阿拉巴马州，皮兹利，R。，萨马拉斯，K.，V·洛德，R。，克拉克，SJ和莫洛伊（PL）（2015年）。从头鉴定人类基因组中差异甲基化的区域。表观遗传染色质8 ：6。

17. R. Pidsley，R. Y Wong，CC，Volta，M.，Lunnon，K.，Mill，J.和Schalkwyk，LC（2013）。一种数据驱动的方法来预处理Illumina 450K甲基化阵列数据。BMC基因组学14 ：293。

18.缅因州里奇（Ritchie，ME），皮普森（B. limma为RNA测序和微阵列研究提供了差异表达分析的能力。核酸研究43 ：e47 - e47。

19.罗伯逊，KD（2005）。DNA甲基化与人类疾病。纳特冯遗传学6 ：597 - 610。

20. Strong，E.，DT，Butcher，DT，Singhania，R.，Mervis，CB，Morris，CA，De Carvalho，D.，Weksberg，R。和Osborne，LR（2015年）。删除或重复7q11.23儿童的对称剂量依赖性DNA甲基化分布图。上午Ĵ人类遗传学97 ：216 - 227。

21. Subramanian，A.，Tamayo，P.，Mootha，VK，Mukherjee，S.，Ebert，BL，Gillette，MA，Paulovich，A.，Pomeroy，SL，Golub，TR，Lander，ES等。（2005）。基因集富集分析：一种基于知识的方法来解释全基因组表达谱。PROC国家科科学院院刊102 ：15545 - 15550。

22. Teh，AL，Pan，H.，Lin，X.，Lim，YI，Patro，CPK，Cheong，CY，Gong，M.，MacIsaac，JL，Kwoh，C.-K.，Meaney，MJ等al。（2016）。在临床样品中进行甲基化组分析的甲基捕获测序与Infinium 450K甲基化阵列的比较。表观遗传学11 ：36 - 48。

23. AE的Teschendorff，F 。的Marabita，M。的Lechner，T。的Bartlett，J。的Tegner，D。的Gomez-Cabrero和S.的Beck（2013）。一种用于纠正Illumina Infinium 450 k DNA甲基化数据中探针设计偏倚的β混合物分位数归一化方法。Bioinforma OXF英格兰29 ：189 - 196。

24. Tche的Triche，DJ的Weisenberger，D。的Van Den Berg，PW的Laird和KD的Siegmund（2013）。Illumina Infinium DNA甲基化BeadArrays的低级处理。核酸研究41 ：e90 - e90。

25. Yuan T.，Horikoshi M.和Li W.（2016）。ggfortify：统一界面，用于可视化常用R包的统计结果。RJ 8 ：474 - 485。