企业官网建设流程全解析

GSVA，全称Gene Set Variation Analysis，也就是基因集变异分析，是我这几年处理组学数据时用得最顺手的方法之一。乍一听这个名字，很多人会觉得它和GSEA差不多，但实际上手之后你就会发现，GSVA的思路完全不同：它不需要预先设定分组，而是给每个样本的每个基因集单独算一个“活性分数”，最后输出一张样本乘基因集的矩阵。这张矩阵能接的下游分析远超你的想象，从差异通路、聚类热图，到生存分析和免疫浸润关联，全部可以基于它展开。这篇文章我就把这套方法从原理到实操、从常见报错到参数选择，完整梳理一遍，希望对正在撸表达矩阵的你有点帮助。

这篇文章适合谁看？如果你手头有芯片或者RNA-seq的表达矩阵，想做通路层面的分析但又不确定该用GSEA还是GSVA；或者你已经跑过GSVA，但总觉得参数是照着教程抄的，不知道为什么这么选，也不知道结果到底靠不靠谱——那这篇就是写给你的。文里的代码都是可以直接跑的，我尽量把每一步背后的逻辑也讲清楚，这样你遇到教程没覆盖到的情况时，也能自己判断怎么处理。

1. 项目概述与核心需求解析

1.1 GSVA到底做了什么，它和其他富集分析的本质区别是什么

先举个例子说明问题域。假设你拿到了一批肿瘤样本和正常样本的表达谱，做完差异基因分析后得到几百个上调基因、几百个下调基因。接下来你想知道这些基因整体上影响了哪些生物学通路，通常的做法是富集分析。传统的富集分析（比如基于超几何检验的GO分析）只看“差异基因列表”里有没有富集到特定通路基因集，它把基因当成等权的二元变量，只在“显著差异”和“不显著”之间做判断，结果很容易受阈值选择影响。

GSEA向前走了一步，它利用所有基因的表达变化排序信息，不再需要硬切阈值，能在两组样本之间判断哪些基因集显著富集。但GSEA有一个隐含前提——你必须先有“两组”样本或者一个连续的“表型”来算出基因的排序。这导致它没法在单个样本层面给出通路活性。

GSVA解决的就是这个痛点。它把每个样本单独拿出来，对该样本内部所有基因按表达量做排序，然后针对每个待检验的基因集计算一个富集分数。这个分数是正数代表该基因集在这个样本中整体高表达，负数代表整体低表达。打分过程完全不依赖样本分组信息，是天然的无监督方法。所以GSVA的输出不是一组显著富集的通路列表，而是一张“基因集×样本”的数值矩阵，本质上相当于把基因层面的表达矩阵升维（或者说降维）到了通路层面。

这一步转换的价值非常大。因为下游你无论做无监督聚类、相关性网络，还是做有监督的差异比较，都可以直接用这张通路活性矩阵。比如你想看不同分子亚型的样本在代谢通路上的差异，可以把GSVA分数拿来跑limma或者t检验；你想看某条通路活性跟患者生存是否相关，可以按分数中位数分组画KM曲线；你想分析通路之间的共活跃关系，直接用这个矩阵算Spearman相关系数就行。

1.2 单样本通路活性这个需求为什么难，GSVA是怎么解决的

要理解GSVA的价值，还得先说清楚为什么“单样本通路活性”本质上是个难题。一个基因集的活性不是简单地把里面的基因表达量平均一下，因为基因之间有复杂的调控关系：有的基因是通路的核心酶，表达量波动很小但活性变化巨大；有的基因只是辅助亚基，表达量高但未必代表通路真的激活。如果直接取均值，你会把很多噪音算进去，而且不同通路基因数量差异很大，也不可比。

GSVA的思路是借用了富集分析里经典的Kolmogorov-Smirnov（K-S）统计量的思想。它对每个样本，先把所有基因按表达量高低排好序，形成一个基因列表，然后看目标基因集里的基因是零散分布在排序列表中，还是整体偏向高表达端或低表达端。如果偏向高表达端，说明该基因集在这个样本里被激活，打分就高；反之就低。这和我们人类去看富集图时的直觉完全一致，只是GSVA把这种直觉量化了。

