# 第 13 章 通讯重编程：生态位间差异分析

> 本章对应脚本：`code/08_differential_comm.R`。
> 前置：第 12 章的 `cc_list`（各生态位的 CellChat 对象列表），以及对应的 scanpy 对象 `adata`。

本章是全书的“高潮”。前两章我们分别拿到了全局通讯网络和生态位内通讯网络；本章要做的是把它们之间的差异“讲成故事”——这就是“通讯重编程”（communication reprogramming）。请注意，“重编程”在这里不是修辞，而是可以操作化的分析对象：生态位之间在通讯数量、强度与通路活性上的系统性差异。读完本章，你将能把“侵袭前沿被重编程为免疫抑制枢纽”这样一句话，变成一张张有数据支撑的图表。

## 本章目标

1. 掌握“通讯重编程”的操作化定义（数量、强度、通路活性三个层面）
2. 学会生态位间差异通讯的比较方法（compareInteractions / 手动汇总）
3. 学会识别差异信号通路（rankNet 对比、通路活性热图、netVisual_diff）
4. 学会“差异配体-受体对”分析，并讲清 log2FC 与方向的含义
5. 建立“三重证据闭环”（空间共定位 + 空间通讯活性 + 表达证据），并用 Sankey 图把重编程故事讲完整

## 1 什么是“通讯重编程”：操作化定义

“重编程”的生物学直觉是：肿瘤在塑造微环境，微环境反过来塑造肿瘤，最终整个组织的通讯格局被系统性改写。要让这个概念可检验，必须把它翻译成可计算的指标。本书的操作化定义是：**通讯重编程 = 生态位之间通讯网络在数量（count）、强度（weight）、通路活性（pathway activity）上的系统性差异**。

- **数量**：显著 L-R 互作的条数，衡量网络“规模”；
- **强度**：通讯概率之和，衡量网络“音量”；
- **通路活性**：每个信号通路的总通讯概率，衡量网络“曲目”。

三者合起来，一个生态位的“通讯指纹”就完整了。图 12-1（见第 12 章）展示的正是第一个层面：不同生态位的通讯数量与强度差异显著。注意，重编程是**系统级差异**，与“某一条 L-R 对在不同生态位强弱不同”这种单点差异不同——前者是格局，后者是细节；本章先看格局（第 2、3 节），再看细节（第 4 节）。这套定义与第 5 章登记的研究假设一一对应：H1（不同生态位的通讯网络存在系统性差异）对应数量与强度层面，H2（侵袭前沿存在免疫抑制性通讯重编程）对应通路活性层面，H3（分区分析能揭示更多空间异质性信号）则贯穿第 12 章与本章的比较逻辑——本章的每一张图，都是在为某条假设提供证据。

## 2 差异通讯数量/强度比较

### 2.1 compareInteractions：官方比较函数

当两个生态位（或两个样本）的细胞类型集合完全一致时，可以先把对象合并，再用官方函数比较：

```r
# 前提：两个对象的细胞类型集合完全相同（如"样本 A vs 样本 B 的侵袭前沿"）
# 若集合不同（本书生态位常有此情况），请直接跳到 2.2 手动汇总
cc_merged <- mergeCellChat(cc_list[c("侵袭前沿", "肿瘤核心")],
                           add.names = c("侵袭前沿", "肿瘤核心"))
# 比较总通讯数量（count）与总强度（weight）；group.1/group.2 填合并时的条件名
compareInteractions(cc_merged, group.1 = "侵袭前沿", group.2 = "肿瘤核心",
                    measure = c("count", "weight"))
```

运行结果解读：输出两张并排柱状图（正是图 12-1 的风格来源），直观显示“前沿 vs 核心”的总通讯数量与总强度差异。如果报“细胞群数量不一致”错误，说明两个对象细胞类型集合不同——这正是 2.2 节手动汇总的用武之地。

### 2.2 手动汇总：更通用的做法

从第 12 章的 `cc_list` 出发，逐生态位统计通讯数量与总强度，一步到位：

