单细胞测序—标准流程代码（3）—marker 基因富集分析_差异基因
过了很久之后才想起来继续整理单细胞测序的标准分析流程。书接上回单细胞测序—标准流程代码(2) — 标记基因与细胞注释，这篇帖子主要关注的是富集分析。
1 marker基因富集分析
主要是step2-anno-go-kegg-reactome.R脚本根据物种中调用com_go_kegg_ReactomePA_human.R或者com_go_kegg_ReactomePA_mice.R脚本。
1.1 step2-anno-go-kegg-reactome.R
代码语言：r复制rm(list=ls())
options(stringsAsFactors = F)
source('scRNA_scripts/lib.R')
getwd()
d='figures-marker_cosg_com_go_kegg/'
dir.create(d)
setwd(d)
load('../check-by-celltype/qc-_marker_cosg.Rdata')
head(marker_cosg)
symbols_list <- as.list(as.data.frame(apply(marker_cosg$names,2,head,100)))
symbols_list
source('../com_go_kegg_ReactomePA_human.R')
#source('../com_go_kegg_ReactomePA_mice.R')
com_go_kegg_ReactomePA_human(symbols_list, pro='SLE' )
setwd('../')代码解释
首先load之前得到的check-by-celltype文件夹的qc-_marker_cosg.Rdata
注：得到qc-_marker_cosg.Rdata的部分代码
代码语言：r复制#存在check-all-markers.R脚本中
marker_cosg <- cosg(
sce.all.int,
groups='all',
assay='RNA',
slot='data',
mu=1,
n_genes_user=100)
save(marker_cosg,file = paste0(pro,'_marker_cosg.Rdata'))
head(marker_cosg)我们先观察下marker_cosg，这是一个list，包含sce对象中每一个细胞分群的前100个的marker基因
接下来
代码语言：r复制symbols_list <- as.list(as.data.frame(apply(marker_cosg$names,2,head,100)))这个代码段将 marker_cosg 对象中的基因符号提取并转换为列表。具体步骤如下：
apply(marker_cosg$names, 2, head, 100)应用于marker_cosg$names的每一列，提取前 100 个基因符号。2表示沿着列操作。
as.data.frame()将apply的结果转换为数据框。as.list()将数据框转换为列表，每个列表元素对应一个细胞类型或分群的前 100 个基因符号。
最终，symbols_list 是一个列表，每个元素包含某个细胞类型或分群的前 100 个基因符号。
同样，观察下symbols_list列表，含有每个分群中的前100的mark基因
接着加载source('../com_go_kegg_ReactomePA_human.R')，调用函数
1.2 com_go_kegg_ReactomePA_human.R
这段代码会出5张关于富集分析的的图，分别是KEGG通路富集、Reactome通路富集、以及GO中BP（生物过程）、CC（细胞组分）、MF（分子功能）的富集分析的气泡图，记录了单细胞各个分群之间marker基因富集结果。
代码语言：r复制com_go_kegg_ReactomePA_human <- function(symbols_list ,pro){
library(clusterProfiler)
library(org.Hs.eg.db)
library(ReactomePA)
library(ggplot2)
library(stringr)
# 首先全部的symbol 需要转为 entrezID
gcSample = lapply(symbols_list, function(y){
y=as.character(na.omit(AnnotationDbi::select(org.Hs.eg.db,
keys = y,
columns = 'ENTREZID',
keytype = 'SYMBOL')[,2])
)
y
})
gcSample
# 第1个注释是 KEGG
xx <- compareCluster(gcSample, fun="enrichKEGG",
organism="hsa", pvalueCutoff=0.05)
dotplot(xx) theme(axis.text.x=element_text(angle=45,hjust = 1))
scale_y_discrete(labels=function(x) str_wrap(x, width=50))
ggsave(paste0(pro,'_kegg.pdf'),width = 10,height = 14)
# 第2个注释是 ReactomePA
xx <- compareCluster(gcSample, fun="enrichPathway",
organism = "human",
pvalueCutoff=0.05)
dotplot(xx) theme(axis.text.x=element_text(angle=45,hjust = 1))
scale_y_discrete(labels=function(x) str_wrap(x, width=50))
ggsave(paste0(pro,'_ReactomePA.pdf'),width = 10,height = 14)
# 然后是GO数据库的BP,CC,MF的独立注释
# Run full GO enrichment test for BP
formula_res <- compareCluster(
gcSample,
fun="enrichGO",
OrgDb="org.Hs.eg.db",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# Run GO enrichment test and merge terms
# that are close to each other to remove result redundancy
lineage1_ego <- simplify(
formula_res,
cutoff=0.5,
by="p.adjust",
select_fun=min
)
pdf(paste0(pro,'_GO_BP_cluster_simplified.pdf') ,width = 15,height = 14)
print(dotplot(lineage1_ego, showCategory=5) theme(axis.text.x=element_text(angle=45,hjust = 1))
scale_y_discrete(labels=function(x) str_wrap(x, width=50)) )
dev.off()
write.csv(lineage1_ego@compareClusterResult,
file=paste0(pro,'_GO_BP_cluster_simplified.csv'))
# Run full GO enrichment test for CC
formula_res <- compareCluster(
gcSample,
fun="enrichGO",
OrgDb="org.Hs.eg.db",
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# Run GO enrichment test and merge terms
# that are close to each other to remove result redundancy
lineage1_ego <- simplify(
formula_res,
cutoff=0.5,
by="p.adjust",
select_fun=min
)
pdf(paste0(pro,'_GO_CC_cluster_simplified.pdf') ,width = 15,height = 14)
print(dotplot(lineage1_ego, showCategory=5) theme(axis.text.x=element_text(angle=45,hjust = 1))
scale_y_discrete(labels=function(x) str_wrap(x, width=50)) )
dev.off()
write.csv(lineage1_ego@compareClusterResult,
file=paste0(pro,'_GO_CC_cluster_simplified.csv'))
# Run full GO enrichment test for MF
formula_res <- compareCluster(
gcSample,
fun="enrichGO",
OrgDb="org.Hs.eg.db",
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# Run GO enrichment test and merge terms
# that are close to each other to remove result redundancy
lineage1_ego <- simplify(
formula_res,
cutoff=0.5,
by="p.adjust",
select_fun=min
)
pdf(paste0(pro,'_GO_MF_cluster_simplified.pdf') ,width = 15,height = 14)
print(dotplot(lineage1_ego, showCategory=5) theme(axis.text.x=element_text(angle=45,hjust = 1))
scale_y_discrete(labels=function(x) str_wrap(x, width=50)) )
dev.off()
write.csv(lineage1_ego@compareClusterResult,
file=paste0(pro,'_GO_MF_cluster_simplified.csv'))
}以这两张图为例
SLE_kegg.pdf
SLE_ReactomePA.pdf
1.3 com_go_kegg_ReactomePA_mice.R
针对mouse物种，与1.2，二者只需要调用一个。
代码语言：r复制com_go_kegg_ReactomePA_human <- function(symbols_list ,pro){
library(clusterProfiler)
library(org.Mm.eg.db)
library(ReactomePA)
library(ggplot2)
library(stringr)
# 首先全部的symbol 需要转为 entrezID
gcSample = lapply(symbols_list, function(y){
y=as.character(na.omit(select(org.Mm.eg.db,
keys = y,
columns = 'ENTREZID',
keytype = 'SYMBOL')[,2])
)
y
})
gcSample
# 第1个注释是 KEGG
xx <- compareCluster(gcSample, fun="enrichKEGG",
organism="mmu", pvalueCutoff=0.05)
dotplot(xx) theme(axis.text.x=element_text(angle=45,hjust = 1))
scale_y_discrete(labels=function(x) str_wrap(x, width=50))
ggsave(paste0(pro,'_kegg.pdf'),width = 10,height = 8)
# 第2个注释是 ReactomePA
xx <- compareCluster(gcSample, fun="enrichPathway",
organism = "mouse",
pvalueCutoff=0.05)
dotplot(xx) theme(axis.text.x=element_text(angle=45,hjust = 1))
scale_y_discrete(labels=function(x) str_wrap(x, width=50))
ggsave(paste0(pro,'_ReactomePA.pdf'),width = 10,height = 8)
# 然后是GO数据库的BP,CC,MF的独立注释
# Run full GO enrichment test for BP
formula_res <- compareCluster(
gcSample,
fun="enrichGO",
OrgDb="org.Mm.eg.db",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# Run GO enrichment test and merge terms
# that are close to each other to remove result redundancy
lineage1_ego <- simplify(
formula_res,
cutoff=0.5,
by="p.adjust",
select_fun=min
)
pdf(paste0(pro,'_GO_BP_cluster_simplified.pdf') ,width = 15,height = 8)
print(dotplot(lineage1_ego, showCategory=5) theme(axis.text.x=element_text(angle=45,hjust = 1))
scale_y_discrete(labels=function(x) str_wrap(x, width=50)) )
dev.off()
write.csv(lineage1_ego@compareClusterResult,
file=paste0(pro,'_GO_BP_cluster_simplified.csv'))
# Run full GO enrichment test for CC
formula_res <- compareCluster(
gcSample,
fun="enrichGO",
OrgDb="org.Mm.eg.db",
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# Run GO enrichment test and merge terms
# that are close to each other to remove result redundancy
lineage1_ego <- simplify(
formula_res,
cutoff=0.5,
by="p.adjust",
select_fun=min
)
pdf(paste0(pro,'_GO_CC_cluster_simplified.pdf') ,width = 15,height = 8)
print(dotplot(lineage1_ego, showCategory=5) theme(axis.text.x=element_text(angle=45,hjust = 1))
scale_y_discrete(labels=function(x) str_wrap(x, width=50)) )
dev.off()
write.csv(lineage1_ego@compareClusterResult,
file=paste0(pro,'_GO_CC_cluster_simplified.csv'))
# Run full GO enrichment test for MF
formula_res <- compareCluster(
gcSample,
fun="enrichGO",
OrgDb="org.Mm.eg.db",
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05
)
# Run GO enrichment test and merge terms
# that are close to each other to remove result redundancy
lineage1_ego <- simplify(
formula_res,
cutoff=0.5,
by="p.adjust",
select_fun=min
)
pdf(paste0(pro,'_GO_MF_cluster_simplified.pdf') ,width = 15,height = 8)
print(dotplot(lineage1_ego, showCategory=5) theme(axis.text.x=element_text(angle=45,hjust = 1))
scale_y_discrete(labels=function(x) str_wrap(x, width=50)) )
dev.off()
write.csv(lineage1_ego@compareClusterResult,
file=paste0(pro,'_GO_MF_cluster_simplified.csv'))
}1.4 一些疑问？
问1：得到的marker_cosg，是一个包含两个元素的列表，一个是names，一个是scores。这个scores存储了哪些信息？
答：在 marker_cosg 列表中，scores 存储了每个基因相对于不同细胞群（groups）的得分信息。这些得分通常用于评估每个基因在特定细胞群中的表达显著性或区分能力。以下是 scores 的具体内容和作用：
scores 的内容
- 得分矩阵：scores通常是一个矩阵或数据框，其中每行对应一个基因，每列对应一个细胞群或组别。
- 得分值：矩阵中的值代表每个基因在不同细胞群中的显著性得分。较高的得分通常表示该基因在该特定细胞群中具有更显著的表达模式或是更具代表性的 marker 基因。
scores 的作用
- 筛选 Marker 基因：通过分析
scores，可以确定哪些基因在特定的细胞群中表现出显著的差异表达，从而筛选出潜在的 marker 基因。 - 功能注释：高得分的基因通常用于进一步的功能注释分析，如 GO、KEGG 或 Reactome 通路富集分析，以探索它们在生物学过程中的角色。
问2:cosg与findallMarker的异同？
答：参考链接：神兵利器——单细胞细胞类群基因marker鉴定新方法：COSG
总体来说COSG执行效果更快，更科学，这里我有一个疑问就是我们采用Seurat 官方的 FindMarkers与FindAllMarkers函数时往往加一个参数only.pos = TRUE，意为仅返回在给定簇中表达上调（正向表达）的基因，而不包括在其他簇中下调的基因。这通常用于识别某个簇中特异性表达的基因。那么执行cosg后返回的marker基因是否也都是上调的呢？
我观察到返回的基因score都是正值，且暂未发现类似于FindAllMarkers中的only.pos参数。因此我猜想cosg默认返回的是上调的marker基因。那么由cosg得到marker基因富集出来的各种通路就都是上调的，而不存在下调。
ps:猜想，不一定正确。
问3:Reactome 通路与kegg的区别？
答：Reactome和KEGG（Kyoto Encyclopedia of Genes and Genomes）是两种常用的生物信息学数据库，它们都提供了关于生物通路的信息，但在内容、数据组织和用途上存在一些区别。
- Reactome:
是一个开放的、基于人类基因的生物通路数据库，但它也包含了其他物种的通路数据，通常通过人类通路的同源映射生成。Reactome 涵盖了广泛的生物过程，包括信号传导、代谢、基因表达、细胞周期、免疫反应等。它不仅关注代谢通路，还深入到细胞内部的分子反应和调控机制。Reactome 的数据是以事件（event）为基础组织的，涵盖了反应（reaction）和通路（pathway）。通路被分为多个层级的子通路，并且可以在不同的细胞上下文或条件下展示。它采用了类似“事件树”的结构，允许用户逐步深入探索从高层次的生物过程到具体的分子反应。更适用于深入研究分子反应和基因调控的机制，尤其是在非代谢通路方面，如信号传导、基因表达和细胞周期等。Reactome 也常用于转录组学和蛋白质组学数据的功能注释，因为它包含了许多非代谢相关的生物过程。
- KEGG:
KEGG 是一个集成的数据库资源，用于理解生物系统（如细胞、器官、生态系统）及其与环境的关系。
它包括代谢通路图（metabolic pathways），以及与代谢密切相关的其他通路，如遗传信息处理、环境信息处理、细胞过程等。KEGG 的通路更侧重于代谢网络及其相关的基因和化学物质的关系。KEGG 的数据主要组织成一系列的通路图，这些通路图展示了代谢物、酶和基因产物之间的关系。KEGG 通路图具有可视化的特点，图中的节点表示化学物质或基因产物，边表示生化反应或调控关系。KEGG 通常用于研究代谢网络、药物代谢和环境适应等问题。它特别适合代谢组学分析，以及代谢通路富集分析。由于 KEGG 的通路图直观且集中于代谢过程，因此在代谢相关研究中应用广泛。
2 差异基因
2.1 step3_deg_then_anno.R
差异基因的选定与可视化主要在step3_deg_then_anno.R这个脚本中
代码语言：r复制rm(list=ls())
options(stringsAsFactors = F)
source('scRNA_scripts/lib.R')
getwd()
d='figures-deg_then_anno/'
dir.create(d)
setwd(d)
sce = readRDS('../2-harmony/sce.all_int.rds')
load('../phe.Rdata')
sce$celltype=phe$celltype
sce$group=sce$stim
pdf("group_celltype_balloon.pdf")
gplots::balloonplot(table(sce$group,sce$celltype ))
dev.off()
Idents(sce) = sce$group
table(Idents(sce))
degs = lapply(unique(sce$celltype), function(x){
FindMarkers(sce[,sce$celltype==x],ident.1 = 'STIM',
ident.2 = 'CTRL')
})
names(degs)=unique(sce$celltype)
x=degs[[1]]
head(x)
do.call(rbind,lapply(degs, function(x){
table(x$avg_log2FC > 1 )
}))
do.call(rbind,lapply(degs, function(x){
table(x$avg_log2FC < -1 )
}))
# 这个差异分析结果有两种可视化方法
# 每个细胞亚群在两个分组差异后，都有上下调基因，可以火山图
# 然后又可以走 批量注释啦
# com_go_kegg_ReactomePA_human.R
# 这里略
# 火山图需要设置阈值，logFC和p值
degs_list=lapply(names(degs),function(i){
#i=names(degs)[1]
x = degs[[i]]
res = x
colnames(res)
res$symbol=rownames(x)
p=EnhancedVolcano::EnhancedVolcano(res,
lab = res$symbol,
x = 'avg_log2FC',
y = 'p_val_adj',pCutoff = 0.05,
FCcutoff =1)
p
ggsave(file = paste0('deg_for_',i,'_by_volcano.pdf' ),
plot = p,width = 8,height = 6)
res$group = i
res
})
degs_allcluster_type_df=do.call(rbind,degs_list)
degs_allcluster_type_df %>%head()
write.csv(degs_allcluster_type_df,file = "./degs_allcluster_type_df.csv")
saveRDS(degs_allcluster_type_df,file = "./allclusters_degs_sce.markers.Rdata")
degs_allcluster_type_df$cluster=degs_allcluster_type_df$group
length(unique(degs_allcluster_type_df$cluster))
degs_allcluster_type_df$gene = degs_allcluster_type_df$symbol
pak::pkg_install('junjunlab/scRNAtoolVis')
p2=scRNAtoolVis::jjVolcano(diffData = degs_allcluster_type_df,
log2FC.cutoff = 2,
adjustP.cutoff = 0.05,
legend.position = c(0.93, 0.99),
topGeneN=0,#top genes to be labeled in plot, default 5.
# cluster.order=seq(0,23,1),
pSize=0.4,
tile.col = c("#EE3B3B", "#FF7F00", "#CD6600", "#8B2323", "#DEB887", "#76EEC6",
"#F0FFFF", "#008B8B", "#FFB90F", "#F5F5DC", "#1F1F1F", "#66CD00",
"#0000FF", "#97FFFF", "#528B8B", "#9400D3", "#EE1289", "#00BFFF",
"#00FF00", "#191970", "#FFFF00", "#4A708B", "#00FF7F", "#8B8B00",
"#FF1493", "#FFA500", "#8B4513"))
p2
ggsave(plot = p2,filename = "allcelltypes_degs_volcano.pdf",width = 13,height = 10)
print(getwd())
setwd('../')2.2 代码解释
代码语言：r复制sce = readRDS('../2-harmony/sce.all_int.rds')
load('../phe.Rdata')
sce$celltype=phe$celltype
sce$group=sce$stim
pdf("group_celltype_balloon.pdf")
gplots::balloonplot(table(sce$group,sce$celltype ))
dev.off()首先sce = readRDS('../2-harmony/sce.all_int.rds'),注意区别与load Rdata，RDS在read的时候，必须要赋值。
其中sce.all_int.rds存储的是单细胞对象降维聚类分群后的结果（还未细胞注释）
后load phe.Rdata，phe中存储的是细胞注释后meta.data。
代码语言：r复制 #step1-run-basic.R
phe=sce.all.int@meta.data
save(phe,file = 'phe.Rdata')再将phe中的celltype信息添加到sce对象中。
因此我觉得这段代码是有优化的空间的，可以直接load，降维聚类分群注释后的seurat对象（sce.all.int），这样的效果和上述三行的效果一致。
最后sce$stim中存储了分组的信息，在meta.data中新增group一列，存储sce$stim信息，主要是为了方便后续的操作。
首先是对数据的有一个总体的把握，呈现下细胞分群与分组（批次）之间的关系。
group_celltype_balloon.pdf
degs = lapply(unique(sce$celltype), function(x){
FindMarkers(sce[,sce$celltype==x], ident.1 = 'STIM', ident.2 = 'CTRL')
})
names(degs)=unique(sce$celltype)- unique(sce$celltype)：sce$celltype表示单细胞对象
sce中的细胞类型信息。unique()函数会返回所有不同的细胞类型，即去重后的细胞类型列表。 - lapply() ：是R中的一个循环函数，作用是对列表中的每个元素应用同一个函数，并返回一个列表。这里，
lapply()遍历每一种独特的细胞类型，并对每种细胞类型执行指定的函数。 function(x) { ... }
function(x)是一个匿名函数，x代表当前的细胞类型。这个函数的目的是对特定细胞类型下的细胞进行差异表达分析。
sce[,sce$celltype==x]
sce[,sce$celltype==x]：这一部分的代码是用来从sce对象中提取属于当前细胞类型x的所有细胞。sce$celltype == x会生成一个逻辑向量，标记哪些细胞属于当前细胞类型x，而sce, sce$celltype==x则根据这个逻辑向量提取这些细胞的数据。
- FindMarkers()：Seurat包中的一个函数，用于执行差异表达分析。它比较两个细胞群体之间的基因表达差异。ident.1 = 'STIM：代表分组1，即STIM条件下的细胞。ident.2 = CTRL：代表分组2，即CTRL条件下的细胞。
- degs：返回一个列表，每个列表元素对应于一个细胞类型，并包含该细胞类型在
STIM与CTRL条件下的差异表达基因结果。FindMarkers()的结果包括基因名、log2倍数变化值（log2 fold change,avg_log2FC）、p值调整值（p_val_adj）等信息。
do.call(rbind,lapply(degs, function(x){
table(x$avg_log2FC > 1 )
}))
do.call(rbind,lapply(degs, function(x){
table(x$avg_log2FC < -1 )
}))这段代码的作用是汇总和统计在差异表达分析中基因表达上调和下调的情况。
lapply(degs, function(x){ ... })：对degs列表中的每一个元素（即每种细胞类型的差异表达结果）应用一个匿名函数。x$avg_log2FC > 1和x$avg_log2FC < -1：x是degs列表中的某个元素，也就是某个细胞类型的差异表达数据。
x$avg_log2FC > 1：这部分代码用于检查基因的log2倍数变化是否大于1，也就是基因在STIM条件下的表达是否相对于CTRL条件显著上调。x$avg_log2FC < -1：这部分代码用于检查基因的log2倍数变化是否小于-1，也就是基因在STIM条件下的表达是否相对于CTRL条件显著下调。
- do.call(rbind, ...)：
rbind()函数用于将多个数据框或矩阵按行绑定在一起。do.call()用于将rbind应用到lapply()生成的结果列表中，将不同细胞类型的统计结果合并为一个矩阵或数据框。
即：
- 第一段代码统计了每个细胞类型中基因表达上调（log2倍数变化 > 1）的基因数量，并将所有细胞类型的统计结果整合为一个矩阵。
- 第二段代码统计了每个细胞类型中基因表达下调（log2倍数变化 < -1）的基因数量，并将所有细胞类型的统计结果整合为一个矩阵。
最终，这两段代码生成了两个矩阵，分别显示了每个细胞类型中基因表达显著上调或下调的基因数量。
执行完之后，每个细胞群只剩下了avg_log2FC < -1和>1的基因。
degs_list=lapply(names(degs),function(i){
#i=names(degs)[1]
...
})
degs_allcluster_type_df=do.call(rbind,degs_list)
write.csv(degs_allcluster_type_df,file = "./degs_allcluster_type_df.csv")
saveRDS(degs_allcluster_type_df,file = "./allclusters_degs_sce.markers.Rdata")- lapply(names(degs), function(i) { ... })
lapply(names(degs), function(i) {...})：对degs列表中的每个细胞类型名称执行指定的函数。i表示当前的细胞类型名称。
- x = degs[i]：提取degs列表中当前细胞类型i的差异表达结果。x现在包含了当前细胞类型的差异表达数据。
- res = x：将x赋值给res，这个步骤并没有改变数据，只是为变量
res赋值，方便后续操作。 - res$symbol = rownames(x):rownames(x)获取
x中的基因名称（通常是数据框的行名）。res$symbol = rownames(x)：将基因名称添加到res数据框中，作为一列，列名为symbol。 EnhancedVolcano()：使用EnhancedVolcano包生成火山图。
- lab = res$symbol：用基因符号（symbol列）作为图中的标签。
- x = 'avg_log2FC：使用avg_log2FC作为x轴变量，表示基因的log2倍数变化。
- y = 'p_val_adj：使用p_val_adj作为y轴变量，表示调整后的p值。
- pCutoff = 0.05：设置p值的阈值，p值小于0.05的基因被认为是显著的。
- FCcutoff = 1：设置log2倍数变化的阈值，大于1或小于-1的基因被认为是显著上调或下调。
- do.call(rbind, degs_list)：将所有细胞类型的差异表达数据按行绑定在一起，生成一个整合了所有细胞类型的差异表达数据框degs_allcluster_type_df。
生成的数据框degs_allcluster_type_df，包含了所有细胞分群的差异基因
这段代码将每个细胞分群都花了一张差异基因的火山图，并将结果保存在allclusters_degs_sce.markers.Rdata中。如
deg_for_B Activated_by_volcano.pdf
deg_for_B_by_volcano.pdf
p2=scRNAtoolVis::jjVolcano(diffData = degs_allcluster_type_df,...)这里安装了junjunlab/scRNAtoolVis包，这是一个用于单细胞RNA测序数据可视化的工具包。
scRNAtoolVis包中的一个函数，用于生成增强版火山图，适用于大规模的差异表达基因数据。
diffData = degs_allcluster_type_df：输入数据，即前面整合的所有细胞类型的差异表达数据框。log2FC.cutoff = 2：设定log2倍数变化的阈值，大于2或小于-2的基因被认为是显著上调或下调。adjustP.cutoff = 0.05：设定调整后的p值的阈值，小于0.05的基因被认为是显著的。
最后生成了
allcelltypes_degs_volcano.pdf
