EnrichedHeatmap 之基本用法

Original JunJunLab 老俊俊的生信笔记 2022-08-15

收录于合集

Do something meaningful,man?

引言

EnrichedHeatmap 也是 顾博士 开发的 R 包，用来可视化基因组富集信号的。

富集热图 是一种特殊类型的热图，它可视化富集特定的目标区域的基因组信号。它被广泛用于可视化，例如，组蛋白修饰如何在转录起始位点富集。

已经有一些工具能够制作这样的热图(例如 ngs.plot、deepTools 或 genomation)。在这里，通过 ComplexHeatmap 包实现了丰富的热图。由于富集热图本质上是一个普通的热图，但是有一些特殊的设置，使用 ComplexHeatmap 的功能，定制热图以及连接到一个热图列表来显示不同数据源之间的对应关系会容易得多。

先放张图养养眼：

基本用法

安装：

BiocManager::install("EnrichedHeatmap")

library(EnrichedHeatmap)

# 加载测试数据
set.seed(123)
load(system.file("extdata", "chr21_test_data.RData", package = "EnrichedHeatmap"))
ls()
[1] "cgi"     "genes"   "H3K4me3" "meth"    "rpkm"

H3K4me3: coverage for H3K4me3 histone modification from the ChIP-seq data。
cgi: CpG islands。
genes: genes。
meth: methylation for CpG sites from WGBS。
rpkm: gene expression from RNASeq。

可视化 H3K4me3 组蛋白修饰是在 TSS 基因周围的富集，首先提取基因的转录起始位点（TSS）信息：

tss = promoters(genes, upstream = 0, downstream = 1)
tss[1:5]
GRanges object with 5 ranges and 0 metadata columns:
                     seqnames    ranges strand
                        <Rle> <IRanges>  <Rle>
   ENSG00000141956.9    chr21  43299591      -
  ENSG00000141959.12    chr21  45719934      +
   ENSG00000142149.4    chr21  33245628      +
  ENSG00000142156.10    chr21  47401651      +
   ENSG00000142166.8    chr21  34696734      +
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

H3K4me3 的覆盖度：

# H3K4me3的覆盖度
H3K4me3[1:5]
GRanges object with 5 ranges and 1 metadata column:
      seqnames          ranges strand |  coverage
         <Rle>       <IRanges>  <Rle> | <integer>
  [1]    chr21 9825468-9825470      * |        10
  [2]    chr21 9825471-9825488      * |        13
  [3]    chr21         9825489      * |        12
  [4]    chr21         9825490      * |        13
  [5]    chr21 9825491-9825493      * |        14
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

可视化主要分为两步：

1、通过归一化到矩阵得到基因组信号和目标之间的关联。
2、使用热图可视化矩阵。

可视化：

mat1 = normalizeToMatrix(H3K4me3, tss, value_column = "coverage",
                         extend = 5000, mean_mode = "w0", w = 50)
mat1
Normalize H3K4me3 to tss:
  Upstream 5000 bp (100 windows)
  Downstream 5000 bp (100 windows)
  Include target regions (width = 1)
  720 target regions
# 查看类型
class(mat1)
[1] "normalizedMatrix" "matrix"

normalizeToMatrix() 将基因组信号（H3K4me3）和目标（tss）之间的关联转换为一个矩阵（实际上 mat1 只是一个带有几个附加属性的普通矩阵）。它首先将扩展的目标区域（向上游和下游的扩展由 extend 参数控制）分割成一个小窗口列表（窗口的宽度由 w 控制），然后将基因组信号重叠到这些小窗口上，计算每个小窗口的值，即与该窗口相交的基因组信号的平均值（该值对应于基因组信号）。

根据不同类型的基因组信号，mean_mode 有几种模式。后面的部分将对此进行解释。

EnrichedHeatmap(mat1, name = "H3K4me3")

默认情况下，行根据目标区域的丰富程度进行 行排序。在热图顶部有一种特定类型的注释，使用 anno_enriched() 绘制为线图。

