前言
DNA 甲基化在許多細胞進程中扮演重要的角色炸枣，例如胚胎發(fā)育搭伤、基因印跡踊挠、X 染色體失活和維持染色體穩(wěn)定性乍桂。
在哺乳動物中，DNA 甲基化很少見效床，其產生位置分布在整個基因組中的確定的 CpG 序列中睹酌，但是卻很少在 CpG 島上發(fā)生甲基化。
CpG 島（CGI）是富含 GC 堿基的短間隔 DNA 序列剩檀。這些 CpG 島通常位于轉錄起始位置憋沿，它們的甲基化會導致基因沉默。
DNA 甲基化會抑制轉錄沪猴，因此辐啄，對 DNA 甲基化的研究對于理解癌癥中調控基因網絡至關重要。所以运嗜，差異甲基化區(qū)域（DMR）的檢測有助于我們研究調控基因網路壶辜。
差異甲基化分析
1. 樣本甲基化均值
我們首先對 DNA 甲基化數(shù)據進行預處理，450k 平臺的 DNA 甲基化數(shù)據有三種探針：
-
cg：CpG位點 -
ch：非CpG位點 -
rs：SNP芯片
最后一種探針可用于識別和跟蹤樣本洗出，應該在差異甲基化分析中被排除士复，所以要刪除 rs 探針图谷。同時為了消除性別的影響翩活，X阱洪、Y 染色體也應該排除在外。最后菠镇，去除包含 NA 值的探針冗荸。
在本示例中，我們分析的是結直腸癌的兩個亞型：
- 結腸癌：
COAD - 直腸癌：
READ
library(TCGAbiolinks)
library(SummarizedExperiment)
library(tidyverse)
#------------------------------------
# 獲取 DNA 同時檢測甲基化和表達的樣本
#------------------------------------
# 結腸和直腸數(shù)據
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 位點
獲取并處理完數(shù)據之后，我們要識別兩組之間的差異甲基化 CpG 位點隘梨。
我們使用的是 TCGAanalyze_DMR 函數(shù)程癌，該函數(shù)首先計算每個探針在分組之間的平均甲基化值（mean of the beta-value）差異，然后轴猎，使用 wilcoxon 檢驗兩組之間的表達差異嵌莉，并使用 BH 方法矯正 P 值。
默認的最小均值差為 0.15捻脖，FDR 值為 0.05
#------- 識別差異甲基化位點 ----------
res <- TCGAanalyze_DMC(
met,
# colData 函數(shù)獲取的矩陣中分組列名
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)
# 獲取臨床數(shù)據
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()
模體分析
在識別出了差異甲基化區(qū)域之后可婶，我們可能希望能從這些區(qū)域中找出某些獨特的特征沿癞。從這些受到累積或缺失 DNA 甲基化作用影響的位點中識別出候選的轉錄因子。
模體（motif）分析的目的是提取隱藏在大部分非功能性基因間序列中的小序列信號矛渴，這些小序列核苷酸信號（6-15bp）可能具有調控基因表達的生物學意義椎扬。這些序列被稱為調控模體。
我們對差異甲基化區(qū)域前后 100 個堿基的序列進行模體分析具温，然后使用 motifStack 包來生成多個具有不同相似性得分的模體
Bioconductor 中的 rGADEM 包為大規(guī)牡两ⅲ基因組序列的數(shù)據分析提供了一種高效的 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)
嗯...很不走運桂躏，運行上述代碼钻趋，并沒有找到模體。
# 獲取模體的數(shù)量
> nMotifs(gadem)
[1] 3
由于我們識別出的位點較少剂习，發(fā)現(xiàn)模體的可能性較低蛮位，為了讓我們能夠繼續(xù)下面的步驟，只能換兩種癌型再試試了
我們使用 GBM 和 LGG 兩種癌型各十個樣本鳞绕，運行上述代碼失仁，首先從 DNA 甲基化均值可以看出這兩種癌型的差異具有顯著性
而從下面的火山圖我們也可以看出差異位點明顯增加了很多
從這些差異甲基化位點來尋找模體，看來是很有希望的了
> nMotifs(gadem)
[1] 3
找到了三個们何，再看看這三個模體出現(xiàn)的次數(shù)
> 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 對照表
在 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"
