简介面向零基础转录组学习者的基因集变异分析GSVA配套资源包含输入数据、分析脚本和结果输出三大模块。脚本已提前测试可一键全选跑通适合想快速上手富集分析的生信新手基础薄弱者可结合配套博文教程边学边练逐句理解代码逻辑。资源共十个文件以六个CSV数据表格为核心涵盖表达矩阵、分组信息、富集评分以及显著性筛选结果另含一个R语言脚本、一个GMT基因集注释文件、一份PDF结果报告和一张PNG可视化图片压缩包整体为106.43MB目录划分清晰便于对照输出文件自查。目前已有206人学习使用是完成转录组下游功能富集分析的实用入门资料能帮助读者省去整理数据和调试代码的时间把精力集中在理解生物学意义上。1. 为什么转录组做完差异分析还不够拿到一组 FPKM 表达矩阵差异基因也筛完了KEGG 和 GO 富集也跑了一轮这时候很多人才意识到一个问题差异富集看的是“哪些通路里有差异基因”但它回答不了“这条通路在样本里到底被激活还是被抑制”。GSVAGene Set Variation Analysis就是用来补这块短板的——它不依赖差异基因列表而是直接对全部基因的表达量做无监督打分算出每个基因集在每个样本里的富集分数。这决定了它特别适合处理那些组间差异不剧烈、但整体表达趋势有偏移的数据比如肿瘤分型、药物处理的时间序列、或者临床样本的亚群比较。这套配套资源是围绕一个完整可跑的 R 脚本r.GSVA.R展开的输入是三份文件表达矩阵data_fpkm.csv、样本分组data_group.csv、基因集c2.all.v2023.2.Hs.symbols.gmt输出则包括 GSVA 分数矩阵、显著性筛选结果、以及可直接用于论文的 PDF/PNG 热图。对于刚接触转录组下游分析的人来说最大的门槛其实不是函数调用而是搞不清 GSVA 和 GSEA 的区别、不知道基因集格式怎么处理、以及拿到gsva_res.csv之后下一步该干嘛。这篇文章就按“原理—数据准备—脚本逐段拆解—结果解读—常见报错”这条线把这套资源讲透。提示文中所涉及的代码均从配套r.GSVA.R的完整逻辑中拆出你可以直接复制到 RStudio 里分段执行也可以全选一键跑通。2. GSVA 算法原理与 GSEA 的边界划分2.1 GSVA 的统计思想从表达矩阵到通路活性矩阵GSVA 的核心是把“基因表达量”转换成“基因集富集分数”。它先对每个样本内部的全部基因表达值做核密度估计再把基因按照表达量排序然后用类似于 Kolmogorov-Smirnov 检验的思路计算每个基因集在这个排序中的富集程度。注意这里的“富集”不是看基因是否属于差异基因而是看基因集内的基因整体上是否偏向表达量排序的某一端——偏向高表达端说明该通路在该样本中被激活偏向低表达端则说明被抑制。GSVA_score(sample_j, gene_set_k) max(ES_kj)其中 ES 是基因集内基因在样本 j 表达排序中的富集得分计算时会对每个基因集独立进行所以最终输出的GSVA_result.csv是“基因集 × 样本”的矩阵每一格都是一个连续的数值正负代表激活方向绝对值大小代表强度。这也解释了为什么 GSVA 不需要预先定义“差异基因”它天然是全局分析。2.2 GSVA 与 GSEA 的根本差异GSEA 在做富集分析时输入的是一份带有统计量的基因排序列表通常是 log2FC 或 t 值它比较的是两个生物学状态之间的差异富集而 GSVA 是对每个样本独立打分输出的分数可以直接用于后续的聚类、相关性分析或者生存分析。一句话概括GSEA 是“组间比较”GSVA 是“样本级打分”。维度GSVAGSEA输入核心表达矩阵 基因集基因排序列表 基因集打分粒度每个样本每个基因集每个对比组合每个基因集是否需要分组不需要但分组有利于后续统计必须有两组及以上输出形式样本 × 通路矩阵富集得分及显著性典型应用分型、相关性、生存分析差异通路筛选理解了这个区别你就明白为什么在这套资源里data_group.csv并不是 GSVA 计算的必需输入——它只在后续做组间显著性筛选时才用到。2.3 ssGSEA 与 GSVA 的适用边界配套资源里只用了 GSVA 函数但在实际分析中很多人会把它和 ssGSEAsingle-sample GSEA混用。ssGSEA 同样是单样本打分但它的算法是直接对基因集内的基因做经验 CDF 累积计算复杂度低于 GSVA而且在免疫浸润分析中更常用。如果你手头的数据是芯片而不是 RNA-seq表达谱经过 log2 归一化后GSVA 和 ssGSEA 的结果差异并不大但如果是 FPKM 这种带有长度归一化的数据GSVA 对表达分布的假设更稳所以我一般建议用 GSVA 而不是直接替换成 ssGSEA。3. 输入数据准备GMT 基因集与表达矩阵的格式陷阱3.1 基因集文件 c2.all.v2023.2.Hs.symbols.gmt 的结构GMT 文件每行代表一个基因集三列以上第一列是基因集名称比如KEGG_OXIDATIVE_PHOSPHORYLATION第二列是描述信息通常是一个网址或简要说明第三列开始才是基因符号。这套资源里用的是 MSigDB 的 c2 集合v2023.2 版本包含 KEGG、Reactome、WikiPathways 等来源的通路条目基因标识是Hs.symbols也就是标准的 Human Gene Symbol。KEGG_OXIDATIVE_PHOSPHORYLATION https://www.gsea-msigdb.org/gsea/msigdb/human/geneset/KEGG_OXIDATIVE_PHOSPHORYLATION.html ATP5MC1 ATP5MC2 ATP5MC3 ...注意这里有个坑如果你的表达矩阵里是 Ensembl ID 或 Entrez ID不能直接和这个 GMT 匹配。常见的做法是先做 ID 转换把 Ensembl 转成 Gene Symbol再删掉空匹配和重复行。提示检查表达矩阵的行名和 GMT 里的基因符号是否同一种 ID 格式这一步错了后面所有结果都是空的。3.2 data_fpkm.csv 的标准化要求FPKM 数据的特征是没有经过跨样本标准化基因的表达量受文库大小和基因长度双重影响。GSVA 函数内部会对每个样本的表达向量做排序所以理论上 FPKM 可以直接进 GSVA但对数化通常能带来更稳健的结果。常见的预处理方式是先做log2(FPKM 1)再对每个样本做 z-score 标准化。# 读取表达矩阵行名为基因列名为样本 expr - read.csv(data_fpkm.csv, row.names 1) # 建议的预处理流程先取对数再按行进行标准化 expr_log - log2(expr 1) expr_z - t(apply(expr_log, 1, function(x) { (x - mean(x)) / sd(x) }))这段代码里apply按行MARGIN 1对每个基因在所有样本中的表达量做 z-score 变换这样做的本质是消除不同基因之间的表达丰度差异让 GSVA 排序比较的是“基因内跨样本的相对变化”而不是绝对表达值。对 FPKM 来说这一步非常重要否则那些高表达基因如看家基因会主导排序方向。3.3 data_group.csv 的写入格式data_group.csv用于后续的分组比较最简单的格式是两列第一列是样本名第二列是分组标签。注意样本名必须和表达矩阵的列名完全一致包括顺序和大小写。如果你直接从 Excel 复制出来很容易在列名里带上不可见字符或者把分组列变成因子导致后面compare环节报错。sample,group Ctrl_1,Control Ctrl_2,Control Treat_1,Treatment Treat_2,Treatment4. 核心脚本 r.GSVA.R 逐段拆解参数与执行逻辑4.1 加载 GSVA 包与基因集解析脚本开头加载的关键包是GSVA和GSEABase。GSVA提供gsva()主函数GSEABase负责把 GMT 文件读成标准的基因集对象GeneSetCollection。这一步没有必要手动逐行解析 GMT直接用getGmt()就行但要注意版本差异——老版本GSEABase读取的GeneSet对象里基因类型默认是SYMBOL如果你的表达矩阵行名也是符号那么不需要额外指定geneIdType。library(GSVA) library(GSEABase) # 读取基因集文件 gene_sets - getGmt(c2.all.v2023.2.Hs.symbols.gmt) # 把 GeneSet 列表转换为 GeneSetCollection 对象 gs_collection - GeneSetCollection(gene_sets)这里getGmt()返回的是一个GeneSet列表GeneSetCollection()把它统一成 GSVA 官方接口需要的格式。如果你读到这一步报错说unable to find an inherited method for function getGmt通常是GSEABase没有成功加载而不是文件路径的问题。4.2 调用 gsva() 主函数参数 scRNAseq 与 kcdf 的含义这是整套资源里最重要的一个函数调用# 执行 GSVA 打分 gsva_res - gsva( expr as.matrix(expr_z), gset.idx.list gs_collection, kcdf Gaussian, method gsva, min.sz 5, max.sz 500, verbose TRUE )expr参数要求是数值矩阵不能是数据框所以这里用as.matrix()做显式转换。矩阵的行名必须和基因集里的基因符号能匹配上GSVA 内部会自动取交集没有交集的基因会被丢弃。kcdf参数这是 GSVA 特有的核密度估计类型。对于 RNA-seq 的 FPKM 数据Gaussian是常规选择。如果你的数据是直接的 counts 矩阵可以考虑Poisson但 FPKM 是连续值用 Poisson 反而不合适。method参数gsva表示使用原始 GSVA 算法设置为ssgsea则计算单样本 GSEA。这里用默认的gsva因为我们要的是完整的通路活性矩阵。min.sz和max.sz基因集大小过滤范围。min.sz 5表示基因数少于 5 的基因集直接丢弃max.sz 500表示大于 500 的也会被过滤掉。前者是为了避免小基因集打分不稳定后者是为了避免过于宽泛的通路把信号稀释掉。c2 集合里有些通路非常大比如某些 Reactome 条目有上千个基因这些通路在转录组数据里往往没有生物学特异性所以过滤掉是合理的。4.3 输出 GSVA_result.csv 与 gsva_res.csv 的区别脚本里会生成两个看起来很像的文件GSVA_result.csv和gsva_res.csv。以我拿到资源后的核对来看前者通常是未经后续筛选的完整打分矩阵后者则是经过分组差异比较后选出的显著基因集打分矩阵。简单说GSVA_result.csv是中间产物行是基因集列是样本gsva_res.csv是用于画图的下游输入。# 写出完整结果 write.csv(gsva_res, GSVA_result.csv) # 后续筛选显著基因集后再次写出 write.csv(gsva_sig, gsva_res.csv)4.4 分组显著性判断从 GSVA 分数到差异通路拿到gsva_res矩阵之后脚本会按data_group.csv的分组信息做差异比较。常见做法是对每个基因集的行向量做 t.test 或 Wilcoxon 检验然后筛选 p 值小于阈值比如 0.05的基因集。group_info - read.csv(data_group.csv, row.names 1) group - factor(group_info$group) # 对每个基因集做 Wilcoxon 秩和检验 pvals - apply(gsva_res, 1, function(row) { wilcox.test(row[group levels(group)[1]], row[group levels(group)[2]])$p.value }) # 多重检验校正 padj - p.adjust(pvals, method BH) # 筛选显著基因集 gsva_sig - gsva_res[padj 0.05, , drop FALSE]这里用 Wilcoxon 而不是 t 检验是因为 GSVA 分数不一定是正态分布尤其在小样本量每组 3~5 个样本的转录组数据里秩和检验更稳健。p.adjust的BH方法控制 FDR适合高维假设检验场景。提示如果每组只有 2 个样本wilcox.test不会报错但 p 值没有区分度建议至少每组 3 个样本再跑显著性筛选。5. 可视化输出从 pheatmap 到论文级热图5.1 01.GSVA_res.pdf / 01.GSVA_res.png 的绘图逻辑脚本用pheatmap绘制聚类热图输出同时存 PDF 和 PNG 两个版本PDF 用于论文投稿不失真PNG 用于快速预览。绘图数据就是gsva_sig矩阵按基因集做聚类样本不聚类保持原有的分组顺序会更便于观察组间差异。library(pheatmap) # 绘制显著基因集的热图 pheatmap( gsva_sig, scale row, cluster_rows TRUE, cluster_cols FALSE, show_rownames TRUE, show_colnames TRUE, color colorRampPalette(c(#4B4BDF, white, #DF4B4B))(100), filename 01.GSVA_res.pdf )scale row会对每个基因集在样本间的打分做标准化这样高分数基因集和低分数基因集能在一个色板下对比。cluster_rows TRUE让相似的基因集聚在一起cluster_cols FALSE保留样本的原始排列顺序方便对照分组信息。5.2 plot_data.csv 在可视化中的角色plot_data.csv是绘图脚本读取的中间数据它的列结构通常是基因集名称、样本名、GSVA 分数、分组标签。这种长格式long format比宽矩阵更适合ggplot2绘图如果你要自己改图建议从gsva_res.csv转换成长格式而不是直接操作宽矩阵。library(tidyr) plot_data - gsva_sig %% as.data.frame() %% rownames_to_column(gene_set) %% pivot_longer(-gene_set, names_to sample, values_to score) # 合并分组信息 plot_data - merge(plot_data, group_info, by sample)pivot_longer把宽矩阵的每个样本列拆成sample列和score列再用merge匹配分组这样后面画箱线图、点图或者做相关性分析都不需要重新整理数据。5.3 分组箱线图的快速加画热图展示整体格局但如果你想在论文里展示某个具体通路在两组间的分数差异箱线图是更直接的证据。library(ggplot2) # 选一个代表性基因集 target_set - gsva_sig[1, ] %% names() %% head(1) plot_data %% filter(gene_set target_set) %% ggplot(aes(x group, y score, fill group)) geom_boxplot() geom_jitter(width 0.2) theme_minimal() labs(title target_set, y GSVA score)geom_jitter用于展示单个样本的分布避免样本量少时箱线图掩盖真实数据点。5.4 GSVA_sig_results.csv 的列含义这个文件是在显著性筛选之后导出的完整结果包含基因集名称、分组比较的 p 值、校正后 p 值、以及两组各自的平均 GSVA 分数。它在功能上等价于DESeq2输出的差异基因表只是这里的“基因”换成了“基因集”。gene_set,pvalue,padj,mean_Control,mean_Treatment KEGG_OXIDATIVE_PHOSPHORYLATION,0.001,0.012,-0.21,0.346. 性能瓶颈用并行化处理大型 GMT 集合6.1 c2 全集合的规模问题c2.all.v2023.2.Hs.symbols.gmt包含约一万个基因集如果只做一次 GSVA 打分普通笔记本也能跑完但耗时可能在 30~60 分钟甚至更久取决于样本数和基因数。如果你不需要全集合建议先用c2.cp.kegg或c2.cp.reactome子集跑起来会快很多。6.2 并行版本的 gsva() 调用GSVA包从 1.40 版本开始支持BPPARAM参数可以用BiocParallel开启多核计算。library(BiocParallel) # 开启 4 核并行计算 param - MulticoreParam(workers 4) gsva_res - gsva( expr as.matrix(expr_z), gset.idx.list gs_collection, kcdf Gaussian, method gsva, min.sz 5, max.sz 500, BPPARAM param, verbose TRUE )MulticoreParam(workers 4)会让 4 个核心同时处理不同的基因集子集加速效果在基因集数量过万时尤其明显。注意 Windows 系统不支持MulticoreParam要改用SnowParam(workers 4)。提示如果verbose TRUE且使用并行计算子进程会打印大量日志建议并行时把verbose设为FALSE。6.3 内存不足的优化策略如果你的表达矩阵比较大比如 2 万基因 × 100 样本gsva()会先生成一个 2 万 × 100 的排序矩阵再对每个基因集做数值计算内存峰值可能达到数 GB。这时可以缩小基因集文件或者把表达矩阵按染色体拆分后分段计算最后再合并结果。7. 常见报错与参数调优实战7.1 Error in getGmt: line 1 did not have 2 or 3 columns这个报错说明 GMT 文件格式有问题第一行只有一列或第二列缺失。检查文件是否被 Excel 编辑过Excel 另存的制表符分隔文件可能把\t吞掉了建议用文本编辑器打开确认第二列存在。7.2 Warning: gene sets overlap detection took too longGSVA 在计算前会检查基因集之间的重叠度c2 全集合中大量通路共享基因这一步可能非常耗时。如果只是想做通路活性打分可以在gsva()之前不执行重叠检测直接估算结果。实际上这一步耗时主要来自交集矩阵的构建可以改用min.sz和max.sz过滤后再跑。7.3 表达矩阵与基因集交集为 0 的排查思路如果gsva_res里全是NA或者行数比基因集数少很多先看基因 ID 格式是否一致。data_fpkm.csv的行名一般是Gene Symbol如TP53如果实际行名是ENSG00000141510需要先做 ID 转换。7.4 GSVA 分数全为 0 或分数差异过小的原因GSVA 分数的动态范围依赖基因集内基因的排序波动。如果你的表达矩阵经过过于激进的标准化比如把每个样本都归一化到均值为 0基因排序会被严重压缩导致每个基因集的 ES 都趋近于 0。解决方法是只做 log2 变换不做 z-score让原始表达丰度信息参与到排序中。7.5 小样本量下 p 值校正的极端情况每组 3 个样本时Wilcoxon 检验的最小 p 值是 0.1BH 校正后几乎没有基因集能通过 0.05 阈值。这种情况我会改用 t 检验每组 n≥3 时仍可用或者直接不设 p 值阈值只按 GSVA 分数差异的绝对值排序选出前后 10 个基因集画热图作为探索性分析结果。脚本里的GSVA_sig_results.csv就是考虑了这种场景里面同时保留了未校正的 p 值和分数均值差方便你自己决定筛选标准。本文还有配套的精品资源点击获取 GSVA转录组通路活性分析：原理、R脚本与可视化实战
简介面向零基础转录组学习者的基因集变异分析GSVA配套资源包含输入数据、分析脚本和结果输出三大模块。脚本已提前测试可一键全选跑通适合想快速上手富集分析的生信新手基础薄弱者可结合配套博文教程边学边练逐句理解代码逻辑。资源共十个文件以六个CSV数据表格为核心涵盖表达矩阵、分组信息、富集评分以及显著性筛选结果另含一个R语言脚本、一个GMT基因集注释文件、一份PDF结果报告和一张PNG可视化图片压缩包整体为106.43MB目录划分清晰便于对照输出文件自查。目前已有206人学习使用是完成转录组下游功能富集分析的实用入门资料能帮助读者省去整理数据和调试代码的时间把精力集中在理解生物学意义上。1. 为什么转录组做完差异分析还不够拿到一组 FPKM 表达矩阵差异基因也筛完了KEGG 和 GO 富集也跑了一轮这时候很多人才意识到一个问题差异富集看的是“哪些通路里有差异基因”但它回答不了“这条通路在样本里到底被激活还是被抑制”。GSVAGene Set Variation Analysis就是用来补这块短板的——它不依赖差异基因列表而是直接对全部基因的表达量做无监督打分算出每个基因集在每个样本里的富集分数。这决定了它特别适合处理那些组间差异不剧烈、但整体表达趋势有偏移的数据比如肿瘤分型、药物处理的时间序列、或者临床样本的亚群比较。这套配套资源是围绕一个完整可跑的 R 脚本r.GSVA.R展开的输入是三份文件表达矩阵data_fpkm.csv、样本分组data_group.csv、基因集c2.all.v2023.2.Hs.symbols.gmt输出则包括 GSVA 分数矩阵、显著性筛选结果、以及可直接用于论文的 PDF/PNG 热图。对于刚接触转录组下游分析的人来说最大的门槛其实不是函数调用而是搞不清 GSVA 和 GSEA 的区别、不知道基因集格式怎么处理、以及拿到gsva_res.csv之后下一步该干嘛。这篇文章就按“原理—数据准备—脚本逐段拆解—结果解读—常见报错”这条线把这套资源讲透。提示文中所涉及的代码均从配套r.GSVA.R的完整逻辑中拆出你可以直接复制到 RStudio 里分段执行也可以全选一键跑通。2. GSVA 算法原理与 GSEA 的边界划分2.1 GSVA 的统计思想从表达矩阵到通路活性矩阵GSVA 的核心是把“基因表达量”转换成“基因集富集分数”。它先对每个样本内部的全部基因表达值做核密度估计再把基因按照表达量排序然后用类似于 Kolmogorov-Smirnov 检验的思路计算每个基因集在这个排序中的富集程度。注意这里的“富集”不是看基因是否属于差异基因而是看基因集内的基因整体上是否偏向表达量排序的某一端——偏向高表达端说明该通路在该样本中被激活偏向低表达端则说明被抑制。GSVA_score(sample_j, gene_set_k) max(ES_kj)其中 ES 是基因集内基因在样本 j 表达排序中的富集得分计算时会对每个基因集独立进行所以最终输出的GSVA_result.csv是“基因集 × 样本”的矩阵每一格都是一个连续的数值正负代表激活方向绝对值大小代表强度。这也解释了为什么 GSVA 不需要预先定义“差异基因”它天然是全局分析。2.2 GSVA 与 GSEA 的根本差异GSEA 在做富集分析时输入的是一份带有统计量的基因排序列表通常是 log2FC 或 t 值它比较的是两个生物学状态之间的差异富集而 GSVA 是对每个样本独立打分输出的分数可以直接用于后续的聚类、相关性分析或者生存分析。一句话概括GSEA 是“组间比较”GSVA 是“样本级打分”。维度GSVAGSEA输入核心表达矩阵 基因集基因排序列表 基因集打分粒度每个样本每个基因集每个对比组合每个基因集是否需要分组不需要但分组有利于后续统计必须有两组及以上输出形式样本 × 通路矩阵富集得分及显著性典型应用分型、相关性、生存分析差异通路筛选理解了这个区别你就明白为什么在这套资源里data_group.csv并不是 GSVA 计算的必需输入——它只在后续做组间显著性筛选时才用到。2.3 ssGSEA 与 GSVA 的适用边界配套资源里只用了 GSVA 函数但在实际分析中很多人会把它和 ssGSEAsingle-sample GSEA混用。ssGSEA 同样是单样本打分但它的算法是直接对基因集内的基因做经验 CDF 累积计算复杂度低于 GSVA而且在免疫浸润分析中更常用。如果你手头的数据是芯片而不是 RNA-seq表达谱经过 log2 归一化后GSVA 和 ssGSEA 的结果差异并不大但如果是 FPKM 这种带有长度归一化的数据GSVA 对表达分布的假设更稳所以我一般建议用 GSVA 而不是直接替换成 ssGSEA。3. 输入数据准备GMT 基因集与表达矩阵的格式陷阱3.1 基因集文件 c2.all.v2023.2.Hs.symbols.gmt 的结构GMT 文件每行代表一个基因集三列以上第一列是基因集名称比如KEGG_OXIDATIVE_PHOSPHORYLATION第二列是描述信息通常是一个网址或简要说明第三列开始才是基因符号。这套资源里用的是 MSigDB 的 c2 集合v2023.2 版本包含 KEGG、Reactome、WikiPathways 等来源的通路条目基因标识是Hs.symbols也就是标准的 Human Gene Symbol。KEGG_OXIDATIVE_PHOSPHORYLATION https://www.gsea-msigdb.org/gsea/msigdb/human/geneset/KEGG_OXIDATIVE_PHOSPHORYLATION.html ATP5MC1 ATP5MC2 ATP5MC3 ...注意这里有个坑如果你的表达矩阵里是 Ensembl ID 或 Entrez ID不能直接和这个 GMT 匹配。常见的做法是先做 ID 转换把 Ensembl 转成 Gene Symbol再删掉空匹配和重复行。提示检查表达矩阵的行名和 GMT 里的基因符号是否同一种 ID 格式这一步错了后面所有结果都是空的。3.2 data_fpkm.csv 的标准化要求FPKM 数据的特征是没有经过跨样本标准化基因的表达量受文库大小和基因长度双重影响。GSVA 函数内部会对每个样本的表达向量做排序所以理论上 FPKM 可以直接进 GSVA但对数化通常能带来更稳健的结果。常见的预处理方式是先做log2(FPKM 1)再对每个样本做 z-score 标准化。# 读取表达矩阵行名为基因列名为样本 expr - read.csv(data_fpkm.csv, row.names 1) # 建议的预处理流程先取对数再按行进行标准化 expr_log - log2(expr 1) expr_z - t(apply(expr_log, 1, function(x) { (x - mean(x)) / sd(x) }))这段代码里apply按行MARGIN 1对每个基因在所有样本中的表达量做 z-score 变换这样做的本质是消除不同基因之间的表达丰度差异让 GSVA 排序比较的是“基因内跨样本的相对变化”而不是绝对表达值。对 FPKM 来说这一步非常重要否则那些高表达基因如看家基因会主导排序方向。3.3 data_group.csv 的写入格式data_group.csv用于后续的分组比较最简单的格式是两列第一列是样本名第二列是分组标签。注意样本名必须和表达矩阵的列名完全一致包括顺序和大小写。如果你直接从 Excel 复制出来很容易在列名里带上不可见字符或者把分组列变成因子导致后面compare环节报错。sample,group Ctrl_1,Control Ctrl_2,Control Treat_1,Treatment Treat_2,Treatment4. 核心脚本 r.GSVA.R 逐段拆解参数与执行逻辑4.1 加载 GSVA 包与基因集解析脚本开头加载的关键包是GSVA和GSEABase。GSVA提供gsva()主函数GSEABase负责把 GMT 文件读成标准的基因集对象GeneSetCollection。这一步没有必要手动逐行解析 GMT直接用getGmt()就行但要注意版本差异——老版本GSEABase读取的GeneSet对象里基因类型默认是SYMBOL如果你的表达矩阵行名也是符号那么不需要额外指定geneIdType。library(GSVA) library(GSEABase) # 读取基因集文件 gene_sets - getGmt(c2.all.v2023.2.Hs.symbols.gmt) # 把 GeneSet 列表转换为 GeneSetCollection 对象 gs_collection - GeneSetCollection(gene_sets)这里getGmt()返回的是一个GeneSet列表GeneSetCollection()把它统一成 GSVA 官方接口需要的格式。如果你读到这一步报错说unable to find an inherited method for function getGmt通常是GSEABase没有成功加载而不是文件路径的问题。4.2 调用 gsva() 主函数参数 scRNAseq 与 kcdf 的含义这是整套资源里最重要的一个函数调用# 执行 GSVA 打分 gsva_res - gsva( expr as.matrix(expr_z), gset.idx.list gs_collection, kcdf Gaussian, method gsva, min.sz 5, max.sz 500, verbose TRUE )expr参数要求是数值矩阵不能是数据框所以这里用as.matrix()做显式转换。矩阵的行名必须和基因集里的基因符号能匹配上GSVA 内部会自动取交集没有交集的基因会被丢弃。kcdf参数这是 GSVA 特有的核密度估计类型。对于 RNA-seq 的 FPKM 数据Gaussian是常规选择。如果你的数据是直接的 counts 矩阵可以考虑Poisson但 FPKM 是连续值用 Poisson 反而不合适。method参数gsva表示使用原始 GSVA 算法设置为ssgsea则计算单样本 GSEA。这里用默认的gsva因为我们要的是完整的通路活性矩阵。min.sz和max.sz基因集大小过滤范围。min.sz 5表示基因数少于 5 的基因集直接丢弃max.sz 500表示大于 500 的也会被过滤掉。前者是为了避免小基因集打分不稳定后者是为了避免过于宽泛的通路把信号稀释掉。c2 集合里有些通路非常大比如某些 Reactome 条目有上千个基因这些通路在转录组数据里往往没有生物学特异性所以过滤掉是合理的。4.3 输出 GSVA_result.csv 与 gsva_res.csv 的区别脚本里会生成两个看起来很像的文件GSVA_result.csv和gsva_res.csv。以我拿到资源后的核对来看前者通常是未经后续筛选的完整打分矩阵后者则是经过分组差异比较后选出的显著基因集打分矩阵。简单说GSVA_result.csv是中间产物行是基因集列是样本gsva_res.csv是用于画图的下游输入。# 写出完整结果 write.csv(gsva_res, GSVA_result.csv) # 后续筛选显著基因集后再次写出 write.csv(gsva_sig, gsva_res.csv)4.4 分组显著性判断从 GSVA 分数到差异通路拿到gsva_res矩阵之后脚本会按data_group.csv的分组信息做差异比较。常见做法是对每个基因集的行向量做 t.test 或 Wilcoxon 检验然后筛选 p 值小于阈值比如 0.05的基因集。group_info - read.csv(data_group.csv, row.names 1) group - factor(group_info$group) # 对每个基因集做 Wilcoxon 秩和检验 pvals - apply(gsva_res, 1, function(row) { wilcox.test(row[group levels(group)[1]], row[group levels(group)[2]])$p.value }) # 多重检验校正 padj - p.adjust(pvals, method BH) # 筛选显著基因集 gsva_sig - gsva_res[padj 0.05, , drop FALSE]这里用 Wilcoxon 而不是 t 检验是因为 GSVA 分数不一定是正态分布尤其在小样本量每组 3~5 个样本的转录组数据里秩和检验更稳健。p.adjust的BH方法控制 FDR适合高维假设检验场景。提示如果每组只有 2 个样本wilcox.test不会报错但 p 值没有区分度建议至少每组 3 个样本再跑显著性筛选。5. 可视化输出从 pheatmap 到论文级热图5.1 01.GSVA_res.pdf / 01.GSVA_res.png 的绘图逻辑脚本用pheatmap绘制聚类热图输出同时存 PDF 和 PNG 两个版本PDF 用于论文投稿不失真PNG 用于快速预览。绘图数据就是gsva_sig矩阵按基因集做聚类样本不聚类保持原有的分组顺序会更便于观察组间差异。library(pheatmap) # 绘制显著基因集的热图 pheatmap( gsva_sig, scale row, cluster_rows TRUE, cluster_cols FALSE, show_rownames TRUE, show_colnames TRUE, color colorRampPalette(c(#4B4BDF, white, #DF4B4B))(100), filename 01.GSVA_res.pdf )scale row会对每个基因集在样本间的打分做标准化这样高分数基因集和低分数基因集能在一个色板下对比。cluster_rows TRUE让相似的基因集聚在一起cluster_cols FALSE保留样本的原始排列顺序方便对照分组信息。5.2 plot_data.csv 在可视化中的角色plot_data.csv是绘图脚本读取的中间数据它的列结构通常是基因集名称、样本名、GSVA 分数、分组标签。这种长格式long format比宽矩阵更适合ggplot2绘图如果你要自己改图建议从gsva_res.csv转换成长格式而不是直接操作宽矩阵。library(tidyr) plot_data - gsva_sig %% as.data.frame() %% rownames_to_column(gene_set) %% pivot_longer(-gene_set, names_to sample, values_to score) # 合并分组信息 plot_data - merge(plot_data, group_info, by sample)pivot_longer把宽矩阵的每个样本列拆成sample列和score列再用merge匹配分组这样后面画箱线图、点图或者做相关性分析都不需要重新整理数据。5.3 分组箱线图的快速加画热图展示整体格局但如果你想在论文里展示某个具体通路在两组间的分数差异箱线图是更直接的证据。library(ggplot2) # 选一个代表性基因集 target_set - gsva_sig[1, ] %% names() %% head(1) plot_data %% filter(gene_set target_set) %% ggplot(aes(x group, y score, fill group)) geom_boxplot() geom_jitter(width 0.2) theme_minimal() labs(title target_set, y GSVA score)geom_jitter用于展示单个样本的分布避免样本量少时箱线图掩盖真实数据点。5.4 GSVA_sig_results.csv 的列含义这个文件是在显著性筛选之后导出的完整结果包含基因集名称、分组比较的 p 值、校正后 p 值、以及两组各自的平均 GSVA 分数。它在功能上等价于DESeq2输出的差异基因表只是这里的“基因”换成了“基因集”。gene_set,pvalue,padj,mean_Control,mean_Treatment KEGG_OXIDATIVE_PHOSPHORYLATION,0.001,0.012,-0.21,0.346. 性能瓶颈用并行化处理大型 GMT 集合6.1 c2 全集合的规模问题c2.all.v2023.2.Hs.symbols.gmt包含约一万个基因集如果只做一次 GSVA 打分普通笔记本也能跑完但耗时可能在 30~60 分钟甚至更久取决于样本数和基因数。如果你不需要全集合建议先用c2.cp.kegg或c2.cp.reactome子集跑起来会快很多。6.2 并行版本的 gsva() 调用GSVA包从 1.40 版本开始支持BPPARAM参数可以用BiocParallel开启多核计算。library(BiocParallel) # 开启 4 核并行计算 param - MulticoreParam(workers 4) gsva_res - gsva( expr as.matrix(expr_z), gset.idx.list gs_collection, kcdf Gaussian, method gsva, min.sz 5, max.sz 500, BPPARAM param, verbose TRUE )MulticoreParam(workers 4)会让 4 个核心同时处理不同的基因集子集加速效果在基因集数量过万时尤其明显。注意 Windows 系统不支持MulticoreParam要改用SnowParam(workers 4)。提示如果verbose TRUE且使用并行计算子进程会打印大量日志建议并行时把verbose设为FALSE。6.3 内存不足的优化策略如果你的表达矩阵比较大比如 2 万基因 × 100 样本gsva()会先生成一个 2 万 × 100 的排序矩阵再对每个基因集做数值计算内存峰值可能达到数 GB。这时可以缩小基因集文件或者把表达矩阵按染色体拆分后分段计算最后再合并结果。7. 常见报错与参数调优实战7.1 Error in getGmt: line 1 did not have 2 or 3 columns这个报错说明 GMT 文件格式有问题第一行只有一列或第二列缺失。检查文件是否被 Excel 编辑过Excel 另存的制表符分隔文件可能把\t吞掉了建议用文本编辑器打开确认第二列存在。7.2 Warning: gene sets overlap detection took too longGSVA 在计算前会检查基因集之间的重叠度c2 全集合中大量通路共享基因这一步可能非常耗时。如果只是想做通路活性打分可以在gsva()之前不执行重叠检测直接估算结果。实际上这一步耗时主要来自交集矩阵的构建可以改用min.sz和max.sz过滤后再跑。7.3 表达矩阵与基因集交集为 0 的排查思路如果gsva_res里全是NA或者行数比基因集数少很多先看基因 ID 格式是否一致。data_fpkm.csv的行名一般是Gene Symbol如TP53如果实际行名是ENSG00000141510需要先做 ID 转换。7.4 GSVA 分数全为 0 或分数差异过小的原因GSVA 分数的动态范围依赖基因集内基因的排序波动。如果你的表达矩阵经过过于激进的标准化比如把每个样本都归一化到均值为 0基因排序会被严重压缩导致每个基因集的 ES 都趋近于 0。解决方法是只做 log2 变换不做 z-score让原始表达丰度信息参与到排序中。7.5 小样本量下 p 值校正的极端情况每组 3 个样本时Wilcoxon 检验的最小 p 值是 0.1BH 校正后几乎没有基因集能通过 0.05 阈值。这种情况我会改用 t 检验每组 n≥3 时仍可用或者直接不设 p 值阈值只按 GSVA 分数差异的绝对值排序选出前后 10 个基因集画热图作为探索性分析结果。脚本里的GSVA_sig_results.csv就是考虑了这种场景里面同时保留了未校正的 p 值和分数均值差方便你自己决定筛选标准。本文还有配套的精品资源点击获取