CIBERSORT免疫浸润分析实战：从原理、数据准备到结果解读

简介：本资源是面向生物信息学零基础学习者的转录组下游分析实战配套材料，聚焦免疫微环境解析中的CIBERSORT算法应用，解决科研人员在肿瘤免疫浸润定量分析中常见的数据输入、R代码运行与结果解读难题。压缩包共13个文件，包含5个CSV格式的表达矩阵与分组数据、3个TXT格式的参考基因集（如LM22）与中间结果、2个已调试通过的R脚本（含主分析与差异检验）、1个PDF格式的可视化报告、1张PNG结果图及1个Rhistory操作记录，整体大小为144.83MB。已有239人下载学习，资源设计兼顾实操性与教学性：R脚本支持一键全选运行，输出涵盖CIBERSORT核心结果、Wilcoxon检验统计表及可视化图表，配套教程提供逐行注释与参数说明，便于新手理解免疫细胞比例估算原理与下游差异分析逻辑。 做转录组下游分析做到免疫浸润这一步的，一般都已经脱离“跑完流程拿到表达矩阵就收工”的阶段了。但你真去查CIBERSORT的资料会发现一个挺尴尬的局面：要么是英文原版文档，术语密集；要么是各种二手教程互相抄，中间缺了很多关键细节。尤其是“零基础”这个定位——很多教程默认你已经会处理表达矩阵、知道什么叫TPM、理解反卷积。可实际情况是，大多数刚接触转录组分析的人连CIBERSORT的输入文件到底长什么样都没概念。这篇我就按自己从零摸爬滚打的经验，把CIBERSORT从原理、数据准备、运行、结果解读到可视化、踩坑修复完整走一遍，配套资源文件怎么用、每个参数该怎么调，一次说清楚。

免疫浸润分析的输入就是你手头最普通的基因表达矩阵，输出却是每种免疫细胞在样本里的相对比例。这个转换过程听起来很黑盒，但核心逻辑其实不复杂。我在刚接触CIBERSORT的时候最大的困惑是：为什么一个基因表达矩阵能算出免疫细胞比例？这些比例到底代表什么？后来搞清楚原理之后，再看那些参数和报错就顺多了。这篇文章的思路就是先讲清楚原理，再带着你一步步跑通流程，最后把我在实际项目中遇到的各种坑和绕过的路全部列出来。

先说清楚适用范围：这篇文章适合手头已经有表达矩阵、想做免疫浸润分析但不太确定从哪下手的初学者；也适合那些已经跑过一遍但结果全都是0、或者不知道怎么看p值的同学。肿瘤样本的转录组数据分析中，免疫浸润比例是很多下游分析（比如生存分析、药物响应预测、肿瘤分型）的基础。CIBERSORT作为引用量极高的免疫反卷积工具，是入手免疫浸润最值得先掌握的一个。

1. 免疫浸润分析到底在做什么——先搞懂这步分析的意义

1.1 肿瘤微环境里的“细胞身份调查”

肿瘤组织不是一个只有肿瘤细胞的均质团块。它里面混着各种免疫细胞（T细胞、B细胞、巨噬细胞、NK细胞等）、基质细胞、成纤维细胞，这些东西组合在一起叫做肿瘤微环境（Tumor Microenvironment, TME）。不同病人肿瘤里的免疫细胞组成差异很大，这直接决定了免疫治疗的效果、预后好坏、甚至对化疗的敏感性。

传统做法是想办法把组织做成切片，用免疫组化或者流式细胞术去数细胞。但这些方法要么只能验证少数几种细胞，要么对样本质量要求极高，批量做几十上百个样本基本不现实。转录组测序给了一个更省事的思路：既然是组织测序，得到的表达信号就是所有细胞表达信号的混合。如果能从这个混合信号里把不同免疫细胞的比例“拆”出来，就可以在不做湿实验的情况下、批量估计免疫细胞组成。这就是免疫反卷积（deconvolution），也是CIBERSORT这类工具存在的意义。

这里要强调一个关键认知：CIBERSORT算出来的比例不是细胞绝对数量，而是“某种免疫细胞在所有免疫细胞中的相对丰度”。因为模型假设样本里的表达信号全部来自免疫细胞，每一行的结果加起来约等于1。理解这一点很关键，因为它决定了你后续怎么解读结果——比如一个样本M2巨噬细胞是0.25，意思是免疫细胞里约25%是M2巨噬细胞，而不是整个组织里M2巨噬细胞占25%。

1.2 CIBERSORT与同类工具的定位区别

做免疫浸润分析的工具有很多，新手经常被绕晕。我把常见的几类放在一起对比，你就明白为什么很多文献首选CIBERSORT，以及它适合什么场景不适什么场景。

CIBERSORT的核心优势有两个：一是用支持向量回归处理高维表达数据，对基因共线性的容忍度高；二是输出带有置换检验p值，告诉你每个样本的估算结果靠不靠谱。这两个特点其实是后来同类工具普遍借鉴的标配——你想，一个没有可靠度指示的结果，拿去做生存分析或者组间差异，万一全是噪声怎么办？

