一、甲基化芯片检测
- 1.1 DNA甲基化
- 1.2 甲基化芯片原理
- 1.3 β值
- 1.4 分析需要考虑的问题
二、甲基化芯片数据分析
- 2.1 Pipeline
- - 2.1.1 450K
  - 2.1.2 EPIC(850k)
- 2.2 数据预处理
- - 2.2.1 数据读取
  - 2.2.2 质控与标准化
- 2.2 甲基化分析
- - 2.2.1 Differential Methylation Probe（差异化CpG位点）
  - 2.2.2 Differential Methylation Region（差异化CpG区域）
  - 2.2.3 Differential Methylation Block（更大范围的差异化region区域）
- 2.4 富集分析
三、参考
四、作业
- 4.1 练习一
- 4.2 练习二
- 4.3 练习三

一、甲基化芯片检测

1.1 DNA甲基化

DNA甲基化是表观遗传学的中最为常见的一种修饰，其主要形式包括：5-甲基胞嘧啶 (5-mC)、少量的N6-甲基腺嘌呤 (N6-mA) 以及7-甲基鸟嘌呤(7-mG)。

目前常说的DNA甲基化一般指CpG岛甲基化，即在DNA甲基化转移酶（DNMTs）的作用下使CpG二核苷酸5’端的胞嘧啶转变为5’甲基胞嘧啶。

CpG岛（CpG islands）:指CpG序列密度相比整个基因组来说是特别高的富集区域，一般位于启动子附近，5’端非翻译区或第一个外显子；一般CpG岛序列长度在500bp以上，GC含量高于55%以及CpG出现比率大于0.65，40%的启动子区域含有CpG岛。
CpG shores指距CpG岛边缘2kb的区域
CpG shelves是指距CpG岛边缘4kb的区域

1.2 甲基化芯片原理

甲基化芯片的原理是基于亚硫酸盐处理后的DNA序列杂交的信号探测。

亚硫酸盐是甲基化探测的“金标准”，不管是芯片或者甲基化测序，都要先对DNA样品进行亚硫酸盐处理，使非甲基化的C变成U，而甲基化的C保持不变，从而在后续的测序或者杂交后区分出来。
450K采用了两种探针对甲基化进行测定，Infinium I采用了两种bead（甲基化M和非甲基化U，如图显示），而Infinium II只有一种bead（即甲基化和非甲基化在一起），这也导致了它们在后续荧光探测的不同，450K采用了两种荧光探测信号（红光和绿光）。

上图左边一列是非甲基化的GpC locus，右边是甲基化的GpC locus，上下分别是Infinium I 和Infinium Ⅱ

1.3 β值

通过计算甲基化（信号A）和未甲基化（信号B）等位基因之间的强度比来确定DNA甲基化水平（β值）。

平均β=信号B /(信号A +信号B + 100)

具体地，β值是由甲基化（M对应于信号A）和未甲基化（U对应于信号B）等位基因的强度计算的，荧光信号的比率β= Max（M，0）/ [Max（ M，0）+ Max（U，0）+ 100]。
因此，β值的范围从0（完全未甲基化）到1（完全甲基化）

具体的β值的意义是：
任何等于或大于0.6的β值都被认为是完全甲基化的。
任何等于或小于0.2的β值被认为是完全未甲基化的。
β值在0.2和0.6之间被认为是部分甲基化的。

1.4 分析需要考虑的问题

背景校正
红光和绿光的校正
控制芯片的使用（illumina450K本身有一些控制芯片，可以用来做质控，如亚硫酸盐处理效率）
探针类型（I型和II型）的校正（不同探针类型产生的数据不同）
位置的校正（芯片上的不同位置产生的数据可能会有偏差）
批次的校正（不同的批次做的数据会有偏差）
探针序列本身是否可靠（有些探针本身位于repeat区或者包含snp等就会影响杂交及最后的结果，应该去除，附上一片参考文献，里边有list可以用来去除不好的探针）

二、甲基化芯片数据分析

图源：ChAMP: updated methylation analysis pipeline for Illumina BeadChips

绿色发光线表示主要的分析步骤，灰色为可选的步骤。黑点表示准备好的甲基化数据。
蓝色表示准备工作，比如Loading, Normalization, Quality Control checks etc.
红色表示产生分析结果：Differentially Methylated Positions (DMPs), Differentially Methylated Regions (DMRs), Differentially methylated Blocks, EpiMod (a method for detecting differentially methylated gene modules derived from FEM package), Pathway Enrichment Results etc.
黄色表示交互界面画图

2.1 Pipeline

2.1.1 450K