```r
niche_summary <- data.frame(
  niche          = names(cc_list),
  n_interactions = sapply(cc_list, function(cc) sum(cc@net$count)),  # 显著互作总条数
  total_strength = sapply(cc_list, function(cc) sum(cc@net$weight)), # 通讯概率之和
  row.names = NULL
)
print(niche_summary)
# 转成长表后画柱状图（对应图 12-1 的"数量"面板）
library(tidyr); library(ggplot2)
niche_long <- pivot_longer(niche_summary,
                           cols = c(n_interactions, total_strength),
                           names_to = "measure", values_to = "value")
ggplot(niche_long, aes(x = niche, y = value, fill = measure)) +
  geom_col(position = "dodge") + theme_minimal() +
  labs(x = "生态位", y = "数值", title = "各生态位通讯数量与强度（模拟数据演示）")
```

运行结果解读：通常看到侵袭前沿的 `n_interactions` 与 `total_strength` 最高。注意这只是描述性比较；要下“显著更高”的结论，需要生物学重复与统计检验（见第 8 节注意事项）。

## 3 差异信号通路识别

### 3.1 rankNet 对比

`rankNet` 给每个生态位输出“通路贡献排名”。把 4 个生态位的排名并列，就能看出“曲目”差异：

```r
# 逐个生态位看通路排名（对比时注意各自排名前列的通路）
for (n in names(cc_list)) {
  cat("==== 生态位：", n, " ====\n")
  print(rankNet(cc_list[[n]]))   # 返回排序后的通路与相对贡献
}
```

运行结果解读：如果侵袭前沿的榜首是 TGFb/MIF，免疫浸润区是 CXCL/CCL，肿瘤核心是 VEGF——“曲目”差异一目了然，这就是重编程的故事素材。

### 3.2 通路活性热图：生态位 × 通路

把每个生态位的通路总概率拼成矩阵画热图，一张图看完全部差异（对应图 13-2，模拟数据演示）：

```r
# 从每个生态位取"通路级总概率"，拼成生态位 × 通路矩阵
pw_list <- lapply(cc_list, function(cc) apply(cc@netP$prob, 3, sum))
all_pw  <- unique(unlist(lapply(pw_list, names)))
pathway_mat <- sapply(pw_list, function(v) v[all_pw])
pathway_mat[is.na(pathway_mat)] <- 0
pathway_mat <- t(pathway_mat)                              # 行 = 生态位，列 = 通路
pathway_mat <- pathway_mat[, colSums(pathway_mat) > 0]     # 去掉全零通路
pheatmap::pheatmap(pathway_mat, cluster_cols = FALSE,
                   main = "生态位 × 通路活性（模拟数据演示）")
```

运行结果解读：热图里能直接看到“前沿高活性”的免疫抑制与侵袭相关通路（TGFb、PD-L1、MIF、GALECTIN），与图 13-2 的风格一致。若提示没有 `pheatmap` 包，先运行 `install.packages("pheatmap")`。

![图 13-2 生态位 × 信号通路活性热图（模拟数据演示）。行是通路（TGFb、PD-L1、MIF、GALECTIN、VEGF、CXCL、ECM、IL10、CCL），列是生态位；可见侵袭前沿在免疫抑制与侵袭相关通路上活性最高。](figures/fig15_pathway_heatmap.png){width=85%}

### 3.3 netVisual_diff：网络级差异

如果走 2.1 的合并路线，还能在网络上直接看差异：

```r
# 在合并对象上可视化两个生态位的网络差异（红色 = 前沿更强，蓝色 = 核心更强）
netVisual_diff(cc_merged, weight.edge = TRUE, measure = "weight")
```

## 4 差异配体-受体对分析：谁在具体地“多说/少说”

通路是“曲目”，L-R 对才是“乐句”。跨生态位比较每条 L-R 对的通讯概率，并计算 log2FC：

