查看原文
其他

单细胞转录组高级分析五:GSEA与GSVA分析

Kinesin 生信会客厅 2022-06-07
上期专题我们介绍了单细胞转录组数据的基础分析,然而那些分析只是揭开了组织异质性的面纱,还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容:多样本批次校正、SCENIC转录调控网络分析、CellPhoneDB细胞通讯分析、inferCNV肿瘤拷贝数变异分析、GSEA/GSVA分析和其他高级分析内容。不足之处请大家批评指正,欢迎添加Kinesin微信交流探讨!(扫描文末二维码)

往期相关文章

单细胞转录组基础分析专题

单细胞转录组高级分析一:多样本合并与批次校正
单细胞转录组高级分析二:转录调控网络分析
单细胞转录组高级分析三:细胞通讯分析
单细胞转录组高级分析四:scRNA数据推断CNV

GSEA与GSVA简介
基因集的概念
GSEA全称Gene Set Enrichment Analysis,GSVA全称Gene Set Variation Analysis,它们都是基于基因集开展的分析,因此我们先要了解基因集的定义。基因集顾名思义就是一些基因的集合,任何一些基因放在一起都可以叫做基因集,但是我们用来分析的基因集要求有一定的生物学意义。最常见的基因集数据库如GO与KEGG,它们一个按照基因本体论将基因分门别类,一个按照代谢通路将相关基因集合在一起。除此之外我们还可以按转录因子调控网络、共表达网络、定义生物状态的marker基因列表等把基因集合成有一定生物学意义的基因集。

MSigDB基因集数据库

GSEA是由Broad研究所开发的一种富集方法,他们在提出该方法的同时还提供了一个基因集数据库——MSigdb。它从位置,功能,代谢途径,靶标结合等多种角度出发,构建出了许多的基因集合,Broad研究所将他们构建的基因集合保存在MSigDB,官网地址如下:
http://software.broadinstitute.org/gsea/msigdb/index.jsp

在MSigDB中,将所有的基因集划分为以下8大类别:

1. H:hallmark gene sets

特征基因集,由定义生物状态和进程的marker基因组成。

2. C1:positional gene sets

位置基因集,包含人类每条染色体上的不同cytoband区域对应的基因集合。

3. C2:curated gene sets

代谢通路基因集,包含KEGG, Reactome, BioCarta数据库,以及文献和专家支持的基因集信息。

4. C3:motif gene sets

靶基因集,包含了miRNA靶基因集和转录因子调控基因集两大类。

5. C4:computational gene sets

计算基因集,计算机软件预测出来的基因集,主要是和癌症相关的基因。

6. C5:GO gene sets

基因本体基因集,包含了Gene Ontology对应的基因集合。

7. C6:oncogenic signatures

癌症扰动基因集,来源于药物处理肿瘤后基因差异表达数据,包含已知条件处理后基因表达量发生变化的基因。

8. C7:immunologic signatures

免疫基因集,包含了免疫系统功能相关的基因集合。


GSEA的分析原理
常规GO/KEGG富集分析需要设定阈值过滤差异基因,阈值太宽富集的结果太多,阈值太严又可能会遗漏一些关键结果。GO/KEGG富集的结果通常还很宽泛,并不能很好地解释生物学现象。有鉴于此,Broad研究所开发了基因集富集分析(GSEA)方法。GSEA使用无监督算法,不用过滤任何基因,配合MSigDB数据库使用,更容易找到解释生物学现象的基因集。其原理如下:

GSEA分析要先将样本做组间对比分析,GSEA自带9种分析方法,分为基因表达值差异分析和相关性分析两大类。对于Case/Contral的实验设计,差异分析方法更为常用,这其中又以默认的信噪比和大家熟悉的差异倍数用的最多。对比分析之后要按结果将基因排序,以差异倍数方法为例,把所有基因按差异倍数(FC)的值降序排列以供后续分析。上图小人脚下的小方块代表排序好的差异基因列表,蓝色之外的其他色块代表属于某个基因集的基因,如黄色属于基因集A,绿色属于基因集B。最下面高低不等的竖条代表与基因列表对应的FC值,红色上调、蓝色下调、黄色没有变化。基因集的富集分析需要经历三步:
  1. 基因集A富集分析时,小人从基因列表的左端走到右端,每经过一个蓝色基因扣分,每遇到一个黄色基因加分,扣分时与FC无关,加分时考虑FC的权重。基因集A最终的富集分数(ES)是小人曾经得过的最高/低分,实际公式比这复杂,但基本理念如此。

  2. 采用置换检验计算基因集A的显著性,即p值。

  3. 基因集A富集分析完成后,按上述同样的方法完成基因集B、C直至所有输入基因集的分析。所有需要富集分析的基因集都计算ES和p值之后,将ES转换为标准富集分数(NES),并计算校正后p值。


听完我的解释之后再看官方的解释可能更容易理解:

