拓冰网站定制 
简介这份资源面向零基础或刚接触转录组下游分析的生物信息学学习者聚焦免疫浸润中的ssGSEA算法实践帮助读者在R语言环境下完成从数据输入到结果可视化的完整流程。压缩包共8个文件约72.61MB包含4个csv数据文件、1个txt注释文件、1个R脚本、1个pdf说明文档和1个png示例图覆盖输入数据、分析代码与输出结果三类内容结构清晰便于按需查阅。其中R脚本已测试可一键全选运行基础较好的读者可直接研读代码逻辑基础薄弱者也能结合配套教程逐步理解算法实现与结果解读。目前已有218人学习下载适合希望快速上手ssGSEA免疫浸润分析、掌握差异比较与可视化输出思路的转录组初学者可作为课程练习或课题预分析的参考材料。1. 零基础也能跑通免疫浸润ssGSEA 到底在算什么转录组测序做完差异分析手里攥着一堆基因名和表达值下一步想回答“这批样本里免疫细胞到底谁多谁少”很多人第一反应是 CIBERSORT。但 CIBERSORT 要跑网页版、要等队列、结果还依赖 LM22 矩阵对零基础的人不算友好。ssGSEA 就不一样了它本质上是单样本基因集富集分析把每个样本单独拎出来算一个免疫细胞特征基因集的富集分数分数越高说明这类细胞在样本里越活跃。配套资源给的是一套 R 语言脚本加示例表达矩阵从读取数据到出免疫浸润热图一条龙。适合刚接触转录组下游、想快速把免疫浸润图做出来的人也适合已经会跑差异分析但卡在可视化这一步的从业者。它不挑物种人和小鼠都能用前提是你手里有基因集文件。2. 把表达矩阵喂给 ssGSEA数据准备与 GSVA 包调用2.1 为什么选 ssGSEA 而不是 CIBERSORTCIBERSORT 做的是解卷积它假设你看到的表达值是多种细胞混合后的结果通过线性方程组反推每种细胞的占比。这个思路很漂亮但有两个硬门槛一是需要 LM22 特征矩阵这个矩阵是固定的换物种或换组织类型就不太灵二是网页版有队列批量样本要等本地版又依赖源码编译。ssGSEA 走的是另一条路它不反推占比而是算富集分数。你把免疫细胞对应的基因集准备好它给每个样本算一个分数这个分数是相对值不是百分比但组间比较完全够用。常见做法是直接从 MSigDB 或文献补充材料里拿免疫细胞基因集整理成 GMT 格式。配套资源里已经带了一份整理好的 GMT省去自己爬文献的功夫。选 ssGSEA 的另一个理由是它对数据分布不敏感。CIBERSORT 要求表达矩阵不能有太多零否则解卷积不稳定。ssGSEA 用的是秩排序零值多也能跑只是分数会偏低。如果你做的是单细胞转录组拟时序或者 bulk RNA-seq 的免疫浸润ssGSEA 的容错率更高。我一般会建议新手先用 ssGSEA 把图跑出来建立信心再去啃 CIBERSORT。2.2 读入表达矩阵与基因集文件配套资源的目录结构很直白一个data文件夹放表达矩阵和基因集一个scripts文件夹放 R 脚本。表达矩阵要求是行名是基因符号列名是样本名值可以是 FPKM、TPM 或 counts但推荐用 TPM因为 TPM 已经做了测序深度和基因长度的校正样本间可比。如果你手里是 FPKM热词里有人问 FPKM 转 TPM 的 R 语言步骤其实就是一个公式先算每个样本的 FPKM 总和再除以总和乘以 1e6。不过 ssGSEA 对 FPKM 和 TPM 的差异不敏感因为它是基于秩的所以不用太纠结。# 加载必要的包 library(GSVA) library(GSEABase) library(limma) # 读取表达矩阵第一列是基因名所以 row.names 1 expr - read.csv(data/expression_matrix.csv, row.names 1, check.names FALSE) # 读取 GMT 基因集文件里面是免疫细胞对应的基因列表 gene_sets - getGmt(data/immune_cell_sets.gmt) # 快速检查一下数据维度 dim(expr) # 预期是基因数 x 样本数 length(gene_sets) # 预期是免疫细胞类型数量这段代码做了三件事读表达矩阵、读 GMT 基因集、检查维度。read.csv里的row.names 1很关键因为表达矩阵第一列通常是基因名不设这个参数 R 会把基因名当成一列数据后面算的时候会报错。check.names FALSE是为了防止 R 把样本名里的特殊字符自动替换成点号比如TCGA-01变成TCGA.01后面画图对不上。getGmt来自 GSEABase 包它把 GMT 文件读成一个 GeneSetCollection 对象GSVA 包能直接识别。如果你没有 GMT 文件也可以用list手动构建但格式要对每个元素是一个字符向量向量里是基因名。2.3 跑 ssGSEA 并理解输出矩阵GSVA 包里的gsva函数是核心method ssgsea指定算法。还有一个参数kcdf需要留意它决定如何把表达值转成累积分布。对于 RNA-seq 的 TPM 或 FPKM推荐kcdf Gaussian因为连续值更接近高斯分布如果是 counts用kcdf Poisson。abs.ranking FALSE是默认值表示保留正负号这样富集分数有正有负正数表示基因集在样本中富集负数表示耗竭。# 运行 ssGSEA ssgsea_scores - gsva( expr as.matrix(expr), gset.idx.list gene_sets, method ssgsea, kcdf Gaussian, abs.ranking FALSE, parallel.sz 4 ) # 查看输出矩阵的前几行几列 ssgsea_scores[1:5, 1:3] # 保存结果方便后续画图 write.csv(ssgsea_scores, results/ssgsea_scores.csv)gsva函数的输出是一个矩阵行是免疫细胞类型列是样本值就是 ssGSEA 分数。parallel.sz 4表示用 4 个核心并行计算样本多的时候能省不少时间。如果你的电脑内存小把parallel.sz设成 1 或 2不然容易爆内存。跑完之后一定要write.csv存下来因为 ssGSEA 计算比较慢尤其是样本数超过 100 的时候存下来下次直接读不用重跑。这里有个血泪经验有一次我跑了 200 个样本没存结果R 会话崩了只能从头再来从那以后我每次跑完都先存盘。3. 免疫浸润结果可视化热图、箱线图与分组比较3.1 用 pheatmap 画免疫细胞热图ssGSEA 分数拿到手第一张图通常是热图行是免疫细胞列是样本颜色深浅代表分数高低。配套资源里用的是pheatmap包因为它比heatmap好看支持列注释能把样本分组标出来。画热图之前要做一件事对分数做标准化。因为不同免疫细胞的分数范围不一样直接画颜色对比不明显。常见做法是按行做 Z-score 标准化也就是每个细胞类型的分数减去该行均值再除以标准差。library(pheatmap) # 按行做 Z-score 标准化 ssgsea_scaled - t(scale(t(ssgsea_scores))) # 读取样本分组信息比如对照组和实验组 group_info - read.csv(data/sample_group.csv, row.names 1) # 画热图 pheatmap( ssgsea_scaled, annotation_col group_info, show_colnames FALSE, cluster_cols TRUE, cluster_rows TRUE, color colorRampPalette(c(navy, white, firebrick3))(100), main Immune Cell Infiltration (ssGSEA) )scale(t(ssgsea_scores))这个写法有点绕先转置让免疫细胞变成列scale默认按列标准化再转置回来就实现了按行标准化。annotation_col接受一个数据框行名是样本名列是分组变量这样热图上方会多一条颜色条标出每个样本属于哪组。show_colnames FALSE是因为样本多的时候列名会挤成一团不如不显示。colorRampPalette定义了一个从深蓝到白到砖红的渐变色比默认配色更符合免疫浸润的直觉红色高蓝色低。如果你想让样本按分组排序而不是聚类把cluster_cols设成 FALSE然后手动指定annotation_col的顺序。3.2 箱线图展示组间差异与统计检验热图看整体模式箱线图看组间差异。把 ssGSEA 分数从矩阵转成长格式用ggplot2画箱线图再叠加上统计检验的显著性标记。配套资源里用的是 Wilcoxon 秩和检验因为免疫浸润分数通常不满足正态分布非参数检验更稳妥。如果你分组超过两组用 Kruskal-Wallis 检验加事后 Dunn 检验。library(ggplot2) library(reshape2) library(ggpubr) # 矩阵转长格式 ssgsea_long - melt(ssgsea_scores, varnames c(CellType, Sample), value.name Score) # 合并分组信息 ssgsea_long$Group - group_info[ssgsea_long$Sample, Group] # 画箱线图加显著性标记 ggplot(ssgsea_long, aes(x CellType, y Score, fill Group)) geom_boxplot(outlier.size 0.5) stat_compare_means(method wilcox.test, label p.signif) theme_bw() theme(axis.text.x element_text(angle 45, hjust 1)) labs(x , y ssGSEA Score, fill Group)melt把宽矩阵变成长表格三列细胞类型、样本、分数。group_info[ssgsea_long$Sample, Group]用样本名去分组表里查对应的组别这一步要求样本名完全一致大小写和连字符都不能差。stat_compare_means来自ggpubr自动在箱线图上加显著性星号p.signif表示用星号代替具体 p 值图更干净。theme(axis.text.x element_text(angle 45, hjust 1))把 x 轴标签旋转 45 度不然免疫细胞名字太长会重叠。如果你觉得箱线图太挤可以只挑差异最显著的几个细胞类型画或者用分面facet_wrap拆成多张小图。3.3 分组比较的统计表与结果导出图是给别人看的表是给自己留底的。跑完统计检验把 p 值、校正后的 FDR、组间均值差整理成一张表方便写文章时直接引用。配套资源里有一个stat_compare脚本输出就是这张表。# 对每个细胞类型做 Wilcoxon 检验 stat_results - sapply(unique(ssgsea_long$CellType), function(ct) { sub_data - subset(ssgsea_long, CellType ct) test - wilcox.test(Score ~ Group, data sub_data) mean_diff - diff(tapply(sub_data$Score, sub_data$Group, mean)) c(p_value test$p.value, mean_diff mean_diff) }) # 转置并校正 p 值 stat_results - as.data.frame(t(stat_results)) stat_results$FDR - p.adjust(stat_results$p_value, method BH) stat_results$CellType - rownames(stat_results) # 导出 write.csv(stat_results, results/immune_stat_results.csv, row.names FALSE)sapply遍历每个细胞类型wilcox.test做检验tapply算两组均值再求差。p.adjust用 Benjamini-Hochberg 方法校正 p 值得到 FDR。这一步很重要因为免疫细胞类型通常有几十种不校正假阳性率会很高。导出的表里有细胞类型、p 值、均值差、FDR按 FDR 排序就能挑出显著差异的细胞类型。我一般会把 FDR 0.05 且均值差绝对值大于 0.1 的细胞类型挑出来作为后续讨论的重点。4. 避坑与排查ssGSEA 跑不动、结果怪、图不对怎么办4.1 报错 “protect(): protection stack overflow”现象跑gsva的时候 R 直接崩提示protect(): protection stack overflow。原因基因集太大或者基因名太多R 的递归保护栈溢出。常见于 GMT 文件里某个基因集包含几千个基因或者表达矩阵基因数超过 2 万。解决在 R 启动时加参数--max-ppsize500000或者用options(expressions 500000)临时调大。如果还不行把基因集拆小每个基因集控制在 500 个基因以内。热词里有人问 “r做转录组t-sne分析时出现错误: protect(): protection stack overflow”其实和 ssGSEA 是同一类问题都是 R 的栈限制调大参数就能解决。4.2 免疫细胞分数全是负数或全是一个值现象跑完 ssGSEA所有样本的分数都是负数或者所有细胞类型的分数几乎一样。原因表达矩阵没有做基因名去重或者基因集里的基因和表达矩阵的基因名不匹配。比如表达矩阵用 Ensembl ID基因集用 Gene Symbol匹配不上ssGSEA 算出来就是随机噪声。解决先检查rownames(expr)和geneIds(gene_sets)的交集数量如果交集很少说明命名系统不一致。用bitr函数做 ID 转换或者手动把 Ensembl ID 转成 Symbol。另外表达矩阵里同一个基因名出现多次gsva会报错用aggregate按基因名取均值去重。4.3 热图聚类把对照组和实验组混在一起现象热图上对照组和实验组的样本没有分开聚类树把两组混在一起。原因免疫浸润分数在两组间差异不大或者样本量太少聚类结果不稳定。解决先看箱线图如果确实没有显著差异热图混在一起是正常的不用强行解释。如果箱线图有差异但热图混可能是标准化方式不对试试按列标准化而不是按行。另外pheatmap默认用欧氏距离对离群样本敏感换成correlation距离或者manhattan距离聚类结果会更稳。4.4 样本名对不上导致分组信息合并失败现象画箱线图时Group列全是 NA或者报错 “undefined columns selected”。原因表达矩阵的列名和分组表的行名不一致比如表达矩阵里是Sample1分组表里是sample1大小写不同。解决用all(colnames(expr) rownames(group_info))检查返回 FALSE 就说明对不上。统一转成大写或小写或者用match函数重新排序。我一般会在读数据之后立刻做一次identical检查不通过就不往下跑省得后面图全错。4.5 ssGSEA 分数和 CIBERSORT 结果对不上现象同一批样本ssGSEA 显示某免疫细胞高CIBERSORT 显示低。原因两种算法的数学基础不同ssGSEA 算的是富集分数CIBERSORT 算的是细胞占比两者没有绝对的可比性。解决不要跨算法比较绝对值只比较组间趋势。如果趋势一致说明结论稳健如果趋势相反检查基因集是否重叠或者用第三种方法如 xCell做验证。常见做法是 ssGSEA 和 CIBERSORT 各跑一遍取交集细胞类型做后续分析这样结论更可靠。5. 进阶技巧把 ssGSEA 分数和临床特征关联起来5.1 生存分析免疫细胞分数与预后ssGSEA 分数不只是画个热图就完了它可以和临床信息结合做生存分析。配套资源里带了一份临床数据包含生存时间和生存状态。把 ssGSEA 分数按中位数分成高分组和低分组用survival包做 Kaplan-Meier 曲线survdiff做 log-rank 检验。如果某个免疫细胞的高分组预后更好说明它可能是保护因素。这一步的关键是选对细胞类型不是所有免疫细胞都和预后相关通常 CD8 T 细胞、NK 细胞、树突状细胞比较有指示意义。library(survival) library(survminer) # 以 CD8 T 细胞为例 cd8_score - ssgsea_scores[CD8_positive_T_cell, ] clinical$CD8_group - ifelse(cd8_score median(cd8_score), High, Low) # 拟合生存曲线 fit - survfit(Surv(OS_time, OS_status) ~ CD8_group, data clinical) # 画图 ggsurvplot(fit, data clinical, pval TRUE, risk.table TRUE, palette c(firebrick3, navy))Surv(OS_time, OS_status)构建生存对象OS_time是生存时间OS_status是生存状态1 表示事件发生0 表示删失。survfit按分组拟合曲线ggsurvplot画图pval TRUE自动加 log-rank p 值。risk.table TRUE在曲线下方加风险人数表这是生存分析的标准配置。如果 p 值小于 0.05说明这个免疫细胞的分组和预后有关。我一般会跑多个细胞类型挑 p 值最小的几个画图不要把所有细胞都画一遍图太多反而没有重点。5.2 免疫检查点基因与 ssGSEA 分数的相关性另一个进阶用法是把 ssGSEA 分数和免疫检查点基因表达做相关性分析。免疫检查点如 PD-1、PD-L1、CTLA-4 是热门靶点如果某个免疫细胞分数和 PD-L1 表达正相关说明这类细胞可能参与了免疫逃逸。做法很简单用cor.test算 Spearman 相关系数再用ggplot2画散点图加拟合线。# 提取 PD-L1 表达值 pdl1_expr - expr[CD274, ] # 算 Spearman 相关 cor_test - cor.test(cd8_score, pdl1_expr, method spearman) # 画散点图 ggplot(data.frame(CD8 cd8_score, PDL1 pdl1_expr), aes(x CD8, y PDL1)) geom_point(alpha 0.6) geom_smooth(method lm, se TRUE, color firebrick3) labs(x CD8 T cell ssGSEA score, y PD-L1 expression) annotate(text, x min(cd8_score), y max(pdl1_expr), label paste0(rho , round(cor_test$estimate, 2), , p , format(cor_test$p.value, digits 2)))cor.test的method spearman指定秩相关因为 ssGSEA 分数和基因表达都不一定正态。geom_smooth(method lm)加线性拟合线和置信区间。annotate把相关系数和 p 值标在图上省得在正文里再写一遍。这里有个细节PD-L1 的基因名是 CD274不是 PD-L1如果你用 PD-L1 去查表达矩阵会找不到。常见做法是查文献确认基因名或者用grep模糊搜索。5.3 批量跑多个基因集并汇总如果你有多个基因集文件比如一个免疫细胞集、一个免疫通路集、一个代谢通路集可以写一个循环批量跑最后把结果合并成一张大表。这样一次分析能覆盖多个维度写文章时素材更丰富。# 定义多个 GMT 文件 gmt_files - c(data/immune_cell_sets.gmt, data/immune_pathway_sets.gmt) # 批量跑 ssGSEA all_scores - lapply(gmt_files, function(f) { gs - getGmt(f) gsva(as.matrix(expr), gs, method ssgsea, kcdf Gaussian) }) # 合并结果 combined_scores - do.call(rbind, all_scores) write.csv(combined_scores, results/combined_ssgsea_scores.csv)lapply遍历每个 GMT 文件do.call(rbind, ...)把结果按行合并。注意不同基因集的细胞类型或通路名不能重复否则合并后行名会冲突。我一般会在合并前给行名加前缀比如Immune_和Pathway_这样一眼就能看出分数来自哪个基因集。跑完这个流程你手里就有一张覆盖多个维度的 ssGSEA 分数表后续画热图、做生存分析、算相关性都够用了。从那以后我每次跑 ssGSEA 都强制走一遍检查先看基因名交集再看样本名是否对齐最后存盘再画图。这三步花不了五分钟但能省下几个小时的返工。希望帮到你。本文还有配套的精品资源点击获取 零基础跑通ssGSEA免疫浸润：从表达矩阵到热图与生存分析
简介这份资源面向零基础或刚接触转录组下游分析的生物信息学学习者聚焦免疫浸润中的ssGSEA算法实践帮助读者在R语言环境下完成从数据输入到结果可视化的完整流程。压缩包共8个文件约72.61MB包含4个csv数据文件、1个txt注释文件、1个R脚本、1个pdf说明文档和1个png示例图覆盖输入数据、分析代码与输出结果三类内容结构清晰便于按需查阅。其中R脚本已测试可一键全选运行基础较好的读者可直接研读代码逻辑基础薄弱者也能结合配套教程逐步理解算法实现与结果解读。目前已有218人学习下载适合希望快速上手ssGSEA免疫浸润分析、掌握差异比较与可视化输出思路的转录组初学者可作为课程练习或课题预分析的参考材料。1. 零基础也能跑通免疫浸润ssGSEA 到底在算什么转录组测序做完差异分析手里攥着一堆基因名和表达值下一步想回答“这批样本里免疫细胞到底谁多谁少”很多人第一反应是 CIBERSORT。但 CIBERSORT 要跑网页版、要等队列、结果还依赖 LM22 矩阵对零基础的人不算友好。ssGSEA 就不一样了它本质上是单样本基因集富集分析把每个样本单独拎出来算一个免疫细胞特征基因集的富集分数分数越高说明这类细胞在样本里越活跃。配套资源给的是一套 R 语言脚本加示例表达矩阵从读取数据到出免疫浸润热图一条龙。适合刚接触转录组下游、想快速把免疫浸润图做出来的人也适合已经会跑差异分析但卡在可视化这一步的从业者。它不挑物种人和小鼠都能用前提是你手里有基因集文件。2. 把表达矩阵喂给 ssGSEA数据准备与 GSVA 包调用2.1 为什么选 ssGSEA 而不是 CIBERSORTCIBERSORT 做的是解卷积它假设你看到的表达值是多种细胞混合后的结果通过线性方程组反推每种细胞的占比。这个思路很漂亮但有两个硬门槛一是需要 LM22 特征矩阵这个矩阵是固定的换物种或换组织类型就不太灵二是网页版有队列批量样本要等本地版又依赖源码编译。ssGSEA 走的是另一条路它不反推占比而是算富集分数。你把免疫细胞对应的基因集准备好它给每个样本算一个分数这个分数是相对值不是百分比但组间比较完全够用。常见做法是直接从 MSigDB 或文献补充材料里拿免疫细胞基因集整理成 GMT 格式。配套资源里已经带了一份整理好的 GMT省去自己爬文献的功夫。选 ssGSEA 的另一个理由是它对数据分布不敏感。CIBERSORT 要求表达矩阵不能有太多零否则解卷积不稳定。ssGSEA 用的是秩排序零值多也能跑只是分数会偏低。如果你做的是单细胞转录组拟时序或者 bulk RNA-seq 的免疫浸润ssGSEA 的容错率更高。我一般会建议新手先用 ssGSEA 把图跑出来建立信心再去啃 CIBERSORT。2.2 读入表达矩阵与基因集文件配套资源的目录结构很直白一个data文件夹放表达矩阵和基因集一个scripts文件夹放 R 脚本。表达矩阵要求是行名是基因符号列名是样本名值可以是 FPKM、TPM 或 counts但推荐用 TPM因为 TPM 已经做了测序深度和基因长度的校正样本间可比。如果你手里是 FPKM热词里有人问 FPKM 转 TPM 的 R 语言步骤其实就是一个公式先算每个样本的 FPKM 总和再除以总和乘以 1e6。不过 ssGSEA 对 FPKM 和 TPM 的差异不敏感因为它是基于秩的所以不用太纠结。# 加载必要的包 library(GSVA) library(GSEABase) library(limma) # 读取表达矩阵第一列是基因名所以 row.names 1 expr - read.csv(data/expression_matrix.csv, row.names 1, check.names FALSE) # 读取 GMT 基因集文件里面是免疫细胞对应的基因列表 gene_sets - getGmt(data/immune_cell_sets.gmt) # 快速检查一下数据维度 dim(expr) # 预期是基因数 x 样本数 length(gene_sets) # 预期是免疫细胞类型数量这段代码做了三件事读表达矩阵、读 GMT 基因集、检查维度。read.csv里的row.names 1很关键因为表达矩阵第一列通常是基因名不设这个参数 R 会把基因名当成一列数据后面算的时候会报错。check.names FALSE是为了防止 R 把样本名里的特殊字符自动替换成点号比如TCGA-01变成TCGA.01后面画图对不上。getGmt来自 GSEABase 包它把 GMT 文件读成一个 GeneSetCollection 对象GSVA 包能直接识别。如果你没有 GMT 文件也可以用list手动构建但格式要对每个元素是一个字符向量向量里是基因名。2.3 跑 ssGSEA 并理解输出矩阵GSVA 包里的gsva函数是核心method ssgsea指定算法。还有一个参数kcdf需要留意它决定如何把表达值转成累积分布。对于 RNA-seq 的 TPM 或 FPKM推荐kcdf Gaussian因为连续值更接近高斯分布如果是 counts用kcdf Poisson。abs.ranking FALSE是默认值表示保留正负号这样富集分数有正有负正数表示基因集在样本中富集负数表示耗竭。# 运行 ssGSEA ssgsea_scores - gsva( expr as.matrix(expr), gset.idx.list gene_sets, method ssgsea, kcdf Gaussian, abs.ranking FALSE, parallel.sz 4 ) # 查看输出矩阵的前几行几列 ssgsea_scores[1:5, 1:3] # 保存结果方便后续画图 write.csv(ssgsea_scores, results/ssgsea_scores.csv)gsva函数的输出是一个矩阵行是免疫细胞类型列是样本值就是 ssGSEA 分数。parallel.sz 4表示用 4 个核心并行计算样本多的时候能省不少时间。如果你的电脑内存小把parallel.sz设成 1 或 2不然容易爆内存。跑完之后一定要write.csv存下来因为 ssGSEA 计算比较慢尤其是样本数超过 100 的时候存下来下次直接读不用重跑。这里有个血泪经验有一次我跑了 200 个样本没存结果R 会话崩了只能从头再来从那以后我每次跑完都先存盘。3. 免疫浸润结果可视化热图、箱线图与分组比较3.1 用 pheatmap 画免疫细胞热图ssGSEA 分数拿到手第一张图通常是热图行是免疫细胞列是样本颜色深浅代表分数高低。配套资源里用的是pheatmap包因为它比heatmap好看支持列注释能把样本分组标出来。画热图之前要做一件事对分数做标准化。因为不同免疫细胞的分数范围不一样直接画颜色对比不明显。常见做法是按行做 Z-score 标准化也就是每个细胞类型的分数减去该行均值再除以标准差。library(pheatmap) # 按行做 Z-score 标准化 ssgsea_scaled - t(scale(t(ssgsea_scores))) # 读取样本分组信息比如对照组和实验组 group_info - read.csv(data/sample_group.csv, row.names 1) # 画热图 pheatmap( ssgsea_scaled, annotation_col group_info, show_colnames FALSE, cluster_cols TRUE, cluster_rows TRUE, color colorRampPalette(c(navy, white, firebrick3))(100), main Immune Cell Infiltration (ssGSEA) )scale(t(ssgsea_scores))这个写法有点绕先转置让免疫细胞变成列scale默认按列标准化再转置回来就实现了按行标准化。annotation_col接受一个数据框行名是样本名列是分组变量这样热图上方会多一条颜色条标出每个样本属于哪组。show_colnames FALSE是因为样本多的时候列名会挤成一团不如不显示。colorRampPalette定义了一个从深蓝到白到砖红的渐变色比默认配色更符合免疫浸润的直觉红色高蓝色低。如果你想让样本按分组排序而不是聚类把cluster_cols设成 FALSE然后手动指定annotation_col的顺序。3.2 箱线图展示组间差异与统计检验热图看整体模式箱线图看组间差异。把 ssGSEA 分数从矩阵转成长格式用ggplot2画箱线图再叠加上统计检验的显著性标记。配套资源里用的是 Wilcoxon 秩和检验因为免疫浸润分数通常不满足正态分布非参数检验更稳妥。如果你分组超过两组用 Kruskal-Wallis 检验加事后 Dunn 检验。library(ggplot2) library(reshape2) library(ggpubr) # 矩阵转长格式 ssgsea_long - melt(ssgsea_scores, varnames c(CellType, Sample), value.name Score) # 合并分组信息 ssgsea_long$Group - group_info[ssgsea_long$Sample, Group] # 画箱线图加显著性标记 ggplot(ssgsea_long, aes(x CellType, y Score, fill Group)) geom_boxplot(outlier.size 0.5) stat_compare_means(method wilcox.test, label p.signif) theme_bw() theme(axis.text.x element_text(angle 45, hjust 1)) labs(x , y ssGSEA Score, fill Group)melt把宽矩阵变成长表格三列细胞类型、样本、分数。group_info[ssgsea_long$Sample, Group]用样本名去分组表里查对应的组别这一步要求样本名完全一致大小写和连字符都不能差。stat_compare_means来自ggpubr自动在箱线图上加显著性星号p.signif表示用星号代替具体 p 值图更干净。theme(axis.text.x element_text(angle 45, hjust 1))把 x 轴标签旋转 45 度不然免疫细胞名字太长会重叠。如果你觉得箱线图太挤可以只挑差异最显著的几个细胞类型画或者用分面facet_wrap拆成多张小图。3.3 分组比较的统计表与结果导出图是给别人看的表是给自己留底的。跑完统计检验把 p 值、校正后的 FDR、组间均值差整理成一张表方便写文章时直接引用。配套资源里有一个stat_compare脚本输出就是这张表。# 对每个细胞类型做 Wilcoxon 检验 stat_results - sapply(unique(ssgsea_long$CellType), function(ct) { sub_data - subset(ssgsea_long, CellType ct) test - wilcox.test(Score ~ Group, data sub_data) mean_diff - diff(tapply(sub_data$Score, sub_data$Group, mean)) c(p_value test$p.value, mean_diff mean_diff) }) # 转置并校正 p 值 stat_results - as.data.frame(t(stat_results)) stat_results$FDR - p.adjust(stat_results$p_value, method BH) stat_results$CellType - rownames(stat_results) # 导出 write.csv(stat_results, results/immune_stat_results.csv, row.names FALSE)sapply遍历每个细胞类型wilcox.test做检验tapply算两组均值再求差。p.adjust用 Benjamini-Hochberg 方法校正 p 值得到 FDR。这一步很重要因为免疫细胞类型通常有几十种不校正假阳性率会很高。导出的表里有细胞类型、p 值、均值差、FDR按 FDR 排序就能挑出显著差异的细胞类型。我一般会把 FDR 0.05 且均值差绝对值大于 0.1 的细胞类型挑出来作为后续讨论的重点。4. 避坑与排查ssGSEA 跑不动、结果怪、图不对怎么办4.1 报错 “protect(): protection stack overflow”现象跑gsva的时候 R 直接崩提示protect(): protection stack overflow。原因基因集太大或者基因名太多R 的递归保护栈溢出。常见于 GMT 文件里某个基因集包含几千个基因或者表达矩阵基因数超过 2 万。解决在 R 启动时加参数--max-ppsize500000或者用options(expressions 500000)临时调大。如果还不行把基因集拆小每个基因集控制在 500 个基因以内。热词里有人问 “r做转录组t-sne分析时出现错误: protect(): protection stack overflow”其实和 ssGSEA 是同一类问题都是 R 的栈限制调大参数就能解决。4.2 免疫细胞分数全是负数或全是一个值现象跑完 ssGSEA所有样本的分数都是负数或者所有细胞类型的分数几乎一样。原因表达矩阵没有做基因名去重或者基因集里的基因和表达矩阵的基因名不匹配。比如表达矩阵用 Ensembl ID基因集用 Gene Symbol匹配不上ssGSEA 算出来就是随机噪声。解决先检查rownames(expr)和geneIds(gene_sets)的交集数量如果交集很少说明命名系统不一致。用bitr函数做 ID 转换或者手动把 Ensembl ID 转成 Symbol。另外表达矩阵里同一个基因名出现多次gsva会报错用aggregate按基因名取均值去重。4.3 热图聚类把对照组和实验组混在一起现象热图上对照组和实验组的样本没有分开聚类树把两组混在一起。原因免疫浸润分数在两组间差异不大或者样本量太少聚类结果不稳定。解决先看箱线图如果确实没有显著差异热图混在一起是正常的不用强行解释。如果箱线图有差异但热图混可能是标准化方式不对试试按列标准化而不是按行。另外pheatmap默认用欧氏距离对离群样本敏感换成correlation距离或者manhattan距离聚类结果会更稳。4.4 样本名对不上导致分组信息合并失败现象画箱线图时Group列全是 NA或者报错 “undefined columns selected”。原因表达矩阵的列名和分组表的行名不一致比如表达矩阵里是Sample1分组表里是sample1大小写不同。解决用all(colnames(expr) rownames(group_info))检查返回 FALSE 就说明对不上。统一转成大写或小写或者用match函数重新排序。我一般会在读数据之后立刻做一次identical检查不通过就不往下跑省得后面图全错。4.5 ssGSEA 分数和 CIBERSORT 结果对不上现象同一批样本ssGSEA 显示某免疫细胞高CIBERSORT 显示低。原因两种算法的数学基础不同ssGSEA 算的是富集分数CIBERSORT 算的是细胞占比两者没有绝对的可比性。解决不要跨算法比较绝对值只比较组间趋势。如果趋势一致说明结论稳健如果趋势相反检查基因集是否重叠或者用第三种方法如 xCell做验证。常见做法是 ssGSEA 和 CIBERSORT 各跑一遍取交集细胞类型做后续分析这样结论更可靠。5. 进阶技巧把 ssGSEA 分数和临床特征关联起来5.1 生存分析免疫细胞分数与预后ssGSEA 分数不只是画个热图就完了它可以和临床信息结合做生存分析。配套资源里带了一份临床数据包含生存时间和生存状态。把 ssGSEA 分数按中位数分成高分组和低分组用survival包做 Kaplan-Meier 曲线survdiff做 log-rank 检验。如果某个免疫细胞的高分组预后更好说明它可能是保护因素。这一步的关键是选对细胞类型不是所有免疫细胞都和预后相关通常 CD8 T 细胞、NK 细胞、树突状细胞比较有指示意义。library(survival) library(survminer) # 以 CD8 T 细胞为例 cd8_score - ssgsea_scores[CD8_positive_T_cell, ] clinical$CD8_group - ifelse(cd8_score median(cd8_score), High, Low) # 拟合生存曲线 fit - survfit(Surv(OS_time, OS_status) ~ CD8_group, data clinical) # 画图 ggsurvplot(fit, data clinical, pval TRUE, risk.table TRUE, palette c(firebrick3, navy))Surv(OS_time, OS_status)构建生存对象OS_time是生存时间OS_status是生存状态1 表示事件发生0 表示删失。survfit按分组拟合曲线ggsurvplot画图pval TRUE自动加 log-rank p 值。risk.table TRUE在曲线下方加风险人数表这是生存分析的标准配置。如果 p 值小于 0.05说明这个免疫细胞的分组和预后有关。我一般会跑多个细胞类型挑 p 值最小的几个画图不要把所有细胞都画一遍图太多反而没有重点。5.2 免疫检查点基因与 ssGSEA 分数的相关性另一个进阶用法是把 ssGSEA 分数和免疫检查点基因表达做相关性分析。免疫检查点如 PD-1、PD-L1、CTLA-4 是热门靶点如果某个免疫细胞分数和 PD-L1 表达正相关说明这类细胞可能参与了免疫逃逸。做法很简单用cor.test算 Spearman 相关系数再用ggplot2画散点图加拟合线。# 提取 PD-L1 表达值 pdl1_expr - expr[CD274, ] # 算 Spearman 相关 cor_test - cor.test(cd8_score, pdl1_expr, method spearman) # 画散点图 ggplot(data.frame(CD8 cd8_score, PDL1 pdl1_expr), aes(x CD8, y PDL1)) geom_point(alpha 0.6) geom_smooth(method lm, se TRUE, color firebrick3) labs(x CD8 T cell ssGSEA score, y PD-L1 expression) annotate(text, x min(cd8_score), y max(pdl1_expr), label paste0(rho , round(cor_test$estimate, 2), , p , format(cor_test$p.value, digits 2)))cor.test的method spearman指定秩相关因为 ssGSEA 分数和基因表达都不一定正态。geom_smooth(method lm)加线性拟合线和置信区间。annotate把相关系数和 p 值标在图上省得在正文里再写一遍。这里有个细节PD-L1 的基因名是 CD274不是 PD-L1如果你用 PD-L1 去查表达矩阵会找不到。常见做法是查文献确认基因名或者用grep模糊搜索。5.3 批量跑多个基因集并汇总如果你有多个基因集文件比如一个免疫细胞集、一个免疫通路集、一个代谢通路集可以写一个循环批量跑最后把结果合并成一张大表。这样一次分析能覆盖多个维度写文章时素材更丰富。# 定义多个 GMT 文件 gmt_files - c(data/immune_cell_sets.gmt, data/immune_pathway_sets.gmt) # 批量跑 ssGSEA all_scores - lapply(gmt_files, function(f) { gs - getGmt(f) gsva(as.matrix(expr), gs, method ssgsea, kcdf Gaussian) }) # 合并结果 combined_scores - do.call(rbind, all_scores) write.csv(combined_scores, results/combined_ssgsea_scores.csv)lapply遍历每个 GMT 文件do.call(rbind, ...)把结果按行合并。注意不同基因集的细胞类型或通路名不能重复否则合并后行名会冲突。我一般会在合并前给行名加前缀比如Immune_和Pathway_这样一眼就能看出分数来自哪个基因集。跑完这个流程你手里就有一张覆盖多个维度的 ssGSEA 分数表后续画热图、做生存分析、算相关性都够用了。从那以后我每次跑 ssGSEA 都强制走一遍检查先看基因名交集再看样本名是否对齐最后存盘再画图。这三步花不了五分钟但能省下几个小时的返工。希望帮到你。本文还有配套的精品资源点击获取