irGSEA：基于秩次的单细胞基因集富集分析整合框架（中文教程）

一、背景

  1. 在单细胞转录组测序中，整合多样本数据进行分析逐渐成为一种趋势。
  2. 越来越多的研究者趋于使用未校正批次效应后的数据进行差异基因分析。
  3. 他们认为，校正批次效应后的数据，会掩盖部分真实的生物学差异。
    但校正批次效应后的数据是否能用于基因集富集分析，以及样本之间的批次效应是否会影响基因富集分析结果仍然是一个争论。
  4. 因此，我们重新审视现有的、主流的基因集富集方法，希望能找到受批次效应影响较小的基因集富集分析方法。
  5. 在这里，我们重新审视了9种常见的功能集打分方法：GSEA、GSVA、PLAGE、Zscore、AddModuleScore、ssGSEA、AUCell、UCell和singscore。

二、审视结果

  1. GSEA：执行前需对所有样本进行分组，然后基于分组去计算得到排序基因列表。这个过程中，我们需要考虑不同分组中样本构成的影响；
  2. GSVA：首先需要对所有样本中每个基因进行累积分布密度函数的核估计。这个过程中，需要考虑所有样本，容易受到样本背景信息的影响；
  3. PLAGE 和 z-score：首先需要对基因表达矩阵执行标准化处理。同样的，这个过程容易受样本构成的影响；
  4. AddModuleScore：Seurat包中的AddModuleScore函数，需要先计算基因集中所有基因的平均值，再根据平均值把表达矩阵切割成若干份，然后从切割后的每一份中随机抽取对照基因（基因集外的基因）作为背景值。因此，在整合不同样本的情况下，即使使用相同基因集为相同细胞打分，也会产生不同的富集评分；
  5. AUCell：基于单个样本中的基因表达排名（gene expression rank）,使用曲线下面积来评估输入基因集是否在单个样本的前5%表达基因内富集；
  6. UCell：基于单个样本的基因表达排名，使用Mann-Whitney U统计量计算单个样本的基因集富集评分；
    7.singscore：基于单个样本的基因表达排名，评估基因集远离中心的程度从而计算基因集富集评分；
    8.ssgsea：基于单个样本的基因表达排名，通过计算单个样本中基因集内和基因集外的经验累积分布函数之间的差值进而生成富集分数。然后，根据基因表达谱最大表达值和最小表达值的差值对富集分数进行标准化。这一步容易受样本构成的影响。

三、原理解析

1.纳入合适的方法：

我们淘汰了所有需要考虑样本组成的基因集富集分析方法，选取了基于单个样本的基因表达排名的基因集分析方法：AUCell、UCell和singscore。同时，我们也部分改进了ssGSEA，取消最后的标准化步骤，使之更贴近于单个样本的基因集富集分析。

2.构建基因集：

为了方便用户获取MSigDB数据库中预先定义好的基因集，我们内置了MSigDB包进行基因集的获取。同时，我们也支持多个物种的基因集获取，以及多种基因格式的表达矩阵的输入。除此之外，我们还允许用户自定义自己的基因集进行打分计算。如果用户的基因集中既包括正向的基因，也包括负向的基因。该基因集将默认只使用UCell包和singscore包进行打分处理。

3.输入对象和数据清洗：

我们允许直接输入单细胞表达矩阵或者Seurat对象。我们内置了Seurat包，可以将多种基因集的富集分数矩阵直接保存到Seurat对象中。同时，我们也支持过滤单细胞表达矩阵中所有细胞表达量为0的基因。当然，用户也可以自定义自己的过滤标准。合适的过滤指标可以改善富集分析的结果。

4.差异分析和综合评估：

为了评估基因集在某个细胞亚群中是否富集，我们通过多种基因集富集方法分别对单个细胞进行打分，并生成多个基因集富集分数矩阵。接着，我们通过wilcox检验计算各个基因集富集分数矩阵中每个细胞亚群差异表达的基因集。
单一的基因集富集分析方法不仅只能反映有限的信息，而且也容易带来误差。我们期待从多个角度解释复杂的生物学问题，并找到生物学问题中的共性部分。简单地为多种基因集富集分析方法的结果取共同交集，不仅容易得到少而保守的结果，而且忽略了富集分析方法中很多的其他信息，例如不同基因集的相对富集程度信息。
单纯地取共同交集无法体现基因集的相对富集程度，而我们希望目标基因集在大部分富集分析方法中都是富集且富集程度没有明显差异。因此，我们通过RobustRankAggreg包中的秩聚合算法(robust rank aggregation, RRA)对差异分析的结果进行综合评估，筛选出在大部分基因集富集分析方法中都显著富集的基因集。

