10X单细胞（10X空间转录组）基因集打分方法汇总

作者: 单细胞空间交响乐 | 来源:发表于2021-10-25 16:30 被阅读0次

hello，大家好，今天周一，做一个简单的汇总，主要针对基因集打分，无论是单细胞还是空间，基因集打分都很普遍，方法也很多，我们来总结一下。

研究背景

在单细胞转录组测序中，整合多样本数据进行分析逐渐成为一种趋势。

越来越多的研究者趋于使用未校正批次效应后的数据进行差异基因分析。

校正批次效应后的数据，会掩盖部分真实的生物学差异。但校正批次效应后的数据是否能用于基因集富集分析，以及样本之间的批次效应是否会影响基因富集分析结果仍然是一个争论。

因此，我们重新审视现有的、主流的基因集富集方法，希望能找到受批次效应影响较小的基因集富集分析方法。

在这里，重新评估了9种常见的功能集打分方法：GSEA、GSVA、PLAGE、Zscore、AddModuleScore、ssGSEA、AUCell、UCell和singscore。

结果汇总

GSEA：执行前需对所有样本进行分组，然后基于分组去计算得到排序基因列表。这个过程中，需要考虑不同分组中样本构成的影响；

GSVA：首先需要对所有样本中每个基因进行累积分布密度函数的核估计。这个过程中，需要考虑所有样本，容易受到样本背景信息的影响；

PLAGE 和 z-score：首先需要对基因表达矩阵执行标准化处理。同样的，这个过程容易受样本构成的影响；

AddModuleScore：Seurat包中的AddModuleScore函数，需要先计算基因集中所有基因的平均值，再根据平均值把表达矩阵切割成若干份，然后从切割后的每一份中随机抽取对照基因（基因集外的基因）作为背景值。因此，在整合不同样本的情况下，即使使用相同基因集为相同细胞打分，也会产生不同的富集评分；

AUCell：基于单个样本中的基因表达排名（gene expression rank）,使用曲线下面积来评估输入基因集是否在单个样本的前5%表达基因内富集；

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，集成了上述工作流程。

示例代码

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()

图片

顺手设置一下数据集的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 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"` </pre>

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图展示了综合评估中每个细胞亚群具有统计学意义差异的基因集的数目，以及不同细胞亚群之间具有交集的差异基因集数目；左边不同颜色的条形图代表不同的细胞亚群；上方的条形图代表具有交集的差异基因集的数目；中间的气泡图单个点代表单个细胞亚群，多个点连线代表多个细胞亚群取交集.

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

图片

④.堆叠条形图

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

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

图片

2）. 局部展示

①.密度散点图

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

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

图片

②.半小提琴图

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

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
#> Picking joint bandwidth of 0.00533

图片

④.密度热图

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

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

图片

生活很好，有你更好