更多示例

參考鏈接

1.為基因表達(dá)矩陣添加更多信息

熱圖非常流行用于可視化基因表達(dá)矩陣。矩陣中的行對(duì)應(yīng)于基因，在表達(dá)熱圖之后可以添加關(guān)于這些基因的更多信息。

在下面的例子中，熱圖可視化了基因的相對(duì)表達(dá)（每個(gè)基因的表達(dá)被縮放）。在右側(cè)，我們將基因的絕對(duì)表達(dá)水平作為單列熱圖?；蜷L(zhǎng)度和基因類(lèi)型（即蛋白質(zhì)編碼或lincRNA）也作為熱圖注釋或熱圖。

在熱圖的最左側(cè)，繪制了彩色矩形 anno_block()以識(shí)別 k 均值聚類(lèi)中的五個(gè)聚類(lèi)。在“基本均值”和“基因類(lèi)型”熱圖之上，還有匯總圖（條形圖和箱形圖）顯示五個(gè)集群中數(shù)據(jù)點(diǎn)的統(tǒng)計(jì)數(shù)據(jù)或分布。

library(ComplexHeatmap)
library(circlize)

expr = readRDS(system.file(package = "ComplexHeatmap", "extdata", "gene_expression.rds"))
mat = as.matrix(expr[, grep("cell", colnames(expr))])
base_mean = rowMeans(mat)
mat_scaled = t(apply(mat, 1, scale))

type = gsub("s\\d+_", "", colnames(mat))
ha = HeatmapAnnotation(type = type, annotation_name_side = "left")

ht_list = Heatmap(mat_scaled, name = "expression", row_km = 5, 
    col = colorRamp2(c(-2, 0, 2), c("green", "white", "red")),
    top_annotation = ha, 
    show_column_names = FALSE, row_title = NULL, show_row_dend = FALSE) +
Heatmap(base_mean, name = "base mean", 
    top_annotation = HeatmapAnnotation(summary = anno_summary(gp = gpar(fill = 2:6), 
        height = unit(2, "cm"))),
    width = unit(15, "mm")) +
rowAnnotation(length = anno_points(expr$length, pch = 16, size = unit(1, "mm"), 
    axis_param = list(at = c(0, 2e5, 4e5, 6e5), 
        labels = c("0kb", "200kb", "400kb", "600kb")),
    width = unit(2, "cm"))) +
Heatmap(expr$type, name = "gene type", 
    top_annotation = HeatmapAnnotation(summary = anno_summary(height = unit(2, "cm"))),
    width = unit(15, "mm"))

ht_list = rowAnnotation(block = anno_block(gp = gpar(fill = 2:6, col = NA)), 
    width = unit(2, "mm")) + ht_list

draw(ht_list)

image

2.麻疹疫苗熱圖

以下代碼重現(xiàn)了此處和此處介紹的熱圖。

mat = readRDS(system.file("extdata", "measles.rds", package = "ComplexHeatmap"))
ha1 = HeatmapAnnotation(
    dist1 = anno_barplot(
        colSums(mat), 
        bar_width = 1, 
        gp = gpar(col = "white", fill = "#FFE200"), 
        border = FALSE,
        axis_param = list(at = c(0, 2e5, 4e5, 6e5, 8e5),
            labels = c("0", "200k", "400k", "600k", "800k")),
        height = unit(2, "cm")
    ), show_annotation_name = FALSE)
ha2 = rowAnnotation(
    dist2 = anno_barplot(
        rowSums(mat), 
        bar_width = 1, 
        gp = gpar(col = "white", fill = "#FFE200"), 
        border = FALSE,
        axis_param = list(at = c(0, 5e5, 1e6, 1.5e6),
            labels = c("0", "500k", "1m", "1.5m")),
        width = unit(2, "cm")
    ), show_annotation_name = FALSE)
year_text = as.numeric(colnames(mat))
year_text[year_text %% 10 != 0] = ""
ha_column = HeatmapAnnotation(
    year = anno_text(year_text, rot = 0, location = unit(1, "npc"), just = "top")
)
col_fun = colorRamp2(c(0, 800, 1000, 127000), c("white", "cornflowerblue", "yellow", "red"))
ht_list = Heatmap(mat, name = "cases", col = col_fun,
    cluster_columns = FALSE, show_row_dend = FALSE, rect_gp = gpar(col= "white"), 
    show_column_names = FALSE,
    row_names_side = "left", row_names_gp = gpar(fontsize = 8),
    column_title = 'Measles cases in US states 1930-2001\nVaccine introduced 1961',
    top_annotation = ha1, bottom_annotation = ha_column,
    heatmap_legend_param = list(at = c(0, 5e4, 1e5, 1.5e5), 
        labels = c("0", "50k", "100k", "150k"))) + ha2
