简介:IOBR免疫浸润分析R包代码,是一套聚焦肿瘤微环境解析的开源分析工具,面向从事生物信息学与肿瘤免疫研究的科研人员。包内集成CIBERSORT、TIMER、xCell等8种主流免疫浸润算法,并收录264个基因集签名,覆盖代谢、m6A等领域,支持PCA、z-score、ssGSEA三种打分方式,可直接用于批量分析与结果可视化。压缩包共18个文件,以R脚本、Rdata数据集为主,辅以CSV表格、PNG结果图、HTML报告及README说明,整体大小9.53MB,结构与函数划分清晰,便于按需调用。其中包含可运行的示例分析流程与模拟数据生成脚本,可帮助读者快速复现热图等图形输出,理解从数据加载到浸润量化、单细胞配受体活性网络分析乃至肿瘤微环境聚类分群的全链路操作。目前已有213人学习下载,适合希望借助现成代码减少踩坑、提升肿瘤微环境分析效率的R语言使用者。
1. IOBR免疫浸润分析R包:为什么我把它放进转录组分析的默认流水线
IOBR这个免疫浸润分析R包,是近两年我在转录组项目里用得最顺手的工具之一。它要解决的事情很具体:从表达矩阵推断肿瘤微环境里有哪几类免疫细胞、比例是多少、跟临床结局有什么关系。过去我得挨个装六七个包,数据格式还得来回转,现在一个接口全部跑完。它的定位不是某种新算法,而是把CIBERSORT、EPIC、quanTIseq、MCPcounter、TIMER、xCell这些成熟方法收拢成统一的输入输出格式。新手跟内置示例走一遍能跑通,熟手直接用自定义表达矩阵替换即可,它适合处理肿瘤转录组、免疫治疗队列,以及需要多种算法互相印证的场景。下面我把从零到一跑通这套免疫浸润分析流程的实操记录拆开讲。
2. IOBR安装与数据准备:依赖坑与表达矩阵的入场检查
2.1 两种安装方式与依赖处理
IOBR目前主要通过Github分发,CRAN上没有收录,所以安装路径跟普通包不太一样。第一种是常规做法:先装好devtools,然后用install_github拉取源码。第二种是离线安装:把源码包下载到本地后,用install.packages("路径/IOBR_xxx.tar.gz", repos=NULL, type="source")装。我一般优先用第一种,因为依赖关系能自动解析,离线方式主要留给内网服务器。
# 方式一:在线安装(推荐) install.packages("devtools") library(devtools) install_github("IOBR/IOBR", ref = "master", dependencies = TRUE) # 方式二:离线安装(适合内网或无外网环境) # install.packages("IOBR_xxx.tar.gz", repos = NULL, type = "source")这里最关键的参数是dependencies=TRUE。IOBR依赖大量Bioconductor包,比如limma、GSEABase、Biobase等,如果依赖缺失,安装过程会报错并提示你先去装BiocManager。常见的坑是在Windows上编译某些依赖包时找不到Rtools,导致安装中断。解决办法是提前装好与当前R版本匹配的Rtools,并在安装后把Rtools的路径写入系统环境变量。
安装完成后验证一下能否正常加载。加载IOBR需要一两分钟,因为它会把内置的signature基因集和参考矩阵载入内存,这个过程卡顿是正常的。
library(IOBR) # 成功后会输出包的版本号和主要函数列表 # 同时会附带提示:建议同时加载 tidyverse 和 survminer 用于后续数据处理和生存分析如果看到关于"没有找到某个依赖包"的报错,比如"there is no package called 'sva'",说明依赖解析没走完。此时不要急着重新安装整个IOBR,先单独装缺失的那个包,再加载一次,速度会快很多。依赖问题解决后,IOBR本身安装基本不会再出错。
2.2 表达矩阵的格式要求与预处理
IOBR的反卷积函数对表达矩阵的格式要求比较固定:必须是基因×样本的结构,基因名作为行名,样本名作为列名;如果方向反了,函数会在内部转置,但有些算法对矩阵类型很敏感,转置后行名可能变成整数,导致后续匹配失效。我踩过这个坑,所以后面养成习惯:任何表达矩阵进入IOBR之前,先强制检查dim()和rownames()。
表达矩阵里的基因名推荐使用官方Symbol,不要用Ensembl ID。原因在于IOBR内置的参考矩阵——比如CIBERSORT的LM22、MCPcounter的marker基因列表——全部基于Symbol。如果用Ensembl ID,函数会把基因名全部转成大写再做匹配,Ensembl ID一般不匹配反而全盘落空。如果你的原始数据是Ensembl ID,我一般会用biomaRt做一次映射,或者用注释包内的对应表转成Symbol再进来。
重复基因名是另一个高频问题。多个探针对应同一个Symbol时,直接去重会丢信息;常见做法是保留表达量最大的那一条,或者按均值合并。我一般用聚合方式:按Symbol分组后取均值,这样既保留基因的总体表达水平,又避免人为拔高。
# 表达矩阵预处理:假设 expr 是 样本×基因 的 data.frame expr <- as.data.frame(expr) expr$gene <- rownames(expr) # 先把基因名列挖出来 # 重复基因名按均值合并 library(dplyr) expr_agg <- expr %>% group_by(gene) %>% summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>% as.data.frame() rownames(expr_agg) <- expr_agg$gene expr_agg$gene <- NULL # 转置为 基因×样本 格式 expr_final <- as.matrix(expr_agg) expr_final <- t(expr_final) expr_final[1:3, 1:3] # 查看前几行几列关于缺失值:IOBR内部不会自动填补NA,如果表达矩阵里存在NA,部分算法会直接报错,另外一些算法会把NA当成0参与计算,导致细胞丰度被低估。最省事的方案是在预处理阶段强制把NA替换为整个矩阵的最小值,或者用该基因在所有样本中的中位数填充。我倾向用中位数,因为最小值的做法会把缺失项抹成零,放大低表达基因的贡献。
表达量是否取log,取决于你用的算法。CIBERSORT在芯片时代通常用不取log的原始表达强度,但放进IOBR后,内部会做标准化,因此输入log2或TPM问题不大。quanTIseq和EPIC对输入类型较敏感,一个期望得到TPM,另一个期望得到read counts或TPM。IOBR的官方示例里用的是log2(TPM+1),我实际用下来这个取值在多数算法下表现稳定,可以作为默认起点。
2.3 用小样本先跑通全流程
环境配好、矩阵整理好之后,不要急着把全部样本灌进去。我的习惯是先取三个样本组成一个小矩阵,跑一遍全流程,确认函数能正常返回结果,再放全量数据。这一步能省掉后面不少排查时间,因为反卷积算法本身很耗时,CIBERSORT一次置换检验就要几十秒,全量数据如果参数不对,跑完发现报错,返工成本极高。
# 取一小部分样本做冒烟测试 test_expr <- expr_final[, 1:3] test_result <- deconvo_tme( eset = test_expr, tumor_type = "pancan", methods = c("epic", "quantiseq", "mcpcounter"), perm = 100, nthreads = 2 ) # 查看返回结果的结构 names(test_result) # 每个算法对应一个数据框 test_result$epic[1:3, 1:5] # 查看 EPIC 输出的前几行这里tumor_type="pancan"表示泛癌模式,CIBERSORT和TIMER的部分参考数据会按泛癌默认参数处理。如果你研究的癌种在函数支持的列表里,比如"nsclc"或"melanoma",换成具体癌种能提高肿瘤纯度相关特征的匹配度。perm控制CIBERSORT的置换次数,冒烟测试时设100就够,正式分析建议1000。nthreads是并行线程数,在本地机器上不要超过物理核心数,否则IO wait反而拖慢速度。
冒烟测试的验收标准很简单:返回列表里的每个数据框行数等于样本数,列名是细胞类型或分数,数值范围大致合理——比例型算法如EPIC每行加起来接近1,MCPcounter则是相对分数,可能出现负值,这都正常。只有全部通过,才进入完整样本的正式运行。
3. IOBR核心操作:批量反卷积与特征基因集打分参数详解
3.1 六种反卷积算法的选型思路
deconvo_tme()支持的六种算法各有脾气,不是随便勾选就能完事。我把它们的输入偏好和输出类型整理成了一张对照表,方便对照选型。
| 算法 | 输入偏好 | 输出类型 | 适合场景 |
|---|---|---|---|
| CIBERSORT | 芯片或RNA-seq均可,需要置换检验 | 各类免疫细胞百分比 | 经典场景,结果直观,但耗时高 |
| EPIC | 推荐TPM,需要肿瘤纯度相关特征 | 细胞比例 | 关注CAF和内皮细胞时效果好 |
| quanTIseq | TPM或read counts | 细胞比例 | 批量RNA-seq数据,免疫细胞分型较细 |
| MCPcounter | 芯片或RNA-seq,标准化表达 | 相对丰度分数 | 做连续变量相关性分析时灵活 |
| TIMER | 基因表达 | 六类免疫细胞丰度 | 快速预览时用,细节偏粗 |
| xCell | 表达矩阵,依赖ssGSEA富集 | 富集分数 | 包含多种基质细胞和免疫细胞 |
选型时我的经验是:如果目的只是报告里画一张免疫细胞组成堆叠图,CIBERSORT加EPIC就够;如果要跟临床数据做相关性分析,加入MCPcounter和xCell能提供更多连续变量;如果样本来自单细胞结合bulk数据,quanTIseq对RNA-seq的兼容性最好。六种全开不是不行,但xCell和CIBERSORT的运行时间会明显拉长,数据量大时建议分两批跑。
3.2 deconvo_tme()的核心参数逐项说明
deconvo_tme()是我用这个包时最常调的函数,名字是deconvo_tme,全称是deconvolute tumor microenvironment。它把六种算法打包进一个统一入口,参数不多但每个都能影响结果。下面是一段正式的全量运行代码,注释里标了每个参数的作用。
full_result <- deconvo_tme( eset = expr_final, # 基因×样本表达矩阵,行名为Symbol tumor_type = "pancan", # 癌种类型,支持 "pancan" 或具体癌种 methods = c("cibersort", "epic", "quantiseq", "mcpcounter", "timer", "xcell"), arrays = FALSE, # TRUE表示芯片数据,FALSE表示RNA-seq perm = 1000, # CIBERSORT置换次数,正式分析用1000 nthreads = 6 # 并行线程数,按机器核心数调 )eset参数接受matrix或data.frame,但内部会先转成matrix再往下走,所以如果你传入的是tibble,部分算法可能因为列名类型问题报错,提前转成matrix更稳妥。tumor_type参数只影响CIBERSORT和TIMER,EPIC和MCPcounter不依赖癌种分类,quanTIseq在泛癌模式下也能正常出结果。
arrays参数需要认真确认。它决定数据是否来自芯片平台。arrays=FALSE时,IOBR会假定你是RNA-seq数据,并对表达矩阵做log2转换和标准化;arrays=TRUE时,内部会走芯片的标准化路径。如果设置错误,最常见的表现是MCPcounter返回全零值,因为参考marker的表达范围跟你的输入量纲完全不匹配。我常用的检查手段是看EPIC输出的总比例是否接近1,如果偏差超过0.2,大概率是arrays参数设反了。
3.3 特征基因集打分:从免疫细胞比例到免疫治疗相关指标
细胞比例只能回答"有哪些免疫细胞",很多人还想知道"这些细胞是否处于激活状态"或"患者对免疫治疗是否可能响应"。IOBR内置了一批特征基因集打分接口,核心思想是把一组已知功能状态的基因表达水平汇总成一个分数。常见做法是采用ssGSEA或GSVA富集打分,IOBR里对应的入口是calculate_sig_score(),支持tibble或matrix输入。
# 计算特征基因集分数,比如IFNG特征、T细胞炎症特征 signature_scores <- calculate_sig_score( eset = expr_final, signatures = IOBR::sig_tme, # 内置的TME特征基因集 method = "ssgsea", # 可选 "ssgsea" 或 "gsva" ssgsea_norm = TRUE # 是否对ssGSEA分数做归一化 )sig_tme是IOBR内置的免疫微环境特征基因集列表,包含IFNG响应、T细胞炎症、基质活化等常见通路。method参数推荐用"ssgsea",因为它在样本量小的时候比GSVA更稳定,且对表达量的绝对值不敏感。ssgsea_norm=TRUE会把分数归一化到0到1之间,便于跨样本比较;如果你后面要做生存分析的cutoff分组,归一化后的分数直接用更顺手。
特征基因集打分的意义在于,它能把细胞比例这种"数量"信息变成"活性"信息。比如CD8 T细胞比例高不代表细胞毒性强,但IFNG特征分数高往往跟细胞毒性相关。把这两类指标结合使用,是免疫浸润分析从描述性走向推断性的一步。
4. 结果整理与下游分析:相关性热图、分组比较和生存曲线
4.1 从列表到整齐数据框
deconvo_tme()返回的结果是一个列表,每个元素对应一种算法。这个结构方便查看,但不方便后续分析。我拿到结果后第一件事是把不同算法的细胞类型列合并成一个大宽表,行是样本,列是"算法_细胞类型"。
library(tidyverse) # 把列表中每个数据框转成样本×细胞类型的格式 clean_results <- function(res_list) { purrr::imap(res_list, function(df, nm) { df %>% rownames_to_column("sample") %>% pivot_longer(-sample, names_to = "cell_type", values_to = "score") %>% mutate(algorithm = nm) }) %>% bind_rows() %>% pivot_wider(id_cols = c(sample, algorithm), names_from = cell_type, values_from = score) } tme_wide <- clean_results(full_result)这里用pivot_longer和pivot_wider做双向变换,目的是把原本每个算法一个宽表的数据,先拉长再合并,最后转成整洁表格。注意算法名作为一列保留,这样后续做分组对比时可以用algorithm列做facet。我在实际项目中会把tme_wide直接写出CSV存档,方便跟临床表型数据用sample列匹配。
4.2 免疫细胞相关性矩阵的计算与可视化
合表之后最常见的分析是看不同免疫细胞之间的相关性。相关性矩阵能告诉你哪些细胞常常共浸润,哪些此消彼长。计算时建议用Spearman相关,因为细胞丰度分数通常不服从正态分布,尤其MCPcounter和xCell的输出偏态明显。
# 选一种算法的结果做相关性矩阵 epic_mat <- clean_results(full_result) %>% filter(algorithm == "epic") %>% column_to_rownames("sample") %>% select(-sample) cor_mat <- cor(epic_mat, method = "spearman", use = "complete.obs") # 可视化 library(corrplot) corrplot( cor_mat, method = "color", type = "upper", order = "hclust", tl.col = "black", tl.cex = 0.7, number.cex = 0.6, addCoef.col = "grey30" )corrplot的order="hclust"会对变量做层次聚类排序,让相关的细胞类型聚在一起,图面可读性比按字母排序好很多。tl.cex控制基因名大小,细胞类型名字较长时建议调到0.7以下,否则会互相重叠。addCoef.col把相关系数数字叠加在格子上,方便直接读数值。
这类图在文章里通常作为肿瘤微环境特征展示,配合下面要写的生存分析一起出现,可以支撑"特定免疫细胞浸润与预后相关"的结论。
4.3 与临床数据整合:生存分析标准流程
免疫浸润分析的项目最终几乎都要落到生存分析:某类细胞丰度高是不是意味着更好的总生存期或无进展生存期。做法是把前面算出的细胞比例或特征分数合并到临床表型表上,然后按中位数或最优cutoff分组画KM曲线。
# 假设 clin 有 sample、OS_time、OS_status 三列 clin_data <- clin %>% left_join(tme_wide %>% select(sample, starts_with("CD8")), by = "sample") # 按中位数分组 clin_data$cd8_group <- ifelse( clin_data$CD8_T_cells >= median(clin_data$CD8_T_cells, na.rm = TRUE), "High", "Low" ) library(survival) library(survminer) fit <- survfit(Surv(OS_time, OS_status) ~ cd8_group, data = clin_data) ggsurvplot( fit, data = clin_data, pval = TRUE, risk.table = TRUE, conf.int = FALSE, palette = c("#E64B35", "#4DBBD5"), legend.title = "CD8 T cells", xlab = "Time (months)" )这里注意合并列名要跟tme_wide里的实际列名一致。列名如果带特殊字符,比如CIBERSORT输出的Naive_CD4_T_cells里有空格,left_join时要用反引号引用。median分组是保守做法,如果样本量超过80,可以试试surv_cutpoint找最优cutoff,但最优cutoff容易过拟合,文章里用它需要额外做内部验证。我在项目里一般先按中位分组出一版结果,如果p值在0.05附近,再用surv_cutpoint确认稳定性。
生存分析的意义是给免疫浸润结果赋予临床价值。某类细胞比例高低本身不能说明问题,但一旦跟生存期挂钩,就变成了可以写进结论的发现。需要注意的是样本量太小时KM曲线不可靠,少于20个事件(死亡或复发)时,我通常不会把生存结果当作主要结论。
5. 避坑与常见问题排查:五个高频报错和一个性能陷阱
5.1 报错"Error in t(eset) : argument is not a matrix"
现象:运行deconvo_tme()时直接报错,提示t()函数参数不是矩阵。原因通常是表达矩阵传入的是data.frame或tibble,且内部某些算法会先做转置操作,data.frame转置后变成matrix,行名反而丢了。解决:进函数之前强制转成matrix,并确认行名存在。
expr_final <- as.matrix(expr_final) storage.mode(expr_final) <- "numeric" class(expr_final) # 必须显示 "matrix" "array"5.2 跑xCell时R会话直接崩溃或内存溢出
现象:六种算法一起跑,前面几种顺利出结果,到xCell时R进程突然被杀掉。原因:xCell内部对数千个基因集做ssGSEA富集,内存占用跟基因数和样本数乘积成正比,全量跑时直接撑爆32G内存。解决:把xCell单独拎出来跑,并且限制样本数,分批次运行后再合并结果。
# 分批次跑 xCell batch1 <- deconvo_tme(eset = expr_final[, 1:20], methods = "xcell", perm = 100, nthreads = 4) batch2 <- deconvo_tme(eset = expr_final[, 21:40], methods = "xcell", perm = 100, nthreads = 4) xcell_full <- cbind(batch1$xcell, batch2$xcell)5.3 安装时报错"ERROR: dependencies 'sva', 'limma' are not available"
现象:install_github时依赖包解析失败,提示多个Bioconductor包缺失。原因:IOBR从Github安装时不会自动配置Bioconductor源,你的R环境里没有设置BiocManager。解决:先手动安装BiocManager并配置Bioconductor源,再重新安装IOBR。
install.packages("BiocManager") BiocManager::install(c("sva", "limma", "GSEABase", "biomaRt")) library(devtools) install_github("IOBR/IOBR", dependencies = TRUE)5.4 CIBERSORT结果出现全零行或负值
现象:cibersort对应的数据框里某些样本全部是0,或者部分细胞类型出现负数。原因:CIBERSORT的置换检验在样本表达谱跟参考矩阵匹配度过低时,会输出全零;负值则来自内部反卷积的约束优化,在肿瘤纯度很高时容易出现。解决:不要直接删除样本,先检查肿瘤类型设置。把tumor_type从"pancan"改成更具体的癌种,匹配度会显著提升。如果已经是具体癌种,考虑用EPIC的结果交叉验证全零样本。
5.5 表达矩阵行名不是gene symbol导致匹配失败
现象:算法运行没有报错,但结果里所有细胞类型的数值都接近常量,不同样本之间几乎无差异。原因:行名是Ensembl ID或探针ID,参考矩阵匹配不到基因,算法退化为在相同基因子集上反复计算。解决:在进入IOBR之前把行名映射成Symbol,映射后用unique()检查是否有重复,重复的按2.2节的方式合并。
5.6 性能陷阱:perm设置过高导致CIBERSORT跑一天
现象:CIBERSORT单独跑一个100样本的队列,perm=1000,跑了十几个小时还没结束。原因:CIBERSORT的置换检验复杂度是O(perm×sample),perm从500提到1000,耗时线性翻倍。解决:先确认CIBERSORT的perm对最终结果产生的影响范围。我对比过同一批数据perm=100和perm=1000的结果,细胞比例的差异通常小于0.01,所以正式分析用perm=500已经足够,先跑一版perm=100的预览结果,再决定是否加高。
这六条踩坑记录的共性是:九成问题出在数据进函数之前,而不是算法本身。IOBR把六种算法封在一个接口下,降低了使用门槛,但也把数据兼容性问题集中到了一起。参数调整前先检查输入格式,能省掉大量返工时间。
6. 一个收尾技巧:为每个样本生成TME报告并顺手做质量检查
反卷积跑完之后,真正让结果可信的是交叉验证和可追溯性。我现在的习惯是写一个批处理脚本,为每个样本输出一页免疫细胞组成堆叠图和一张细胞比例表,同时在表格末尾附上质量检查列,帮助快速定位异常样本。
# 为每个样本生成免疫细胞组成堆叠图 library(ggplot2) tme_long <- clean_results(full_result) %>% filter(algorithm == "epic") %>% pivot_longer(-c(sample, algorithm), names_to = "cell_type", values_to = "score") sample_list <- unique(tme_long$sample) for (s in sample_list) { p <- tme_long %>% filter(sample == s) %>% mutate(cell_type = fct_reorder(cell_type, score)) %>% ggplot(aes(x = sample, y = score, fill = cell_type)) + geom_bar(stat = "identity", width = 0.6) + coord_flip() + theme_minimal(base_size = 11) + labs(title = paste0("TME composition: ", s), x = NULL, y = "Estimated proportion") ggsave(paste0("TME_report_", s, ".pdf"), p, width = 8, height = 4.5) }这个循环从每个样本的EPIC结果里生成堆叠条形图,文件名带样本ID,方便跟临床数据一一对应。fct_reorder让细胞类型按分数排序,图面更有层次。每张图都是一个检查点:如果某张图的堆叠柱异常短,说明该样本的免疫细胞比例整体偏低,需要回头确认表达矩阵是否存在批次效应。
质量检查方面,我通常在合表之后加一个判断:比例型算法每一行的总和应该落在0.95到1.05之间,超出这个区间就标记为"check"。这个检查不针对绝对精度,而是快速识别那些输入表达量异常导致总比例偏移的样本。对于MCPcounter和xCell这种非比例型分数,我则检查是否有超过3倍四分位距的极端值,有的话单独画箱线图确认。
从那以后,我每次跑反卷积都会强制走一遍这套流程:先小样本冒烟测试,再检查行名和矩阵类型,然后全量运行,最后逐样本看图。这套流程看起来繁琐,但它帮我挡掉了好几次因为表达矩阵格式问题而产出的无效结果。如果你正准备在自己的转录组数据上做免疫浸润分析,这份IOBR免疫浸润分析R包源码值得下载一份,按上面的顺序跑一遍,能少走不少弯路,希望帮到你。
本文还有配套的精品资源,点击获取