各个单细胞亚群的差异基因数量投射到umap图

生信技能树jimmy

发布于 2022-06-13 09:51:54

2.3K0

发布于 2022-06-13 09:51:54

文章被收录于专栏：单细胞天地单细胞天地

收到一个有意思的求助，希望可以把各个单细胞亚群的差异基因数量投射到umap图，如下所示：

各个单细胞亚群的差异基因数量投射到umap图

我简单读了一下文章，其实就降维聚类分群后，每个单细胞亚群在两个分组简单的做一下差异分析，有多少个单细胞亚群就做多少次差异分析，差异分析的上下调基因数量就是umap图里面的每个细胞的颜色情况。

我在 JCI Insight 2022 https://doi.org/10.1172/jci.insight.152616 文章里面也看到了类似的图：

因为每个单细胞亚群都有一个差异分析结果，所以也就是会有一个火山图等等，就跟常规表达量数据分析类似，公众号推文在：

单细胞的降维聚类分群标准分析这里我就不再赘述，也可以看基础10讲：

我们以大家熟知的pbmc3k数据集为例。大家先安装这个数据集对应的包，并且对它进行降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释 ，而且每个亚群找高表达量基因，都存储为Rdata文件。

# 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数据集并没有分组，所以我们没办法做差异分析。如果你一定要知道如何对每个单细胞亚群都在两个分组做一下差异分析并且统计上下调基因数量，也可以看前些天我们在《单细胞天地》的推文笔记：各个单细胞亚群独立在两个分组做差异分析

其实就是每个单细胞亚群都跑一下如下所示的示例代码：

# 假设是两个分组：
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图也是可以的：

> 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里面，如下所示：

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图

其实这个图呢，就是把前面的表格内容：

> 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

进行了简单的可视化，并没有太多的其它意义。

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2022-06-08，如有侵权请联系 cloudcommunity@tencent.com 删除

数据处理