单细胞GSVA分析该用什么数据？

经常有人问我单细胞GSVA分析应该用Seurat对象中的哪个数据，因为我此前的推文《单细胞转录组高级分析五：GSEA与GSVA分析》用的counts数据，后面有一篇推文《非人物种的GSEA&GSVA分析》用的是data数据。还有人推荐用scale.data，据说运行起来比counts数据快不少。很多人对此迷惑了，到底该用那个数据呢？

测试数据来源

测试数据来源于10x genomics官网的示例数据集。

数据下载链接：

https://cf.10xgenomics.com/samples/cell-exp/3.0.2/5k_pbmc_v3_nextgem/5k_pbmc_v3_nextgem_filtered_feature_bc_matrix.h5

保存文件名：pbmc.h5

运行时间对比

library(Seurat)
library(tidyverse)
library(GSVA)
library(msigdbr)

## 创建Seurat对象
#sc <- Read10X_h5("/home/james5woo/project/2020/2009_18pkgs/garnett/garnett_demo/pbmc.h5")
sc <- Read10X_h5("pbmc.h5")
sc <- sc[,sample(colnames(sc),1000)]  #将细胞数量随机抽样到1000个细胞
sc <- CreateSeuratObject(sc, project = "pbmc", min.cells = 10) #按较高标准过滤基因，主要是提高gsva的计算速度
dim(sc)
#[1] 12564  1000
sc <- NormalizeData(sc) %>% FindVariableFeatures() %>% ScaleData()

## 准备gsva分析的基因集(hallmark)
genesets <- msigdbr(species = "Homo sapiens", category = "H") %>% select("gs_name","gene_symbol") %>% as.data.frame()
genesets <- split(genesets$gene_symbol, genesets$gs_name)

## 测试counts数据
expr <- as.matrix(sc@assays$RNA@counts)
system.time({res.counts = gsva(expr, genesets, method="ssgsea", parallel.sz=10)}) 
#   用户    系统    流逝 
#196.838 533.948 207.484 

## 测试data数据
expr <- as.matrix(sc@assays$RNA@data)
system.time({res.data = gsva(expr, genesets, method="ssgsea", parallel.sz=10)}) 
#   用户    系统    流逝 
#195.250 501.749 197.244

## 测试scaledata数据
expr <- as.matrix(sc@assays$[e-mail])
system.time({res.scaledata = gsva(expr, genesets, method="ssgsea", parallel.sz=10)}) 
#  用户    系统    流逝 
#63.273 414.530 144.077

## 再次测试scaledata数据
sc2 <- ScaleData(sc, features = rownames(sc))
expr <- as.matrix(sc2@assays$[e-mail])
system.time({res.scaledata2 = gsva(expr, genesets, method="ssgsea", parallel.sz=10)}) 
#   用户    系统    流逝 
#235.689 721.989 272.767

从上面记录的运行时间来看，scale.data数据CPU的计算时间（63.273s）相比counts和data数据快了不少，这是因为默认的scale.data数据只有2000个基因。对所有基因scale之后再次测试scale.data数据，CPU的计算时间变成了235.689s，相比counts数据和data数据没有任何优势。

小结：scale.data数据并不能加快GSVA的运行时间。

分析结果对比

为了客观地对比不同数据运行GSVA之后的差异，我用pearson相关性热图给大家展示。**从左上到右下的对角线代表相同细胞用不同数据运行GSVA分析后结果的相关性。**为了节省计算时间，我只取了前100个细胞计算相关性。

1、counts数据与data数据的结果对比，相同细胞的pearson相关性等于1，说明运行结果完全相同。

2、counts数据与2000个基因的scale.data数据的结果对比，相同细胞的pearson相关性很低。

3、counts数据与所有基因的scale.data数据的结果对比，相同细胞的pearson相关性不高。

小结：GSVA分析使用counts数据和data数据没有差别，但是使用scale.data数据会影响结果。

减少基因数量可行吗？

写这篇推文时我突发奇想：使用高变基因来做GSVA分析可行吗？

# 筛选高变基因运行GSVA
EXPR <- as.matrix(sc@assays$RNA@counts)
for(i in seq(2000, 8000, 2000)){
  HVGs <- FindVariableFeatures(sc, nfeatures = i) %>% VariableFeatures()
  expr <- EXPR[HVGs,]
  res.hvg = gsva(expr, genesets, method="ssgsea", parallel.sz=10)
  coorda <- psych::corr.test(as.data.frame(res.counts[,1:100]), as.data.frame(res.hvg[,1:100]), method="pearson")
  pheatmap::pheatmap(coorda$r, cluster_rows = F, cluster_cols = F, display_numbers = F, show_rownames = F, 
                     show_colnames = F, filename = paste0("AllGenes_vs_",i,"HVGs.png"), width = 8, height = 7)
}

原始表达矩阵与2000高变基因表达矩阵的GSVA结果相关性

原始表达矩阵与4000高变基因表达矩阵的GSVA结果相关性

原始表达矩阵与6000高变基因表达矩阵的GSVA结果相关性

原始表达矩阵与8000高变基因表达矩阵的GSVA结果相关性

小结：不能使用高变基因的表达矩阵代替原始表达矩阵做GSVA分析。

交流探讨：如果您阅读此文有所疑惑，或有不同见解，亦或其他单细胞需求，可以点击阅读原文联系。

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注