再论热图 (heatmap)

~ admin

之前我在R绘图基础（四）热图 heatmap非常粗浅地讲述了一些关于热图的知识，利用示例代码让大家了解了如何使用R的众多工具来绘制热图。通过转载的情况来看，那应该是本博客最受欢迎的一篇日志了。但是那已经是五年前的一篇老黄历了。现在在R中，又涌现了许多优秀的heatmap的工具。应一些读者的要求，我再论热图，希望能介绍一下关于heatmap的一些新的包。因为本人能力有限，难免挂一漏万，如有遗漏，希望大家留言补充。

首先要肯定华人在R数据图像化方面巨大贡献，尤其是在heatmap的高级应用中，华人应该是主要力量。本文介绍的几个package中，大多数都为华人写就。

首先, 我们先生成一个随机数据。

set.seed(1)
mat <- matrix(runif(36), nrow=6, 
              dimnames=list(paste("row", letters[1:6], sep="_"), 
                            paste("col", LETTERS[1:6], sep="_")))

2D热图

2D热图是为了同交互热图相区分而起的名字。除了之前介绍的软件包以外，现在在R当中出现了很多新的软件包，比如heatmap3, ComplexHeatmap, fields。在这里，我首先将这些包都跑一遍，让大家看看每个包在生成图象上的不同，最后重点介绍ComplexHeatmap包。

我先按字母排序来试用每一个包。

ComplexHeatmap

library(ComplexHeatmap)
Heatmap(mat)

plot of chunk unnamed-chunk-3

fields

fields的特点就是图象中的每一个格的位置都是可以计算的，由此可以产生很多灵活性。下面的例子就是产生一个heatmap table.

library(fields)
image.plot(mat)
# calculate the lowest and highest value
min.z <- min(mat)
max.z <- max(mat)
# 黄色和绿色范围内文字要用黑色
z.yellows <- min.z + (max.z - min.z)/64*c(20,45) 
# 每行每列加上数值
for(i in 1:nrow(mat)){
  for(j in 1:ncol(mat)){
    if((mat[i,j] > z.yellows[1])&(mat[i,j] < z.yellows[2])){
      text((i-1)/(nrow(mat)-1),
           (j-1)/(ncol(mat)-1),
           formatC(mat[i,j], width=4), 
           col="black", cex = 0.8)
    }else{
      text((i-1)/(nrow(mat)-1),
           (j-1)/(ncol(mat)-1),
           formatC(mat[i,j], width=4), 
           col="white", cex = 0.8)     
    }
  }
}

plot of chunk unnamed-chunk-4

heatmap3

在heatmap, heatmap.plus以及gplots::heatmap.2之后，赵同学写了heatmap3，并且在BioC2015上做了介绍。heatmap3的特点是可以对heatmap的每一个单元格做出特殊的标记，并且对行和列加上2D的基本注释图像。

library(heatmap3)
colorCell<-data.frame(row=c(2, 4, 6),
                      col=c(1, 3, 5),
                      color=c("cyan","black","blue"),
                      stringsAsFactors=F)
highlightCell <- data.frame(row=c(2, 4, 6),
                            col=c(1, 3, 5),
                            color=c("black","white","white"),
                            lwd=1:3, 
                            stringsAsFactors=F)
ColSideAnn <- data.frame(Information=rnorm(6),
                         Group=c(rep("Ctr", 2),
                                 rep(c("TrtA", "TrtB"),2)))
row.names(ColSideAnn) <- colnames(mat)
RowSideColors <- colorRampPalette(c("chartreuse4",
                                    "white", 
                                    "firebrick"))(nrow(mat))
heatmap3(mat,
         colorCell=colorCell,
         highlightCell=highlightCell,
         ColSideCut=1.2,
         ColSideAnn=ColSideAnn,
         ColSideFun=function(x) showAnn(x),
         ColSideWidth=0.8,
         RowSideColors=RowSideColors,
         col=colorRampPalette(c("#0F7D12", 
                                "#FFFD38", 
                                "#FC0D1B"))(10),
         RowAxisColors=1,
         legendfun=function() 
           showLegend(legend=c("High", "Low"),
                      col=c("#FC0D1B","#0F7D12")))

