再论热图 (heatmap) 3

“`{r echo=FALSE,warning=FALSE,message=FALSE}
library(heatmap3)
library(Heatplus)
library(ComplexHeatmap)
library(fields)
library(d3heatmap)
library(heatmaply)
library(rbokeh)
“`

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

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

首先, 我们先生成一个随机数据。
“`{r}
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 ```{r} library(ComplexHeatmap) Heatmap(mat) ``` ### fields fields的特点就是图象中的每一个格的位置都是可以计算的，由此可以产生很多灵活性。下面的例子就是产生一个heatmap table. ```{r} 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) } } } ``` ### heatmap3 在heatmap, heatmap.plus以及gplots::heatmap.2之后，赵同学写了heatmap3，并且在BioC2015上做了介绍。heatmap3的特点是可以对heatmap的每一个单元格做出特殊的标记，并且对行和列加上2D的基本注释图像。 ```{r} 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"))) ``` ### HeatPlus HeatPlus的特点是和eSet数据类型结合得比较紧密，并且可以对行和列进行多种不同数据类型的注释，也就是说对行或者列可以绘制2D的基本注释图像，比如说曲线，点图，以及histograms等。 ```{r,warning=FALSE} 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)))) ``` ## 交互热图 interactive heatmap 对于交互热图，基本上都是基于Web浏览器的。当图像生成之后，可以使用网页浏览器对数据交互体验。 ### d3heatmap d3heatmap是大名顶顶的交互式热图包。 ```{r, eval=FALSE} library(d3heatmap) d3heatmap(mat) ``` ### heatmaply heatmaply是基于plotly的热图包。[plotly.js](https://plot.ly/javascript/) 是一个基于javascript的强大的交互绘图包。它推出了多种不同语言的API。 ```{r, eval=FALSE} library(heatmaply) heatmaply(mat) ``` ### Bokeh [Bokeh](http://bokeh.pydata.org/en/latest/)和plotly.js一样都是javascript的包。它也有多种接口。但是在风格上更加简洁。而且写作风格与ggplot2一致。不足之处在于，它开发的程度不足。 ```{r, eval=FALSE} library(rbokeh) (figure(padding_factor = 0) %>% ly_image(mat))
“`

## 详解ComplexHeatmap

接下来就是本文的重点了，如何使用ComplexHeatmap。ComplexHeatmap的作者给出了很多非常实用的实例：[bioinformatics](http://bioinformatics.oxfordjournals.org/content/early/2016/05/20/bioinformatics.btw313)以及[Bioconductor](https://bioconductor.org/packages/release/bioc/vignettes/ComplexHeatmap/inst/doc/s1.introduction.html)。本文的目的在于推广ComplexHeatmap的使用。

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


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

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

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

“`{r}
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) ``` 也可以在画图的过程中对图例进行调整。 ```{r} draw(ht_list, heatmap_legend_side = "left", show_annotation_legend = FALSE) ``` 从上面的步骤中，我们可以看到，使用Heatmap可以实例化一个Heatmap，使用HeatmapAnnotation可以实例化一个HeatmapAnnotation。这些实例可以通过'+'号连接起来。画图使用draw或者print命令(注意，因为我们在interactive模式下是自动调用show命令的，所以系统会直接画图，如果不是在interactive模式下，一定要使用draw或者pirnt命令来出图)。 当然annotation可以是多种图型，其中包括points, boxplot, barplot, density_line, violin, histogram, 甚至heatmap。 ```{r, fig.height=8} 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")) ``` 对于rowAnnotation，以上各种类型的注释，可以使用row_前缀，比如 * `row_anno_points()` * `row_anno_barplot()` * `row_anno_boxplot()` * `row_anno_histogram()` * `row_anno_density()` ```{r, fig.width=8, fig.height=8} 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")) }) ``` 如果想只对个别的行或者列标示名称，其它的不标注，如何实现呢？ ```{r} 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)) ``` ### OncoPrint cBioPortal给我们带来了lollipop plot以及OncoPrint来表示基因中的突变情况。这里的`oncoPrint`就是对OncoPrint的R包装。 ```{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 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) ``` 上面的例子有两个重要的概念，一个是`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设置一个全局变量, 画两次。也许有更好的办法，希望同学们留言交流。 ```{r} 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)
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) ``` 下面就是一个实际的例子，其实也就是ComplexHeatmap中的实例。 ```{r} 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] unique(unlist(strsplit(mat, ";"))) ``` 我们可以看到，这里有三种不同的突变类型，“MUT”，“AMP”，“HOMDEL”。我们要对它们分别写出`alter_fun`。 ```{r} 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)) } ) ``` 同时对我们使用到的颜色进行定义，它们将用在上方和右方的统计中。 ```{r} col = c("MUT" = "#008000", "AMP" = "red", "HOMDEL" = "blue") ``` 使用`remove_empty_columns`来去除没有突变的样品。 使用`column_order`以及`row_order`来按需要设定样品顺序。 ```{r} 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") ``` =====EndrmdSourceCode=====

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

library(ComplexHeatmap)
Heatmap(mat)

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)     
    }
  }
}

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")))

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))))

library(d3heatmap)
d3heatmap(mat)

library(heatmaply)
heatmaply(mat)

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

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)

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

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"))

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"))
})

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))

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)

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) 

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) 

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"

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")

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")