那CIBERSORT的短板也明显：特征矩阵LM22是基于芯片数据构建的，而且只针对人的免疫细胞，细胞类型固定为22种。你用RNA-seq数据跑，原则上要做相应的数据转换和参数调整；你拿小鼠样本跑，虽然可以通过同源基因转换硬跑，但结果可信度打折扣。这些细节我后面都会展开说。

2. CIBERSORT的算法原理与LM22——理解核心概念才能用好结果

2.1 基因表达反卷积：把混合信号拆成单组分

要理解CIBERSORT到底在干什么，可以拿一个生活化的例子类比。假设你录了一段音频，里面同时有钢琴、吉他和人声，你想知道这三种声音各占多大比例。如果你有一份“每种乐器的标准音色特征库”，就可以拿混合音频去对照，反推出每种乐器的强度。CIBERSORT做的事情本质上一样：组织测序的表达矩阵是“混合音频”，LM22特征矩阵就是“每种免疫细胞的标准表达特征库”，最后通过算法反推出每种细胞在混合信号里的“强度”，也就是比例。

数学上，CIBERSORT用支持向量回归来求解这个问题。为什么不是简单的线性回归？因为免疫细胞比例有两个天然约束：比例不能为负、总和要接近1。线性回归解出来可能会出现负值，这在生物学上完全说不通。支持向量回归可以加入这些约束条件，并且在高维基因空间中做回归时更稳健——你想，输入的特征不是两三个基因，而是LM22里的几百个标记基因，这么多基因之间还有共线性，普通回归很容易被个别离群基因带偏。

CIBERSORT的求解流程大致是：把某个样本的基因表达向量作为因变量，LM22里各种细胞类型的表达特征作为自变量，用支持向量回归拟合，得到一个回归系数向量，再经过归一化处理，就是各种免疫细胞的比例估计。这个过程对每个样本独立进行，所以样本之间互不影响。

这里有个实操上很重要的理解：CIBERSORT实际上“混合矩阵里的基因必须和LM22基因对齐”。你的表达矩阵里基因名再全，如果和LM22匹配不上的话，那些基因就被丢弃了。很多新手第一次跑出来全是0，八成就是这个匹配环节出了问题。所以数据准备阶段最关键的任务，就是保证基因名格式和LM22一致。

2.2 LM22特征矩阵的适用范围与限制

LM22是CIBERSORT作者基于公开的免疫细胞表达数据集构建的特征矩阵，包含了22种人类免疫细胞亚型，每种细胞用一组特征基因来代表，总共涉及547个基因。这22种细胞包括：初始B细胞、记忆B细胞、浆细胞、CD8 T细胞、CD4初始T细胞、CD4记忆静息T细胞、CD4记忆激活T细胞、滤泡辅助T细胞、调节性T细胞（Tregs）、γδ T细胞、静息NK细胞、激活NK细胞、单核细胞、M0巨噬细胞、M1巨噬细胞、M2巨噬细胞、静息树突状细胞、激活树突状细胞、静息肥大细胞、激活肥大细胞、嗜酸性粒细胞、中性粒细胞。

这个矩阵的适用范围有三个必须记住的限制。第一，它是基于人源细胞构建的，理论上只适用于人的转录组数据。第二，它基于芯片表达谱构建，处理RNA-seq数据时需要谨慎，尤其是不要直接套用芯片时代的归一化参数。第三，它只包含免疫细胞类型，实体瘤组织里大量存在的肿瘤细胞、基质细胞并不在模型内。这就导致一个潜在问题：如果你的组织样本里免疫细胞占比很低，而不是很纯的免疫细胞环境，CIBERSORT会把这些非免疫细胞的表达信号强行“分配”到22种免疫细胞头上，结果可能失真。

在实操上怎么处理这个限制？我自己的习惯是：实体瘤样本尽量别单独依赖CIBERSORT一个结果，至少再用一种不依赖LM22的方法（比如ssGSEA或xCell）做交叉验证；如果样本是纯化的免疫细胞群或者免疫细胞高度富集的微环境，CIBERSORT的参考价值就大得多。另外，后来的CIBERSORTx允许用户自定义特征矩阵，可以把肿瘤细胞、基质细胞也纳入模型，这才是正解。

2.3 为什么要有permutation和p值

CIBERSORT输出结果里有一列P-value，这一列特别容易被误解。它不是组间差异分析的p值，而是“这个样本的反卷积结果是否可信”的置换检验p值。具体逻辑是：为了评估你估算出的免疫细胞比例是否显著优于随机结果，算法会把样本的表达值随机打乱，重新用同样的流程跑一遍反卷积，重复很多次（permutation次数，通常设为1000次）。如果真实数据得到的拟合效果显著好于随机打乱后的拟合效果，p值就小；如果打乱后也能得到类似的拟合效果，说明你得到的结果很可能是噪声，p值就大。

