视频小教程_如何画出没有教程的图?

(本次操作有配套的视频教程,在果子学生信公众号回复“果子对角线”自行获取,和代码一起以project的形式分享)

果子唠嗑

这是一篇引起血战的帖子,一开始发表在Y叔的公众号上,然后就出现了
ggplot2版本,excel版本,complexheatmap版本,base plot版本。

这篇帖子实现的略显复杂,但是,这个帖子是我做事风格的见证:

用有限的工具,尽量去做无限的事情。就像钢琴的88键,键有数,律无穷。

比如,我只熟悉R语言,那么批量读取文件,数据调整,绘图都只用R语言完成,尽管有很多更好的工具。
比如,我只会ggplot2,所以当出现以下需求的时候,我都尽量用ggplot2去完成,因为没有选择。

这篇帖子最终不是教大家如何画一张图,画一张图没有意义,而是展示R语言里面的一些基本技能,包括

  • 如何创建数据框?
  • 如何融合多个数据框?
  • 如何实现数据的宽长转换?
  • 如何调整x轴label的顺序?
  • 多张图如何拼接?
  • ggplot2怎么用?
  • 多边形怎么画?
  • 如何把大需求拆解成可完成的小步骤?
  • R语言中批量操作神器lapply怎么用?

既然有这么多技能,估计初学者也消化不了,所以我就变废为宝,录了一个视频教程,来讲解这些技能。
教程以R语言project的形式发布,里面包含了原始数据,代码,生成的图片,以及手把手教程。
有需要的在果子学生信公众号回复"果子对角线"自助获取。

以下是正文

最近有人在群里提问,下面的热图该怎么画。

这张热图,在每一个格子里面用对角线一分为二,呈现了两类信息。图片来自于这篇文章的图1D。

通常,遇到这种需求,我都会直接检索,看看有没有现成的R包可用。可以提取检索词检索一下,这里主要信息应该是"对角线 热图 R语言",
查阅字典后,关键词就变成了英文,“diagonal heatmap r”。在浏览器中检索后得到以下结果。

然后有人给出了使用base plot画图的方案

d=data.frame(p1=rep(LETTERS[1:2],times=2,each=2),
             p2=rep(LETTERS[4:5],times=4),
             value=c(.2,.4,.56,.32,.7,.16,.12,.71),
             group=rep(c("A1","B1"),each=4))

x=as.numeric(d$p1)
y=as.numeric(d$p2)

plot(1,xlim=c(1,length(unique(x))+1),ylim=c(1,length(unique(y))+1),
  type="n",bty="n",xaxt="n",yaxt="n",xlab="",ylab="")

for(i in 1:nrow(d)) {
  if(d$group[i]=="A1") polygon(x[i]+c(0,1,1),y[i]+c(0,0,1),col=gray(d$value[i]))
  if(d$group[i]=="B1") polygon(x[i]+c(0,1,0),y[i]+c(0,1,1),col=gray(d$value[i]))
}

axis(1,at=sort(unique(x))+.5,labels=levels(d$p1),lty=0)
axis(2,at=sort(unique(y))+.5,labels=levels(d$p2),lty=0)

复制到R语言后运行得到如下结果

这意味着,只要用基本功扎实,复现原图是没有问题的。
大鱼海棠应该就是用的这个方法。
base plot版本的对角线热图
但是,我不会base plot,我只会用ggplot2, 检索了一圈也没有现成的解决方案,所以我就准备自己画。

自己画我就开始拆分这张图,一个热图的每一个格子由上下两个三角形构成,我可以先画上三角形,再画下三角形,然后给他们分别配色就可以。
(果子:此时,我还不会写图层,尽管有这样的思想,但是没有趁手的工具)

三角形可以通过ggplot2中内置的geom_polygon来画。

## 下三角数据
cone1 = data.frame(x = c(1,2,2),
                   y = c(1,1,2))
library(dplyr)
library(ggplot2)
## 下三角画图
cone1 %>% 
  ggplot(aes(x=x, y=y)) +
  geom_polygon(fill="#213c18")

上三角也是一样的

## 上三角数据
cone2 = data.frame(x = c(1,1,2),
                   y = c(1,2,2)) 
## 上三角画图
cone2 %>% 
  ggplot(aes(x=x, y=y)) +
  geom_polygon(fill = "#668c6f")

通过叠加图层,我们可以把两个三角形放在一起,这就是热图的一个小格子

cone1 %>% 
  ggplot(aes(x=x, y=y)) +
  geom_polygon(fill="#213c18")+
  geom_polygon(data=cone2, fill = "#668c6f")

解决了基本的元素,我就可以批量产生数据。
下面写了一个函数,批量算出给定行和列的三角形数据。
这里用上了,我很喜欢的lapply和do.call函数

