EnrichedHeatmap 教程

本文是《ATAC-Seq 分析流程》的补充，简单展示用 EnrichedHeatmap 生成 ATAC-Seq 在指定区域的信号富集情况，当然也适用于 CHIP-Seq 或组蛋白修饰的信号展示。

EnrichedHeatmap 是基于 ComplexHeatmap 的定制化热图，所以如果需要进行个性化的调整，先学习熟悉 ComplexHeatmap 比较适合。

先导入需要的 R 包。

library(TxDb.Hsapiens.UCSC.hg38.knownGene, quietly = TRUE)
library(EnrichedHeatmap, quietly = TRUE)
library(rtracklayer, quietly = TRUE)

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
bwPath <- file.path("/Example/WT1_ATAC.bw")

读取 bigwig 文件。

bws <- rtracklayer::import(bwPath, format = "BIGWIG")
> bws
GRanges object with 76967086 ranges and 1 metadata column:
               seqnames      ranges strand |     score
                  <Rle>   <IRanges>  <Rle> | <numeric>
         [1]       chr1     1-10040      * |   0.00000
         [2]       chr1 10041-10050      * |   1.90888
         [3]       chr1 10051-10060      * |   3.18146
         [4]       chr1 10061-10070      * |   3.81775
         [5]       chr1 10071-10110      * |   5.72663
         ...        ...         ...    ... .       ...
  [76967082] KI270757.1 70811-70900      * |  0.000000
  [76967083] KI270757.1 70901-70910      * |  0.636292
  [76967084] KI270757.1 70911-70980      * |  1.272580
  [76967085] KI270757.1 70981-70990      * |  0.636292
  [76967086] KI270757.1 70991-71251      * |  0.000000
  -------
  seqinfo: 194 sequences from an unspecified genome

将启动子区域数据转换成矩阵。

tss <- promoters(genes(txdb), upstream = 3000, downstream = 3000)
tssMat <- normalizeToMatrix(bws, tss, extend = 0, value_column = "score", k = 400, background = 0, mean_mode = "w0")

# 设置 k 为 400 所以有 400 windows
> tssMat
Normalize bws to tss:
  Upstream 0 bp (0 window)
  Downstream 0 bp (0 window)
  Include target regions (400 windows)
  25750 target regions

这里 value_column 是你要画图的信号值列，函数将用该列值转换成画图矩阵；k 参数是划分窗口的数目，将决定画图矩阵的列数，一定要设置到合适；mean_mode 参数有以下几种。

Following illustrates different settings for mean_mode (note there is one signal region overlapping with other signals):

      40      50     20     values in signal regions
    ++++++   +++    +++++   signal regions
           30               values in signal region
         ++++++             signal region
      =================     a window (17bp), there are 4bp not overlapping to any signal regions.
        4  6  3      3      overlap

    absolute: (40 + 30 + 50 + 20)/4
    weighted: (40*4 + 30*6 + 50*3 + 20*3)/(4 + 6 + 3 + 3)
    w0:       (40*4 + 30*6 + 50*3 + 20*3)/(4 + 6 + 3 + 3 + 4)
    coverage: (40*4 + 30*6 + 50*3 + 20*3)/17 

将数据转换成矩阵后用 EnrichedHeatmap 画图。

EnrichedHeatmap(tssMat)

效果如下。

数据表现比较差

参考
Make Enriched Heatmaps
Using EnrichedHeatmap for visualization of NGS experiments
An Introduction to the GenomicRanges Package