五分钟学会甲基化芯片处理,快上车!!!
甲基化分析系列之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 plot
champ.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)差异分析
####差异甲基化区域,DMR
myDMR <- 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和DMR
head(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.7777
group_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
我们看看基因拷贝数的相关图片。
我们先看看对照组的样本的拷贝数情况(横轴是染色体区域,纵轴是拷贝数情况):
我们再看看肿瘤组的样本拷贝数情况:
大家再运行完上述流程之后,所在的文件夹自动就会生成几个文件夹,里面保存着我们分析流程中自动产生的文件。如图:
使用AI拼一拼图就成为我们文章中的一个Figure啦。
好啦,我是风间琉璃,我们下期再聊聊另外一个包——methylationArrayAnalysis”。
回复“风间琉璃甲基化04”得到本次分享的代码。
—
END—
撰文丨风间琉璃
排版丨四金兄
值班 | 小雪球
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
最新评论
推荐文章
作者最新文章
你可能感兴趣的文章
Copyright Disclaimer: The copyright of contents (including texts, images, videos and audios) posted above belong to the User who shared or the third-party website which the User shared from. If you found your copyright have been infringed, please send a DMCA takedown notice to [email protected]. For more detail of the source, please click on the button "Read Original Post" below. For other communications, please send to [email protected].
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。