plot of chunk unnamed-chunk-5

HeatPlus

HeatPlus的特点是和eSet数据类型结合得比较紧密，并且可以对行和列进行多种不同数据类型的注释，也就是说对行或者列可以绘制2D的基本注释图像，比如说曲线，点图，以及histograms等。

library(Heatplus)
annoRow <- data.frame(means=rowMeans(mat), 
                      logic=rowMeans(mat)>.5)
annoCol <- data.frame(sums=colSums(mat), 
                      logic=sample(c(T, F), ncol(mat), 
                                   replace = TRUE))
plot(annHeatmap2(mat, 
                ann=list(Col=list(data=annoCol),
                         Row=list(data=annoRow))))

plot of chunk unnamed-chunk-6

交互热图 interactive heatmap

对于交互热图，基本上都是基于Web浏览器的。当图像生成之后，可以使用网页浏览器对数据交互体验。

d3heatmap

d3heatmap是大名顶顶的交互式热图包。

library(d3heatmap)
d3heatmap(mat)

heatmaply

heatmaply是基于plotly的热图包。plotly.js 是一个基于javascript的强大的交互绘图包。它推出了多种不同语言的API。

library(heatmaply)
heatmaply(mat)

Bokeh

Bokeh和plotly.js一样都是javascript的包。它也有多种接口。但是在风格上更加简洁。而且写作风格与ggplot2一致。不足之处在于，它开发的程度不足。

library(rbokeh)
(figure(padding_factor = 0) %>% ly_image(mat))

详解ComplexHeatmap

接下来就是本文的重点了，如何使用ComplexHeatmap。ComplexHeatmap的作者给出了很多非常实用的实例：bioinformatics以及Bioconductor。本文的目的在于推广ComplexHeatmap的使用。

ComplexHeatmap设计了两个重要的类，一个是Heatmap，一个是HeatmapAnnotation。多个Heatmap类还可以组成HeatmapList，HeatmapList可以将多个heatmap肩并肩平行放置。

ComplexHeatmap的布局设计

在绘图上ComplexHeatmap使用了与ggplot2类似的加法运算符(‘+’)，但是其背后使用的是grid包。因此作者特意保留了所有的grid中的viewport，所以可以对每个部分进行多层修饰。但是如果用户没有使用过grid的话，可能会比较难理解。

下面我们举两个例子，来讲解如何使用ComplexHeatmap来绘制复杂的热图。

Heatmap和annotation

我尽量使用ComplexHeatmap中原有的例子来讲解。

ha_column <- 
  HeatmapAnnotation(
    df = data.frame(## data.frame, 行数与mat的列数一致。
      type1 = sample(c("a", "b"), 
                     ncol(mat), 
                     replace = TRUE)),
    col = list(type1 = c("a" = "red", "b" = "blue")) 
    ## 颜色list, 名字与df中的列字一致，
    ##           元素名与df列中的出现的元素一致。
)
ha_row <- rowAnnotation( ## 对行进行注释
  df = data.frame(type2 = sample(c("A", "B"),
                                 nrow(mat), 
                                 replace = TRUE)),
  col = list(type2 = c("A" = "green", "B" = "cyan")),
  width = unit(1, "cm") ##指定宽度，必须是unit，参见grid::unit
)

## 生成两个Heatmap实例。这里我们使用了相同的数据。
## 在实际应用中，它应该是行数一致的不同数据。
ht1 <- Heatmap(mat, name = "ht1", #name将出现的图例中
               column_title = "Heatmap 1", #title将出现在标题中
               top_annotation = ha_column) #ha_column将出现在热图上方（也可以是下方）
ht2 <- Heatmap(mat, name = "ht2", 
               column_title = "Heatmap 2")
ht_list <- ht1 + ht2 + ha_row
draw(ht_list)