draw(ht_list, ht_gap = unit(3, "mm"))
decorate_heatmap_body("cases", {
    i = which(colnames(mat) == "1961")
    x = i/ncol(mat)
    grid.lines(c(x, x), c(0, 1), gp = gpar(lwd = 2, lty = 2))
    grid.text("Vaccine introduced", x, unit(1, "npc") + unit(5, "mm"))
})

image

3.從單細(xì)胞 RNASeq 中可視化細(xì)胞異質(zhì)性

在這個(gè)例子中，小鼠 T 細(xì)胞的單細(xì)胞 RNA-Seq 數(shù)據(jù)被可視化以顯示細(xì)胞的異質(zhì)性。 ( mouse_scRNAseq_corrected.txt) 數(shù)據(jù)來(lái)自Buettner 等人，2015 年，補(bǔ)充數(shù)據(jù) 1mouse_scRNAseq_corrected.txt，表“細(xì)胞周期校正基因表達(dá)”。你可以到 這里。

在以下代碼中，刪除了重復(fù)的基因。

expr = read.table("data/mouse_scRNAseq_corrected.txt", sep = "\t", header = TRUE)
expr = expr[!duplicated(expr[[1]]), ]
rownames(expr) = expr[[1]]
expr = expr[-1]
expr = as.matrix(expr)

過(guò)濾掉一半以上細(xì)胞中未表達(dá)的基因。

expr = expr[apply(expr, 1, function(x) sum(x > 0)/length(x) > 0.5), , drop = FALSE]

該get_correlated_variable_rows()函數(shù)在這里被定義。它提取在細(xì)胞之間可變表達(dá)與其他基因相關(guān)的特征基因。

get_correlated_variable_genes = function(mat, n = nrow(mat), cor_cutoff = 0, n_cutoff = 0) {
    ind = order(apply(mat, 1, function(x) {
            q = quantile(x, c(0.1, 0.9))
            x = x[x < q[1] & x > q[2]]
            var(x)/mean(x)
        }), decreasing = TRUE)[1:n]
    mat2 = mat[ind, , drop = FALSE]
    dt = cor(t(mat2), method = "spearman")
    diag(dt) = 0
    dt[abs(dt) < cor_cutoff] = 0
    dt[dt < 0] = -1
    dt[dt > 0] = 1

    i = colSums(abs(dt)) > n_cutoff

    mat3 = mat2[i, ,drop = FALSE]
    return(mat3)
}

特征基因被定義為一個(gè)基因列表，其中每個(gè)基因的絕對(duì)相關(guān)性超過(guò) 0.5 并且有 20 多個(gè)相關(guān)基因。

mat2包含按基因縮放的表達(dá)值，這意味著它包含每個(gè)基因跨細(xì)胞的相對(duì)表達(dá)。由于單細(xì)胞 RNASeq 數(shù)據(jù)高度可變且異常值頻繁出現(xiàn)，基因表達(dá)僅在第10和第90分位數(shù)內(nèi)縮放。

mat = get_correlated_variable_genes(expr, cor_cutoff = 0.5, n_cutoff = 20)
mat2 = t(apply(mat, 1, function(x) {
    q10 = quantile(x, 0.1)
    q90 = quantile(x, 0.9)
    x[x < q10] = q10
    x[x > q90] = q90
    scale(x)
}))
colnames(mat2) = colnames(mat)

加載細(xì)胞周期基因和核糖核蛋白基因。細(xì)胞周期基因列表來(lái)自Buettner et al., 2015 , 補(bǔ)充表 1, “ Union of Cyclebase and GO gene”。核糖核蛋白基因來(lái)自 GO:0030529。基因列表存儲(chǔ)在 mouse_cell_cycle_gene.rds和mouse_ribonucleoprotein.rds中。這兩個(gè)文件可以在這里和 這里找到 。

cc = readRDS("data/mouse_cell_cycle_gene.rds")
ccl = rownames(mat) %in% cc
cc_gene = rownames(mat)[ccl]

rp = readRDS("data/mouse_ribonucleoprotein.rds")
rpl = rownames(mat) %in% rp

由于縮放每個(gè)基因的表達(dá)值，一個(gè)基因相對(duì)于其他基因的表達(dá)水平已經(jīng)丟失，我們將基本平均值計(jì)算為所有樣本中基因的平均表達(dá)?；酒骄悼捎糜诒容^基因之間的表達(dá)水平。

