單細(xì)胞富集分析系列：
1. 背景
-
Gene Set Enrichment Analysis (GSEA)用來評估一個(gè)預(yù)先定義的基因集的基因在與表型相關(guān)度排序的基因表中的分布趨勢那伐，從而判斷其對表型的貢獻(xiàn)瑰妄。其輸入數(shù)據(jù)包含兩部分杯巨，一是已知功能的基因集 (可以是GO注釋泣矛、MsigDB的注釋或其它符合格式的基因集定義)，一是表達(dá)矩陣(也可以是排序好的列表)齿风，軟件會(huì)對基因根據(jù)其于表型的關(guān)聯(lián)度(可以理解為表達(dá)值的變化)從大到小排序药薯，然后判斷基因集內(nèi)每條注釋下的基因是否富集于表型相關(guān)度排序后基因表的上部或下部，從而判斷此基因集內(nèi)基因的協(xié)同變化對表型變化的影響救斑。 -
Gene Set Variation Analysis(GSVA)被稱為基因集變異分析童本，是一種非參數(shù)的無監(jiān)督分析方法，主要用來評估芯片和轉(zhuǎn)錄組的基因集富集結(jié)果脸候。主要是通過將基因在不同樣品間的表達(dá)量矩陣轉(zhuǎn)化成基因集在樣品間的表達(dá)量矩陣穷娱，從而來評估不同的代謝通路在不同樣品間是否富集绑蔫。GSVA和GSEA兩者的輸入數(shù)據(jù)相似，包括已知功能的基因集和表達(dá)矩陣 鄙煤。 - GSEA和GSVA與GO等富集分析的不同：
一般的富集分析（GO和KEGG）僅考慮差異基因晾匠，沒有考慮基因表達(dá)量 茶袒，這容易遺漏部分差異表達(dá)不顯著卻有重要生物學(xué)意義的基因梯刚，忽略一些基因的生物特性、基因調(diào)控網(wǎng)絡(luò)之間的關(guān)系及基因功能和意義等有價(jià)值的信息薪寓。一刀切的閾值亡资，對于發(fā)現(xiàn)真正的生物學(xué)效應(yīng)，許多時(shí)候是一種障礙向叉。而GSEA與GSVA使用的不是差異基因集而是全部基因锥腻。從基因集的富集角度出發(fā)，理論上更容易發(fā)現(xiàn)細(xì)微變化對生物通路的影響母谎，尤其是差異倍數(shù)不太大的基因集瘦黑。
參考：基因功能富集方法和基因注釋數(shù)據(jù)庫介紹
2. 原理
2.1 GSEA
(1) 背景基因排序：將全部基因按照某種指標(biāo)（差異分析p值，表型相關(guān)性奇唤，表達(dá)量等）進(jìn)行排序幸斥，比如log2FC排序。
(2) 目標(biāo)基因富集：將某個(gè)特定類型的基因在排序表中標(biāo)出咬扇，目標(biāo)基因可以是某個(gè)通路或GO terms的基因等甲葬。
(3) 計(jì)算富集分?jǐn)?shù)：使用加權(quán)法，計(jì)算ES值變化懈贺。對位于中部（與性狀相關(guān)性低）的部分采用較小的權(quán)值经窖，所以越集中在兩端，與表型的相關(guān)性越高梭灿。ES曲線最大值為富集分?jǐn)?shù)（Enrichment Score）画侣。
(4) Permutation test：對基因集的ES值進(jìn)行顯著性檢驗(yàn)及多重假設(shè)檢驗(yàn)，從而計(jì)算出顯著富集的基因集堡妒。
2.2 GSVA
GSEA雖然是一種強(qiáng)大的富集分析工具配乱，但是它的應(yīng)用場景通常局限于Case/Control的實(shí)驗(yàn)設(shè)計(jì)。對于表型（分組）復(fù)雜的大樣本量研究涕蚤，比如scRNAseq和TCGA這樣的項(xiàng)目宪卿，分析起來就顯得困難。因此万栅，Broad研究所又開發(fā)了GSVA算法來拓展基因集分析的應(yīng)用佑钾。GSVA不需要預(yù)先進(jìn)行樣本之間的差異分析，它依據(jù)表達(dá)矩陣就可以計(jì)算每個(gè)樣本中特定基因集的變異分?jǐn)?shù)烦粒。
在輸入后基因表達(dá)矩陣和特定基因集以后：
(1) 算法會(huì)對表達(dá)數(shù)據(jù)進(jìn)行核密度估計(jì)休溶；
(2) 基于第一步的結(jié)果對樣本進(jìn)行表達(dá)水平排序代赁；
(3) 對于每一個(gè)基因集進(jìn)行類似K-S檢驗(yàn)的秩統(tǒng)計(jì)量計(jì)算；
(4) 獲取GSVA富集分?jǐn)?shù)兽掰。最終輸出為以每個(gè)基因集對應(yīng)每個(gè)樣本的數(shù)據(jù)矩陣芭碍。
簡單的說，輸入以基因?yàn)樾械谋磉_(dá)矩陣和基因集數(shù)據(jù)庫給GSVA孽尽，它就輸出以基因集名稱為行的變異分?jǐn)?shù)矩陣窖壕。
2.3 總結(jié)：
GSEA和GSVA都是基于對基因的某一個(gè)值的排序來進(jìn)行富集分析速勇。
GSEA主要是用case和control之間的差異倍數(shù)或信噪比來進(jìn)行排序（一次處理兩個(gè)樣本）。GSVA則不需要做對比坎拐，而是對每個(gè)樣本或單個(gè)細(xì)胞按基因的表達(dá)量進(jìn)行單獨(dú)排序烦磁，然后將富集分?jǐn)?shù)的值做個(gè)標(biāo)準(zhǔn)化，再在不同的樣本間對比哼勇。
3. 分析
library(Seurat)
library(msigdbr)
library(GSVA)
library(tidyverse)
library(clusterProfiler)
library(patchwork)
rm(list=ls())
setwd("Function")
dir.create("GSEA")
dir.create("GSVA")
3.1 GSEA
- 讀入注釋好的數(shù)據(jù)都伪，進(jìn)行差異分析，得到差異表達(dá)矩陣
scRNA <- readRDS("scRNA.classified.rds")
# 樣本分組
scRNA$group <- recode(scRNA$orig.ident,
"HNC01PBMC" = "hnc_pbmc",
"HNC01TIL" = "hnc_til",
"HNC10PBMC" = "hnc_pbmc",
"HNC10TIL" = "hnc_til",
"HNC20PBMC" = "hnc_pbmc",
"HNC20TIL" = "hnc_til",
"PBMC1" = "PBMC",
"PBMC2" = "PBMC",
"Tonsil1" = "Tonsil",
"Tonsil2" = "Tonsil")
scRNA$celltype <- scRNA$SingleR
# 按細(xì)胞類型將對象分割
scRNA.list <- SplitObject(scRNA, split.by = "celltype")
# 數(shù)據(jù)標(biāo)準(zhǔn)化
scRNA.list <- lapply(scRNA.list, function(sco){
DefaultAssay(sco) <- "RNA"
sco <- NormalizeData(sco)
return(sco)
})
# 檢查數(shù)據(jù)
for(i in names(scRNA.list)){
print(i)
print(table(scRNA.list[[i]][["group"]]))
}
# 剔除NKT細(xì)胞（數(shù)目太少）
scRNA.list <- scRNA.list[1:6]
# 按細(xì)胞類型配對差異分析
library(future)
options(future.globals.maxSize = 20 * 1024^3)
plan(multisession, workers = 8)
deg.list <- lapply(scRNA.list, FUN = function(x) {
FindMarkers(x, ident.1 = "hnc_til", ident.2 = "hnc_pbmc", assay = "RNA",
group.by = "group", logfc.threshold = 0, min.pct = 0)})
saveRDS(deg.list, file = "deg.list.rds")
View(deg.list)
- GSEA分析（KEGG數(shù)據(jù)集）
此前寫過選取單個(gè)cluster的差異基因進(jìn)行GSEA分析猴蹂，見：單細(xì)胞GSEA分析流程院溺，在這里演示批量操作。
filenames = c("CD8_TCs","DCs","Monocytes", "CD4_TCs", "NKs", "BCs")
#??filenames順序一定要和deg.list中順序一致磅轻，否則最后導(dǎo)出文件會(huì) 出錯(cuò)珍逸。
dir.create("GSEA/kegg")
setwd("./GSEA/kegg")
# 設(shè)置基因集
genesets = msigdbr(species = "Homo sapiens", category = "C2") #msigdbr提供多個(gè)物種的基因集數(shù)據(jù)
# View(msigdbr_collections()) #查看msigdbr包中所有的基因集
unique(genesets$gs_subcat) # 有多個(gè)數(shù)據(jù)庫來源的基因集可選，這里選用KEGG
genesets <- subset(genesets, gs_subcat=="CP:KEGG", select = c("gs_name", "gene_symbol"))
unique(genesets$gs_name) #查看有多少條通路（186個(gè)）
# GSEA分析
deg.list <- readRDS("deg.list.rds")
res.list <- lapply(deg.list, FUN = function(x){
x = x[order(x$avg_log2FC, decreasing = T),]
genelist <- structure(x$avg_log2FC, names = rownames(x))
res <- GSEA(genelist, TERM2GENE = genesets, eps = 0)
return(res) #按差異倍數(shù)進(jìn)行排序并返回GSEA結(jié)果
})
# 導(dǎo)出表格
for(i in seq_along(res.list)){
res <- data.frame(res.list[[i]])
write.csv(res, paste0(filenames[i], ".csv"), row.names = F)
}
# 導(dǎo)出圖形
for(i in seq_along(res.list)){
res <- res.list[[i]]
for(j in seq_along(res@result$ID)){
p <- gseaplot(res, geneSetID = j, title = res@result$ID[j], by = "runningScore")
filename <- paste0(filenames[i], "_", res@result$ID[j], '.pdf')
ggsave(filename = filename, p, width = 8, height = 4)
}
} #每個(gè)細(xì)胞類型針對每個(gè)通路注盈，都會(huì)出一個(gè)圖
setwd("~/project/Function")
3.2 分組平均水平GSVA
單細(xì)胞表達(dá)矩陣上有一萬多個(gè)基因，但實(shí)際上在單個(gè)細(xì)胞水平上叙赚，很多基因的表達(dá)都是0老客，因此在做單個(gè)細(xì)胞水平上的GSEA僚饭，排序結(jié)果會(huì)有問題。但是按分組排序的話胧砰，基本上每個(gè)基因都有一個(gè)表達(dá)量的值（非零）鳍鸵，非零的基因會(huì)比單細(xì)胞矩陣要多得多，得出的排序結(jié)果會(huì)更準(zhǔn)尉间。
- 數(shù)據(jù)準(zhǔn)備
rm(list = ls())
scRNA <- readRDS("scRNA.classified.rds")
DefaultAssay(scRNA) <- "RNA"
scRNA <- NormalizeData(scRNA)
- 基因集準(zhǔn)備（Hallmark基因集）
setwd("Function")
dir.create("GSVA/avg_hallmark")
setwd("./Function/GSVA/avg_hallmark")
# 選擇基因集
genesets <- msigdbr(species = "Homo sapiens", category = "H")
genesets <- subset(genesets, select = c("gs_name","gene_symbol")) %>% as.data.frame()
genesets <- split(genesets$gene_symbol, genesets$gs_name)
- 提取分組平均表達(dá)矩陣
#Idents設(shè)置按什么來取組的表達(dá)均值（計(jì)算case control之間的均值也可以）
Idents(scRNA) <- "SingleR"
expr <- AverageExpression(scRNA, assays = "RNA", slot = "data")[[1]]
expr <- expr[rowSums(expr)>0,] #選取非零基因
expr <- as.matrix(expr)
- GSVA富集分析
# gsva默認(rèn)開啟全部線程計(jì)算
gsva.res <- gsva(expr, genesets, method="ssgsea")
saveRDS(gsva.res, "gsva.res.rds")
gsva.df <- data.frame(Genesets=rownames(gsva.res), gsva.res, check.names = F)
write.csv(gsva.df, "gsva_res.csv", row.names = F)
pheatmap::pheatmap(gsva.res, show_colnames = T, scale = "row")
另：如果做KEGG基因集的變異分析偿乖，和Hallmark基因集一樣，只是在做基因集準(zhǔn)備的時(shí)候換成下面的代碼即可乌妒。 ??注意切換路徑
genesets <- msigdbr(species = "Homo sapiens", category = "C2")
genesets <- subset(genesets, gs_subcat=="CP:KEGG", select = c("gs_name", "gene_symbol")) %>% as.data.frame()
genesets <- split(genesets$gene_symbol, genesets$gs_name)
3.3 單細(xì)胞水平GSVA
單細(xì)胞水平GSVA汹想，在細(xì)胞數(shù)比較多（比如說幾萬個(gè)細(xì)胞）的情況下，再選用pathway比較多的數(shù)據(jù)集比如KEGG撤蚊，運(yùn)行會(huì)很慢，可能需要數(shù)小時(shí)到1-2天损话。
- 數(shù)據(jù)準(zhǔn)備
rm(list = ls())
scRNA <- readRDS("scRNA.classified.rds")
DefaultAssay(scRNA) <- "RNA"
scRNA <- NormalizeData(scRNA)
- 基因集準(zhǔn)備（KEGG基因集）
setwd("Function")
dir.create("GSVA/kegg")
setwd("./GSVA/kegg")
# 選擇基因集
genesets <- msigdbr(species = "Homo sapiens", category = "C2")
genesets <- subset(genesets, gs_subcat=="CP:KEGG", select = c("gs_name", "gene_symbol")) %>% as.data.frame()
genesets <- split(genesets$gene_symbol, genesets$gs_name)
- 提取單細(xì)胞表達(dá)矩陣
expr <- as.matrix(scRNA@assays$RNA@data)
#set.seed(seed = 123)
#expr <- expr[,sample(colnames(expr), 1000)]
expr <- expr[rowSums(expr)>0,] #選取非零基因
- GSVA富集分析
# gsva默認(rèn)開啟全部線程計(jì)算
gsva.res <- gsva(expr, genesets, method="ssgsea", parallel.sz=8)
saveRDS(gsva.res, "gsva.res.rds")
gsva.df <- data.frame(Genesets=rownames(gsva.res), gsva.res, check.names = F)
write.csv(gsva.df, "gsva_res.csv", row.names = F)
- 指定基因集展示
??注意基因集的名字中把“_”改為“-”侦啸，不然會(huì)顯示找不到pathway
scRNA[["gsva"]] <- CreateAssayObject(gsva.res)
DefaultAssay(scRNA) <- "gsva"
VlnPlot(scRNA, features = "KEGG-ABC-TRANSPORTERS",
group.by = "SingleR", slot = "counts", pt.size = 0)
pheatmap::pheatmap(gsva.res, show_colnames = F, scale = "row")
