☰
Cibersort免疫浸润分析全流程:反卷积原理、参数设置与可视化
2026/10/3 1:17:16 网站建设 项目流程

这个工具我在几次项目里来回用过,也帮人远程排查过不少报错。很多刚接触免疫浸润分析的朋友,拿到表达矩阵就直接往网页版上传,结果要么格式报错,要么跑出来结果全是0,要么根本不知道输出表里那一堆英文列名是什么意思。这篇教程我按自己的实操流程写,尽量把原理、输入格式、本地运行、结果筛选和可视化一次讲透,帮你少走几个月的弯路。

先说我自己的使用结论:如果你手里是常规的bulk转录组数据(芯片或RNA-seq),想快速知道样本里22种免疫细胞亚群的大致构成,Cibersort依然是目前最稳妥、审稿人最容易接受的方案之一。它输出的是相对比例,不是绝对数量,理解这一点很重要。下面我从算法核心开始讲。

1. 免疫细胞反卷积的本质:Cibersort到底在算什么

1.1 一锅汤里的食材构成:从混合表达谱反推细胞比例

组织测序拿到的那张表达矩阵,本质上是一锅煮好的汤。肿瘤组织里有肿瘤细胞、成纤维细胞、各种免疫细胞,测序仪把所有细胞的RNA混在一起读出平均值,你根本分不清哪个基因来自哪种细胞。反卷积要解决的问题,就是像品尝一锅汤判断里面放了多少种食材一样,根据特征基因的表达信号,估算出每种免疫细胞在混合样本中的占比。

当你输入一个肿瘤样本的表达谱,Cibersort会告诉你:CD8 T细胞占12%,巨噬细胞M0占8%,静息NK细胞占3%……这些数字就是免疫浸润分析最常被引用的结果。整个计算过程是一种解混合问题,本质上是已知混合物中某些标志物(marker基因)的信号强度,反推各组分比例。

1.2 LM22特征矩阵:22种免疫细胞的标准指纹库

Cibersort的核心依赖一个固定的签名矩阵,叫做LM22。这个矩阵由论文作者经过大量纯化细胞表达谱训练得到,包含22种成熟免疫细胞亚群的547个特征基因的表达特征,可以理解为每个细胞类型都有一个标准指纹库。

LM22里的22种细胞具体包括:naive和memory B细胞、浆细胞、CD8 T细胞、naive/静息记忆/激活记忆CD4 T细胞、滤泡辅助T细胞、调节性T细胞、γδ T细胞、静息和激活NK细胞、单核细胞、M0/M1/M2巨噬细胞、静息和激活树突状细胞、静息和激活肥大细胞、嗜酸性粒细胞、中性粒细胞。

实际使用中特别要注意,LM22矩阵不建议随意修改或替换。有些同学会把里面的基因名按自己的习惯改格式,结果反而破坏了签名矩阵的基因对关系,跑出来的比例明显异常。如果你想用自己数据构建新的签名矩阵,那是Cibersortx做的事情,原版Cibersort没这个功能。

1.3 SVR回归+1000次排列检验:算法核心与p值来源

Cibersort的算法核心是支持向量回归(SVR),它和常见的线性回归思路不同,更擅长处理高维特征。在分析单个样本时,算法会提取该样本中LM22包含的特征基因表达量,对22种免疫细胞分别做回归拟合,得到一个初始估计值,再通过归一化把结果换算成相对比例(所有细胞类型加起来通常接近1)。

算法在每个样本上独立建模,样本之间互不影响,这是Cibersort的一个特点。另一个特点是它内置了蒙特卡洛排列检验,官方推荐perm参数设为1000,也就是针对每个样本随机打乱基因标签后重新跑1000次反卷积,生成一个“随机结果”的分布,再比较你的真实结果落在分布的什么位置,最终得到一个p值。

这个p值很容易被误解。它不是某个细胞类型比例差异的显著性,而是该样本反卷积整体结果相对于随机情况的可信程度。p值越小,说明你的样本结果越不像随机碰出来的,结果越值得信任。后面我会专门讲怎么用这个p值做质量筛选。

2. 输入数据准备:90%的报错都发生在这一关

2.1 官方表达矩阵格式逐列拆解

