Escape v1.5.1套件-單細胞基因集通路分析(ssGSEA)範例

前言

Gene Set Enrichment Analysis (GSEA) 是一種常用於功能基因組學(Functional genomics)研究的分析方法。它的目的是評估在特定條件下,例如正常組和病變組樣本,特定基因集(gene set)是否在基因表達數據中具有統計上的顯著差異。這種方法可以幫助分析人員了解哪些可能的通路(pathway),會因為特定條件而可能被啟動或抑制。而ssGSEA(Single-sample Gene Set Enrichment Analysis)是GSEA方法的擴展,與傳統的GSEA方法相比,ssGSEA是針對個別樣本進行分析,而不是對比兩個不同條件之間的差異,主要用於評估在單個樣本中特定基因集的活性程度,因此適合用於單細胞數據的基因集通路分析。

Escape套件是一個分析ssGSEA的R package,由Nicholas Borcherding團隊所開發,可對樣本進行通路分析、可視化以及顯著分析,基本上只要一個package就能將ssGSEA分析從頭做到尾,不需要再去考慮如何解決下游的結果呈現等問題,推薦給大家使用。

下載以及安裝Escape套件

首先,你必須在電腦上安裝R語言和RStudio,然後打開RStudio視窗後輸入安裝Escape套件的指令:

# 下載devtools套件
install.packages("devtools") 
# 下載Escape套件
devtools::install_github("ncborcherding/escape")

範例數據

在這次的範例中,我們將分析來自10X Genomics的周邊血(PBMC)的兩組單細胞數據集,分別包含了2700和5025顆單細胞的數據集,原始數據的.rds儲存檔可以在我的雲端下載。如果對檔案有疑慮,可以到10X Genomics數據庫自行下載,然後按照單細胞數據分析:批次效應與Seurat Integration分析詳解這篇文章進行分析。

數據集下載位置:

3k PBMCs from a Healthy Donor

5k Peripheral Blood Mononuclear Cells (PBMCs) from a Healthy Donor (v3 chemistry)

1. 加載套件和載入分析好的PBMC數據集

加載所需套件進入R環境,並將pbmc3k5K_final.rds數據載入R中,這裡使用的指令是readRDS()函數,代碼中提供的路徑只是一個示例路徑,需要根據電腦實際文件所在的路徑進行修改

# 加載套件
library(Seurat)
library(dplyr)
library(patchwork)
library(escape)

# 需用電腦中的正確路徑,須將"\"更改為"/"
pbmc.combined <- readRDS(file = "C:/Users/Administrator/Desktop/pbmc3k5K_final.rds")

# 檢查輸入結果
DimPlot(pbmc.combined, reduction = "umap", group.by = "orig.ident")
DimPlot(pbmc.combined, reduction = "umap")
DimPlot(pbmc.combined, reduction = "umap",label = TRUE, split.by = "orig.ident")

2. 選擇想要分析的基因集(gene set)

Molecular Signature Database(MSigDB)是一個收錄了數以萬計基因集的資料庫,目前有人類和小鼠兩種資料庫,可以用於GO、KEGG、GSEA、ssGSEA等分析。

escape2

人類基因集資料庫共分為九個子資料庫,分別為:

  • H:hallmark gene sets (癌症)特徵基因集合,共有50組基因集,最常使用。
  • C1:positional gene sets,共有300組基因集。
  • C2:curated gene sets,共有6495組基因集。
  • C3:regulatory target gene sets,共有3713組基因集。
  • C4:computational gene sets,共有858組基因集。
  • C5:ontology gene sets,包括biological process(BP),cellular component(CC)以及molecular function(MF)三個部分,共有15937組基因集,最常使用。
  • C6:oncogenic signatures gene sets,共有189組基因集。
  • C7:immunologic signatures gene sets,共有5219組基因集,最常使用。
  • C8:cell type signature gene sets,共有830組基因集。

Escape套件很貼心地將這個資料庫歸納進套件內,所以我們只要使用getGeneSets()函數,便可以從MSigDB資料庫中選擇我們想要進行分析的基因集,通過設定library()函數可以指定想要的基因集子資料庫,指令如下:

# 從MSigDB資料庫提取想要分析的基因集,以hallmark gene sets為例子
GS <- getGeneSets(library = "H")

3. 進行ssGSEA分析

# 進行ssGSEA分析
ES <- enrichIt(obj = pbmc.combined, gene.sets = GS, groups = 1000, cores = 16)

# 將分析結果加入pbmc Seurat 資料中
pbmc.combined <- AddMetaData(pbmc.combined, ES)

4. 可視化結果(Visualizations)

接著便是將分析好的結果進行可視化,需要先下載dittoSeq套件,方法如下:

# 下載dittoSeq套件
BiocManager::install("dittoSeq")

4.1 熱圖 (Heatmap)

# 載入dittoSeq套件
library(dittoSeq)
# 載入調色板
colorblind_vector <- colorRampPalette(c("#FF4B20", "#FFB433", "#C6FDEC", "#7AC5FF", "#0348A6"))
# 熱圖(Heatmap)
pbmc.combined@meta.data$active.idents <- pbmc.combined@active.ident
dittoHeatmap(pbmc.combined, genes = NULL, metas = names(ES),
             heatmap.colors = rev(colorblind_vector(50)),
             annot.by = c("active.idents", "orig.ident"),
             fontsize = 7)
Escape-1

