TCGA 数据分析实战 —— 差异甲基化区域模体分析

前言

DNA 甲基化在许多细胞进程中扮演重要的角色，例如胚胎发育、基因印迹、X 染色体失活和维持染色体稳定性。

在哺乳动物中，DNA 甲基化很少见，其产生位置分布在整个基因组中的确定的 CpG 序列中，但是却很少在 CpG 岛上发生甲基化。

CpG 岛（CGI）是富含 GC 碱基的短间隔 DNA 序列。这些 CpG 岛通常位于转录起始位置，它们的甲基化会导致基因沉默。

DNA 甲基化会抑制转录，因此，对 DNA 甲基化的研究对于理解癌症中调控基因网络至关重要。所以，差异甲基化区域（DMR）的检测有助于我们研究调控基因网路。

差异甲基化分析

样本甲基化均值

我们首先对 DNA 甲基化数据进行预处理，450k 平台的 DNA 甲基化数据有三种探针：

• cg：CpG 位点
• ch：非 CpG 位点
• rs：SNP 芯片

最后一种探针可用于识别和跟踪样本，应该在差异甲基化分析中被排除，所以要删除 rs 探针。同时为了消除性别的影响，X、Y 染色体也应该排除在外。最后，去除包含 NA 值的探针。

在本示例中，我们分析的是非小细胞肺癌的两个亚型：

• 肺腺癌：LUAD
• 肺鳞癌：LUSC

library(TCGAbiolinks)
library(SummarizedExperiment)
library(tidyverse)

#------------------------------------
# 获取 DNA 同时检测甲基化和表达的样本
#------------------------------------
# 肺腺癌和肺鳞癌
luad.samples <- matchedMetExp("TCGA-LUAD", n = 10)
lusc.samples <- matchedMetExp("TCGA-LUSC", n = 10)
samples <- c(luad.samples, lusc.samples)

query <- GDCquery(
    project = c("TCGA-LUAD", "TCGA-LUSC"),
    data.category = "DNA Methylation", 
    data.type = "Methylation Beta Value",
    platform = "Illumina Human Methylation 450",
    barcode = samples
)
GDCdownload(query)
met <- GDCprepare(query, save = FALSE)

# 删除包含 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
)

整体看，两组样本之间的 P 值并不显著

识别差异甲基化 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-LUAD",
    group2 = "TCGA-LUSC",
    p.cut = 0.05,
    diffmean.cut = 0.15,
    save = FALSE,
    legend = "State",
    plot.filename = "~/Downloads/metvolcano.png",
    cores =  1
)

从这个火山图的结果中，可以看出有挺多显著差异的甲基化位点。

最后，我们使用热图来显示这些探针在所有样本中的 DNA 甲基化水平

#--------------------------
# DNA Methylation heatmap
#--------------------------
library(ComplexHeatmap)

# 获取临床数据
luad_clin <- GDCquery_clinic(project = "TCGA-LUAD", type = "Clinical")
lusc_clin <- GDCquery_clinic(project = "TCGA-LUSC", type = "Clinical")

use_cols <- c("bcr_patient_barcode", "tissue_or_organ_of_origin","gender","vital_status","race")
clinical <- luad_clin %>% 
    dplyr::select(use_cols) %>%
    add_row(dplyr::select(lusc_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("LUAD" = "grey", "LUSC" = "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)[,-7],  # 一个样本多次取样，先随机删除一个
    name = "DNA methylation",
    col = matlab::jet.colors(21),
    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 包来生成多个具有不同相似性得分的模体

首先获取 JASPAR2018 中已知的 motifs

#---------------------------
#          motif 分析
#---------------------------
library(JASPAR2018)
library(TFBSTools)

# 获取 JASPAR2018 中的 motif
opts <- list()
opts["species"] <- 9606
opts["collection"] <- "CORE"
out  <- getMatrixSet(JASPAR2018, opts)


if (!isTRUE(all.equal(TFBSTools::name(out), names(out)))) 
    names(out) <- paste(names(out), TFBSTools::name(out), sep = "_")
motifs <- out

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 = "*"
)

我们使用 Bioconductor 中的 motifmatchr 包来识别 motif。

library(motifmatchr)

# 获取 score
motif_scores <- matchMotifs(motifs, sequence, genome = "hg38", 
                         out = "scores")
# 获取富集到的 motif 位置
motif_pos <- matchMotifs(motifs, sequence, genome = "hg38",
                         out = "positions")
motif_pos
# GRangesList object of length 452:
# $MA0025.1_NFIL3
# GRanges object with 22 ranges and 1 metadata column:
#        seqnames              ranges strand |     score
#           <Rle>           <IRanges>  <Rle> | <numeric>
#    [1]     chr1   42462465-42462475      - |   9.03835
#    [2]     chr1   47533338-47533348      - |   9.54086
#    [3]     chr1   47533338-47533348      - |   9.54086
#    [4]     chr1   94898657-94898667      - |   8.73465
#    [5]     chr2   42814041-42814051      + |  11.43626

其中 MA0025.1_NFIL3 为 motif 的名字，GRanges 类型的属性为富集到该 motif 的位置

绘制该 motif

library(motifStack)

pfm  <- new("pfm", mat = pcm2pfm(motifs$MA0025.1_NFIL3@profileMatrix), 
            name = "MA0025.1_NFIL3")
plotMotifLogo(pfm)

一个序列的 logo 图展示的是多序列比对之后可能的碱基堆积标记，logo 的高度与序列的保守程度有关，碱基或氨基酸的保守程度越高，字母越高大，每个位置的碱基根据频率从高到低进行堆叠显示，可以从序列的顶端读取一致性序列