base_mean = rowMeans(mat)

現(xiàn)在可以使用以下信息：

1. 縮放表達(dá)式, mat2,
2. 基本均值, base_mean,
3. 基因是否為核糖核蛋白基因rpl，
4. 基因是否是細(xì)胞周期基因ccl，
5. 細(xì)胞周期基因的符號(hào), cc_gene,

在下一步中，我們可以將信息放在一起并將其可視化為熱圖列表。最后添加了一個(gè)基因-基因相關(guān)熱圖，并定義為 main_heatmap，這意味著所有熱圖/行注釋的行順序都基于此相關(guān)矩陣的聚類(lèi)。

對(duì)于表達(dá)水平相對(duì)較高的細(xì)胞周期基因（大于所有基因的 25% 分位數(shù)），基因名稱(chēng)以文本標(biāo)簽表示。在第一個(gè)熱圖中，列樹(shù)狀圖下面有兩種不同的顏色，這些顏色基于通過(guò)層次聚類(lèi)得出的兩個(gè)主要組，以突出顯示兩個(gè)亞群。

library(GetoptLong)
ht_list = Heatmap(mat2, col = colorRamp2(c(-1.5, 0, 1.5), c("blue", "white", "red")), 
    name = "scaled_expr", column_title = qq("relative expression for @{nrow(mat)} genes"),
    show_column_names = FALSE, width = unit(8, "cm"),
    heatmap_legend_param = list(title = "Scaled expr")) +
    Heatmap(base_mean, name = "base_expr", width = unit(5, "mm"),
        heatmap_legend_param = list(title = "Base expr")) +
    Heatmap(rpl + 0, name = "ribonucleoprotein", col = c("0" = "white", "1" = "purple"), 
        show_heatmap_legend = FALSE, width = unit(5, "mm")) +
    Heatmap(ccl + 0, name = "cell_cycle", col = c("0" = "white", "1" = "red"), 
        show_heatmap_legend = FALSE, width = unit(5, "mm")) +
    rowAnnotation(link = anno_mark(at = which(ccl & base_mean > quantile(base_mean, 0.25)), 
        labels = rownames(mat)[ccl & base_mean > quantile(base_mean, 0.25)], 
        labels_gp = gpar(fontsize = 10), padding = unit(1, "mm"))) +
    Heatmap(cor(t(mat2)), name = "cor", 
        col = colorRamp2(c(-1, 0, 1), c("green", "white", "red")), 
        show_row_names = FALSE, show_column_names = FALSE, row_dend_side = "right", 
        show_column_dend = FALSE, column_title = "pairwise correlation between genes",
        heatmap_legend_param = list(title = "Correlation"))
ht_list = draw(ht_list, main_heatmap = "cor")
decorate_column_dend("scaled_expr", {
    tree = column_dend(ht_list)$scaled_expr
    ind = cutree(as.hclust(tree), k = 2)[order.dendrogram(tree)]

    first_index = function(l) which(l)[1]
    last_index = function(l) { x = which(l); x[length(x)] }
    x1 = c(first_index(ind == 1), first_index(ind == 2)) - 1
    x2 = c(last_index(ind == 1), last_index(ind == 2))
    grid.rect(x = x1/length(ind), width = (x2 - x1)/length(ind), just = "left",
        default.units = "npc", gp = gpar(fill = c("#FF000040", "#00FF0040"), col = NA))
})

image

熱圖清楚地表明細(xì)胞被分成兩個(gè)亞群。第一個(gè)熱圖中左側(cè)的群體表現(xiàn)出細(xì)胞周期基因子集的高表達(dá)（細(xì)胞周期基因在“ cell_cycle ”熱圖中表示）。然而，這些基因的整體表達(dá)水平相對(duì)較低（參見(jiàn)“ base_expr ”熱圖）。右側(cè)的群體在其他特征基因中具有更高的表達(dá)。有趣的是，在該亞群中高表達(dá)的特征基因富含編碼核糖核蛋白的基因（參見(jiàn)“核糖核蛋白”熱圖）。核糖核蛋白基因的一個(gè)子集顯示出強(qiáng)烈的共表達(dá)（參見(jiàn)相關(guān)熱圖）和整體高表達(dá)水平（“ base_expr”熱圖）。

4.甲基化、表達(dá)值和其他基因組特征之間的相關(guān)性

