TCGA 数据分析实战 —— 差异甲基化区域模体分析
前言
DNA 甲基化在许多细胞进程中扮演重要的角色，例如胚胎发育、基因印迹、X 染色体失活和维持染色体稳定性。
在哺乳动物中，DNA 甲基化很少见，其产生位置分布在整个基因组中的确定的 CpG 序列中，但是却很少在 CpG 岛上发生甲基化。
CpG 岛（CGI）是富含 GC 碱基的短间隔 DNA 序列。这些 CpG 岛通常位于转录起始位置，它们的甲基化会导致基因沉默。
DNA 甲基化会抑制转录，因此，对 DNA 甲基化的研究对于理解癌症中调控基因网络至关重要。所以，差异甲基化区域（DMR）的检测有助于我们研究调控基因网路。
差异甲基化分析
1. 样本甲基化均值
我们首先对 DNA 甲基化数据进行预处理，450k 平台的 DNA 甲基化数据有三种探针：
-
cg：CpG位点 -
ch：非CpG位点 -
rs：SNP芯片
最后一种探针可用于识别和跟踪样本，应该在差异甲基化分析中被排除，所以要删除 rs 探针。同时为了消除性别的影响，X、Y 染色体也应该排除在外。最后，去除包含 NA 值的探针。
在本示例中，我们分析的是结直肠癌的两个亚型：
- 结肠癌：
COAD - 直肠癌：
READ
library(TCGAbiolinks)
library(SummarizedExperiment)
library(tidyverse)
#------------------------------------
# 获取 DNA 同时检测甲基化和表达的样本
#------------------------------------
# 结肠和直肠数据
coad.samples <- matchedMetExp("TCGA-COAD", n = 10)
read.samples <- matchedMetExp("TCGA-READ", n = 10)
samples <- c(coad.samples, read.samples)
#------------------------------------
# 1 - Methylation
# -----------------------------------
query <- GDCquery(
project = c("TCGA-COAD", "TCGA-READ"),
data.category = "DNA methylation",
platform = "Illumina Human Methylation 450",
legacy = TRUE,
barcode = samples
)
GDCdownload(query)
met <- GDCprepare(query, save = FALSE)
# 我们以 chr9 为例
met <- subset(met,subset = as.character(seqnames(met)) %in% c("chr9"))
# 删除包含 NA 值的探针
met <- subset(met,subset = (rowSums(is.na(assay(met))) == 0))
# 去除重复样本
met <- met[, substr(colnames(met), 14, 16) != "01B"]
我们根据 project_id 对样本进行分组，并使用 T 检验计算样本甲基化值之间的差异
TCGAvisualize_meanMethylation(
met,
groupCol = "project_id",
group.legend = "Groups",
filename = NULL,
print.pvalue = TRUE
)
==================== DATA Summary ====================
groups Mean Median Max Min
1 TCGA-COAD 0.4670816 0.4682460 0.5076270 0.4222167
2 TCGA-READ 0.4730700 0.4758822 0.4904898 0.4498025
整体看，两组样本之间的 P 值并不显著
2. 识别差异甲基化 CpG 位点
获取并处理完数据之后，我们要识别两组之间的差异甲基化 CpG 位点。
我们使用的是 TCGAanalyze_DMR 函数，该函数首先计算每个探针在分组之间的平均甲基化值（mean of the beta-value）差异，然后，使用 wilcoxon 检验两组之间的表达差异，并使用 BH 方法矫正 P 值。
默认的最小均值差为 0.15，FDR 值为 0.05
#------- 识别差异甲基化位点 ----------
res <- TCGAanalyze_DMC(
met,
# colData 函数获取的矩阵中分组列名
groupCol = "project_id",
group1 = "TCGA-COAD",
group2 = "TCGA-READ",
p.cut = 0.05,
diffmean.cut = 0.15,
save = FALSE,
legend = "State",
plot.filename = "~/Downloads/COAD_READ_metvolcano.png",
cores = 1
)
从这个火山图的结果中，我们可以看到有 6 个低甲基化位点
最后，我们使用热图来显示这些探针在所有样本中的 DNA 甲基化水平
#--------------------------
# DNA Methylation heatmap
#--------------------------
library(ComplexHeatmap)
# 获取临床数据
coad_clin <- GDCquery_clinic(project = "TCGA-COAD", type = "Clinical")
read_clin <- GDCquery_clinic(project = "TCGA-READ", type = "Clinical")
use_cols <- c("bcr_patient_barcode", "disease","gender","vital_status","race")
clinical <- coad_clin %>%
dplyr::select(use_cols) %>%
add_row(dplyr::select(read_clin, use_cols)) %>%
subset(bcr_patient_barcode %in% substr(samples, 1, 12))
# 获取 Hypermethylated 和 Hypomethylated 的探针
sig_met <- filter(res, status != "Not Significant")
# 获取探针的表达值
res_data <- subset(met,subset = (rownames(met) %in% rownames(sig_met)))
# 添加注释
ta <- HeatmapAnnotation(
df = clinical[, c("disease", "gender", "vital_status", "race")],
col = list(
disease = c("COAD" = "grey", "READ" = "black"),
gender = c("male" = "blue", "female" = "pink")
))
ra <- rowAnnotation(
df = sig_met$status,
col = list(
"status" =
c("Hypomethylated" = "orange",
"Hypermethylated" = "darkgreen")
),
width = unit(1, "cm")
)
heatmap <- Heatmap(
assay(res_data),
name = "DNA methylation",
col = matlab::jet.colors(200),
show_row_names = FALSE,
cluster_rows = TRUE,
cluster_columns = FALSE,
show_column_names = FALSE,
bottom_annotation = ta,
column_title = "DNA Methylation"
)
# 保存结果
png("~/Downloads/heatmap.png",width = 600, height = 400)
draw(heatmap, annotation_legend_side = "bottom")
dev.off()
模体分析
在识别出了差异甲基化区域之后，我们可能希望能从这些区域中找出某些独特的特征。从这些受到累积或缺失 DNA 甲基化作用影响的位点中识别出候选的转录因子。
模体（motif）分析的目的是提取隐藏在大部分非功能性基因间序列中的小序列信号，这些小序列核苷酸信号（6-15bp）可能具有调控基因表达的生物学意义。这些序列被称为调控模体。
我们对差异甲基化区域前后 100 个碱基的序列进行模体分析，然后使用 motifStack 包来生成多个具有不同相似性得分的模体
Bioconductor 中的 rGADEM 包为大规模基因组序列的数据分析提供了一种高效的 de novo 模体识别算法。
#---------------------------
# motif 分析
#---------------------------
library(rGADEM)
library(GenomicRanges)
library(BSgenome.Hsapiens.UCSC.hg19)
library(motifStack)
probes <- rowRanges(res_data)
# 获取差异的探针并设置 200bp 大小的窗口
sequence <- GRanges(
seqnames = as.character(seqnames(probes)),
ranges = IRanges(start = start(ranges(probes)) - 100,
end = start(ranges(probes)) + 100),
strand = "*"
)
# 识别模体
gadem <- GADEM(sequence, verbose = FALSE, genome = Hsapiens)
嗯...很不走运，运行上述代码，并没有找到模体。
# 获取模体的数量
> nMotifs(gadem)
[1] 3
由于我们识别出的位点较少，发现模体的可能性较低，为了让我们能够继续下面的步骤，只能换两种癌型再试试了
我们使用 GBM 和 LGG 两种癌型各十个样本，运行上述代码，首先从 DNA 甲基化均值可以看出这两种癌型的差异具有显著性
而从下面的火山图我们也可以看出差异位点明显增加了很多
从这些差异甲基化位点来寻找模体，看来是很有希望的了
> nMotifs(gadem)
[1] 3
找到了三个，再看看这三个模体出现的次数
> nOccurrences(gadem)
[1] 187 110 143
查看模体的一致性序列
> consensus(gadem)
[1] "nrrAGrrArArrrrrnrArr" "nCTTCCTsn" "sCCCnsmnCCCn"
展示模体 logo 图
pwm <- getPWM(gadem)
pfm <- new("pfm",mat=pwm[[1]],name="Novel Site 1")
plotMotifLogo(pfm)
一个序列的 logo 图展示的是多序列比对之后可能的碱基堆积标记，logo 的高度与序列的保守程度有关，碱基或氨基酸的保守程度越高，字母越高大，每个位置的碱基根据频率从高到低进行堆叠显示，可以从序列的顶端读取一致性序列
序列的 IUPAC codes 对照表
DNA
Protein
在 rGADEM 获取到结果之后，还可以使用 MotIV 进行模体配对分析
library(MotIV)
# motif 配对分析
> analysis.jaspar <- motifMatch(pwm)
Ungapped Alignment
Scores read
Database read
Motif matches : 5
> summary(analysis.jaspar)
Number of input motifs : 3
Input motifs names : TGGGGTGAGGGGAAGGCGGA TCTTCCTGG CCCCGCCTCCCC
Number of matches per motif: 5
Matches repartition :
SPI1 SPIB ELF5 EWSR1-FLI1 FEV IRF1 IRF2 Klf4 Pax4
2 2 1 1 1 1 1 1 1
PLAG1 SP1 Stat3 TFAP2A
1 1 1 1
Arguments used :
-metric name : PCC
-alignment : SWU
绘制配对分析的结果
plot(analysis.jaspar, ncol=2, top=5, rev=FALSE, main="", bysim=TRUE, cex=0.4)
查看比对结果的质量
> alignment <- viewAlignments(analysis.jaspar)
> print(alignment[[1]])
EWSR1-FLI1 IRF1 SPIB SPI1
seq "NRGAGRNANARRRRNNNARN-" "NRGAGRNANARRRRNNNARN" "NRGAGRNANARRRRNNNARN" "NRGAGRNANARRRRNNNARN"
match "---GGAAGGAAGGAAGGAAGG" "----NAAANYGAAACC----" "-------ASMGGAA------" "---AGGAAGT----------"
evalue "3.0277e-06" "1.4371e-03" "6.4823e-03" "1.5616e-02"
IRF2
seq "NRGAGRNANARRRRNNNARN-"
match "---SGAAAGYGAAASYNWWNN"
evalue "1.9415e-02"
