Loading...
墨滴

weiiioyo

2021/12/10  阅读:46  主题:默认主题

单细胞文章绘图——热图

最近在读一篇文章1 ,看到了以下热图(图1),这张热图解决了基因数目过多会基因名称重叠以及细胞类型重叠的问题。于是,我“牺牲”了夜生活对此图进行了剖析和复刻。

图1
图1

剖析如下:

  1. 右半边热图,可使用ggplot2中的geom_tile()完成,然后黑色虚线可通过geom_hline()geom_vline()添加。横向虚线对基因进行了分组,纵向虚线对细胞类型进行了分组。
  2. 细胞类型的标签分上下两部分,代表这个热图有双x轴。双x轴的添加,要求x轴必须是连续型函数,不能直接是类别因子。通用scale_x_continous()dup_aixs添加。
  3. 左半部分是基因的注释。基因的背景色块可通过geom_tile()添加,然后在色块上用geom_text()添加基因名称。
  4. 每个色块的基因分成三列,那么每个基因的横坐标就是他们所在的列号(1,2,3),基因的纵坐标就是他们的行号(1,2,3,4,5,...)。
  5. 基因色块中间有空隙,可以通过geom_hline添加白色伪装成空隙。
  6. 最后将左右部分的图拼接起来,可以通过patchwork这个包实现。

首先,我们来绘制右半部分的热图。在此步骤中,我遇到了如何添加双x轴的难题(因为菜),因为scale_x_discrete中是没有sec.axis参数的,通过 "ggplot2 dual discrete axis"关键词检索找到了解决方法。找到答案后,开始整活!

ps: 我这里没有用图片文章的数据,只是用的手头上现有的数据,downsample到5000个细胞做的图。

library(Seurat)
library(ggplot2)
library(tidyverse)
library(patchwork)
# 我这里导入的是 GSE131907中的5000个细胞。
obj=readRDS('obj_5000.rds')
# 多虑掉细胞类型不明确的细胞
obj=subset(obj,subset=Minor.CellType!="Undetermined")
Idents(obj)=obj$Minor.CellType
# 这里,我已经给每个基因编好了横纵坐标
gene=read.table('gene.txt',header=T,sep="\t",quote='',stringsAsFactors = F)
gene
# group为基因的分组信息,用来画热图中的横向黑色虚线以及基因注释色块中的空隙
#        gene group x  y
#1       IL7R    g1 1 37
#2       TRAC    g1 1 36
#3       CD3E    g1 1 35
#4     GIMAP7    g1 1 34
#5        LTB    g1 1 33
#6       CCR7    g1 1 32
#7      TRBC2    g1 1 31
#8       CD3D    g1 1 30
# 利用Seurat中AverageExpression函数为每个Minor.CellType中的细胞类型求取相应基因表达量的平均值
vag_exp=AverageExpression(obj,assays = "RNA",features = gene$gene,group.by = "Minor.CellType",slot="data")
# 提取Minor.CellType以及他们对应的分组,他们的分组信息存储在Major.CellType里
celltype=unique(FetchData(obj,vars = c("Major.CellType","Minor.CellType")))
# Major.CellType用来画热图中的纵向虚线
celltype$Major.CellType=factor(celltype$Major.CellType,levels=c("T/NK cells","B lymphocytes","Myeloid cells","MAST cells""Endothelial cells","Epithelial cells","Fibroblasts"))
celltype=celltype[order(celltype$Major.CellType),]
celltype$Minor.CellType=factor(celltype$Minor.CellType,levels = unique(celltype$Minor.CellType))
gene=gene[!duplicated(gene$gene),]
gene$gene=factor(gene$gene,levels = gene$gene)

# 求z score,以及通过因子设置固定热图的基因顺序和细胞类型的标签的顺序
dat=base::apply(vag_exp$RNA,1,function(x) (x-mean(x))/sd(x))%>%t()%>%as.data.frame()%>%rownames_to_column('Gene')%>%reshape2::melt()
dat$Gene=factor(dat$Gene,levels = rev(unique(dat$Gene)))
dat$variable=factor(dat$variable,levels = levels(celltype$Minor.CellType))
dat=dat[order(dat$Gene),]
# 热图颜色
heatmap_color=RColorBrewer::brewer.pal(name = "RdBu",n=11)
pal=rev(colorRampPalette(heatmap_color)(100))

# 双x轴的标签
label1=levels(celltype$Minor.CellType)[seq(1,length(levels(celltype$Minor.CellType)),by=2)]
label2=levels(celltype$Minor.CellType)[seq(2,length(levels(celltype$Minor.CellType)),by=2)]

dat[dat$value>3,'value']=3
dat[dat$value < -2.5,'value']=-2.5
p1=ggplot(dat,aes(as.numeric(variable),Gene,fill=value))+
geom_tile()+
scale_y_discrete(expand = c(0,0))+
scale_x_continuous(expand = c(0,0),
                       breaks = seq(1,length(levels(celltype$Minor.CellType)),by=2),
                       labels = label1,
                      sec.axis =dup_axis(breaks = seq(2,length(levels(celltype$Minor.CellType)),by=2),
                                        labels = label2))+
    scale_fill_gradientn(colors=pal,limits=c(-2.5,3),name="Z Score")+
    # 累加每个group中的成员数量,每次的累加值就是每条虚线所在的位置
    geom_hline(yintercept=as.numeric(cumsum(rev(table(gene$group)[-1]))+.5),linetype=2)+
    geom_vline(xintercept=as.numeric(cumsum(table(celltype$Major.CellType))+.5),linetype=2)+
    theme(text = element_text(face="bold"),
          axis.title = element_blank(),
          axis.ticks=element_blank(),
         axis.text.y=element_blank(),
         axis.text.x.top=element_text(angle=90,hjust=0,vjust=.5),
         axis.text.x.bottom=element_text(angle=90,hjust=1,vjust=.5))

好了,右半部分的热图我们已经完成,如图2所示。

图2
图2

接下来,剩下左半部分的基因注释了。继续整!

p2=ggplot(gene,aes(x,y,fill=group))+
geom_tile()+
geom_text(aes(label=gene),family="serif",fontface="italic")+
scale_y_continuous(expand = c(0,0))+
scale_fill_brewer(palette = "Pastel2")+
theme(text=element_text(),
      axis.text = element_blank(),
      axis.title =element_blank(),
      axis.ticks = element_blank(),
      legend.position="none")+
scale_x_continuous(expand = c(0,0))+geom_hline(yintercept = as.numeric(cumsum(rev(table(gene$group)[-1]/3)))+.5,color="white")

左半部分如图3所示:

图3
图3

左右两部分的图都画完成后,我们就可以拼图啦。

pic.heatmap=p2+p1+plot_layout(ncol = 2, widths  = c(13))
pic.heatmap

拼图效果如图4: 图4

大家发现没有,两个图之间的空隙有些大。通过关键词"patchwork change distance between plots"检索可获得答案

pic.heatmap & theme(plot.margin = margin(0,0,0,0))

最终的效果图(图5)如下:

图5
图5

总结:

  1. 复刻过程中,学会用关键字检索解决问题
  2. 如果换了数据集和基因,基因坐标得重新计算。需要添加一个自动计算基因位置的函数功能。

此文使用的数据集和基因集可在github中获得:https://github.com/weiiioyo/singlecell-heatmap

参考文献

  1. Single-cell analyses reveal key immune cell subsets associated with response to PD-L1 blockade in triple-negative breast cancer

weiiioyo

2021/12/10  阅读:46  主题:默认主题

作者介绍

weiiioyo