```r
library(ggplot2)
# 提取每个生态位"每条 L-R 对的全网络平均通讯概率"
pair_prob <- function(cc) apply(cc@net$prob, 3, mean)   # 3D 数组按第三维（L-R 对）求平均
prob_front <- pair_prob(cc_list[["侵袭前沿"]])
prob_core  <- pair_prob(cc_list[["肿瘤核心"]])
pairs <- union(names(prob_front), names(prob_core))
df <- data.frame(
  pair  = pairs,
  front = prob_front[pairs],
  core  = prob_core[pairs]
)
df[is.na(df)] <- 0
# log2FC > 0：前沿更强（上调）；log2FC < 0：核心更强（前沿下调）
df$log2FC <- log2((df$front + 1e-4) / (df$core + 1e-4))   # +1e-4 防止除零
df$direction <- ifelse(df$log2FC > 0.5, "前沿上调",
                ifelse(df$log2FC < -0.5, "前沿下调", "无显著差异"))
head(df[order(-abs(df$log2FC)), ], 10)
# 画成点图（对应图 13-1 的左面板风格）
top15 <- head(df[order(-abs(df$log2FC)), ], 15)
ggplot(top15, aes(x = log2FC, y = reorder(pair, log2FC), color = direction)) +
  geom_point(size = 3) + geom_vline(xintercept = 0, linetype = 2) +
  scale_color_manual(values = c("前沿上调" = "#C44E52", "前沿下调" = "#4C72B0",
                                "无显著差异" = "grey50")) +
  theme_minimal() + labs(x = "log2FC（侵袭前沿 / 肿瘤核心）", y = NULL,
                         title = "生态位间差异配体-受体对（模拟数据演示）")
```

**log2FC 与方向的含义（务必讲清楚）**：log2FC 是对数尺度的比值——log2FC = 1 表示前沿的平均通讯概率是核心的 2 倍，log2FC = -1 表示只有一半。这里的“方向”指的是“谁更强”：红色（前沿上调）表示该 L-R 对在侵袭前沿的通讯概率系统性高于肿瘤核心，例如 PD-L1（CD274-PDCD1）、MIF-CD74、TGFB1-TGFBR2；蓝色（前沿下调）则相反，例如只活跃于核心的血管生成相关对。图 13-1 就是这种分析的标准呈现：左面板气泡图对比两生态位的通讯概率（实心 = 前沿、空心 = 核心，红色 = 前沿上调、蓝色 = 前沿下调），右面板是 Δ 通讯概率条形热图。

![图 13-1 侵袭前沿与肿瘤核心的差异通讯（模拟数据演示）。左：气泡图，实心圆代表侵袭前沿、空心圆代表肿瘤核心，红色表示前沿上调、蓝色表示前沿下调；右：Δ 通讯概率条形热图。](figures/fig14_diff_comm.png){width=90%}

注意：这里的 log2FC 是“均值之比”的描述性指标，本身没有 p 值。严谨做法是把 CellChat 给出的每条 L-R 对的排列检验 p 值（存在 `cc@net$pval`）纳入筛选（比如只保留至少在一个生态位显著的对），再结合重复样本做正式统计检验（见第 8 节）。

## 5 空间活性验证：三重证据闭环

CellChat 的推断本质是“表达层面的共现”，它不知道配体细胞与受体细胞是否真的相邻。本章的可信度建立在三重证据闭环上（对应课题设计中的“空间共定位证据增强”）：

- **证据 1（空间共定位）**：Squidpy 邻域富集，检验“配体细胞类型与受体细胞类型在空间上相邻”；
- **证据 2（空间通讯活性）**：COMMOT 把空间距离纳入通讯推断，输出每个 spot 的通讯活性；
- **证据 3（表达证据）**：在目标生态位内直接检查配体/受体的表达。

