免疫浸润分析（CIBERSORT/xCell/MCP-counter/TIMER）¶

# ===== CIBERSORT/CIBERSORTx R实现 =====
# 方法1：使用IOBR包（整合多种方法）
# remotes::install_github("IOBR/IOBR")  # 从GitHub安装，不在Bioconductor上
library(IOBR)

# 准备表达矩阵（基因symbol为行名，TPM/FPKM值）
expr_matrix <- as.matrix(read.csv("expression_tpm.csv", row.names = 1))

# CIBERSORT分析
cibersort_result <- deconvo_tme(
  eset = expr_matrix,
  method = "cibersort",
  arrays = FALSE,           # FALSE=RNA-seq, TRUE=microarray
  perm = 1000              # 置换次数
)
head(cibersort_result)

# 方法2：使用immunedeconv包
# 注意：immunedeconv 不在 CRAN 上，需从 GitHub 安装
# remotes::install_github("omnideconv/immunedeconv")
library(immunedeconv)

# CIBERSORT（需要先在 https://cibersort.stanford.edu 注册获取学术许可，
#   下载 CIBERSORT.R 脚本和 LM22.txt 参考矩阵，然后设置路径：）
# set_cibersort_binary("/path/to/CIBERSORT.R")
# set_cibersort_mat("/path/to/LM22.txt")
# cibersort_res <- deconvolute(expr_matrix, "cibersort")

# CIBERSORT absolute mode
# cibersort_abs <- deconvolute(expr_matrix, "cibersort_abs")

# 方法3：手动实现CIBERSORT（使用e1071包的SVR）
library(e1071)
library(preprocessCore)

# 读取LM22参考矩阵
lm22 <- read.table("LM22.txt", header = TRUE, sep = "\t", row.names = 1)

# 找共有基因
common_genes <- intersect(rownames(expr_matrix), rownames(lm22))
expr_sub <- expr_matrix[common_genes, ]
lm22_sub <- lm22[common_genes, ]

# 分位数归一化
expr_norm <- normalize.quantiles(as.matrix(expr_sub))
rownames(expr_norm) <- common_genes
colnames(expr_norm) <- colnames(expr_sub)

# 对每个样本做SVR反卷积
cibersort_proportions <- matrix(NA, ncol(expr_norm), ncol(lm22_sub))
colnames(cibersort_proportions) <- colnames(lm22_sub)
rownames(cibersort_proportions) <- colnames(expr_norm)

for (i in 1:ncol(expr_norm)) {
  # ν-SVR
  svr_model <- svm(x = as.matrix(lm22_sub), y = expr_norm[, i],
                   type = "nu-regression", kernel = "linear", nu = 0.5)
  coefs <- t(svr_model$coefs) %*% svr_model$SV
  coefs[coefs < 0] <- 0  # 非负约束
  coefs <- coefs / sum(coefs)  # 归一化到和为1
  cibersort_proportions[i, ] <- coefs
}

# ===== xCell分析 =====
# 注意：xCell 不在 CRAN 上，需从 GitHub 安装
# devtools::install_github("dviraran/xCell")
# 也可使用升级版 xCell2（已正式收录于Bioconductor）：
# BiocManager::install("xCell2")
# 或从GitHub安装开发版：devtools::install_github("AlmogAngel/xCell2")
library(xCell)

# 运行xCell（输入为TPM或归一化后的微阵列数据）
xcell_result <- xCellAnalysis(expr_matrix)
dim(xcell_result)  # 67种细胞/特征 × 样本数

# 查看结果
head(xcell_result[, 1:5])
# 结果包含64种细胞类型 + ImmuneScore + StromaScore + MicroenvironmentScore

# 提取特定细胞类型
immune_cells <- c("B-cells", "CD4+ T-cells", "CD8+ T-cells",
                  "NK cells", "Macrophages M1", "Macrophages M2",
                  "Tregs", "Dendritic cells")
xcell_immune <- xcell_result[immune_cells, ]

# 使用immunedeconv统一接口
xcell_res <- deconvolute(expr_matrix, "xcell")

