拓十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

零基础跑通ssGSEA免疫浸润分析:从FPKM转TPM到生存分析全流程

零基础跑通ssGSEA免疫浸润分析:从FPKM转TPM到生存分析全流程 简介这份资源面向零基础或刚接触转录组下游分析的生物信息学学习者聚焦免疫浸润中的ssGSEA算法实践帮助读者在R语言环境下完成从数据输入到结果可视化的完整流程。压缩包共8个文件包含4个csv数据表、1个txt注释文件、1个R脚本、1个pdf说明文档和1个png示例图整体约72.61MB覆盖输入数据、可一键运行的R脚本与输出结果三部分脚本已测试可全选后直接跑通。已有218人学习下载适合基础薄弱者搭配教程逐步理解也方便有经验者直接研读代码。通过这份资源读者能获得完整的免疫浸润分析素材、可复用的ssGSEA脚本、差异比较与Wilcoxon检验结果示例以及配套图表输出便于快速复现流程并迁移到自己的转录组项目中。1. 零基础也能跑通免疫浸润ssGSEA 到底在算什么你手里有一份转录组表达矩阵可能是肿瘤 vs 癌旁也可能是处理组 vs 对照组差异分析做完了富集也跑完了但审稿人或者老板追问一句「免疫细胞浸润情况怎么样」很多人就卡住了。ssGSEAsingle-sample Gene Set Enrichment Analysis就是解决这个问题的常用手段它不需要你有原始测序数据只要有一个基因表达矩阵和一套免疫细胞标记基因集就能给每个样本算出一组免疫细胞浸润评分。和 CIBERSORT 这类需要解卷积的方法比ssGSEA 不要求绝对定量、对批次不敏感、零基础也能在 R 里几行代码跑通所以成了转录组下游免疫浸润分析里出镜率最高的算法之一。这篇笔记就按「准备数据 → 跑 ssGSEA → 出图 → 排错」的顺序把每一步的参数和坑讲清楚让你拿到自己的矩阵就能复现。2. 跑 ssGSEA 之前表达矩阵和基因集的准备2.1 表达矩阵要满足的三个硬条件ssGSEA 对输入矩阵其实不挑但有几个前提不满足结果会直接失真。第一矩阵必须是基因 × 样本行名是基因符号SYMBOL不是 Ensembl ID也不是探针 ID第二值必须是归一化后的表达量FPKM、TPM、CPM 或者 log2 转换后的都可以但绝不能直接塞原始 counts因为 ssGSEA 内部做的是秩排序counts 的尺度差异会让排序结果偏离真实生物学信号第三矩阵里不能有大量 NA 和全零行免疫基因本身表达量偏低如果过滤太狠标记基因被删掉评分就会系统性偏低。这里顺带说一个高频热搜问题转录组测序 FPKM 值换算成 TPM 的 R 语言步骤。很多人手里只有 FPKM想换成 TPM 再跑 ssGSEA。换算逻辑是先把 FPKM 按基因长度还原成 reads 数再按样本总 reads 归一化# fpkm_matrix: 行是基因列是样本值是对应 FPKM # gene_length: 命名向量名字与 rownames(fpkm_matrix) 对应单位是 bp fpkm_to_tpm - function(fpkm_matrix, gene_length) { # 第一步FPKM 反推 reads count除以基因长度再乘 1e3 rate - fpkm_matrix / gene_length[rownames(fpkm_matrix)] * 1e3 # 第二步每个样本的 rate 求和作为该样本的归一化因子 scaling - colSums(rate, na.rm TRUE) # 第三步rate 除以 scaling 再乘 1e6得到 TPM tpm - t(t(rate) / scaling) * 1e6 return(tpm) }逻辑说明FPKM 和 TPM 的差别就在归一化顺序FPKM 先按总 reads 归一化再除基因长度TPM 先除基因长度再按总和归一化所以换算必须知道基因长度。参数上gene_length的单位要和 FPKM 计算时一致一般用 exon 长度之和如果长度缺失可以从 GTF 文件用GenomicFeatures提取。跑完 ssGSEA 前建议对 TPM 做一次log2(tpm 1)让分布更接近正态后续做组间比较时 t 检验更稳。2.2 免疫细胞基因集从哪来怎么选ssGSEA 的核心是基因集。免疫浸润最常用的是Charoentong 等 2017 年发表的 28 种免疫细胞标记集覆盖了活化 CD8 T 细胞、巨噬细胞、MDSC、Treg 等做肿瘤免疫微环境基本够用。另一个常见来源是Bindea 等 2013 的免疫细胞标记细胞类型更细但部分基因集只有十几个基因稳定性差一些。选基因集的原则是基因数少于 10 个的慎用因为 ssGSEA 的富集评分对基因集大小敏感同一批分析里不要混用两套基因集否则评分不可比。拿到基因集后格式整理成 list每个元素是一个细胞类型对应的基因符号向量# 假设 gene_set 是 data.frame两列cell_type 和 gene_symbol library(dplyr) immune_sets - gene_set %% group_by(cell_type) %% summarise(genes list(unique(gene_symbol))) %% tibble::deframe() # 检查每个基因集大小过滤掉过小的 immune_sets - immune_sets[sapply(immune_sets, length) 10]逻辑说明deframe把两列 data.frame 转成命名 list名字是细胞类型值是基因向量这正是GSVA::gsva()需要的格式。参数上length 10是经验阈值低于这个数评分波动大如果你的矩阵基因数本来就少可以放宽到 5但要在文章里说明。3. 用 GSVA 包跑 ssGSEA参数怎么设、结果怎么读3.1 最小可运行命令与四个关键参数R 里跑 ssGSEA 最稳的是GSVA包底层调用的是单样本 GSEA。最小命令如下library(GSVA) library(GSEABase) # expr: 基因 × 样本的表达矩阵已 log2 转换 # immune_sets: 上一步得到的 list ssgsea_score - gsva( expr as.matrix(expr), gset.idx.list immune_sets, method ssgsea, kcdf Gaussian, abs.ranking FALSE, min.sz 10, max.sz 500, parallel.sz 4, mx.diff TRUE, tau 0.25, ssgsea.norm TRUE, verbose TRUE )逻辑说明method ssgsea指定算法kcdf Gaussian表示输入是连续表达值log2 后的 TPM/FPKM 用这个如果输入是 counts 则用Poissonabs.ranking FALSE保留基因排序方向免疫浸润分析必须保留tau 0.25是 ssGSEA 的权重参数默认 0.25调大更强调高表达基因一般不动ssgsea.norm TRUE会对每个样本的评分做归一化让不同样本可比这个一定要开。parallel.sz按你机器核数设4 到 8 都行设太大反而因为内存拷贝变慢。跑完得到一个细胞类型 × 样本的矩阵值就是每个样本每种免疫细胞的富集评分。评分本身没有绝对意义只有样本间相对高低和组间差异有意义。3.2 结果矩阵的标准化与分组比较原始 ssGSEA 评分在不同细胞类型之间尺度不同做热图前一般按行做 Z-score 标准化# ssgsea_score: 细胞类型 × 样本 z_score - t(scale(t(ssgsea_score))) # 分组比较group 是长度等于样本数的因子向量 library(limma) group - factor(group, levels c(Control, Treat)) design - model.matrix(~ 0 group) colnames(design) - levels(group) fit - lmFit(ssgsea_score, design) contrast - makeContrasts(Treat - Control, levels design) fit2 - contrasts.fit(fit, contrast) fit2 - eBayes(fit2) diff_immune - topTable(fit2, number Inf, adjust.method BH, sort.by P)逻辑说明scale(t(...))先转置再按行标准化因为scale默认按列操作转置两次就实现了按细胞类型标准化。limma 这里用的是评分矩阵直接做线性模型比逐细胞类型跑 t 检验更稳能借用整体方差信息。参数上adjust.method BH是 FDR 校正免疫浸润分析里细胞类型多必须校正sort.by P按 P 值排序方便挑显著差异的细胞类型。注意ssGSEA 评分做组间比较时如果两组样本量都小于 5limma 的方差估计会不稳这时候建议改用 Wilcoxon 秩和检验并在结果里标注检验方法。4. 出图与验证热图、箱线图和 GSVA 打分的关系4.1 差异免疫细胞热图与箱线图热图用pheatmap输入 Z-score 矩阵按分组加注释条library(pheatmap) annotation_col - data.frame(Group group) rownames(annotation_col) - colnames(z_score) pheatmap( z_score, annotation_col annotation_col, cluster_rows TRUE, cluster_cols TRUE, show_colnames FALSE, color colorRampPalette(c(navy, white, firebrick3))(100), border_color NA, fontsize_row 9, main Immune infiltration (ssGSEA Z-score) )逻辑说明cluster_cols TRUE让样本按免疫浸润模式聚类常常能看出分组结构show_colnames FALSE在样本多的时候避免标签重叠。颜色用 navy-white-firebrick3 是免疫浸润文章里最常见的配色低浸润蓝、高浸润红审稿人一眼能看懂。箱线图挑显著差异的细胞类型画library(ggplot2) library(reshape2) diff_cells - rownames(diff_immune)[diff_immune$P.Value 0.05] plot_data - melt(ssgsea_score[diff_cells, , drop FALSE]) colnames(plot_data) - c(CellType, Sample, Score) plot_data$Group - group[match(plot_data$Sample, colnames(ssgsea_score))] ggplot(plot_data, aes(x CellType, y Score, fill Group)) geom_boxplot(outlier.size 0.5) theme_bw() theme(axis.text.x element_text(angle 45, hjust 1)) labs(y ssGSEA score, x )逻辑说明melt把宽矩阵转成长表match按样本名对齐分组避免顺序错乱。参数上outlier.size调小是因为免疫评分偶尔有极端值点太大会盖住箱体。4.2 用 GSVA 打分和免疫检查点基因做相关性验证单跑 ssGSEA 容易被质疑「评分是不是真的反映免疫浸润」一个常见验证是把关键免疫检查点基因PDCD1、CD274、CTLA4、LAG3的表达和对应细胞评分做相关性checkpoint_genes - c(PDCD1, CD274, CTLA4, LAG3) checkpoint_expr - expr[intersect(checkpoint_genes, rownames(expr)), , drop FALSE] cor_result - sapply(rownames(ssgsea_score), function(cell) { sapply(rownames(checkpoint_expr), function(g) { cor(ssgsea_score[cell, ], checkpoint_expr[g, ], method spearman) }) })逻辑说明Spearman 相关对非线性单调关系也稳免疫评分和检查点表达通常不是严格线性。参数上method spearman比 Pearson 更适合评分数据结果矩阵行是检查点基因、列是细胞类型正相关说明评分方向合理。如果出现大面积负相关先检查表达矩阵是否做了 log 转换、基因集方向是否搞反。5. 避坑与排查ssGSEA 跑崩、结果反常的 5 个现场现象一gsva()报错protect(): protection stack overflow。原因通常是基因集 list 嵌套太深或者矩阵太大R 的递归保护栈溢出。解决在 R 启动前执行ulimit -s 65536Linux/Mac或者在 R 里options(expressions 500000)再把parallel.sz降到 1 先跑通确认不是并行导致的内存问题。现象二所有样本所有细胞评分几乎一样热图一片白。原因多半是输入矩阵没有做 log 转换或者基因名是 Ensembl ID 而基因集是 SYMBOL匹配率极低。解决检查sum(rownames(expr) %in% unlist(immune_sets))如果低于基因集总基因数的 30%说明 ID 类型不匹配需要用clusterProfiler::bitr或org.Hs.eg.db转换。现象三评分能跑出来但组间比较全不显著。原因可能是样本量太小、或者分组本身免疫浸润差异就小。解决先看 PCA 或热图聚类确认两组在整体表达上是否分开如果整体就分不开免疫浸润不显著是正常结果不要硬调参数凑显著性。现象四kcdf参数设错导致评分异常。输入是 log2 TPM 却设了kcdf Poisson评分会整体偏移。解决连续值一律Gaussian原始 counts 才用Poisson不确定就先跑一小部分样本对比两种设置的评分相关性低于 0.9 说明设错了。现象五换一套基因集后结论反转。不同来源的免疫细胞标记基因集覆盖的基因不同评分不可直接跨研究比较。解决同一篇文章里固定一套基因集方法部分写清楚来源和版本如果要做敏感性分析可以两套都跑但结论只以主分析为准另一套放补充材料。6. 进阶技巧把 ssGSEA 评分接到生存分析和亚型分型上ssGSEA 评分真正的价值不在热图本身而在于它能作为连续变量接下游分析。一个我常用的套路是挑出差异显著的免疫细胞评分按中位数把样本分成高/低浸润两组做 Kaplan-Meier 生存曲线。这里有个细节连续变量按中位数二分会损失信息更稳的做法是用survminer::surv_cutpoint找最优截断值library(survival) library(survminer) # surv_data: 包含 OS.time, OS.status, 以及某个细胞评分 surv_data$score - ssgsea_score[Activated.CD8.T.cell, ] cut - surv_cutpoint( surv_data, time OS.time, event OS.status, variables score, minprop 0.2 ) surv_data$group - surv_cutgroup(surv_data, cut, variables score) fit - survfit(Surv(OS.time, OS.status) ~ group, data surv_data) ggsurvplot(fit, data surv_data, pval TRUE, risk.table TRUE)逻辑说明surv_cutpoint通过最大化 log-rank 统计量找截断值minprop 0.2保证每组至少 20% 样本避免切出极端分组。参数上variables可以一次传多个细胞评分函数会分别找截断值。这一步做完ssGSEA 就从一张热图变成了有临床意义的结论。另一个进阶方向是把多个免疫细胞评分做无监督聚类识别免疫亚型。我一般用ConsensusClusterPlus对 Z-score 后的评分矩阵聚类k 从 2 试到 6看 CDF 曲线拐点定亚型数。这里最容易翻车的地方是聚类前一定要确认评分矩阵没有批次效应否则聚出来的亚型可能只是批次差异。我自己的习惯是跑完 ssGSEA 先存一份原始评分矩阵和一份 Z-score 矩阵所有下游分析都从这两份文件出发中间不再重新跑gsva()这样结果可复现也省得每次调图都等计算。免疫浸润分析不难难的是每一步都留下可追溯的中间文件希望帮到你。本文还有配套的精品资源点击获取
返回列表