论文
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、生物信息学入门学习资料及自己的学习笔记!