A GSEA overview illustrating the method. (A) An expression dataset sorted by correlation with phenotype, the corresponding heat map, and the ‘‘gene tags,’’ i.e., location of genes from a set S within the sorted list. (B) Plot of the running sum for S in the dataset, including the location of the maximum enrichment score (ES) and the leading-edge subset.


GSVA的分析原理
GSEA虽然是一种强大的富集分析工具,但是它的应用场景通常局限于Case/Control的实验设计。对于表型(分组)复杂的大样本量研究,例如TCGA和单细胞图谱这样的项目,分析起来就非常困难。因此,Broad研究所在GSEA发布8年之后,又开发了GSVA算法来拓展基因集分析的应用。GSVA不需要预先进行样本之间的差异分析,它依据表达矩阵就可以计算每个样本中特定基因集的变异分数。简单的说,输入以基因为行的表达矩阵和基因集数据库给GSVA,它就输出以基因集名称为行的变异分数矩阵,如下图所示:

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

GSEA的安装与使用
GSEA官方提供很多版本,但是不提倡R语言版,下载网址如下:
https://www.gsea-msigdb.org/gsea/downloads.jsp


准备分析文件
单细胞GSEA分析的难点是输入文件的准备,大家在R语言中运行下列代码,就可以得到符合GSEA要求的特定格式文件。
##R语言准备gsea输入文件library(Seurat)library(tidyverse)library(GSEA)dir.create("GSEA")   dir.create("GSEA/input")dir.create("GSEA/output")scRNA <- readRDS("scRNA.rds")#提取HNC01TIL和Tonsil2的T细胞做富集分析tmp <- subset(tmp, subset = (tmp$orig.ident=='HNC01TIL'|tmp$orig.ident=='Tonsil2')&tmp$celltype_Monaco=='T cells')sub.cells <- rownames(tmp)scRNAsub <- subset(scRNA, cells=sub.cells)expr <- GetAssayData(scRNAsub, slot = 'counts')expr <- data.frame(NAME=rownames(expr), Description=rep('na', nrow(expr)), expr, stringsAsFactors=F)write('#1.2', "GSEA/input/expr.gct", ncolumns=1)write(c(nrow(expr),(ncol(expr)-2)), "GSEA/input/expr.gct", ncolumns=2, append=T, sep='\t')write.table(expr, "GSEA/input/expr.gct", row.names=F, sep='\t', append=T, quote=F)line.1 <- c((ncol(expr)-2), 2, 1)tmp <- table(scRNAsub@meta.data$orig.ident)line.2 <- c("#", names(tmp))line.3 <- c(rep(names(tmp)[1],tmp[1]), rep(names(tmp)[2],tmp[2]))write(line.1, 'GSEA/input/group.cls', ncolumns=length(line.1), append=T, sep='\t')write(line.2, 'GSEA/input/group.cls', ncolumns=length(line.2), append=T, sep='\t')write(line.3, 'GSEA/input/group.cls', ncolumns=length(line.3), append=T, sep='\t')
准备好的文件格式:
表达数据gct文件内容节选

分组cls文件内容节选

得到上面的gct和cls文件,再加上基因集数据库gmt文件,就可以进行GSEA分析了。建议在java版本上运行,官方的R包太简陋了,java版上的好多功能都不支持。


R语言版代码与结果展示

##安装官方R包GSEAlibrary(devtools)install_github("GSEA-MSigDB/GSEA_R")#gs.db参数对应的是基因集数据文件,这里只对kegg收录的基因集做了分析GSEA::GSEA("GSEA/input/expr.gct", "GSEA/input/group.cls",            gs.db="~/database/GSEA/c2.cp.kegg.v7.1.symbols.gmt", output.directory="GSEA/output/")

分析得到的结果很多,分为结果总览,分组概况,以及各个基因集的富集情况详细。这里展示NES分数最高的通路的富集结果:

上图对应数据


java版分析代码

提前安装好java11,下载图中所示GSEA版本,解压即可
#以下命令在linux环境中运行,windows环境用gsea-cli.bat代替gsea-cli.shbash gsea-cli.sh GSEA -res ~/project/2020/2007_10xDemo2/GSEA/input/expr.gct \ -cls ~/project/2020/2007_10xDemo2/GSEA/input/group.cls#HNC01TIL_versus_Tonsil2 \ -gmx ~/project/2020/2007_10xDemo2/GSEA/input/c2.cp.kegg.v7.1.symbols.gmt \ -collapse false \ -out ~/project/2020/2007_10xDemo2/GSEA/out/

图形界面版GSEA操作
表达数据和分组文件准备好之后,图形界面版操作没有难度,简单介绍一下:

这是开始界面,先加载数据

点击箭头所指按钮之后,就会弹出对话框让你选择文件,分次输入我们之前准备好的表达数据文件(*.gct)和分组文件(*.cls)

选好文件之后,文件目录会出现在左右两边箭头所指的地方。如果右边没有,双击左边即可。我这儿左边有三个文件是因为加载了本地保存的基因集数据库文件。

