850K甲基化芯片数据的分析

作者是生信技能树组建的表观遗传学学习小组的小组长，前面已经发过一个：

学员分享-Chip-seq 实战分析流程

本文是看到生信技能树有个450K甲基化芯片数据处理传送门，我呢，恰好不久前用一个集成度很高的ChAMP包分析过850K的甲基化芯片数据。所以，就想着把自己的笔记整理下，可以和更多的小伙伴学习交流，还有个原因可能是因为这是四月份打算学生信时，接手的第一个任务，曲曲折折好几个月才跑通流程，遇到的坑也比较多，想记录下来。

我之前分析时是参考ChAMP包的源文档，非常详细的整个流程的介绍，但是，在笔记快整理完时突然发现作者的博客也写过一篇介绍的文章，博客里写的不像源文档很官方，这里面有很多作者很直白的解释和补充，还有作者一些很深刻的思考。看了之后发现自己对很多分析理解的还不是很深刻。所以如果想学甲基化芯片数据分析的小伙伴可以以官方源文档和作者的博客为主，这篇笔记仅仅作为额外的参考吧。

---

Illumina甲基化芯片目前仍是很多实验室做甲基化项目的首选，尤其是对于大样本研究而言，其性价比相当高；目前在临床上应用还是很广的。这种芯片的发展主要经历了27K、450K以及850K（27K，450K，850K指能测到的CpG甲基化位点），目前积累的数据主要是450K芯片的，之后850K可能会成为甲基化芯片的主流。楼主之前写过一篇450K芯片预处理的帖子，其中很详细介绍了这种芯片的基础知识以及流程图和代码，大家可以先看看。芯片的处理流程一般就是：数据读入——数据过滤——数据校正——下游分析。

数据处理一种时基于GenomeStudio（illumina开发的软件），但是只对于小样本，另一种基于R的各种package，如lumi、minfi、wateRmelon、ChAMP等。

与测序相比，芯片的处理可能对计算资源的要求不算高，主要使用的工具就是R，但是R的使用比较耗内存，尤其是处理大批量数据的时候。

Step1: 基础知识的补充

在正式分析前，我结合作业先将有关甲基化和芯片的基础知识整理了一下。

Illumina 甲基化芯片的原理及探针的设计（I型探针和II型探针)

原理：简而言之，基于亚硫酸盐处理后的DNA序列杂交的信号探测。亚硫酸盐是甲基化探测的“金标准”，不管是芯片或者甲基化测序，都要先对DNA样品进行亚硫酸盐处理，使非甲基化的C变成U，而甲基化的C保持不变，从而在后续的测序或者杂交后区分出来。

450K和850K采用了两种探针Infinium Ⅰ 和Infinium Ⅱ对甲基化进行测定，Infinium I采用了两种bead（甲基化M和非甲基化U，如图显示），而II只有一种bead（即甲基化和非甲基化在一起），这也导致了它们在后续荧光探测的不同，450K采用了两种荧光探测信号（红光和绿光）（图1）。

甲基化概述：

DNA甲基化被认为是表观遗传调控的一种方式，如Cytosine methylation (5-mC)是研究最多的，被认为是哺乳动物中常见的甲基化方式, 最近有一些研究也发现了其他形式的甲基化，如2016年Nature上发表了一篇关于鼠的胚胎干细胞的m6A（N6-methyladenine）形式的甲基化。DAN甲基化被认为对基因表达，染色质重塑，细胞分化，疾病等都有重要影响（图2）。

甲基化的检测方法：

目前甲基化检测的方法可以概括为三种：芯片、测序、免疫沉淀。具体选择何种方法主要还是根据实验目的和实验室条件了。但目前来说，甲基化芯片技术从覆盖度，检测灵敏度和价格综合考虑，还是性价比相对高的（图3）。

关于甲基化芯片常见的Glossary:

CpG island: Defned as regions 500 bp, 55% GC and expected/observed CpG ratio of 0.65. 40% of gene promoters contain islands.

CpG shelves: ~4Kb from islands.

CpG shores: ~2Kb from islands, 75% of tissuespecifc differentially methylated regions found in shores. Methylation in shores shows higher correlation with gene expression than CpG islands.

Differentially methylated regions (DMR): Cell-, tissue-, and condition- specifc differences in methylation.

Enhancer: A short region of DNA that can activate transcription and is often regulated by methylation.

Hypermethylation: Most cytosines are methylated. Hypomethylation: Most cytosines do not have 5-mC. Euchromatin and active gene promoters are hypomethylated.

Beta value:通常的甲基化衡量方法被称为“Beta”值; 等于甲基化百分比，并定义为“Meth”除以“Meth + Unmeth”。

CGI: CpG island 即甲基化岛。

---