实际操作中，我建议把P-value、Correlation、RMSE三个指标一起看。P-value 小于0.05说明反卷积结果比随机好；Correlation 是真实表达谱与拟合表达谱的相关性，越高越好，通常大于0.8比较理想；RMSE是拟合误差的均方根，越小越好。这三个指标刻画的是同一个问题：“用LM22和估算出的细胞比例，能不能很好地重建出你观测到的表达谱？”如果不能，那估算出的细胞比例就不太可信。

常见错误是有人跑完结果直接把22列比例拿去分析，完全不管p值。如果一个样本p值是0.36，那这个样本的免疫浸润比例基本就是个随机结果，放进下游分析会污染整个统计。正确做法是：先过滤掉p值大于0.05（有人更严格用0.01）的样本，或者至少做一个敏感性分析，证明剔除低质量样本后主要结论不变。

3. 零基础上手第一步：数据准备是成败关键

3.1 输入表达矩阵长什么样

CIBERSORT的输入文件是纯文本格式的基因表达矩阵，行是基因，列是样本。它不认Excel格式，也不认R的数据框格式——原版脚本要求传文件路径，脚本内部用read.table读取。矩阵长这样：

GeneSymbol Sample1 Sample2 Sample3 TP53 152.3 98.4 203.7 EGFR 45.2 67.8 12.9 CD8A 0.3 12.5 78.2 ...

第一行是表头，第一列是基因符号（Gene Symbol），不能有重复；列名是样本名，尽量不要有特殊字符（空格、括号、中文都会出问题）。第一列基因名最好只有基因名称这一个字段——很多人从Ensembl或者NCBI下载注释文件，直接把带版本号的或者带染色体位置的基因名塞进第一列，这就不行。LM22用的是标准的官方基因Symbol，比如CD8A、GZMB、CD274这种。

如果你的表达矩阵行名是Ensembl ID，跑之前必须转成Symbol。转换工具有很多：biomaRt、clusterProfiler的bitr函数、org.Hs.eg.db包都可以。转换之后要注意：同一个Symbol可能对应多个Ensembl ID，需要去重。去重策略我一般取表达量均值最大的那条（代表主要转录本），或者直接取平均值。网上有人取最大值，有人取平均，我个人倾向取平均，因为这样对多转录本基因更公平。

还有一个小但非常重要的问题：所有表达值必须是数值，不能有NA，不能有负数。有缺失值的行直接删除或者补0（但补0会影响后面的归一化，所以首选删除）。有负值说明你的数据经过了log2处理且有负的fold change，这种情况在芯片数据里偶尔出现，需要小心处理——不是说不能跑，但要意识到这些负值会影响后续的归一化和相关性计算。

3.2 环境准备：R版本与依赖包

CIBERSORT原版脚本是用R写的，所以先确保你的电脑有R环境。我推荐R 4.2以上版本，Windows、Mac、Linux都行。需要用到的包有三个：e1071（支持向量机）、parallel（并行计算，加快置换检验）、preprocessCore（分位数归一化，从Bioconductor安装）。

依赖包的安装其实是个容易卡住的点，尤其是preprocessCore。e1071和parallel直接用install.packages就能装，preprocessCore是Bioconductor的包，在Windows上用install.packages装不了。正确安装方式是：

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("preprocessCore", update = FALSE, ask = FALSE)

Linux服务器上如果编译报错，通常需要先装系统依赖，比如libcurl、libxml2这些。Windows上则需要配置Rtools。第一次装这个包的时候我在这上面耗了整整一下午，后来发现一个更省事的方法：直接用conda建一个R环境，然后conda install -c bioconda r-preprocesscore，编译问题全没了。如果你是在服务器上跑，这条路确实省心。

3.3 官方源代码怎么跑（原版CIBERSORT.R）

官方提供了两套运行方式：一个是斯坦福实验室维护的网页版工具，图形界面，上传表达矩阵就可以跑，适合不求甚解快速出结果的场景；另一个是原版R脚本CIBERSORT.R，需要下载到本地，然后source进R会话里调用。R脚本的优势是可重复、可批量、可集成到自动化流程，而且能自己控制permutation次数和归一化参数。这里我重点讲R脚本方式，因为网页版反而有很多限制，比如样本数上限、无法自定义参数，对后续批量分析不友好。

拿到CIBERSORT.R和LM22.txt两个文件后，把它们放在你的工作目录下。标准调用方式：

source("CIBERSORT.R") result <- CIBERSORT("LM22.txt", "my_expression.txt", perm = 1000, QN = TRUE)

这个函数有三个核心参数：perm是置换检验次数，建议1000，要更严格的检验可以设2000，但计算时间会成倍增加；QN是是否做分位数归一化，TRUE或FALSE，这个参数选错是很多人踩坑的重灾区，我单独拿出来讲。函数会先做基因匹配，然后逐样本反卷积，最后返回一个数据框（实际上是矩阵），每一行是一个样本，前22列是免疫细胞比例，最后三列是P-value、Correlation、RMSE。