myLoad <- cham.load(testDir)
# Or you may separate about code as champ.import(testDir) + champ.filter()
CpG.GUI()
champ.QC() # Alternatively: QC.GUI()
myNorm <- champ.norm()
champ.SVD()
# If Batch detected, run champ.runCombat() here.
myDMP <- champ.DMP()
DMP.GUI()
myDMR <- champ.DMR()
DMR.GUI()
myBlock <- champ.Block()
Block.GUI()
myGSEA <- champ.GSEA()
myEpiMod <- champ.EpiMod()
myCNA <- champ.CNA()# If DataSet is Blood samples, run champ.refbase() here.
myRefbase <- champ.refbase()

2.1.2 EPIC(850k)

myLoad <- champ.load(directory = testDir,arraytype="EPIC")
# We simulated EPIC data from beta value instead of .idat file,
# but user may use above code to read .idat files directly.
# Here we we started with myLoad.CpG.GUI(arraytype="EPIC")
champ.QC() # Alternatively QC.GUI(arraytype="EPIC")
myNorm <- champ.norm(arraytype="EPIC")
champ.SVD()
# If Batch detected, run champ.runCombat() here.This data is not suitable.
myDMP <- champ.DMP(arraytype="EPIC")
DMP.GUI()
myDMR <- champ.DMR()
DMR.GUI()
myDMR <- champ.DMR(arraytype="EPIC")
DMR.GUI(arraytype="EPIC")
myBlock <- champ.Block(arraytype="EPIC")
Block.GUI(arraytype="EPIC") # For this simulation data, not Differential Methylation Block is detected.
myGSEA <- champ.GSEA(arraytype="EPIC")
myEpiMod <- champ.EpiMod(arraytype="EPIC")# champ.CNA(arraytype="EPIC")
# champ.CNA() function call for intensity data, which is not included in our Simulation data.

ps. 要做好装包装半天的准备！

2.2 数据预处理

2.2.1 数据读取

本次实验数据来源(GSE108785)：

文献：2019-Front. Neurosci-Genome-Wide DNA Methylation Analysis Reveals Epigenetic Pattern of SH2B1 in Chinese Monozygotic Twins Discordant for Autism Spectrum Disorder

数据：Query DataSets for GSE108785

.idat files 为原始芯片文件，还需pd file (Sample_Sheet.csv)文件（表型，编号等），整理至一个目录下：

sample_sheet按照下载的数据进行填写：

强迫症还可以把下载得到的文件名前面的GSM给去掉（bash 删11个字符）：

for file in `ls *.idat.gz`;do mv $file `echo $file|sed 's/.\{11\}//'`;done;

然后利用champ.load()读取数据：
testDir为文件存放的路径

testDir <- "/public/workspace/stu18230130/MyCourse/BioChip/lab6_ChAMP/MethyData/"
myLoad <- champ.load(testDir,arraytype="450K")
str(myLoad)## 查看beta value
myLoad$beta[1:6,]

champ.load()其实等于champ.import() + champ.filter()，读数据和过滤一步到位。
过滤标准：

过滤掉detection p-value大于0.01的探针
过滤掉在至少5%样本中bead count小于3的探针
过滤掉非GpC位点的探针
过滤掉所有SNP相关的探针
过滤掉multi-hit探针，即映射到多个位置的
过滤掉X和Y染色体上的探针

ChAMP包还提供了CpG.GUI()函数用于展示distribution of CpGs on chromosome, CpG island, TSS reagions. e.g，其实是几张简单展示myLoad$beta分布的图片，但是界面很友好

2.2.2 质控与标准化

champ.QC()和QC.GUI函数用于检查看看过滤后的数据是否符合要求，能否进行下一步的分析。

champ.QC()函数结果三张图

mdsPlot (Multidimensional Scaling Plot)，基于前1000个最易变化的位点查看样本的相似度，用颜色标记不同的样本分组,主要看看不同分组样本是否分开；
densityPlot，查看每个样本的beta值分布，有严重偏离的样本预示着质量较差（如亚硫酸盐处理不完全等；
dendrogram，样本的聚类图

Normalization
type-I 和 type-II 两种探针化学反应不同，设计也不同，导致分布区域也不同。这两种探针检测出的差异可能是因为探针所在位置不平衡导致的生物学差异引起的（如CpG位置的差异引起的）。
因此，针对 type-II probe bias的矫正是必要的。champ.norm() 函数可以实现这个功能。针对type-II 探针有4种标准化的方法：BMIQ, SWAN, PBC 和 unctionalNormliazation。

myNorm<-champ.norm(arraytype="450K")
write.csv(myNorm,file="./Normalization Data.csv",quote=F,row.names = T)

SVD
The singular value decomposition method (SVD) 用来用于评估数据集中变量的主要成分。这些显著性位点可能与我们感兴趣的生物学现象相关联，也可能与技术相关，如批次效应或群体效应。

