hello，大家好，今天周一，做一個(gè)簡(jiǎn)單的匯總岩梳，主要針對(duì)基因集打分啃匿，無(wú)論是單細(xì)胞還是空間，基因集打分都很普遍雷蹂，方法也很多伟端，我們來(lái)總結(jié)一下。

研究背景

在單細(xì)胞轉(zhuǎn)錄組測(cè)序中匪煌，整合多樣本數(shù)據(jù)進(jìn)行分析逐漸成為一種趨勢(shì)责蝠。

越來(lái)越多的研究者趨于使用未校正批次效應(yīng)后的數(shù)據(jù)進(jìn)行差異基因分析。

校正批次效應(yīng)后的數(shù)據(jù)萎庭，會(huì)掩蓋部分真實(shí)的生物學(xué)差異霜医。但校正批次效應(yīng)后的數(shù)據(jù)是否能用于基因集富集分析，以及樣本之間的批次效應(yīng)是否會(huì)影響基因富集分析結(jié)果仍然是一個(gè)爭(zhēng)論驳规。

因此肴敛，我們重新審視現(xiàn)有的、主流的基因集富集方法，希望能找到受批次效應(yīng)影響較小的基因集富集分析方法医男。

在這里砸狞，重新評(píng)估了9種常見(jiàn)的功能集打分方法：GSEA、GSVA镀梭、PLAGE刀森、Zscore、AddModuleScore报账、ssGSEA研底、AUCell、UCell和singscore笙什。

結(jié)果匯總

GSEA：執(zhí)行前需對(duì)所有樣本進(jìn)行分組飘哨，然后基于分組去計(jì)算得到排序基因列表。這個(gè)過(guò)程中琐凭，需要考慮不同分組中樣本構(gòu)成的影響芽隆；

GSVA：首先需要對(duì)所有樣本中每個(gè)基因進(jìn)行累積分布密度函數(shù)的核估計(jì)。這個(gè)過(guò)程中统屈，需要考慮所有樣本胚吁，容易受到樣本背景信息的影響；

PLAGE 和 z-score：首先需要對(duì)基因表達(dá)矩陣執(zhí)行標(biāo)準(zhǔn)化處理愁憔。同樣的腕扶，這個(gè)過(guò)程容易受樣本構(gòu)成的影響；

AddModuleScore：Seurat包中的AddModuleScore函數(shù)吨掌，需要先計(jì)算基因集中所有基因的平均值半抱，再根據(jù)平均值把表達(dá)矩陣切割成若干份，然后從切割后的每一份中隨機(jī)抽取對(duì)照基因（基因集外的基因）作為背景值膜宋。因此窿侈，在整合不同樣本的情況下，即使使用相同基因集為相同細(xì)胞打分秋茫，也會(huì)產(chǎn)生不同的富集評(píng)分史简；

AUCell：基于單個(gè)樣本中的基因表達(dá)排名（gene expression rank）,使用曲線下面積來(lái)評(píng)估輸入基因集是否在單個(gè)樣本的前5%表達(dá)基因內(nèi)富集；

UCell：基于單個(gè)樣本的基因表達(dá)排名肛著，使用Mann-Whitney U統(tǒng)計(jì)量計(jì)算單個(gè)樣本的基因集富集評(píng)分圆兵；7.singscore：基于單個(gè)樣本的基因表達(dá)排名，評(píng)估基因集遠(yuǎn)離中心的程度從而計(jì)算基因集富集評(píng)分枢贿；8.ssgsea：基于單個(gè)樣本的基因表達(dá)排名殉农，通過(guò)計(jì)算單個(gè)樣本中基因集內(nèi)和基因集外的經(jīng)驗(yàn)累積分布函數(shù)之間的差值進(jìn)而生成富集分?jǐn)?shù)。然后局荚，根據(jù)基因表達(dá)譜最大表達(dá)值和最小表達(dá)值的差值對(duì)富集分?jǐn)?shù)進(jìn)行標(biāo)準(zhǔn)化超凳。這一步容易受樣本構(gòu)成的影響。

分析原理

原理解析

1.納入合適的方法：

淘汰了所有需要考慮樣本組成的基因集富集分析方法，選取了基于單個(gè)樣本的基因表達(dá)排名的基因集分析方法：AUCell聪建、UCell和singscore钙畔。同時(shí)，也部分改進(jìn)了ssGSEA金麸，取消最后的標(biāo)準(zhǔn)化步驟擎析，使之更貼近于單個(gè)樣本的基因集富集分析。

2.構(gòu)建基因集：

為了方便獲取MSigDB數(shù)據(jù)庫(kù)中預(yù)先定義好的基因集挥下，內(nèi)置MSigDB包進(jìn)行基因集的獲取揍魂。同時(shí)，也支持多個(gè)物種的基因集獲取棚瘟，以及多種基因格式的表達(dá)矩陣的輸入现斋。除此之外，還允許自定義自己的基因集進(jìn)行打分計(jì)算偎蘸。如果基因集中既包括正向的基因庄蹋，也包括負(fù)向的基因。該基因集將默認(rèn)只使用UCell包和singscore包進(jìn)行打分處理迷雪。

