探索生信之美，解构每一段代码的故事
更多生信干货教程
私信我回复“数据库”即可免费领取哦
大家好呀，我是风间琉璃，前面我们花了许多期的时间介绍甲基化分析入门的教程，包括minfi包、CHAMP包的使用方法。大家可能觉得甲基化分析也就这样了吧，谬矣！ 今天我带大家甲基化分析登堂入室之分析——ELMER包（甲基化+转录因子梦幻联动），我将用2-3期的内容进行介绍，并在最后一期用一篇高分文献进行实例讲解，让大家看一下这一个包怎们在顶刊文献中进行应用的 。小声地说一句，如果是肿瘤的学员有福啦，本次分析全程使用TCGA数据进行分析，所以你们直接上手，毫无阻碍。另外非肿瘤的小伙伴也不要失望，我也会给大家介绍不是TCGA数据如何进行构建对象以及下游分析。好的，那我们废话不多说，正式开始吧。
〇、背景知识
首先我们需要补充和一些基础的甲基化与转录因子调控基因转录知识。DNA链上的CG位点能够被甲基化酶的修饰下在胞嘧啶上添加甲基从而发生甲基化，如果一DN*片A**段富含CG序列的密度比平均密度高1020倍，GC含量大于50%，长度大于200bp的区域，称为CpG岛（CpG island）。如果某一区域高度甲基化,一般认为在顺式调控元件部位（包括启动子、增强子等部位）甲基化，将影响DNA的结构，从而阻遏该部位基因的转录。
转录因子（transcription factor，TF）是与基因5’上游序列结合从而影响基因转录的蛋白。而DNA双链上与转录因子结合的位点称之为转录因子结合位点（transcription factor binding site, TFBS）。这个位点一般是实验分析出来的，但是有时候TFBS的碱基序列并不会一直不变，因此引入了转录因子结合基序的概念（transcription factor binding motif, TFBM）代表一个TF的结合特异性，通常通过汇总一系列结合位点的保守和可变位点而来。
〇、设计思路
背景知识介绍完毕，我们再来看看ELMER包的设计者的设计这个包的思路。由于甲基化能够调控基因表达，因此筛选出差异甲基化探针的碱基序列，从而找到富集的motif（转录因子结合基序），从而进一步反推出与这些motif相互作用的转录因子，从而构建出TF—甲基化位点——靶基因网络。
这就是ELMER包的分析思路。这一点我们通过ELMER的英文全称也可以看出，ELMER（Enhancer Linking by Methylation/Expression Relationships），顾名思义，增强子相关甲基化和表达相关网络。
说完思路我们再说一说操作流程。基于同一批样本基因的甲基化和基因的表达水平从而重构建出转录因子-基因调控网络。ELMER通过将处于顺式调控模块（通常是增强子）上甲基化位点水平变化从而将上游的主要调控转录因子和下游的靶基因进行联系。ELMER通过构建MAE（S4）对象，构建出一共1639个转录因子与潜在结合motif对应的靶基因调控网络。
ELMER的主要分析流程包含：
1.识别HM450k或者EPIC甲基化芯片的远端探针（distal probe）
2.识别正常vs.对照组间甲基化水平差异
3.识别差异甲基化探针的靶基因
4.确定具有差异甲基化差异探针以及靶基因相关探针富集的motif
5.确认与富集motif上甲基化相关的调控转录因子
解释一下，distal probe是距离转录起始位点位点上游大于2kb的甲基化探针，默认这一区域的甲基化探针区域为增强子区域。
好啦，前言我们就说这么多啦，我们直接上手吧~
〇、安装包
# #首先安装包
# if (!requireNamespace("BiocManager", quietly=TRUE))
# install.packages("BiocManager")
# BiocManager::install("ELMER")
# #安装依赖包
# BiocManager::install("sesameData")
# BiocManager::install("GenomicInteractions")
# #*载下**ELMER的数据
# devtools::install_github(repo = "tiagochst/ELMER.data")
#加载包
library(MultiAssayExperiment)
library(ELMER)
## Loading required package: ELMER.data
library(ELMER.data)
library(sesameData)
一、数据导入
首先我们需要明确我们输入的数据，由于我们ELMER是根据基因的转录组数据和甲基化数据共同推断出转录因子-靶基因调控网络，并且ELMER需要通过输入的转录组和甲基化数据从而构建出Multi Assay Experiment对象，从而进行下游分析。所以ELMER需要以下经过处理后的数据
1.以探针为行，以样本为列的β值甲基化矩阵
2.以基因为为行，以样本为列并且标准化处理后的转录组矩阵（FPKM等）
3.甲基化样本名和转录组样本名对应的转化矩阵
4.样本本身的临床信息（metadata）
其中编号1和2的文件是必需的，后面两个文件我们可以不用。这里我们以TCGA项目中肺鳞癌（LUSC）项目为例，为减少运算负荷，我们只选择1号染色体的甲基化探针和基因进行分析。除此之外我们的探针也要提前做出筛选，我们需要获得远端探针（distal probe）的id。
distal.probes <- get.feature.probe(genome = "hg19",
met.platform = "450K",
rm.chr = paste0("chr",c(2:22,"X","Y")))
## Returning distal probes: 15561
#这里同样可以获取promoter区域内的探针
# promoter.probes <- get.feature.probe(genome = "hg19",
# met.platform = "450K",
# TSS.range = list(upstream = 2000, downstream = 2000),
# promoter = T,
# rm.chr = paste0("chr",c(2:22,"X","Y")))
#展示一下，有点长
head(distal.probes)
## GRanges object with 6 ranges and 52 metadata columns:
## seqnames ranges strand | address_A address_B channel
## <Rle> <IRanges> <Rle> | <integer> <integer> <character>
## cg00168193 chr1 790667-790668 - | 34713309 <NA> Both
## cg08258224 chr1 800083-800084 - | 23712478 <NA> Both
## cg13938959 chr1 834183-834184 + | 38601418 <NA> Both
## cg12445832 chr1 834295-834296 - | 70726486 <NA> Both
## cg23999112 chr1 834356-834357 - | 49747457 <NA> Both
## cg11527153 chr1 837536-837537 + | 47662324 <NA> Both
## designType nextBase nextBaseRef probeType orientation
## <character> <character> <character> <character> <character>
## cg00168193 II G/A C cg down
## cg08258224 II G/A C cg down
## cg13938959 II G/A C cg up
## cg12445832 II G/A C cg down
## cg23999112 II G/A C cg down
## cg11527153 II G/A C cg up
## probeCpGcnt context35 probeBeg probeEnd ProbeSeq_A
## <integer> <integer> <integer> <numeric> <character>
## cg00168193 1 2 790668 790717 TTCATATACACATATACTCR..
## cg08258224 0 1 800084 800133 CTAACATACAATCCAAAAAC..
## cg13938959 1 2 834134 834183 CATTAAAAAACAAAAATCAC..
## cg12445832 2 4 834296 834345 TTCCCTCRCAATTCCRAAAA..
## cg23999112 1 4 834357 834406 AAAAAAATACCCCCAACTCC..
## cg11527153 1 2 837487 837536 ACACCTCATCTTAAAACCAA..
## ProbeSeq_B gene gene_HGNC chrm_A beg_A
## <character> <character> <character> <character> <integer>
## cg00168193 LINC01128 LINC01128 chr1 790668
## cg08258224 <NA> <NA> chr1 800084
## cg13938959 <NA> <NA> chr1 834134
## cg12445832 <NA> <NA> chr1 834296
## cg23999112 <NA> <NA> chr1 834357
## cg11527153 <NA> <NA> chr1 837487
## flag_A mapQ_A cigar_A NM_A chrm_B beg_B
## <integer> <integer> <character> <integer> <character> <integer>
## cg00168193 16 53 50M 0 <NA> <NA>
## cg08258224 16 60 50M 0 <NA> <NA>
## cg13938959 0 60 50M 0 <NA> <NA>
## cg12445832 16 60 50M 0 <NA> <NA>
## cg23999112 16 60 50M 0 <NA> <NA>
## cg11527153 0 60 50M 0 <NA> <NA>
## flag_B mapQ_B cigar_B NM_B wDecoy_chrm_A
## <integer> <integer> <character> <integer> <character>
## cg00168193 <NA> <NA> <NA> <NA> chr1
## cg08258224 <NA> <NA> <NA> <NA> chr1
## cg13938959 <NA> <NA> <NA> <NA> chr1
## cg12445832 <NA> <NA> <NA> <NA> chr1
## cg23999112 <NA> <NA> <NA> <NA> chr1
## cg11527153 <NA> <NA> <NA> <NA> chr1
## wDecoy_beg_A wDecoy_flag_A wDecoy_mapQ_A wDecoy_cigar_A
## <integer> <integer> <integer> <character>
## cg00168193 790668 16 53 50M
## cg08258224 800084 16 60 50M
## cg13938959 834134 0 60 50M
## cg12445832 834296 16 60 50M
## cg23999112 834357 16 60 50M
## cg11527153 837487 0 60 50M
## wDecoy_NM_A wDecoy_chrm_B wDecoy_beg_B wDecoy_flag_B wDecoy_mapQ_B
## <integer> <character> <integer> <integer> <integer>
## cg00168193 0 <NA> <NA> <NA> <NA>
## cg08258224 0 <NA> <NA> <NA> <NA>
## cg13938959 0 <NA> <NA> <NA> <NA>
## cg12445832 0 <NA> <NA> <NA> <NA>
## cg23999112 0 <NA> <NA> <NA> <NA>
## cg11527153 0 <NA> <NA> <NA> <NA>
## wDecoy_cigar_B wDecoy_NM_B posMatch MASK_mapping
## <character> <integer> <logical> <logical>
## cg00168193 <NA> <NA> TRUE FALSE
## cg08258224 <NA> <NA> TRUE FALSE
## cg13938959 <NA> <NA> TRUE FALSE
## cg12445832 <NA> <NA> TRUE FALSE
## cg23999112 <NA> <NA> TRUE FALSE
## cg11527153 <NA> <NA> TRUE FALSE
## MASK_typeINextBaseSwitch MASK_rmsk15 MASK_sub40_copy
## <logical> <logical> <logical>
## cg00168193 FALSE TRUE FALSE
## cg08258224 FALSE FALSE FALSE
## cg13938959 FALSE FALSE FALSE
## cg12445832 FALSE FALSE FALSE
## cg23999112 FALSE FALSE FALSE
## cg11527153 FALSE FALSE FALSE
## MASK_sub35_copy MASK_sub30_copy MASK_sub25_copy MASK_snp5_common
## <logical> <logical> <logical> <logical>
## cg00168193 FALSE FALSE FALSE TRUE
## cg08258224 FALSE FALSE FALSE FALSE
## cg13938959 FALSE FALSE FALSE FALSE
## cg12445832 FALSE FALSE FALSE FALSE
## cg23999112 FALSE FALSE FALSE TRUE
## cg11527153 FALSE FALSE TRUE FALSE
## MASK_snp5_GMAF1p MASK_extBase MASK_general
## <logical> <logical> <logical>
## cg00168193 FALSE FALSE FALSE
## cg08258224 FALSE FALSE FALSE
## cg13938959 FALSE FALSE FALSE
## cg12445832 FALSE FALSE FALSE
## cg23999112 FALSE FALSE FALSE
## cg11527153 FALSE FALSE FALSE
## -------
## seqinfo: 26 sequences from an unspecified genome; no seqlengths
得到1号染色体的distal probe的ID后，我们就需要准备前面提到的甲基化和转录组数据啦。从而创建MAE对象，这里同样导入我们在ELMER.data中的数据
data(LUSC_RNA_refined,package = "ELMER.data")
data(LUSC_meth_refined,package = "ELMER.data")
#展示一下转录组数据
GeneExp[1:5,1:5]
## TCGA-22-5472-01A-01R-1635-07 TCGA-22-5489-01A-01R-1635-07
## ENSG00000188984 0.0000000 0.000000
## ENSG00000204518 0.4303923 0.000000
## ENSG00000108270 10.0817831 10.717673
## ENSG00000198691 6.4462711 6.386644
## ENSG00000135776 8.5929182 9.333097
## TCGA-22-5491-11A-01R-1858-07 TCGA-56-5898-01A-11R-1635-07
## ENSG00000188984 0.000000 0.5233612
## ENSG00000204518 0.000000 0.0000000
## ENSG00000108270 10.180863 10.1595162
## ENSG00000198691 5.755627 4.8354795
## ENSG00000135776 8.558358 8.7772810
## TCGA-90-6837-01A-11R-1949-07
## ENSG00000188984 0.000000
## ENSG00000204518 1.167294
## ENSG00000108270 9.975092
## ENSG00000198691 3.152329
## ENSG00000135776 8.604804
#展示一下甲基化数据
Meth[1:5,1:5]
## TCGA-43-3394-11A-01D-1551-05 TCGA-43-3920-11B-01D-1551-05
## cg00045114 0.8190894 0.8073763
## cg00050294 0.8423084 0.8241138
## cg00066722 0.9101127 0.9162212
## cg00093522 0.8751903 0.8864599
## cg00107046 0.3326016 0.3445508
## TCGA-56-8305-01A-11D-2294-05 TCGA-56-8307-01A-11D-2294-05
## cg00045114 0.8907009 0.8483227
## cg00050294 0.5597787 0.3488952
## cg00066722 0.7228874 0.6238963
## cg00093522 0.8050060 0.8194921
## cg00107046 0.4312738 0.4328108
## TCGA-56-8308-01A-11D-2294-05
## cg00045114 0.7612094
## cg00050294 0.3908054
## cg00066722 0.7727631
## cg00093522 0.7507631
## cg00107046 0.4260053
#开始创建MAE对象
mae <- createMAE(exp = GeneExp,
met = Meth,
save = TRUE,
linearize.exp = TRUE,
save.filename = "mae.rda",#同时保存数据为mae.rda
filter.probes = distal.probes,
met.platform = "450K",
genome = "hg19",
TCGA = TRUE)
## =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
## Creating a SummarizedExperiment from gene expression input
## =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
## Creating a SummarizedExperiment from DNA methylation input
## Checking if samples have both DNA methylation and Gene expression and if they are in the same order...
## Starting to add information to samples
## => Add clinical information to samples
## => Adding TCGA molecular information from marker papers
## => Information will have prefix 'paper_'
## lusc subtype information from:doi:10.1038/nature11404
## Creating MultiAssayExperiment
## MAE saved as: mae.rda
创建出MAE对象后我们可以对这个对象进行适当了解，看一下我们如何获取MEA对象中相应的数据。
#获取甲基化表达谱数据
getMet(mae) %>% assay() %>% as.data.frame() %>% .[1:5,1:5]
## TCGA-43-3394-11A-01D-1551-05 TCGA-56-8305-01A-11D-2294-05
## cg23625715 0.01539862 0.1097107
## cg26222311 0.02163782 0.1563133
## cg18432292 0.40672158 0.3114181
## cg24517686 0.65124883 0.4567180
## cg23926866 0.44623986 0.4382566
## TCGA-56-8307-01A-11D-2294-05 TCGA-56-8308-01A-11D-2294-05
## cg23625715 0.02174192 0.02259196
## cg26222311 0.03046606 0.02755373
## cg18432292 0.18142269 0.25889453
## cg24517686 0.41442389 0.47301394
## cg23926866 0.29348669 0.32135912
## TCGA-56-8309-01A-11D-2294-05
## cg23625715 0.01918069
## cg26222311 0.02575283
## cg18432292 0.37044511
## cg24517686 0.49564970
## cg23926866 0.46165625
#获取MEA对象转录组数据
getExp(mae) %>% assay() %>% as.data.frame() %>% .[1:5,1:5]
## TCGA-43-3394-11A-01R-1758-07 TCGA-56-8305-01A-11R-2296-07
## ENSG00000188984 0.0000000 0.000000
## ENSG00000204518 0.2796211 0.000000
## ENSG00000108270 3.4657077 3.526386
## ENSG00000198691 2.5867196 2.047156
## ENSG00000135776 3.2674645 3.320399
## TCGA-56-8307-01A-11R-2296-07 TCGA-56-8308-01A-11R-2296-07
## ENSG00000188984 0.000000 0.000000
## ENSG00000204518 0.000000 0.000000
## ENSG00000108270 3.548234 3.544454
## ENSG00000198691 2.116421 1.800524
## ENSG00000135776 3.400905 3.327408
## TCGA-56-8309-01A-11R-2296-07
## ENSG00000188984 0.000000
## ENSG00000204518 1.653326
## ENSG00000108270 3.494669
## ENSG00000198691 2.508659
## ENSG00000135776 3.248861
#获取MEA对象的metadata信息
colData(mae) %>% as.data.frame() %>% .[1:5,1:5]
## patient sample shortLetterCode
## TCGA-43-3394-11A TCGA-43-3394 TCGA-43-3394-11A NT
## TCGA-56-8305-01A TCGA-56-8305 TCGA-56-8305-01A TP
## TCGA-56-8307-01A TCGA-56-8307 TCGA-56-8307-01A TP
## TCGA-56-8308-01A TCGA-56-8308 TCGA-56-8308-01A TP
## TCGA-56-8309-01A TCGA-56-8309 TCGA-56-8309-01A TP
## definition sample_submitter_id
## TCGA-43-3394-11A Solid Tissue Normal TCGA-43-3394-11A
## TCGA-56-8305-01A Primary solid Tumor TCGA-56-8305-01A
## TCGA-56-8307-01A Primary solid Tumor TCGA-56-8307-01A
## TCGA-56-8308-01A Primary solid Tumor TCGA-56-8308-01A
## TCGA-56-8309-01A Primary solid Tumor TCGA-56-8309-01A
二、筛选出差异甲基化探针
我们通过比较病变组vs.正常组样本中所有远端探针（distal probe）的甲基化水平，并筛选出adj.P value<0.01，并且两组之间的Δβ＞0.3的甲基化位点，并且是在肿瘤组织低甲基化水平的位点。
除此之外，ELMER提供了两种方式进行筛选远端探针方式。
（1）supervised，即我们定义好的疾病组和正常组中所有样本进行分组，从而进行差异分析筛选出在疾病组低甲基化远端探针。
（2）unsupervised，即对甲基化探针在所有样本中甲基化水平进行排序，选取最高甲基化水平20%以及最低甲基化水平20%进行差异分析，从而得到在疾病组低甲基化远端探针。这一类做法的原因是：疾病（尤其是肿瘤）具有明确的异质性，通过提取20%的子集能够进一步找到疾病*特中**定的分子亚型。
我们这儿以unsupervised方式为例。
sig.diff <- get.diff.meth(data = mae,
group.col = "definition",
group1 = "Primary solid Tumor",
group2 = "Solid Tissue Normal",
minSubgroupFrac = 0.2,
sig.dif = 0.3,
diff.dir = "hypo", # 在group 1为低甲基化水平的差异探针
cores = 1,
dir.out ="result",
pvalue = 0.01)
## ELMER will search for probes hypomethylated in group Primary solid Tumor (n:226) compared to Solid Tissue Normal (n:8)
## Saving results
## Warning: The `path` argument of `write_csv()` is deprecated as of readr 1.4.0.
## Please use the `file` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Number of relevant probes found: 710
#展示一下结果
head(sig.diff)
## probe pvalue
## cg25340711 cg25340711 3.898829e-05
## cg24852892 cg24852892 5.418625e-24
## cg06358191 cg06358191 2.617048e-09
## cg06646682 cg06646682 8.147636e-18
## cg15844005 cg15844005 2.104992e-30
## cg06205922 cg06205922 2.302866e-22
## Primary.solid.Tumor_Minus_Solid.Tissue.Normal adjust.p
## cg25340711 -0.3193081 5.509361e-05
## cg24852892 -0.4664680 4.835534e-23
## cg06358191 -0.3832813 6.310122e-09
## cg06646682 -0.4730666 4.041818e-17
## cg15844005 -0.3818866 4.019066e-29
## cg06205922 -0.3646617 1.742537e-21
筛选出来了！！！！，我们可以看一下result中volcano plot。
效果不错，对吧？我们的文章又多了一张图了。
总结一下，本期现简单介绍了转录因子以及甲基化调控基因表达的潜在机制，并初步介绍ELMER probe处理流程以及怎么使用ELMER包进行甲基化位点的差异分析。下一次我们将揭开ELMER最激动人心的部分，靶基因的推断以及上游转录因子的预测。
好啦，我是风间琉璃，咱们下期再见~
—END—