在以下示例中，數(shù)據(jù)是根據(jù)未發(fā)表分析中發(fā)現(xiàn)的模式隨機(jī)生成的。首先我們加載數(shù)據(jù)。meth.rds可以在這里找到。

res_list = readRDS("data/meth.rds")
type = res_list$type
mat_meth = res_list$mat_meth
mat_expr = res_list$mat_expr
direction = res_list$direction
cor_pvalue = res_list$cor_pvalue
gene_type = res_list$gene_type
anno_gene = res_list$anno_gene
dist = res_list$dist
anno_enhancer = res_list$anno_enhancer

不同的信息來(lái)源和相應(yīng)的變量是：

1. type：顯示樣本是腫瘤還是正常的標(biāo)簽。
2. mat_meth：一個(gè)矩陣，其中的行對(duì)應(yīng)于差異甲基化區(qū)域 (DMR)。矩陣中的值是每個(gè)樣品中 DMR 中的平均甲基化水平。
3. mat_expr：一個(gè)矩陣，其中的行對(duì)應(yīng)于與 DMR 相關(guān)的基因（即最接近 DMR 的基因）。矩陣中的值是每個(gè)樣本中每個(gè)基因的表達(dá)水平。表達(dá)針對(duì)樣本中的每個(gè)基因進(jìn)行縮放。
4. direction：甲基化變化的方向（hyper 表示腫瘤樣本中較高的甲基化，hypo 表示腫瘤樣本中較低的甲基化）。
5. cor_pvalue：甲基化和相關(guān)基因表達(dá)之間相關(guān)性檢驗(yàn)的 p 值。
6. gene_type：基因的類(lèi)型（例如蛋白質(zhì)編碼基因或lincRNA）。
7. anno_gene: 基因模型的注釋?zhuān)ɑ蜷g、基因內(nèi)或 TSS）。
8. dist：相關(guān)基因的 DMR 到 TSS 的距離。
9. anno_enhancer：與增強(qiáng)子重疊的 DMR 的分?jǐn)?shù)。

數(shù)據(jù)僅包括相關(guān)基因的甲基化和表達(dá)呈負(fù)相關(guān)的 DMR。

首先計(jì)算甲基化矩陣的列聚類(lèi)，以便可以將表達(dá)矩陣中的列調(diào)整為與甲基化矩陣中的列順序相同。

column_tree = hclust(dist(t(mat_meth)))
column_order = column_tree$order

library(RColorBrewer)
meth_col_fun = colorRamp2(c(0, 0.5, 1), c("blue", "white", "red"))
direction_col = c("hyper" = "red", "hypo" = "blue")
expr_col_fun = colorRamp2(c(-2, 0, 2), c("green", "white", "red"))
pvalue_col_fun = colorRamp2(c(0, 2, 4), c("white", "white", "red"))
gene_type_col = structure(brewer.pal(length(unique(gene_type)), "Set3"), 
    names = unique(gene_type))
anno_gene_col = structure(brewer.pal(length(unique(anno_gene)), "Set1"), 
    names = unique(anno_gene))
dist_col_fun = colorRamp2(c(0, 10000), c("black", "white"))
enhancer_col_fun = colorRamp2(c(0, 1), c("white", "orange"))

我們首先定義兩列注釋?zhuān)缓笾谱鲝?fù)雜的熱圖。

ht_opt(
    legend_title_gp = gpar(fontsize = 8, fontface = "bold"), 
    legend_labels_gp = gpar(fontsize = 8), 
    heatmap_column_names_gp = gpar(fontsize = 8),
    heatmap_column_title_gp = gpar(fontsize = 10),
    heatmap_row_title_gp = gpar(fontsize = 8)
)

ha = HeatmapAnnotation(type = type, 
    col = list(type = c("Tumor" = "pink", "Control" = "royalblue")),
    annotation_name_side = "left")
ha2 = HeatmapAnnotation(type = type, 
    col = list(type = c("Tumor" = "pink", "Control" = "royalblue")), 
    show_legend = FALSE)

ht_list = Heatmap(mat_meth, name = "methylation", col = meth_col_fun,
    column_order= column_order,
    top_annotation = ha, column_title = "Methylation") +
    Heatmap(direction, name = "direction", col = direction_col) +
    Heatmap(mat_expr[, column_tree$order], name = "expression", 
        col = expr_col_fun, 
        column_order = column_order, 
        top_annotation = ha2, column_title = "Expression") +
    Heatmap(cor_pvalue, name = "-log10(cor_p)", col = pvalue_col_fun) +
    Heatmap(gene_type, name = "gene type", col = gene_type_col) +
    Heatmap(anno_gene, name = "anno_gene", col = anno_gene_col) +
    Heatmap(dist, name = "dist_tss", col = dist_col_fun) +
    Heatmap(anno_enhancer, name = "anno_enhancer", col = enhancer_col_fun, 
        cluster_columns = FALSE, column_title = "Enhancer")