champ.SVD()`

champ.SVD()函数将把pd文件中的所有协变量和表型数据纳入进行分析。但是对于分类变量和数字变量处理方法是不一样的。分类变量要转换成“factor” or “character”类型，数字变量转换成数字类型。
champ.SVD()分析时会把协变量打印在屏幕上，结果是热图，保存为SVDsummary.pdf文件。黑色表示最显著的p值。如果发现技术因素有影响，就需要用ComBat等方法重新标准化数据，包括variation related to the beadchip, position and/or plate。

图(a)是自己选的6个样本的绘制结果，可以看出它并没有combat，图(b)是用自带的测试数据绘制的，不是很复杂，看不出来。图©是用GSE40279的656个样本绘制的，其中年龄是数字变量，其他都为分类变量，可以看出存在很明显的批次效应及混杂因素，需进行处理：

#矫正批次效应，这步非常耗时。
myCombat <- champ.runCombat(beta=myNorm,pd=myLoad$pd,batchname=c("Slide"))
#查看矫正后 的结果
champ.SVD()

2.2 甲基化分析

2.2.1 Differential Methylation Probe（差异化CpG位点）

DMP（Differential Methylation Probe），用于找出差异的甲基化位点，然后ChAMP包已经将分析过程（主要还是基于Limma包的linear model，这个包是专门用于分析芯片差异表达的包）；然后根据你的分组信息，获得差异甲基化位点。champ.DMP()返回的是list。

champ.DMP(beta = myNorm,
pheno = myLoad$pd$Sample_Group,
compare.group = NULL,
adjPVal = 0.05,
adjust.method = "BH",
arraytype = "450K")

myDMP <- champ.DMP(arraytype="450K",adjPVal = 0.05,compare.group=c("C","T"))
## 查看结果
head(myDMP[[1]])
## 保存DMP
write.csv(myDMP,file="./DMP_analysis_result.csv",quote=F,row.names = F)

2.2.2 Differential Methylation Region（差异化CpG区域）

DMRs主要指一连串的CpG都会出现很明显的差异，champ.DMR()函数计算并返回一个数据框，包括：detected DMRs, with their length, clusters, number of CpGs annotated.
函数包含三种算法Bumphunter, ProbeLasso and DMRcate.。Bumphunter比较可靠，精确度可以有90%以上，ProbeLasso有75%左右。Bumphunter 算法首先将所有的探针分成几小类，然后用随机permutation方法评估候选的DMRs.

champ.DMR(beta=myNorm,
pheno=myLoad$pd$Sample_Group,
compare.group=NULL,
arraytype="450K",
method = "Bumphunter",
minProbes=7,
adjPvalDmr=0.05,
cores=3)

myDMR <- champ.DMR(arraytype="450K",adjPvalDmr=0.05,cores=8,minProbes=5)
write.csv(myDMR,file="./DMR_analysis_result.csv",quote=F,row.names = F)

2.2.3 Differential Methylation Block（更大范围的差异化region区域）

在Block-finder 功能中，champ.Block()函数首先在全基因组范围上计算small clusters (regions) ，然后对于每个cluster，计算平均值和位置，将每个区域压缩为一个单元。

champ.Block(beta=myNorm,
pheno=myLoad$pd$Sample_Group,
arraytype="450K",
maxClusterGap=250000,
B=500,
bpSpan=250000,
minNum=10,   #Threshold to filtering Blocks with too few probes in it
cores=3)

myBlock <- champ.Block(minNum=5,cores=8)
write.csv(myBlock$Block,file="./BLOCK_analysis_result.csv",quote=F,row.names = F)

2.4 富集分析

当我们得到差异的探针或者差异的甲基化区域之后，通常都会分析这些差异区域对应的基因是否在特定功能上有富集。在ChAMP中，通过champ.GSEA函数来实现功能富集分析。

champ.GSEA()默认对差异CpG位点和差异甲基化区域对应的基因做富集分析，采用的方式是Fisher exact test，分析的是Gene Set 来自MSigDB。

champ.GSEA(beta=myNorm,
DMP=myDMP[[1]],
DMR=myDMR,
CpGlist=NULL,
Genelist=NULL,
pheno=myLoad$pd$Sample_Group,
method="fisher",
arraytype="450K",
Rplot=TRUE,
adjPval=0.05,
cores=1)

myGSEA <- champ.GSEA()

默认对DMP和DMR对应的基因都是富集分析，所以结果是一个长度为2的列表，第一个列表是DMP富集分析的结果，第二个列表是DMR富集分析的结果，每个富集结果是一个data.frame对象。(由于数据原因，这里只富集出了DMR，与设定的阈值有关)

如果要做GO/KEGG 富集分析，同时又要考虑CpG位点个数，就必须考虑missMethyl包。

最近还新开发了基于经典贝叶斯的无偏GSEA方法

champ.ebGSEA(beta=myNorm, pheno=myLoad$pd$Sample_Group, minN=5, adjPval=0.05, arraytype="450K", cores=1)

三、参考

ChAMP 包分析甲基化数据
Package ‘ChAMP’-Bioconductor Maunal
2017-Bioinformatics-ChAMP: updated methylation analysis pipeline for Illumina BeadChips
R: 使用CHAMP包进行甲基化数据的差异分析(QC, CNV, DMP, DMR等)
Infinium Methylation Assay Overview
Illumina HumanMethylation450 BeadChip (甲基化450k芯片) 预处理初探

四、作业

4.1 练习一

1）提取chr11上DMR的数据：grep()/which()
2）提取chr1的数据：which()
3）将DMR结果写入到本地文件

利用which提取chr11上的DMR

 chr11_dat <- myDMR[[1]][which(myDMR[[1]]$seqnames == "chr11"),]

2. 利用which提取chr1的数据，并将其整合到List

4.2 练习二

对标准化后的beta用wilcox rank test 寻找差异甲基化位点，并和champ.DMP结果比较

对标准化后的数据进行秩和检验，输出p value向量

apply_P <- apply(norm.dat, 1, test)test <- function(x){res <- wilcox.test(x[1:4], x[5:8])res$p.value
}

需要补充信息，不然无法和champ.DMR()的结果比较

利用注释数据hm450.manifest.hg19，选取其中有特征的几列

将其与apply_P进行合并操作：

 res <- merge(probe_info,apply_P,by.x = "probeID",by.y = 0)

按P值升序取前30个与DMR结果进行比较：

champ.DMR()按照p值升序排序后结果：

好像不是很匹配，我觉得应该是因为做秩和检验的样本量太少了，统计学效能不够，在做wilcoxin的时候出现了一个warning，好像说的也是这个问题~

4.3 练习三

自己下载一套illumina甲基化芯片
1）该数据参考基因组是？请写出追溯流程
2）完成上述处理流程

参考基因组：在读入数据后会生成“hm450.manifest.hg19"注释文件，可以看出参考基因组是hg19
处理流程即本文~

【Bioinfo Blog 013】【R Code 011】——甲基化芯片数据分析（ChAMP包）相关推荐

甲基化系列 3. 甲基化芯片数据分析完整版（ChAMP）
点击关注,桓峰基因桓峰基因生物信息分析,SCI文章撰写及生物信息基础知识学习:R语言学习,perl基础编程,linux系统命令,Python遇见更好的你 104篇原创内容公众号桓峰基因的教程不 ...
【Bioinfo Blog 006】【R Code 005】——GEO表达谱数据质控
目录一.数据下载及质控 1.1 GEO数据下载 1.1.1 GEOquery包安装 1.1.2 .cel数据下载 1.2 读取.cel文件 1.2.1 Affy包安装 1.2.2 利用ReadAff ...
r语言怎么把txt数据变成一个Rdata格式_甲基化芯片数据下载如何读入到R里面
数据是一切的开始,前面我们介绍了一些背景知识,主要是理解什么是DNA甲基化,为什么要检测它,以及芯片和测序两个方向的DNA甲基化检测技术.具体介绍在:甲基化的一些基础知识,也了解了甲基化芯片的一般分析 ...
amt630a芯片中文资料_甲基化芯片学习记录
今天是生信星球陪你的第539天大神一句话,菜鸟跑半年.我不是大神,但我可以缩短你走弯路的半年~ 就像歌儿唱的那样,如果你不知道该往哪儿走,就留在这学点生信好不好~ 这里有豆豆和花花的学习历程,从新 ...
国家生物信息中心开发DNA甲基化芯片数据标准化方法—GMQN
过去十年来,由于DNA甲基化芯片技术的不断发展以及测序成本的快速下降,DNA甲基化芯片数据呈现爆发式增长.这些数据是表观基因组关联研究(Epigenome-Wide Association Studi ...
UA MATH571A 回归分析概念与R code总结
UA MATH571A 回归分析概念与R code总结 Simple Linear Regression Multivariate Linear Regression Part 0 Basic R ...
用R和BioConductor进行基因芯片数据分析(四)：芯片内归一化
接前一篇: 用R和BioConductor进行基因芯片数据分析(三):计算median 归一化是从normalization翻译过来的.归一化的目的是使各次/组测量或各种实验条件下的测量可以相互比较, ...
甲基化芯片入门学习-基础知识（一）
基本概念梳理什么是DNA甲基化 DNA甲基化是表观遗传学的中最为常见的一种修饰,其主要形式包括:5-甲基胞嘧啶 (5-mC).少量的N6-甲基腺嘌呤 (N6-mA) 以及7-甲基鸟嘌呤(7-mG). ...
R code execution error处理
R code execution error 解决方法: Ctrl + Shift + F10 to restart your R session

【Bioinfo Blog 013】【R Code 011】——甲基化芯片数据分析（ChAMP包）

目录