这里要说一个很多人容易误解的点：GSVA并不仅仅是在单个样本内部做基因排序，它其实还考虑了基因在不同样本间表达分布的差异。官方算法里有一个核密度估计步骤，用的是所有样本的数据，然后再回到单个样本计算基因集的富集分数。这也是它和后来的ssGSEA（single-sample GSEA）不一样的地方。ssGSEA更纯粹，只看单个样本内部的排序；GSVA多了一步跨样本的分布估计，所以在样本量足够大的情况下，GSVA对表达值的分布形态更敏感，对数据质量的要求也更高。

1.3 典型的应用场景与适合项目的地方

就我实际用下来的经验，GSVA最常出现的地方有三类。

第一类是肿瘤分子分型与微环境分析。比如用TCGA的转录组数据，把Hallmark基因集和免疫相关基因集都跑一遍GSVA，然后在不同免疫亚型或者不同TMB分组的样本之间比较通路活性的差异，往往能找出几条解释性很强的通路，比如干扰素响应、EMT（上皮间质转化）、炎症反应等。有了GSVA分数，这个过程能非常干净地衔接后续的聚类和差异分析。

第二类是通路层面的共活性网络构建。基因与基因之间的调控网络很难直接推断，但在通路层面，如果两条通路在多个样本中呈现稳定一致的活性变化，那它们背后很可能有共同的调控机制。用GSVA分数矩阵做WGCNA（加权基因共表达网络分析）或者相关网络，比直接用基因表达矩阵做更稳，也更贴近生物学功能层面。

第三类是功能筛选和药敏关联。在药物处理组的转录组数据上跑GSVA，可以快速看出哪些通路被显著激活或抑制；再把GSVA分数和药物IC50做相关性分析，能帮助筛选潜在的药敏标志通路。我自己做过一个项目，就是用GSVA分数矩阵替代原来的基因表达矩阵，跑了一轮随机森林筛选，结果找到的通路标志物在独立验证集里比基因标志物稳定得多。

2. 核心原理与关键参数拆解

2.1 GSVA的算法过程，用非数学的方式理解它

我不打算堆公式，但如果你要调好参数，算法里几个关键节点还是必须搞清楚的。GSVA对每个基因集在每个样本中打分的完整流程大致是这么走的。

第一步，对表达矩阵做标准化处理。这里的标准化不是我们常说的TMM或CPM归一化，而是针对每个基因，在所有样本的表达值上做一个排序变换，让不同基因的表达量分布可比。这一步的意义在于：不同基因本身的表达本底差异很大，有的基因平均表达几百，有的只有几，如果直接在原始值上比较，高表达基因会主导结果。

第二步，用核密度估计（kernel density estimation）拟合每个基因的表达分布。算法里有两种核函数，一种是对称的高斯核（Gaussian），适用于连续型数据，比如芯片的log2信号值或者经过标准化处理的RNA-seq数据；另一种是泊松核（Poisson），适用于原始的RNA-seq整数计数数据，因为计数数据本质上是离散的泊松分布。

第三步，对每个样本，把该样本的所有基因按表达量升序排列。然后对某一个基因集，检查这个基因集中的基因在排序列表中的位置分布。GSVA用两个经验累积分布函数的差值——基因集内部基因的分布和基因集外部基因的分布——来量化这种富集倾向。如果基因集中的基因显著集中在高表达端，那么这两个分布之间的差异会很大，打分就高。

这个打分的过程会同时处理所有样本和所有基因集，所以GSVA的时间复杂度不低。特别是当表达矩阵有几万个基因、基因集几千个、样本几百个时，计算量会明显上升。好在GSVA包支持并行计算，后面实操部分我会给出parallel.sz参数的设置建议。

还有一个细节需要注意：GSVA对基因集大小有敏感度。如果基因集太小（比如少于10个基因），K-S统计量的方差会增大，分数容易虚高或虚低；如果基因集太大（比如超过500个基因），它又会变得过于平滑，区分度下降。所以构建基因集时要适当过滤，最小基因数和最大基因数可以通过后续代码控制。