draw(ht_list, row_km = 2, row_split = direction,
    column_title = "Comprehensive correspondence between methylation, expression and other genomic features", 
    column_title_gp = gpar(fontsize = 12, fontface = "bold"), 
    merge_legends = TRUE, heatmap_legend_side = "bottom")

image

ht_opt(RESET = TRUE)

復(fù)雜的熱圖顯示高度甲基化的 DMR 富含基因間和基因內(nèi)區(qū)域，很少與增強(qiáng)子重疊。相比之下，低甲基化的 DMR 富含轉(zhuǎn)錄起始位點(diǎn) (TSS) 和增強(qiáng)子。

5. 使用復(fù)雜注釋可視化甲基化概況

在這個(gè)例子中，Strum et al., 2012 中的圖 1 重新實(shí)現(xiàn)了一些調(diào)整。

需要先加載包。

library(matrixStats)
library(GenomicRanges)

甲基化圖譜可以從GEO 數(shù)據(jù)庫(kù)下載。該GEOquery 包用于從GEO檢索數(shù)據(jù)。

library(GEOquery)
gset = getGEO("GSE36278")

已通過(guò) Illumina HumanMethylation450 BeadChip 陣列測(cè)量甲基化圖譜。我們通過(guò)IlluminaHumanMethylation450kanno.ilmn12.hg19包加載探針數(shù)據(jù)。

調(diào)整矩陣中的行名稱(chēng)與探針相同。

library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
data(Locations)

mat = exprs(gset[[1]])
colnames(mat) = phenoData(gset[[1]])@data$title
mat = mat[rownames(Locations), ] 

probe包含探針的位置以及 CpG 位點(diǎn)是否與 SNP 重疊的信息。在這里，我們?nèi)コ诵匀旧w上的探針和SNP 重疊的探針。

data(SNPs.137CommonSingle)
data(Islands.UCSC)
l = Locations$chr %in% paste0("chr", 1:22) & is.na(SNPs.137CommonSingle$Probe_rs)
mat = mat[l, ]

相應(yīng)地獲取探針位置的子集和對(duì) CpG 島的注釋。

cgi = Islands.UCSC$Relation_to_Island[l]
loc = Locations[l, ]

將基質(zhì)分為腫瘤樣品基質(zhì)和正常樣品基質(zhì)。還要修改腫瘤樣本的列名以和表型數(shù)據(jù)一致。

mat1 = as.matrix(mat[, grep("GBM", colnames(mat))])   # tumor samples
mat2 = as.matrix(mat[, grep("CTRL", colnames(mat))])  # normal samples
colnames(mat1) = gsub("GBM", "dkfz", colnames(mat1))

表型數(shù)據(jù)來(lái)自Sturm 等人，2012 年，補(bǔ)充表 S1，可在此處找到。

表型數(shù)據(jù)的行被調(diào)整為與甲基化矩陣的列相同。

phenotype = read.table("data/450K_annotation.txt", header = TRUE, sep = "\t", 
    row.names = 1, check.names = FALSE, comment.char = "", stringsAsFactors = FALSE)
phenotype = phenotype[colnames(mat1), ]

請(qǐng)注意，我們僅使用來(lái)自 DKFZ 的 136 個(gè)樣本，而在Sturm et al., 2012 中，使用了額外的 74 個(gè) TCGA 樣本。

提取腫瘤樣本中甲基化變化最大的前 8000 個(gè)探針，并相應(yīng)地對(duì)其他信息進(jìn)行子集化。

ind = order(rowVars(mat1, na.rm = TRUE), decreasing = TRUE)[1:8000]
m1 = mat1[ind, ]
m2 = mat2[ind, ]
cgi2 = cgi[ind]
cgi2 = ifelse(grepl("Shore", cgi2), "Shore", cgi2)
cgi2 = ifelse(grepl("Shelf", cgi2), "Shelf", cgi2)
loc = loc[ind, ]

對(duì)于每個(gè)探頭，找到到最近的 TSS 的距離。pc_tx_tss.bed包含來(lái)自蛋白質(zhì)編碼基因的 TSS 位置。