3.輸入對(duì)象和數(shù)據(jù)清洗：

允許直接輸入單細(xì)胞表達(dá)矩陣或者Seurat對(duì)象限书。可以內(nèi)置Seurat包章咧，可以將多種基因集的富集分?jǐn)?shù)矩陣直接保存到Seurat對(duì)象中倦西。同時(shí)，也支持過(guò)濾單細(xì)胞表達(dá)矩陣中所有細(xì)胞表達(dá)量為0的基因赁严。當(dāng)然扰柠，也可以自定義自己的過(guò)濾標(biāo)準(zhǔn)。合適的過(guò)濾指標(biāo)可以改善富集分析的結(jié)果疼约。

4.差異分析和綜合評(píng)估：

為了評(píng)估基因集在某個(gè)細(xì)胞亞群中是否富集卤档，通過(guò)對(duì)多種基因集富集方法分別對(duì)單個(gè)細(xì)胞進(jìn)行打分，并生成多個(gè)基因集富集分?jǐn)?shù)矩陣忆谓。接著裆装，通過(guò)wilcox檢驗(yàn)計(jì)算各個(gè)基因集富集分?jǐn)?shù)矩陣中每個(gè)細(xì)胞亞群差異表達(dá)的基因集踱承。單一的基因集富集分析方法不僅只能反映有限的信息倡缠，而且也容易帶來(lái)誤差。期待從多個(gè)角度解釋復(fù)雜的生物學(xué)問(wèn)題茎活，并找到生物學(xué)問(wèn)題中的共性部分昙沦。簡(jiǎn)單地為多種基因集富集分析方法的結(jié)果取共同交集，不僅容易得到少而保守的結(jié)果载荔，而且忽略了富集分析方法中很多的其他信息盾饮，例如不同基因集的相對(duì)富集程度信息。單純地取共同交集無(wú)法體現(xiàn)基因集的相對(duì)富集程度，而目標(biāo)基因集在大部分富集分析方法中都是富集且富集程度沒(méi)有明顯差異丘损。因此普办，通過(guò)RobustRankAggreg包中的秩聚合算法(robust rank aggregation, RRA)對(duì)差異分析的結(jié)果進(jìn)行綜合評(píng)估，篩選出在大部分基因集富集分析方法中都顯著富集的基因集徘钥。

5.可視化展示：

內(nèi)置了多種可視化的函數(shù)衔蹲，不僅允許通過(guò)熱圖、氣泡圖呈础、柱狀圖和upset圖展示它們綜合結(jié)果舆驶，而且允許通過(guò)密度散點(diǎn)圖、半小提琴圖而钞、山巒圖和密度熱圖展示目標(biāo)基因集在具體富集分析方法中的表達(dá)水平和數(shù)據(jù)分布沙廉。其中，熱圖臼节、upset圖撬陵、密度熱圖主要由ComplexHeatmap包生成；氣泡圖和柱狀圖主要由ggplot2包网缝、ggtree包袱结、aplot包生成；密度散點(diǎn)圖主要由Seurat包和Nebulosa包生成途凫，半小提琴圖主要由Seurat包和gghalves包生成垢夹；山巒圖主要由Seurat包和ggridges包生成。最后维费，為了方便將可視化結(jié)果與其他的ggplot2對(duì)象進(jìn)行拼圖操作果元，也可以通過(guò)ggplotify包把輸出結(jié)果轉(zhuǎn)換為ggplot2對(duì)象。最后犀盟，為了方便而晒，開(kāi)發(fā)了R包irGSEA，集成了上述工作流程阅畴。

示例代碼

1.安裝前的準(zhǔn)備

在安裝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.加載示例數(shù)據(jù)

通過(guò)SeuratData包加載示例數(shù)據(jù)集（注釋好的PBMC數(shù)據(jù)集）

# devtools::install_github('satijalab/seurat-data')
library(SeuratData)
# view all available datasets
View(AvailableData())
# download 3k PBMCs from 10X Genomics
# 示例數(shù)據(jù)有點(diǎn)大倡怎，建議在網(wǎng)速好的時(shí)候下載，如果失敗了贱枣，請(qǐng)重啟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()

圖片

順手設(shè)置一下數(shù)據(jù)集的ident

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

3.加載R包

這一步出錯(cuò)的話监署，要看一下前面的包有沒(méi)有裝好

library(UCell)
library(irGSEA)

4.計(jì)算富集分?jǐn)?shù)

當(dāng)你的ncore設(shè)置大于1的時(shí)候，發(fā)生下面的錯(cuò)誤：Error (Valid ‘mctype’: ‘snow’ or ‘doMC’)纽哥，你應(yīng)該檢查一下你的AUCell 版本钠乏，確保版本大于等于1.14 。如果為了方便春塌，直接把ncore設(shè)置為1也是可以的晓避，只是運(yùn)行速度會(huì)稍微慢一點(diǎn)簇捍。

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