EnrichedHeatmap() 内部调用 Heatmap() 并返回一个 Heatmap 类对象，因此 Heatmap() 的参数可以直接应用于 EnrichedHeatmap()。用户可以到 ComplexHeatmap 包获取更全面的帮助。

颜色：

与普通的热图类似，设置颜色的最简单方法是提供一个颜色向量：

EnrichedHeatmap(mat1, col = c("white", "red"), name = "H3K4me3")

你可能想知道为什么颜色看起来这么浅。原因是 H3K4me3 的覆盖值中 存在一些极值，导致 mat1 中出现极值。

查看数据分布：

quantile(H3K4me3$coverage, c(0, 0.25, 0.5, 0.75, 0.99, 1))
  0%  25%  50%  75%  99% 100%
  10   18   29   43   87  293

quantile(mat1, c(0, 0.25, 0.5, 0.75, 0.99, 1))
       0%       25%       50%       75%       99%      100%
  0.00000   0.00000   0.00000   0.00000  38.92176 174.78000

如果指定了颜色向量，则从最小到最大的顺序值将映射到颜色，而其他值则进行线性插值。要摆脱这种极端值，有两种方法。第一个是指定保留选项，该选项在下限和上限处修剪极值。（在下文中，这意味着只修剪大于第 99 个百分位数的值）

mat1_trim = normalizeToMatrix(H3K4me3, tss, value_column = "coverage",
    extend = 5000, mean_mode = "w0", w = 50, keep = c(0, 0.99))

# keep参数
EnrichedHeatmap(mat1_trim, col = c("white", "red"), name = "H3K4me3")

第二种方法是定义颜色映射函数，该函数只将颜色映射到小于第 99 个百分位数的值，大于第 99 个百分位数的值使用与第 99 个百分位数相同的颜色。

library(circlize)
col_fun = colorRamp2(quantile(mat1, c(0, 0.99)), c("white", "red"))

EnrichedHeatmap(mat1, col = col_fun, name = "H3K4me3")

按行分割：

通过指定 row_split 选项按向量或数据框拆分行：

mat1 = mat1_trim
EnrichedHeatmap(mat1, col = col_fun, name = "H3K4me3",
                row_split = sample(c("A", "B"), length(genes), replace = TRUE),
                column_title = "Enrichment of H3K4me3")

通过指定 row_km 选项，通过 k 均值聚类 拆分行：

set.seed(123)
EnrichedHeatmap(mat1, col = col_fun, name = "H3K4me3",
                row_km = 3,
                column_title = "Enrichment of H3K4me3",
                row_title_rot = 0)

当行被分割时，注释的图形参数可以是一个长度和行相等的向量：

# 改变线颜色和线形
set.seed(123)
EnrichedHeatmap(mat1, col = col_fun, name = "H3K4me3",
                row_km = 3,
                top_annotation = HeatmapAnnotation(enriched = anno_enriched(gp = gpar(col = 2:4, lty = 1:3))),
                column_title = "Enrichment of H3K4me3", row_title_rot = 0)

对行进行聚类：

默认情况下，show_row_dend 处于关闭状态，因此无需在此处指定它。更多关于行聚类的选项可以在 Heatmap() 的帮助页面中找到：

EnrichedHeatmap(mat1, col = col_fun,
                name = "H3K4me3",
                cluster_rows = TRUE,
                column_title = "Enrichment of H3K4me3")

延伸目标区域：

可以通过延伸单个值或长度为 2 的向量来控制对上游和下游的扩展：

# upstream 1kb, downstream 2kb
mat12 = normalizeToMatrix(H3K4me3, tss,
                          value_column = "coverage",
                          extend = c(1000, 2000),
                          mean_mode = "w0", w = 50)

EnrichedHeatmap(mat12, name = "H3K4me3", col = col_fun)

上游或下游也可以设置为 0：