2.2 四个核心参数的选择逻辑

GSVA的核心函数gsva()里有几个参数，实操中你最需要关心的是kcdf、method、abs.ranking和mx.diff。我逐个说。

kcdf是核密度估计的核函数类型，默认是"Gaussian"。如果你输入的是RNA-seq的原始count数据，我强烈建议把它改成"Poisson"。原因是count数据是离散整数，用连续型的高斯核去拟合会失真。如果你输入的是经过log2转换或TMM标准化的表达矩阵，那就用默认的高斯核没问题。判断标准很简单：你矩阵里如果还能看到类似“257”、“0”、“18"这种整数，那多半是原始count；如果看起来都是负数和带小数的值，像"5.62”、“-1.03”这种，那就是log后的连续值。

method是基因集打分所用的算法，可选"gsva"、"ssgsea"、"plage"和"zscore"。默认是"gsva"，用的是前面介绍的K-S类统计量。"ssgsea"是单样本GSEA方法，它对单个样本内部基因排序后计算富集分数，速度更快，更适用于样本间批次差异较大、或者样本量比较少的情况。我自己在单细胞数据上很少直接用GSVA（太慢），如果非要用也会选ssgsea。"plage"和"zscore"分别是基于基因载荷和Z分数的简化打分方法，速度最快但信息量损失较大，一般做初步筛选时用。

abs.ranking参数控制的是排序时是否取表达值的绝对值。默认FALSE，表示按表达量的原始值排序。如果设为TRUE，它会先对表达值取绝对值再排序，这时正负表达变化的基因都会被同等对待。这个参数在芯片数据里应用比较多，因为芯片数据有上下调倍数变化的概念；在RNA-seq标准化数据中一般用默认值就行。

mx.diff控制的是打分时是否对K-S统计量做最大化处理。默认TRUE，推荐保持默认，这样分数范围更宽，区分度更好。如果你发现结果中大部分基因集的分数都在零附近挤成一团，可以检查一下是不是这里被改了。

2.3 输入数据的前期处理和基因集准备

GSVA对输入表达矩阵有一些硬性要求，踩过坑的人都知道。首先，表达矩阵必须是数值矩阵，行名是基因符号或Entrez ID，列名是样本名，中间不能有NA值。如果矩阵里有NA，GSVA会直接报错或者输出NaN。我习惯在跑GSVA之前先做一遍过滤：删除在所有样本中表达量都为零的基因，删除NA比例超过5%的基因。RNA-seq数据最好先做低表达过滤，这不是GSVA明确要求的，但能明显减少核密度估计时出现极端值的概率。

其次，基因名的格式需要和你的基因集保持统一。如果表达矩阵用基因符号，那么基因集最好也用基因符号；反之如果基因集是Entrez ID，表达矩阵的行名也得是Entrez ID。很多人在这里翻车：从MSigDB下载的GMT文件标的都是Entrez ID，而表达矩阵却是基因符号，结果一跑，几百个基因集全部匹配不上，输出的矩阵几乎全是零。

基因集准备方面，我推荐用msigdbr这个R包，它能直接以数据框形式返回MSigDB里各个集合的基因列表，而且自动帮你标好基因符号和Entrez ID两种格式，比去官网手动下载GMT文件再解析方便得多。常用的基因集有Hallmark（50条核心通路）、KEGG（代谢和信号通路）、GO的BP（生物学过程）等。如果你做的是免疫微环境分析，还可以用文献里整理的免疫基因集列表，做成list对象直接喂给GSVA。

3. 实操过程与核心环节实现

3.1 环境准备与R包安装

在动手跑GSVA之前，先把环境搭好。GSVA是一个Bioconductor包，推荐直接用BiocManager安装。

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("GSVA") BiocManager::install("msigdbr")

除了这两个核心包，建议把limma（做差异通路分析）、pheatmap（热图）、survival和survminer（生存分析）也一起装好。这些都是后续分析会用到的基础工具。