trianle <- function(a,b,type="up"){
  ## 单个下三角数据的函数
  trianle_down <- function(a,b){
    data.frame(x=c(0,1,1)+a,
               y= c(0,0,1)+b,
               group=paste0(a,"_",b),
               stringsAsFactors = F)
  }
  ## 单个上三角数据的函数
  trianle_up <- function(a,b){
    data.frame(x=c(0,0,1)+a,
               y= c(0,1,1)+b,
               group=paste0(a,"_",b),
               stringsAsFactors = F)
  }
  ### 批量产生上三角的数据
  if(type=="up"){
    data <- do.call(rbind,lapply(1:b, function(i){
      do.call(rbind,lapply(1:a,trianle_up,i))
    }))
  }
  ### 批量产生下三角的数据
  if(type=="down"){
    data <- do.call(rbind,lapply(1:b, function(i){
      do.call(rbind,lapply(1:a,trianle_down,i))
    }))
  }
  return(data)
}

因为我们的测试数据是33行,20列,所以用这个函数来产生对应个数的三角形

updata <- trianle(33,20,"up")
downdata <- trianle(33,20,"down")

updata的数据是这样的,每一个框线里面都是一个三角形的数据,group那一列是为了geom_polygon画图需要。

可以作图来测试一下上三角形

library(ggplot2)
updata %>% 
  ggplot(aes(x=x, y=y)) +
  geom_polygon(aes(group=group,fill=group))+
  theme(legend.position = "none")

下三角也是一样的

downdata %>% 
  ggplot(aes(x=x, y=y)) +
  geom_polygon(aes(group=group,fill=group))+
  theme(legend.position = "none")

现在我比较安心了,只要调整原始数据的格式,然后把他跟三角形的数据匹配起来就行。匹配就依靠group那一列

读入第一个数据, 行是基因,列是癌症类型

dd <- data.table::fread("easy_input_amp.txt",data.table = F)


接下来宽数据变长
这里用到的技能我也很喜欢,可以参考这个帖子
我喜欢的gather快要被淘汰了,迎面走来了更好用的宽长转换工具

library(dplyr)
library(tidyr)
data1 <- dd %>% 
  ## 数据变长
  pivot_longer(-1,
               names_to = "type",
               values_to = "exp") %>% 
  rename(gene=V1) 

现在的数据是这个样子的

如果想要改变未来热图基因的顺序,可以使用因子的水平来控制
这个技能在这里讲过
因子(factor)就像贤内助,让你始终分清主次,拨开云雾。

data1$gene <- factor(data1$gene,
                     levels = c("ALKBH5","FTO",
                                "HNRNPA2B1","HNRNPC", "IGF2BP1", "IGF2BP2","IGF2BP3","RBMX","YTHDC1","YTHDC2","YTHDF1","YTHDF2","YTHDF3",
                                "METTL14","METTL3","RBM15","RBM15B","VIRMA","WTAP","ZC3H13"),
                     ordered = F)

接下来就是给这个数据增加能够和三角形group匹配的列,
这一步的操作有trick,比如,因子转为数值,是需要先变成字符串再变成字符的,但是我们这里直接as.numeric可以获取他的levels数值

data1 <- data1 %>% 
  mutate(a=as.numeric(as.factor(type)),
         b=as.numeric(gene)) %>% 
  mutate(group = paste0(a,"_",b)) %>% 
  as.data.frame() %>% 
  inner_join(updata,by="group")

最终获取了完整的作图数据

有了数据,作图就很简单了

library(ggplot2)
data1 %>% 
  ggplot(aes(x=x, y=y)) +
  geom_polygon(aes(group=group,fill=exp))+
  theme(legend.position = "none")

第二个数据也是同样的处理

dd <- data.table::fread("easy_input_del.txt",data.table = F)
data2 <- dd %>% 
  pivot_longer(-1,
               names_to = "type",
               values_to = "del") %>% 
  rename(gene=V1)

data2$gene <- factor(data2$gene,
                     levels = c("ALKBH5","FTO",
                                "HNRNPA2B1","HNRNPC", "IGF2BP1", "IGF2BP2","IGF2BP3","RBMX","YTHDC1","YTHDC2","YTHDF1","YTHDF2","YTHDF3",
                                "METTL14","METTL3","RBM15","RBM15B","VIRMA","WTAP","ZC3H13"),
                     ordered = F)
data2 <- data2 %>% 
  mutate(a=as.numeric(as.factor(type)),
         b=as.numeric(as.factor(gene))) %>% 
  mutate(group = paste0(a,"_",b)) %>% 
  as.data.frame() %>% 
  inner_join(downdata,by="group")

data2 %>% 
  ggplot(aes(x=x, y=y)) +
  geom_polygon(aes(group=group,fill=del))+
  theme(legend.position = "none")