# 返回一個(gè)Seurat對(duì)象，富集分?jǐn)?shù)矩陣存放在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 < 0.05 are used to integrated through RRA.
# Among them, Gene sets with p value < 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íng)價(jià)中具體基因集在每個(gè)細(xì)胞亞群是否具有統(tǒng)計(jì)學(xué)意義差異俏拱；其中暑塑，淺藍(lán)色的格子無(wú)統(tǒng)計(jì)學(xué)差異，紅色的格子具有統(tǒng)計(jì)學(xué)差異锅必。格子中的星號(hào)越多梯投，格子的P值越小况毅；左邊的聚類樹(shù)代表不同基因集在不同細(xì)胞亞群中表達(dá)模式的相似性分蓖；上方的條形圖分別代表不同的細(xì)胞亞群，以及差異基因集在細(xì)胞亞群中是呈現(xiàn)上調(diào)還是下調(diào)趨勢(shì)尔许；你還可以把method從'RRA"換成“ssgsea”么鹤，展示特定基因集富集分析方法中差異上調(diào)或差異下調(diào)的基因集；

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

圖片

②.氣泡圖

氣泡圖展示了綜合評(píng)價(jià)中具體基因集在每個(gè)細(xì)胞亞群是否具有統(tǒng)計(jì)學(xué)意義差異味廊；其中蒸甜，淺藍(lán)色的點(diǎn)無(wú)統(tǒng)計(jì)學(xué)差異，紅色的點(diǎn)具有統(tǒng)計(jì)學(xué)差異余佛。點(diǎn)越大柠新，P值越小辉巡；左邊的聚類樹(shù)代表不同基因集在不同細(xì)胞亞群中表達(dá)模式的相似性恨憎；上方的條形圖分別代表不同的細(xì)胞亞群，以及差異基因集在目標(biāo)細(xì)胞亞群中是呈現(xiàn)上調(diào)還是下調(diào)趨勢(shì)郊楣；你還可以把method從'RRA"換成“ssgsea”憔恳，展示特定基因集富集分析方法中差異上調(diào)或差異下調(diào)的基因集；如果運(yùn)行中提示“ error (argument “caller_env” is missing, with no default)”净蚤，請(qǐng)卸載了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圖展示了綜合評(píng)估中每個(gè)細(xì)胞亞群具有統(tǒng)計(jì)學(xué)意義差異的基因集的數(shù)目，以及不同細(xì)胞亞群之間具有交集的差異基因集數(shù)目今瀑；左邊不同顏色的條形圖代表不同的細(xì)胞亞群程梦；上方的條形圖代表具有交集的差異基因集的數(shù)目；中間的氣泡圖單個(gè)點(diǎn)代表單個(gè)細(xì)胞亞群橘荠，多個(gè)點(diǎn)連線代表多個(gè)細(xì)胞亞群取交集.

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

圖片

④.堆疊條形圖

堆疊柱狀圖具體展示每種基因集富集分析方法中每種細(xì)胞亞群中上調(diào)屿附、下調(diào)和沒(méi)有統(tǒng)計(jì)學(xué)差異的基因集數(shù)目；上方的條形代表每個(gè)亞群中不同方法中差異的基因數(shù)目砾医，紅色代表上調(diào)的差異基因集拿撩，藍(lán)色代表下調(diào)的差異基因集衣厘；中間的柱形圖代表每個(gè)亞群中不同方法中上調(diào)如蚜、下調(diào)和沒(méi)有統(tǒng)計(jì)學(xué)意義的基因集的比例压恒；

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

圖片

2）. 局部展示

①.密度散點(diǎn)圖

密度散點(diǎn)圖將基因集的富集分?jǐn)?shù)和細(xì)胞亞群在低維空間的投影結(jié)合起來(lái)，展示了特定基因集在空間上的表達(dá)水平错邦。其中探赫，顏色越黃，代表富集分?jǐn)?shù)越高撬呢；

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

圖片

②.半小提琴圖

半小提琴圖同時(shí)以小提琴圖（左邊）和箱線圖（右邊）進(jìn)行展示伦吠。不同顏色代表不同的細(xì)胞亞群；

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

圖片

③.山巒圖

山巒圖中上方的核密度曲線展示了數(shù)據(jù)的主要分布魂拦，下方的條形編碼圖展示了細(xì)胞亞群具體的數(shù)量毛仪。不同顏色代表不同的細(xì)胞亞群，而橫坐標(biāo)代表不同的表達(dá)水平芯勘；

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

圖片

④.密度熱圖

密度熱圖展示了具體差異基因在不同細(xì)胞亞群中的表達(dá)和分布水平箱靴。顏色越紅，代表富集分?jǐn)?shù)越高荷愕；

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

圖片

生活很好衡怀，有你更好