单细胞之富集分析-2：批量GSEA和GSVA分析

2021-05-10 本文已影响0人 Hayley笔记

单细胞富集分析系列：

单细胞之富集分析-1：单细胞GSEA分析流程

1. 背景

Gene Set Enrichment Analysis (GSEA)用来评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋势，从而判断其对表型的贡献。其输入数据包含两部分，一是已知功能的基因集 (可以是GO注释、MsigDB的注释或其它符合格式的基因集定义)，一是表达矩阵(也可以是排序好的列表)，软件会对基因根据其于表型的关联度(可以理解为表达值的变化)从大到小排序，然后判断基因集内每条注释下的基因是否富集于表型相关度排序后基因表的上部或下部，从而判断此基因集内基因的协同变化对表型变化的影响。
Gene Set Variation Analysis(GSVA)被称为基因集变异分析，是一种非参数的无监督分析方法，主要用来评估芯片和转录组的基因集富集结果。主要是通过将基因在不同样品间的表达量矩阵转化成基因集在样品间的表达量矩阵，从而来评估不同的代谢通路在不同样品间是否富集。GSVA和GSEA两者的输入数据相似，包括已知功能的基因集和表达矩阵。
GSEA和GSVA与GO等富集分析的不同：
一般的富集分析（GO和KEGG）仅考虑差异基因，没有考虑基因表达量，这容易遗漏部分差异表达不显著却有重要生物学意义的基因，忽略一些基因的生物特性、基因调控网络之间的关系及基因功能和意义等有价值的信息。一刀切的阈值，对于发现真正的生物学效应，许多时候是一种障碍。而GSEA与GSVA使用的不是差异基因集而是全部基因。从基因集的富集角度出发，理论上更容易发现细微变化对生物通路的影响，尤其是差异倍数不太大的基因集。
参考：基因功能富集方法和基因注释数据库介绍

2. 原理

2.1 GSEA

(1) 背景基因排序：将全部基因按照某种指标（差异分析p值，表型相关性，表达量等）进行排序，比如log2FC排序。
(2) 目标基因富集：将某个特定类型的基因在排序表中标出，目标基因可以是某个通路或GO terms的基因等。
(3) 计算富集分数：使用加权法，计算ES值变化。对位于中部（与性状相关性低）的部分采用较小的权值，所以越集中在两端，与表型的相关性越高。ES曲线最大值为富集分数（Enrichment Score）。
(4) Permutation test：对基因集的ES值进行显著性检验及多重假设检验，从而计算出显著富集的基因集。

2.2 GSVA

GSEA虽然是一种强大的富集分析工具，但是它的应用场景通常局限于Case/Control的实验设计。对于表型（分组）复杂的大样本量研究，比如scRNAseq和TCGA这样的项目，分析起来就显得困难。因此，Broad研究所又开发了GSVA算法来拓展基因集分析的应用。GSVA不需要预先进行样本之间的差异分析，它依据表达矩阵就可以计算每个样本中特定基因集的变异分数。
在输入后基因表达矩阵和特定基因集以后：
(1) 算法会对表达数据进行核密度估计；
(2) 基于第一步的结果对样本进行表达水平排序；
(3) 对于每一个基因集进行类似K-S检验的秩统计量计算；
(4) 获取GSVA富集分数。最终输出为以每个基因集对应每个样本的数据矩阵。
简单的说，输入以基因为行的表达矩阵和基因集数据库给GSVA，它就输出以基因集名称为行的变异分数矩阵。

左侧输入基因表达矩阵和基因集数据库，中间是GSVA算法原理，右侧是输出的基因集变异分数矩阵。基因集变异分数可以理解为基因集内所有基因的综合表达值。

2.3 总结：

GSEA和GSVA都是基于对基因的某一个值的排序来进行富集分析。
GSEA主要是用case和control之间的差异倍数或信噪比来进行排序（一次处理两个样本）。GSVA则不需要做对比，而是对每个样本或单个细胞按基因的表达量进行单独排序，然后将富集分数的值做个标准化，再在不同的样本间对比。

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

读入注释好的数据，进行差异分析，得到差异表达矩阵

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