gr = GRanges(loc[, 1], ranges = IRanges(loc[, 2], loc[, 2]+1))
tss = read.table("data/pc_tx_tss.bed", stringsAsFactors = FALSE)
tss = GRanges(tss[[1]], ranges = IRanges(tss[, 2], tss[, 3]))

tss_dist = distanceToNearest(gr, tss)
tss_dist = tss_dist@elementMetadata$distance

因?yàn)?code>NA矩陣中有一些( sum(is.na(m1))/length(m1)= 0.0011967) 會(huì)破壞cor()函數(shù)，所以我們替換NA為中間甲基化 (0.5)。請(qǐng)注意，盡管ComplexHeatmap允許NA在矩陣中，但刪除NA會(huì)加速聚類(lèi)。

m1[is.na(m1)] = 0.5
m2[is.na(m2)] = 0.5

以下注釋將添加到甲基化矩陣的列中：

1. 年齡
2. DKFZ 亞型分類(lèi)
3. TCGA 亞型分類(lèi)
4. 基于表達(dá)譜的 TCGA 亞型分類(lèi)
5. IDH1突變
6. H3F3A突變
7. TP53突變
8. chr7增益
9. chr10 損失
10. CDKN2A 缺失
11. EGFR擴(kuò)增
12. PDGFRA 擴(kuò)增

在下面的代碼中，我們?cè)?code>ha變量中定義了列注釋。我們還自定義注釋的顏色、圖例和高度。

mutation_col = structure(names = c("MUT", "WT", "G34R", "G34V", "K27M"), 
    c("black", "white", "#4DAF4A", "#4DAF4A", "#377EB8"))
cnv_col = c("gain" = "#E41A1C", "loss" = "#377EB8", "amp" = "#E41A1C", 
    "del" = "#377EB8", "normal" = "white")
ha = HeatmapAnnotation(
    age = anno_points(phenotype[[13]], 
        gp = gpar(col = ifelse(phenotype[[13]] > 20, "black", "red")), 
        height = unit(3, "cm")),
    dkfz_cluster = phenotype[[1]],
    tcga_cluster = phenotype[[2]],
    tcga_expr = phenotype[[3]],
    IDH1 = phenotype[[5]],
    H3F3A = phenotype[[4]],
    TP53 = phenotype[[6]],
    chr7_gain = ifelse(phenotype[[7]] == 1, "gain", "normal"),
    chr10_loss = ifelse(phenotype[[8]] == 1, "loss", "normal"),
    CDKN2A_del = ifelse(phenotype[[9]] == 1, "del", "normal"),
    EGFR_amp = ifelse(phenotype[[10]] == 1, "amp", "normal"),
    PDGFRA_amp = ifelse(phenotype[[11]] == 1, "amp", "normal"),
    col = list(dkfz_cluster = structure(names = c("IDH", "K27", "G34", "RTK I PDGFRA", 
            "Mesenchymal", "RTK II Classic"), brewer.pal(6, "Set1")),
        tcga_cluster = structure(names = c("G-CIMP+", "Cluster #2", "Cluster #3"), 
            brewer.pal(3, "Set1")),
        tcga_expr = structure(names = c("Proneural", "Classical", "Mesenchymal"), 
            c("#377EB8", "#FFFF33", "#FF7F00")),
        IDH1 = mutation_col,
        H3F3A = mutation_col,
        TP53 = mutation_col,
        chr7_gain = cnv_col,
        chr10_loss = cnv_col,
        CDKN2A_del = cnv_col,
        EGFR_amp = cnv_col,
        PDGFRA_amp = cnv_col),
    na_col = "grey", border = TRUE,
    show_legend = c(TRUE, TRUE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE),
    show_annotation_name = FALSE,
    annotation_legend_param = list(
        dkfz_cluster = list(title = "DKFZ Methylation"),
        tcga_cluster = list(title = "TCGA Methylation"),
        tcga_expr = list(title = "TCGA Expression"),
        H3F3A = list(title = "Mutations"),
        chr7_gain = list(title = "CNV"))
)

在最后的圖中，添加了四個(gè)熱圖。從左到右，有

1. 腫瘤樣本中甲基化的熱圖
2. 正常樣本中的甲基化
3. 到最近的 TSS 的距離
4. CpG 島 (CGI) 注釋。

熱圖根據(jù) CGI 注釋按行拆分。

繪制熱圖后，decorate_*()函數(shù)會(huì)添加附加圖形，例如注釋標(biāo)簽 。