因为手头的数据是850K的甲基化数据，之前也只接触过ChAMP包，所以这里就以ChAMP包介绍850K甲基化数据分析。ChAMP包是一个集成度很高的包，它包括450K和EPIC(即通常所说的850K)两套分析流程，完整的包括了数据的载入，标准化，矫正，差异甲基化和富集分析等功能（图4）。

Step2:计算机资源的准备

作业1 安装好R软件及相应的包，下载R包的说明书，整理它们的官网链接。

R的使用真的很耗内存，我有28个样本（14个control, 14个case), 之前4G内存的电脑，本地分析总时半路电脑就卡死了。所以最好配置高一点，或者在服务器上下载安装R和Rstudio（这里最好安装Rstudio, 因为ChAMP包中有很多的GUI图形功能，Rstudio可以更好实现，或者含有X11功能的linux系统）。

软件的安装：

R和Rstudio 的本地安装很简单，直接到官网下载，只要注意安装时的路径不要有中文，Rstudio安装前要先安装R。

服务器版本的Rstudio安装好后，在网页地址栏输入访问地址：<服务器IP:8787，用户名和密码为Linux用户的用户名和密码。

具体安装方法可以参考生信宝典陈老师的一篇文章http://www.biotrainee.com/thread-1808-1-1.html。

下载R包：

下载ChAMP 包，官网给出了很详细的流程说明（https://bioconductor.org/packages/release/bioc/vignettes/ChAMP/inst/doc/ChAMP.html）。

1. source("https://bioconductor.org/biocLite.R")

2. biocLite("ChAMP")

NOTE: ChAMP有很多依赖包，安装时，若报错有哪个包没有，就继续安装 biocLite("YourErrorPackage")，可能3-4次就可以安装成功。

导入ChAMP包并测试：

导入ChAMP包后，根据是450K的数据或者是850K的数据，导入测试数据集，走一下分析流程，检测包是否正常工作，更重要的是看该包的文档，理解每一步流程的意义。该包的文档很详细，建议大家看原文档，下面给出的啰啰嗦嗦的介绍基本上都来自官网的文档说明（https://bioconductor.org/packages/release/bioc/vignettes/ChAMP/inst/doc/ChAMP.html）。

1. library("ChAMP")

2. #450K的数据导入：

3. testDir=system.file("extdata",package="ChAMPdata")

4. myLoad <- champ.load(testDir,arraytype="450K")

5. #850K的数据

6. data(EPICSimData)

Step 3: 数据读入

数据的读入这里可能坑最多，首先450K 和850K甲基化芯片的原始数据格式都是IDAT, 因为数组是用两种不同的颜色来测量的，所以每个样本都有两个文件，通常是扩展名Grn.idat和Red.idat。数据在载入时还需要一个SampleSheet.csv文件（图5）（也称做pd file）, 这个文件很重要，它包含了样本的信息，可以对照测试数据的csv文件和自己的csv文件，对信息不全的地方进行补充。尤其要注意SampleGroup 这一列信息是否有，这一列信息代表你想比较的表型类型，比如癌和癌旁。另一个我遇到过的一个隐形坑在Sentrix_ID，这一列数因为数字串很长，在Excel中可能以科学计数法显示，然后本来是长数字串后两位不一样的数字串都变为一样的，在读入时就会报重复字符的错误，所以这里一定要核查下长数字串的信息，如果有错误，自己重新输入时以文档格式输入，或者前面加右单引'。csv文件准备好后，将csv文件与所有样本的芯片数据（即IDAT文件）放在一个文件下，然后就可以正常读入了。

图5 Sample_Sheet.csv fiel

1. library("ChAMP")

2. myLoad <- champ.load("F:/850K Methylation Chip/biotree_850K/methy_rawData",arraytype = "EPIC")

3. save(myLoad,file="myLoad.rda")

champ.load()包含了 champ.import() 和champ.filter()，这里会自动过滤p值0.01; probes beadcount <3 in at least 5% of samples;NoCG;probes with SNPs; MultiHit; probes located on X,Y chromosome。

在读入数据之后，最好保存，后续重复读入时会加快速度。

Step 4: 质控和标准化

CpG overview:

质控前可以先看看CpG的分布，包括在染色体上的分布；CpG岛附近的 open sea, shelf,shore (参考图2，理解具体意思) ; UTR,TSS; I 型探针和II探针上的分布（图6），这个信息对后续DMP的分析有帮助。

1. CpG.GUI(arraytype="EPIC")

质控:

然后进行质控，有两种方式：champ.QC() 和 QC.GUI()。champ.QC会产生三种类型的图（点图，beta 分布图，聚类图）以pdf格式输出，QC.GUI产生5个图，多了一个I型、II型探针图和热图（图7）。所有的GUI功能都比较耗内存，且产生的是网页交互式的图片，每幅图的右上角给的都有保存按钮，要注意的是保存时文件名要加上.png的后缀（图7）。

