跟着NatureMetabolism学作图：R语言ggplot2画热图展示基因表达量

论文

Independent phenotypic plasticity axes define distinct obesity sub-types

https://www.nature.com/articles/s42255-022-00629-2#Sec15

s42255-022-00629-2.pdf

论文中没有公开代码，但是所有作图数据都公开了，我们可以试着用论文中提供的数据模仿论文中的图

今天的推文重复一下论文中的Fig3a 热图展示差异表达基因的表达量

论文中提供的数据没有上调下调的分组，这里我就随便选择数据了，两个热图之间的空白通过拼图的方式来实现

部分示例数据

总共553个基因，前300个标记为上调，后253个标记为下调

读取数据

library(readr)
dat<-read_delim("data/20220921/fig3a.txt",delim = "\t")
dim(dat)
colnames(dat)

dat01<-dat[1:300,]
dim(dat01)
dat02<-dat[301:553,]
dim(dat02)

第一个热图

library(ggplot2)
library(latex2exp)
library(tidyverse)

library(RColorBrewer)
cols<-colorRampPalette(rev(brewer.pal(n = 7, name =
                                        "RdYlBu")))(100)

dat01 %>% 
  pivot_longer(-gene_id,
               names_to = "treatment",
               values_to = "expre") -> dat01.1

dat01.1$treatment<-factor(dat01.1$treatment,
                          levels = colnames(dat01))

ggplot(data=dat01.1,
       aes(x=treatment,y=gene_id))+
  geom_tile(aes(fill=expre))+
  theme(axis.text = element_blank(),
        axis.ticks = element_blank())+
  labs(x=NULL,y=TeX(r"(\overset${\textit{Nnat}{^+}{^/}{^-}{^p}-\textit{Heavy}}{upregulated}$)"))+
  scale_fill_gradientn(colours = cols,
                       breaks=c(-1.6,-0.6,0.5,1.4,2.3),
                       labels=c(-2,-1,0,1,2),
                       name=NULL)+
  guides(fill=guide_colorbar(barheight = 25))+
  theme(panel.background = element_blank())

第二个热图

dat02 %>% 
  pivot_longer(-gene_id,
               names_to = "treatment",
               values_to = "expre") -> dat02.1

dat02.1$treatment<-factor(dat02.1$treatment,
                          levels = colnames(dat01))

ggplot(data=dat02.1,
       aes(x=treatment,y=gene_id))+
  geom_tile(aes(fill=expre))+
  theme(axis.text = element_blank(),
        axis.ticks = element_blank())+
  labs(x=NULL,y=TeX(r"(\overset${\textit{Nnat}{^+}{^/}{^-}{^p}-\textit{Heavy}}{downregulated}$)"))+
  scale_fill_gradientn(colours = cols,
                       breaks=c(-1.6,-0.6,0.5,1.4,2.3),
                       labels=c(-2,-1,0,1,2),
                       name=NULL)+
  guides(fill=guide_colorbar(barheight = 25))+
  theme(panel.background = element_blank(),
        legend.position = "none")

顶部表示分组的信息，用柱形图来展示

ggplot()+
  geom_rect(aes(xmin=0.5,xmax=3.5,ymin=0.5,ymax=1.5),
            fill="#c1bebf")+
  geom_rect(aes(xmin=3.5,xmax=6.5,ymin=0.5,ymax=1.5),
            fill="#71cddd")+
  geom_rect(aes(xmin=6.5,xmax=9.5,ymin=0.5,ymax=1.5),
            fill="#3a54a3")+
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        panel.background = element_blank(),
        axis.title = element_blank())+
  geom_text(aes(x=2,y=1.6),vjust=1,
            label=TeX(r"(\overset${WT}$)"))+
  geom_text(aes(x=5,y=1.6),vjust=0.7,
            label=TeX(r"(\overset${\textit{Nnat}{^+}{^/}{^-}{^p}}{\textit{Light}}$)"))+
  geom_text(aes(x=8,y=1.6),vjust=0.7,
            label=TeX(r"(\overset${\textit{Nnat}{^+}{^/}{^-}{^p}}{\textit{Heavy}}$)"))

本来是想用拼图的方式把3个图组合到一起的，但是遇到了点问题，暂时不知道是什么原因了，可以把第三个图的代码叠加到第一个图上

p1+
  geom_rect(aes(xmin=0.5,xmax=3.5,ymin=305.5,ymax=310.5),
            fill="#c1bebf",
            inherit.aes = FALSE)+
  geom_rect(aes(xmin=3.5,xmax=6.5,ymin=305.5,ymax=310.5),
            fill="#71cddd",
            inherit.aes = FALSE)+
  geom_rect(aes(xmin=6.5,xmax=9.5,ymin=305.5,ymax=310.5),
            fill="#3a54a3",
            inherit.aes = FALSE)+
  theme(plot.margin = unit(c(1.5,1,1,1),'lines'))+
  annotate(geom="text",x=2,y=335.6,vjust=1,
            label=TeX(r"(\overset${WT}$)"))+
  annotate(geom="text",x=5,y=335.6,vjust=0.7,
            label=TeX(r"(\overset${\textit{Nnat}{^+}{^/}{^-}{^p}}{\textit{Light}}$)"))+
  annotate(geom="text",x=8,y=335.6,vjust=0.7,
            label=TeX(r"(\overset${\textit{Nnat}{^+}{^/}{^-}{^p}}{\textit{Heavy}}$)"))+
  coord_cartesian(clip = "off") -> p1.1
p1.1

最后再拼图

library(patchwork)
p1.1/p2+
  plot_layout(heights = c(4,2))

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！