装好之后加载：

library(GSVA) library(msigdbr) library(limma) library(pheatmap)

我这里用一套模拟数据来演示完整流程。实际项目中你替换成自己的表达矩阵就行。

3.2 标准GSVA运行流程

先构造一个演示用的表达矩阵。假设我有20000个基因、40个样本，其中前20个是肿瘤组，后20个是正常组。我随机生成一部分差异基因来模拟生物学变异。

set.seed(123) n_genes <- 20000 n_samples <- 40 expr_matrix <- matrix(rnbinom(n_genes * n_samples, size = 10, prob = 0.8), nrow = n_genes, ncol = n_samples) rownames(expr_matrix) <- paste0("GENE", 1:n_genes) colnames(expr_matrix) <- paste0("SAMPLE", 1:n_samples) # 模拟肿瘤样本中某些基因上调 diff_genes <- sample(1:n_genes, 2000) expr_matrix[diff_genes, 1:20] <- expr_matrix[diff_genes, 1:20] * 3

实际的RNA-seq计数数据应该先做归一化。这里我用edgeR的CPM加log2转换做演示：

library(edgeR) # 注意：真实项目中应该用DGEList对象原始count来做归一化 dge <- DGEList(counts = expr_matrix) dge <- calcNormFactors(dge) cpm_log <- cpm(dge, log = TRUE)

GSVA包也支持直接输入原始count并指定kcdf="Poisson"。我自己的经验是：如果数据样本量比较大（超过50），用原始count加Poisson核效果不错；样本量小的时候，还是先做CPM-log转换再用高斯核更稳。两种方式都可以，关键是保持同一个项目里前后一致。

接着获取基因集。这里以Hallmark基因集为例：

hallmark <- msigdbr(species = "Homo sapiens", category = "H") hallmark_list <- split(hallmark$gene_symbol, hallmark$gs_name)

msigdbr返回的category="H"就是50条Hallmark通路。split()直接按通路名转换成list，每个元素就是一个基因集，这正是GSVA需要的基因集格式。如果你要用KEGG，把category改成"C2"，并且subcategory设为"CP:KEGG"。

真正跑GSVA就一行：

gsva_res <- gsva(expr = cpm_log, gset.idx.list = hallmark_list, method = "gsva", kcdf = "Gaussian", abs.ranking = FALSE, parallel.sz = 4, verbose = TRUE)

parallel.sz=4表示用4个并行线程，如果你的机器核心多，可以设大一点。verbose=TRUE会打印运行进度，样本量和基因集多的时候很有用。

跑完后的gsva_res是一个矩阵，行名是通路名，列名是样本名，每个单元格是该通路在该样本中的活性分数。我拿到这个矩阵后的第一件事，永远是检查分布：

summary(apply(gsva_res, 1, sd))

如果很多通路的标准差趋近于零，说明基因集和表达矩阵匹配度太低，或者表达矩阵本身区分度太差。正常的GSVA结果中，大部分通路的标准差应该在0.3以上（对数尺度下）。这一步检查能帮你避免拿着一个全是零方差的结果往下走。

3.3 结果解读与常用可视化

GSVA矩阵的可视化最常用的是热图。做热图之前建议先对通路做聚类，不然50条通路乱排在一起根本看不出模式。

# 做一个简单的分组注释 sample_anno <- data.frame(Group = rep(c("Tumor", "Normal"), each = 20)) rownames(sample_anno) <- colnames(gsva_res) # 筛选变化较大的通路，例如方差前25条 top_var <- order(apply(gsva_res, 1, var), decreasing = TRUE)[1:25] gsva_heatmap <- gsva_res[top_var, ] pheatmap(gsva_heatmap, annotation_col = sample_anno, show_colnames = FALSE, scale = "row", color = colorRampPalette(c("#4575B4", "white", "#D73027"))(100))