关于QN参数，你需要记住这条分界线：芯片数据用TRUE，RNA-seq数据用FALSE。为什么？芯片数据的探针信号强度在不同芯片间存在系统性分布差异，分位数归一化可以校正这种技术偏差，让不同芯片之间可比。而RNA-seq数据（尤其是TPM、FPKM这类经过文库大小和基因长度校正后的数据）本身的分布形态与芯片完全不同，再做分位数归一化会把真实生物学差异也抹掉一部分，容易出现大量样本p值变大、细胞比例分布异常的情况。我实测过同一个RNA-seq数据集，QN=TRUE跑出来有一半样本p值大于0.05，QN=FALSE跑出来大部分样本都小于0.05。差距就是这么明显。

4. 从表达矩阵到免疫浸润结果：完整运行过程

4.1 用示例数据跑通的完整流程

为了让你能跟着完整操作一遍，我准备了一个模拟的演示过程。假设你有一个表达矩阵文件exp_TPM.txt，列是样本名（10个肿瘤样本Tumor_01~Tumor_10，10个正常样本Normal_01~Normal_10），行是基因Symbol。第一步先读进来看看数据概况：

# 读取表达矩阵 expr <- read.table("exp_TPM.txt", header = TRUE, row.names = 1, sep = "\t", check.names = FALSE) # 查看维度 dim(expr) # 应该类似 [1] 18000 20 # 查看前几行 head(expr[, 1:4]) # 检查是否有NA sum(is.na(expr)) # 删除有缺失值的行 expr <- na.omit(expr) # 删除全为0的行（完全没表达信息的基因没有意义） expr <- expr[rowSums(expr > 0) > 0, ] # 处理重复基因名：取平均值 expr <- aggregate(expr, by = list(rownames(expr)), FUN = mean) rownames(expr) <- expr$Group.1 expr <- expr[, -1] # 写出CIBERSORT要求的纯文本格式 write.table(expr, "exp_TPM_clean.txt", sep = "\t", quote = FALSE, col.names = TRUE, row.names = TRUE)

这里有几个操作细节值得说一下。read.table里的check.names=FALSE非常重要——如果不设这个参数，R会自动把样本名里的横杠“-”变成点“.”，后面你拿结果去跟临床数据匹配时会莫名其妙对不上。write.table里quote=FALSE也别忘了，否则生成的文本文件里基因名和样本名都会被加上引号，CIBERSORT读入后可能报错。

接下来正式运行CIBERSORT：

source("CIBERSORT.R") set.seed(123) # 固定随机种子，保证置换检验结果可重复 result <- CIBERSORT("LM22.txt", "exp_TPM_clean.txt", perm = 1000, QN = FALSE) # 查看结果结构 dim(result) # 20行 × 25列（22种细胞 + 3个统计量） head(result[, 22:25]) # 保存结果 write.csv(result, "cibersort_result.csv")

这个set.seed很多人会忽略。由于置换检验涉及随机抽样，如果不固定随机种子，每次运行的结果会有一点点不同（细胞比例的主结果基本稳定，但p值会有轻微波动）。为了可重复性，强烈建议设置种子。

4.2 结果文件里每一列代表什么

CIBERSORT返回的对象有25列：前22列就是LM22定义的22种免疫细胞亚型，列名直接是细胞类型名称（比如B cells naive、T cells CD8、Macrophages M0等）；第23列是P-value，第24列是Correlation，第25列是RMSE。

第一次拿到结果的人最常犯的错，就是看到22列比例后直接做标准化或者归一化。注意，CIBERSORT输出的比例已经是归一化好的——每一行的22列加起来应该接近1（因为数值精度可能略微偏离，但基本在0.99到1.00之间）。你不需要再做任何转换，直接可以用于下游分析。

关于三列质量控制指标，我给一个可操作的筛选建议：

我的习惯是先看P-value，剔除不合格样本；再看Correlation，如果整体偏低就要怀疑LM22是否适用于这批数据。如果大部分样本的Correlation都在0.6左右，即使p值显著，结果解释也要非常谨慎——很有可能你的样本组织构成和LM22假设相差太远。

4.3 低质量样本过滤与细胞比例矩阵整理

跑完CIBERSORT后，第一步不是画图，而是过滤。筛选p < 0.05的样本：

# 过滤低质量样本 pass <- result[result[, "P-value"] < 0.05, ] cat("通过质量控制的样本数:", nrow(pass), "/", nrow(result), "\n") # 只保留22种细胞比例 prop_matrix <- pass[, 1:22]

如果过滤后样本数太少（比如只剩3个），下游分析基本没法做。这时候要回头检查数据：基因匹配率是否太低？QN参数是否选错？表达矩阵是否经过不恰当的标准化？有时候样本量少是正常的，比如你本来就只有10个样本，过滤掉2个还剩8个，可以接受。但如果50个样本过滤完只剩10个，一定是数据准备环节出了问题。