# 可视化
library(pheatmap)
annotation_col <- data.frame(
  Group = sample_info$Group,
  row.names = colnames(xcell_immune)
)
pheatmap(xcell_immune,
         annotation_col = annotation_col,
         scale = "row",
         clustering_distance_rows = "correlation",
         main = "xCell Immune Infiltration")

# ===== MCP-counter =====
# 注意：MCPcounter 不在 CRAN 上，需从 GitHub 安装
# devtools::install_github("ebecht/MCPcounter", ref = "master", subdir = "Source")
library(MCPcounter)

# 运行MCP-counter
mcp_result <- MCPcounter.estimate(
  expression = expr_matrix,
  featuresType = "HUGO_symbols",  # 基因名类型
  # 需联网下载；或提前下载到本地后改用本地路径
  genes = read.table(curl("http://raw.githubusercontent.com/ebecht/MCPcounter/master/Signatures/genes.txt"),
                     sep = "\t", stringsAsFactors = FALSE, header = TRUE,
                     colClasses = "character", check.names = FALSE)
)

dim(mcp_result)  # 10种细胞类型 × 样本数
rownames(mcp_result)
# [1] "T cells"          "CD8 T cells"      "Cytotoxic lymphocytes"
# [4] "B lineage"        "NK cells"         "Monocytic lineage"
# [7] "Myeloid dendritic cells" "Neutrophils"  "Endothelial cells"
# [10] "Fibroblasts"

# 或使用immunedeconv
mcp_res <- deconvolute(expr_matrix, "mcp_counter")

# 可视化
library(ggplot2)
library(tidyr)

mcp_long <- as.data.frame(t(mcp_result))
mcp_long$Sample <- rownames(mcp_long)
mcp_long$Group <- sample_info$Group
mcp_long <- pivot_longer(mcp_long, cols = -c(Sample, Group),
                          names_to = "CellType", values_to = "Score")

ggplot(mcp_long, aes(x = CellType, y = Score, fill = Group)) +
  geom_boxplot() +
  coord_flip() +
  theme_classic() +
  labs(title = "MCP-counter Immune Infiltration", y = "MCP-counter Score")

# ===== TIMER (R实现) =====
# 使用immunedeconv包
timer_res <- deconvolute(expr_matrix, "timer",
                          indications = rep("brca", ncol(expr_matrix)))
# TIMER需要指定肿瘤类型(indication)

# ===== EPIC =====
# install.packages("remotes")
# remotes::install_github("GfellerLab/EPIC")  # 从GitHub安装，不在CRAN上
library(EPIC)
epic_res <- EPIC(bulk = expr_matrix)
epic_proportions <- epic_res$cellFractions

# ===== quanTIseq =====
quantiseq_res <- deconvolute(expr_matrix, "quantiseq",
                              tumor = TRUE,  # 是否含肿瘤细胞
                              arrays = FALSE)

# ===== 多方法整合比较（immunedeconv一站式） =====
library(immunedeconv)

# 运行所有方法
methods <- c("xcell", "mcp_counter", "epic", "quantiseq")
all_results <- list()
for (method in methods) {
  all_results[[method]] <- deconvolute(expr_matrix, method)
}

# 比较不同方法对同一细胞类型的估计
# 例如CD8 T cells
cd8_comparison <- data.frame(
  xCell = as.numeric(all_results$xcell["CD8+ T-cells", -1]),
  MCP = as.numeric(all_results$mcp_counter["CD8+ T-cells", -1]),
  EPIC = as.numeric(all_results$epic["CD8+ T-cells", -1]),
  quanTIseq = as.numeric(all_results$quantiseq["CD8+ T-cells", -1])
)
# 计算方法间相关性
cor(cd8_comparison, use = "complete.obs", method = "spearman")

# ===== ESTIMATE =====
library(estimate)

# 准备输入文件
write.table(data.frame(Gene = rownames(expr_matrix), expr_matrix),
            "expr_for_estimate.txt", sep = "\t", quote = FALSE, row.names = FALSE)

# 过滤共同基因
filterCommonGenes(input.f = "expr_for_estimate.txt",
                  output.f = "expr_filtered.gct", id = "GeneSymbol")

# 计算ESTIMATE分数
estimateScore(input.ds = "expr_filtered.gct",
              output.ds = "estimate_scores.gct",
              platform = "illumina")

