甲基化芯片入门学习-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转载请注明出处

文章目录
  1. 1. ChAMP包的安装
  2. 2. 读入数据
  3. 3. 数据过滤
  4. 4. 数据质控QC
  5. 5. 数据标准化
  6. 6. SVD plot
  7. 7. 批次效应校正
  8. 8. 甲基化探针差异分析
  9. 9. 其他分析类型
  10. 10. Summary
|