```python
# ============ 第 13 章 Python 代码：空间活性验证 ============
import scanpy as sc
import squidpy as sq
import commot as ct

# 前置：adata 是 scanpy 对象（第 9–10 章结果），
# adata.obs["celltype"] 为细胞类型，adata.obs["niche"] 为生态位

# ---- 证据 1：Squidpy 邻域富集（配体细胞与受体细胞是否相邻）----
sq.gr.spatial_neighbors(adata, coord_type="generic", n_neighs=6)
sq.gr.nhood_enrichment(adata, cluster_key="celltype")
z = adata.uns["celltype_nhood_enrichment"]["zscore"]   # z > 0 = 倾向共定位
print(z.loc["CAF", "肿瘤细胞"])    # 把 CAF/肿瘤细胞换成你数据里的细胞类型名
sq.pl.nhood_enrichment(adata, cluster_key="celltype")  # 全矩阵热图（风格参见第 10 章）

# ---- 证据 2：COMMOT 空间通讯活性（仓库 https://github.com/zcang/COMMOT）----
# dis_thr：通讯距离阈值（μm）；pathway_sum：聚合到通路级
ct.tl.spatial_communication(adata, database_name="CellChat", dis_thr=250,
                            heteromeric=True, pathway_sum=True)
# 绘制指定通路的空间活性（不同 COMMOT 版本绘图 API 略有差异，以官方文档为准）
ct.pl.plot_cell_communication(adata, cell_type_key="celltype",
                              signaling_type="pathway", communication_key="TGFb")

# ---- 证据 3：表达证据（只看向侵袭前沿生态位）----
sc.pl.dotplot(adata[adata.obs["niche"] == "侵袭前沿"],
              var_names={"配体": ["TGFB1", "MIF", "CD274"],
                         "受体": ["TGFBR2", "CD74", "PDCD1"]},
              groupby="celltype")
```

运行结果解读：三重证据的判读规则是“相互印证”——CellChat 说“前沿有 TGFb 通讯”，Squidpy 说“CAF 与肿瘤细胞确实相邻”，COMMOT 说“TGFb 活性在空间上确实落在前沿”，dotplot 说“配体与受体确实在前沿表达”。四者一致，结论才站得住；任何一环缺失（比如共定位不显著），都要在论文里如实说明，或把该信号降级为“候选”。黄金标准仍是实验验证（中和抗体、基因敲除），计算证据只能提高置信度。

## 6 重编程故事化：Sankey 流量图

差异分析会产出几十上百条信号，直接堆在论文里没人读得懂。Sankey 图（图 13-3，模拟数据演示）把“来源生态位 → 重编程通路 → 靶细胞类型”串成一条河流，带宽 ∝ 通讯活性，让“谁通过哪条通路影响了谁”一目了然：

```python
# ============ Sankey 流量图（plotly；数据来自第 4–5 节汇总表）============
import plotly.graph_objects as go

# 手工整理三元组：来源生态位 → 通路 → 靶细胞（数值 = 归一化通讯活性，模拟数据演示）
links = [
    ("侵袭前沿",   "TGFb",  0.35), ("侵袭前沿",   "MIF",   0.28), ("侵袭前沿",   "PD-L1", 0.22),
    ("免疫浸润区", "CXCL",  0.30), ("免疫浸润区", "CCL",   0.25),
    ("肿瘤核心",   "VEGF",  0.40),
]
nodes = sorted({s for s, _, _ in links} | {t for _, t, _ in links})
idx   = {n: i for i, n in enumerate(nodes)}
fig = go.Figure(go.Sankey(
    node=dict(label=nodes, pad=15, thickness=20),
    link=dict(source=[idx[s] for s, _, _ in links],
              target=[idx[t] for _, t, _ in links],
              value=[v for _, _, v in links])))
fig.update_layout(title="通讯重编程信号流（模拟数据演示）")
fig.show()
```

![图 13-3 通讯重编程信号流 Sankey 图（模拟数据演示）。左列是来源生态位，中列是重编程信号通路，右列是靶细胞类型；带宽与通讯活性成正比，直观呈现“哪个生态位通过哪条通路影响哪些靶细胞”。](figures/fig16_reprogramming.png){width=90%}

解读要领：Sankey 的“流向”必须严格忠于数据——来源生态位是“发送者细胞所在生态位”，靶细胞是“接收者细胞类型”，中间是通路。不要为了画面好看而添加数据里不存在的连线；每条连线的数值应能在第 4、5 节的结果里找到出处。本章产出的这些“素材图”（图 13-1 至图 13-3）到第 14 章会按论文 Figure 的版式重新编排（版式示例见 fig17_paper_figure）。

## 7 结果解读示例：一段示范性 Results 文字

