收到一个有意思的求助，希望可以把各个单细胞亚群的差异基因数量投射到umap图 ，如下所示：
各个单细胞亚群的差异基因数量投射到umap图
我简单读了一下文章，其实就降维聚类分群后，每个单细胞亚群在两个分组简单的做一下差异分析，有多少个单细胞亚群就做多少次差异分析，差异分析的上下调基因数量就是umap图里面的每个细胞的颜色情况。
我在 JCI Insight 2022 https://doi.org/10.1172/jci.insight.152616 文章里面也看到了类似的图：
因为每个单细胞亚群都有一个差异分析结果，所以也就是会有一个火山图等等，就跟常规表达量数据分析类似，公众号推文在：
- 解读GEO数据存放规律及下载，一文就够
- 解读SRA数据库规律一文就够
- 从GEO数据库下载得到表达矩阵 一文就够
- GSEA分析一文就够（单机版 R语言版）
- 根据分组信息做差异分析- 这个一文不够的
- 差异分析得到的结果注释一文就够
单细胞的降维聚类分群标准分析这里我就不再赘述，也可以看基础10讲：
- 01. 上游分析流程
- 02.课题多少个样品，测序数据量如何
- 03. 过滤不合格细胞和基因（数据质控很重要）
- 04. 过滤线粒体核糖体基因
- 05. 去除细胞效应和基因效应
- 06.单细胞转录组数据的降维聚类分群
- 07.单细胞转录组数据处理之细胞亚群注释
- 08.把拿到的亚群进行更细致的分群
- 09.单细胞转录组数据处理之细胞亚群比例比较
我们以大家熟知的pbmc3k数据集为例。大家先安装这个数据集对应的包，并且对它进行降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释 ，而且每个亚群找高表达量基因，都存储为Rdata文件。
代码语言：javascript复制# 0.安装R包 ----
# InstallData("pbmc3k")
library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")
data("pbmc3k")
sce <- pbmc3k.final
library(Seurat)
table(Idents(sce))
DimPlot(sce,label = T)
library(future)
# check the current active plan
plan()
plan("multiprocess", workers = 4)
plan()
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE,
min.pct = 0.25,
thresh.use = 0.25)
DT::datatable(sce.markers)
pro='markers'
write.csv(sce.markers,file=paste0(pro,'_sce.markers.csv'))
save(sce.markers,file = paste0(pro, 'sce.markers.Rdata'))
这个时候，因为它pbmc3k数据集并没有分组，所以我们没办法做差异分析。如果你一定要知道如何对每个单细胞亚群都在两个分组做一下差异分析并且统计上下调基因数量，也可以看前些天我们在《单细胞天地》的推文笔记：各个单细胞亚群独立在两个分组做差异分析
其实就是每个单细胞亚群都跑一下如下所示的示例代码 ：
代码语言：javascript复制# 假设是两个分组：
Idents(sce) = sce$group
table(Idents(sce))
deg = FindMarkers(sce,ident.1 = 'group1',
ident.2 = 'group1')
head(deg[order(deg$p_val),])
table(Idents(sce))
library(EnhancedVolcano)
res=deg
head(res)
EnhancedVolcano(res,
lab = rownames(res),
x = 'avg_log2FC',
y = 'p_val_adj')
虽然我们没办法跑差异分析，但是统计每个单细胞亚群的标记基因，也是可以的啊！
我们把每个单细胞亚群的标记基因数量投射到umap图也是可以的：
代码语言：javascript复制> as.data.frame(table(sce.markers$cluster))
Var1 Freq
1 Naive CD4 T 162
2 Memory CD4 T 176
3 CD14 Mono 391
4 B 147
5 CD8 T 162
6 FCGR3A Mono 608
7 NK 364
8 DC 633
9 Platelet 242
代码比较简单，把每个单细胞亚群的标记基因数量首先加入到seurat对象的meta.data里面，如下所示：
代码语言：javascript复制sce$celltype = Idents(sce)
df = as.data.frame(table(sce.markers$cluster))
colnames(df) = c('celltype','NofDEGs')
sce$NofDEGs = df[match(sce$celltype,df$celltype),2]
library(patchwork)
DimPlot(sce,label = T) FeaturePlot(sce,'NofDEGs')
如下所示：
把每个单细胞亚群的标记基因数量投射到umap图
其实这个图呢，就是把前面的表格内容：
代码语言：javascript复制> as.data.frame(table(sce.markers$cluster))
Var1 Freq
1 Naive CD4 T 162
2 Memory CD4 T 176
3 CD14 Mono 391
4 B 147
5 CD8 T 162
6 FCGR3A Mono 608
7 NK 364
8 DC 633
9 Platelet 242
进行了简单的可视化，并没有太多的其它意义。
