单细胞转录组基因表达差异的显著性分析

之前算基因表达量已知都用的averageexpression,但是这样就没有办法标注出两个组之间的差异是否具有显著性。如果把一个细胞当做普通转录组的一个样本，那么按理来说不同组之间有N个样本，但是之前一直也没有考虑过这样的问题。

#因为我想合并不同的sample, 所以需要更改orig.ident
BCs$orig.ident<-factor(BCs$orig.ident,
                       levels = c('A1','A2','B1','B2'),
                       labels = c('A','A','B','B'))

1. 提取自定义基因集的表打矩阵

ls1<- GetAssayData(object = sample1, slot = "counts")
ls2 <- GetAssayData(object = sample2, slot = "counts")
gene <- c("")
gene <- as.vector(gene)
ls1_plot <- ls1[gene,]
ls2_plot <- ls2[gene,]
ls1_plot<-as.data.frame(ls1_plot)
ls2_plot<-as.data.frame(ls2_plot)
rm(ls1,ls2,sample1,sample2)

2. 其实实质上是获得两个文件
①表达量数据（列为样本，行为基因）
②注释信息（第一列是样本，第二列是组别）

#生成注释信息
sample1 <-colnames(ls1_plot)
sample1 <-as.data.frame(sample1)
sample1$Type <- c(rep("",times =))
sample2<-colnames(ls2_plot)
sample2 <-as.data.frame(sample2)
sample2$Type <- c(rep("",times = ))
colnames(sample1) <-c("Sample1","Type")
colnames(sample2) <-c("Sample2","Type")
annotate <- rbind(sample1,sample2)
##合并表达量数据，列为样本，行为基因
expression <- cbind(ls1_plot,ls2_plot)
rm(sample1,sample2,ls1_plot,ls2_plot)

3.接下来就是通过ggplot2进行绘图了,这里列出的是箱线图，其他的类似

library(RColorBrewer)
library(ggpubr)
library(ggplot2)
library(cowplot)

expression <- t(expression)#这里需要给表达矩阵转置，因为要符合ggplot作图方式
Exp_plot <- as.data.frame(expression)
Exp_plot$sam=annotate$Type

plist2<-list()#创建一个空列表，用来存储循环的产出
for (i in 1:length(gene)){
  bar_tmp<-Exp_plot[,c(gene[i],"sam")]#循环提取每个基因表达信息
  colnames(bar_tmp)<-c("Expression","sam")#统一命名
  my_comparisons1 <- list(c("Asymptomatic", "Mild")) #设置比较组
  #...(这里可以添加很多comparsions的分类)
  p1<-ggboxplot(bar_tmp,#ggboxplot画箱线图
                 x="sam",#x轴为组别
                 y="Expression",#y轴为表达量
                 color="sam",#用样本分组填充
                 fill=NULL,
                 add = "jitter",#添加散点
                 bxp.errorbar.width = 0.6,
                 width = 0.4,
                 size=0.01,
                 font.label = list(size=30), 
                 palette = col)+theme(panel.background =element_blank())
 p1<-p1+theme(axis.line=element_line(colour="black"))+theme(axis.title.x = element_blank())#坐标轴修饰
 p1<-p1+theme(axis.title.y = element_blank())+theme(axis.text.x = element_text(size = 15,angle = 45,vjust = 1,hjust = 1))#横坐标文字设置
 p1<-p1+theme(axis.text.y = element_text(size = 15))+ggtitle(gene[i])+theme(plot.title = element_text(hjust = 0.5,size=15,face="bold"))#标题设置
 p1<-p1+theme(legend.position = "NA")#（因为有组图，横坐标分组了，所以不需要设置legend）
 p1<-p1+stat_compare_means(method="t.test",hide.ns = F,
                              comparisons =c(my_comparisons1,my_comparisons2,my_comparisons3,my_comparisons4,my_comparisons5,my_comparisons6),
                              label="p.signif")#显著性检验用t检验，添加不同比较组。详情可以查看stat_compare_means函数帮助信息
plist2[[i]]<-pb1 #将画好的图储存于plist2列表，并不断赋值循环直到结束
}

参考链接：https://www.jianshu.com/p/1bfc545d7cfa