ARTICLE DETAIL

资讯详情

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

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

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支持向量回归(SVR)反卷积基因表达矩阵22种免疫细胞比例经典引用量高需LM22特征矩阵支持置换检验CIBERSORTx改进版反卷积表达矩阵支持构建自定义特征矩阵自定义细胞类型比例官方网页版支持RNA-seq批次校正更灵活ssGSEA单样本基因集富集分析表达矩阵细胞类型富集分数相对值不能直接算比例但适合分组比较xCellssGSEA扩展表达矩阵64种免疫/基质细胞分数覆盖细胞类型多适合做广谱筛查MCPcounter标记基因统计表达矩阵8种免疫/基质细胞分数简单快速对芯片和测序都稳TIMER线性回归校正表达矩阵6种免疫细胞浸润TCGA数据友好但细胞类型少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参数你需要记住这条分界线芯片数据用TRUERNA-seq数据用FALSE。为什么芯片数据的探针信号强度在不同芯片间存在系统性分布差异分位数归一化可以校正这种技术偏差让不同芯片之间可比。而RNA-seq数据尤其是TPM、FPKM这类经过文库大小和基因长度校正后的数据本身的分布形态与芯片完全不同再做分位数归一化会把真实生物学差异也抹掉一部分容易出现大量样本p值变大、细胞比例分布异常的情况。我实测过同一个RNA-seq数据集QNTRUE跑出来有一半样本p值大于0.05QNFALSE跑出来大部分样本都小于0.05。差距就是这么明显。4. 从表达矩阵到免疫浸润结果完整运行过程4.1 用示例数据跑通的完整流程为了让你能跟着完整操作一遍我准备了一个模拟的演示过程。假设你有一个表达矩阵文件exp_TPM.txt列是样本名10个肿瘤样本Tumor_01~Tumor_1010个正常样本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.namesFALSE非常重要——如果不设这个参数R会自动把样本名里的横杠“-”变成点“.”后面你拿结果去跟临床数据匹配时会莫名其妙对不上。write.table里quoteFALSE也别忘了否则生成的文本文件里基因名和样本名都会被加上引号CIBERSORT读入后可能报错。接下来正式运行CIBERSORTsource(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 0.05反卷积结果显著优于随机置换Correlation 0.8用估算比例重建的表达谱与真实表达谱高度相关RMSE越小越好重建表达谱与真实表达谱的误差我的习惯是先看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个一定是数据准备环节出了问题。这里分享一个我的个人习惯在正式分析之前先跑一遍perm100做快速测试看看整体的p值分布情况。如果连perm100都有一大半样本p值大于0.05那就别急着上perm1000先排查数据问题。perm1000比perm100慢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数据用了QNTRUE出现大批p0.05这个坑我在前面已经预告过再强调一遍因为太典型了。有一次我帮同事跑一批TCGA的RNA-seq数据她用的教程里写着QNTRUE因为教程作者主要处理芯片数据结果跑出来一半样本p值大于0.05Correlation也很低。我当时查了好久最后把QN改成FALSE一切都正常了。为什么RNA-seq不能随便用分位数归一化核心原因是RNA-seq数据的分布形态和芯片数据差异很大。芯片数据的信号强度有上限背景噪声水平大致一致不同芯片间的分位数分布差异主要来自技术因素RNA-seq的count或者TPM数据是稀疏的、高度偏态的而且不同样本的文库组成差异本身就携带生物学信息。分位数归一化会强制所有样本的分布完全一致这等于把一部分真实的生物学差异也抹掉了。所以用RNA-seq数据时QNFALSE用芯片数据时QNTRUE这是一个经验法则。当然也有特殊情况RNA-seq数据经过TMM或DESeq2的vst/rlog标准化后分布已经相对一致这时候有人还是会用QNTRUE做额外校正。我的建议是不要依赖特殊情况对新手来说RNA-seq就选FALSE简单安全。6.3 样本名带特殊字符导致列错乱有次我从某平台下载数据样本名长这样“TCGA-XX-XXXX-01A-11R-XXXX-07”。read.table读进来没问题但CIBERSORT内部处理时样本名带有“-”有时候会出现解析问题。更麻烦的是遇到样本名里有空格、括号、中文R会默默地把列名改成合法形式空格变点、中文变乱码。等跑完拿结果去跟临床数据merge的时候才发现对不上号。解决方法是读入表达矩阵时务必用check.namesFALSE写文本文件时保持列名原样。如果样本名确实包含特殊字符最好在读入之后统一改一次名比如把横杠换成下划线然后再写出去。这样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——只要QNFALSE脚本内部不会调用分位数归一化函数那个包其实用不到。这句话很关键如果你用RNA-seq数据且QNFALSEpreprocessCore缺失不影响运行但如果你要用芯片数据QNTRUE那必须装好它。我之前在一个Linux服务器上遇到过一个更奇怪的问题BiocManager::install报错提示“package ‘preprocessCore’ is not available for this version of R”。后来发现是服务器的R版本太老Bioconductor的旧版本仓库已经不维护了。最后我升级了R版本才解决。所以如果你也遇到这种提示先检查R版本是不是太旧。6.5 运行慢的真相perm1000没你想的那么轻松CIBERSORT的置换检验非常耗时间。perm1000意味着每个样本要做1000次SVR训练样本越多越慢。40个样本perm1000在个人电脑上可能要跑15到30分钟这还得看CPU性能。如果你有几百个样本那基本要按小时算了。有个提高效率的办法明确分配的核心数。CIBERSORT.R脚本默认调用detectCores()来检测可用核心但有时会检测到很多逻辑核心导致并行开销反而拖慢速度。我一般会在脚本里手动设置# CIBERSORT.R中修改并行核心数 num_cores - 4 # 最多不要超过物理核心数 cl - makeCluster(num_cores)如果你的机器有8核16线程设4到8都行。如果是在共享服务器上跑注意别把全部核心占满别人还要跑任务呢。不过我也要提醒一句perm值直接影响p值精度。perm100和perm1000跑出来的主结果细胞比例几乎一致但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数据、示例输出结果。新手拿到这些文件后我建议按照这样一个顺序来学习先用示例数据跑通一遍确认环境和代码都没问题然后替换成自己的表达矩阵先跑perm100快速测试确认结果合理后再用perm1000正式运行最后再做可视化和下游分析。我自己在实际操作中形成的一些个人习惯分享出来供你参考。第一每次都固定set.seed保证可重复性第二任何一步操作前先备份原始数据尤其是做基因名转换和去重时原始矩阵一旦被覆盖就很难恢复第三跑完正式结果后第一时间检查质量控制指标而不是急着画图——画一个精美的图出来才发现数据不合格那才是真的浪费时间第四所有中间文件表达矩阵、匹配基因列表、过滤后样本列表都保留后续如果审稿人质疑结果你可以快速定位问题。如果让我给一个刚入门的同学推荐最快的上手路径我的建议是不要先纠结算法公式先拿着示例数据跑通一遍再把自己的数据换进去遇到报错再回头查资料。CIBERSORT这个工具最大的特点就是“只要数据格式对、参数选对结果就是稳的”反过来数据准备稍微偷懒一点后面问题重重。把这篇文章里的几个关键检查点做完你的结果大概率不会太差。踩过那么多坑之后我最大的体会是在整个免疫浸润分析流程里真正的瓶颈极少是算法本身而是你对输入数据、参数选择和输出质量控制的敬畏程度。数据干净、流程严谨好结果是水到渠成的事数据粗糙、参数乱调再牛的算法也救不了你。本文还有配套的精品资源点击获取
返回列表