如果不喜歡調色板的顏色,可以選擇不加或自行更換調色板的顏色:

dittoHeatmap(pbmc.combined, genes = NULL, metas = names(ES),
             annot.by = c("active.idents", "orig.ident"),
             fontsize = 7)
Escape-3

可自行挑選需要的基因集:

dittoHeatmap(pbmc.combined, genes = NULL, 
             metas = c("HALLMARK_APOPTOSIS", "HALLMARK_DNA_REPAIR", "HALLMARK_P53_PATHWAY"), 
             heatmap.colors = rev(colorblind_vector(50)),
             annot.by = c("active.idents", "orig.ident"),
             fontsize = 7)
Escape-4

4.2 小提琴圖(Violin Plot)

dittoPlot(pbmc.combined, "HALLMARK_DNA_REPAIR", group.by = "orig.ident") + scale_fill_manual(values = colorblind_vector(2))
Escape-5

多個基因集通路小提琴圖指令:

multi_dittoPlot(pbmc.combined, vars = c("HALLMARK_WNT_BETA_CATENIN_SIGNALING", "HALLMARK_XENOBIOTIC_METABOLISM","HALLMARK_UV_RESPONSE_UP"), 
                group.by = "orig.ident", plots = c("jitter", "vlnplot", "boxplot"), 
                ylab = "Enrichment Scores", 
                theme = theme_classic() + theme(plot.title = element_text(size = 10)))
Escape-6

4.3 密集富集圖(Hex Density Enrichment Plots)

如果認為兩個通路之間具有一定的相關性,可以使用密集富集圖來做進一步的觀察:

# 下載所需套件
BiocManager::install("hexbin")
# 密集富集圖可視化
dittoScatterHex(pbmc.combined, x.var = "HALLMARK_DNA_REPAIR", y.var = "HALLMARK_MTORC1_SIGNALING", do.contour = TRUE) +             
  theme_classic() + scale_fill_gradientn(colors = rev(colorblind_vector(11))) +
  geom_vline(xintercept = 0, lty=2) + geom_hline(yintercept = 0, lty=2) 
Escape-7

可以將pbmc3k和pbmc5K的結果以split.by()函數分開:

dittoScatterHex(pbmc.combined,
                x.var = "HALLMARK_DNA_REPAIR", y.var = "HALLMARK_MTORC1_SIGNALING", do.contour = TRUE, split.by = "orig.ident") + 
  theme_classic() +
  scale_fill_gradientn(colors = rev(colorblind_vector(11))) +
  geom_vline(xintercept = 0, lty=2) + 
  geom_hline(yintercept = 0, lty=2) 
Escape-8

4.4  山峰圖(Ridge Plot)

ES2 <- data.frame(pbmc.combined[[]], Idents(pbmc.combined))
colnames(ES2)[ncol(ES2)] <- "cluster"
ridgeEnrichment(ES2, gene.set = "HALLMARK_DNA_REPAIR", group = "orig.ident", add.rug = TRUE)
Escape-9
ridgeEnrichment(ES2, gene.set = "HALLMARK_DNA_REPAIR", group = "seurat_clusters", facet = "orig.ident", add.rug = TRUE)
Escape-10

4.5 分組小提琴圖(Split Violin Plot)

splitEnrichment(ES2, split = "orig.ident", gene.set = "HALLMARK_DNA_REPAIR")
splitEnrichment(ES2, x.axis = "cluster", split = "orig.ident", gene.set = "HALLMARK_DNA_REPAIR")
Escape-11

5. 擴展分析(Expanded Analysis)- PCA分析

使用performPCA()函數在基因集富集度分數上執行PCA分析,可以幫助分析人員理解基因集富集度分析結果中不同基因集之間的變異性和關聯性。

# 執行PCA分析
PCA <- performPCA(enriched = ES2, groups = c("orig.ident", "cluster"))
pcaEnrichment(PCA, PCx = "PC1", PCy = "PC2", contours = FALSE, facet = "cluster") 
masterPCAPlot(ES2, gene.sets = names(GS), PCx = "PC1", PCy = "PC2", top.contribution = 10)
Escape-12

6. 顯著性分析(Significance Testing)

Escape套件還提供了顯著性分析,一共有五種方法,分別為T.test、linear regression(LR)、 Wilcoxon Rank Sum Test(Wilcoxon)、ANOVA 以及 Kruskal-Wallis(KW),分析人員可以依照需求選擇適合的分析方法,這裡以T.test為例子:

# 顯著性分析
output <- getSignificance(ES2, group = "orig.ident", fit = "T.test") 
View(output)
Escape-13

結語

Escape套件是我非常常用的一個分析ssGSEA的R套件,只要跟著教學步驟做,基本上就能得到需要的結果,並且還附帶了可視化方式,所以我個人滿推薦這個套件的。得到ssGSEA的結果後,接下來就是要針對結果做生物學上的解釋了,沒有其它捷徑,靠的就是各位強大的生物學背景以及知識!

參考文獻

Borcherding, N., Vishwakarma, A., Voigt, A.P. et al. Mapping the immune environment in clear cell renal carcinoma by single-cell genomics. Commun Biol 4, 122 (2021).

您的分享,我非常感激!!!
MillionQuesn
MillionQuesn

一個旅居台灣的異鄉人,分享突然靈光一現的亮點。

文章: 46

發佈留言

發佈留言必須填寫的電子郵件地址不會公開。 必填欄位標示為 *

error: Alert: Content selection is disabled!!