“`{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=====
厉害了,一个热图被玩出各种花样来了,佩服你们
您好,我想问R语言的学习方法。
因为我在网上搜索并没有得到比较好的建议,
比如找到了某个包的document,但其中的说明并不能满足我的画图需要,
或并不了解其中的含义。
希望您给我一些建议,谢谢!
熟能生巧。只有多做练习,才是学习的捷径。