两组单细胞样品的不同亚群比例差异的火山图展现
到目前为止，单细胞转录组费用仍然是居高不下，所以绝大部分情况下大家做两个分组，每个组内也就是三五个样品而已。这样的话两个分组之间的不同单细胞亚群的比例差异其实往往是需要最后使用流式细胞等价格相对低廉的实验技术去扩大样品队列去验证一下。
而不同单细胞样品的不同亚群比例差异，前面我们介绍过：展示细胞比例变化之balloonplot和马赛克图，以及 展示细胞比例变化之桑基图，但它们通常并没有分组比较。最近看到了2020发表在NC杂志的文章：《Integrated single cell analysis of blood and cerebrospinal fluid leukocytes in multiple sclerosis.》，其数据集是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138266 ，可以看到是2个分组， 每个分组是11个样品，可以两个分组进行单细胞比例比较了。
首先，仍然是经典的降维聚类分群和标记基因对亚群进行命名，如下所示：
这些基因大家基本上都是可以背诵下来了，然后，可以根据样品的分组拆开看单细胞亚群比例差异：
如果肉眼看，基本上也可以判断出来NK1这个细胞亚群在CSF分组里面基本上没有了，而Mono2相反，本来是在blood里面基本上没有，但是在CSF这个疾病分组里面比例还蛮高的。但是肉眼看不清楚其它并不很明显的细胞亚群，所以有了右边的火山图展现两个分组的单细胞亚群比例变化。
下面我们来演示一下这样的火山图如何绘制，其实最重要的反而是数据如何获得！我们这里只能说选择模拟数据，如下所示的代码：
library(dplyr)
library(ggplot2)
library(dplyr)
set.seed(1)
n=260000
phe = data.frame(
celltype= sample(LETTERS ,n,replace = T) ,
orig.ident = sample(paste0(c('case','control'),c(1:10,1:10)),n,replace = T)
)
head(phe)
table(phe$celltype)
table(phe$orig.ident)
phe$group = ifelse(grepl('case',phe$orig.ident) , 'case','control')
table(phe$group)
head(phe)
如下所示：
> table(phe$group)
case control
130094 129906
> dim(phe)
[1] 260000 3
> head(phe)
celltype orig.ident group
1 Y case3 case
2 D case1 case
3 G control8 control
4 A control8 control
5 B case1 case
6 W case9 case
我们的两个分组，case和control都是10个单细胞样品，我们的单细胞亚群是26个，总共是26万细胞。这个模拟的结果，就是大家对单细胞数据集进行降维聚类分群后的，参考前面的例子：人人都能学会的单细胞聚类分群注释 ，自己拿到这样的 phe变量。
接下来就可以对这个phe进行系列统计啦，代码如下所示是：
pro='group'
tb <- data.frame(table(phe$celltype,phe$orig.ident,
phe[,pro]))
tb=tb[which(tb$Freq != 0),]
tb=tb[,c(1,3,4)]
head(tb)
tb$Total <- apply(tb,1,function(x)sum(tb[tb$Var1 == x[1],3]))
tb<- tb %>% mutate(Percentage = round(Freq/Total,3) * 100)
tb=tb[,c(1,2,5)]
tb$Var1=as.factor(tb$Var1)
tb$Var3=as.factor(tb$Var3)
head(tb)
df= do.call(rbind,
lapply(split(tb,tb$Var1), function(x){
# x= split(tb,tb$Var1)[[1]]
tmp = t.test(x$Percentage ~ x$Var3)
return(c(tmp$p.value, tmp$estimate[1]-tmp$estimate[2]))
}))
head(df)
可以看到，这26个单细胞亚群，在2个分组的比例差异都不大，而且基本上不太可能有统计学显著，因为我是随机模拟的数据，并不是真正的单细胞数据分析实战。
> df
mean in group case
A 0.58496895 0.10
B 0.22798735 -0.18
C 0.52050726 0.14
D 0.09569494 0.20
E 0.61076971 0.08
F 0.66543654 -0.12
G 0.67605420 0.06
H 0.32041576 0.28
I 0.49641609 -0.10
J 0.57703994 0.14
K 0.91762160 0.02
L 0.12634692 -0.28
M 0.39033329 0.18
N 0.70214693 0.08
O 0.30149845 -0.14
P 0.86727855 -0.02
Q 1.00000000 0.00
R 0.75319312 -0.08
S 0.45515512 0.14
T 0.45181182 -0.16
U 0.87092289 -0.04
V 0.22497991 -0.24
W 0.29969641 0.16
X 0.41407380 0.12
Y 0.54543081 -0.08
Z 0.36185885 0.18
下面就需要对这个统计结果进行火山图绘制啦，如下所示的基础函数：
pdf("p.pdf",width = 5.5,height = 5.5)
plot(df[,2],-log10(df[,1]),pch = 19)
abline(v = 0,lty = 2)
text(df[,2],-log10(df[,1]),rownames(df),cex = 0.6,pos = 2,col = "red")
dev.off()
以及高级绘图函数：
library(ggrepel)
library(ggplot2)
colnames(df) = c("pval","Difference")
df = as.data.frame(df)
df$threshold = factor(ifelse(df$Difference > 0 ,'Down','UP'))
ggplot(df,aes(x=Difference,y=-log10(pval),color=threshold))+
geom_point()+
geom_text_repel(
aes(label = rownames(df)),
size = 4,
segment.color = "black", show.legend = FALSE ) + #添加细胞群的名字
theme_bw()+#修改图片背景
theme(legend.title = element_blank()) + #不显示图例标题
ylab('-log10(pval)')+ #修改y轴名称
xlab('Difference')+ #修改x轴名称
geom_vline(xintercept=c(0),lty=3,col="black",lwd=0.5)
ggsave("p1.pdf",width = 5,height = 4)
效果如下所示：
学徒作业
就是根据我们的单细胞流程，对这个数据集是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138266 ，进行降维聚类分群，计算不同单细胞亚群比例，在两个分组绘制这样的火山图。
写在文末
我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。