159_肿瘤微环境评估¶

# 安装
# install.packages("estimate")
# 或从Bioconductor/GitHub
library(estimate)

# 1. 输入: 基因表达矩阵 (基因×样本, 行名为基因symbol)
# ESTIMATE需要特定格式的GCT文件
write.table(expression_matrix, "expression.txt", sep="\t", quote=FALSE)

# 使用ESTIMATE函数
# 过滤共有基因
filterCommonGenes(
  input.f = "expression.txt",
  output.f = "expression_filtered.gct",
  id = "GeneSymbol"
)

# 计算三个评分
estimateScore(
  input.ds = "expression_filtered.gct",
  output.ds = "estimate_scores.gct",
  platform = "illumina"  # "affymetrix" 或 "illumina"
)

# 读取结果
scores <- read.delim("estimate_scores.gct", skip = 2)
# 包含: StromalScore, ImmuneScore, ESTIMATEScore, TumorPurity

# 可视化
library(ggplot2)
scores_long <- tidyr::pivot_longer(scores, cols = c("StromalScore", "ImmuneScore"),
                                    names_to = "Score_Type", values_to = "Score")
ggplot(scores_long, aes(x = Score_Type, y = Score)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "ESTIMATE Scores")

# 按肿瘤纯度分层
scores$PurityGroup <- ifelse(scores$TumorPurity > median(scores$TumorPurity),
                              "High_Purity", "Low_Purity")

# TIMER2.0在线平台: http://timer.cistrome.org/
# 或使用immunedeconv R包调用

# 安装immunedeconv (统一接口)
# devtools::install_github("omnideconv/immunedeconv")
library(immunedeconv)

# 输入: TPM表达矩阵 (基因symbol × 样本)
expr_tpm <- as.matrix(expression_tpm)

# === TIMER方法 ===
timer_result <- deconvolute(expr_tpm, method = "timer",
                             indications = rep("blca", ncol(expr_tpm)))  # TCGA肿瘤类型缩写

# === quanTIseq ===
quantiseq_result <- deconvolute(expr_tpm, method = "quantiseq",
                                 tumor = TRUE)

# === CIBERSORT (需要注册token) ===
# cibersort_result <- deconvolute_cibersort(expr_tpm, 
#                                            arrays = FALSE)

# === EPIC ===
epic_result <- deconvolute(expr_tpm, method = "epic")

# === 综合多种方法 ===
methods <- c("timer", "quantiseq", "mcp_counter", "epic", "xcell")
all_results <- list()
for (m in methods) {
  tryCatch({
    if (m == "timer") {
      all_results[[m]] <- deconvolute(expr_tpm, method = m,
                                       indications = rep("blca", ncol(expr_tpm)))
    } else {
      all_results[[m]] <- deconvolute(expr_tpm, method = m)
    }
  }, error = function(e) cat(m, "failed:", e$message, "\n"))
}

# 合并T细胞估计并比较
library(ggplot2)
tcell_comparison <- data.frame(
  Sample = colnames(expr_tpm),
  EPIC_Tcell = as.numeric(all_results$epic["T cell CD4+", ]),
  quanTIseq_Tcell = as.numeric(all_results$quantiseq["T cell CD4+", ]),
  MCP_Tcell = as.numeric(all_results$mcp_counter["T cell", ])
)

# 安装
# devtools::install_github("ebecht/MCPcounter", ref="master", subdir="Source")
library(MCPcounter)

# MCPcounter估计10种细胞类型的丰度
# 输入: 基因表达矩阵(基因×样本)
mcp_result <- MCPcounter.estimate(
  expression = expr_tpm,
  featuresType = "HUGO_symbols"  # 或 "ENSEMBL_ID"
)

# 输出: 10种细胞类型 × N样本的评分矩阵
# T cells, CD8 T cells, Cytotoxic lymphocytes, B lineage
# NK cells, Monocytic lineage, Myeloid dendritic cells
# Neutrophils, Endothelial cells, Fibroblasts
print(rownames(mcp_result))

# 热图可视化
library(pheatmap)
pheatmap(mcp_result,
         scale = "row",
         clustering_method = "ward.D2",
         annotation_col = sample_annotation,
         color = colorRampPalette(c("blue", "white", "red"))(100),
         main = "MCPcounter Cell Type Scores")

# 与临床变量关联分析
library(survival)
library(survminer)