我见过太多人在输入格式上栽跟头,其实官方对格式的要求非常明确。Cibersort要求的表达矩阵是:第一列必须是基因名,表头必须是GeneSymbol这个单词,注意大小写;第二列开始每一列是一个样本,列名是样本ID;每一行是一个基因。

文件要求是tab分隔的纯文本,UTF-8编码更保险。很多人在Excel里整理数据后直接另存为txt,结果Excel自动加了引号或者改变了分隔符,上传后解析全乱。最稳妥的做法是用R或Python脚本读写,不要用Excel做最后一步。

有个细节容易被忽略:表达矩阵里不要有额外的注释列、基因描述列、或者行名列(即第一列表头不能是空的)。如果第一列表头写成了Gene或者Gene_id,Cibersort会找不到基因名列,直接报错。

2.2 基因名格式问题:大小写、版本号、Ensembl ID与重复基因

基因名格式是重灾区。LM22里用的是gene symbol,比如CD3D、CD8A这种。如果你手里的表达矩阵是Ensembl ID格式,比如ENSG00000167286,直接拿去跑,匹配到的特征基因可能只有几十个,剩下全部对不上,结果自然不可信。

我的做法是在输入之前先做一步基因名转换。R里用org.Hs.eg.db包可以完成从Ensembl ID到symbol的映射,也可以用clusterProfiler里的转换函数。处理代码大致是这样:

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("org.Hs.eg.db") BiocManager::install("clusterProfiler") library(clusterProfiler) library(org.Hs.eg.db) expr <- read.table("expression_ensembl.txt", header = TRUE, row.names = 1, sep = "\t", check.names = FALSE) # 假设行名是Ensembl ID,去掉版本号 rownames(expr) <- gsub("\\..*", "", rownames(expr)) gene_map <- bitr(rownames(expr), fromType = "ENSEMBL", toType = "SYMBOL", OrgDb = org.Hs.eg.db) expr_symbol <- merge(gene_map, expr, by.x = "ENSEMBL", by.y = "row.names")

另外,有些数据库导出的基因symbol后面带着版本号,比如TP53.1,需要先去掉。还有重复基因名问题:不同探针可能对应同一个基因,直接去重的话会随机丢掉表达量更高的那个,建议按行平均值排序后保留每个基因平均表达量最高的那一行,再执行去重:

expr <- expr[order(rowMeans(expr), decreasing = TRUE), ] expr <- expr[!duplicated(rownames(expr)), ] write.table(expr, "expression_clean.txt", sep = "\t", quote = FALSE)

2.3 RNA-seq与芯片数据在参数上的差异化处理

Cibersort最初是基于microarray数据训练和验证的,但这几年大家手里的数据基本是RNA-seq,因此参数设置上要区分处理。

对于芯片数据,用RMA或GC-RMA标准化后的表达矩阵即可,官方推荐QN参数保持TRUE,因为QN会执行分位数标准化,可以在一定程度上消除不同平台间的系统性差异,提升跨平台可比性。

对于RNA-seq数据,建议把表达量转换为TPM形式,不要使用raw count,不建议再取log2。运行时的QN参数我一般设置为FALSE。分位数标准化在这里可能会过度压缩表达值差异,导致结果偏向某种细胞类型,比例分布非常不自然。这一点是很多教程没有讲清楚的,实测影响很大。

2.4 缺失值、样本量和批次效应的一般建议

表达矩阵中尽量不要有NA。如果有少量缺失,可以先做填补,或者删除缺失严重的基因。缺失值过多会直接影响特征基因的匹配和回归过程。

样本量方面,Cibersort对单个样本逐样本计算,理论上少到几个样本也能跑,但从统计角度看,样本量太少时p值很容易极端,组间差异也缺乏说服力。我建议至少保证每个分组有5个以上样本,总数尽量不少于20个,否则后续差异比较和生存分析很难出可靠结论。

还有一个经常被忽视的点:不要在跑Cibersort之前对表达矩阵做跨样本的批次效应校正,比如用limma的removeBatchEffect或者ComBat。反卷积算法依赖的是基因间的相对表达特征,批次校正会改变这些关系,反而可能破坏免疫信号。批次效应的影响可以放在结果解读阶段讨论,不要预处理阶段就消除掉。

3. 本地运行全流程:从脚本下载到结果输出的完整操作

3.1 准备两个核心文件

