百万级单细胞GSVA如何提速？

生信技能树

发布于 2025-01-01 11:10:36

16900

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

那么，当我们遇到大数据量的时候，如何加速单细胞的GSVA分析呢？

下面以 seurat 官网的 pbm3k 数据为例，进行演示：

首先是加载这个经典的数据集：

library(gplots)
library(ggplot2) 
library(clusterProfiler)
library(org.Hs.eg.db)
library(GSVA) # BiocManager::install('GSVA')
library(GSEABase)
# install.packages('devtools')
# devtools::install_github('satijalab/seurat-data')
library(SeuratData) #加载seurat数据集  
getOption('timeout')
options(timeout=10000)
# InstallData("pbmc3k")  
data("pbmc3k") 

sce <- pbmc3k.final  
table(Idents(sce)) 
p1 <- DimPlot(sce,label = T)
sce$celltype <- Idents(sce)

sce
table(Idents(sce))

现在这个sce是一个seurat对象，并且Idents为注释的细胞亚类：

做GSVA还需要一个基因集，去msigdb上薅一个下来：

library(msigdbr)
all_gene_sets <- msigdbr(species = "Homo sapiens", category='C7')
length(unique(table(all_gene_sets$gs_name)))
tail(sort(table(all_gene_sets$gs_name)))

# 转换成gsva分析需要的对象
gs <- split(all_gene_sets$gene_symbol,all_gene_sets$gs_name)
gs <- lapply(gs, unique)
gsc <- GeneSetCollection(mapply(function(geneIds, keggId) {
  GeneSet(geneIds, geneIdType=EntrezIdentifier(),
          collectionType=KEGGCollection(keggId),
          setName=keggId)
}, gs, names(gs)))
gsc
geneset <- gsc
geneset # 是一个 GeneSetCollection 对象

常规的gsva分析是在细胞水平做

尽管可以使用parallel.sz进行加速，但是数据量大的依然非常耗费计算机资源，耗费时间：

# 非常耗费计算机资源，耗费时间哦
es.max <- gsva(as.matrix( sce@assays$RNA$scale.data), geneset, mx.diff=FALSE, verbose=FALSE, parallel.sz=8)
save(es.max,file = 'gsva_all_pbmc.Rdata')

提速一：使用亚群水平进行GSVA分析

这样从原来的单个细胞变成注释后的亚群，样本量直接降级。使用 AverageExpression 函数提取每个细胞亚群的表达矩阵：

# v4 是 AverageExpression
av <- AggregateExpression(sce , group.by = "celltype", assays = "RNA") 
av <- av[[1]]
dim(av)
head(av) # 可以看到是整数矩阵

X <- av 
## X修改为as.matrix(X)
# 这里为 gsva 老版本代码
es.max <- gsva(as.matrix(X), geneset, mx.diff=FALSE, verbose=FALSE,  parallel.sz=8)
head(es.max)
dim(es.max)
save(es.max,file = 'gsva_celltype_pbmc.Rdata')

提速二：对单细胞进行了downsample处理

#抽样100个细胞测试展示
set.seed(1234)
table(Idents(sce))
sce_down <- subset(sce, downsample=100)
table(Idents(sce_down))

# 新版 gsva 软件
y <- as.matrix(sce_down@assays$RNA$scale.data)
gsvapar <- gsvaParam(y, geneset, maxDiff=TRUE,kcdf="Gaussian")
es.max <- gsva(as.matrix(sce_down@assays$RNA$scale.data), geneset, mx.diff=FALSE, verbose=FALSE,  parallel.sz=8)
head(es.max)
dim(es.max)
# [1]  50 751
save(es.max,file = 'gsva_celltype100_pbmc.Rdata')

降采样后，细胞亚群中细胞数大于100的取100个，小于100的为原来的细胞数：

得到为单个细胞水平的打分，单细细胞数已经是降采样后的：

抽样还可以解决的另外一个问题是计算机资源限制，现在的单细胞数据集很容易达到几十万甚至上百万的细胞数量，就容易出现下面的报错：

拿到了每个细胞的每个通路的打分后的可视化

这个时候，可以直接参考Seurat标准流程的可视化方法学，不过它默认出图并不好看，详见以前我们做的投票：可视化单细胞亚群的标记基因的5个方法，下面的5个基础函数相信大家都是已经烂熟于心了：

VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
FeaturePlot(pbmc, features = c("MS4A1", "CD79A"))
RidgePlot(pbmc, features = c("MS4A1", "CD79A"), ncol = 1)
DotPlot(pbmc, features = unique(features)) + RotatedAxis()
DoHeatmap(subset(pbmc, downsample = 100), features = features, size = 3)

也可以把打分矩阵（每个细胞或者单细胞亚群在每个通路）单独拿出来，使用热图的方法学去可视化，如下所示：

library(pheatmap)

# note：降采样后的的打分对每个亚群取gsva打分均值
score <- as.data.frame(t(es.max))
score$cellid <- rownames(score)
head(sce_down@meta.data)
sce_down$cellid <- rownames(sce_down@meta.data)
score <- merge(sce_down@meta.data[,c("cellid","celltype")],score, by="cellid")
head(score[,1:10])

score_ave <- limma::avereps(score[,-c(1,2)],ID=score$celltype)
score_ave <- t(score_ave)
head(score_ave)

pheatmap(score_ave[1:50, ],show_rownames = F)

随便找一篇高分文章，你可以看到里面大部分都是以细胞亚群水平进行GSVA打分展示：

如 2020年 Cell Research 发表的《Single-cell transcriptomics reveals regulators underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma》：