单细胞基因集打分实操之irGSEA

原创 Immugent 生信宝库

Immugent在推文一文带你了解单细胞数据基因集打分的所有算法中，介绍了11种常用的对单细胞数据进行基因集打分的算法，其中提到了一个非常给力的R包：irGSEA。它内置了"AUCell", "UCell", "singscore", "ssgsea"四种算法，并且可以将多种基因集的富集分数矩阵直接保存到Seurat对象中，那么这期推文Immugent就给大家用示例数据来演示一下这个包的常见用法。

首先就是包的安装问题，可以参考我前一篇文章的代码。此外，Immugent已经咨询过包的作者，最好是用最新版(4.1) 的R上进行安装，但是小编用R version 4.0.5 (2021-03-31)也装上这个包了，只是在调用其中一个函数时出现了以下问题：

原因是这个包链接到的"GSVA"包函数是最新版的1.40.1；而4.05版本的R只能装低版本的“GSVA"包，能装上的最高版本是1.38.2的。

所以在低版本上即使能装上irGSEA包，其中的 “ssgsea” 算法是不能使用的，但是从专业角度看，"AUCell", "UCell", "singscore"更适用于单细胞数据，所以使用前3种算法就足够使用的了。

好了，下面开始展示具体的用法，为了每位伙伴都能很好的复现，小编就用Seuratdata的数据进行演示。

#-------------------------------1.准备数据---------------------------
library(SeuratData)
# view all available datasets
View(AvailableData())
# download 3k PBMCs from 10X Genomics
InstallData("pbmc3k")

library(Seurat)
# loading dataset
data("pbmc3k.final")
pbmc3k.final <- UpdateSeuratObject(pbmc3k.final)
pbmc3k.final
# plot
DimPlot(pbmc3k.final, reduction = "umap",
        group.by = "seurat_annotations",label = T) + NoLegend()

我们可以看到这是来源于人外周血的单细胞数据，包含各种免疫细胞，其主要分为三大块：B细胞，T细胞和髓系淋巴细胞。下面就是要使用上述的四种方法计算基因集的富集分数了，这里我们使用人的Hallmark数据集来进行展示。

#-------------------------------2.计算富集分数---------------------------
library(UCell)
library(irGSEA)

Idents(pbmc3k.final) <- pbmc3k.final$seurat_annotations #定义要分析的细胞亚群
pbmc3k.final <- irGSEA.score(object = pbmc3k.final, assay = "RNA",
                             slot = "data", seeds = 123, ncores = 1,
                             min.cells = 3, min.feature = 0,
                             custom = F, geneset = NULL, msigdb = T,
                             species = "Homo sapiens", category = "H",  
                             subcategory = NULL, geneid = "symbol",
                             method = c("AUCell", "UCell", "singscore",
                                        "ssgsea"),
                             aucell.MaxRank = NULL, ucell.MaxRank = NULL,
                             kcdf = 'Gaussian')

# 返回一个Seurat对象，富集分数矩阵存放在RNA外的assay中
Seurat::Assays(pbmc3k.final)

这个包最方便的就是可以通过一个函数直接使用4种算法，同时计算目的基因集的富集分数，而且函数还被作者进行优化，使得大大节约了所需计算时间。此外，因为内置了”msigdb“包，我们可以方便的计算人和其它模式物种的各种常见基因集，如GO, KEGG, Hallmark等。

在计算出所有细胞的通路富集分数后，接下来就需要找出差异的通路进行重点分析了。

#-------------------------------3.整合差异基因集---------------------------
result.dge <- irGSEA.integrate(object = pbmc3k.final,
                               group.by = "seurat_annotations",
                               metadata = NULL, col.name = NULL,
                               method = c("AUCell","UCell","singscore",
                                          "ssgsea"))

class(result.dge)
#> [1] "list"

为了评估基因集是否在某个细胞亚群中富集，irGSEA通过Wilcox检验计算富集分数矩阵中的差异基因集（差异基因的过滤标准为为校正后P值小于0.05）。然后，又通过RobustRankAggreg包中的秩聚合算法(RRA)对差异分析的结果进行综合评估，并筛选出在大部分基因集富集分析方法中显著富集的基因集（综合评估的过滤标准为P值小于0.05）。但是小编感觉这里只能对所有的细胞亚群同时进行统计学检验，而实际应用中我们往往只关注某两群细胞，因此如果未来能开发出像Seurat那样可以做到指定两群细胞进行统计学差异检验那就更完美了。

这个包的另一个亮点就是可以做出很多好看的图，下面就简单做出几个图展示一下

#-------------------------------4.可视化展示---------------------------
irGSEA.barplot.plot <- irGSEA.barplot(object = result.dge,
                                      method = c("AUCell", "UCell", "singscore",
                                                 "ssgsea"))
irGSEA.barplot.plot

halfvlnplot <- irGSEA.halfvlnplot(object = pbmc3k.final,
                                  method = "UCell",
                                  show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")
halfvlnplot

ridgeplot <- irGSEA.ridgeplot(object = pbmc3k.final,
                              method = "UCell",
                              show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")
ridgeplot

上面说到我们可以对经典数据集进行分析，但现实操作中我们往往想要对自己构建的目的基因集来进行富集分析，那么下面小编就进行展示如果构建自己的目的基因集；此外，值得一提的是"UCell", "singscore"函数可以通过正反向基因同时进行计算富集分数，那我们接下来就来看一下效果如何

#-------------------------------5.个性化分析---------------------------
markers <- list()
markers$TcellN <- c("CD3E+","CD3G+","ITGAE-", "ITGAM-")
markers$TcellY <- c("CD3E+","CD3G+","CD19-","CD79A-","CD79B-","FGFBP2-", "KLRF1-","ITGAE-", "ITGAM-")
pbmc3k.final <- irGSEA.score(object = pbmc3k.final, assay = "RNA", slot = "data", seeds = 123, ncores = 1,msigdb=F, custom = T, geneset = markers, method = c( "UCell", "singscore"), kcdf = 'Gaussian')

scatterplot <- irGSEA.density.scatterplot(object = pbmc3k.final,
                                          method = "singscore",
                                          show.geneset = c("TcellN","TcellY")
                                          reduction = "umap")
scatterplot

咋一看用这两种不同的基因集做出的富集分数很类似，但是细细看来就有不一样的发现。如右边那群髓系淋巴细胞，在没有排除B细胞的marker基因时的得分低于右侧排除后的，这是因为髓系淋巴细胞不表达B细胞的marker，用这种方式做出的综合评分会适当加分，而且这种影响随着基因数目增加而增加；此外，虽然T细胞群和其它群比较颜色上没有差异，但是打分范围右侧高于左侧，也可以认为，髓系的打分没有变，是B细胞群变得更低了，这其实是更有利于发现差异的。

好了，说到这，关于单细胞基因集打分的故事就说完了。接下来小编将会介绍另一个基于bulk数据进行打分的强大R包：IOBR，敬请期待！

继续滑动看下一个