plot of chunk unnamed-chunk-10

也可以在画图的过程中对图例进行调整。

draw(ht_list, heatmap_legend_side = "left", 
     show_annotation_legend = FALSE)

plot of chunk unnamed-chunk-11

从上面的步骤中，我们可以看到，使用Heatmap可以实例化一个Heatmap，使用HeatmapAnnotation可以实例化一个HeatmapAnnotation。这些实例可以通过’+’号连接起来。画图使用draw或者print命令(注意，因为我们在interactive模式下是自动调用show命令的，所以系统会直接画图，如果不是在interactive模式下，一定要使用draw或者pirnt命令来出图)。

当然annotation可以是多种图型，其中包括points, boxplot, barplot, density_line, violin, histogram, 甚至heatmap。

ha_mix_top <- HeatmapAnnotation(
  points = anno_points(runif(ncol(mat))), ## points
  barplot = anno_barplot(runif(ncol(mat)), axis=TRUE), ## barplot
  histogram = anno_histogram(mat, 
                             gp = gpar(fill = 1:6)), ## histogram
  density_line = anno_density(mat, type = "line", 
                              gp = gpar(col = 1:6)), ## density
  violin = anno_density(mat, type = "violin", gp = gpar(fill = 6:1)), ## violin
  box = anno_boxplot(mat), ## boxplot
  heatmap = anno_density(mat, type = "heatmap") ## heatmap
)
Heatmap(mat, name = "mat", 
        top_annotation = ha_mix_top, 
        top_annotation_height = unit(8, "cm"))

plot of chunk unnamed-chunk-12

对于rowAnnotation，以上各种类型的注释，可以使用row_前缀，比如 * row_anno_points() * row_anno_barplot() * row_anno_boxplot() * row_anno_histogram() * row_anno_density()

ha1 <- rowAnnotation(b1 = row_anno_boxplot(mat, axis = TRUE),
    p1 = row_anno_points(runif(nrow(mat)), axis = TRUE),
    width = unit(4, "cm"))
ha2 <- rowAnnotation(b2 = row_anno_boxplot(mat, axis = TRUE),
    p2 = row_anno_points(rowMeans(mat), axis = TRUE), 
    width = unit(4, "cm"))
Heatmap(mat, name = "foo", top_annotation = ha_mix_top, 
        top_annotation_height = unit(8, "cm")) + ha1 + ha2

#在现有注释的基础上，可以使用`decorate_annotation`来对图像近一步修饰。

decorate_annotation("barplot", ## 第一个参数是要找的annotation的名字
                    {##第二个参数是要运行的命令。
  grid.text("barplot", unit(-10, "mm"), just="bottom", rot=90)
  grid.lines(c(0, 1), unit(c(.5, .5), "native"), ##这里的native的意思是指实际坐标系。 
             gp=gpar(lty=2, col="blue"))
})

decorate_annotation("p2", {
  grid.lines(unit(c(.5, .5), "native"), c(0, 1),
             gp=gpar(col="red"))
})

plot of chunk unnamed-chunk-13

如果想只对个别的行或者列标示名称，其它的不标注，如何实现呢？

mat = matrix(rnorm(10000), nr = 1000)
rownames(mat) = sprintf("%.2f", rowMeans(mat))
subset = sample(1000, 20)
labels = rownames(mat)[subset]
Heatmap(mat, show_row_names = FALSE, show_row_dend = FALSE, show_column_dend = FALSE) + 
rowAnnotation(link = row_anno_link(at = subset, labels = labels),
  width = unit(1, "cm") + max_text_width(labels))

plot of chunk unnamed-chunk-14

OncoPrint

cBioPortal给我们带来了lollipop plot以及OncoPrint来表示基因中的突变情况。这里的oncoPrint就是对OncoPrint的R包装。