mat12 = normalizeToMatrix(H3K4me3, tss,
                          value_column = "coverage",
                          extend = c(0, 2000),
                          mean_mode = "w0", w = 50)

EnrichedHeatmap(mat12, name = "H3K4me3", col = col_fun)

有时信号是在目标的上游还是下游并不重要，用户只想显示与目标的相对距离。如果 flip_upstream 设置为 TRUE，则归一化矩阵中的上游部分被翻转并添加到下游部分。仅当目标是单点目标或目标被排除在归一化矩阵中时才允许翻转（通过设置 include_target = FALSE）。如果上游和下游的扩展长度不相等，则较小的扩展长度用于最终矩阵：

mat_f = normalizeToMatrix(H3K4me3, tss,
                          value_column = "coverage",
                          extend = 5000,
                          mean_mode = "w0", w = 50,
                          flip_upstream = TRUE)
mat_f
Normalize H3K4me3 to tss:
  Extension 5000 bp (100 window)
    upstream is flipped to downstream.
  Include target regions (width = 1)
  720 target regions

EnrichedHeatmap(mat_f, name = "H3K4me3", col = col_fun)

富集注释的坐标轴：

注释的轴由 anno_enriched() 中的 axis_param 控制。所有可以设置的参数都可以在 ComplexHeatmap::default_axis_param() 中找到：

EnrichedHeatmap(mat1, col = col_fun, name = "H3K4me3",
    top_annotation = HeatmapAnnotation(
        enriched = anno_enriched(
            ylim = c(0, 10),
            axis_param = list(
                at = c(0, 5, 10),
                labels = c("zero", "five", "ten"),
                side = "left",
                facing = "outside"
    )))
)

均值模式：

当将基因组信号标准化为目标区域时，目标的上游和下游（如果包括在内，也包括目标区域本身）被分成小窗口。然后将基因组信号重叠到每个窗口并计算每个窗口的平均信号。当一个窗口没有被基因组信号的区域完全覆盖时，应该应用适当的平均方法来总结窗口中的值。

根据不同的场景，EnrichedHeatmap 提供了 三个指标 进行平均。

重叠模型如下图所示。底部的 红线代表小窗口。顶部的黑线是与窗口重叠的基因组信号区域。粗线表示信号区域和窗口之间的相交部分：

下面说明了 mean_mode 的不同设置（注意有一个信号区域与其他信号区域重叠）：

平滑化：

生成矩阵时，可以通过将 smooth 设置为 TRUE 来平滑化。稍后将演示平滑也有助于估算 NA 值。

由于平滑可能会改变原始数据范围，因此这里的颜色映射函数 col_fun 确保调色板仍然与未平滑的调色板相同。

背景对应于没有信号重叠的区域。适当的值取决于特定的场景。在这里，由于可视化了 ChIP-Seq 数据的覆盖范围，因此将 0 分配给没有 H3K4me3 信号的区域是合理的。

在以下示例中，由于富集热图也是热图，我们可以通过+连接两个 heamtaps：

mat1_smoothed = normalizeToMatrix(H3K4me3, tss,
                                  value_column = "coverage",
                                  extend = 5000, mean_mode = "w0",
                                  w = 50, background = 0, smooth = TRUE)

EnrichedHeatmap(mat1_smoothed, col = col_fun,
                name = "H3K4me3_smoothed",
                column_title = "smoothed") +
EnrichedHeatmap(mat1, col = col_fun, name = "H3K4me3",
                column_title = "unsmoothed")

在上面的图中，可能会觉得左侧的热图比右侧的未平滑热图更好。下面，我们将证明平滑可以显着改善 甲基化数据集的富集模式。

以下热图可视化了 TSS 上低甲基化区域的富集。灰色代表没有 CpG 位点的窗口（注意我们将 NA 设置为背景，灰色是 ComplexHeatmap 的 NA 值的默认颜色）：