文件加载完毕之后,点击箭头所指运行会出现右边的对话框,注意一下红框标注的选项。Collapse/Remap项选择“No_Collapse”(Collapse适用芯片数据),Permutation项选择“phenotype”(更适合单细胞数据),Chip项可不选或选择有“Human_Symbol_with_”关键字的项目(选择之后结果中基因有简短注释)。

Basic fields点开之后选择保存结果的目录,点击“Run”就开始运行了。java版(图形界面的也是java版)运行之后结果与R语言版一致,不过分析报告是网页版,点击超链接就可以查看细节,使用起来非常方便。

图形也更美观


GSVA的安装与使用
安装GSVA
#从bioconductor安装if (!requireNamespace("BiocManager", quietly = TRUE))    install.packages("BiocManager")BiocManager::install("GSVA")#或者从github安装library(devtools)install_github("rcastelo/GSVA")
准备输入数据
GSVA要求输入表达矩阵和基因集列表。表达矩阵从seurat对象导入即可,以列表的形式提供基因集数据有一定难度,我编写了一个函数解决此问题。
library(Seurat)library(GSVA)library(tidyverse)##创建gmt文件转list函数gmt2list <- function(gmtfile){ sets <- as.list(read_lines(gmtfile)) for(i in 1:length(sets)){ tmp = str_split(sets[[i]], '\t') n = length(tmp[[1]]) names(sets)[i] = tmp[[1]][1] sets[[i]] = tmp[[1]][3:n] rm(tmp, n) } return(sets)}#读取基因集数据库s.sets = gmt2list("GSEA/input/c2.cp.kegg.v7.1.symbols.gmt")#读取表达矩阵scRNA <- readRDS("scRNA.rds")# 随机提取1000个细胞演示GSVA,非常规操作# tmp <- sample(colnames(scRNA),1000) %>% sort()# scRNA <- scRNA[,tmp]expr <- as.matrix(scRNA@assays$RNA@counts)meta <- scRNA@meta.data[,c("seurat_clusters", "celltype_Monaco")]es.matrix = gsva(expr, s.sets, kcdf="Poisson")write.table(es.matrix, 'GSVA/gsva.xls', row.names=T, col.names=NA, sep='\t')

结果可视化
GSVA的分析结果相当于表达矩阵,用热图来可视化比较常用,为了体现单细胞的特点我也用FeaturePlot展示了几个基因集。
library(pheatmap)library(patchwork)#绘制热图pheatmap(es.matrix, show_rownames=1, show_colnames=0, annotation_col=meta, fontsize_row=5, filename='GSVA/gsva_demo.png', width=15, height=12)#挑选感兴趣的基因集绘制featureplotes <- data.frame(t(es.matrix),stringsAsFactors=F)scRNA <- AddMetaData(scRNA, es)p1 <- FeaturePlot(scRNA, features = "KEGG_PRIMARY_BILE_ACID_BIOSYNTHESIS", reduction = 'tsne')p2 <- FeaturePlot(scRNA, features = "KEGG_ETHER_LIPID_METABOLISM", reduction = 'tsne')p3 <- FeaturePlot(scRNA, features = "KEGG_RIBOSOME", reduction = 'tsne')p4 <- FeaturePlot(scRNA, features = "KEGG_ASTHMA", reduction = 'tsne')plotc = (p1|p2)/(p3|p4)ggsave('GSVA/featureplot_demo.png', plotc, width = 10, height = 8)

随意挑了几个基因集展示gsva富集分数

GSVA分析给我们的不是阶段性结果,而是观察数据的一个新的视角,让我们以基因集为单位研究数据。基于基因集的富集矩阵,我们可以做差异分析、表型相关性分析、生存分析,也可以基于基因集的功能来鉴定细胞类型。

References
https://www.gsea-msigdb.org/gsea/doc/GSEAUserGuideFrame.html
https://www.gsea-msigdb.org/gsea/msigdb/collection_details.jsp
Laurila E , Mootha V K , Lindgren C M , et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes[J]. 2003, 34(3):267-273.
Subramanian, Tamayo, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. 2005, PNAS 102, 15545-15550.
Hnzelmann S , Castelo R , Guinney J . GSVA: gene set variation analysis for microarray and RNA-Seq data[J]. Bmc Bioinformatics, 2013, 14(1):7-7.

获取帮助

本教程的目的在于把单细胞分析流程串起来,适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释,官网有详细的教程和权威的解释,翻译好的中文教程也有多个版本,有兴趣的可以搜索一下。如果您学习本教程有一定困难,可以分享此篇文章到朋友圈,截图微信发给Kinesin(文末二维码),我会抽时间组群直播讲解单细胞数据分析的全过程。本专题所用的软件、数据、代码脚本和中间结果,我打包放在了百度云上,需要的朋友添加Kinesin微信索取。


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存