ChAMP包处理甲基化芯片数据

学习笔记

一、甲基化芯片
二、数据处理流程
三、详细步骤

一、甲基化芯片

甲基化定义
简介

探针类型
芯片注释信息

#查看450k芯片
library(IlluminaHumanMethylation450kmanifest)
show(IlluminaHumanMethylation450kmanifest)
#查看850k芯片信息使用另一个R包 IlluminaHumanMethylationEPICkmanifest
#minfi包有特定函数getManifest()可以查看指定对象对应芯片的注释信息

芯片注释文件

二、数据处理流程

#full pipeline in a step-wise process
library(ChAMP)
#-----使用原始数据:包含8个样本，4个肿瘤样本，4个对照样本
testdir<-system.file("extdata",package="ChAMPdata")
# 返回extdata所在的绝对路径(extdata中包含了testDataSet测试数据)
myLoad<-champ.load(testdir)
champ.QC()
myNorm<-champ.norm()
champ.SVD()
batchNorm<-champ.runCombat()
myDMP<-champ.DMP()
myDMR<-champ.DMR()
myDMB<-champ.Block()
myRefBase<-champ.refbase()
champ.CNA()
#找不到champ.MVP()和champ.reffree(),但是data(testDataSet)后有"myMVP"和"myRefFree"，说明这两个函数应该是旧版本的，但现在已经不使用了
#默认针对450K芯片的数据进行处理，如果要对850K芯片处理，在load、norm及后续识别甲基化的函数中设置参数arraytype="EPIC"

三、详细步骤

数据准备

1.使用ChAMPdata包中的数据（具体操作见后续）

2.数据库下载，以GEO为例说明数据准备

ChAMP包读取数据不仅要.idat原始文件，还需要一个pd(phenotype)文件，名为sample_sheet.csv（当然sample_sheet这个名字不重要，重要的是这是一个.csv文件）
这个文件记录了各样本信息，具体形式如下

(数据 GSE109904)可以看到没有sample_sheet.csv文件
（自认为该数据不适合用于甲基化分析）
GEO数据库中没有记录样本信息的pd文件，因此你需要根据各样本信息自行创建一个sample_sheet.csv

编辑.csv文件：

其中比较重要的属性列是Sample_Name、Sample_Group、Slide(Sentrix_ID)、Array(Sentrix_position)

注：使用的是GSE109904中4个使用850K芯片测得的样本，Sample_Name具体是什么不重要，只要能区分各个样本就行，如,可以命名为A/B/C/D

一步跑完整个流程