1. #champ.QC()

2. QC.GUI(arraytype="EPIC")

图7 QC Overview

标准化:

champ.norm 提供了四种方法：BMIQ, SWAN1, PBC2 and FunctionalNormliazation4。默认的方法是BMIQ, 且BMIQ对850K的标准化方法更好一点，所以这里我选择的是BMIQ的标准化方法，没有尝试其他的标准化方法。

1. myNorm <- champ.norm(arraytype="EPIC")

2. QC.GUI(myNorm,arraytype="EPIC")

3. save(myNorm,file="myNorm.rda")

SVD plot 和批次效应：

SVD(singular value decomposition) 这里用于评估数据集中变量的主要成分。这种成分可能确实是你感兴趣的生物因素，也可能是技术来源的一些变量成分（称为批次效应）（图8）。如果存在批次效应，就进行批次效应的矫正，矫正完之后可以再看看SVD plot。

1. champ.SVD()

图8 SVD Plot

Step 5: 差异甲基化分析（DMP & DMR & DMB)

差异分析是多数研究都要分析的，这里包括三种方法：DMP，DMR，DMB。DMP代表找出Differential Methylation Probe（差异化CpG位点），DMR代表找出Differential Methylation Region（差异化CpG区域），Block代表Differential Methylation Block（更大范围的差异化region区域）

简单来说，DMP是找出一个一个的差异甲基化CpG位点，DMR就是一个连续不断都比较长的差异片段，科学家们觉得，这样的连续差异片段，对于基因的影响会更加明显，只找这样的片段，可以使得计算生物学的打击精度更为准确，也可以让最终找出来的结论数据更少，便于实验人员筛选。另外一个类似的东西就是DMB，那个东西出现的原因是，有的科学家觉得，DMR这样的区域还不够显著，DNA上的甲基化出现变化，可能是绵延几千位点的！而且只会在基因以外的区域，但是这些基因以外的区域发生变化，却会导致基因的表达发生变化。你可以想象成，北京周边的河北在大炼钢铁，然后北京也跟着雾霾了，大概就是这意思。

DMP,DMR,DMB的结果都是基于的shiny的交互页面，左栏上方是 P-value 和 abs(logFC) ，可以选择想看的值，然后点submit, 右栏可以生成差异甲基化表，热图，feature&cgi, 左栏下方还有基因，CpG按钮，选择你想看的结果，submit, 右栏就会生成相应gene,CpG结果（图9）。

1. myDMP <- champ.DMP(arraytype="EPIC")

2. save(myDMP,file="myDMP.rda")

3. DMP.GUI()

4. myDMR <- champ.DMR(arraytype = "EPIC",method="DMRcate",cores=1)

5. save(myDMR,file="myDMR.rda")

6. DMR.GUI(arraytype="EPIC")

7. #myBlock <- champ.Block(arraytype = "EPIC")

8. #Block.GUI(arraytype="EPIC",compare.group=c("PrEC_cells","LNCaP_cells"))

图9 DMP Overview

Step 6: 基因富集和网络分析（GSEA & EpiMod）

差异甲基化分析后，你可能想知道DMP,DMR中涉及到的基因是否可以富集到某个生物功能或通路，GSEA(Gene Set Enrichment Analysis)和EpiMod（Differential Methylated Interaction Hotspots）提供了可以寻找作用通路网络中的疾病关联小网络的功能 （图 10）。

1. myGSEA <- champ.GSEA(arraytype = "EPIC")

2. save(myGSEA,file="myGSEA.rda")

3. myEpiMod <- champ.EpiMod(arraytype="EPIC")

4. save(myEpiMod,file="myEpiMod.rda")

Step 7: 拷贝数变异分析（CNA)

拷贝数变异，也就是有些基因片段被复制的此处过多或者过少，从而导致某些疾病。但是这个函数作者觉得有点粗糙，精度还不够。我试着跑了一下，时间超长（图11）。

1. myCNA <- champ.CNA(control = F,arraytype = "EPIC")

2. save(myCNA,file=myCNA)

图11 Frequency Plot of Cancer Sample

---

小结：如果用ChAMP包对450K或850K甲基化数据进行分析时，一是最好有个配置高一点的电脑；二是初始数据导入时，注意csv文件的格式，且要和IDAT文件放在一个文件下；其余的流程很少会遇到bug, 但最关键的是理解每一步的意义，能够根据分析的结果挖掘出想要的东西。

ps: 这次作业提供的公共数据，有IDAT文件，也有个csv文件，但是这里的csv文件和我的csv文件差别很大，不是很明白这里的csv文件是什么，有什么作用。

有不当之处欢迎指正和补充。

另外，如果还有想去南京面基的朋友，可以看看这个哈：江苏省生信专委会实战培训--非编码RNA专场  一起吃个饭咯