五分钟学会甲基化芯片处理，快上车！！！

甲基化分析系列之CHAMP包的使用（二）

大家好，我是风间琉璃，这里是每周四的甲基化生信分析专栏。接着上一期我们讲解了CHAMP的数据处理流程以及数据导入，今天我们将继续讲解之后的分析步骤。
传送到上期戳这里：
手把手教你甲基化分析——甲基化CHAMP包的使用（一）

步骤回顾：

1.数据导入（champ.import）,包括甲基化位点筛选（champ.filter）和缺失值的插补（champ.impute）
2.数据质控（champ.QC）
3.数据标准化（champ.norm）
4.奇异值分析了解批次效应（champ.SVD）
5.批次校正（champ.runCombat）和细胞类型异质性校正（champ.refbase）
6.甲基化差异位点分析（champ.DMP）和差异区域分析（champ.DMR）
7.甲基化位点/区域的通路富集分析（champ.GSEA）

第二步：数据质控

数据质控在CHAMP中很容易实现，只需要一个命令就可达成。

#####数据质控#####champ.QC()###QC质控###可视化#CpG.GUI(CpG=rownames(myLoad$beta),arraytype="450K")#QC.GUI(beta=myLoad$beta,arraytype="450K")

在champ.QC中，我们可以得到三个图，分别是MDS图、密度图和树状图。

我们发现在MDS图中，对照组（C）与肿瘤组（T）是明显分开的。四个对照组样本聚在一起，而肿瘤组样本是分得很开的，这可能与肿瘤具有极高的异质性相关。

在密度图发现对照组和肿瘤组的beta值在0.1、0.9是最富集的，这和我们的生物学知识背景也相同。

树状图我们可以看到，对照组四个都聚在一起，但是肿瘤组的样本间差异就很大。

后面的两行被我注释的代码就是CHAMP包特有工具，使用shiny构建的GUI，我们可以根据用鼠标点点点从而得到自己想看到的图。进行初步的分析之后，这八个样本认为是都是符合质量控制的，那我们就进入下一步——标准化。

第三步：标准化

其实在CHAMP包中的标准化流程也非常简单，只需要使用一个命令champ.norm即可以得到标准化后的甲基化数据。

#######标准化，这里使用的是β，并允许有5个现成一起计算###默认使用的是BMIQ方法，还有PBC、SWAN等方法myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)#标准化后得到新的数据##在看看标准化后得champ.QC(myNorm)##还是GUI的方式进行可视化#QC.GUI(myNorm)

我们发现好像校正后的图片并没有什么变化，可能是我们的样本本身就比较好吧。接下来就到了我们的第四步。

第四步：奇异值分析了解批次

这一步也只需要一个命令。

######SVD plotchamp.SVD(beta=myNorm,pd=myLoad$pd)

## Sample_Group Sample_Well Slide Array## [1,] 0.02092134 0.3286227 0.2227875 0.4365574## [2,] 0.14891467 0.3286227 0.1097210 0.2903377## [3,] 0.77282999 0.5545449 0.6998607 0.7256622

虽然代码很简单，但是理解起来可能并不容易。这里的结果提供了两张图。

第一张图展示的是碎石图：

这一张图是进行SVD反卷积（SVD deconvolution）之后得到的每一个成分（componnent）能够解释变异的百分比，我们发现成分一很明显非常重要，因为它占据总百分比的60%以上。

第二张图是相关性的热图：

我们进一步看我们的每一个成分究竟与我们提供的pdata文件中的什么表型有关。我们发现我们的PC主要还是与分组有关，而和其他的几乎没有关系，再一次证明了我们的分组情况非常好。我们看一看下面这张：

发现其实并没有这个数据分析的结果这么好了对吧？主成分和许多表型有关，包括种族、性别、上级的板次等等。经过奇异值分析之后，如果我们很明显的看到具有批次，也就是和我们上样的板次（plate）有关，这时候我们需要进行批次校正了，这一点其实也是CHAMP非常好的一点。从另外的一个角度讲，如果我们需要整合多个甲基化芯片数据，使用CHAMP的SVD分析+批次校正不就可以达到多样本整合分析的目的了吗？

第五步：批次矫正

在这里因为我们的数据并不存在明显的批次，但是为了演示，我们还是校正了上样的板次（plate）。在进行校正之前，我们还需要将数据的pdata的批次转化为因子格式。

####combat####在此之前需要把我们得批次改成因子格式myLoad$pd$Slide=as.factor(myLoad$pd$Slide)###很耗内存耐心等一会儿myCombat <- champ.runCombat(beta=myNorm,pd=myLoad$pd, variablename="Sample_Group",batchname="Slide")## ~Sample_Group## <environment: 0x000000005cc63698>## Found 1 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.

进行完上述步骤之后，就到了我们CHAMP分析的重头戏了。

第六步：差异分析