# 按细胞类型将对象分割
scRNA.list <- SplitObject(scRNA, split.by = "celltype")
# 数据标准化
scRNA.list <- lapply(scRNA.list, function(sco){
  DefaultAssay(sco) <- "RNA"
  sco <- NormalizeData(sco)
  return(sco)
})
# 检查数据
for(i in names(scRNA.list)){
  print(i)
  print(table(scRNA.list[[i]][["group"]]))
}
# 剔除NKT细胞（数目太少）
scRNA.list <- scRNA.list[1:6]
# 按细胞类型配对差异分析
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)

使用FindMarkers对hnc_til和hnc_pbmc两组的6个细胞群的18818个基因做了差异分析

GSEA分析（KEGG数据集）
此前写过选取单个cluster的差异基因进行GSEA分析，见：单细胞GSEA分析流程，在这里演示批量操作。

filenames = c("CD8_TCs","DCs","Monocytes", "CD4_TCs", "NKs", "BCs")
#⚠️filenames顺序一定要和deg.list中顺序一致，否则最后导出文件会 出错。
dir.create("GSEA/kegg")
setwd("./GSEA/kegg")

# 设置基因集
genesets = msigdbr(species = "Homo sapiens", category = "C2") #msigdbr提供多个物种的基因集数据
# View(msigdbr_collections()) #查看msigdbr包中所有的基因集
unique(genesets$gs_subcat)  # 有多个数据库来源的基因集可选，这里选用KEGG
genesets <- subset(genesets, gs_subcat=="CP:KEGG", select = c("gs_name", "gene_symbol"))
unique(genesets$gs_name) #查看有多少条通路（186个）
# 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) #按差异倍数进行排序并返回GSEA结果
})

res.list的结构。输入进行比较的通路有186个，输出的每个组的通路各不相同（CD8T有29个，B只有1个）。

# 导出表格
for(i in seq_along(res.list)){
  res <- data.frame(res.list[[i]])
  write.csv(res, paste0(filenames[i], ".csv"), row.names = F)
}

运行完上面的代码，工作目录下每个细胞都会有一个GSEA差异分析矩阵。矩阵的横轴是通路，纵轴是富集分数、NES值、p值、q值以及富集到每条通路上的基因等

# 导出图形
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)
  }
} #每个细胞类型针对每个通路，都会出一个图
setwd("~/project/Function")

3.2 分组平均水平GSVA

单细胞表达矩阵上有一万多个基因，但实际上在单个细胞水平上，很多基因的表达都是0，因此在做单个细胞水平上的GSEA，排序结果会有问题。但是按分组排序的话，基本上每个基因都有一个表达量的值（非零），非零的基因会比单细胞矩阵要多得多，得出的排序结果会更准。

数据准备

rm(list = ls())
scRNA <- readRDS("scRNA.classified.rds")
DefaultAssay(scRNA) <- "RNA"
scRNA <- NormalizeData(scRNA)

基因集准备（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)

提取分组平均表达矩阵

#Idents设置按什么来取组的表达均值（计算case control之间的均值也可以）
Idents(scRNA) <- "SingleR" 
expr <- AverageExpression(scRNA, assays = "RNA", slot = "data")[[1]]
expr <- expr[rowSums(expr)>0,]  #选取非零基因
expr <- as.matrix(expr)

GSVA富集分析

# gsva默认开启全部线程计算
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)

gsva.df

pheatmap::pheatmap(gsva.res, show_colnames = T, scale = "row")

另：如果做KEGG基因集的变异分析，和Hallmark基因集一样，只是在做基因集准备的时候换成下面的代码即可。 ⚠️注意切换路径

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 单细胞水平GSVA

单细胞水平GSVA，在细胞数比较多（比如说几万个细胞）的情况下，再选用pathway比较多的数据集比如KEGG，运行会很慢，可能需要数小时到1-2天。

数据准备

rm(list = ls())
scRNA <- readRDS("scRNA.classified.rds")
DefaultAssay(scRNA) <- "RNA"
scRNA <- NormalizeData(scRNA)

基因集准备（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)

提取单细胞表达矩阵

expr <- as.matrix(scRNA@assays$RNA@data)
#set.seed(seed = 123)
#expr <- expr[,sample(colnames(expr), 1000)] 
expr <- expr[rowSums(expr)>0,]  #选取非零基因

GSVA富集分析

# gsva默认开启全部线程计算
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)

指定基因集展示
⚠️注意基因集的名字中把“_”改为“-”，不然会显示找不到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")