5.可视化展示：

我们内置了多种可视化的函数，不仅允许用户通过热图、气泡图、柱状图和upset图展示它们综合结果，而且允许用户通过密度散点图、半小提琴图、山峦图和密度热图展示目标基因集在具体富集分析方法中的表达水平和数据分布。其中，热图、upset图、密度热图主要由ComplexHeatmap包生成；气泡图和柱状图主要由ggplot2包、ggtree包、aplot包生成；密度散点图主要由Seurat包和Nebulosa包生成，半小提琴图主要由Seurat包和gghalves包生成；山峦图主要由Seurat包和ggridges包生成。最后，为了方便用户将可视化结果与其他的ggplot2对象进行拼图操作，我们也通过ggplotify包把输出结果转换为ggplot2对象。最后，为了方便用户的实现，我们开发了R包irGSEA (https://github.com/chuiqin/irGSEA/)，集成了上述工作流程。

四、使用教程

1.安装前的准备

在安装irGSEA 需要查看自己是否安装了以下这些包

# install packages from CRAN
cran.packages <- c("msigdbr", "dplyr", "purrr", "stringr","magrittr",
                   "RobustRankAggreg", "tibble", "reshape2", "ggsci",
                   "tidyr", "aplot", "ggfun", "ggplotify", "ggridges", 
                   "gghalves", "Seurat", "SeuratObject", "methods", 
                   "devtools", "BiocManager","data.table","doParallel",
                   "doRNG")
if (!requireNamespace(cran.packages, quietly = TRUE)) { 
    install.packages(cran.packages, ask = F, update = F)
}

# install packages from Bioconductor
bioconductor.packages <- c("GSEABase", "AUCell", "SummarizedExperiment", 
                           "singscore", "GSVA", "ComplexHeatmap", "ggtree", 
                           "Nebulosa")
if (!requireNamespace(bioconductor.packages, quietly = TRUE)) { 
    BiocManager::install(bioconductor.packages, ask = F, update = F)
}
# install packages from Github
if (!requireNamespace("UCell", quietly = TRUE)) { 
    devtools::install_github("carmonalab/UCell")
}
if (!requireNamespace("irGSEA", quietly = TRUE)) { 
    devtools::install_github("chuiqin/irGSEA")
}

2.加载示例数据

我们通过SeuratData包加载示例数据集（注释好的PBMC数据集）

# devtools::install_github('satijalab/seurat-data')
library(SeuratData)
# view all available datasets
View(AvailableData())
# download 3k PBMCs from 10X Genomics
# 示例数据有点大，建议在网速好的时候下载，如果失败了，请重启R后多尝试几遍
InstallData("pbmc3k")
# the details of pbmc3k.final
?pbmc3k.final

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

顺手设置一下数据集的ident

# set cluster to idents
Idents(pbmc3k.final) <- pbmc3k.final$seurat_annotations

3.加载R包

这一步出错的话，要看一下前面的包有没有装好

library(UCell)
library(irGSEA)

4.计算富集分数

当你的ncore设置大于1的时候，发生下面的错误：Error (Valid ‘mctype’: ‘snow’ or ‘doMC’)，你应该检查一下你的AUCell 版本，确保版本大于等于1.14 。如果你比较懒，那你直接把ncore设置为1也是可以的，只是运行速度会稍微慢一点。

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')
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper strucutre
#> Ensuring feature names don't have underscores or pipes
#> Object representation is consistent with the most current Seurat version
#> Calculate AUCell scores
#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')

#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')
#> Finish calculate AUCell scores
#> Calculate UCell scores
#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')

#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')
#> Finish calculate UCell scores
#> Calculate singscore scores
#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')

#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')
#> Finish calculate singscore scores
#> Calculate ssgsea scores
#> Warning in .filterFeatures(expr, method): 1 genes with constant expression
#> values throuhgout the samples.

#> Warning in .filterFeatures(expr, method): Feature names cannot have underscores
#> ('_'), replacing with dashes ('-')

#> Warning in .filterFeatures(expr, method): Feature names cannot have underscores
#> ('_'), replacing with dashes ('-')
#> Finish calculate ssgsea scores

# 返回一个Seurat对象，富集分数矩阵存放在RNA外的assay中
Seurat::Assays(pbmc3k.final)
#> [1] "RNA"       "AUCell"    "UCell"     "singscore" "ssgsea"

5.整合差异基因集

# Wlicox test is perform to all enrichment score matrixes and gene sets
# with adjusted p value &lt; 0.05 are used to integrated through RRA.
# Among them, Gene sets with p value &lt; 0.05 are statistically
# significant and common differential in all gene sets enrichment analysis
# methods. All results are saved in a list.
result.dge <- irGSEA.integrate(object = pbmc3k.final, 
                               group.by = "seurat_annotations",
                               metadata = NULL, col.name = NULL,
                               method = c("AUCell","UCell","singscore",
                                          "ssgsea"))
#> Calculate differential gene set : AUCell
#> Calculate differential gene set : UCell
#> Calculate differential gene set : singscore
#> Calculate differential gene set : ssgsea
class(result.dge)
#> [1] "list"

6.可视化展示

1）. 全局展示

①.热图

热图展示了综合评价中具体基因集在每个细胞亚群是否具有统计学意义差异；
其中，浅蓝色的格子无统计学差异，红色的格子具有统计学差异。格子中的星号越多，格子的P值越小；
左边的聚类树代表不同基因集在不同细胞亚群中表达模式的相似性；
上方的条形图分别代表不同的细胞亚群，以及差异基因集在细胞亚群中是呈现上调还是下调趋势；
你还可以把method从'RRA"换成“ssgsea”，展示特定基因集富集分析方法中差异上调或差异下调的基因集；

irGSEA.heatmap.plot <- irGSEA.heatmap(object = result.dge, 
                                      method = "RRA",
                                      top = 50, 
                                      show.geneset = NULL)
irGSEA.heatmap.plot
热图

②.气泡图

气泡图展示了综合评价中具体基因集在每个细胞亚群是否具有统计学意义差异；
其中，浅蓝色的点无统计学差异，红色的点具有统计学差异。点越大，P值越小；
左边的聚类树代表不同基因集在不同细胞亚群中表达模式的相似性；
上方的条形图分别代表不同的细胞亚群，以及差异基因集在目标细胞亚群中是呈现上调还是下调趋势；
你还可以把method从'RRA"换成“ssgsea”，展示特定基因集富集分析方法中差异上调或差异下调的基因集；
如果运行中提示“ error (argument “caller_env” is missing, with no default)”，请卸载了ggtree包，并重新安装remotes::install_github(”YuLab-SMU/ggtree“)

irGSEA.bubble.plot <- irGSEA.bubble(object = result.dge, 
                                    method = "RRA", 
                                    top = 50)
irGSEA.bubble.plot
气泡图

③.upset plot

upset图展示了综合评估中每个细胞亚群具有统计学意义差异的基因集的数目，以及不同细胞亚群之间具有交集的差异基因集数目；
左边不同颜色的条形图代表不同的细胞亚群；
上方的条形图代表具有交集的差异基因集的数目；
中间的气泡图单个点代表单个细胞亚群，多个点连线代表多个细胞亚群取交集.；
如果运行代码发生以下warning“the condition has length > 1 and
only the first element will be used”。这个警告可忽略。

irGSEA.upset.plot <- irGSEA.upset(object = result.dge, 
                                  method = "RRA")
#> Warning in if (as.character(ta_call[[1]]) == "upset_top_annotation") {: the
#> condition has length > 1 and only the first element will be used
irGSEA.upset.plot
upset图

④.堆叠条形图

堆叠柱状图具体展示每种基因集富集分析方法中每种细胞亚群中上调、下调和没有统计学差异的基因集数目；
上方的条形代表每个亚群中不同方法中差异的基因数目，红色代表上调的差异基因集，蓝色代表下调的差异基因集；
中间的柱形图代表每个亚群中不同方法中上调、下调和没有统计学意义的基因集的比例；

irGSEA.barplot.plot <- irGSEA.barplot(object = result.dge,
                                      method = c("AUCell", "UCell", "singscore",
                                                 "ssgsea"))
irGSEA.barplot.plot
image.png

2）. 局部展示

①.密度散点图

密度散点图将基因集的富集分数和细胞亚群在低维空间的投影结合起来，展示了特定基因集在空间上的表达水平。其中，颜色越黄，代表富集分数越高；

scatterplot <- irGSEA.density.scatterplot(object = pbmc3k.final,
                             method = "UCell",
                             show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE",
                             reduction = "umap")
scatterplot
image.png

②.半小提琴图

半小提琴图同时以小提琴图（左边）和箱线图（右边）进行展示。不同颜色代表不同的细胞亚群；

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

③.山峦图

山峦图中上方的核密度曲线展示了数据的主要分布，下方的条形编码图展示了细胞亚群具体的数量。不同颜色代表不同的细胞亚群，而横坐标代表不同的表达水平；

ridgeplot <- irGSEA.ridgeplot(object = pbmc3k.final,
                              method = "UCell",
                              show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")
ridgeplot
#> Picking joint bandwidth of 0.00533
image.png

④.密度热图

密度热图展示了具体差异基因在不同细胞亚群中的表达和分布水平。颜色越红，代表富集分数越高；

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

五、拓展：自己如何做一个R包

如果这篇推文只到教程结束，那一切都将索然无味。
我也没有必要在前面水那么多字数的背景，并详细介绍每一步我是怎么做的。
因为我的单细胞R包框架是可以被模仿的，并且只需要以下3步：
1.针对一个生物学问题收集多种解决方法；
2.对多种解决方法的结果进行一个综合评估；
3.对综合评估的结果进行可视化展示；
举个例子，你可以做一个计算细胞亚群差异基因的单细胞R包：
1.收集多种单细胞差异基因集计算方法（例如wilicox,AUC,MAST等）；
2.分别计算，每种方法中每种细胞亚群中差异表达的基因，用一种综合评估的方法（例如RRA算法）评估这些差异分析方法的结果，找出在大部分方法中都具有差异的基因；
3.寻找合适的可视化方法（例如热图、散点图、气泡图等）可视化展示综合评估的结果；
再举个例子，你可以做一个计算TCGA中不同分子亚型免疫浸润差异的R包：
1.收集多种免疫浸润的计算方法（例如cibersortX,estimate,xcell,ssgsea,等）；
2.分别计算，每种方法中每种分子亚型中差异表达的免疫细胞群，用一种综合评估的方法（例如RRA算法）评估这些差异分析方法的结果，找出在大部分方法中都具有差异的基因；
3.寻找合适的可视化方法（例如热图、散点图、气泡图等）可视化展示综合评估的结果；
再再再举个例子，收集多种单细胞转录因子调控网络的计算方法（例如SCENIC，DoRothEA等），收集多种转录组差异分析的方法（例如limma-voom，edgeR, DESeq2等）...然后按照上文介绍的流程，如法炮制即可。
简单总结一下这个套路：多种方法 + 综合评估 + 可视化 = R包
如果你已经厌倦了GEO/TCGA的数据挖掘，可以尝试一下开发一两个有趣的R包。
当然，只有框架还是不够，你还需要有一定的R语言基础，以及R包开发的经验。
关于R包开发，推荐你们去看生信技能树在B站的视频，以及Hadley Wickham写的《R包开发》这本书，相信会让你减少很多摸索的时间。在此，我也要感谢一下生信技能树组织的《R包开发》课程，以及永和大神的指点。

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 206,378评论 6赞 481
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 88,356评论 2赞 382
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 152,702评论 0赞 342
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 55,259评论 1赞 279
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 64,263评论 5赞 371
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,036评论 1赞 285
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,349评论 3赞 400
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,979评论 0赞 259
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 43,469评论 1赞 300
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,938评论 2赞 323
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,059评论 1赞 333
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,703评论 4赞 323
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,257评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,262评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,485评论 1赞 262
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,501评论 2赞 354
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,792评论 2赞 345

推荐阅读更多精彩内容