######差异分析#####myDMP <- champ.DMP(beta = myCombat,pheno=myLoad$pd$Sample_Group)## Contrasts## Levels pT-pC## pC -1## pT 1#查看我们的分析结果head(myDMP[[1]])## logFC AveExpr t P.Value adj.P.Val B## cg26970841 0.5183941 0.2937171 58.52077 1.516600e-10 3.340464e-05 12.11173## cg06649610 -0.2799447 0.6084112 -56.64793 1.896502e-10 3.340464e-05 12.04163## cg20097929 -0.3245835 0.7689923 -50.98097 3.912603e-10 3.340464e-05 11.79103## cg16448058 0.4092639 0.4251169 50.43682 4.211893e-10 3.340464e-05 11.76343## cg04461197 -0.2983230 0.6822998 -47.94310 5.966500e-10 3.340464e-05 11.62759## cg06356785 -0.4377791 0.6456765 -47.41555 6.437495e-10 3.340464e-05 11.59674## C_AVG T_AVG deltaBeta CHR MAPINFO Strand Type gene## cg26970841 0.03452007 0.5529142 0.5183941 16 85932666 R I IRF8## cg06649610 0.74838353 0.4684389 -0.2799447 11 65405268 F II SIPA1## cg20097929 0.93128403 0.6067005 -0.3245835 10 83464779 R II ## cg16448058 0.22048495 0.6297489 0.4092639 10 102880841 R I TLX1NB## cg04461197 0.83146136 0.5331383 -0.2983230 3 172635656 R II SPATA16## cg06356785 0.86456604 0.4267869 -0.4377791 21 45844767 R II TRPM2## feature cgi feat.cgi UCSC_CpG_Islands_Name DHS## cg26970841 TSS200 island TSS200-island chr16:85932121-85932942 TRUE## cg06649610 TSS1500 shelf TSS1500-shelf chr11:65408344-65408631 NA## cg20097929 IGR opensea IGR-opensea NA## cg16448058 5'UTR shelf 5'UTR-shelf chr10:102882977-102883551 NA## cg04461197 Body opensea Body-opensea NA## cg06356785 Body shore Body-shore chr21:45844992-45845693 NA## Enhancer Phantom Probe_SNPs Probe_SNPs_10## cg26970841 NA ## cg06649610 TRUE rs931127 ## cg20097929 TRUE ## cg16448058 NA ## cg04461197 TRUE ## cg06356785 NA####还是GUI的可视化#DMP.GUI(DMP=myDMP[[1]],beta=myCombat,pheno=myLoad$pd$Sample_Group)

# myDMP是一个列表格式, 其中myDMP[[1]], myDMP[[2]], myDMP[[3]]等列表中的数据框，因为我们只分析两组，所以就只有一个数据框存储差异甲基化位点信息

####羟甲基化（Hydroxymethylation Analysis）##羟甲基化一般用不到，所以注释掉#myDMP <- champ.DMP(beta=myCombat, # pheno=myLoad$pd$Sample_Group, compare.group=c("oxBS", "BS"))#提取羟甲基化位点信息#hmc <- myDMP[[1]][myDMP[[1]]$deltaBeta>0,]

进行完甲基化位点的差异分析之后，常规到了我们的下一步，甲基化区域（region）甲基化模块（block）差异分析

####差异甲基化区域，DMRmyDMR <- champ.DMR(beta=myCombat,pheno=myLoad$pd$Sample_Group,method="Bumphunter")#查看差异甲基化区域head(myDMR$BumphunterDMR)## seqnames start end width strand value area cluster## DMR_1 chr12 115134148 115136308 2160 * 3.050049 152.50247 51386## DMR_2 chr7 27224343 27226329 1986 * 2.664508 85.26425 179533## DMR_3 chr6 32115964 32117401 1437 * 2.541977 83.88523 167665## DMR_4 chr7 45961289 45962236 947 * 3.322501 79.74003 180997## DMR_5 chr17 46655289 46656093 804 * 4.971320 79.54112 88509## DMR_6 chr6 30094947 30095802 855 * 2.644733 74.05253 166933## indexStart indexEnd L clusterL p.value fwer p.valueArea fwerArea## DMR_1 24376 24425 50 50 0 0 4.998001e-05 0.012## DMR_2 84217 84248 32 35 0 0 9.996002e-05 0.016## DMR_3 77405 77437 33 36 0 0 9.996002e-05 0.016## DMR_4 84867 84890 24 29 0 0 1.499400e-04 0.020## DMR_5 38745 38760 16 16 0 0 1.499400e-04 0.020## DMR_6 74017 74044 28 28 0 0 2.165800e-04 0.028##还是GUI可视化#DMR.GUI(DMR=myDMR)#这一步使用DMR.GUI可视化的时候可能会非常耗时。####差异甲基化模块，DMBlock###速度比较慢，多点耐心# myBlock <- champ.Block(beta=myCombat,pheno=myLoad$pd$Sample_Group,# arraytype="450K")# #查看差异甲基化模块# head(myBlock$Block)##还是GUI可视化#Block.GUI(Block=myBlock,beta=myCombat,pheno=myLoad$pd$Sample_Group,# runDMP=TRUE,compare.group=NULL,arraytype="450K")