本地运行需要两个文件:CIBERSORT.R和LM22.txt。这两个文件可以在Cibersort官网注册后从下载页面获取。LM22.txt就是前面说的22种免疫细胞特征基因矩阵,CIBERSORT.R是官方提供的R脚本。

下载后我习惯把这两个文件连同表达矩阵放在同一个工作目录下。文件路径最好全部用英文,不要带中文和空格,Windows环境下因为路径问题导致的source失败和读取失败我已经见过太多次了。

3.2 运行Cibersort函数的完整R代码

本地运行的核心代码其实非常简洁。先把工作目录设置好,然后source脚本,再调用CIBERSORT函数:

setwd("D:/project/cibersort") source("CIBERSORT.R") res <- CIBERSORT("LM22.txt", "expression_clean.txt", perm = 1000, QN = TRUE)

第一次跑的时候建议先设置QN = TRUE或者FALSE之前先确认自己的数据类型是芯片还是测序。个人习惯是芯片数据用TRUE,RNA-seq数据用FALSE,这一点也可以先小规模试跑对比一下比例分布是否合理。

运行结束后把结果保存到文件:

write.table(res, "Cibersort_result.txt", sep = "\t", quote = FALSE, col.names = NA)

如果之前没有安装过需要依赖的包,运行过程中可能会提示缺包,按提示安装即可。另外,CIBERSORT.R脚本使用了一些Matrix包里的功能,确保R版本在4.0以上,旧版本在部分函数调用上会报错。

3.3 perm与QN参数选择说明

参数方面,perm控制排列检验次数,官方推荐1000,这是保证p值稳定性的通用选择。如果想先跑通流程、看结果大概长相,可以把perm设为0,这样运行速度非常快,但输出矩阵里没有p值。正式分析时一定要用1000,否则审稿人问到p值来源你解释不清楚。

QN参数的含义是是否执行跨平台分位数标准化。microarray默认TRUE,RNA-seq建议FALSE,这个前面已经提过。如果你不确定手里的RNA-seq数据应该怎么设置,可以先分别用TRUE和FALSE跑一遍,比较22种细胞类型比例的分布。如果某一种细胞类型在所有样本中占了80%以上,多半是参数设置有问题。

运行时间方面,几十个样本配合perm=1000大概需要几十分钟到一小时,取决于机器性能。如果样本上百个,建议留出足够时间,或者分批运行再合并结果。

3.4 输出文件结构

跑完之后res是一个矩阵,行是每个样本,前22列对应LM22里的22种免疫细胞类型,每一行的22个值加起来接近1。最后三列分别是P-value、Correlation和RMSE。

这个输出结构就是你后续做差异比较、可视化和生存分析的原材料。需要注意列名里包含空格,比如“T cells CD8”,在R里取列时要用反引号或字符串索引,否则会报错。

4. 结果可信度先过三关:p值、相关性、RMSE怎么看

4.1 三个指标分别说明什么

很多同学拿到输出表后直接开始画图,这是不严谨的。Cibersort在每个样本后面附了三个质量指标,应该是结果解读的第一道关卡。

P-value来自前面说的蒙特卡洛排列检验。它衡量的是当前样本的反卷积结果和随机打乱后的结果之间差异的显著性,p值越小,说明结果越不像随机造成的。我的筛选习惯是优先保留p值小于0.05的样本,p值大于0.1的样本结果会打一个比较明显的问号。

Correlation是拟合相关性,越接近1说明表达谱特征与LM22签名矩阵的整体拟合程度越好,结果越可信。

RMSE是均方根误差,越小说明拟合误差越小。

实际操作中我不会机械地用单一阈值卡死,而是三个指标放在一起看:相关性在0.8以上、RMSE相对较低、p值小于0.05,基本可以放心使用;如果p值合格但相关性偏低,结果解读时就要谨慎。

4.2 哪些样本结果最不可信

根据我的经验,有几类样本的结果特别容易出现质量指标异常。第一类是纯化程度很高的肿瘤样本,肿瘤细胞占比过高、免疫细胞信号太微弱,反卷积结果容易被噪声主导。第二类是RNA质量较差的样本,基因检出率低,匹配到的特征基因数量骤减。第三类是小样本量队列里的个别离群样本,p值和表达分布容易出现极端值。

遇到质量指标不合格的样本,我一般不会直接把整列删掉,而是根据分析目的决定。如果做的是组间差异比较,可以剔除掉明显离群的样本;如果做的是相关性分析,保留但标注异常值也是可接受的。关键是整个过程要在论文方法部分写清楚,包括筛选标准和剔除数量。