# 以CD8 T细胞为例
clinical$CD8_score <- mcp_result["CD8 T cells", rownames(clinical)]
clinical$CD8_group <- ifelse(clinical$CD8_score > median(clinical$CD8_score),
                              "High", "Low")

fit <- survfit(Surv(OS_time, OS_event) ~ CD8_group, data = clinical)
ggsurvplot(fit, data = clinical, pval = TRUE, risk.table = TRUE,
           title = "OS by CD8 T cell Infiltration")

# 安装: devtools::install_github("GfellerLab/EPIC")
library(EPIC)

# EPIC估计细胞比例(和为1)，包含"otherCells"(未识别成分)
epic_result <- EPIC(
  bulk = expr_tpm,           # TPM矩阵
  reference = NULL,          # 使用默认参考(TRef/BRef)
  mRNA_cell = NULL,          # 默认mRNA/cell估计
  sigGenes = NULL            # 默认签名基因
)

# 结果包含:
# $cellFractions: 细胞比例矩阵
# $fit.gof: 拟合优度
cell_fractions <- epic_result$cellFractions
print(head(cell_fractions))
# 细胞类型: B cells, CAFs, CD4+ T cells, CD8+ T cells, 
#           Endothelial, Macrophages, NK cells, otherCells

# 堆叠柱状图
library(tidyr)
library(ggplot2)
fractions_df <- as.data.frame(cell_fractions)
fractions_df$Sample <- rownames(fractions_df)
fractions_long <- pivot_longer(fractions_df, -Sample, 
                               names_to = "CellType", values_to = "Fraction")

ggplot(fractions_long, aes(x = Sample, y = Fraction, fill = CellType)) +
  geom_bar(stat = "identity") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) +
  labs(title = "EPIC Cell Type Fractions") +
  scale_fill_brewer(palette = "Set3")

# === 整合多种方法的免疫评分 ===
# 标准化各方法的评分到0-1范围
scale_01 <- function(x) (x - min(x, na.rm=T)) / (max(x, na.rm=T) - min(x, na.rm=T))

immune_scores <- data.frame(
  ESTIMATE_immune = scale_01(estimate_scores$ImmuneScore),
  MCP_CD8 = scale_01(mcp_result["CD8 T cells", ]),
  EPIC_CD8 = scale_01(epic_result$cellFractions[, "CD8+ T cells"]),
  row.names = colnames(expr_tpm)
)

# 相关性检验
cor_matrix <- cor(immune_scores, use = "complete.obs", method = "spearman")
print(round(cor_matrix, 3))

# === 免疫亚型分类 ===
# 基于免疫浸润模式将肿瘤分为免疫亚型
library(ConsensusClusterPlus)

# 使用多种免疫评分做聚类
immune_matrix <- t(scale(t(rbind(mcp_result, estimate_immune = estimate_scores$ImmuneScore))))

cc_results <- ConsensusClusterPlus(
  d = immune_matrix,
  maxK = 6,
  reps = 100,
  pItem = 0.8,
  clusterAlg = "hc",
  distance = "pearson",
  seed = 42
)

# 选择最优K
immune_subtypes <- cc_results[[3]]$consensusClass  # 假设K=3最优

# === 与免疫治疗响应关联 ===
# 比较不同免疫亚型的PD-L1表达和TMB
boxplot_data <- data.frame(
  Subtype = factor(immune_subtypes),
  PDL1 = expr_tpm["CD274", ],  # PD-L1 = CD274
  TMB = sample_tmb
)

ggplot(boxplot_data, aes(x = Subtype, y = PDL1, fill = Subtype)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "PD-L1 Expression by Immune Subtype")

# === R包安装 ===
Rscript -e '
install.packages("estimate")
devtools::install_github("omnideconv/immunedeconv")
devtools::install_github("ebecht/MCPcounter", ref="master", subdir="Source")
devtools::install_github("GfellerLab/EPIC")
BiocManager::install("GSVA")
'

# === TIMER2.0在线分析 ===
# 访问 http://timer.cistrome.org/
# 上传TPM表达矩阵
# 选择肿瘤类型和分析模块

# === 批量分析脚本 ===
Rscript tme_analysis.R --expr expression_tpm.csv \
  --clinical clinical.csv \
  --methods "ESTIMATE,MCPcounter,EPIC,quanTIseq" \
  --output results/