简介转录组测序揭示的是组织整体表达信号如何从中解析免疫微环境组成是肿瘤研究、免疫治疗响应预测及预后分析中的高频需求。免疫浸润分析旨在从bulk转录组数据中推断不同免疫细胞的相对丰度或活性而ssGSEA单样本基因集富集分析凭借原理透明、无需特征矩阵、单样本可算的优势成为零基础学习者首选算法。它通过对每个样本内部基因排序结合免疫细胞特征基因集计算富集分数输出样本-细胞类型分数矩阵进而支持热图、组间差异比较、相关性及生存分析等下游可视化。本文从R环境配置、表达矩阵整理、基因集选择到核心代码实操系统梳理了ssGSEA的完整流程并总结了基因ID转换、重复基因处理、表达量类型选择等常见踩坑细节助你快速掌握免疫浸润分析并产出可直接用于论文的图表。 转录组测序做完了差异基因筛出来上百个KEGG富集到几条通路你以为下游分析到此为止结果导师或者审稿人轻飘飘问一句样本的免疫微环境有没有差异——这是多数人接触免疫浸润的真实场景。作为转录组下游分析中最常用也最容易被追问的分析模块免疫浸润的任务是从bulk转录组数据里把组织整体表达信号拆解出不同免疫细胞的相对贡献。而在众多算法里ssGSEAsingle-sample Gene Set Enrichment Analysis单样本基因集富集分析是最适合零基础学习者跨过第一道门槛的算法原理直观、代码量小、结果好解释。这篇教程把我从配环境到出图踩过的坑全部整理出来配套代码和基因集资源一起提供跟着走完就能出免疫浸润热图和分组比较箱线图。1. 为什么转录组下游分析里免疫浸润永远绕不开1.1 你真正要回答的生物学问题先回到最根本的问题免疫浸润分析到底在算什么在肿瘤研究、自身免疫病、感染性疾病里组织的免疫微环境很大程度上决定了疾病走向。肿瘤组织里CD8 T细胞多不多、Treg多不多、巨噬细胞偏向M1还是M2这些信息直接和预后、免疫治疗响应相关。但bulk转录组测序有个天然限制——测出来的是整个组织里所有细胞的平均表达信号。你没法像单细胞测序那样直接数出样本里有几个CD8 T细胞。这时候就需要算法来做信号分离借助免疫细胞的特异性基因表达特征从混合转录组数据里推断这类免疫细胞的相对丰度或活性。要理解这个逻辑可以打个比方。你想知道一个合唱团里女高音声部的声音大不大但手里只有整场录音的混音文件。ssGSEA做的事情就是先列出一份女高音成员名单免疫细胞特征基因集然后去混音文件里看名单上的这些人是不是集中在音量最大的那一层。如果是就说明女高音声部在这场合唱里确实突出。这个思路用在生物学上就是把CD8 T细胞特征基因在样本里是否整体高表达作为判断CD8 T细胞浸润水平的依据。所以免疫浸润分析不是只看一两个marker基因而是看一组细胞类型特异基因的综合信号这比单基因判断稳健得多。1.2 主流免疫浸润算法对比为什么ssGSEA适合入门免疫浸润的算法从原理上可以分成两大类一类是基于反卷积的定量算法代表性的是CIBERSORT一类是基于基因集评分的富集类算法代表性的是ssGSEA、xCell、MCPcounter。两类算法解决的问题不完全一样输出结果的解释方式也不同。算法核心原理输出含义是否需要特征矩阵典型特点CIBERSORT线性支持向量回归反卷积22种免疫细胞比例需要LM22特征矩阵精度高但特征矩阵构建复杂比例和需为100%CIBERSORTxCIBERSORT升级版细胞比例/绝对丰度可自定义特征矩阵支持单细胞数据构建矩阵功能更强xCell基于ssGSEA的思想多种数据集校正64种细胞富集分数不需要细胞类型覆盖广但算法是一个黑盒MCPcounter基于标记基因表达计数10种免疫/基质细胞丰度不需要简单快速但覆盖细胞类型少ssGSEA单样本基因集富集分析任意基因集的富集分数不需要原理透明基因集可自定义入门门槛最低CIBERSORT在文献里出现频率最高但它需要LM22特征矩阵而且对输入数据的处理和p值解读有不少细节新手第一次跑很容易被绕晕。xCell虽然方便但内部做了很多数据校正改起来不灵活。相比之下ssGSEA有四个明确的优势原理透明就是排序加富集打分每一步都可以解释清楚。不需要特征矩阵只要有一份免疫细胞基因集就能跑。单个样本也能算每个样本独立打分不依赖其他样本的分布。基因集完全可以自定义不局限于免疫细胞换成任何感兴趣的基因集都能跑。对于零基础入门转录组下游分析我始终建议第一个免疫浸润算法就学ssGSEA因为它能让你把算法在做什么这件事想明白。想明白之后再去学CIBERSORT会轻松很多。2. ssGSEA的算法逻辑一次理清富集分数怎么来的2.1 从GSEA到ssGSEA从组间比较到单样本打分很多初学者看到ssGSEA这个缩写第一反应是是不是和GSEA长得很像。确实如此。GSEA基因集富集分析是经典的通路富集算法它的经典场景是对比两组样本先根据基因在两组间的差异程度给所有基因排个序然后看某个通路的基因是不是集中在排序列表的头部或尾部。如果是说明这个通路在两组之间有显著差异。但GSEA有个前提——它需要两组样本或者说需要基于样本分组来计算。这就带来一个实际问题如果我只有一批样本想对每个样本单独评估某个通路的活跃程度GSEA就不太合适了。ssGSEA就是在这个需求下提出的。Barbie等人在2009年发表ssGSEA算法时核心改动只有一个把组间比较换成每个样本内部比较。它先在单个样本内部把所有基因按表达量排序然后直接计算给定基因集在这个样本的排序列表里是否富集。这样每个样本都能得到一个独立的富集分数。放到免疫浸润场景里每个免疫细胞类型对应一个基因集每个样本跑一遍所有免疫细胞基因集的富集分析最后就得到一张样本 × 免疫细胞的分数矩阵。这个矩阵里的数值就是后续所有可视化和统计分析的原材料。2.2 ssGSEA分数的计算过程拆解ssGSEA的具体计算过程可以拆成四步我尽量讲得比教材更好懂。第一步对单个样本的所有基因按表达值从高到低排序得到一个排序列表L。表达量最高的基因排在第一位。第二步选定一个免疫细胞基因集S检查S里的基因在L中的位置分布。这里的关键是如果S里的基因大量集中在L的高位即这些基因在样本里表达量整体偏高说明该免疫细胞类型在这个样本里比较活跃。第三步计算富集分数。严谨一点说ssGSEA会计算基因集内基因和基因集外基因两个经验累积分布函数的差值然后沿着排序列表做一次类似Kolmogorov-Smirnov检验的随机游走把累积差值最大的偏移量作为该基因集的富集分数。这就是为什么你会在一些教程里看到ssGSEA分数有时候也叫enrichment score。加权版本还会把基因的表达值作为权重表达越高的基因贡献越大。第四步归一化。因为不同基因集的大小不同直接比较会不公平。ssGSEA会把原始富集分数除以基因集大小和排序列表长度造成的最大可能值最终得到在-1到1之间或类似范围的归一化分数。这四步走完一个样本的一个免疫细胞类型分数就出来了。不用被经验累积分布函数这种名词吓到你只需要记住最终效果一个基因集的基因如果在样本里整体高表达分数就高整体低表达分数就低。2.3 为什么基于排序的设计让ssGSEA这么稳ssGSEA有一个容易被低估的优点——它用的是表达值的排序而不是原始值。这个设计带来两个实际好处一是对数据尺度不敏感。同样是测序数据A实验室用TPMB实验室用FPKM绝对值会有差异。但基因之间的相对高低顺序通常不会因为定量方式不同而剧烈变化。所以ssGSEA在不同平台、不同定量方式的数据上都有一定的可比性。二是对极端值不敏感。RNA-seq数据里偶尔会出现个别基因表达量爆炸的情况比如某些核糖体蛋白或线粒体基因。如果用原始值计算这么一两个离群点就可能主导整个结果。但排序之后这些极端值最多排在第一位不会对后面的积累分数产生不成比例的影响。当然稳不等于万能。ssGSEA得到的分数是相对富集程度不是真实的细胞比例。它不能告诉你样本里CD8 T细胞精确占所有细胞的百分之几只能告诉你这批样本里CD8 T细胞的基因程序相对活跃。这一点在后面写文章、解释结果的时候一定要记住。3. 跑通ssGSEA前的环境与数据准备3.1 R环境与R包安装ssGSEA的实操基本都在R里完成。如果你用的是R 4.2及以上版本建议直接装新版GSVA包它是目前执行ssGSEA最方便的工具。需要安装的包分成两类核心分析包和可视化辅助包。# 安装BiocManager如果还没有 if (!require(BiocManager, quietly TRUE)) install.packages(BiocManager) # 核心分析包 BiocManager::install(c(GSVA, GSEABase)) # 可视化与数据处理辅助包 install.packages(c(pheatmap, ggplot2, ggpubr, reshape2, corrplot)) # 生存分析可选如果要做预后分析 BiocManager::install(c(survival, survminer))其中GSVA包从2023年底发布的1.50版本开始推荐使用新的参数对象写法ssgseaParam旧写法虽然还能用但会给出过时警告。这篇教程会以新版写法为主同时把旧版写法放在代码注释里供参考。3.2 表达矩阵的整理规范表达矩阵是ssGSEA的输入核心。它必须满足几个基本条件第一行名是基因列名是样本。不要反了这是新手最容易犯的错误。GSVA对基因在行、样本在列的格式要求很严格。第二矩阵里只能是数值不能有缺失值。如果有NA需要先处理否则会在运算时报错或者算出NA分数。第三行名建议使用基因Symbol。后面要对齐免疫细胞基因集基因集几乎都是Symbol格式。如果你的表达矩阵行名是Ensembl ID需要先做ID转换具体方法在第6部分踩坑实录里详细讲。第四表达量建议使用log2(TPM1)或log2(CPM1)。如果你用的数据是芯片数据通常已经做了log2处理可以直接使用。如果是从TCGA下载的Count数据建议先转换成TPM再做log2变换。count数据直接跑ssGSEA的问题在于不同样本测序深度差异很大排序会被高表达的看家基因主导。第五处理重复基因。如果有多行是同一个基因Symbol比如不同转录本需要先合并成一行否则GSVA会报错或者产生不可靠的结果。这个在第6.2节会给出具体处理策略。3.3 免疫细胞基因集的来源与选择ssGSEA所需的免疫细胞基因集最常用的是文献发表的基因集格式一般是GMT。以下的几个选择在文献中出现频率最高Bindea et al. 2013Immunity24种免疫细胞基因集经典老牌。Charoentong et al. 2017Cell Reports28种免疫细胞基因集覆盖更细包括各类T细胞亚群、NK细胞、树突状细胞等是目前最常用的选择之一。Thorsson et al. 2018ImmunityTCGA免疫亚型相关的特征基因集适合做泛癌分析。LM22CIBERSORT配套22种免疫细胞特征矩阵虽然是CIBERSORT的配套但也可以提取基因集部分用于ssGSEA。我个人的建议是入门阶段先用Charoentong 2017的28种免疫细胞基因集覆盖全面注释清楚教程配套资源里也包含了这份GMT文件。拿到基因集之后先别急着跑。要做一步命中率检查用intersect()看一下基因集里的基因在表达矩阵里实际命中多少。如果某个免疫细胞基因集在表达矩阵里只命中三五个基因这个分数基本不可靠需要在结果解读时剔除或者更换基因集。4. 核心代码实操从表达矩阵到免疫浸润矩阵4.1 载入表达矩阵与基因集下面的代码假定你已经把表达矩阵保存为tab分隔的txt文件免疫细胞基因集保存为GMT格式文件。# ---------- 载入表达矩阵 ---------- expr - read.table(expression_matrix.txt, header TRUE, row.names 1, sep \t, check.names FALSE) expr - as.matrix(expr) # 过滤掉含NA的行 expr - expr[apply(expr, 1, function(x) sum(is.na(x))) 0, ] # 过滤掉全为0或全低表达的行 expr - expr[apply(expr, 1, function(x) max(x) 0), ] # 查看矩阵维度确认基因行、样本列 dim(expr) expr[1:5, 1:5]# ---------- 读取免疫细胞基因集 ---------- library(GSEABase) geneSets - getGmt(immune_cell_signatures.gmt) # 转成list格式ssGSEA要求基因集是list of character vectors gset_list - lapply(geneSets, geneIds) names(gset_list) - sapply(geneSets, setName) # 查看基因集名称和每个基因集的基因数量 sapply(gset_list, length)载入完成后检查基因集是否与表达矩阵匹配# 统计每个基因集在表达矩阵中命中的基因数 hits - lapply(gset_list, function(gs) intersect(gs, rownames(expr))) sapply(hits, length)如果某个基因集命中数明显偏少比如少于10个建议直接删除这个基因集或者换更合适的版本。基因集命中太少时算出的富集分数很容易被一两个基因主导不稳定。4.2 运行ssGSEA新旧版GSVA兼顾GSVA包新版写法基于参数对象运行更规范报错信息也更友好。library(GSVA) # ---------- 新版GSVA ( 1.50) 推荐写法 ---------- ssgsea_param - ssgseaParam(exprData expr, geneSets gset_list, norm TRUE) ssgsea_score - gsva(ssgsea_param)如果你在别人代码里看到的是下面这种老写法也可以运行但建议尽早迁移到新写法# ---------- 旧版GSVA写法仅供参考不推荐在新环境使用 ---------- # ssgsea_score - gsva(as.matrix(expr), gset_list, # method ssgsea, # kcdf Gaussian, # verbose TRUE)有人会问norm参数是什么含义。norm TRUE表示对基因集大小做归一化让不同细胞类型之间的分数具有可比性。如果不归一化基因集大的细胞类型天然会获得更高分数这会给后续比较带来偏差。默认保持TRUE即可。4.3 输出结果的结构与含义ssGSEA跑完之后ssgsea_score是一个矩阵行是免疫细胞类型列是样本值是该免疫细胞类型在该样本中的富集分数比如ssgsea_score[CD8_T_cells, Sample01]就是Sample01里CD8 T细胞相关基因程序的富集分数。看到这里可以顺手把结果保存下来# 保存免疫浸润分数 write.csv(ssgsea_score, file ssgsea_immune_scores.csv)这个矩阵就是后续所有可视化和统计分析的入口。你可能会注意到分数有正有负这是归一化后的正常现象。正负号不代表细胞存在或不存在只代表富集程度高于或低于平均水平。5. 结果可视化与下游组学分析的衔接5.1 免疫浸润热图与样本分层观察拿到了免疫浸润分数矩阵第一个最直观的呈现方式就是热图。library(pheatmap) # 构造样本分组注释示例前10个样本Tumor后10个样本Normal annotation_col - data.frame( Group factor(c(rep(Tumor, 10), rep(Normal, 10))) ) rownames(annotation_col) - colnames(ssgsea_score) # 画热图对行做标准化 pheatmap(ssgsea_score, scale row, annotation_col annotation_col, show_colnames FALSE, color colorRampPalette(c(#4B0082, white, #FF4500))(100), main Immune Cell Infiltration (ssGSEA))这里解释一下scale row的作用。ssGSEA分数在不同免疫细胞类型之间的绝对水平不同有的细胞类型整体分数高有的整体分数低。如果直接用原始分数画热图颜色深浅会被整体水平主导看不出样本之间的差异。按行标准化之后每一行内部比较能更清晰地看到哪些样本在某种免疫细胞上偏高、哪些偏低。我在实际项目里观察热图的第一眼永远不是看某个具体值而是看样本是否按分组聚成两个明显簇。如果Tumor样本和Normal样本在热图上犬牙交错说明免疫浸润模式在两组之间差异不大如果形成明显的两个簇那后面做差异分析的把握就大很多。5.2 组间差异比较与统计检验热图适合全局观察但统计分析还是需要组间比较。最常见的做法是把某一种或多种免疫细胞的分数按分组画箱线图并做显著性检验。library(reshape2) library(ggplot2) library(ggpubr) # 把矩阵转成长格式 score_m - melt(ssgsea_score, varnames c(CellType, Sample)) score_m$Group - ifelse(grepl(^T, score_m$Sample), Tumor, Normal) # 以CD8 T细胞为例 df_cd8 - subset(score_m, CellType CD8_T_cells) ggboxplot(df_cd8, x Group, y value, color Group, palette jco, add jitter) stat_compare_means(method wilcox.test, label p.format)这里用Wilcoxon秩和检验是比较稳妥的因为ssGSEA分数通常不完全符合正态分布。样本量小的时候也可以用它不要求正态性。如果是三组及以上比较用kruskal.test这也是非参数方法不需要正态假设。画图的时候注意一个小技巧如果样本量很小散点和箱线图叠加容易看不清分布。这时候可以只用散点加中位数线或者增加alpha 0.5让点半透明。5.3 免疫细胞相关性分析免疫细胞之间经常存在协同或拮抗关系比如CD8 T细胞和Treg在肿瘤微环境里往往呈负相关因为Treg会抑制CD8 T细胞的活性。相关性分析可以帮助你发现这些关系。library(corrplot) # 计算免疫细胞之间的Spearman相关系数 cor_mat - cor(t(ssgsea_score), method spearman) # 相关性热图 corrplot(cor_mat, method color, type upper, tl.cex 0.8, tl.col black, addCoef.col grey30, number.cex 0.6)这里用Spearman而不是Pearson是因为ssGSEA分数之间的关系不一定是线性的基于排序的Spearman更稳健。当你看到某对细胞类型相关系数大于0.6或者小于-0.6通常在生物学上值得拿来讨论。5.4 从浸润分数到生存分析在肿瘤研究里免疫浸润最常做的下游分析之一就是和生存数据结合评估某类免疫细胞的浸润水平是否影响预后。library(survival) library(survminer) # 假设你已经有了临床数据clin包含样本ID、OS.time、OS.event # OS.time单位是月OS.event 0删失 1死亡 # 转置ssGSEA分数样本变成行 ssgsea_df - as.data.frame(t(ssgsea_score)) # 以CD8 T细胞中位数为界分为高低两组 ssgsea_df$CD8_Group - ifelse( ssgsea_df$CD8_T_cells median(ssgsea_df$CD8_T_cells), High, Low ) ssgsea_df$SampleID - rownames(ssgsea_df) # 合并临床信息 surv_data - merge(ssgsea_df, clin[, c(SampleID, OS.time, OS.event)], by SampleID) # 生存曲线与log-rank检验 fit - survfit(Surv(OS.time, OS.event) ~ CD8_Group, data surv_data) ggsurvplot(fit, data surv_data, pval TRUE, risk.table TRUE, palette c(#E64B35, #4DBBD5), legend.title CD8 T cell)这里有一个我反复强调的经验以中位数分组是一种常见做法但它只有在样本量足够、且样本整体分布没有明显偏移时才靠谱。如果你的样本数少于30或者高低分组后某组只剩个位数样本生存分析结果的可靠性就要打问号。这种情况下可以考虑用连续变量做单因素Cox回归输出风险比而不是简单的高低分组。6. 新手最容易翻车的四个细节6.1 基因名格式与ID转换这是我见过最多人卡住的地方。ssGSEA本质上是在做表达矩阵的基因名和基因集的基因名之间的匹配。只要两边格式不一致命中率就会掉得厉害结果直接失效。常见的格式问题有表达矩阵行名是Ensembl ID基因集是Symbol表达矩阵行名带版本号比如ENSG00000141510.17基因名大小写不一致小鼠基因名是Cd8a人的基因集一般是大写CD8A基因名里带了_、-等特殊字符。比较稳妥的做法是提前统一成Symbol并且检查一遍。library(org.Hs.eg.db) # 如果你的行名是Ensembl ID ids - mapIds(org.Hs.eg.db, keys rownames(expr), keytype ENSEMBL, column SYMBOL) # 去掉无法映射的基因 expr_symbol - expr[!is.na(ids), ] rownames(expr_symbol) - ids[!is.na(ids)] # 再去掉仍然为空的行 expr_symbol - expr_symbol[rownames(expr_symbol) ! , ]映射完成之后再跑一次基因集命中率检查。如果命中率没有显著提升就要看看基因集文件本身有没有问题。6.2 重复基因的处理策略表达矩阵里出现重复基因名是另一个高频问题。芯片数据里同一个基因可能有多个探针测序数据里同一个基因可能对应多个转录本ID转换完Symbol之后就出现了重复行。GSVA遇到重复行名时不同的GSVA版本报错方式不一样但总之不能直接跑。处理策略一般有两种取最大值或取平均值。# 检查重复基因 table(duplicated(rownames(expr))) # 用aggregate对重复基因取最大值 expr_dedup - aggregate(expr, by list(gene rownames(expr)), FUN max) rownames(expr_dedup) - expr_dedup$gene expr_dedup$gene - NULL expr_dedup - as.matrix(expr_dedup)我个人的偏好是取最大值而不是平均值。原因很简单芯片或转录本级别的表达量取平均会把真正的高表达信号稀释掉而取最大值能保留这个基因在该样本中最强的信号。当然如果你的目的是降噪而不是保留信号取平均值也可以但一致性原则是——一个数据集里只能统一用一种方式不能一部分基因取max另一部分取mean。6.3 表达量类型的选择TPM、CPM还是countssGSEA基于排序理论上比基于绝对值的算法对表达量类型更不敏感。但排序不等于不需要预处理。如果用count数据跑会有两个问题第一不同样本的测序深度不同。测序深度高的样本所有基因的count值整体偏高排序时高表达基因占优这会让样本间的差异部分被测序深度因素主导。第二基因长度不同。count值本身没有纠正基因长度长基因通常count数更高。这在基因集比较时可能会引入偏差——如果某个免疫细胞基因集恰好包含较多长基因分数就会虚高。所以我的建议是RNA-seq数据至少转成CPM再log2最好是转成TPM再log2。Count → CPM的转换很简单按每个样本的总count数缩放即可。TPM的转换公式稍微复杂一点建议参考R包edgeR的cpm()函数或DESeq2的vst()函数辅助处理。芯片数据则不需要纠结这个问题因为芯片数据在标准化阶段已经做了校正通常直接用matrix series文件里的表达值就行。6.4 基因集版本与参数设置的影响同样一份表达矩阵用不同的免疫细胞基因集会得到不同的结果。这不是bug而是基因集定义差异带来的必然结果。Charoentong 2017和Bindea 2013对CD8 T细胞的基因集定义不完全相同所以你在写文章时必须明确说明用的是哪个版本、多少种免疫细胞否则审稿人无法复现。GSVA的ssgseaParam里还有两个参数值得留意ssgsea_param - ssgseaParam(exprData expr, geneSets gset_list, norm TRUE, minSize 10, maxSize 500)minSize设定基因集包含基因数的下限低于这个值的基因集不会被计算maxSize设定上限。默认情况下minSize可能设得很低我建议手动设成10。基因集太小的时候富集分数的方差极大一个基因的改变就能把分数从正拉到负这种结果不稳健。还有一个很容易被忽略的细节如果表达矩阵里有大量基因在多个样本里表达量为0建议先做一步低表达过滤只保留在部分样本中有一定表达量的基因。全零基因行参与排序会让排序列表里出现大面积的并列最低值虽然不是致命错误但会影响富集分数的敏感性。7. 配套资源的使用说明与后续拓展7.1 配套资源里有什么这套教程的配套资源包含四个部分示例表达矩阵模拟的1000基因 × 20样本数据其中10个Tumor、10个Normal。模拟数据特意保留了一些真实数据中常见的特性比如部分基因表达为0、个别高表达离群基因方便你体验完整流程。免疫细胞基因集Charoentong 2017的28种免疫细胞GMT文件覆盖T细胞亚群、B细胞、NK细胞、巨噬细胞亚型等可以直接用于ssGSEA分析。完整R脚本从读入数据、运行ssGSEA到出热图和箱线图的完整脚本脚本里的路径和你下载到的资源目录做了对应。结果示例图我跑这批示例数据得到的免疫浸润热图、CD8 T细胞分组比较图供你核对结果是否一致。拿到资源后推荐的操作路径是这样新建一个RStudio Project把资源文件夹放进去先跑示例数据确认结果和示例图一致然后再把自己的表达矩阵按同样的格式放进去替换。7.2 拿到素材后的完整操作路径第一步用RStudio打开项目目录下的ssGSEA_tutorial.R不要一开始就全选运行。逐段运行每运行一小段就停下来看清楚数据结构。我特别建议在运行ssGSEA之前先执行str(gset_list)和dim(expr)确认两边格式没问题再继续往下。第二步跑通示例数据确认输出文件ssgsea_immune_scores.csv的行列数和预期一致28种免疫细胞 × 20个样本。第三步替换成自己的数据。把自己的表达矩阵命名为expression_matrix.txt保持同样的格式覆盖示例文件重新运行。如果自己的数据太大比如超过2万个基因跑起来会比较慢这是正常的。第四步根据实际分组修改可视化代码里的分组变量。示例代码里写死了前10个样本Tumor后10个样本Normal你替换数据后要改成自己的分组信息。这一步骤人最容易搞混建议在R里先用colnames(expr)确认样本顺序别凭感觉分。7.3 从ssGSEA到更多免疫浸润算法的拓展跑通ssGSEA只是入门后面有几条很自然的进阶路线一条是把ssGSEA和CIBERSORTx结合起来。ssGSEA看相对富集程度CIBERSORTx输出细胞比例两者结合使用一个看趋势一个看幅度结论互相印证。具体的做法是先跑CIBERSORTx可以在线服务拿到22种免疫细胞比例矩阵再和ssGSEA分数矩阵做相关性分析看哪些细胞类型在两个算法下结论一致。另一条是把免疫浸润和WGCNA加权基因共表达网络分析结合。先跑WGCNA得到基因模块计算模块特征基因与免疫浸润分数的相关性筛选出与特定免疫细胞浸润显著相关的模块再去模块里找hub基因。这个方向在肿瘤免疫相关论文里很常见也是把免疫浸润从描述性分析变成机制探索的常见路径。还有一条是结合公共数据做扩展。TCGA的转录组数据可以直接下载很多肿瘤类型都有几百个样本跑一遍ssGSEA之后你可以做泛癌的免疫浸润模式比较、做免疫亚型分类、做免疫评分与药物敏感性之间的关联。TCGA免疫浸润生存分析这个组合做出来的内容足够撑起一篇生信论文的核心分析部分。最后说点我自己的体会。我带过的学生里跑通代码的一天就能出图真正花时间的往往是结果解释。ssGSEA给的是一个相对富集分数不是真实细胞比例写文章时不要写成CD8 T细胞含量显著升高更稳妥的表述是CD8 T细胞相关基因程序富集分数升高或CD8 T细胞浸润水平升高。审稿人看的不是你图多漂亮而是你有没有准确理解自己这个分数的生物学含义。这也是为什么我在这篇教程里花了大量篇幅讲算法原理和参数含义——代码是最容易复制的东西但理解算法为什么这么做才是你以后不被各种报错难住的关键。本文还有配套的精品资源点击获取