ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
122 stars 30 forks source link

转录组+甲基化来点新花样:ELMER包 #3052

Closed ixxmu closed 1 year ago

ixxmu commented 1 year ago

https://mp.weixin.qq.com/s/COMXQ1vlW_CjK5mCWSgJ3A

ixxmu commented 1 year ago

转录组+甲基化来点新花样:ELMER包 by YuYuFiSH

ELMER 多组学R包

如果你手上有甲基化和转录组两种数据,有太多分析的角度能辅助你完成一个漂亮的故事了

自我介绍

ELMER(Enhancer Linking by Methylation/Expression Relationships),增强子相关甲基化和表达相关网络

关键是以甲基化位点作为网络的中心枢纽,在下游筛选出差异甲基化位点对应的碱基序列与对应的靶基因,而上游找到富集的motif(转录因子结合基序即有规律有特殊功能的一串碱基序列),再反推出与这些motif相互作用的转录因子,从而构建出TF—甲基化位点——靶基因网络

主要流程

1.识别远端探针(distal probe):HM450k或者EPIC甲基化芯片

2.识别组间甲基化水平差异

3.识别差异甲基化探针的靶基因

4.确定差异甲基化差异探针以及靶基因相关探针富集的motif

5.确认与motif上甲基化相关的调控转录因子

R包第一步可以理解为识别关键位置的甲基化探针,这里指增强子(enhancer)或者启动子区域(promoter)的探针。普遍认为在顺式调控元件部位(包括启动子、增强子等部位)高度甲基化,将影响DNA的结构,从而阻遏该部位基因的转录

distal probe:距离转录起始位点位点上游大于2kb的甲基化探针,默认这一区域的甲基化探针区域为增强子区域。

最新2019年版本

新版本再次升级

🍥 帮助文档:https://academic.oup.com/bioinformatics/article/35/11/1974/5145133?login=false

http://bioconductor.org/packages/release/bioc/vignettes/ELMER/inst/doc/analysis_data_input.html

http://bioconductor.org/packages/release/bioc/manuals/ELMER/man/ELMER.pdf

🍥 工作流程workshop:

https://benbermanlab.com/assets/code/ELMER%20workshop.html

https://rpubs.com/tiagochst/elmer-data-workshop-2019

🍥 文章TCGA-BRCA完整报告:

https://media.githubusercontent.com/media/tiagochst/ELMER_supplemental/master/supplemental_files.zip

目录

  • 获得远端探针:get.feature.probe
  • 载入甲基化,转录组,临床数据
  • 创建MAE对象:createMAE
  • 差异分析:get.diff.meth
  • 寻找探针上下游基因:GetNearGenes
  • 识别probe-gene对:get.pair
  • 筛选motif:get.enriched.motif
  • 确定转录因子TF:get.TFs
  • 各种可视化:
    • 散点图:scatter.plot
    • 染色体示意图:schematic.plot
    • motif富集森林图:motif.enrichment.plot
    • 转录因子排序图:TF.rank.plot
    • 热图:heatmapPairs
  • 直接下载TCGA数据:TCGA.pipe

领域展开----

1.R包安装

两种方法:从github中安装 or Bioconductor安装

devtools::install_github(repo = "tiagochst/ELMER.data")
devtools::install_github(repo = "tiagochst/ELMER")

if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("ELMER")
BiocManager::install("sesameData")

library(ELMER)
library(sesameData)
library(DT)

2. 获得远端探针

这里以TCGA项目中肺鳞癌(LUSC)项目为例,只选择1号染色体的甲基化探针和基因作为例子进行分析。

首先,探针要提前做出筛选,怎么获得增强子和启动子位置上的呢?r包这里已经根据数据库GENCODE和UCSC分好类,不过如果需要启动子上的探针,需要用参数TSS.range提前定义好启动子区域,默认是TSS上2000bp内,不过也可以自定义,毕竟没有硬性标准

  • 获得远端探针:get.feature.probe