这里分享一个我的个人习惯：在正式分析之前，先跑一遍perm=100做快速测试，看看整体的p值分布情况。如果连perm=100都有一大半样本p值大于0.05，那就别急着上perm=1000，先排查数据问题。perm=1000比perm=100慢10倍，40个样本可能要跑20多分钟，等跑完发现结果全不合格再回头改数据更浪费时间。

4.4 从22种细胞类型到生物学结论

22种细胞类型很多，但在实际解读时不需要逐个看。我的做法是先归纳成几大类：T细胞（CD8 T、CD4 T各亚群、Treg、γδ T）、B细胞（naive、memory、plasma）、NK细胞（resting、activated）、髓系细胞（单核细胞、巨噬细胞M0/M1/M2、树突状细胞、肥大细胞）、粒细胞（中性粒细胞、嗜酸性粒细胞）。然后在分组比较时先看大类趋势，再下钻到具体亚群。

举例来说：肿瘤样本相比正常样本，通常CD8 T细胞和M1巨噬细胞有变化，Treg和M2巨噬细胞会增多（反映免疫抑制微环境）。这些生物学趋势是判断结果是否合理的一个“直觉校验”。如果跑出来的结果完全反常识，比如正常组织中CD8 T细胞比例异常高、所有肿瘤样本Treg都接近0，那大概率是数据或者参数出了问题，而不是生物学真的这样。

5. 跑通之后的结果可视化——堆叠图、热图、箱线图实战

5.1 免疫浸润比例的堆叠条形图

可视化的第一步，通常是用堆叠条形图展示每个样本的22种细胞组成。这个图能直观地看出不同组之间免疫构成的总体差异。

library(ggplot2) library(tidyr) # 准备绘图数据 plot_df <- as.data.frame(prop_matrix) plot_df$sample <- rownames(plot_df) plot_df$group <- ifelse(grepl("Tumor", plot_df$sample), "Tumor", "Normal") # 宽表转长表 plot_long <- pivot_longer(plot_df, cols = all_of(colnames(prop_matrix)), names_to = "cell_type", values_to = "proportion") # 按组排序样本 plot_df$sample <- factor(plot_df$sample, levels = c(paste0("Normal_", 1:10), paste0("Tumor_", 1:10))) ggplot(plot_long, aes(x = sample, y = proportion, fill = cell_type)) + geom_bar(stat = "identity", width = 0.8) + scale_y_continuous(labels = scales::percent) + labs(x = "Sample", y = "Relative Proportion", fill = "Cell Type") + theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8), legend.position = "bottom") + guides(fill = guide_legend(ncol = 3))

这个图出来之后，先检查一个基本特征：每个柱子（样本）的22段加起来应该等于100%。如果有的柱子的总高度明显低于1，说明有样本的质量指标不过关（虽然你已经过滤了，但还是核对一下）。如果所有柱子都正常，接下来就按照分组观察：Tumor组和Normal组在哪些细胞组分上有明显差异。堆叠图的视觉效果比较强，但精确定量还是得靠箱线图。

5.2 免疫细胞相关性热图

第二个常用可视化是免疫细胞两两之间的相关性热图。这个图能帮你看出哪些细胞倾向于共同出现（正相关）或互相排斥（负相关）。生物学上，Treg和CD8 T细胞经常呈负相关，因为Treg抑制CD8 T细胞活性；M1和M2巨噬细胞往往此消彼长。

library(corrplot) # 计算Spearman相关矩阵 cor_mat <- cor(prop_matrix, method = "spearman") # 绘制热图 corrplot(cor_mat, method = "color", type = "upper", tl.col = "black", tl.cex = 0.8, col = colorRampPalette(c("#4575B4", "white", "#D73027"))(200), addCoef.col = "black", number.cex = 0.6)

相关性热图还可以结合聚类分析，把细胞类型分成几个模块。如果你的数据里有明显的免疫抑制微环境信号，通常能看到Treg、M2巨噬细胞聚在一个模块里，与CD8 T细胞模块负相关。这种发现往往是文章里的一个小亮点，比单纯列比例更有信息量。

5.3 分组差异的箱线图与统计检验

后续分析中最常见的需求是：两组（比如肿瘤vs正常，或者响应vs不响应）之间，哪些免疫细胞比例有显著差异。这时候用箱线图加Wilcoxon检验（两组比较）。

library(ggpubr) # 以CD8 T细胞为例 cd8_df <- data.frame( value = prop_matrix[, "T cells CD8"], group = ifelse(grepl("Tumor", rownames(prop_matrix)), "Tumor", "Normal") ) p <- ggboxplot(cd8_df, x = "group", y = "value", fill = "group", palette = c("#00AFBB", "#E7B800"), add = "jitter", shape = "group") + stat_compare_means(method = "wilcox.test", label = "p.format") + labs(y = "CD8 T cell Proportion", title = "CD8 T cell: Tumor vs Normal") print(p) ggsave("CD8_boxplot.pdf", p, width = 4, height = 5)