col_fun = colorRamp2(c(0, 0.5, 1), c("#377EB8", "white", "#E41A1C"))
ht_list = Heatmap(m1, col = col_fun, name = "Methylation",
    clustering_distance_columns = "spearman",
    show_row_dend = FALSE, show_column_dend = FALSE,
    show_column_names = FALSE,
    bottom_annotation = ha, column_title = qq("GBM samples (n = @{ncol(m1)})"),
    row_split = factor(cgi2, levels = c("Island", "Shore", "Shelf", "OpenSea")), 
    row_title_gp = gpar(col = "#FFFFFF00")) + 
Heatmap(m2, col = col_fun, show_column_names = FALSE, 
    show_column_dend = FALSE, column_title = "Controls",
    show_heatmap_legend = FALSE, width = unit(1, "cm")) +
Heatmap(tss_dist, name = "tss_dist", col = colorRamp2(c(0, 2e5), c("white", "black")), 
    width = unit(5, "mm"),
    heatmap_legend_param = list(at = c(0, 1e5, 2e5), labels = c("0kb", "100kb", "200kb"))) + 
Heatmap(cgi2, name = "CGI", show_row_names = FALSE, width = unit(5, "mm"),
    col = structure(names = c("Island", "Shore", "Shelf", "OpenSea"), c("red", "blue", "green", "#CCCCCC")))
draw(ht_list, row_title = paste0("DNA methylation probes (n = ", nrow(m1), ")"),
    annotation_legend_side = "left", heatmap_legend_side = "left")

annotation_titles = c(dkfz_cluster = "DKFZ Methylation",
    tcga_cluster = "TCGA Methylation",
    tcga_expr = "TCGA Expression",
    IDH1 = "IDH1",
    H3F3A = "H3F3A",
    TP53 = "TP53",
    chr7_gain = "Chr7 gain",
    chr10_loss = "Chr10 loss",
    CDKN2A_del = "Chr10 loss",
    EGFR_amp = "EGFR amp",
    PDGFRA_amp = "PDGFRA amp")
for(an in names(annotation_titles)) {
    decorate_annotation(an, {
        grid.text(annotation_titles[an], unit(-2, "mm"), just = "right")
        grid.rect(gp = gpar(fill = NA, col = "black"))
    })
}
decorate_annotation("age", {
    grid.text("Age", unit(8, "mm"), just = "right")
    grid.rect(gp = gpar(fill = NA, col = "black"))
    grid.lines(unit(c(0, 1), "npc"), unit(c(20, 20), "native"), gp = gpar(lty = 2))
})
decorate_annotation("IDH1", {
    grid.lines(unit(c(-40, 0), "mm"), unit(c(1, 1), "npc"))
})
decorate_annotation("chr7_gain", {
    grid.lines(unit(c(-40, 0), "mm"), unit(c(1, 1), "npc"))
})

image

6. 添加多個(gè)箱線圖到單行

注釋函數(shù)anno_boxplot()只為一行繪制一個(gè)箱線圖。當(dāng)多個(gè)熱圖連接在一起，或者已經(jīng)定義了列的組時(shí)，對(duì)于每一行，我們要在熱圖或列組之間進(jìn)行比較，因此，需要為每一行繪制多個(gè)箱線圖。

在以下示例中，我們演示了如何實(shí)現(xiàn)為單行繪制多個(gè)箱線圖的注釋函數(shù)。該grid.boxplot()函數(shù)來(lái)自ComplexHeatmap包，它可以很容易地在網(wǎng)格系統(tǒng)下繪制箱線圖。

m1 = matrix(sort(rnorm(100)), 10, byrow = TRUE)
m2 = matrix(sort(rnorm(100), decreasing = TRUE), 10, byrow = TRUE)

ht_list = Heatmap(m1, name = "m1") + Heatmap(m2, name = "m2")

rg = range(c(m1, m2))
rg[1] = rg[1] - (rg[2] - rg[1])* 0.02
rg[2] = rg[2] + (rg[2] - rg[1])* 0.02
anno_multiple_boxplot = function(index) {
    nr = length(index)
    pushViewport(viewport(xscale = rg, yscale = c(0.5, nr + 0.5)))
    for(i in seq_along(index)) {
        grid.rect(y = nr-i+1, height = 1, default.units = "native")
        grid.boxplot(m1[ index[i], ], pos = nr-i+1 + 0.2, box_width = 0.3, 
            gp = gpar(fill = "red"), direction = "horizontal")
        grid.boxplot(m2[ index[i], ], pos = nr-i+1 - 0.2, box_width = 0.3, 
            gp = gpar(fill = "green"), direction = "horizontal")
    }
    grid.xaxis()
    popViewport()
}