注意这里我用了scale="row"做行内标准化。这个非常关键：GSVA打分本身是相对值，不同通路之间的绝对分数不能直接比大小，但同一通路在不同样本之间的相对高低是可靠的。行内标准化之后，热图展示的就是“这条通路在哪些样本里相对更活跃”，生物学解释更直观。

除了热图，我还经常画分组箱线图和小提琴图。比如想比较肿瘤和正常样本中EMT通路活性的差异：

library(ggplot2) library(reshape2) emt_scores <- data.frame(Group = sample_anno$Group, Score = as.numeric(gsva_res["HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", ])) ggplot(emt_scores, aes(x = Group, y = Score, fill = Group)) + geom_boxplot() + geom_jitter(width = 0.2, size = 0.5) + theme_bw() + labs(title = "EMT pathway activity", y = "GSVA score")

从这种图上你能很直观地看到通路在不同组间的活性差异，这也是很多高分论文里标准的GSVA结果展示形式。

3.4 从GSVA分数到差异通路与生存分析，一套代码串起来

GSVA矩阵拿到手之后，最常做的下游分析是差异通路分析和生存分析。先说差异通路分析。

用limma对GSVA矩阵做差异分析，和做差异基因分析的流程完全一致，只是输入从基因表达矩阵换成了通路活性矩阵：

design <- model.matrix(~ 0 + factor(sample_anno$Group)) colnames(design) <- c("Normal", "Tumor") fit <- lmFit(gsva_res, design) contrast_mat <- makeContrasts(Tumor - Normal, levels = design) fit2 <- contrasts.fit(fit, contrast_mat) fit2 <- eBayes(fit2) diff_pathways <- topTable(fit2, coef = 1, number = Inf, adjust.method = "BH")

diff_pathways里每一行是一条通路，logFC表示肿瘤相对于正常组织的通路活性变化。为了减少假阳性，一般把adj.P.Val < 0.05且abs(logFC) > 0.2作为筛选阈值。注意这里logFC的阈值不能套用差异基因时的1或2，因为GSVA分数的变化范围本身只有±3左右，0.5的差异已经算明显了。

生存分析也很常见。以EMT通路为例，先按分数中位数把样本分成高低两组，再看生存差异：

library(survival) library(survminer) # 假设你有一个生存数据框surv_df，包含time和status列，行名是样本名 score <- gsva_res["HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", ] median_score <- median(score) group <- ifelse(score >= median_score, "High", "Low") surv_df <- merge(surv_df, data.frame(Group = group), by = "row.names") fit_km <- survfit(Surv(time, status) ~ Group, data = surv_df) ggsurvplot(fit_km, data = surv_df, pval = TRUE, risk.table = TRUE, title = "EMT pathway activity and survival")

这种基于通路的生存分析比单基因的生存分析稳定得多，因为通路活性是众多基因的综合反映，对单个基因的表达噪声不敏感。我在实际项目里做过对比，同一个数据集中，通路水平生存分析的重现性比单基因水平明显好，尤其在独立验证集里。

如果想进一步提升可解释性，可以把GSVA的分数矩阵和临床变量（年龄、分期、性别）做相关性矩阵，用corrplot画出相关热图。这一步往往能发现一些很有价值的关联，比如某条代谢通路活性和TNM分期显著正相关，那就是一个可以深挖的切入点。

4. 常见问题与排查技巧实录

4.1 典型报错与解决方案速查表

GSVA虽然用起来简单，但跑的过程中报错和警告是家常便饭。下面整理一份我遇到过的典型问题速查表，按出现频率排序。

最坑的一个问题我单独拎出来说：基因名重复。表达矩阵里如果行名有重复基因名，GSVA不会主动报错，但结果会悄悄地不可靠。因为算法里要对每个基因做排序，重复的基因名会导致部分信息被覆盖或错位。我一般先用limma包的avereps函数按基因名取均值合并去重，这一步在跑GSVA之前一定要做。

# 假设expr_matrix行名有重复基因符号 expr_matrix <- avereps(expr_matrix, ID = rownames(expr_matrix))

4.2 数据质量层面的避坑心得

