# ============================================================
# 08_differential_comm.R —— 生态位间差异通讯 / 通讯重编程（对应教程第 13 章）
# 用法：Rscript code/08_differential_comm.R（或 RStudio 逐段运行）
# 前置：已运行 07_niche_cellchat.R，读取 data/processed/cellchat_by_niche.rds
# 核心问题：侵袭前沿 vs 肿瘤核心（示例），哪些配体-受体对/通路被"重编程"？
# ============================================================

library(CellChat)
library(dplyr)

cellchat_list <- readRDS("data/processed/cellchat_by_niche.rds")
niches <- names(cellchat_list)
cat("可用生态位：", paste(niches, collapse = ", "), "\n")

## 1. 选两个生态位做差异比较（可改）
ni_a <- "侵袭前沿"   # 感兴趣生态位
ni_b <- "肿瘤核心"   # 对照生态位
stopifnot(ni_a %in% niches, ni_b %in% niches)
cc_a <- cellchat_list[[ni_a]]
cc_b <- cellchat_list[[ni_b]]

## 2. 提取配体-受体水平通讯表
df_a <- subsetCommunication(cc_a)   # source, target, ligand, receptor, pathway_name, prob, pval
df_b <- subsetCommunication(cc_b)

## 3. 合并并计算差异（log2FC，方向 = A 相对 B）
df <- full_join(
  df_a %>% select(source, target, ligand, receptor, pathway_name, prob) %>%
    rename(prob_a = prob),
  df_b %>% select(source, target, ligand, receptor, pathway_name, prob) %>%
    rename(prob_b = prob),
  by = c("source", "target", "ligand", "receptor", "pathway_name")
) %>%
  mutate(prob_a = replace_na(prob_a, 0),
         prob_b = replace_na(prob_b, 0),
         log2FC = log2((prob_a + 1e-3) / (prob_b + 1e-3)),
         direction = ifelse(log2FC > 1, paste0(ni_a, "上调"),
                    ifelse(log2FC < -1, paste0(ni_b, "上调"), "无显著差异"))) %>%
  arrange(desc(abs(log2FC)))

print(head(df, 20))
write.csv(df, "results/differential_communication.csv", row.names = FALSE)

## 4. 通路活性比较（重编程热点通路）
path_a <- df_a %>% group_by(pathway_name) %>% summarise(strength_a = sum(prob))
path_b <- df_b %>% group_by(pathway_name) %>% summarise(strength_b = sum(prob))
path_df <- full_join(path_a, path_b, by = "pathway_name") %>%
  mutate(strength_a = replace_na(strength_a, 0),
         strength_b = replace_na(strength_b, 0),
         log2FC = log2((strength_a + 1e-3) / (strength_b + 1e-3))) %>%
  arrange(desc(abs(log2FC)))
print(path_df)
write.csv(path_df, "results/differential_pathways.csv", row.names = FALSE)

## 5. 可视化
# 各生态位通路强度对比（水平条形）
barplot(t(as.matrix(path_df[, c("strength_a", "strength_b")])),
        beside = TRUE, names.arg = path_df$pathway_name, las = 2,
        col = c("#C0392B", "#2980B9"), legend.text = c(ni_a, ni_b),
        main = "通路通讯强度对比", ylab = "通讯强度")

# 差异气泡图（选定通路）
netVisual_bubble(cc_a, signaling = head(path_df$pathway_name, 6),
                 title = paste0(ni_a, " 通讯气泡图"))
netVisual_bubble(cc_b, signaling = head(path_df$pathway_name, 6),
                 title = paste0(ni_b, " 通讯气泡图"))

cat("差异通讯分析完成！关键输出：\n")
cat("  - results/differential_communication.csv（L-R 水平差异）\n")
cat("  - results/differential_pathways.csv（通路水平差异）\n")
cat("重编程候选 = 在", ni_a, "显著上调且与免疫抑制/侵袭相关的通路。\n")