mat <- read.table(textConnection(
",s1,s2,s3
g1,snv1;indel,snv1;snv2,indel
g2,,snv2;indel,snv1
g3,snv1,,indel;snv1;snv2"), 
row.names = 1, header = TRUE, 
sep = ",", stringsAsFactors = FALSE)
mat <- as.matrix(mat)
mat
##    s1           s2           s3               
## g1 "snv1;indel" "snv1;snv2"  "indel"          
## g2 ""           "snv2;indel" "snv1"           
## g3 "snv1"       ""           "indel;snv1;snv2"
library(ComplexHeatmap)
col <- c(snv1 = "#008000", snv2 = "#008000", indel = "blue")
oncoPrint(mat, get_type = function(x) strsplit(x, ";")[[1]],
    alter_fun = list(
      indel = function(x, y, w, h) 
        grid.rect(x, y, w*0.9, h*0.9, 
                  gp = gpar(fill = col["indel"], 
                            col = NA)) , #must be draw first
        snv1 = function(x, y, w, h) 
          grid.rect(x, y-.15*h, w*0.9, h*0.2, 
                    gp = gpar(fill = col["snv1"], 
                              col = NA)),
        snv2 = function(x, y, w, h) 
          grid.rect(x, y+.15*h, w*0.9, h*0.2, 
                    gp = gpar(fill = col["snv2"], 
                              col = NA))
    ), col = col)

plot of chunk unnamed-chunk-15

上面的例子有两个重要的概念，一个是get_type，一个是alter_fun。两者都接受一个函数，其中get_type是如何将mat中的每个单元如何分解成突变类型，alter_fun是如何对每一种突变类型进行绘图。 get_type接受一个参数x，这个x将会是mat中的每一个值。比如这里的mat[g1, s1]就是“snv1;indel”。而alter_fun将会是一个list of function，list中的每个元素名称都将是可能出现在get_type得到的值中。list中的每一个函数有四个参数，x,y,w,h，它们分别指heatmap中对应mat中的每个单元格的x, y的坐标以及单元格的高度和宽度。

现在出一个问题：如果我们希望所有的snv都合并考虑，当只有一个snv的时候，画在cell的中间，如果有两个的时候，对称的画两条。也就是说，无论是snv1还是snv2，如果只出现其中一个，都画在正中间。这时应该怎么做呢？

我的想法是对每一个cell设置一个全局变量, 画两次。也许有更好的办法，希望同学们留言交流。

rm(list=ls(pattern="cell_.*_.*"))
max <- 0
snv <- function(x, y, w, h) {
  globalV <- paste("cell", x, y, sep="_")
  cnt <- 1
  if(exists(globalV)){
    cnt <- get(globalV, envir=globalenv())+1
  }
  assign(globalV, cnt, envir=globalenv())
  if(cnt>max) assign("max", cnt, envir=globalenv())
  grid.rect(x, y, w*0.9, h*0.2, 
            gp = gpar(fill = col["snv1"], 
                      col = NA))
}
oncoPrint(mat, get_type = function(x) strsplit(x, ";")[[1]],
    alter_fun = list(
      indel = function(x, y, w, h) 
        grid.rect(x, y, w*0.9, h*0.9, 
                  gp = gpar(fill = col["indel"], 
                            col = NA)) , #must be draw first
        snv1 = snv,
        snv2 = snv
    ), col = col)

plot of chunk unnamed-chunk-16

H.factor <- 1 / ( 2 * max - 1)
rm(list=ls(pattern="cell.cnt_.*_.*"))
snv <- function(x, y, w, h) {
  globalV <- paste("cell", x, y, sep="_")
  total <- get(globalV, envir=globalenv())
  localV <- paste("cell.cnt", x, y, sep="_")
  if(exists(localV)){
    cnt <- get(localV, envir=globalenv())+1
  }else{
    cnt <- 1
  }
  assign(localV, cnt, envir=globalenv())
  pos.y <- seq(from=0, to=1, length.out=total+1)
  pos.y <- (pos.y[-length(pos.y)] + pos.y[-1]) / 2
  pos.y <- (pos.y[cnt] - .5 ) * h
  grid.rect(x, y + pos.y, w*0.9, h*H.factor, 
            gp = gpar(fill = col["snv1"], 
                      col = NA))
}
oncoPrint(mat, get_type = function(x) strsplit(x, ";")[[1]],
    alter_fun = list(
      indel = function(x, y, w, h) 
        grid.rect(x, y, w*0.9, h*0.9, 
                  gp = gpar(fill = col["indel"], 
                            col = NA)) , #must be draw first
        snv1 = snv,
        snv2 = snv
    ), col = col)

plot of chunk unnamed-chunk-16

下面就是一个实际的例子，其实也就是ComplexHeatmap中的实例。

mat <- read.delim(file.path(system.file("extdata", "tcga_lung_adenocarcinoma_provisional_ras_raf_mek_jnk_signalling.txt",
                                       package = "ComplexHeatmap")), 
                 stringsAsFactors=FALSE, row.names=1)
mat[is.na(mat)] <- "  "
mat <-  mat[, -ncol(mat)]
mat = t(as.matrix(mat))
mat[1:3, 1:3]
##      TCGA-05-4384-01 TCGA-05-4390-01 TCGA-05-4425-01
## KRAS "  "            "MUT;"          "  "           
## HRAS "  "            "  "            "  "           
## BRAF "  "            "  "            "  "
unique(unlist(strsplit(mat, ";")))
## [1] "  "     "MUT"    "AMP"    "HOMDEL"

我们可以看到，这里有三种不同的突变类型，“MUT”，“AMP”，“HOMDEL”。我们要对它们分别写出alter_fun。

alter_fun <- list(
    background = function(x, y, w, h) {
        grid.rect(x, y, w-unit(0.5, "mm"), h-unit(0.5, "mm"), gp = gpar(fill = "#CCCCCC", col = NA))
    },
    HOMDEL = function(x, y, w, h) {
        grid.rect(x, y, w-unit(0.5, "mm"), h-unit(0.5, "mm"), gp = gpar(fill = "blue", col = NA))
    },
    AMP = function(x, y, w, h) {
        grid.rect(x, y, w-unit(0.5, "mm"), h-unit(0.5, "mm"), gp = gpar(fill = "red", col = NA))
    },
    MUT = function(x, y, w, h) {
        grid.rect(x, y, w-unit(0.5, "mm"), h*0.33, gp = gpar(fill = "#008000", col = NA))
    }
)

同时对我们使用到的颜色进行定义，它们将用在上方和右方的统计中。

col = c("MUT" = "#008000", "AMP" = "red", "HOMDEL" = "blue")

使用remove_empty_columns来去除没有突变的样品。 使用column_order以及row_order来按需要设定样品顺序。

sample_order <- scan(system.file("extdata", "sample_order.txt",
                                package = "ComplexHeatmap"), 
                    what = "character")
ht <- oncoPrint(mat, get_type = function(x) strsplit(x, ";")[[1]],
                alter_fun = alter_fun, col = col, 
                row_order = NULL, column_order = sample_order,
                remove_empty_columns = TRUE,
                column_title = "OncoPrint for TCGA Lung Adenocarcinoma, genes in Ras Raf MEK JNK signalling",
                heatmap_legend_param = 
                  list(title = "Alternations", 
                       at = c("AMP", "HOMDEL", "MUT"), 
                       labels = c("Amplification", 
                                  "Deep deletion", 
                                  "Mutation"),
                       nrow = 1, title_position = "leftcenter"))
draw(ht, heatmap_legend_side = "bottom")

plot of chunk unnamed-chunk-20

发布者：admin

3 thoughts on “再论热图 (heatmap)”

  2. 您好，我想问R语言的学习方法。
    因为我在网上搜索并没有得到比较好的建议，
    比如找到了某个包的document，但其中的说明并不能满足我的画图需要，
    或并不了解其中的含义。
    希望您给我一些建议，谢谢！

    回复

发表评论

电子邮件地址不会被公开。 必填项已用*标注