但有个重要的统计学细节：22种细胞同时做差异检验，必须做多重假设检验校正。常用的有BH校正（Benjamini-Hochberg），也就是FDR。如果不校正，22次检验里随机可能就有1-2个假阳性，你拿来写进文章容易被审稿人质疑。实际操作中，我先把22种细胞的p值都算出来，然后统一做BH校正，再画图。

# 对所有22种细胞做循环检验 p_values <- sapply(colnames(prop_matrix), function(cell) { df <- data.frame( value = prop_matrix[, cell], group = ifelse(grepl("Tumor", rownames(prop_matrix)), "Tumor", "Normal") ) wilcox.test(value ~ group, data = df)$p.value }) # BH校正 p_adjust <- p.adjust(p_values, method = "BH") # 查看显著差异的细胞 significant <- names(p_adjust[p_adjust < 0.05]) print(significant)

顺便说一句，箱线图加散点（jitter）是常规操作，但样本量很小时（比如每组只有3个样本），散点比箱线图更能传达样本分布的信息。我见过太多只画箱线图不画点的图，样本量小的时候箱子根本没有统计意义，还是把每个样本的点标出来实在。

5.4 与临床信息/其他指标的关联分析（简要）

免疫浸润比例的终极价值不在描述本身，而在于关联分析。最经典的是三联分析：免疫细胞比例与肿瘤分期、分级的关系；免疫细胞比例与生存预后的关系（KM曲线）；免疫细胞比例与免疫治疗响应或其他分子标志物的关系。

生存分析的切入点一般是：用某种免疫细胞比例的中位数把样本分成高、低两组，然后做Kaplan-Meier曲线和log-rank检验，或者进一步用Cox回归校正临床变量。这个分析用survival和survminer包就能完成：

library(survival) library(survminer) # 假设clinical里有OS.time（生存时间）和OS.event（生存状态） # data_merge是细胞比例和临床信息合并后的数据框 data_merge$cd8_group <- ifelse(data_merge$`T cells CD8` > median(data_merge$`T cells CD8`), "High", "Low") fit <- survfit(Surv(OS.time, OS.event) ~ cd8_group, data = data_merge) ggsurvplot(fit, data = data_merge, pval = TRUE, risk.table = TRUE)

这个延伸方向一篇足够写一大章，这里点到为止。但记住一个原则：细胞比例与临床指标关联分析之前，一定要确保前面的质量控制是严格的，否则你的“发现”很可能只是技术误差的体现。

6. 全程踩坑记录——那些让我多花一星期的错误

6.1 不检查基因名格式，结果全为0

这是我第一次跑CIBERSORT时踩的坑。当时我用的表达矩阵从某个数据库下载，行名是Ensembl ID加版本号（比如ENSG00000141510.17），完全没转换成Symbol就直接跑了。结果CIBERSORT输出一大堆0，个别细胞类型全是0，我当时还以为算法出了问题。后来发现是基因匹配环节几乎全部失败——CIBERSORT内部只保留了与LM22基因名完全匹配的行，Ensembl ID和Symbol根本对不上。

排查方法其实很简单：运行CIBERSORT之后，控制台会输出匹配到的基因数量。如果你发现匹配到的基因不到100个，而LM22本身有547个基因，说明你的基因名格式有严重问题。我现在的习惯是：在跑CIBERSORT之前先单独写一段代码，检查表达矩阵的行名与LM22基因的交集数量：

lm22 <- read.table("LM22.txt", header = TRUE, row.names = 1, check.names = FALSE) matched <- intersect(rownames(expr), rownames(lm22)) cat("匹配到的基因数量:", length(matched), "/", nrow(lm22), "\n")

如果匹配率低于80%，就别跑了，先处理基因名转换。用biomaRt转换Ensembl ID到Symbol是最常见的做法，clusterProfiler的bitr函数也很方便。

6.2 对RNA-seq数据用了QN=TRUE，出现大批p>0.05

这个坑我在前面已经预告过，再强调一遍因为太典型了。有一次我帮同事跑一批TCGA的RNA-seq数据，她用的教程里写着QN=TRUE（因为教程作者主要处理芯片数据），结果跑出来一半样本p值大于0.05，Correlation也很低。我当时查了好久，最后把QN改成FALSE，一切都正常了。

为什么RNA-seq不能随便用分位数归一化？核心原因是RNA-seq数据的分布形态和芯片数据差异很大。芯片数据的信号强度有上限，背景噪声水平大致一致，不同芯片间的分位数分布差异主要来自技术因素；RNA-seq的count或者TPM数据是稀疏的、高度偏态的，而且不同样本的文库组成差异本身就携带生物学信息。分位数归一化会强制所有样本的分布完全一致，这等于把一部分真实的生物学差异也抹掉了。所以用RNA-seq数据时QN=FALSE，用芯片数据时QN=TRUE，这是一个经验法则。

当然也有特殊情况：RNA-seq数据经过TMM或DESeq2的vst/rlog标准化后，分布已经相对一致，这时候有人还是会用QN=TRUE做额外校正。我的建议是不要依赖特殊情况，对新手来说，RNA-seq就选FALSE，简单安全。