meth[1:5]
GRanges object with 5 ranges and 1 metadata column:
      seqnames    ranges strand |      meth
         <Rle> <IRanges>  <Rle> | <numeric>
  [1]    chr21   9432427      * |  0.267104
  [2]    chr21   9432428      * |  0.267107
  [3]    chr21   9432964      * |  0.272710
  [4]    chr21   9432965      * |  0.272735
  [5]    chr21   9433315      * |  0.285115
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

mat2 = normalizeToMatrix(meth, tss, value_column = "meth",
                         mean_mode = "absolute",
                         extend = 5000, w = 50,
                         background = NA)

meth_col_fun = colorRamp2(c(0, 0.5, 1), c("blue", "white", "red"))

EnrichedHeatmap(mat2, col = meth_col_fun,
                name = "methylation",
                column_title = "methylation near TSS")

当将 CpG 位置重叠到分割的目标区域时，可能在某些窗口中没有 CpG 位点，特别是对于仅包含 100000 个 CpG 位点的 meth 在 21 号染色体中随机采样。这些不包含 CpG 位点的窗口的值可以通过平滑进行估算。虽然将甲基化值分配给非 CpG 窗口似乎不合适，但它会大大改善可视化：

mat2 = normalizeToMatrix(meth, tss, value_column = "meth",
                         mean_mode = "absolute",
                         extend = 5000, w = 50,
                         background = NA,
                         smooth = TRUE)
EnrichedHeatmap(mat2, col = meth_col_fun,
                name = "methylation",
                column_title = "methylation near TSS")

为了进行平滑，默认情况下，locfit() 首先应用于原始矩阵中的每一行。如果失败，则随后应用 loess() 平滑。如果两种平滑方法都失败，则会出现警告并保留原始值。

用户可以通过 smooth_fun 参数提供自己的平滑功能。这个自定义函数接受一个数值向量（可能包含 NA 值）并返回一个具有相同长度的向量。如果平滑失败，函数应该调用 stop() 来抛出错误，以便 normalizeToMatrix() 可以捕获平滑失败的行数。查看 default_smooth_fun() 的源代码以获取示例。

目标为范围区域：

在 H3K4me3 的例子中，目标区域是单点。目标也可以是 宽度大于 1 的区域。以下热图可视化了 CpG 岛上低甲基化的富集：

mat3 = normalizeToMatrix(meth, cgi, value_column = "meth",
                         mean_mode = "absolute",
                         extend = 5000, w = 50,
                         background = NA,
                         smooth = TRUE,
                         target_ratio = 0.3)
EnrichedHeatmap(mat3, col = meth_col_fun,
                name = "methylation",
                axis_name_rot = 90,
                column_title = "methylation near CGI")

热图上显示的 目标区域的宽度 可以由 target_ratio 控制，它相对于完整热图的宽度。

目标区域也被分成小窗口。由于目标区域的宽度不等，每个目标被分成 k 个相等的窗口，其中 (k = (n_1 +n_2)*r/(1-r)) 其中 n_1 在上游窗口中的数量，n_2 是下游窗口的数量，r 是矩阵中目标列的比率。normalizeToMatrix() 中有一个 k 参数，但它仅在目标没有上游或下游时使用。

当基因组目标是区域时，可以在热图中排除上游和或下游：

mat3 = normalizeToMatrix(meth, cgi, value_column = "meth",
                         mean_mode = "absolute",
                         extend = c(0, 5000), w = 50,
                         background = NA, smooth = TRUE,
                         target_ratio = 0.5)
EnrichedHeatmap(mat3, col = meth_col_fun,
                name = "methylation",
                column_title = "methylation near CGI")

当没有上游和下游时，热图中的列数由 k 参数控制：

mat3 = normalizeToMatrix(meth, cgi, value_column = "meth",
                         mean_mode = "absolute",
                         extend = 0, k = 20, background = NA,
                         smooth = TRUE, target_ratio = 1)
EnrichedHeatmap(mat3, col = meth_col_fun,
                name = "methylation",
                column_title = "methylation near CGI")