ht_list = ht_list + rowAnnotation(boxplot = anno_multiple_boxplot, width = unit(4, "cm"), 
    show_annotation_name = FALSE)
lgd = Legend(labels = c("m1", "m2"), title = "boxplots",
    legend_gp = gpar(fill = c("red", "green")))
draw(ht_list, padding = unit(c(20, 2, 2, 2), "mm"), heatmap_legend_list = list(lgd))

image

其他技巧

參考鏈接

1.為不同維度的不同熱圖設(shè)置相同的單元格大小

假設(shè)您有一個(gè) heatmaps/oncoPrints 列表要保存為不同的文件（例如 png 或 pdf ），您可能想要做的一件事是使熱圖中每個(gè)網(wǎng)格/單元格的大小在熱圖中相同，因此，您需要根據(jù)熱圖中的行數(shù)或列數(shù)計(jì)算png/pdf文件的大小。在ComplexHeatmap生成的熱圖中，所有熱圖組件都具有絕對(duì)大小，并且只有熱圖主體的大?。ɑ騿卧竦拇笮。┦强筛牡模ɑ蛘邠Q句話說(shuō)，如果您更改最終圖形設(shè)備的大小，例如通過(guò)拖動(dòng)圖形窗口，如果您繪制時(shí)，僅調(diào)整熱圖主體的大?。@意味著整個(gè)圖的大小與熱圖中的行數(shù)或列數(shù)呈線性相關(guān)。這意味著我們實(shí)際上可以擬合一個(gè)線性模型y = a*x + b，其中y是整個(gè)圖的高度，x是行數(shù)。

在下面的例子中，我們簡(jiǎn)單地演示了如何在熱圖中建立繪圖高度和行數(shù)之間的關(guān)系。我們首先定義一個(gè)函數(shù)，該函數(shù)生成一個(gè)具有特定行數(shù)的 10 列矩陣。請(qǐng)注意，矩陣中的值在本演示中并不重要。

random_mat = function(nr) {
    m = matrix(rnorm(10*nr), nc = 10)
    colnames(m) = letters[1:10]
    return(m)
}

由于是絕對(duì)線性的關(guān)系，我們只需要測(cè)試兩個(gè)具有不同行數(shù)的熱圖，其中單行的高度為unit(5, "mm")。在熱圖中，還有列標(biāo)題、列樹(shù)狀圖、列注釋和列名。

下面的代碼有幾點(diǎn)需要注意：

1. 熱圖對(duì)象應(yīng)該由draw()返回，因?yàn)闊釄D的布局是在執(zhí)行后才計(jì)算的draw()。
2. component_height()返回一個(gè)單位向量，這些單位對(duì)應(yīng)于熱圖中從上到下的所有熱圖組件的高度。（component_width()返回?zé)釄D組件的寬度）。
3. 在計(jì)算 時(shí)ht_height，我們添加unit(4, "mm")，因?yàn)樵谧罱K圖的頂部和底部，有2mm白色邊框。
4. ht_height需要轉(zhuǎn)換為cm或inch單位。

在下文中，y包含以inch單位測(cè)量的值。

y = NULL
for(nr in c(10, 20)) {
    ht = draw(Heatmap(random_mat(nr), height = unit(5, "mm")*nr, 
        column_title = "foo", # one line text
        top_annotation = HeatmapAnnotation(bar = 1:10)))
    ht_height = sum(component_height(ht)) + unit(4, "mm")
    ht_height = convertHeight(ht_height, "inch", valueOnly = TRUE)
    y = c(y, ht_height)
}

然后我們可以擬合y和行數(shù)之間的線性關(guān)系：

x = c(10, 20)
lm(y ~ x)

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      1.2222       0.1969

這意味著行數(shù)x和繪圖高度y之間的關(guān)系是：y = 0.1969*x + 1.3150。

您可以通過(guò)以下代碼測(cè)試不同行的熱圖的單行高度是否相同。請(qǐng)注意，所有熱圖配置都應(yīng)與您準(zhǔn)備的y相同。

for(nr in c(10, 20)) {
    png(paste0("test_heatmap_nr_", nr, ".png"), width = 5, height = 0.1969*nr + 1.3150, 
        units = "in", res = 100)
    draw(Heatmap(random_mat(nr), height = unit(5, "mm")*nr, 
        column_title = "foo", # column title can be any one-line string
        top_annotation = HeatmapAnnotation(bar = 1:10)))
    dev.off()
}

image

image