ixxmu
commented
2 years ago CNS图表复现：Boxplot比较单细胞亚群的组织分布 by 生信随笔
对单细胞数据进行亚群注释之后，我们往往想比较某亚群，例如CD8Tex，是倾向于分布在实验组还是对照组，例如癌组织，癌旁组织，转移癌组织，淋巴组织？这时候有很多策略去做这种多组间的比较。如我在《CNS图表复现：OR指数比较单细胞亚群的组织偏好》所述：
第一种，是《Therapy-Induced Evolution of Human Lung Cancer Revealed by Single-Cell RNA Sequencing》这篇Cell的做法，这篇文章有三种处理组（TN，RD和PD），或许是因为用的Smart-seq2测序，每个样本得到的细胞数量确实不多，因此作者简单粗暴的把同一组内的样本细胞加和计算亚群细胞的频率，进行统计，绘图如下：
这种做法其实是存在一定问题的，它画不出所有样本分布的散点图，因为本质上作者把同一组的亚群看成一个样本。我在单细胞思考：Cell作者一定是对的吗？此文中复现了作者的处理思路和图表。
我也在如上推文中给出了第二种多组间亚群比较的策略，即按照每个样本的细胞总和进行百分比的校正，然后比较频率。这也是绝大多数文章采取的策略，例如《Pan-cancer single-cell landscape of tumor-infiltrating T cells》张泽民团队比较多癌种某单细胞亚群的频率分布差异，一个散点代表一个样本：
第三种策略也是出现在《Pan-cancer single-cell landscape of tumor-infiltrating T cells》张泽民老师的这篇文章中，即Fig1F的这幅图，利用OR比值比的统计学方法，比较血液，正常组织和肿瘤组织，各单细胞亚群的分布差异，详见《CNS图表复现：OR指数比较单细胞亚群的组织偏好》。
第四种策略也是张泽民团队经常使用的一个统计方法，Ro/e，这个指标是观察到的细胞数与期望细胞数的比值，用于量化每个亚群对组织的偏好程度
有一位锲而不舍的读者想让我复现张泽民老师团队的Boxplot画法，其实张泽民团队公开的代码就有这个图的画法，我在此基础上做了一些调整，见下文。
二. 图表复现
library(data.table)
library(readr)
library(dplyr)
library(plyr)
library(ggpubr)
library(ggplot2)
library(patchwork)
# 本次画图的技巧支持多线程运行
RhpcBLASctl::omp_set_num_threads(2)
doParallel::registerDoParallel(cores = 2)
#### 1.加载数据
meta.tb = read_rds("./test_data.rds")
#### 2.数据预处理
### 2.1 （可选）这里选择CD8T亚群进行可视化，这里的loc的L分组作者原文是不进行画图展示的，因此过滤掉
dat.plot = data.table(filter(meta.tb,stype == "CD8"&loc != c("L")))
### 2.2 绘图需要指定celltype标签(celltype)、分组信息(Group)、以及患者编号(SampleID)
dat.plot$Celltype = dat.plot$meta.cluster
dat.plot$Group = factor(dat.plot$loc,levels = c("P","N","T"))
dat.plot$SampleID = dat.plot$patient.uid
### 2.3 统计细胞出现的频数和频率
dat.plot <- dat.plot[,.(N=.N),by=c("cancerType","SampleID","Celltype","Group")]
sum.data <- dat.plot[,{.(NTotal = sum(.SD$N))},by=c("cancerType","SampleID","Group")]
dat.plot = left_join(dat.plot,sum.data)
dat.plot$freq = dat.plot$N / dat.plot$NTotal *100
head(dat.plot)
#### 3.可视化
### 3.1 指定输出的目录
out.prefix.loc.cmp = "./Outplot/"
### 3.2 指定比较组以生成p值
compraisions.list=list(c("P","N"),c("N","T"),c("P","T"))
### 3.3.设置配色
g.colSet = list(group = c("P" = "#e75945","N" = "#f39b7f","T" = "#00a087","L" = "#BC3C29FF"))
### 3.4 ggboxplot可视化（换成ggplot2也行），同时lapply支持多线程运行
p.list.perMcls <- lapply(unique(sort(dat.plot$Celltype)),function(mcls){
text.size = 10
text.angle = 45
text.hjust = 1
legend.position = "none"
p <- ggboxplot(dat.plot[Celltype== mcls ,],x="Group",y="freq",
color = "Group", legend="none",title=mcls,
#fill = "loc",alpha=0.8,
#add = "none",outlier.shape=NA) +
xlab="",ylab="frequency",
add = "jitter",outlier.shape=NA) +
scale_color_manual(values=g.colSet$group) +
#scale_color_manual(values=c("P"="#E41A1C","N"="#377EB8","T"="#4DAF4A")) +
#scale_fill_manual(values=c("P"="#E41A1C","N"="#377EB8","T"="#4DAF4A")) +
stat_compare_means(label="p.format",comparisons=compraisions.list) +
coord_cartesian(clip="off") +
theme(plot.title = element_text(size = text.size+2,color="black",hjust = 0.5),
axis.ticks = element_line(color = "black"),
axis.title = element_text(size = text.size,color ="black"),
axis.text = element_text(size=text.size,color = "black"),
axis.text.x = element_text(angle = text.angle, hjust = text.hjust ), #,vjust = 0.5
#axis.line = element_line(color = "black"),
#axis.ticks = element_line(color = "black"),
#panel.grid.minor.y = element_blank(),
#panel.grid.minor.x = element_blank(),
# panel.grid=element_blank(), # 去网格线
legend.position = legend.position,
legend.text = element_text(size= text.size),
legend.title= element_text(size= text.size),
# panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),
strip.background = element_rect(color="black",size= 1, linetype="solid") # fill="#FC4E07",
)
# ggsave打开即可保存图片
#ggsave(sprintf("%s.compostion.PNT.%s.png",out.prefix.loc.cmp,mcls),width=2,height=4.5)
# ggsave(sprintf("%sepi.compostion.PNT.%s.pdf",out.prefix.loc.cmp,make.names(mcls)),
# width=5,height= 7.5 ,useDingbats=F,units = "cm")
return(p)
})
p.prop = wrap_plots(p.list.perMcls,ncol = 6)
p.prop
ggsave(p.prop,filename = paste0(out.prefix.loc.cmp,"Propotion_boxplot.pdf"),width = 22,height = 12,units = "cm")
测试数据及全文markdown代码后台回复：
https://mp.weixin.qq.com/s/eWYGQGbW4XXLIakUkA7fag