GSVA对输入数据的质量敏感程度比我想象中高，这也是很多新手跑不出理想结果的原因。首要问题是样本量。GSVA核密度估计需要一定的样本量才能对每个基因的表达分布做出稳定拟合。我自己的经验是样本量少于15个时，GSVA结果会明显不稳定，换一个样本子集跑，通路排序会变；样本量超过30个，结果就相当稳定了。所以如果你的项目只有8个样本，我会建议你谨慎使用GSVA，或者至少要认识到结果的方差会偏大。

批次效应也是一个隐藏杀手。因为GSVA会利用所有样本的表达分布来调整打分，如果样本来自两个批次，表达分布整体偏移，GSVA分数可能会被批次因素主导而不是生物学差异主导。跑GSVA之前最好先用removeBatchEffect()把已知的批次协变量去掉，或者用sva包估计隐藏的批次因子。

对照组的选择也会影响结果。GSVA打分是相对于当前样本集分布的，如果你把一组纯肿瘤样本放进GSVA，得到的分数反映的是这条通路在这个样本里相对于这批肿瘤样本的高低，而不是绝对活性。所以不同数据集之间跑出来的GSVA分数不能直接比较，这点在泛癌分析时尤其要注意——每个癌种单独跑GSVA，然后看通路活性模式是否一致，而不能把多个癌种的表达矩阵合在一起跑GSVA再互相比较，除非你有充分的批次校正理由。

4.3 结果可靠性评估的土办法

我不太建议拿到GSVA结果就直接往下游放。在真正开始差异分析之前，我会先用几个“土办法”快速校验结果的可靠性。

第一是抽一个生物学已知明确的方向来验证。比如免疫相关数据中，肿瘤样本的干扰素响应通路活性应该普遍升高；如果是抗感染实验，炎症相关通路应该在处理组升高。如果这些常识性的预期都没有出现，那要么基因集选错了，要么表达数据有问题，要么匹配率太低。先验证一条已知通路，比检查一百个统计指标都直接。

第二是看生物学一致方案的重复性。如果你同一组实验有生物学重复样本，理论上这些重复样本的GSVA分数应该比较接近。我会算一下同一组内样本间的平均相关系数，如果低于0.5，就要警惕数据质量问题。有一条简单标准：生物学重复的样本在GSVA分数热图里应该聚在一起，如果连重复样本都四散分布，这个数据跑下游意义不大。

第三是换了参数重新跑一次。比如分别用method="gsva"和method="ssgsea"各跑一遍，然后看主要结论是否一致。两种方法在算法上有差异，但对真正强的通路活性信号，结果应该大同小异。如果某个关键通路的活跃方向在不同方法下截然相反，那大概率是数据本身的问题，而不是算法差异。

这些土办法虽然不“论文级严谨”，但在项目早期能帮你快速止损，避免在一个错误的数据基础上花几周时间做下游分析。

5. 一点个人体会与扩展建议

跑GSVA这几年，我最大的体会是：它不是一个“跑完就完事”的工具，而是整个通路层面分析的第一棒。GSVA分数矩阵的价值在于它把几百个基因的复杂变化压缩成了几十个通路的简洁画像，让后续的分析变量大大减少，又保留了生物学可解释性。我自己现在的标准流程是：拿到表达矩阵后先跑GSVA，用热图快速扫描样本分组模式下哪些通路有明显变化，锁定三四条关键通路，再回到基因层面去看这些通路内部哪些核心节点在驱动信号变化。这个“通路→基因”的层级分析思路，比一开始就在几万个基因里大海捞针要高效得多。

最后再分享一个小技巧。GSVA算出的通路活性分数，还可以拿来做主成分分析（PCA）。用GSVA分数矩阵代替基因矩阵做PCA，往往能得到比基因层面更清晰的分组分离。因为通路层面的信号已经是去噪浓缩后的结果。我最近一个项目里，基因表达PCA只能勉强区分亚型，但GSVA分数的PCA第一主成分就把两个亚型干净地分开了。这算是一个性价比很高的加分项，你在自己的项目里不妨试试。