尧图网站建设 
简介面向零基础转录组学习者提供CIBERSORT免疫浸润分析的一体化配套资源涵盖输入数据、R脚本与输出结果适合想用免疫浸润算法解析表达矩阵、绘制细胞比例可视化图表的读者也适合结合博文教程逐步上手实践。包内共13个文件以csv数据表、txt结果文件、R运行脚本为主并附PDF与PNG可视化图片整体约144.83MB。csv与txt分别存放表达矩阵、样本分组信息及CIBERSORT运行结果R脚本已实测可一键全选后跑通便于直接复现代分析流程。目前已有240人学习资源附有总目录跳转与配套教程链接可按需选择知识点学习。通过输入数据、完整脚本和输出结果三部分读者能快速理解免疫浸润分析的核心步骤节省环境配置与排错时间尤其适合零基础入门者对照学习。1. CIBERSORT不是黑匣子免疫浸润结果是怎么算出来的拿到一批肿瘤样本的转录组表达矩阵后最常见的下游问题不是「哪些基因差异表达」而是「肿瘤微环境里到底有哪些免疫细胞、各占多少」。CIBERSORT 是目前被引最多的免疫浸润推断工具之一它不靠显微镜也不靠流式分选而是从批量转录组数据里用数学方法把 22 种免疫细胞的相对比例拆出来。这份配套资源包含完整可跑的 R 脚本、LM22 特征矩阵和示例数据属于「零基础能直接复现」的类型。它面向的是已经拿到表达矩阵、想做免疫浸润却无处下手的生信新手以及想快速验证一批样本免疫组成的科研人员。全文不涉及算法黑匣子的空谈直接告诉你每一步怎么选参数、怎么避坑。2. 免疫浸润分析的选型逻辑为什么 CIBERSORT 适合零基础上手2.1 免疫浸润分析的两条路线去卷积与富集打分免疫浸润分析在方法学上大体分两条路线。第一条是富集打分路线代表工具是 ssGSEA 和 xCell。这类方法把免疫细胞的标记基因集合看成一个个「特征库」然后拿样本的表达谱去和特征库做富集分析输出一个相对活性分数。它不要求所有标记基因的表达值都在一个量纲下对输入数据的格式宽容度较高。第二条是去卷积路线代表工具就是 CIBERSORT。它的逻辑更接近解方程把样本的表达谱看成 22 种免疫细胞表达谱的线性混合目标是求出一组比例系数让这组系数乘以各细胞类型的特征表达矩阵后能尽可能还原观察到的样本表达谱。对刚入门的人来说这两条路线最直观的差别在输出结果上。ssGSEA 给的是一个分数这个分数只在同一个数据集内部做相对比较才有意义而 CIBERSORT 输出的是百分比所有细胞类型的比例加起来接近 1。你从 CIBERSORT 结果里能直接说「这个样本里 CD8 T 细胞约占 12%」这是很多临床研究者习惯的表述方式。也正因为这个输出形式直观CIBERSORT 在肿瘤免疫相关的转录组分析里出场率极高。不过它的代价是输入矩阵的格式要求更严格。基因名必须是标准的 Gene Symbol表达值最好转换成 TPM而且对数据质量非常敏感——这些坑后面会一条条展开。我自己的习惯是如果只想快速看看趋势用 ssGSEA如果要做分组比较且最终图表要呈现比例关系就上 CIBERSORT。2.2 CIBERSORT 的核心机制ν-SVR 与 LM22 特征矩阵要跑好 CIBERSORT至少要理解它的两个核心组件。第一个是 LM22 特征矩阵。这是一个包含 22 种免疫细胞亚型、由 547 个标记基因构成的「标准表达谱」。这 547 个基因是作者从纯化细胞样本的芯片数据里筛出来的每种细胞类型在该矩阵里对应一组相对特异的表达模式。第二个是 ν-SVR全称是支持向量回归。CIBERSORT 没有用传统的最小二乘拟合而是用 ν-SVR 来做回归求解。为什么这么做因为免疫细胞的表达谱之间有很强的相关性普通线性回归在多重共线性下会给出非常不稳定的解ν-SVR 通过引入一个不敏感损失带允许一定的拟合误差换来的是比例解在小扰动下的稳定性。这两个组件决定了 CIBERSORT 的输入约束。LM22 的特征基因只有 547 个所以你的表达矩阵里必须能匹配到足够多的这些基因。如果匹配上的基因数太少算法会直接提示结果不可靠。另外LM22 建于 Affymetrix 芯片平台它记录的是绝对荧光信号强度这与 RNA-seq 的计数单位有天然差异。所以原版论文和后续实践都建议把 RNA-seq 数据先转换成 TPM 再输入而不是直接用 FPKM 或 read count。这一点很关键我在 3.1 节会给出具体的转换代码。很多新手在这里翻车——表达矩阵的基因名倒是匹配上了但因为量纲差异跑出来的比例分布明显偏离预期。提示CIBERSORT 的名字也意味着它的输出本质是「相对丰度」。22 种细胞的比例总和为 1某一种细胞的比例升高必然伴随另一种相对下降。这个数学约束会在免疫学解读时带来一些反直觉的结论后面第 5 章会专门讨论。2.3 输入数据的工程要求基因名与表达量单位把技术细节讲透前先给你一个我在实际数据上反复踩出来的结论CIBERSORT 对输入矩阵的要求可以归纳成「六个字」——去重、改名、转单位。去重指的是表达矩阵中不能有重复的基因名。很多从 GEO 下载的矩阵或自己比对产生的 count 矩阵里一个基因 Symbol 对应多行转录本或探针如果不处理就直接跑CIBERSORT 会报错让你无从下手。改名指的是确认你的矩阵行名是 Gene Symbol 而不是 Ensembl ID 或 Entrez ID。LM22 用的是 Gene Symbol如果你的矩阵是 Ensembl ID哪怕基因完全重合也匹配不上。转单位就是前面强调的 RNA-seq 数据要转 TPM。这三个要求里最容易忽略的是第一个。我见过不少同学拿到矩阵后直接丢进 CIBERSORT报错了才发现基因有重复。处理重复基因的标准做法是按基因取均值聚合而不是简单去重保留第一行。保留第一行会丢表达信息而且在后续可视化里会出现奇怪的离群样本。聚合操作我一般用两行aggregate就能完成具体代码在 3.1 节。你可以把这一小节当成一个检查清单——先把这三件事做完再谈跑算法否则后面的每一步都在为前面的马虎买单。3. 跑通第一次 CIBERSORT表达矩阵、LM22 与三个关键参数3.1 数据准备读入表达矩阵并做基因去重聚合配套资源里自带一个示例表达矩阵expression_matrix.txt为了让你理解格式我先说明它的结构第一列是 Gene Symbol其余列是样本名单元格里的数值是表达量。日常工作中你的矩阵可能是从 featureCounts、Salmon 或者 GEO 下载后整理而来的但进入 CIBERSORT 前的处理逻辑是一样的。下面这段代码完成读入与去重聚合# 读入表达矩阵第一列基因名作为行名 expr - read.table(expression_matrix.txt, header TRUE, row.names 1, sep \t, check.names FALSE) # 把基因名从行名提取为显式列方便 aggregate 按基因分组 expr$gene - rownames(expr) # 对同一基因的多行表达值取均值得到唯一基因行 expr_agg - aggregate(. ~ gene, data expr, FUN mean) # 恢复基因名为行名删除多余的 gene 列 rownames(expr_agg) - expr_agg$gene expr_agg$gene - NULL # 保留基因名与 LM22 特征矩阵的交集并输出基因匹配数 lm22 - read.table(LM22.txt, header TRUE, row.names 1, sep \t, check.names FALSE) common_genes - intersect(rownames(expr_agg), rownames(lm22)) cat(匹配到 LM22 的基因数, length(common_genes), \n) # 用交集基因重新过滤表达矩阵和 LM22保证两者基因顺序一致 expr_final - as.matrix(expr_agg[common_genes, ]) lm22_final - as.matrix(lm22[common_genes, ])代码逻辑上分成三块。第一块是读入并检查格式check.names FALSE是为了防止 R 把样本名里的横线或空格改成点号这个细节在后续匹配样本时很关键。第二块用aggregate按基因分组取均值它处理的是多对一映射比直接duplicated删除更稳妥因为取均值保留了该基因在多个转录本上的整体表达水平。第三块是求交集并过滤。这里我建议你务必打印length(common_genes)看一眼如果这个数字低于 300说明你的基因名版本或 ID 类型有问题建议回头检查转换步骤。我遇到过一次匹配数只有 200 出头的情况原因是表达矩阵用的是版本较旧的 Gene Symbol和 LM22 的命名有出入换用最新的注释文件后就恢复正常了。3.2 运行 CIBERSORT脚本调用与 perm、QN 的取舍数据准备完毕后进入正式运行环节。配套资源里的核心文件是cibersort.R它是一个封装好的函数库。使用时先source加载函数再调用CIBERSORT函数。下面是最小可运行的调用方式# 加载 CIBERSORT 官方函数 source(cibersort.R) # 调用主函数输入为 LM22 文件路径、表达矩阵、排列次数、是否做分位数归一化 res - CIBERSORT(sig_matrix LM22.txt, mixture expr_final, perm 1000, QN TRUE) # 查看结果的前几行和前几列确认输出格式 res[1:5, 1:5]这里的参数有三个需要重点解释。第一个是perm它代表置换次数CIBERSORT 用它来估计每个样本结果的显著性 P 值。默认值是 100但实践中 100 次置换得到的 P 值波动很大同一个样本跑两次可能一个显著一个不显著。我在正式分析里一般设 1000这会增加一些运行时间但对结果的稳定性提升非常明显。第二个是QN布尔值代表是否做分位数归一化。如果输入数据来自芯片平台设TRUE是合理选择如果输入的是 RNA-seq 的 TPM 数据是否设TRUE存在争议。原因是分位数归一化会把不同样本的表达分布强行拉齐这在芯片数据里是标准操作但对于 TPM 数据部分方法学文章认为它会引入虚假信号。我的经验是如果你的样本分组来自同一批次测序设TRUE问题不大如果样本来自不同平台或不同批次建议设FALSE并先自己用limma的normalizeBetweenArrays做过一次归一化再输入。第三个参数mixture它接收的是矩阵对象而不是文件路径所以必须先读入并处理成上一步的expr_final。3.3 结果文件解析三列统计量不是摆设运行结束后res是一个行为样本、列为细胞类型的矩阵但它的列并不全是比例值。最后三列依次是P-value、Correlation、RMSE。这三个统计量很重要很多人只看比例忽略了它们。P-value来自置换检验它回答的问题是「这个样本的免疫细胞比例估计是否显著优于随机」。如果某样本的 P 值大于 0.05说明该样本的 deconvolution 结果并不可靠后续做差异分析时建议谨慎纳入或直接剔除。Correlation表示用估计出的细胞比例重建表达谱后与原始表达谱的相关性。这个值越高说明拟合越好。RMSE是均方根误差越低越好。一个常见的筛查习惯是剔除Correlation小于 0.8 或P-value大于 0.05 的样本再进入下游分析和绘图。# 提取细胞比例部分去掉最后三列统计量 props - res[, 1:22] # 检查有多少样本 P 值大于 0.05决定是否需要剔除 pval_col - res[, P-value] cat(P值大于0.05的样本数, sum(pval_col 0.05), \n) # 如果只是看分组趋势可以先把 P 大于 0.05 的样本标记出来不急着删 res - as.data.frame(res) res$sample_id - rownames(res)这里我补充一个实际经验如果你跑出来的所有样本 P 值都小于 0.05但比例分布很不合理比如所有样本的巨噬细胞都是 0那大概率不是算法问题而是表达矩阵里对应 LM22 基因的表达值缺失严重。检验方法很简单把表达矩阵里那部分巨噬细胞标记基因单独提取出来看表达量分布。如果大量基因表达值都是 0说明你的上游定量流程或基因注释有问题跟 CIBERSORT 参数无关。先查数据再怀疑算法这是避免浪费半天时间的原则。4. 把结果画到能投稿堆叠条形图、热图与小提琴图4.1 堆叠条形图一眼看出每个样本的免疫组成CIBERSORT 输出的比例矩阵本身是一张表格但审稿人更愿意看图。最常用的第一张图是堆叠条形图每个样本一根柱子柱子的不同颜色段代表不同免疫细胞的占比。这张图能直观展示样本间的免疫组成差异。画图首选ggplot2但需要先把比例矩阵从宽格式转换成长格式这一步对 ggplot2 来说是必须的library(ggplot2) library(reshape2) # props 是第 3 节提取的 22 列比例矩阵先加一列样本名 props_df - as.data.frame(props) props_df$sample_id - rownames(props_df) # 宽表转长表每一行变成一个样本-细胞类型的组合 props_long - melt(props_df, id.vars sample_id, variable.name cell_type, value.name proportion) # 按样本排序并指定细胞类型在堆叠条中的顺序 props_long$sample_id - factor(props_long$sample_id, levels rownames(props_df)) # 堆叠条形图positionfill 让每根柱子高度都为 1 ggplot(props_long, aes(x sample_id, y proportion, fill cell_type)) geom_bar(stat identity, position fill, width 0.8) theme_minimal() theme(axis.text.x element_text(angle 60, hjust 1)) labs(x 样本, y 细胞比例, fill 细胞类型)长格式的转换是画图前的固定动作melt把 22 列比例值堆叠成一列id.vars指定保留样本名列。position fill让每根柱子等高展示的是组成结构而不是绝对含量如果你想让柱子的高度反映某个细胞类型分数的高低改用position stack即可。横轴样本标签比较多时可以旋转 60 度避免重叠。如果你希望图里每个样本按照某个细胞比例升序排列把levels改成按该细胞比例排序后的样本名向量就行。再有就是配色22 种细胞类型手动配色很麻烦直接用scale_fill_manual配一套有色差的色板别用默认的 hue 循环打印出来分不清相近色。4.2 热图与小提琴图分组比较的标准配置堆叠条形图适合展示单样本组成但当你需要比较肿瘤组和对照组时热图和小提琴图是更通用的选择。热图适合展示所有样本、所有细胞类型表达矩阵的整体格局通常配合样本聚类来观察是否存在亚型分组。小提琴图则适合对某一细胞类型做两组差异比较。下面的代码分别给出热图和小提琴图的最小实现library(pheatmap) # 热图对细胞比例做行标准化聚类展示格局 pheatmap(as.matrix(props), cluster_cols TRUE, cluster_rows TRUE, scale row, border_color NA, color colorRampPalette(c(navy, white, firebrick3))(100)) # 小提琴图以 CD8_T_cells 为例做两组比较 library(ggplot2) # 构造示例分组信息实际分析中替换为你的临床分组表 group_info - data.frame( sample_id rownames(props), group c(rep(Tumor, 10), rep(Normal, 10)) ) # 取目标细胞比例并合并分组 plot_df - data.frame( sample_id rownames(props), cd8 props[, T cells CD8], group group_info$group[match(rownames(props), group_info$sample_id)] ) ggplot(plot_df, aes(x group, y cd8, fill group)) geom_violin(trim FALSE) geom_boxplot(width 0.15, outlier.shape NA) geom_jitter(width 0.1, size 1.5, alpha 0.6) theme_classic() labs(x 分组, y CD8 T 细胞比例)热图的scale row会按细胞类型做标准化这会让颜色反映的是每种细胞在不同样本间的相对高低而不是原始比例。如果你的目的是直接比较比例大小把scale去掉用原始比例值配色更符合直觉。小提琴图里我刻意叠加了箱线和散点这个组合是生信论文里最常见的呈现方式。group_info里我假设前 10 个样本是肿瘤、后 10 个是正常实际使用时请按照你的样本注释表格改务必注意样本名与rownames(props)严格对应。另外要留意T cells CD8这个列名的写法不同版本的 LM22 命名略有差异小写、空格、下划线的细微差别会让取列时报错。稳妥做法是先colnames(props)看一眼实际列名再写代码。提示小提琴图背后的统计检验没有写进代码但正式分析时需要用wilcox.test(cd8 ~ group, data plot_df)算一下 P 值再标到图上。样本量小的时候 P 值容易不显著这是正常现象不要为了出图强行套 t 检验。5. 避坑记录从报错到错误解读的五个高频问题5.1 报错与冲突R 版本、包依赖和函数缺失现象source(cibersort.R)之后调用函数报错could not find function ...或者 R 直接提示加载某个包失败常见于 R 4.2 及以上版本。原因CIBERSORT 的 R 脚本发布年份较早内部依赖了部分在后续 R 版本中被移除的包和函数例如gplots中的某些绘图函数和preprocessCore的编译依赖。新装的环境里这些包要么装不上要么函数名变了。解决我一般会先按顺序装齐三件套preprocessCore从 Bioconductor 安装e1071和gplots从 CRAN 安装。如果仍报错不要盲目改源码先看错误信息定位到具体函数再决定是降级安装老版本包还是绕道走。一个比较省事的方案是直接在项目目录下新建一个renv环境固定 R 版本为 4.1.x。这块没有秘诀就是依赖管理问题。5.2 结果异常全零、基因匹配数过低与 P 值全为 1现象跑出来的比例矩阵中某个或某几个细胞类型在所有样本中全为 0或者匹配基因数不到 200另外有一类情况是所有样本 P 值都是 1。原因全为 0 通常不是算法问题而是表达矩阵里这一类细胞的标记基因几乎不表达或者基因名版本与 LM22 对不上。匹配数过低基本可以断定是 ID 类型问题最常见的是 Ensembl ID 直接当 Symbol 用。P 值全为 1 则是因为perm参数太小或表达矩阵中有极端离群值导致 SVR 求解不稳定。解决先做两步排查。第一步打印交集基因数如果低于 300优先怀疑 ID 映射用biomaRt或本地注释文件做转换。第二步检查表达矩阵是否存在全零行或包含巨大离群值的样本这类行会在 SVR 求解时干扰权重计算。我习惯在进入 CIBERSORT 之前按表达量过滤一次删除在所有样本中表达量都低于某个阈值的基因。至于 P 值问题把perm提到 1000并且检查是否有某列表达值明显超出其他列几个数量级有的话先做 log2 变换再跑。5.3 数据来源混用把 FPKM 直接灌进 CIBERSORT现象用 TCGA 的 FPKM 矩阵直接跑结果和已发表文献里的免疫浸润分布对不上CD4 和 CD8 T 细胞比例总是偏低或偏高得很离谱。原因CIBERSORT 的 LM22 模型是在芯片绝对信号强度上训练的RNA-seq 数据输入前需要把 FPKM 转换成 TPM。FPKM 与 TPM 虽然在很多样本里高度相关但 FPKM 的样本间可比性不如 TPM尤其在基因长度偏好明显的转录组数据里直接用 FPKM 会让 SVR 的拟合残差变大。解决跑之前统一做一次 FPKM 到 TPM 的转换公式是 TPM FPKM / sum(FPKM) * 1e6。这个转换有现成函数自己写三行 R 代码也能完成。第 6 章我会给出可直接复制的实现。另外如果你的上游是 STAR 定量出的 read count建议先用DESeq2的vst或edgeR的cpm做标准化再转 TPM不要直接拿原始 count 进来。这里的核心认知是CIBERSORT 需要的是「相对丰度」语义的表达值任何破坏样本间可比性的操作都会污染结果。5.4 置换检验的玄学perm 次数与结果可复现性现象同一份数据、同一个参数跑两次 CIBERSORT 得到的比例矩阵并不完全一致某些细胞类型的比例在小数点后第二位浮动P 值甚至会在 0.04 和 0.06 之间来回跳。原因perm置换次数本质上是蒙特卡洛过程次数越少随机性越大。默认值 100 次对大多数样本够用但边缘显著的样本就会非常不稳定。解决方法只有一句话——把perm调到 1000 或更高。代价是运行时间线性增长一个 50 样本的矩阵在普通笔记本上大约跑几十秒到几分钟。从那以后我每次做正式分析都固定用perm 1000并且在论文方法部分明确写清楚。如果你跑的是大规模队列比如几百个样本可以先把perm设 100 跑一遍过滤掉明显可靠的样本再对边缘样本重跑高 perm这样能省不少时间。5.5 解读陷阱相对丰度不是绝对含量现象分析报告中写「A 组的 CD8 T 细胞比例显著高于 B 组」但结合病理切片或流式数据却发现实际绝对数量没有差异甚至方向相反。原因CIBERSORT 输出的是相对丰度所有细胞比例总和固定为 1。某一种免疫细胞比例升高不代表它的绝对数量增加也可能是其他细胞类型减少导致的「分母变小」。这是去卷积方法天然带有的数学约束不是算法 bug。解决在论文表述里避免使用「浸润程度显著升高」这类暗示绝对量的措辞改用「相对比例变化」或「在免疫细胞组成中占比发生变化」。如果必须与绝对数量挂钩建议联合病理切片的免疫组化染色或流式定量数据做验证。这个坑在投稿时最容易被审稿人抓住要格外注意。我在一次项目中被审稿人质疑过这点后来把结果描述改写为「相对比例改变」并且补充了一组 IHC 验证数据才过关。注意如果你看到某个样本的 22 种细胞比例之和不是 1比如 0.93 或 1.07这不算 bug。CIBERSORT 做了归一化处理但小数位舍入可能导致轻微偏差。偏差过大超过 0.1时通常意味着表达矩阵里存在异常值。6. 进阶用法FPKM 转 TPM、批量运行与结果落地当单次运行已经熟练之后下一步值得做的事有三件把表达量单位转换脚本化、让 CIBERSORT 在多个数据集或多种参数下批量运行、以及把结果与下游分析衔接起来。关于 FPKM 转 TPM下面是固定套路# fpkm_matrix 是行为基因、列为样本的 FPKM 矩阵 fpkm_to_tpm - function(fpkm_matrix) { # 对每个样本每个基因的 FPKM 除以该样本 FPKM 总和再乘 1e6 tpm - t(t(fpkm_matrix) / colSums(fpkm_matrix) * 1e6) return(tpm) } tpm_matrix - fpkm_to_tpm(fpkm_matrix)这个函数内部做了两次转置第一次t(fpkm_matrix)把样本放到行上colSums后分母变成每个样本的总 FPKM外层再转置回来保持基因在行。函数式写法是一劳永逸的。我一般在拿到任何来源的表达矩阵时都会先确认它的单位然后第一时间转成 TPM再存入后续分析流程。批量运行的核心是循环封装。你可以把第 3 节的整个流程包成一个函数输入是表达矩阵文件路径和分组信息输出是比例矩阵和对应图。注意每次循环里都需要重新做基因交集过滤和source加载或者把cibersort.R的依赖包提前加载好放在循环外。批量运行时另一个实用技巧是把结果统一写入一个汇总目录文件名带上数据集名称和perm值防止自己把不同参数的结果搞混。至于与下游分析的衔接我推荐把 CIBERSORT 比例矩阵保存为 CSV 后直接读入 GSVA 的结果一起做相关分析——免疫浸润比例与通路活性矩阵的相关性热图是肿瘤转录组论文里出现频率很高的组合图。做法很简单cor_mat - cor(props, gsva_scores, method spearman)然后用pheatmap画相关性热图p.adjust做多重校正后筛选显著配对。最后想说的是整个过程中最值得养成的一个习惯是拿到任何表达矩阵后不要急着跑 CIBERSORT先在纸上列出基因名版本、表达量单位、是否有重复基因、LM22 匹配数这四项检查清单。从那次被审稿人质疑相对丰度解读之后我每做一个免疫浸润项目都会强制走一遍这个检查流程再顺手把 FPKM 转 TPM 的脚本跑一遍确认结果分布没有离群。这份配套资源里的示例数据和脚本足够你把整个流程从上游矩阵一路跑到出图但数据和参数替换成你自己的项目时上面的坑一个都省不掉。希望帮到你。本文还有配套的精品资源点击获取 CIBERSORT免疫浸润分析：从转录组表达矩阵到细胞比例推断
简介面向零基础转录组学习者提供CIBERSORT免疫浸润分析的一体化配套资源涵盖输入数据、R脚本与输出结果适合想用免疫浸润算法解析表达矩阵、绘制细胞比例可视化图表的读者也适合结合博文教程逐步上手实践。包内共13个文件以csv数据表、txt结果文件、R运行脚本为主并附PDF与PNG可视化图片整体约144.83MB。csv与txt分别存放表达矩阵、样本分组信息及CIBERSORT运行结果R脚本已实测可一键全选后跑通便于直接复现代分析流程。目前已有240人学习资源附有总目录跳转与配套教程链接可按需选择知识点学习。通过输入数据、完整脚本和输出结果三部分读者能快速理解免疫浸润分析的核心步骤节省环境配置与排错时间尤其适合零基础入门者对照学习。1. CIBERSORT不是黑匣子免疫浸润结果是怎么算出来的拿到一批肿瘤样本的转录组表达矩阵后最常见的下游问题不是「哪些基因差异表达」而是「肿瘤微环境里到底有哪些免疫细胞、各占多少」。CIBERSORT 是目前被引最多的免疫浸润推断工具之一它不靠显微镜也不靠流式分选而是从批量转录组数据里用数学方法把 22 种免疫细胞的相对比例拆出来。这份配套资源包含完整可跑的 R 脚本、LM22 特征矩阵和示例数据属于「零基础能直接复现」的类型。它面向的是已经拿到表达矩阵、想做免疫浸润却无处下手的生信新手以及想快速验证一批样本免疫组成的科研人员。全文不涉及算法黑匣子的空谈直接告诉你每一步怎么选参数、怎么避坑。2. 免疫浸润分析的选型逻辑为什么 CIBERSORT 适合零基础上手2.1 免疫浸润分析的两条路线去卷积与富集打分免疫浸润分析在方法学上大体分两条路线。第一条是富集打分路线代表工具是 ssGSEA 和 xCell。这类方法把免疫细胞的标记基因集合看成一个个「特征库」然后拿样本的表达谱去和特征库做富集分析输出一个相对活性分数。它不要求所有标记基因的表达值都在一个量纲下对输入数据的格式宽容度较高。第二条是去卷积路线代表工具就是 CIBERSORT。它的逻辑更接近解方程把样本的表达谱看成 22 种免疫细胞表达谱的线性混合目标是求出一组比例系数让这组系数乘以各细胞类型的特征表达矩阵后能尽可能还原观察到的样本表达谱。对刚入门的人来说这两条路线最直观的差别在输出结果上。ssGSEA 给的是一个分数这个分数只在同一个数据集内部做相对比较才有意义而 CIBERSORT 输出的是百分比所有细胞类型的比例加起来接近 1。你从 CIBERSORT 结果里能直接说「这个样本里 CD8 T 细胞约占 12%」这是很多临床研究者习惯的表述方式。也正因为这个输出形式直观CIBERSORT 在肿瘤免疫相关的转录组分析里出场率极高。不过它的代价是输入矩阵的格式要求更严格。基因名必须是标准的 Gene Symbol表达值最好转换成 TPM而且对数据质量非常敏感——这些坑后面会一条条展开。我自己的习惯是如果只想快速看看趋势用 ssGSEA如果要做分组比较且最终图表要呈现比例关系就上 CIBERSORT。2.2 CIBERSORT 的核心机制ν-SVR 与 LM22 特征矩阵要跑好 CIBERSORT至少要理解它的两个核心组件。第一个是 LM22 特征矩阵。这是一个包含 22 种免疫细胞亚型、由 547 个标记基因构成的「标准表达谱」。这 547 个基因是作者从纯化细胞样本的芯片数据里筛出来的每种细胞类型在该矩阵里对应一组相对特异的表达模式。第二个是 ν-SVR全称是支持向量回归。CIBERSORT 没有用传统的最小二乘拟合而是用 ν-SVR 来做回归求解。为什么这么做因为免疫细胞的表达谱之间有很强的相关性普通线性回归在多重共线性下会给出非常不稳定的解ν-SVR 通过引入一个不敏感损失带允许一定的拟合误差换来的是比例解在小扰动下的稳定性。这两个组件决定了 CIBERSORT 的输入约束。LM22 的特征基因只有 547 个所以你的表达矩阵里必须能匹配到足够多的这些基因。如果匹配上的基因数太少算法会直接提示结果不可靠。另外LM22 建于 Affymetrix 芯片平台它记录的是绝对荧光信号强度这与 RNA-seq 的计数单位有天然差异。所以原版论文和后续实践都建议把 RNA-seq 数据先转换成 TPM 再输入而不是直接用 FPKM 或 read count。这一点很关键我在 3.1 节会给出具体的转换代码。很多新手在这里翻车——表达矩阵的基因名倒是匹配上了但因为量纲差异跑出来的比例分布明显偏离预期。提示CIBERSORT 的名字也意味着它的输出本质是「相对丰度」。22 种细胞的比例总和为 1某一种细胞的比例升高必然伴随另一种相对下降。这个数学约束会在免疫学解读时带来一些反直觉的结论后面第 5 章会专门讨论。2.3 输入数据的工程要求基因名与表达量单位把技术细节讲透前先给你一个我在实际数据上反复踩出来的结论CIBERSORT 对输入矩阵的要求可以归纳成「六个字」——去重、改名、转单位。去重指的是表达矩阵中不能有重复的基因名。很多从 GEO 下载的矩阵或自己比对产生的 count 矩阵里一个基因 Symbol 对应多行转录本或探针如果不处理就直接跑CIBERSORT 会报错让你无从下手。改名指的是确认你的矩阵行名是 Gene Symbol 而不是 Ensembl ID 或 Entrez ID。LM22 用的是 Gene Symbol如果你的矩阵是 Ensembl ID哪怕基因完全重合也匹配不上。转单位就是前面强调的 RNA-seq 数据要转 TPM。这三个要求里最容易忽略的是第一个。我见过不少同学拿到矩阵后直接丢进 CIBERSORT报错了才发现基因有重复。处理重复基因的标准做法是按基因取均值聚合而不是简单去重保留第一行。保留第一行会丢表达信息而且在后续可视化里会出现奇怪的离群样本。聚合操作我一般用两行aggregate就能完成具体代码在 3.1 节。你可以把这一小节当成一个检查清单——先把这三件事做完再谈跑算法否则后面的每一步都在为前面的马虎买单。3. 跑通第一次 CIBERSORT表达矩阵、LM22 与三个关键参数3.1 数据准备读入表达矩阵并做基因去重聚合配套资源里自带一个示例表达矩阵expression_matrix.txt为了让你理解格式我先说明它的结构第一列是 Gene Symbol其余列是样本名单元格里的数值是表达量。日常工作中你的矩阵可能是从 featureCounts、Salmon 或者 GEO 下载后整理而来的但进入 CIBERSORT 前的处理逻辑是一样的。下面这段代码完成读入与去重聚合# 读入表达矩阵第一列基因名作为行名 expr - read.table(expression_matrix.txt, header TRUE, row.names 1, sep \t, check.names FALSE) # 把基因名从行名提取为显式列方便 aggregate 按基因分组 expr$gene - rownames(expr) # 对同一基因的多行表达值取均值得到唯一基因行 expr_agg - aggregate(. ~ gene, data expr, FUN mean) # 恢复基因名为行名删除多余的 gene 列 rownames(expr_agg) - expr_agg$gene expr_agg$gene - NULL # 保留基因名与 LM22 特征矩阵的交集并输出基因匹配数 lm22 - read.table(LM22.txt, header TRUE, row.names 1, sep \t, check.names FALSE) common_genes - intersect(rownames(expr_agg), rownames(lm22)) cat(匹配到 LM22 的基因数, length(common_genes), \n) # 用交集基因重新过滤表达矩阵和 LM22保证两者基因顺序一致 expr_final - as.matrix(expr_agg[common_genes, ]) lm22_final - as.matrix(lm22[common_genes, ])代码逻辑上分成三块。第一块是读入并检查格式check.names FALSE是为了防止 R 把样本名里的横线或空格改成点号这个细节在后续匹配样本时很关键。第二块用aggregate按基因分组取均值它处理的是多对一映射比直接duplicated删除更稳妥因为取均值保留了该基因在多个转录本上的整体表达水平。第三块是求交集并过滤。这里我建议你务必打印length(common_genes)看一眼如果这个数字低于 300说明你的基因名版本或 ID 类型有问题建议回头检查转换步骤。我遇到过一次匹配数只有 200 出头的情况原因是表达矩阵用的是版本较旧的 Gene Symbol和 LM22 的命名有出入换用最新的注释文件后就恢复正常了。3.2 运行 CIBERSORT脚本调用与 perm、QN 的取舍数据准备完毕后进入正式运行环节。配套资源里的核心文件是cibersort.R它是一个封装好的函数库。使用时先source加载函数再调用CIBERSORT函数。下面是最小可运行的调用方式# 加载 CIBERSORT 官方函数 source(cibersort.R) # 调用主函数输入为 LM22 文件路径、表达矩阵、排列次数、是否做分位数归一化 res - CIBERSORT(sig_matrix LM22.txt, mixture expr_final, perm 1000, QN TRUE) # 查看结果的前几行和前几列确认输出格式 res[1:5, 1:5]这里的参数有三个需要重点解释。第一个是perm它代表置换次数CIBERSORT 用它来估计每个样本结果的显著性 P 值。默认值是 100但实践中 100 次置换得到的 P 值波动很大同一个样本跑两次可能一个显著一个不显著。我在正式分析里一般设 1000这会增加一些运行时间但对结果的稳定性提升非常明显。第二个是QN布尔值代表是否做分位数归一化。如果输入数据来自芯片平台设TRUE是合理选择如果输入的是 RNA-seq 的 TPM 数据是否设TRUE存在争议。原因是分位数归一化会把不同样本的表达分布强行拉齐这在芯片数据里是标准操作但对于 TPM 数据部分方法学文章认为它会引入虚假信号。我的经验是如果你的样本分组来自同一批次测序设TRUE问题不大如果样本来自不同平台或不同批次建议设FALSE并先自己用limma的normalizeBetweenArrays做过一次归一化再输入。第三个参数mixture它接收的是矩阵对象而不是文件路径所以必须先读入并处理成上一步的expr_final。3.3 结果文件解析三列统计量不是摆设运行结束后res是一个行为样本、列为细胞类型的矩阵但它的列并不全是比例值。最后三列依次是P-value、Correlation、RMSE。这三个统计量很重要很多人只看比例忽略了它们。P-value来自置换检验它回答的问题是「这个样本的免疫细胞比例估计是否显著优于随机」。如果某样本的 P 值大于 0.05说明该样本的 deconvolution 结果并不可靠后续做差异分析时建议谨慎纳入或直接剔除。Correlation表示用估计出的细胞比例重建表达谱后与原始表达谱的相关性。这个值越高说明拟合越好。RMSE是均方根误差越低越好。一个常见的筛查习惯是剔除Correlation小于 0.8 或P-value大于 0.05 的样本再进入下游分析和绘图。# 提取细胞比例部分去掉最后三列统计量 props - res[, 1:22] # 检查有多少样本 P 值大于 0.05决定是否需要剔除 pval_col - res[, P-value] cat(P值大于0.05的样本数, sum(pval_col 0.05), \n) # 如果只是看分组趋势可以先把 P 大于 0.05 的样本标记出来不急着删 res - as.data.frame(res) res$sample_id - rownames(res)这里我补充一个实际经验如果你跑出来的所有样本 P 值都小于 0.05但比例分布很不合理比如所有样本的巨噬细胞都是 0那大概率不是算法问题而是表达矩阵里对应 LM22 基因的表达值缺失严重。检验方法很简单把表达矩阵里那部分巨噬细胞标记基因单独提取出来看表达量分布。如果大量基因表达值都是 0说明你的上游定量流程或基因注释有问题跟 CIBERSORT 参数无关。先查数据再怀疑算法这是避免浪费半天时间的原则。4. 把结果画到能投稿堆叠条形图、热图与小提琴图4.1 堆叠条形图一眼看出每个样本的免疫组成CIBERSORT 输出的比例矩阵本身是一张表格但审稿人更愿意看图。最常用的第一张图是堆叠条形图每个样本一根柱子柱子的不同颜色段代表不同免疫细胞的占比。这张图能直观展示样本间的免疫组成差异。画图首选ggplot2但需要先把比例矩阵从宽格式转换成长格式这一步对 ggplot2 来说是必须的library(ggplot2) library(reshape2) # props 是第 3 节提取的 22 列比例矩阵先加一列样本名 props_df - as.data.frame(props) props_df$sample_id - rownames(props_df) # 宽表转长表每一行变成一个样本-细胞类型的组合 props_long - melt(props_df, id.vars sample_id, variable.name cell_type, value.name proportion) # 按样本排序并指定细胞类型在堆叠条中的顺序 props_long$sample_id - factor(props_long$sample_id, levels rownames(props_df)) # 堆叠条形图positionfill 让每根柱子高度都为 1 ggplot(props_long, aes(x sample_id, y proportion, fill cell_type)) geom_bar(stat identity, position fill, width 0.8) theme_minimal() theme(axis.text.x element_text(angle 60, hjust 1)) labs(x 样本, y 细胞比例, fill 细胞类型)长格式的转换是画图前的固定动作melt把 22 列比例值堆叠成一列id.vars指定保留样本名列。position fill让每根柱子等高展示的是组成结构而不是绝对含量如果你想让柱子的高度反映某个细胞类型分数的高低改用position stack即可。横轴样本标签比较多时可以旋转 60 度避免重叠。如果你希望图里每个样本按照某个细胞比例升序排列把levels改成按该细胞比例排序后的样本名向量就行。再有就是配色22 种细胞类型手动配色很麻烦直接用scale_fill_manual配一套有色差的色板别用默认的 hue 循环打印出来分不清相近色。4.2 热图与小提琴图分组比较的标准配置堆叠条形图适合展示单样本组成但当你需要比较肿瘤组和对照组时热图和小提琴图是更通用的选择。热图适合展示所有样本、所有细胞类型表达矩阵的整体格局通常配合样本聚类来观察是否存在亚型分组。小提琴图则适合对某一细胞类型做两组差异比较。下面的代码分别给出热图和小提琴图的最小实现library(pheatmap) # 热图对细胞比例做行标准化聚类展示格局 pheatmap(as.matrix(props), cluster_cols TRUE, cluster_rows TRUE, scale row, border_color NA, color colorRampPalette(c(navy, white, firebrick3))(100)) # 小提琴图以 CD8_T_cells 为例做两组比较 library(ggplot2) # 构造示例分组信息实际分析中替换为你的临床分组表 group_info - data.frame( sample_id rownames(props), group c(rep(Tumor, 10), rep(Normal, 10)) ) # 取目标细胞比例并合并分组 plot_df - data.frame( sample_id rownames(props), cd8 props[, T cells CD8], group group_info$group[match(rownames(props), group_info$sample_id)] ) ggplot(plot_df, aes(x group, y cd8, fill group)) geom_violin(trim FALSE) geom_boxplot(width 0.15, outlier.shape NA) geom_jitter(width 0.1, size 1.5, alpha 0.6) theme_classic() labs(x 分组, y CD8 T 细胞比例)热图的scale row会按细胞类型做标准化这会让颜色反映的是每种细胞在不同样本间的相对高低而不是原始比例。如果你的目的是直接比较比例大小把scale去掉用原始比例值配色更符合直觉。小提琴图里我刻意叠加了箱线和散点这个组合是生信论文里最常见的呈现方式。group_info里我假设前 10 个样本是肿瘤、后 10 个是正常实际使用时请按照你的样本注释表格改务必注意样本名与rownames(props)严格对应。另外要留意T cells CD8这个列名的写法不同版本的 LM22 命名略有差异小写、空格、下划线的细微差别会让取列时报错。稳妥做法是先colnames(props)看一眼实际列名再写代码。提示小提琴图背后的统计检验没有写进代码但正式分析时需要用wilcox.test(cd8 ~ group, data plot_df)算一下 P 值再标到图上。样本量小的时候 P 值容易不显著这是正常现象不要为了出图强行套 t 检验。5. 避坑记录从报错到错误解读的五个高频问题5.1 报错与冲突R 版本、包依赖和函数缺失现象source(cibersort.R)之后调用函数报错could not find function ...或者 R 直接提示加载某个包失败常见于 R 4.2 及以上版本。原因CIBERSORT 的 R 脚本发布年份较早内部依赖了部分在后续 R 版本中被移除的包和函数例如gplots中的某些绘图函数和preprocessCore的编译依赖。新装的环境里这些包要么装不上要么函数名变了。解决我一般会先按顺序装齐三件套preprocessCore从 Bioconductor 安装e1071和gplots从 CRAN 安装。如果仍报错不要盲目改源码先看错误信息定位到具体函数再决定是降级安装老版本包还是绕道走。一个比较省事的方案是直接在项目目录下新建一个renv环境固定 R 版本为 4.1.x。这块没有秘诀就是依赖管理问题。5.2 结果异常全零、基因匹配数过低与 P 值全为 1现象跑出来的比例矩阵中某个或某几个细胞类型在所有样本中全为 0或者匹配基因数不到 200另外有一类情况是所有样本 P 值都是 1。原因全为 0 通常不是算法问题而是表达矩阵里这一类细胞的标记基因几乎不表达或者基因名版本与 LM22 对不上。匹配数过低基本可以断定是 ID 类型问题最常见的是 Ensembl ID 直接当 Symbol 用。P 值全为 1 则是因为perm参数太小或表达矩阵中有极端离群值导致 SVR 求解不稳定。解决先做两步排查。第一步打印交集基因数如果低于 300优先怀疑 ID 映射用biomaRt或本地注释文件做转换。第二步检查表达矩阵是否存在全零行或包含巨大离群值的样本这类行会在 SVR 求解时干扰权重计算。我习惯在进入 CIBERSORT 之前按表达量过滤一次删除在所有样本中表达量都低于某个阈值的基因。至于 P 值问题把perm提到 1000并且检查是否有某列表达值明显超出其他列几个数量级有的话先做 log2 变换再跑。5.3 数据来源混用把 FPKM 直接灌进 CIBERSORT现象用 TCGA 的 FPKM 矩阵直接跑结果和已发表文献里的免疫浸润分布对不上CD4 和 CD8 T 细胞比例总是偏低或偏高得很离谱。原因CIBERSORT 的 LM22 模型是在芯片绝对信号强度上训练的RNA-seq 数据输入前需要把 FPKM 转换成 TPM。FPKM 与 TPM 虽然在很多样本里高度相关但 FPKM 的样本间可比性不如 TPM尤其在基因长度偏好明显的转录组数据里直接用 FPKM 会让 SVR 的拟合残差变大。解决跑之前统一做一次 FPKM 到 TPM 的转换公式是 TPM FPKM / sum(FPKM) * 1e6。这个转换有现成函数自己写三行 R 代码也能完成。第 6 章我会给出可直接复制的实现。另外如果你的上游是 STAR 定量出的 read count建议先用DESeq2的vst或edgeR的cpm做标准化再转 TPM不要直接拿原始 count 进来。这里的核心认知是CIBERSORT 需要的是「相对丰度」语义的表达值任何破坏样本间可比性的操作都会污染结果。5.4 置换检验的玄学perm 次数与结果可复现性现象同一份数据、同一个参数跑两次 CIBERSORT 得到的比例矩阵并不完全一致某些细胞类型的比例在小数点后第二位浮动P 值甚至会在 0.04 和 0.06 之间来回跳。原因perm置换次数本质上是蒙特卡洛过程次数越少随机性越大。默认值 100 次对大多数样本够用但边缘显著的样本就会非常不稳定。解决方法只有一句话——把perm调到 1000 或更高。代价是运行时间线性增长一个 50 样本的矩阵在普通笔记本上大约跑几十秒到几分钟。从那以后我每次做正式分析都固定用perm 1000并且在论文方法部分明确写清楚。如果你跑的是大规模队列比如几百个样本可以先把perm设 100 跑一遍过滤掉明显可靠的样本再对边缘样本重跑高 perm这样能省不少时间。5.5 解读陷阱相对丰度不是绝对含量现象分析报告中写「A 组的 CD8 T 细胞比例显著高于 B 组」但结合病理切片或流式数据却发现实际绝对数量没有差异甚至方向相反。原因CIBERSORT 输出的是相对丰度所有细胞比例总和固定为 1。某一种免疫细胞比例升高不代表它的绝对数量增加也可能是其他细胞类型减少导致的「分母变小」。这是去卷积方法天然带有的数学约束不是算法 bug。解决在论文表述里避免使用「浸润程度显著升高」这类暗示绝对量的措辞改用「相对比例变化」或「在免疫细胞组成中占比发生变化」。如果必须与绝对数量挂钩建议联合病理切片的免疫组化染色或流式定量数据做验证。这个坑在投稿时最容易被审稿人抓住要格外注意。我在一次项目中被审稿人质疑过这点后来把结果描述改写为「相对比例改变」并且补充了一组 IHC 验证数据才过关。注意如果你看到某个样本的 22 种细胞比例之和不是 1比如 0.93 或 1.07这不算 bug。CIBERSORT 做了归一化处理但小数位舍入可能导致轻微偏差。偏差过大超过 0.1时通常意味着表达矩阵里存在异常值。6. 进阶用法FPKM 转 TPM、批量运行与结果落地当单次运行已经熟练之后下一步值得做的事有三件把表达量单位转换脚本化、让 CIBERSORT 在多个数据集或多种参数下批量运行、以及把结果与下游分析衔接起来。关于 FPKM 转 TPM下面是固定套路# fpkm_matrix 是行为基因、列为样本的 FPKM 矩阵 fpkm_to_tpm - function(fpkm_matrix) { # 对每个样本每个基因的 FPKM 除以该样本 FPKM 总和再乘 1e6 tpm - t(t(fpkm_matrix) / colSums(fpkm_matrix) * 1e6) return(tpm) } tpm_matrix - fpkm_to_tpm(fpkm_matrix)这个函数内部做了两次转置第一次t(fpkm_matrix)把样本放到行上colSums后分母变成每个样本的总 FPKM外层再转置回来保持基因在行。函数式写法是一劳永逸的。我一般在拿到任何来源的表达矩阵时都会先确认它的单位然后第一时间转成 TPM再存入后续分析流程。批量运行的核心是循环封装。你可以把第 3 节的整个流程包成一个函数输入是表达矩阵文件路径和分组信息输出是比例矩阵和对应图。注意每次循环里都需要重新做基因交集过滤和source加载或者把cibersort.R的依赖包提前加载好放在循环外。批量运行时另一个实用技巧是把结果统一写入一个汇总目录文件名带上数据集名称和perm值防止自己把不同参数的结果搞混。至于与下游分析的衔接我推荐把 CIBERSORT 比例矩阵保存为 CSV 后直接读入 GSVA 的结果一起做相关分析——免疫浸润比例与通路活性矩阵的相关性热图是肿瘤转录组论文里出现频率很高的组合图。做法很简单cor_mat - cor(props, gsva_scores, method spearman)然后用pheatmap画相关性热图p.adjust做多重校正后筛选显著配对。最后想说的是整个过程中最值得养成的一个习惯是拿到任何表达矩阵后不要急着跑 CIBERSORT先在纸上列出基因名版本、表达量单位、是否有重复基因、LM22 匹配数这四项检查清单。从那次被审稿人质疑相对丰度解读之后我每做一个免疫浸润项目都会强制走一遍这个检查流程再顺手把 FPKM 转 TPM 的脚本跑一遍确认结果分布没有离群。这份配套资源里的示例数据和脚本足够你把整个流程从上游矩阵一路跑到出图但数据和参数替换成你自己的项目时上面的坑一个都省不掉。希望帮到你。本文还有配套的精品资源点击获取