# 读取结果
estimate_scores <- read.table("estimate_scores.gct", skip = 2, header = TRUE, sep = "\t")
rownames(estimate_scores) <- estimate_scores[, 1]
estimate_scores <- estimate_scores[, -c(1, 2)]

# 提取三个分数
stromal_score <- as.numeric(estimate_scores["StromalScore", ])
immune_score <- as.numeric(estimate_scores["ImmuneScore", ])
estimate_total <- as.numeric(estimate_scores["ESTIMATEScore", ])

# 计算肿瘤纯度
tumor_purity <- cos(0.6049872018 + 0.0001467884 * estimate_total)

# 可视化
plot_data <- data.frame(
  ImmuneScore = immune_score,
  StromalScore = stromal_score,
  TumorPurity = tumor_purity,
  Group = sample_info$Group
)

ggplot(plot_data, aes(x = Group, y = ImmuneScore, fill = Group)) +
  geom_boxplot() +
  stat_compare_means() +
  theme_classic() +
  labs(title = "ESTIMATE Immune Score")

# ===== 生存分析 =====
library(survival)
library(survminer)

# 对每种细胞类型做KM生存分析
cell_types <- rownames(cibersort_proportions_t)  # 转置后行为细胞类型
surv_results <- list()

for (ct in cell_types) {
  scores <- cibersort_proportions[, ct]
  group <- ifelse(scores > median(scores), "High", "Low")

  surv_data <- data.frame(
    time = clinical$OS_time,
    status = clinical$OS_status,
    group = factor(group, levels = c("Low", "High"))
  )

  fit <- survfit(Surv(time, status) ~ group, data = surv_data)
  cox <- coxph(Surv(time, status) ~ group, data = surv_data)
  p <- summary(cox)$coefficients[, "Pr(>|z|)"]
  hr <- exp(coef(cox))

  surv_results[[ct]] <- list(fit = fit, p = p, hr = hr)
}

# 筛选显著的细胞类型
sig_cells <- names(which(sapply(surv_results, function(x) x$p) < 0.05))

# 绘制显著的KM曲线
for (ct in sig_cells) {
  print(ggsurvplot(surv_results[[ct]]$fit,
                   data = data.frame(time = clinical$OS_time,
                                     status = clinical$OS_status),
                   pval = TRUE, risk.table = TRUE,
                   title = ct))
}

# ===== 免疫浸润与免疫治疗反应 =====
# 比较responder vs non-responder的免疫浸润差异
responders <- sample_info$Response == "Responder"
for (ct in colnames(cibersort_proportions)) {
  p <- wilcox.test(cibersort_proportions[responders, ct],
                   cibersort_proportions[!responders, ct])$p.value
  if (p < 0.05) {
    cat(ct, ": P =", p, "\n")
  }
}

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

# 使用多种免疫细胞评分做聚类
immune_matrix <- t(cibersort_proportions)  # 行=细胞类型, 列=样本
results_cc <- ConsensusClusterPlus(
  d = immune_matrix,
  maxK = 6,
  reps = 1000,
  pItem = 0.8,
  pFeature = 1,
  clusterAlg = "km",
  distance = "euclidean",
  seed = 42
)

# 选择最佳K值（通常看CDF和delta area）
# 假设K=3最佳
immune_subtypes <- results_cc[[3]]$consensusClass
sample_info$ImmuneSubtype <- paste0("IS", immune_subtypes)

# 比较免疫亚型的生存差异
fit_subtype <- survfit(Surv(OS_time, OS_status) ~ ImmuneSubtype, data = sample_info)
ggsurvplot(fit_subtype, pval = TRUE, risk.table = TRUE)

# ===== 综合可视化 =====
library(ggplot2)
library(tidyr)
library(ComplexHeatmap)

# 1. 堆积柱状图（CIBERSORT比例）
cib_long <- as.data.frame(cibersort_proportions)
cib_long$Sample <- rownames(cib_long)
cib_long <- pivot_longer(cib_long, -Sample, names_to = "CellType", values_to = "Proportion")

# 按某细胞排序
sample_order <- rownames(cibersort_proportions)[order(cibersort_proportions[, "T cells CD8"])]
cib_long$Sample <- factor(cib_long$Sample, levels = sample_order)