6.3 样本名带特殊字符导致列错乱

有次我从某平台下载数据，样本名长这样：“TCGA-XX-XXXX-01A-11R-XXXX-07”。read.table读进来没问题，但CIBERSORT内部处理时样本名带有“-”有时候会出现解析问题。更麻烦的是遇到样本名里有空格、括号、中文，R会默默地把列名改成合法形式（空格变点、中文变乱码）。等跑完拿结果去跟临床数据merge的时候才发现对不上号。

解决方法是读入表达矩阵时务必用check.names=FALSE，写文本文件时保持列名原样。如果样本名确实包含特殊字符，最好在读入之后统一改一次名，比如把横杠换成下划线，然后再写出去。这样CIBERSORT运行时样本名干净，后续匹配也不会出错。

6.4 preprocessCore包安装失败

Windows用户安装preprocessCore的痛，估计每个跑过CIBERSORT的人都懂。这个包从Bioconductor安装，但在Windows上经常需要编译，而编译又依赖Rtools。有一个情况特别坑：你装了R 4.3，但Rtools还是老版本，编译直接报错。

我的几个可行方案：第一种，确保Rtools版本和R版本匹配，然后重新install；第二种，用conda创建环境安装r-preprocesscore，在Windows上可以用WSL或者直接装Linux虚拟机，但这对新手太折腾；第三种，如果只是临时用，可以先不装preprocessCore也能跑CIBERSORT——只要QN=FALSE，脚本内部不会调用分位数归一化函数，那个包其实用不到。这句话很关键：如果你用RNA-seq数据且QN=FALSE，preprocessCore缺失不影响运行；但如果你要用芯片数据QN=TRUE，那必须装好它。

我之前在一个Linux服务器上遇到过一个更奇怪的问题：BiocManager::install报错提示“package ‘preprocessCore’ is not available for this version of R”。后来发现是服务器的R版本太老，Bioconductor的旧版本仓库已经不维护了。最后我升级了R版本才解决。所以如果你也遇到这种提示，先检查R版本是不是太旧。

6.5 运行慢的真相：perm=1000没你想的那么轻松

CIBERSORT的置换检验非常耗时间。perm=1000意味着每个样本要做1000次SVR训练，样本越多越慢。40个样本perm=1000在个人电脑上可能要跑15到30分钟，这还得看CPU性能。如果你有几百个样本，那基本要按小时算了。

有个提高效率的办法：明确分配的核心数。CIBERSORT.R脚本默认调用detectCores()来检测可用核心，但有时会检测到很多逻辑核心，导致并行开销反而拖慢速度。我一般会在脚本里手动设置：

# CIBERSORT.R中修改并行核心数 num_cores <- 4 # 最多不要超过物理核心数 cl <- makeCluster(num_cores)

如果你的机器有8核16线程，设4到8都行。如果是在共享服务器上跑，注意别把全部核心占满，别人还要跑任务呢。

不过我也要提醒一句：perm值直接影响p值精度。perm=100和perm=1000跑出来的主结果（细胞比例）几乎一致，但p值会有差异，perm太小时p值分辨率很低。正式分析还是用1000，这是文献共识。

6.6 LM22里的细胞注释与文献中的命名差异

最后一个小坑是命名问题。LM22里的细胞类型名称是缩写形式，比如“T cells CD4 naive”“Macrophages M0”，写文章时通常需要展开成完整的生物学名称。很多新手直接把列名复制到文章里，读起来很别扭。我建议在图表和文章中统一用规范名称：比如“CD4+ naive T cells”“M0 macrophages”。最保险的做法是参考已发表文献里的命名方式，保持全文一致。

还有一个容易被审稿人质疑的点：LM22里的细胞类型注释实际上是基于“亚群相似性”的估算，比如“T cells CD4 memory activated”并不完全等同于流式分选出来的那一群细胞。所以在方法部分描述时，直接写“CIBERSORT algorithm with the LM22 signature matrix was used to estimate the relative proportions of 22 immune cell subtypes”，不要过度承诺精确的细胞身份。

7. 进阶方向与替代工具——别让CIBERSORT成为你的终点

7.1 CIBERSORTx与更灵活的特征矩阵

如果你做的是RNA-seq数据，并且对CIBERSORT默认的LM22特征矩阵不满意，升级到CIBERSORTx是更优的选择。CIBERSORTx是原团队推出的增强版，有几个显著的改进：支持bulk RNA-seq数据的分位数归一化优化、支持批次效应校正、最关键的是支持用户用单细胞RNA-seq数据构建自定义特征矩阵。

这意味着什么？如果你手头正好有同一批样本的单细胞数据，可以先从单细胞数据里鉴定细胞类型，构建出这个组织类型特有的特征矩阵，再反卷积回去做bulk数据的免疫浸润估算。相比通用的LM22，这种“定制化特征矩阵”的准确度会明显提高，因为LM22毕竟是基于外周血和免疫细胞亚群构建的，对特定实体瘤组织未必完全适配。

