上期专题我们介绍了单细胞转录组数据的基础分析，然而那些分析只是揭开了组织异质性的面纱，还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容：多样本批次校正、转录因子分析、细胞通讯分析、基因集变异分析和更全面的基因集富集分析。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨!
GSEA全称Gene Set Enrichment Analysis，GSVA全称Gene Set Variation Analysis，它们都是基于基因集开展的分析，因此我们先要了解基因集的定义。基因集顾名思义就是一些基因的集合，任何一些基因放在一起都可以叫做基因集，但是我们用来分析的基因集要求有一定的生物学意义。最常见的基因集数据库如GO与KEGG，它们一个按照基因本体论将基因分门别类，一个按照代谢通路将相关基因集合在一起。除此之外我们还可以按转录因子调控网络、共表达网络、定义生物状态的marker基因列表等把基因集合成有一定生物学意义的基因集。MSigDB基因集数据库 GSEA是由Broad研究所开发的一种富集方法，他们在提出该方法的同时还提供了一个基因集数据库——MSigdb。它从位置，功能，代谢途径，靶标结合等多种角度出发，构建出了许多的基因集合，Broad研究所将他们构建的基因集合保存在MSigDB，官网地址如下：http://software./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免疫基因集，包含了免疫系统功能相关的基因集合。 常规GO/KEGG富集分析需要设定阈值过滤差异基因，阈值太宽富集的结果太多，阈值太严又可能会遗漏一些关键结果。GO/KEGG富集的结果通常还很宽泛，并不能很好地解释生物学现象。有鉴于此，Broad研究所开发了基因集富集分析(GSEA)方法。GSEA使用无监督算法，不用过滤任何基因，配合MSigDB数据库使用，更容易找到解释生物学现象的基因集。其原理如下：GSEA分析要先将样本做组间对比分析，GSEA自带9种分析方法，分为基因表达值差异分析和相关性分析两大类。对于Case/Contral的实验设计，差异分析方法更为常用，这其中又以默认的信噪比和大家熟悉的差异倍数用的最多。对比分析之后要按结果将基因排序，以差异倍数方法为例，把所有基因按差异倍数(FC)的值降序排列以供后续分析。上图小人脚下的小方块代表排序好的差异基因列表，蓝色之外的其他色块代表属于某个基因集的基因，如黄色属于基因集A，绿色属于基因集B。最下面高低不等的竖条代表与基因列表对应的FC值，红色上调、蓝色下调、黄色没有变化。基因集的富集分析需要经历三步：基因集A富集分析时，小人从基因列表的左端走到右端，每经过一个蓝色基因扣分，每遇到一个黄色基因加分，扣分时与FC无关，加分时考虑FC的权重。基因集A最终的富集分数(ES)是小人曾经得过的最高/低分，实际公式比这复杂，但基本理念如此。 采用置换检验计算基因集A的显著性，即p值。 基因集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. GSEA虽然是一种强大的富集分析工具，但是它的应用场景通常局限于Case/Control的实验设计。对于表型（分组）复杂的大样本量研究，例如TCGA和单细胞图谱这样的项目，分析起来就非常困难。因此，Broad研究所在GSEA发布8年之后，又开发了GSVA算法来拓展基因集分析的应用。GSVA不需要预先进行样本之间的差异分析，它依据表达矩阵就可以计算每个样本中特定基因集的变异分数。简单的说，输入以基因为行的表达矩阵和基因集数据库给GSVA，它就输出以基因集名称为行的变异分数矩阵，如下图所示：左侧输入基因表达矩阵和基因集数据库，中间是GSVA算法原理，右侧是输出的基因集变异分数矩阵。基因集变异分数可以理解为基因集内所有基因的综合表达值。
GSEA官方提供很多版本，但是不提倡R语言版，下载网址如下：
https://www./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')
分组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分数最高的通路的富集结果： 上图对应数据 提前安装好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/
表达数据和分组文件准备好之后，图形界面版操作没有难度，简单介绍一下：
点击箭头所指按钮之后，就会弹出对话框让你选择文件，分次输入我们之前准备好的表达数据文件(*.gct)和分组文件(*.cls)选好文件之后，文件目录会出现在左右两边箭头所指的地方。如果右边没有，双击左边即可。我这儿左边有三个文件是因为加载了本地保存的基因集数据库文件。
文件加载完毕之后，点击箭头所指运行会出现右边的对话框，注意一下红框标注的选项。Collapse/Remap项选择“No_Collapse”（Collapse适用芯片数据），Permutation项选择“phenotype”（更适合单细胞数据），Chip项可不选或选择有“Human_Symbol_with_”关键字的项目（选择之后结果中基因有简短注释）。
Basic fields点开之后选择保存结果的目录，点击“Run”就开始运行了。java版（图形界面的也是java版）运行之后结果与R语言版一致，不过分析报告是网页版，点击超链接就可以查看细节，使用起来非常方便。
#从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分析给我们的不是阶段性结果，而是观察数据的一个新的视角，让我们以基因集为单位研究数据。基于基因集的富集矩阵，我们可以做差异分析、表型相关性分析、生存分析，也可以基于基因集的功能来鉴定细胞类型。https://www./gsea/doc/GSEAUserGuideFrame.htmlhttps://www./gsea/msigdb/collection_details.jspLaurila 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语言基础的朋友参考。分析原理和代码我没有详细解释，官网有详细的教程和权威的解释，翻译好的中文教程也有多个版本，有兴趣的可以搜索一下。如果您学习本教程有一定困难，可以点击 “阅读原文” 找到作者联系方式，获取帮助。
|