现在只要把两个图叠加在一起就行了,但是这里面有个问题。这两个图是分别用两个不同的变量来控制颜色的。
但是如果画在一起,ggplot2是不允许的,但是我记得Y叔发过的大部分帖子,我记得他曾经解决过这种情况,就去翻了一下。
答案在这个帖子里面
连续型和离散型数据一起画热图,外加分开配色和图例

这是我们常见的使用gheatmap做图的代码,为了实现不同数据类型的热图放在一起,让我们祭出ggnewscale包,只需要用new_scall_fill(),就可以把后面的fill映射和前面的分离开,然后我们就可以愉快地,再加一条gheatmap再画一个热图。

这样我的问题就完美解决了。
安装和加载 R包

options("repos"=c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
if(!require("ggnewscale")) install.packages("ggnewscale",update = F,ask = F)
library(ggnewscale)

分别抽离出基因名称和亚型名称,用来加在图上

## x轴文字
xlabels <- data1 %>% 
  filter(y==1) %>% 
  distinct(type,.keep_all = T) %>% 
  arrange(x) %>% 
  pull(type)
## y轴文字
ylabels <- data1 %>% 
  filter(x==1) %>% 
  distinct(gene,.keep_all = T) %>% 
  arrange(y) %>% 
  pull(gene) %>% 
  as.character()

指定三个颜色,来自小丫画图的144号作品

red  <- "#AB221F"
blue <- "#3878C1"
nake <- "#FFFADD"

画图开始

p1 <- ggplot(data1,aes(x=x, y=y)) +
  ## 画出上三角
  geom_polygon(data = data1,aes(x=x, y=y, group=group,fill=exp),color="black") +
  ## 上三角用表达值来配色
  scale_fill_gradientn(colours = c(nake,red)) +
  ## 神技能,清空aes
  new_scale_fill() +
  ## 画出下三角
  geom_polygon(data = data2,aes(x=x, y=y, group=group,fill=del),color="black") +
  ## 下三角用表达值来配色
  scale_fill_gradientn(colours = c(nake,blue))+
  ## 调整x轴
  scale_x_continuous(limits = c(1, NA), expand = c(0,0),breaks = c(1:33)+0.5,labels=xlabels) +
  ## 调整y轴
  scale_y_continuous(limits = c(1, NA), expand = c(0,0),breaks = c(1:20)+0.5,labels=ylabels) +
  ## 定主题
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.x = element_text(angle = 90, hjust = 1))+
  ## legend置于底部
  theme(legend.position="bottom")+
  theme(plot.margin = margin(0.5,0.01,0.5,0.01, "cm"))
p1

现在来画侧边的标识,也是自己设计四边形的数据就行了

data <- data.frame(
  x = rep(c(1,1,2,2),3),
  y = c(1,3,3,1,3,14,14,3,14,21,21,14),
  type= rep(c("Eraser","Reader","Writer"),each=4)
)

p2 <- ggplot(data,aes(x=x, y=y)) +
  geom_polygon(aes(x=x, y=y, group=type,fill=type),color="black",alpha = 0.5)+
  scale_fill_manual(values = c(red,nake,blue))+
  ## 打标签
  geom_text(data=data.frame(x=1.5,y=c(2,8.5,17.5)),
            aes(label=c("Eraser","Reader","Writer")),
            angle = 90,
            size=4)+
  scale_x_continuous(limits = c(1, NA), expand = c(0,0))+
  scale_y_continuous(limits = c(1, NA), expand = c(0,0)) +
  theme(axis.title=element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())+
  theme(legend.position="none")+
  theme(plot.margin = margin(0.5,0.1,0.5,0.01, "cm"))
p2

现在我们用cowplot把他们拼接起来

library(cowplot)
plot_grid(p1,p2,align = "h", axis = "tb",nrow = 1, rel_widths = c(33, 1))

如果用patchwork就更加简单

options("repos"=c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
if(!require("patchwork")) install.packages("patchwork",update = F,ask = F)
library(patchwork)
p1 + p2 +
  plot_layout(widths = c(33, 1))

虽然解决了问题,但是我觉得还是不够优雅,我还在学习图层中。Y叔有一篇帖子是让我夜不能寐的
扪心自问,meme几何?
这一年多来,很多次都想去回答这里面的问题,但是次次无功而返。好在近期洲更给了初步的回答,
答: 扪心自问,meme几何?

洲更说他对于图层顿悟了,就跟我问Y叔,他给出的答案是一样的,顿悟。但是我还没有,我还在找寻属于我自己的答案。学习图层,对我个人而言并不是为了画图,我就那么点数据,以后的知识以及足够去呈现了,我只是想从高手的视角去看待作图这件事,就跟Y叔给我启发最大的话那样:

图就是数据,数据就是图。

果子备注

这还没玩,下一篇帖子就讲一下我在画图上的突破。