# Selection of probes within biofeatures
# get distal probes that are 2kb away from TSS on chromosome 1
distal.probes <- get.feature.probe(
  genome = "hg19"# hg38 (default) or hg19
  met.platform = "450K"#  EPIC or 450K (default)
  TSS.range = list(upstream = 2000, downstream = 2000), # define promoter regions定义启动子区域
  promoter = FALSE, # 是否选择启动子
  rm.chr = paste0("chr",c(2:22,"X","Y")) # 移除的染色体

length(distal.probes)# [1] 15561
head(distal.probes)

3. 输入数据(TCGA格式)

得到1号染色体的distal probe后,我们就需要

  • 载入甲基化,转录组数据
  • 创建MAE对象: createMAE
# Selection of probes within biofeatures
library(MultiAssayExperiment)
library(ELMER.data)
# 转录组数据
data(LUSC_RNA_refined,package = "ELMER.data")
# 甲基化数据
data(LUSC_meth_refined,package = "ELMER.data")

dim(GeneExp) #3842  234
dim(Meth) #1728  268

可以看到,两种类型的数据样本数量不完全一致

> GeneExp[1:3,1:3]
                TCGA-22-5472-01A-01R-1635-07 TCGA-22-5489-01A-01R-1635-07
ENSG00000188984                    0.0000000                      0.00000
ENSG00000204518                    0.4303923                      0.00000
ENSG00000108270                   10.0817831                     10.71767
                TCGA-22-5491-11A-01R-1858-07
ENSG00000188984                      0.00000
ENSG00000204518                      0.00000
ENSG00000108270                     10.18086
> Meth[1:3,1:3]
           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
           TCGA-56-8305-01A-11D-2294-05
cg00045114                    0.8907009
cg00050294                    0.5597787
cg00066722                    0.7228874


  • 创建MAE对象:createMAE
mae <- createMAE(
  exp = GeneExp, 
  met = Meth,
  save = TRUE,
  linearize.exp = TRUE, #  log2(exp + 1) 
  save.filename = "mae.rda"# 保存的文件
  filter.probes = distal.probes, # promoter regions or distal feature regions
  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对象将数据进行整合,提取同时拥有两种数据的234个样本。可以根据下面提示的函数探索下MAE对象

> mae
A MultiAssayExperiment object of 2 listed
 experiments with user-defined names and respective classes.
 Containing an ExperimentList class object of length 2:
 [1] DNA methylation: RangedSummarizedExperiment with 1642 rows and 234 columns
 [2] Gene expression: RangedSummarizedExperiment with 3808 rows and 234 columns
Functionality:
 experiments() - obtain the ExperimentList instance
 colData() - the primary/phenotype DataFrame
 sampleMap() - the sample coordination DataFrame
 `$`, `[`, `[[` - extract colData columns, subset, or experiment
 *Format() - convert into a long or wide DataFrame
 assays() - convert ExperimentList to a SimpleList of matrices
 exportClass() - save data to flat files
# 解析mae对象
#获取甲基化表达谱数据
getMet(mae) %>% assay()
#获取MEA对象转录组数据
getExp(mae) %>% assay()
#获取MEA对象的metadata信息
colData(mae) %>% as.data.frame()

3.2 输入数据(非TCGA格式)

当然,如果你有自己的测序数据,也可以构建属于自己的MAE对象咯

需要准备:

1.以探针为行,以样本为列的β值甲基化矩阵

2.以基因为行,以样本为列并且标准化处理后的转录组矩阵(FPKM等)

3.甲基化样本名和转录组样本名对应的转化矩阵

4.样本本身的临床信息(metadata)

> head(exp)
                Sample1 Sample2 Sample3 Sample4 Sample5
ENSG00000073282       0       0       0       0       0
ENSG00000078900       0       0       0       0       0
ENSG00000141510       0       0       0       0       0
> head(met)
           Sample1 Sample2 Sample3 Sample4 Sample5
cg26928153       0       0       0       0       0
cg16269199       0       0       0       0       0
cg13869341       0       0       0       0       0
> head(colData) # 后续可以增加多几列临床信息
         sample
Sample1 Sample1
Sample2 Sample2
Sample3 Sample3
Sample4 Sample4
Sample5 Sample5
> sampleMap # 指定样本类型,防止有些样本名不一致,此数据集三个列名固定不变
             assay primary colname
1  DNA methylation Sample1 Sample1
2  DNA methylation Sample2 Sample2
3  DNA methylation Sample3 Sample3
4  DNA methylation Sample4 Sample4
5  DNA methylation Sample5 Sample5
6  Gene expression Sample1 Sample1
7  Gene expression Sample2 Sample2
8  Gene expression Sample3 Sample3
9  Gene expression Sample4 Sample4
10 Gene expression Sample5 Sample5

这种情况下,你就照着葫芦画瓢,模仿文档给出的数据格式构造即可

distal.probes <- get.feature.probe(
  genome = "hg19"
  met.platform = "EPIC"
)

mae <- createMAE(
  exp = exp, 
  met = met,
  save = TRUE,
  filter.probes = distal.probes,
  colData = colData,
  sampleMap = sampleMap,
  linearize.exp = TRUE,
  save.filename = "mae.rda",
  met.platform = "EPIC",
  genome = "hg19",
  TCGA = FALSE
)

4. 筛选差异甲基化探针

  • 差异分析:get.diff.meth

ELMER提供了两种方式进行筛选差异探针

(1)supervised,即使用未配对的单侧t检验,将组1中所有样本的的DNA甲基化水平与组2中的所有样本进行比较

(2)unsupervised,即根据每个探针的甲基化水平分别对两个组的样本先进行排序,每组中选择最低20%甲基化水平的样本用于识别此探针是否低甲基化,高甲基化探针的鉴定正好相反。(原因是:疾病(尤其是肿瘤)具有明确的异质性,通过提取20%的子集能够进一步找到疾病中特定的分子亚型)

因此,在无监督模式下,分析每个探针的样本可能基于总样本的不同子集,可以表示来自多个分子亚型的探针集的结果。而在监督模式下,所有测试都基于所有样本。

这里以unsupervised方式为例,通过比较不同组间样本中所有远端探针(distal probe)的甲基化水平,筛选出adj.P value<0.01,Δβ>0.3的甲基化位点,并且是在肿瘤组织低甲基化水平的位点

### 差异分析----
table(colData(mae)[,"definition"])
# Primary solid Tumor Solid Tissue Normal 
# 226                   8 
sig.diff <- get.diff.meth(data = mae, 
                          group.col = "definition"# 组别类型对应的列名
                          group1 =  "Primary solid Tumor"#  if direction is hyper, get probes hypermethylated in group 1 compared to group 2
                          group2 = "Solid Tissue Normal",
                          minSubgroupFrac = 0.2,  # 默认指定排序后的20%甲基化水平的样本进行比较
                          sig.dif = 0.3, # 组间最小beta值差异
                          diff.dir = "hypo"# "hypo", "hyper" or "both"
                          cores = 1, # 并行核数
                          dir.out ="result",
                          mode = "unsupervised",
                          pvalue = 0.01)

head(sig.diff)
dir(path = "result", pattern = "getMethdiff"

这里要注意参数diff.dir = "hypo",意味着在group1为低甲基化水平的探针

"hypo" which is only selecting hypomethylated probes (one tailed test); "hyper" which is only selecting hypermethylated probes (one tailed test); or "both" which are probes differenly methylated (two tailed test).

输出:所有和显著差异位点的csv结果,以及一张火山图

> dim(sig.diff) #710   4
> head(sig.diff)
       probe       pvalue Primary.solid.Tumor_Minus_Solid.Tissue.Normal
1 cg25340711 3.898829e-05                                    -0.3193081
2 cg24852892 5.418625e-24                                    -0.4664680
3 cg06358191 2.617048e-09                                    -0.3832813
4 cg06646682 8.147636e-18                                    -0.4730666
5 cg15844005 2.104992e-30                                    -0.3818866
6 cg06205922 2.302866e-22                                    -0.3646617
      adjust.p
1 5.509361e-05
2 4.835534e-23
3 6.310122e-09
4 4.041818e-17
5 4.019066e-29
6 1.742537e-21

> dir(path = "result", pattern = "getMethdiff")  
[1] "getMethdiff.hypo.probes.csv"            
[2] "getMethdiff.hypo.probes.significant.csv"

5. 识别甲基化-基因对

这一步,我们将把远端探针的甲基化水平和靶基因的表达水平进行相关性分析,从而构建出甲基化-基因对(probe-gene)。这个包具体的做法是分别找到差异甲基化远端探针的上游和下游最近的10个基因,分析其余对应探针甲基化水平是否存在负相关,从而筛选出来的潜在probe-gene对。

同样提供了两类方式:

(1)supervised,即通过上游分析疾病组vs.正常组差异远端探针,在所有样本中进行非参数检,找出同时在疾病组高表达的基因

(2)unsupervised,根据基因的distal probe甲基化水平进行排序,提取出甲基化水平分别在前20%和后20%的样本,分别为M和U组,通过U检验方式比较两组间基因表达水平,如果在M组的表达水平显著低于U组,并通过迭代的方式计算出所有远端探针和靶基因的相关性P value,从而筛选出具有意义的关系对。


  • 寻找探针上下游基因:GetNearGenes
### 识别probe-gene对-----
# 加载我们之前的数据
mae <- get(load("mae.rda"))
#加载我们筛选得到的差异甲基化远端探针(distal probe)
sig.diff <- read.csv("result/getMethdiff.hypo.probes.significant.csv")
#找到差异甲基化远端探针(distal probe)上下游10个基因
nearGenes <- GetNearGenes(data = mae, 
                          probes = sig.diff$probe
                          numFlankingGenes = 20) # 分别是上游10个和下游10基因
## Searching for the 20 near genes
## Identifying gene position for each probe
head(nearGenes)


  • 识别probe-gene对:get.pair
#因为前面使用unsupervised,我们这而也用unsupervised
hypo.pair <- get.pair(data = mae,
                      group.col = "definition",
                      group1 =  "Primary solid Tumor",
                      group2 = "Solid Tissue Normal",
                      nearGenes = nearGenes,
                      minSubgroupFrac = 0.4,# 默认创建M和U组样本的比例
                      mode = "unsupervised",
                      permu.dir = "result/permu",
                      permu.size = 100, # 一般要设置迭代次数为10000,处于演示目的减少次数
                      raw.pvalue = 0.05,   
                      Pe = 0.01, # P值,默认0.001
                      filter.probes = TRUE, 
                      filter.percentage = 0.05, # 最小考虑的样本比例
                      filter.portion = 0.3, # 甲基化水平截断点
                      dir.out = "result",#结果同样输出到result目录
                      cores = 1,
                      label = "hypo")
# Selecting U (unmethylated) and M (methylated) groups. Each groups has 20% of samples
# -------------------
#   * Filtering probes
# -------------------
#   For more information see function preAssociationProbeFiltering
# Making sure we have at least 5% of beta values lesser than 0.3 and 5% of beta values greater 0.3.
# Removing 592 probes out of 1642
# Calculating Pp (probe - gene) for all nearby genes
# |============================================|100%                      Completed after 5 s 
# File created: result/getPair.hypo.all.pairs.statistic.csv             474.73MB/s
# Calculating Pr (random probe - gene). Permutating  100 probes for 747 genes
# Calculate empirical P value.              
head(hypo.pair)
dim(hypo.pair) # [1] 167  7
> head(nearGenes)
# A tibble: 6 × 5
# Groups:   ID [1]
  ID         GeneID          Symbol  Distance Side 
  <chr>      <chr>           <chr>      <dbl> <chr>
1 cg00045114 ENSG00000076258 FMO4    -1362935 L10  
2 cg00045114 ENSG00000117523 PRRC2C  -1111508 L9   
3 cg00045114 ENSG00000034971 MYOC    -1052335 L8   
4 cg00045114 ENSG00000117533 VAMP4    -962771 L7   
5 cg00045114 ENSG00000010165 METTL13  -890995 L6   
6 cg00045114 ENSG00000197959 DNM3     -286552 L5 

> head(hypo.pair)
                                Probe          GeneID Symbol Distance
cg20701183.ENSG00000143013 cg20701183 ENSG00000143013   LMO4        0
cg19403323.ENSG00000143469 cg19403323 ENSG00000143469  SYT14    87476
cg12213388.ENSG00000143674 cg12213388 ENSG00000143674   MLK4  -984182
cg26607897.ENSG00000143199 cg26607897 ENSG00000143199 ADCY10   267784
cg10574861.ENSG00000143013 cg10574861 ENSG00000143013   LMO4        0
cg26607897.ENSG00000143147 cg26607897 ENSG00000143147 GPR161   543156
                           Sides        Raw.p         Pe
cg20701183.ENSG00000143013    L1 7.453984e-14 0.00990099
cg19403323.ENSG00000143469    R1 1.671937e-12 0.00990099
cg12213388.ENSG00000143674    L4 2.527644e-12 0.00990099
cg26607897.ENSG00000143199    R3 4.593610e-12 0.00990099
cg10574861.ENSG00000143013    L1 4.770162e-12 0.00990099
cg26607897.ENSG00000143147    R6 8.048248e-12 0.00990099

> dir(path = "result", pattern = "getPair"
[1] "getPair.hypo.all.pairs.statistic.csv"                  
[2] "getPair.hypo.pairs.significant.csv"                    
[3] "getPair.hypo.pairs.statistic.with.empirical.pvalue.csv"

6. 筛选motif

提取探针上下游250bp区域的碱基序列,从而找到富集的motif。接下来再通过转录因子(TF)结合motif数据库(HOMER数据库),从而预测出motif结合的转录因子

  • 筛选motif:get.enriched.motif
### 识别motif-----

enriched.motif <- get.enriched.motif(data = mae,
                                     probes = hypo.pair$Probe
                                     dir.out = "result"
                                     label = "hypo",
                                     min.incidence = 1, #  the minimum incidence of motif
                                     lower.OR = 1.1) # motif最低95%置信区间OR阈值
## Loading object: Probes.motif.hg19.450K

names(enriched.motif)
enriched.motif[[1]]
# 识别出六种motif
> names(enriched.motif)
[1] "HXA9_HUMAN.H11MO.0.B"  "FOXH1_HUMAN.H11MO.0.A" "PO5F1_HUMAN.H11MO.1.A"
[4] "HXB13_HUMAN.H11MO.0.A" "IRF9_HUMAN.H11MO.0.C"  "SOX10_HUMAN.H11MO.1.A"

# 展示第一个motif包含的甲基化探针
> enriched.motif[[1]]
 [1] "cg13067635" "cg24873093" "cg25210796" "cg18345456" "cg14094320"
 [6] "cg17901924" "cg22736624" "cg12213388" "cg13305336" "cg26607897"
[11] "cg11718886" "cg24593832" "cg25729466" "cg24446417"

# 输出的结果文件
> dir(path = "result", pattern = "getMotif"
[1] "getMotif.hypo.enriched.motifs.rda"  "getMotif.hypo.motif.enrichment.csv"

# 输出的图片
> dir(path = "result", pattern = "motif.enrichment.pdf"
[1] "hypo.quality.A-DS_with_summary.motif.enrichment.pdf"
[2] "hypo.quality.A-DS.motif.enrichment.pdf"

7. 确定转录因子TF

这一步,ELMER包将上一步我们通过富集分析找到富集的motif以及对应motif的转录因子进行分析,筛选出对应转录水平发生改变的转录因子。

同样地,提供了两种方式:

(1)unsupervised,即将存在相同motif的远端探针(distal probe)根据甲基化水平分为高甲基化组(一般取前20%)和低甲基化组(后20%),比较其两组间对应TF的表达值是否存在差异。

(2)Supervised,即对疾病组vs.对照组相同motif对应TF表达值进行分析。进而筛选出前5% P value最小的TFs,并视为潜在的上游调控转录因子。

  • 确定转录因子TF:get.TFs
## 确定转录因子TF-----
TF <- get.TFs(data = mae, 
              group.col = "definition",
              group1 =  "Primary solid Tumor",
              group2 = "Solid Tissue Normal",
              mode = "unsupervised",
              enriched.motif = enriched.motif,
              dir.out = "result"
              cores = 1, 
              label = "hypo")
head(TF)
dir(path = "result", pattern = "getTF"
> dir(path = "result", pattern = "getTF"
[1] "getTF.hypo.significant.TFs.with.motif.summary.csv"
[2] "getTF.hypo.TFs.with.motif.pvalue.rda" 


总结一下,我们通过同一批样本的甲基化数据和转录组数据,首先找到具有差异的远端探针(distal probe),并根据远端探针和对应基因的表达水平构建出probe-gene对。并将probe-gene对中的probe进一步进行motif富集分析,找出潜在富集motif,最后根据富集的motif推测出上游的调控TFs

8. 各种可视化

  • 散点图:scatter.plot
  • 染色体示意图:schematic.plot
  • motif富集森林图:motif.enrichment.plot
  • 转录因子排序图:TF.rank.plot
  • 热图:heatmapPairs
### 各种可视化----
##### 1.散点图----
#加载数据
mae <- get(load("mae.rda"))
load("result/getMotif.hypo.enriched.motifs.rda")
#探针和上下游分别10个基因的散点图
scatter.plot(data = mae,
             byProbe = list(probe = c("cg19403323"), 
                            numFlankingGenes = 20), 
             category = "definition"
             lm = TRUE, # 画出线条
             save = FALSE)
## Searching for the 20 near genes
## Identifying gene position for each probe

#只展现一对probe-gene散点图
scatter.plot(data = mae,
             byPair = list(probe = c("cg19403323"), gene = c("ENSG00000143469")), 
             category = "definition", save = TRUE, lm_line = TRUE)


#特定转录因子表达水平以及特定motif对应探针甲基化水平散点图
scatter.plot(data = mae,
             byTF = list(TF = c("SOX2"),
                         probe = enriched.motif[[names(enriched.motif)[1]]]), 
             category = "definition",
             save = F, #保存为PDF版本
             lm_line = TRUE)
##### 2.染色体示意图----
# 使用genome browser的方式展示远端探针(distal probe)甲基化水平和附近20个基因在染色体上的位置。
library(GenomicInteractions)
#加载probe-gene对
pair <- read.csv("result/getPair.hypo.pairs.significant.csv")
#展现探针附近的gene
schematic.plot(pair = pair, 
               data = mae,
               group.col = "definition",
               byProbe = pair$Probe[1],
               save = FALSE)
## Searching for the 20 near genes
## Identifying gene position for each probe

#展示基因附近的probe
schematic.plot(pair = pair, 
               data = mae,   
               group.col = "definition"
               byGene = pair$GeneID[1],
               save = FALSE)

##### 3.motif富集森林图-----
#初步展示
motif.enrichment.plot(motif.enrichment = "result/getMotif.hypo.motif.enrichment.csv"
                      significant = list(OR = 1.5,lowerOR = 1.3), 
                      label = "hypo"
                      save = FALSE)

#展现OR值和探针数
motif.enrichment.plot(motif.enrichment = "result/getMotif.hypo.motif.enrichment.csv"
                      significant = list(OR = 1.5,lowerOR = 1.3), 
                      label = "hypo"
                      summary = TRUE,
                      save = FALSE)
##### 4.转录因子排序图-----
load("result/getTF.hypo.TFs.with.motif.pvalue.rda")
motif <- colnames(TF.meth.cor)[1]
TF.rank.plot(motif.pvalue = TF.meth.cor, 
             motif = motif,
             save = FALSE)
## Retrieving TFClass family classification from ELMER.data.
## Retrieving TFClass subfamily classification from ELMER.data.
## |                                                    |  0%                      |====================================================|100% ~0 s remaining       |====================================================|100%                      Completed after 0 s
## $HXA9_HUMAN.H11MO.0.B
## Warning: ggrepel: 69 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
##### 5. 热图-----
heatmapPairs(data = mae, 
             group.col = "definition",
             group1 = "Primary solid Tumor"
             annotation.col = c("years_smoked","gender"),
             group2 = "Solid Tissue Normal",
             pairs = pair,
             filename =  NULL)

9. 直接下载TCGA数据

当然,可以直接使用TCGA.pipe,联合r包TCGAbiolinks直接下载TCGA的数据分析

The following command will do distal DNA methylation analysis and predict putative target genes, motif analysis and identify regulatory transcription factors.

TCGA.pipe("LUSC",
          wd = "./ELMER.example",
          cores = parallel::detectCores()/2,
          mode = "unsupervised",
          permu.size = 300,
          Pe = 0.01,
          analysis = c("distal.probes","diffMeth","pair","motif","TF.search"),
          diff.dir = "hypo",
          rm.chr = paste0("chr",c("X","Y")))
  • 还可以直接在线分析
### 在线版本----
BiocManager::install("TCGAbiolinksGUI")
library(TCGAbiolinksGUI)
TCGAbiolinksGUI()

10. N个好奇心

10.1 Supervised versus Unsupervised mode?

可以看到前面很多功能,作者都提供两种模式,监督 or 非监督模式,那在实际工作中采用哪个比较多呢?

在第一版ELMER v. 1时,作者推荐当两组样本均为同质性(即处理过的vs未处理的,A分子亚型vs B分子亚型)时,应采用监督模式;当至少有一组样本为异质性(即肿瘤样本)时,应采用无监督模式

而最新版,作者使用ELMER v. 2.4.3分析TCGA-BRCA数据,探索增强子-基因对(enhancer-gene)时,发现监督模式有更强的统计效能,并且能识别更多的增强子-基因对

This comparison suggests that while Unsupervised mode can serve as a useful exploratory tool when sample subtype is unknown a priori, Supervised mode offers greater statistical power when sample subtype is pre-defined.

而对于产生的假阳性率情况,作者建议两种方法都可以尝试下

For heterogeneous patient samples composed of multiple subtypes, it thus appears that Unsupervised and Supervised analyses can offer complementary merits, with Unsupervised analysis displaying a higher false negative rate, but a lower false positive rate. It is recommended to run both Supervised and Unsupervised analyses, as we demonstrated here, to gain maximum insight.

10.2 能分析非肿瘤数据吗?

虽然文章中测评用的都是肿瘤数据,并且能搭配TCGABiolinks一起使用。我觉得还是可以用的,很多数据挖掘类型R包都可以通用,可能会有局限,可以搜索下有没有类似文章

we previously developed an R/Bioconductor package ELMER (Enhancer Linking by Methylation/Expression Relationships) to infer regulatory element landscapes and GRNs from cancer methylomes