解析细胞功能基因集变异——遇见GSVA|生信开发实战

企业   2024-12-05 17:00   浙江  

什么是GSVA?

GSVA,即Gene Set Variation Analysis,基因集变异分析,是一种非参数的无监督分析方法,主要用来评估芯片和转录组的基因集富集结果。与GSEA不同,GSVA不需要预先对样本进行分组,可以计算每个样本中特定基因集的富集分数。GSVA主要是通过将基因在不同样品间的表达量矩阵转化成基因集在样品间的表达量矩阵,从而来评估不同的代谢通路在不同样品间是否富集。

不同于GSEA之处的在于,对于不同的数据类型(只支持log表达值或原始的read counts),假设了不同的累积密度函数(cumulative density function,CDF),针对芯片数据,假设为正态分布密度函数;针对RNA-seq数据,假设为泊松分布密度函数,而且,GSVA是为每个样本的每个基因计算对应的CDF值,然后根据该值对基因进行排序,这样,每个样本都有一个从大到小排序的基因列表。

GSVA的应用组学多样,bulk层面的转录组、代谢组、蛋白组等均可使用,但单细胞组学需考虑数据的稀疏程度或细胞数量,特别是后者,以细胞为单位进行分析时,GSVA的分析时间和内存占用量会急剧上升。

 

基本原理

GSVA算法接受的输入为基因表达矩阵(经过log2标准化的芯片数据或者转录组测序的reads count数据)以及特定基因集。第一步,算法会对表达数据进行核密度估计;第二步,基于第一步的结果对样本进行表达水平排序;第三步,对于每一个基因集进行类似K-S检验的秩统计量计算;第四步,获取GSVA富集分数(ES值)。最终输出为以每个基因集对应每个样本的数据矩阵。

下面以转录组测序中的基因表达谱的GSVA分析进行介绍:

文件准备

由于分组是非必须的,因此GSVA相对于GSEA无需分组文件。GSVA需要以下两个文件:

  - 表达谱文件,无论是基因的reads count,还是基因标准化的表达量,如FPKM、TPM值等,需要准备表达量矩阵;

  - 基因集文件,GSEA官方提供的人(Homo sapiens,hsa)的基因集,可以在http://www.gsea-msigdb.org/gsea/login.jsp下载,GSEA提供基因ID(entrez)和基因名(symbols)两种形式,可以基于基因表达矩阵的形式选择对应格式的基因集(一般选择symbols即可)。需要注意的是基因集中的基因名形式需要和表达谱中一致。

除了hsa外,其他物种官方并没有提供基因集文件,需要使用者自己制作基因集文件。考虑到GSEA分析需求的广泛增长,可以利用基因的GO、KEGG注释生成gmt文件并用于GSEA富集分析。gmt格式是多列注释文件,第一列是基因所属基因集的名字,可以是通路名字,也可以是自己定义的任何名字。第二列,官方提供的格式是URL,可以是任意字符串(如NA)。第三列之后是基因集内基因的名字,有几个写几列。列与列之间都是tab键分割。GSVA的gmt文件和GSEA是一致的。下为基因集示例图,每一行代表一个基因集(前面数字编号为行编号):

分析代码

##安装GSVA### BiocManager::install("GSVA",force =TRUE)

##载入分析、数据处理和可视化包##library(GSVA)library(GSEABase)library(pheatmap)library(ggplot2)library(dplyr)library(tidyr)

##读取表达量矩阵,以标准化的基因表达量矩阵为例##exp <- read.table('expression_matrix_fpkm.txt', header = T, sep = "\t")

##对表达量矩阵的相同基因名的基因进行处理,比如取均值或中位数##exp.mean <- aggregate(x = exp[, 2:ncol(exp)], by = list(exp$gene_name), FUN = mean)rownames(exp.mean) = exp.mean[,1]exp.mean <- exp.mean[,-1]

##对表达量取log,GSVA支持log的表达值##exp.log <- log2(exp.mean+1)

##读取基因集gmt文件##geneset <- getGmt('gene_sets.gmt')

##进行GSVA分析,使用标准化的表达量时,kcdf为高斯分布,即正态分布##es <- gsva(as.matrix(exp.log), geneset, min.sz=10, max.sz=500, verbose=TRUE, kcdf="Gaussian")

##输出GSVA分析矩阵##write.csv(es, "gsva_result.csv")

##进行热图绘制##plot<-pheatmap(es, cluster_cols = T, cluster_rows = T, scale = "row", cutree_rows = NA, cutree_cols = NA, color=colorRampPalette(c("green","black","red"))(100), legend=T, border_color=NA, cellwidth = 20, cellheight = 6, fontsize = 6, fontsize_row = 3, fontsize_col = 9 )ggsave("gsea.heatmap.pdf",plot,width=30,height=30,units="in")
结果说明

GSVA分析运行结束后,核心结果为各样本的基因集打分情况,即以下矩阵文件——矩阵第一列为基因集,如果是KEGG数据库,则为通路名,后面几列为基因集在各个样本中表达的富集分数,与第一列一一对应。需要注意的是,相同通路的分数可以在不同样本中比较,富集分数越高,代表该通路在对应样本中的功能更活跃。但是同样本不同通路的分数之间没有比较意义。

获得结果矩阵后,一方面可以在组间进行差异分析,比较具有显著变化的基因集,以箱线图或柱状图的形式进行展示,另一方面,也可以直接以热图形式直观展示各通路在不同样本中的富集情况,效果如下:

以上为GSVA常用的分析与展示方法,该工具已经整合到联川云平台中(https://www.omicstudio.cn/tool/122),欢迎各位使用!



相关阅读


国自然2025:单细胞+bulk转录组如何开展研究?|生信开发实战


单细胞评分分析R包Ucell


国自然2025:空间组学和Xenium原位技术如何实现细胞分割?


单细胞细胞通讯分析之CellChat v2


本文系联川生物公众号原创文章,未经授权禁止转载,侵权必究!

扫描下方二维码





点分享


点点赞


点在看


联川生物
一个提供科研入门学习资源、经验的平台。 分享前沿测序技术资讯、实用生信绘图技巧及工具。 发布高质量的科研论文精度、精炼科研思路。 我们的目标是持续提供“干货”,滋润您的科研生涯。
 最新文章