ARTICLE DETAIL

资讯详情

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

Cibersort免疫浸润分析全流程:反卷积原理、参数设置与可视化

Cibersort免疫浸润分析全流程:反卷积原理、参数设置与可视化 这个工具我在几次项目里来回用过也帮人远程排查过不少报错。很多刚接触免疫浸润分析的朋友拿到表达矩阵就直接往网页版上传结果要么格式报错要么跑出来结果全是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_idCibersort会找不到基因名列直接报错。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之前先确认自己的数据类型是芯片还是测序。个人习惯是芯片数据用TRUERNA-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默认TRUERNA-seq建议FALSE这个前面已经提过。如果你不确定手里的RNA-seq数据应该怎么设置可以先分别用TRUE和FALSE跑一遍比较22种细胞类型比例的分布。如果某一种细胞类型在所有样本中占了80%以上多半是参数设置有问题。运行时间方面几十个样本配合perm1000大概需要几十分钟到一小时取决于机器性能。如果样本上百个建议留出足够时间或者分批运行再合并结果。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参数。默认参数QNTRUE听起来很合理但在RNA-seq数据上跑出来的结果非常夸张某一种巨噬细胞亚群在几乎每个样本里都占60%以上。改成QNFALSE之后比例分布才回到正常范围。第三类是Windows下使用旧版本R运行脚本时报障碍。当时我用的R版本是3.6脚本调用某个矩阵运算函数时会报错。解决办法很简单升级到R4.x并重新安装依赖包。这里建议遇到不明报错时优先检查R版本和依赖包而不是怀疑数据有问题。6.2 常见报错速查表下面这张表是我自己总结的高频报错和排查方向照着检查基本能解决大部分问题报错或异常可能原因排查方向输出结果全部为0或极低基因名不匹配、表达值类型错误检查基因名格式计算匹配到的LM22基因数报错提示找不到GeneSymbol第一列表头不是GeneSymbol修改表头确认第一列是基因名报错Error in read.table文件含引号或特殊字符分隔用标准tab分隔文本避免Excel直接另存报错subscript out of bounds表达矩阵含重复基因名或缺失值去重、处理缺失值后再运行运行时间过长样本多且perm1000确认不是死循环可适当降低perm并提前批量测试RNA-seq结果过度集中于某种细胞QN参数设置不当调QNFALSE并检查表达量是否为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表达量QNFALSEperm1000先把结果跑出来然后第一时间检查p值分布和CD8 T细胞等关键细胞类型是否与已有文献一致。只有这些基础检查通过了再去做那些花哨的可视化和进阶分析。
返回列表