#-----一步跑完整个流程-----#
champ.process(testdir)
# Usage:
champ.process(runload=TRUE,#---是否运行champ.load()，运行后会加载一系列相关的变量directory = getwd(),#---.idat file所在路径filters=c("XY","DetP","Beads","NoCG","SNP","MultiHit"),#---数据预处理过程中要过滤掉的内容runimpute=TRUE,#---是否保留或填补由过滤导致的NAimputemethod="Combine",#---填补空缺值使用的方法runQC=TRUE,#---是否运行champ.QC()，用于可视化查看样本的相关信息，如是否有异常样本QCplots=c("mdsPlot","densityPlot","dendrogram"),#---指定使用champ.QC()要产生的图像runnorm=TRUE,#---是否运行champ.norm()，用于校正type Ⅱ探针normalizationmethod="BMIQ",#---指定运行champ.norm()时要使用的方法，除此之外还有PBC、SWAN、FunctionalNormalization方法runSVD=TRUE,#---是否运行champ.SVD()，可视化检验变异组分的显著性runCombat=TRUE,#---是否运行champ.runCombat()，用于处理批次效应batchname=c("Slide"),#---champ.runCombat处理批次效应时的批次标签runDMP=TRUE,#---是否运行champ.DMR()runDMR=TRUE,#---是否运行champ.DMR()DMRmethod="Bumphunter",#---运行champ.DMR时使用的方法，除此之外还有ProbeLasso方法runBlock=TRUE,#---是否运行champ.Block()，更大范围的差异甲基化区域runCNA=TRUE,#---是否运行champ.CNA()，用于CNA分析runRefBase=FALSE,#---是否运行champ.refbase(),用于细胞异质性处理resultsDir="./CHAMP_RESULT/",arraytype="450K",PDFplot=TRUE,Rplot=TRUE,saveStepresults=TRUE)

建议分步跑流程，一旦process出错

step by step跑完整个流程

1.数据读取和筛选

library(ChAMP)
testdir<-system.file("extdata",package="ChAMPdata")# 一、
#-----数据读取-------------------------------------
myImport<-champ.import(testdir)
# Usage:
champ.import(directory = getwd(),#---.idat file所在路径offset=100,#---beta.value <- M / (M + U + offset)，防止M和U太小而使分母为0arraytype="450K")#---芯片类型，可设置为"EPIC"用于加载850K芯片的数据#返回一个列表包含8个元素："beta"，"M"，"pd"，"intensity"，"detP"，"beadcount"，"Meth"，"UnMeth"
#beta等是一个矩阵，列代表样本，行代表探针
#pd记录样本信息，即自己编写的sample_sheet.csv文件，重要的属性列"Slide","Array"
#新增一个列表变量Anno："Annotation"，"ControlProbe"#-----数据过滤--------------------------------------
myFilter<-champ.filter(beta=myImport$beta,pd=myImport$pd)
# Usage:
champ.filter(beta=myImport$beta,#---原始数据中的β矩阵，甲基化程度:<未-部分-完全> → <0.2 - 0.6>M=NULL,#---logit(β)=log(Meth/Unmeth)pd=myImport$pd,#---pd(phenotype data)，来自Sample_sheet.csv文件，记录了每个样本的信息，如表型、编号等intensity=NULL,#---强度值矩阵Meth=NULL,#---甲基化矩阵UnMeth=NULL,#---非甲基化矩阵detP=NULL,#---Detected P value identifies failed positions defined as both the methylated and unmethylated channel reporting background signal levels.每一个芯片上的数据的可信度#---detP计算：https://blog.csdn.net/Joshua_HIT/article/details/72868459beadcount=NULL,#---Beadcount information for Green and Red Channal, need for filterBeads.(default = NULL) form ChAMP#---number of beads used to summarize the red and green channel from minfiautoimpute=TRUE,#---参见champ.process()filterDetP=TRUE,#---是否过滤detP较大的探针detPcut=0.01,#---设定的detP阈值，值越小越goodSampleCutoff=0.1,#---每个样本的failed detP value 占比超过SampleCutoff=0.1时该sample会被丢弃ProbeCutoff=0,#---在移除failed sample后，每个探针中failed detP value 占比超过ProbeCutoff时该探针会被丢弃filterBeads=TRUE,#---每个探针结合了序列的bead少于3个时会被设为NAbeadCutoff=0.05,#---设定的beadcount frequency阈值filterNoCG = TRUE,#---是否过滤非CpG位点的探针filterSNPs = TRUE,#---是否过滤SNP位点附近的探针filterMultiHit = TRUE,#---是否过滤映射到基因组多个位置的探针filterXY = TRUE,#---是否过滤检测性染色体上的位点的探针population = NULL,#---筛选特定的人群，如"AFR"(Africa),人群列表见http://www.internationalgenome.org/category/population/ fixOutlier = TRUE,#---是否修改异常值，当该参数为T时，若beta值小于0，则会用最小正数代替该值，若beta值大于1，则会用最大整数代替该值arraytype = "450K")
#返回一个列表，包含两个元素："beta"，"pd"
#新增3个变量："hm450.manifest.hg19"，"multi.hit"，"probe.features"#-----数据读取与数据过滤------------------------------
myLoad<-champ.load(testdir)
# Usage:
champ.load()参数见champ.import()和champ.filter()
# champ.load()==champ.import()+champ.filter()
#返回一个列表包含3个元素："beta","intensity","pd"
#新增4个变量
#champ.load()利用minfi包来读入数据

注：每一步处理都有过程的解析，见下图

注：在此使用champ.load()进行后续数据分析

读入数据后有6种数据过滤:

① filterDetP、detPcut：默认过滤掉<0.01的探针，使用detPcut来修改阈值

② filterBeads、beadCutoff：探针与DNA序列杂交，beadcounts就是与探针结合的reads数，若结合得较少，则认为该种探针信号不可靠，实际处理中，默认如果这个探针在5%的样本中，beadcount<3这个探针就会被过滤掉，使用beadCutoff来修改阈值

③ filterNoCG：过滤非CpG位点的探针，对甲基化的C而言，有CpG\CHG\CHH等甲基化的C，若不想过滤这些探针可以设置filterNoCG=FALSE

④ filterSNPs：过滤SNP附近的CpG位点的探针，SNP位点附近的CpG位点不稳定，设filterSNPs=FALSE不进行过滤

⑤ filterMultiHit：过滤非特异性的探针，Nordlund等研究人员使用bwa将探针序列与基因组进行比对，若唯一比对，说明探针特异性好，根据比对结果，去除比对到基因组上多个位置的探针，修改参数filterMultiHit=FALSE不进行过滤

⑥ filterXY：过滤性染色体上的探针，性染色体上的差异CpG是不准确的，无法解释这个差异是否是实验造成，若性别相同时，可以不用过滤，修改参数filterXY=FALSE不进行过滤

2.查看数据过滤后的情况

#-----可视化查看质量控制情况---------------------------
champ.QC()
# Usage:
champ.QC(beta = myLoad$beta,#---产生的β矩阵pheno=myLoad$pd$Sample_Group,#---样本分类信息mdsPlot=TRUE,densityPlot=TRUE,dendrogram=TRUE,PDFplot=TRUE,Rplot=TRUE,Feature.sel="None",#---绘制dendrogram时所使用的计算方法resultsDir="./CHAMP_QCimages/")
QC.GUI()
#QC.GUI(beta=myLoad$beta,pheno=myLoad$pd$Sample_Group,arraytype="450K")
#QC.GUI() 函数也可以画图，但是比较耗内存
#包括5张图:除了上述三种图还包括type-I&II densityPlot,top 1000 most variable CpG’s heatmap.
#type-I&II densityPlot图可以帮助查看两个探针的标准化状态
#Top variable CpGs’ heatmap将前1000个差异最大的位点和状态表示出来

图像解读：C为control组，T为tumor组

3.Type2 probe矫正

两种探针设计不同、分布区域不同，最主要的还是type2会动态衰减导致信号读取有差异，因而对 type2 探针的矫正很重要，对探针的矫正即消除两探针因技术差异带来的误差

#-----type Ⅱ探针的矫正-----#
myNorm<-champ.norm(plotBMIQ=T)
# Usage:
champ.norm(beta=myLoad$beta,rgSet=myLoad$rgSet,mset=myLoad$mset,resultsDir=".CHAMP_Normalization/",method="BMIQ",plotBMIQ=FALSE,arraytype="450K",cores=3)#---运行程序时使用的电脑核数，用于提高计算效率
#矫正beta值，矫正后的值存放在myNorm变量中
#reSet是Red和Green Channel，具体可参见minfi包：read.metharray()和read.metharray.exp()
#标准化的方法有BMIQ、PBC、SWAN、FunctionalNormalize，SWAN需要rgSet和mSet,F....需要rgSet,PBC和BMIQ仅需要beta值

矫正后的beta矩阵

每个样本都生成3个图(plotBMIQ=T,默认不会产生plot)关于探针矫正前后的拟合图（type1,type2,type1+type2）

图像解读：

4.变异组分检验

SVD plot:识别最显著的变异成分，这些成分可能是我们感兴趣的生物因素，也可能是技术差异带来的误差

如果从 .idat导入原始文件，设置champ.SVD()函数的RGEffect=TRUE ，芯片上18个内置的对照探针会作为协变量被纳入研究进行分析
champ.SVD()函数将把pd文件中的所有协变量(如,position,batch,well)和表型数据(如,gender,age)纳入进行分析。可以用cbind()函数将自己的协变量与myLoad$pd合并进行分析，你也可以自己写一个.txt文件用于存储这些表型数据，但第一列一定为pd文件的Sample_Name

对于分类变量和数字变量处理方法是不一样的。分类变量要转换成“factor” or “character”类型，数字变量转换成数字类型

#-----变异组分显著性检验-----#
champ.SVD()
# Usage:
champ.SVD(beta = myNorm,rgSet=NULL,pd=myLoad$pd,RGEffect=FALSE,#---If Green and Red color control probes would be calculatedPDFplot=TRUE,Rplot=TRUE,resultsDir="./CHAMP_SVDimages/")
#产生2个图像：Singular Value Decomposition Analysis (SVD)；Scree Plot

图像解读：

颜色越深，p-value越小变异越显著，表明变异成分越大

从图像可以看出变异主要来自样本间(control组和tumor组)，即我们关注的生物学因素，但变异显著性不大，经过过滤后得到的数据比较clear，适合后续分析。
如果发现在array或slide等有较显著的p-value，则需要检查实验设计以及使用其他normalization methods来进行后续的批次效应矫正

scree plot：各主成分可解释的变异占总变异的百分比

！！interesting：champ.QC()查看样本质量；champ.norm()查看探针质量；champ.SVD()主要查看混杂因素

5.批次效应矫正

champ.runCombat() 函数自动把Sample_Group作为协变量矫正，现在又加入了另一个参数variablename用来加入自定义的协变量进行矫正

#------批次效应矫正------#
mycombat<-champ.runCombat()
# Usage:
champ.runCombat(beta=myNorm,pd=myLoad$pd,variablename="Sample_Group",#---要进行矫正的协变量，参考sva包的ComBat函数batchname=c("Slide"),#---批次，来自pd文件logitTrans=TRUE)#---对beta是否进行logit转化，利用原始beta进行矫正时设为T，利用M值进行计算时设为F

ComBat如果直接用beta值进行矫正，输出值可能不在0-1之间，所以在计算之前会进行一个logit转化成M-value，若用M值进行矫正则没必要进行转化，因此设置logitTrans=F

beta与M-value的关系：

有时候批次效应和变异会混杂在一起，如果矫正了批次效应，变异也会消失，因此在接下来的分析不使用批次效应矫正后的数据

6.甲基化分析

DMP代表Differential Methylation Probe（差异甲基化CpG位点）
DMR代表Differential Methylation Region（差异化CpG区域）
Block代表Differential Methylation Block（更大范围的差异化region区域）

#甲基化分析的函数可参考limma包
# 差异甲基化位点分析
myDMP<-champ.DMP()
# Usage:
champ.DMP(beta = myNorm,pheno = myLoad$pd$Sample_Group,compare.group = NULL,#---当group>2时设定该参数为一个列表，其中每个元素为两个Sample_Group，进而指定要进行比较的group，否则每两个group之间都会进行比较adjPVal = 0.05,#---设定的p-value检验水平adjust.method = "BH",#---p-value矫正，使用的method参考函数p.adjust()arraytype = "450K")
# 直接通过基于类别group的探针比较计算得到DMP# 可视化查看分布
DMP.GUI(DMP=myDMP[[1]],beta=myNorm,pheno=myLoad$pd$Sample_Group)
#使用GUI需要shiny包

生成的myDMP
myDMP的元素C_to_T的含义：C(control组) vs T(tumor组)

logFC = deltaBeta
XXX_AVG = mean beta of the group XXX
MAPINFO = position of CpG

# 差异甲基化区域分析
myDMR<-champ.DMR()
# Usage:
champ.DMR(beta=myNorm,pheno=myLoad$pd$Sample_Group,compare.group=NULL,arraytype="450K",method = "Bumphunter",minProbes=7,#---DMR包含很多CpG位点,每个位点都对应一个探针，minProbes设置每个DMR包含的最小探针数adjPvalDmr=0.05,#---DMR的p-value是用Stouffler's method计算得到，再进行矫正后得到adjPvalDmrcores=3,## following parameters are specifically for Bumphunter method.maxGap=300,cutoff=NULL,pickCutoff=TRUE,smooth=TRUE,smoothFunction=loessByCluster,useWeights=FALSE,permutations=NULL,B=250,nullMethod="bootstrap",## following parameters are specifically for probe ProbeLasso method.meanLassoRadius=375,minDmrSep=1000,minDmrSize=50,adjPvalProbe=0.05,Rplot=T,PDFplot=T,resultsDir="./CHAMP_ProbeLasso/",## following parameters are specifically for DMRcate method.rmSNPCH=T,fdr=0.05,dist=2,mafcut=0.05,lambda=1000,C=2)# 可视化查看分布
DMR.GUI(DMR=myDMR)

生成的myDMR by Bumphunter

For Bumphunter algotithm, it will replying on no previous output. Firstly, Bumphunter all cluster all probes into small clusters. The size and distances between CpGs could be speciﬁed by users. Then Bumphunter will ﬁlter all clusters contain too few CpGs in it, the default threshold is 7, but can be speciﬁed by user as well. After ﬁltering, Bumphunter will select candidate DMRs based on probes’ diﬀerential t value between two groups of samples, and their location in one cluster. Finally random permutation technic will return signiﬁcance of each candidate DMR and p value will be returned. For more information, user may turn to bumphunter package.
The result of bumphunter algotithm is a dataframe contain all detected DMRs, with their length, clusters, number of CpGs.e.g annotated. Also, Bumphunter algotithm will return the annotation of all DMR contained CpGs

Bumphunter可能产生的plot（使用当前2020年12月的champ.DMR()不会产生plot）

图像解读：
每列代表一个probe/CpG位点，因为使用的数据有八个sample，所以每一列有八个点
不同的颜色表示样本所属不同group，红色为tumor group，黑色为control group
两条线分别为两个group的mean beta
蓝色轮廓的点所在的整个区域为DMR（甲基化稳定 by 两线间隔基本不变）

生成的myDMR2 by ProbeLasso

For ProbeLasso algotithm, because the ProbeLasso(lasso is defined as dynamic window) DMR Hunter takes into account probe spacing, the ﬁrst step to identify DMRs is to calculate a dataset-speciﬁc record of probe spacing.As such, each dataset is bestowed a unique set of probes, which underscores the need for an experiment-speciﬁc catalogue of probe-spacing. Once the dataset has been deﬁned, the nearest neighbouring probe is calculated for all available probes; this data is then partitioned into one of 28 diﬀerent categories comprising information on the genomic feature (1st exon, 30UTR, 50UTR, body, intergenic region, TSS200, or TSS1500) and that features relation (if any) to a nearby CpG island (island, shore, shelf, or not associated).The user then speciﬁes a maximum (or minimum) lasso size (bp) and the Probe Lasso DMR Hunter determines which feature/CGI category fulﬁls this criterion ﬁrst and what quantile of the relevant feature/CGI category it corresponds to. This quantile is applied to all feature/CGI categories to deﬁne feature/CGI category-speciﬁc lasso sizes. Next, operating solely on the signiﬁcant probes in the dataset, an appropriatelysized lasso is thrown around each probe; if the lasso captures a user-speciﬁed number of signiﬁcant probes (including itself), that probe (and the probes captured by its lasso) is set aside, in essence, a ”mini-DMR”. Next, ProbeLasso DMR Hunter attempts to close-up gaps between neighbouring miniDMRs whose lasso boundaries are either overlapping or less than a userspeciﬁed distance apart (e.g., 1000bp), eﬀectively deﬁning a ”ﬁnal DMR”.
Finally, ProbeLasso DMR Hunter pulls back all probes within the DMR coordinates and returns them as data frame containing genomic annotation and association statistics of each MVP. Additionally, a p-value for each DMR is calculated using Stouﬄer’s method. This method combines the p-values of individual probes within a DMR by weighting them according to the underlying correlation structure of methylation scores between probes; pvalues of probes whose methylation scores are correlated are down-weighted, while independent probes are not penalised.

ProbeLasso可能产生的plot（使用当前2020年12月的champ.DMR()不会产生box plot和cumulative quantile plot，只会生成最后一个dot plot，但格式会有一点点不同）

a box plot of probe spacing for the input data as a function of genomic feature/CGI relation

Fig. A study speciﬁc plot showing the number of probes found in each feature. Note the diﬀerent spacing on the 450k or EPIC array between neighbouring probes for diﬀerent features (1st exon, 30UTR, 50UTR, gene body, intergenic region, transcription start sites) and their relation to the nearest CpG island (in the island, in a shore, in a shelf, not associated). Intergenic regions, 30UTRs and gene bodies with no CGI relation are very sparsely spaced; whereas probes in the 1st exon and TSSs are very closely spaced.

a cumulative quantile plot illustrating how the user-speciﬁed parameters aﬀect the sizes of all windows employed for DMR-calling

Fig. Deﬁning the lasso for each feature type. This illustrates when someone has speciﬁed a minimum lasso of 10bp and how this sets the quantile for all feature/CGI relation-speciﬁc lassos to ˜48%. These lassos are then cast out (centred on each probe) to capture a minimum number of signiﬁcant probes. Probes that fulﬁl this are set aside.

a dot plot with size-scaled dots to further illustrate the resulting window sizes for each genomic feature that follows from the user-speciﬁed parameters

Fig. This image shows the radius size for each feature

使用ChAMP中的数据产生plot如下

#Differential Methylation Block
champ.Block()
# Usage:
champ.Block(beta=myNorm,pheno=myLoad$pd$Sample_Group,arraytype="450K",maxClusterGap=250000,B=500,bpSpan=250000,minNum=10,cores=3)

除此之外，还可以进行细胞异质性处理champ.refbase()、CNA分析champ.CNA()、GSEA富集分析champ.GSEA()，具体使用参考帮助文档

日常偷懒：（如果只是简单分析甲基化情况可以仅进行以下步骤，函数都使用默认）
myLoad<-champ.load(路径名)
myNorm<-champ.norm()
myDMP<-champ.DMP()
myDMR<-champ.DMR()
注意，预处理步骤的变量名必须为myLoad myNorm，这样就可以不用在后续函数调用中添加参数beta=…,pd=…之类的，直接函数名后跟()就可以执行