下面是一段“结果-解读”示范（数值对应本书模拟数据演示）：

> 与肿瘤核心相比，侵袭前沿生态位的显著通讯数量与总强度均更高（图 12-1）。通路活性分析显示，侵袭前沿在 TGFb、PD-L1、MIF 与 GALECTIN 通路上活性最高，而肿瘤核心以 VEGF 通路为主（图 13-2）。配体-受体对水平上，TGFB1-TGFBR2、CD274-PDCD1、MIF-CD74 在侵袭前沿显著上调（log2FC > 1；图 13-1），且配体细胞（CAF、肿瘤细胞）与受体细胞（肿瘤细胞、免疫细胞）在空间上共定位（Squidpy 邻域富集 z > 2），COMMOT 显示 TGFb 通讯活性集中于侵袭前沿。这些结果提示，侵袭前沿存在以 TGF-β、PD-L1 与 MIF 为特征的免疫抑制性通讯重编程：TGF-β 驱动肿瘤 EMT 与免疫逃逸，PD-L1/PD-1 直接抑制 T 细胞杀伤，MIF 招募并极化免疫抑制细胞——三者协同将侵袭前沿塑造为“免疫抑制 + 侵袭”的枢纽（图 13-3）。

注意这段文字的两处纪律：一是每个“提示”背后都对应一张图与一个统计量；二是结论句用的是“提示”而非“证明”——第 8 节的注意事项就是为这句话兜底。

## 8 注意事项

1. **差异 ≠ 因果**。生态位间的通讯差异是相关性发现：可能是通讯驱动了生态位（肿瘤重编程微环境），也可能是生态位环境（缺氧、ECM）驱动了通讯，还可能两者互相强化。要谈因果，需要扰动实验（中和抗体、敲除、类器官共培养）或至少时序数据。
2. **多生态位比较要做统计校正**。我们进行了成百上千次比较（通路 × 生态位、L-R 对 × 生态位），多重检验会膨胀假阳性率，默认使用 BH-FDR（Benjamini-Hochberg）校正，并在方法部分写明校正方法与阈值。
3. **重复样本是底线**。单个切片只能支持“探索性假说”；任何“重编程”结论都应以多个生物学重复（多个患者/多个切片）为准，最好在独立队列中验证。
4. **组成差异的混淆**。生态位间通讯差异可能只是细胞组成差异的投影（CAF 多的地方当然 CAF 信号强），不一定是“通讯程序被重编程”。严谨做法是结合去卷积丰度做组成匹配或校正，并在讨论中说明这一局限。
5. **分辨率与混合信号**。spot 的混合信号会让“发送者/接收者”标签带噪声；高分辨率平台（Visium HD、Xenium）或单细胞-空间联合分析可以缓解。

## 本章小结

本章完成了从“差异”到“重编程”的升级：用数量、强度、通路活性三个层面操作化定义重编程；用 compareInteractions / 手动汇总比较差异；用通路活性热图与 netVisual_diff 识别差异通路；用 log2FC 讲清差异 L-R 对的方向；用 Squidpy + COMMOT + 表达证据建立三重证据闭环；最后用 Sankey 图把故事讲完整。至此，全书分析主线（数据 → 生态位 → 通讯 → 重编程）已经跑通，下一章进入论文级图表与写作环节。

## 练习与思考

1. “通讯重编程”的操作化定义包含哪三个层面？为什么说它是“系统级差异”而不是单点差异？
2. （动手题）用 `cc_list` 计算侵袭前沿 vs 肿瘤核心的差异 L-R 对（log2FC），列出上调最强的 5 对，并查它们在 CellChatDB 中属于哪些通路。
3. （动手题）在侵袭前沿内运行 COMMOT（dis_thr = 250）并绘制 TGFb 空间活性，与 CellChat 的 TGFb 网络结论对照，写一段“三重证据”判读。
4. （思考题）如果侵袭前沿的 CAF 比例远高于肿瘤核心，TGFb 的差异还能算“重编程”证据吗？你会怎么排除组成混淆？
5. （思考题）为什么单切片数据的“重编程”结论不能直接写进论文的结论部分？需要补什么？