使用CIBERSORTx需要去官网注册账号，把表达矩阵上传到服务器运行。要注意：数据上传涉及隐私，如果你处理的是未发表的数据，务必确认自己的权限。

7.2 immunedeconv包：多算法横向对比

R里有个immunedeconv包，把CIBERSORT、EPIC、QUANTISEQ、TIMER、MCPcounter等多套免疫反卷积方法统一封装在一起，输入同一个表达矩阵，一键输出多种算法的结果。这个包最有用的场景是验证稳健性。

审稿人常问的一个问题：你只用CIBERSORT一种算法，结果可靠吗？如果你能用两三种算法得到一致的结论，说服力会强很多。immunedeconv就是干这个的：

# 安装 # install.packages("immunedeconv") library(immunedeconv) # 一次性计算多种方法的结果 res_mcp <- deconvolute(expr, method = "mcp_counter") res_quantiseq <- deconvolute(expr, method = "quantiseq")

注意一点：immunedeconv内嵌的CIBERSORT方法输出没有置换检验p值，它调用的是另一种封装方式。所以如果你需要p值做样本过滤，还是用原版脚本；如果你只是想多算法交叉验证，用immunedeconv就够了。

7.3 多算法交叉验证的基本思路

做交叉验证时，有一件事要先说清楚：不同算法输出的细胞类型粒度不一样。CIBERSORT输出22种，MCPcounter输出8种，TIMER输出6种，直接比数值没意义。应该比的是“趋势一致性”：比如CD8 T细胞富集的方向是否一致、巨噬细胞M2相关的免疫抑制信号是否在同样的分组中出现。

实际操作上，我会做一个简单的相关性表：把CIBERSORT的细胞大类（比如CD8 T细胞）与MCPcounter的对应细胞类型（比如CD8 T cell）在所有样本中的比例/分数做Spearman相关。如果相关系数高、方向一致，说明不同算法得到的结论是互相印证的。这种分析在文章里虽然只占一小段，但对增强结果的可信度非常有帮助。

7.4 单细胞时代还需要CIBERSORT吗

单细胞测序已经很普及了，有人会问：既然单细胞能直接数细胞，为什么还要CIBERSORT做反卷积？这个问题值得认真回答。单细胞测序成本高、技术门槛高，大规模临床队列（几百上千例）做单细胞不现实。但转录组测序（RNA-seq或者芯片）在临床队列里几乎成了标配数据。CIBERSORT的作用，就是用少量单细胞数据或者公共单细胞资源作为参考，去推断大规模转录组队列里的免疫浸润组成。本质上它充当了一座桥——从高分辨率的单细胞参考，推到大样本的bulk数据上。

所以CIBERSORT不会因为单细胞技术的发展而失去意义，反而因为单细胞参考矩阵的丰富而变得更强大。CIBERSORTx自定义特征矩阵的路径，实际上就是把单细胞数据转化为免疫浸润参考资源的标准方法。

8. 配套资源与个人经验

说了这么多，最后把配套资源的使用方式交代清楚。通常情况下，一套完整的CIBERSORT配套资源包含以下文件：CIBERSORT.R（主脚本）、LM22.txt（特征矩阵）、示例表达矩阵（demo数据）、示例输出结果。新手拿到这些文件后，我建议按照这样一个顺序来学习：先用示例数据跑通一遍，确认环境和代码都没问题；然后替换成自己的表达矩阵，先跑perm=100快速测试；确认结果合理后再用perm=1000正式运行；最后再做可视化和下游分析。

我自己在实际操作中形成的一些个人习惯，分享出来供你参考。第一，每次都固定set.seed，保证可重复性；第二，任何一步操作前先备份原始数据，尤其是做基因名转换和去重时，原始矩阵一旦被覆盖就很难恢复；第三，跑完正式结果后第一时间检查质量控制指标，而不是急着画图——画一个精美的图出来才发现数据不合格，那才是真的浪费时间；第四，所有中间文件（表达矩阵、匹配基因列表、过滤后样本列表）都保留，后续如果审稿人质疑结果，你可以快速定位问题。

如果让我给一个刚入门的同学推荐最快的上手路径，我的建议是：不要先纠结算法公式，先拿着示例数据跑通一遍，再把自己的数据换进去，遇到报错再回头查资料。CIBERSORT这个工具最大的特点就是“只要数据格式对、参数选对，结果就是稳的”；反过来，数据准备稍微偷懒一点，后面问题重重。把这篇文章里的几个关键检查点做完，你的结果大概率不会太差。踩过那么多坑之后，我最大的体会是：在整个免疫浸润分析流程里，真正的瓶颈极少是算法本身，而是你对输入数据、参数选择和输出质量控制的敬畏程度。数据干净、流程严谨，好结果是水到渠成的事；数据粗糙、参数乱调，再牛的算法也救不了你。

本文还有配套的精品资源，点击获取