4.3 从结果里提取分组信息的常规操作

做完质量筛选后,下一步是把22列比例数据拆出来,加上你自己定义的分组信息(比如肿瘤组vs对照组、高复发组vs低复发组),整理成一张分析用表。这一步看似简单,但到后面画图时会节省大量时间。

# 去掉质量指标列,只保留22种细胞 cell_prop <- res[, 1:22] # 加上样本ID和分组 analysis_df <- data.frame( Sample = rownames(cell_prop), Group = group_info$Group, # 按样本名匹配分组 cell_prop ) write.table(analysis_df, "analysis_table.txt", sep = "\t", quote = FALSE, row.names = FALSE)

5. 从结果表格到论文级可视化:四类常用图表的R实现

5.1 堆叠柱状图:展示每个样本的免疫构成

免疫浸润分析最常用的总览图是堆叠柱状图,每个柱子代表一个样本,柱子里不同颜色代表不同免疫细胞类型的占比。这种图适合快速展示所有样本的整体免疫景观变化。

library(ggplot2) library(reshape2) cell_prop <- res[, 1:22] cell_prop$Sample <- rownames(cell_prop) plot_data <- melt(cell_prop, id.vars = "Sample") ggplot(plot_data, aes(x = Sample, y = value, fill = variable)) + geom_bar(stat = "identity", width = 0.8) + labs(x = "", y = "Relative percentage", fill = "Cell type") + theme_bw(base_size = 12) + theme(axis.text.x = element_text(angle = 45, hjust = 1))

样本名带下划线或太长时,可能需要调整旋转角度和图片宽度。

5.2 分组箱线图与差异检验:找到真正有变化的细胞类型

总览图之后,读者更关心的是不同分组之间哪些免疫细胞有显著差异。箱线图加差异检验是论文里出镜率最高的图。

library(ggpubr) plot_df <- data.frame( Group = group_info$Group, CD8T = res[, "T cells CD8"], Macrophage_M0 = res[, "Macrophages M0"] ) ggboxplot(plot_df, x = "Group", y = "CD8T", color = "Group", palette = "jco", add = "jitter") + stat_compare_means(method = "wilcox.test")

如果分组数量超过两组,可以考虑用Kruskal-Wallis检验,组内两两比较再加注释。要注意的是,由于一次比较22种细胞类型,存在多重检验的问题,我这里不会强行推荐FDR校正,但至少有意识地报告p值未经校正或说明校正方式,审稿人会更容易接受。

5.3 免疫细胞占比热图:看细胞之间的相关关系

热图在免疫浸润分析中的作用有两种。一种是按样本聚类展示22种细胞在不同样本间的丰度模式;另一种是计算细胞类型之间的相关性矩阵。

第一种热图用pheatmap直接画:

library(pheatmap) pheatmap(cell_prop, scale = "column", clustering_method = "ward.D2", show_colnames = FALSE, fontsize_row = 8)

对列做scale之后,每一列都变成了均值为0、标准差为1的Z-score,颜色表达的是相对高低而不是原始比例。

第二种相关性热图更多人用来观察免疫细胞间的协同或排斥关系,比如CD8 T细胞与M1巨噬细胞在多个研究中呈正相关,这就是有生物学意义的发现线索。

5.4 与临床特征或生存数据关联的简单扩展

当免疫浸润比例算完之后,最常见的下游延伸就是把某一种细胞比例的中位数作为分组依据,做生存曲线。例如筛选出预后相关的免疫细胞类型,然后按CD8 T细胞比例高低把患者分为两组,再用survival包做KM曲线和log-rank检验。

这一步不需要Cibersort参与,只需要把res里某一列的比例和临床生存数据按样本ID匹配起来就可以。需要注意中位分组只是探索性分析,不能替代多因素cox回归,投稿时方法部分写清楚是单因素探索即可。

6. 实战踩坑记录与工具选型建议:Cibersort的边界在哪里

6.1 我遇到过的三类典型问题

第一类是基因名格式问题。之前处理一批公开的RNA-seq数据,下载下来全是Ensembl ID,直接把原始表达矩阵拿去跑Cibersort,结果超过一半样本的p值是1,我一开始以为数据有问题,检查后发现问题出在基因匹配率只有不到10%。转换完基因名后,结果立刻变得正常。

