保姆级教程：用R语言estimate包给TCGA数据算免疫评分和肿瘤纯度（附完整代码）

零基础实战：用R语言ESTIMATE包解析TCGA肿瘤微环境

肿瘤微环境分析已成为癌症研究的核心方向之一。想象一下，你手头有一批TCGA的转录组数据，如何从中挖掘出肿瘤纯度、免疫细胞浸润程度这些关键指标？来自MD Anderson的ESTIMATE算法为我们提供了一种高效解决方案。不同于需要复杂湿实验的病理分析，这个方法仅需基因表达数据就能推算出基质细胞、免疫细胞比例及肿瘤纯度，特别适合没有湿实验条件的生物信息学研究者。

1. 环境准备与数据获取

1.1 安装必要的R包

首先确保你的R版本≥3.6.0。打开RStudio，依次执行以下命令安装核心依赖：

# 设置CRAN镜像加速安装 options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) install.packages(c("utils", "BiocManager")) BiocManager::install("estimate")

常见问题处理：

• 若遇到dependency 'xxx' is not available错误，尝试先单独安装缺失依赖
• 网络超时可添加timeout = 600参数延长下载时限

1.2 获取TCGA表达矩阵

推荐从UCSC Xena浏览器下载已标准化的FPKM数据：

1. 访问 https://xenabrowser.net/datapages/
2. 选择TCGA TARGET GTEx→Gene Expression→FPKM
3. 筛选目标癌种后下载TSV文件

典型数据结构示例：

注意：确保第一列为基因ID，后续列为样本ID，不要包含额外注释行

2. 数据预处理关键步骤

2.1 表达量log2转换

原始FPKM值通常呈现右偏分布，需进行对数转换：

# 读取下载的表达矩阵 expr_data <- read.delim("TCGA_BRCA_FPKM.tsv", row.names=1) # 过滤低表达基因（至少10%样本表达量>1） keep <- rowSums(expr_data > 1) >= ncol(expr_data)*0.1 expr_filtered <- expr_data[keep, ] # log2转换并处理零值 expr_log2 <- log2(expr_filtered + 1)

转换前后分布对比：

2.2 基因ID转换

ESTIMATE要求使用Hugo Symbol标识基因，常用biomaRt进行转换：

library(biomaRt) ensembl <- useMart("ensembl", dataset="hsapiens_gene_ensembl") gene_symbols <- getBM( attributes=c('ensembl_gene_id', 'hgnc_symbol'), filters='ensembl_gene_id', values=rownames(expr_log2), mart=ensembl)

匹配失败的可能原因：

• 旧版ENSEMBL ID已弃用
• 非编码基因缺乏标准命名
• 物种信息不匹配

3. ESTIMATE分析全流程

3.1 生成GCT格式文件

使用estimate包的专用函数进行格式转换：

library(estimate) write.table(expr_log2, "expr_log2.txt", sep="\t", quote=F) filterCommonGenes( input.f="expr_log2.txt", output.f="expr_estimate.gct", id="GeneSymbol")

成功转换后的GCT文件头两行示例：

#1.2 10182 50

3.2 计算三大评分

根据平台类型选择参数（RNA-seq选择illumina）：

estimateScore( "expr_estimate.gct", "estimate_scores.txt", platform="illumina")

输出文件包含四列关键指标：

1. StromalScore - 基质细胞浸润程度
2. ImmuneScore - 免疫细胞浸润水平
3. ESTIMATEScore - 综合评分
4. TumorPurity - 肿瘤纯度估计值

3.3 结果可视化

建议用ggplot2绘制评分分布：

library(ggplot2) scores <- read.delim("estimate_scores.txt", row.names=1) ggplot(scores, aes(x=TumorPurity)) + geom_histogram(bins=30, fill="#69b3a2") + labs(title="Tumor Purity Distribution", x="Purity", y="Count")

4. 实战问题排查指南

4.1 常见报错解决方案

4.2 结果验证技巧

1. 交叉验证：用相同数据运行CIBERSORT比较免疫评分趋势
2. 生物学合理性检查：高纯度样本应显示低基质/免疫评分
3. 极端值排查：检查ESTIMATEScore超出[-2000,2000]范围的样本

4.3 性能优化建议

对于大型数据集（>500样本）：

• 使用doParallel并行处理
• 预先过滤低表达基因减少计算量
• 考虑分批次运行后合并结果

library(doParallel) registerDoParallel(cores=4) # 后续estimateScore会自动利用多核

5. 高级应用场景

5.1 多癌种联合分析

合并不同TCGA项目数据时需注意：

• 批次效应矫正（推荐ComBat-seq）
• 平台差异处理（统一转为TPM/FPKM）
• 样本均衡性考虑（各癌种样本量匹配）

5.2 与临床数据关联

将ESTIMATE结果与临床信息合并分析：

clinical <- read.csv("TCGA_clinical.csv") merged_data <- merge(scores, clinical, by.x="row.names", by.y="patient_id") # 生存分析示例 library(survival) coxph(Surv(OS_time, OS_status) ~ ImmuneScore, data=merged_data)

5.3 自定义特征基因集

进阶用户可修改源码中的基因列表：

1. 下载estimate包源码
2. 修改data/stromal_signature.rda和immune_signature.rda
3. 重新编译安装

重要提示：修改特征基因需严格验证，建议保留原始版本备份

6. 结果解读与报告生成

6.1 评分生物学意义

6.2 自动化报告模板

推荐使用R Markdown生成可交互报告：

--- title: "ESTIMATE Analysis Report" output: html_document --- ```{r setup} library(estimate) scores <- read.delim("estimate_scores.txt")

Key Findings

Median tumor purity:r median(scores$TumorPurity)

library(plotly) plot_ly(scores, x=~StromalScore, y=~ImmuneScore, color=~TumorPurity)

### 6.3 数据存储建议 建立标准化存储结构：

/project /raw_data # 原始TCGA文件 /processed # 转换后的矩阵 /results # ESTIMATE输出 /scripts # 分析代码 /reports # 最终报告

实际项目中，我们通常会遇到不同癌种间的评分分布差异。比如在乳腺癌数据中经常观察到基质评分与ER状态的相关性，而肺癌数据中免疫评分往往与PD-L1表达呈正相关。这些模式验证了ESTIMATE结果的生物学合理性，也提示我们可以进一步探索特定癌种的微环境特征。