甲基化芯片入门学习-ChAMP包(二)

DNA甲基化芯片分析有不少R包实现,如:minfi、lumi以及ChAMP等,我只粗略看过minfi和ChAMP,发现ChAMP的功能更加齐全以及使用也较为简单,并且其也集成了minfi包的部分功能,所以下面以ChAMP包作为学习对象

ChAMP包的安装
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("ChAMP")

如果是墙内,最好加上options选择中科大镜像,因为有好多依赖包需要安装

ChAMP包也依赖了minfi包的一些功能,还似乎安装了好多配套注释文件的包,省的以后需要时再去找了,后续需要的话直接加载即可了

读入数据

这里有个问题,由于ChAMP包以及minfi包读入文件时,除了需要.idat文件外,还需要sample_sheet.csv文件,里面记录了样本信息、样本分组信息、每个GSM文件对应哪个样本等信息。但是之前给的那个测试数据我在GEO上下载的RAW文件里没有该样本信息文件。后来也尝试了选择其他GSE号的数据,也未发现有该文件,最终暂时放弃尝试了,决定选择ChAMP自带的450K芯片的测试数据(理论上illumina下机后除了有.idat文件外,还有sample_sheet.csv文件的)

先下载测试数据

testDir=system.file("extdata",package="ChAMPdata")
myLoad <- champ.load(testDir,arraytype="450K")

数据可以在R包library/ChAMPdata/extdata文件夹下找到,里面除了.idat文件外,还有一个lung_test_set.csv样本信息文件;所以我猜想,尽管有些GSE数据没有样本信息.csv文件,我们可以根据lung_test_set.csv格式模仿写一个(Header部分可有可无,但Data部分几个关键的列必须要有),前提当然要知道样本分组信息(如果不知道的话,就编一个吧。。。),模板如下图所示:

sample_sheet_csv

除了上述文件外,还有一个文件比较重要的是450K注释文件,如果是GEO下载的数据化,就是XXXHumanMethylation450XXX.csv文件,里面有详细的每个探针的数据,比如探针的ID、不同探针Type的序列,在染色体上的定位信息以及注释信息等

数据过滤

准备好上述文件后,接下来则是读入数据,使用ChAMP包的champ.load函数,按照其说法,champ.load函数其实做了两个步骤,第一是读入数据(champ.import函数),第二是过滤数据(champ.filter函数),所以我们需要了解champ.load以什么标准过滤:

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

以上标准均是champ.load的默认参数,如果是850K芯片还需加上arraytype=”850K”,因为默认是450K芯片。。

myLoad <- champ.load("E:/R-3.4.1/library/ChAMPdata/extdata")

读入数据时会出现很多信息,其中可以看到有613个control probes,这个探针主要用于计算芯片数据的可信度Detect P value,具体解释可参考DNA Methylation中的Detect P value解析

Fetching NEGATIVE ControlProbe. Totally, there are 613 control probes in Annotation. Your data set contains 613 control probes.

还有不少过滤过程的信息,基本上是根据你设定的过滤参数来的

做完上述步骤,ChAMP包还提供了CpG.GUI()函数用于展示distribution of CpGs on chromosome, CpG island, TSS reagions. e.g,其实是几张简单展示myLoad$beta分布的图片,但是界面很友好,也难怪依赖包里要shiny包了

CpG.GUI(CpG = rownames(myLoad$beta))

对于

数据质控QC

ChAMP包不愧是集成度很高的DNA甲基化分析的R包,啥功能都有。。champ.QC()QC.GUI函数用于检查看看过滤后的数据是否符合要求,能否进行下一步的分析。

champ.QC()函数结果三张图,如有mdsPlot (Multidimensional Scaling Plot),主要看看不同分组样本是否分开;densityPlot,每个样本的beta值的分布图,主要看看有无异常的样本;dendrogram,样本的聚类图

champ.QC(beta = myLoad$beta, pheno = myLoad$pd$Sample_Name)

上述代码我的电脑会报错:LAPACK routines cannot be loaded,网上查了下,似乎是电脑的原因,具体图片样式可看R包说明文档,主要也是根据myLoad$beta数据以及myLoad$pd$Sample_Name样本信息来绘制图片的

QC.GUI(beta=myLoad$beta)

这个也是报同样的错误,但ChAMP包的结果可视化做的确实蛮好的

数据标准化

由于Illumina beadarrays是由两套探针,使用不同的杂交方法,所以如果上述QC结果中发现type-I探针和type-II探针有较大偏差,则需进行校正;简单的说就是消除测序芯片上的两种测序技术之间的差异。ChAMP包标准化方法有BMIQ、SWAN、PBC以及FunctionalNormliazation,我这里就先选择默认的BMIQ,后续如果有需求再看看其他标准化方法

myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)
SVD plot

SVD(singular value decomposition) 只为检测变异组分的显著性,主要用于查看变异的主要成分是生物学处理等影响的,还是技术因素所造成的,如果是后者则需要进行后续的批次校正。从下面代码也可看出其主要根据myNorm(校正后的探针的beta)和pd(样本信息)来计算变异程度

champ.SVD(beta=myNorm,pd=myLoad$pd)
批次效应校正

如果上面SVD Plot显示红色块并不在Sample_Group行,而是较多的在Array行,那么则需要对这次数据进行批次校正,使用champ.runCombat函数,batchname则是根据上述SVD plot出现的红块位置而定,决定哪个batch factor需要被校正;这个过程比较耗时,做完后最好再用champ.SVD检查下

myCombat <- champ.runCombat(beta=myNorm,pd=myLoad$pd,batchname=c("Slide"))
甲基化探针差异分析

比较常见的差异甲基化分析是DMP(Differential Methylation Probe),用于找出差异的甲基化位点,然后ChAMP包已经将分析过程(主要还是基于Limma包,这个包是专门用于分析芯片差异表达的包);然后根据你的分组信息,获得差异甲基化位点;最后用DMP.GUI查看下结果

myDMP <- champ.DMP(beta = myNorm,pheno=myLoad$pd$Sample_Group)
DMP.GUI(DMP=myDMP[[1]],beta=myNorm,pheno=myLoad$pd$Sample_Group)

除了可以分析差异甲基化位点外,还可以分析DMR(Differential Methylation Regions),用于找出差异甲基化片段?用的函数也很好记champ.DMR(),该函数现在只支持两组样本进行差异分析

myDMR <- champ.DMR(beta=myNorm,pheno=myLoad$pd$Sample_Group,method="Bumphunter")

这里需要注意的是method的选择,主要有Bumphunter, ProbeLasso and DMRcate三者可以选择,区分的话

  1. Bumphunter,groups all probes into small clusters (or regions), then applies random permutation method to estimate candidate DMRs,不用依赖上述差异位点分析的结果
  2. ProbeLasso,the final data frame is distilled from a limma output of probes and their association statistics,必须要有champ.DMP()的结果作为输入
  3. DMRcate,a new method recently incoporated into ChAMP,pass the square of the moderated t statistic calculated on each 450K probe to DMR-finding function(这个还是看文档吧)

除了DMR外还有一个DMB(Differential Methylation Blocks),也是用来识别大片段的差异甲基化

myBlock <- champ.Block(beta=myNorm,pheno=myLoad$pd$Sample_Group,arraytype="450K")
其他分析类型

除了上述比较常规差异甲基化分析的外,还有GSEA富集分析(Gene Set Enrichment Analysis)、网络分析(Differential Methylated Interaction Hotspots)和拷贝数变异(Copy Number Variation),ChAMP包都集成到一个个函数里了,用起来确实比较方便

Summary

本来这个笔记应该是实际数据操作的过程,自从数据质控QC开始报错后,后面只要用到LAPACK的函数均报错了,找了好久都没找到在windows下的解决办法(只能这过段时间在服务器上跑跑看了)。所以这笔记变成了ChAMP包文档笔记了,顺着文档粗略的扫了一遍,没细读,以后需要的时候再补了

最后上述的BUG还是解决了,原来报错的是因为LAPACK routines没有被加载,但是真正的原因则是因为 maximal number of DLLs reached,出现这个的原因是ChAMP包加载很多DLLs,而R默认的最大DLLs是100,如果超过100个了,那么就不让继续加载新的DLL。

因此解决办法是更改R的默认设置,让R_MAX_NUM_DLLS变量的值设为100+,之前查的解决办法都是在linux或者mac上的,今天终于让我找到在windows下的一种解决办法,就是在设置电脑环境变量中新增一个变量R_MAX_NUM_DLLS,并设置值为200,然后重新启动R,然后上述的那些ChAMP包都能正常工作了。可以好好看看整个流程的结果了

参考资料:

The Chip Analysis Methylation Pipeline
计算表观学公众号:Illumina 850K(EPIC)甲基化芯片及其分析工具ChAMP
生信技能树公众号:850K甲基化芯片数据的分析
Shiny和Plotly实现可交互DNA甲基化分析包ChAMP
DNA Methylation中的Detect P value解析
读取DNA甲基化IDAT文件
集识别差异甲基化和拷贝数变异于一体的R包——ChAMP
bioconductor的minfi处理甲基化数据
Illumina HumanMethylation450 BeadChip (甲基化450k芯片) 预处理初探
450K甲基化芯片数据处理传送门

本文出自于http://www.bioinfo-scrounger.com转载请注明出处