ggplot(cib_long, aes(x = Sample, y = Proportion, fill = CellType)) +
  geom_bar(stat = "identity") +
  theme_classic() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  labs(title = "CIBERSORT Immune Cell Composition")

# 2. 相关性热图（细胞类型间相关性）
cell_cor <- cor(cibersort_proportions, method = "spearman")
pheatmap(cell_cor,
         color = colorRampPalette(c("blue", "white", "red"))(100),
         main = "Cell Type Correlation")

# 3. 小提琴图（分组比较）
library(ggpubr)
for (ct in c("T cells CD8", "Macrophages M2", "T cells regulatory Tregs")) {
  p <- ggviolin(data.frame(score = cibersort_proportions[, ct],
                           group = sample_info$Group),
                x = "group", y = "score", fill = "group",
                add = "boxplot", add.params = list(fill = "white")) +
    stat_compare_means() +
    labs(title = ct, y = "CIBERSORT Proportion")
  print(p)
}

# 4. 棒棒糖图（多方法结果汇总）
# 比较Disease vs Control, 所有方法所有细胞

# ============================================
# 免疫浸润综合分析流程
# ============================================

# 1. 环境准备
install.packages(c("ggpubr", "pheatmap", "survival", "survminer"))
# immunedeconv 和 IOBR 需从 GitHub 安装（不在 CRAN/Bioconductor 上）
# install.packages("remotes")  # 如未安装
remotes::install_github("omnideconv/immunedeconv")
remotes::install_github("IOBR/IOBR")
# 注意：使用 CIBERSORT 方法需额外在 https://cibersort.stanford.edu 注册获取许可
library(immunedeconv); library(IOBR); library(ggpubr)
library(pheatmap); library(survival); library(survminer)

# 2. 读取数据（TPM矩阵，行=基因symbol，列=样本）
expr <- as.matrix(read.csv("expr_tpm.csv", row.names = 1))
clinical <- read.csv("clinical_data.csv")

# 3. 多方法免疫浸润分析
# CIBERSORT (通过IOBR)
cib <- deconvo_tme(expr, method = "cibersort", arrays = FALSE, perm = 100)
# xCell
xcell <- xCellAnalysis(expr)
# MCP-counter
mcp <- MCPcounter::MCPcounter.estimate(expr, featuresType = "HUGO_symbols")
# EPIC
epic <- EPIC::EPIC(bulk = expr)$cellFractions
# quanTIseq
quant <- deconvolute(expr, "quantiseq", tumor = TRUE)

# 4. 组间差异分析
group <- clinical$Group
for (ct in colnames(cib)[-c(1:2)]) {
  p <- wilcox.test(cib[group == "Tumor", ct], cib[group == "Normal", ct])$p.value
  if (p < 0.05) cat(ct, "P =", formatC(p, digits = 3), "\n")
}

# 5. 生存分析
for (ct in colnames(cib)[-c(1:2)]) {
  med <- median(cib[, ct])
  grp <- ifelse(cib[, ct] > med, "High", "Low")
  cox <- coxph(Surv(clinical$OS_time, clinical$OS_status) ~ grp)
  p <- summary(cox)$coefficients[, 5]
  if (p < 0.05) cat(ct, ": HR =", round(exp(coef(cox)), 2), ", P =", round(p, 4), "\n")
}

# 6. 可视化输出
pdf("immune_infiltration_results.pdf", width = 12, height = 8)
pheatmap(xcell[c("CD8+ T-cells", "CD4+ T-cells", "B-cells", "NK cells",
                  "Macrophages M1", "Macrophages M2", "Tregs"), ],
         scale = "row", annotation_col = data.frame(Group = group))
dev.off()

你需要什么：
├── 22种免疫细胞详细比例 → CIBERSORT
├── 64种细胞类型全面覆盖 → xCell
├── 快速简单的绝对丰度估计 → MCP-counter
├── 统计严格的比例估计 → quanTIseq / EPIC
├── 肿瘤纯度和整体免疫状态 → ESTIMATE
├── 多方法比较 → TIMER3 在线平台（15种算法）
└── 自定义细胞标记 → ssGSEA / GSVA