最后，甲基化同样可以进行功能富集分析，我们来试试吧。

第七步：差异甲基化位点/区域功能富集分析

###GSEA分析，输入差异分析后的输出结果，包括DMP和DMR(DMR不是必须的)myGSEA <- champ.GSEA(beta=myCombat,DMP=myDMP[[1]], DMR=myDMR, arraytype="450K",adjPval=0.05, method="fisher")#查看DMP和DMRhead(myGSEA$DMP)## Gene_List## MEISSNER_BRAIN_HCP_WITH_H3K4ME3_AND_H3K27ME3 MEISSNER_BRAIN_HCP_WITH_H3K4ME3_AND_H3K27ME3## BENPORATH_ES_WITH_H3K27ME3 BENPORATH_ES_WITH_H3K27ME3## GGGAGGRR_V$MAZ_Q6 GGGAGGRR_V$MAZ_Q6## CAGGTG_V$E12_Q6 CAGGTG_V$E12_Q6## BENPORATH_SUZ12_TARGETS BENPORATH_SUZ12_TARGETS## AACTTT_UNKNOWN AACTTT_UNKNOWN## nList nRep fRep nOVLAP## MEISSNER_BRAIN_HCP_WITH_H3K4ME3_AND_H3K27ME3 1069 1031 0.9644528 926## BENPORATH_ES_WITH_H3K27ME3 1118 1020 0.9123435 915## GGGAGGRR_V$MAZ_Q6 2274 2046 0.8997361 1753## CAGGTG_V$E12_Q6 2485 2266 0.9118712 1929## BENPORATH_SUZ12_TARGETS 1038 900 0.8670520 803## AACTTT_UNKNOWN 1890 1676 0.8867725 1444## OR P.value adjPval## MEISSNER_BRAIN_HCP_WITH_H3K4ME3_AND_H3K27ME3 2.621172 2.721428e-25 2.265861e-21## BENPORATH_ES_WITH_H3K27ME3 2.588050 1.462291e-24 6.087515e-21## GGGAGGRR_V$MAZ_Q6 1.803739 1.186352e-21 3.292522e-18## CAGGTG_V$E12_Q6 1.729581 6.452669e-21 1.343123e-17## BENPORATH_SUZ12_TARGETS 2.444107 5.146664e-20 8.570224e-17## AACTTT_UNKNOWN 1.861984 7.465872e-20 1.036014e-16###展示DMR的功能富集分析结果#head(myGSEA$DMR)####Empirical Bayes GSEA method###这一步非常漫长，不要轻易尝试#myebayGSEA <- champ.ebGSEA(beta=myCombat,pheno=myLoad$pd$Sample_Group,# arraytype="450K")# head(myebayGSEA[[1]])# head(myebayGSEA[[2]])

彩蛋部分，我们还可以通过甲基化的情况分析基因拷贝数的程度变化，是不是没想到呢？

我们来看一看：

#####Copy Number Variation###仍然比较慢。。。。myCNA <- champ.CNA(intensity=myLoad$beta,pheno=myLoad$pd$Sample_Group)# ##默认保存为pdf图片sample_CNA=myCNA[[1]]####这个可以使用DNAcopy包进行可视化，每个样本可以进行绘图####同样可以导出seg文件，并进一步在IGV软件中进行可视化# ###看看一下数据长什么样head(sample_CNA[[1]])## ID chrom loc.start loc.end num.mark seg.mean## 1 C1.qn 1 69591 1334575 920 -0.1060## 2 C1.qn 1 1334791 1342720 17 -1.4858## 3 C1.qn 1 1342795 2322553 1092 0.0133## 4 C1.qn 1 2322827 2323870 9 -1.5377## 5 C1.qn 1 2324205 2344183 54 0.2686## 6 C1.qn 1 2344280 2344699 7 -1.7777group_CNA=myCNA[[2]]head(group_CNA[[1]])## ID chrom loc.start loc.end num.mark seg.mean## 1 C1.C.qn 1 69591 1406845 1011 -0.1272## 2 C1.C.qn 1 1406955 1407183 4 -2.7196## 3 C1.C.qn 1 1408104 2322553 1014 0.0308## 4 C1.C.qn 1 2322827 2323870 9 -1.5377## 5 C1.C.qn 1 2324205 2344183 54 0.2686## 6 C1.C.qn 1 2344280 2344699 7 -1.7777

我们看看基因拷贝数的相关图片。

我们先看看对照组的样本的拷贝数情况（横轴是染色体区域，纵轴是拷贝数情况）：