第二类是RNA-seq数据没有调整QN参数。默认参数QN=TRUE听起来很合理,但在RNA-seq数据上跑出来的结果非常夸张,某一种巨噬细胞亚群在几乎每个样本里都占60%以上。改成QN=FALSE之后比例分布才回到正常范围。

第三类是Windows下使用旧版本R运行脚本时报障碍。当时我用的R版本是3.6,脚本调用某个矩阵运算函数时会报错。解决办法很简单,升级到R4.x并重新安装依赖包。这里建议遇到不明报错时优先检查R版本和依赖包,而不是怀疑数据有问题。

6.2 常见报错速查表

下面这张表是我自己总结的高频报错和排查方向,照着检查基本能解决大部分问题:

报错或异常可能原因排查方向
输出结果全部为0或极低基因名不匹配、表达值类型错误检查基因名格式,计算匹配到的LM22基因数
报错提示找不到GeneSymbol第一列表头不是GeneSymbol修改表头,确认第一列是基因名
报错Error in read.table文件含引号或特殊字符分隔用标准tab分隔文本,避免Excel直接另存
报错subscript out of bounds表达矩阵含重复基因名或缺失值去重、处理缺失值后再运行
运行时间过长样本多且perm=1000确认不是死循环,可适当降低perm并提前批量测试
RNA-seq结果过度集中于某种细胞QN参数设置不当调QN=FALSE并检查表达量是否为TPM

运行前可以先用一小段代码检查LM22基因在表达矩阵中的匹配率:

lm22 <- read.table("LM22.txt", header = TRUE, row.names = 1, check.names = FALSE) expr <- read.table("expression_clean.txt", header = TRUE, row.names = 1, check.names = FALSE) overlap <- length(intersect(rownames(lm22), rownames(expr))) cat("matched genes:", overlap, "/", nrow(lm22), "\n")

如果匹配上的基因数量低于400,建议回头认真检查基因名格式,这基本可以解释大部分结果异常。

6.3 Cibersortx、ssGSEA、xCell和MCPcounter如何选

Cibersort是免疫浸润领域的老牌工具,但不是唯一选择。根据不同场景,有些工具可能更合适。

如果手里的数据是单细胞转录组,或者你想研究LM22之外的细胞类型,比如特定组织里的罕见细胞亚群,可以考虑Cibersortx。它是Cibersort的升级版,支持基于单细胞数据构建自定义签名矩阵,也支持绝对模式输出,但需要在官网注册使用,且免费额度有限。

如果只是想快速比较组间免疫细胞活性,不追求精确比例,ssGSEA或xCell会更省事。ssGSEA基于基因集富集分数,输出的是活性得分而不是细胞比例,interpretation上更接近通路活性。xCell在此基础上扩展到了60多种细胞类型,覆盖面更广。

MCPcounter则是基于marker基因表达均值的打分方法,简单快速,适合做大样本筛选。TIMER对TCGA数据友好,内置了一些可视化模块,但对独立数据集的支持相对受限。

工具输出类型特色适合场景
Cibersort22种免疫细胞相对比例算法经典,审稿接受度高bulk转录组的免疫比例估算
Cibersortx自定义细胞类型比例/绝对丰度支持单细胞构建签名矩阵需要自定义细胞类型或绝对丰度
ssGSEA基因集活性得分快速,无需反卷积组间活性比较
xCell64种细胞丰度得分覆盖细胞类型多需要更广细胞谱
MCPcounter相对丰度得分简单稳健大样本初筛
TIMER6种免疫细胞浸润深度整合TCGATCGA数据快速探索

从我个人的使用经验来看,除非有明确的场景升级需求,否则常规免疫浸润分析直接用Cibersort就够了。它输出的是相对比例,适合做样本间横向比较,但不适合跨数据集直接对比绝对值。跑完Cibersort之后,建议再看一下匹配率和质量指标,画图之前先把不可信的样本剔除掉,这样后续分析才站得住脚。

如果你手里正好是RNA-seq数据,我最后的建议是:用TPM表达量,QN=FALSE,perm=1000,先把结果跑出来,然后第一时间检查p值分布和CD8 T细胞等关键细胞类型是否与已有文献一致。只有这些基础检查通过了,再去做那些花哨的可视化和进阶分析。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询