3. Harmony算法原理：基于软聚类迭代矫正嵌入、R语言实现机制

Harmony的核心思想，说白了就是“让数据先混在一起，再根据混血程度做微调”。

这和Seurat CCA那种“先找锚点、再匹配邻居”的硬对齐思路完全不同。我个人更喜欢Harmony的这种柔性方案——它不强行让细胞对号入座，而是让每个细胞带着“可能属于多个批次”的概率去矫正，最后迭代出一个全局平衡的嵌入空间。

为什么这样做有效？今天我们就拆开这个黑盒子，看看里面到底是什么。

3.1 算法核心：软聚类 + 迭代矫正

Harmony的流程可以拆成三步走：

1. 初始化主成分 (PC) 嵌入：先用标准的PCA将高维表达数据降维，得到一个初始的低维空间。
2. 软聚类分配 (Soft Clustering)：在这个空间里，Harmony不是把细胞硬塞进某个确定的簇，而是给每个细胞一个“属于每个簇的概率向量”。
3. 迭代矫正 (Correction)：根据每个簇的批次特异性偏移，调整细胞的嵌入坐标。然后重新聚类、重新矫正，直到收敛。

我用一个简单的SVG图来展示这个循环逻辑，你看一眼就明白了。

3.2 软聚类：为什么不是 K-Means？

你可能要问：直接跑个 K-Means 硬聚类不就行了？何必搞什么“软”的？

嗯，这里有个关键区别。硬聚类要求每个细胞只能属于一个簇——比如“这个细胞是cluster 1的”。但在批次矫正的场景里，很多细胞其实处于边界状态，强行划分会丢失信息，甚至引入错误。

Harmony 用的是软 K-Means 变体，每个细胞对每个簇都有一个隶属度（0到1之间），且所有隶属度加起来等于1。我来给你看看它的数学形式：

# 软聚类更新公式（概念示意，非完整R代码）
# 每个细胞的隶属度 w_ik 由高斯核函数计算

w_ik = exp(-||z_i - μ_k||^2 / 2σ^2)
w_ik = w_ik / sum(w_ij, j=1..K)

# 其中：
# z_i 是细胞i的当前嵌入坐标
# μ_k 是簇k的质心
# σ  是软聚类宽度（sigma参数）
# K  是聚类数（theta参数，默认自动估计）

我在处理一个10x和Smart-seq2的整合项目时，发现边界区（比如中间态T细胞）的细胞，硬聚类会直接把它们推到最近的批次簇里，矫正完反而产生伪差异。换成Harmony的软聚类后，这些细胞以混合身份参与了矫正，结果就自然多了。

3.3 迭代矫正：怎么移动细胞坐标？

矫正这一步才是Harmony的精髓。它不像CCA那样做“全局线性变换”，而是做了个局部加权平均。具体说：

1. 对于每个簇，计算该簇内各批次的 质心偏移（centroid shift）。
2. 每个细胞的矫正向量，就是它隶属的所有簇的质心偏移的加权平均。权重就是软聚类隶属度。
3. 矫正完成后，整个嵌入空间会变得更“混匀”。然后基于新空间重新做软聚类，再矫正。重复。

你可以把每一次迭代理解成一次“微调”——每次只挪一小步，不搞一步到位。这样做的好处是：不会把批次结构彻底抹掉，只是让不同批次的同类细胞靠得更近。

3.4 R语言实现机制：RunHarmony 到底干了什么？

用 harmony 包时，你只需要调用一行 RunHarmony()，但背后机制你需要了解。我拿一个典型调用来拆解：

library(harmony)

# 假设你已经有了 Seurat 对象，并跑过 NormalizeData、FindVariableFeatures、ScaleData、RunPCA
# 标准的 Harmony 调用

harmony_obj <- RunHarmony(
  object = seurat_obj,
  group.by.vars = c("orig.ident", "batch"),  # 需要矫正的批次列
  reduction = "pca",                          # 基于 PCA 嵌入
  dims.use = 1:30,                            # 使用前30个PC
  sigma = 0.1,                                # 软聚类宽度
  theta = NULL,                               # 聚类数，NULL为自动估计
  lambda = 1,                                 # 正则化强度
  max.iter.harmony = 10,                      # 最大迭代次数
  verbose = TRUE                              # 打印收敛日志
)

# 返回的 harmony_obj 是一个 S4 对象
# 它的主要插槽：
#  - H 矩阵：最终矫正后的细胞嵌入（细胞 x 维度）
#  - loss_history：每次迭代的损失值
#  - cluster_assignment：每个细胞最终的硬聚类分配
#  - 其他元数据

RunHarmony 内部执行的动作，我帮你整理成一张清单，你可以对照着看：

步骤	内部操作	关键参数
1. 初始化	从 Seurat 对象中提取 PCA 嵌入矩阵 Z，复制一份作为当前嵌入	dims.use
2. 软聚类	运行软 K-Means，计算每个细胞的隶属度矩阵 W 和簇质心矩阵 M	sigma, theta
3. 计算偏移	对每个簇，计算各批次的质心偏移量（相对于该簇的全局质心）	group.by.vars
4. 矫正嵌入	Z_new = Z - W · Δ（加权质心偏移的负梯度方向）	lambda
5. 判断收敛	计算新旧嵌入的欧氏距离变化，如果变化 < 1e-4 则停止	max.iter.harmony
6. 输出	返回 Harmony 对象，并将矫正嵌入写回 Seurat 对象的 reducctions 槽	—

3.5 你需要注意的几个“坑”

在实际项目中，我几乎每次都会遇到下面这几个问题：

• 矫正过度（Over-correction）：如果 sigma 设得太小（比如 < 0.05）或 max.iter 设得太大，Harmony 会过度混匀不同批次的细胞，导致原本的生物差异被抹除。我一般会在整合后跑一个 RunUMAP 并按细胞类型着色，如果不同批次中同一细胞类型完美重叠甚至分不开，那基本上就是矫正过度了。
• 聚类数选择：theta 参数控制预期有多少个“软簇”。你可能会觉得“我已知有6种细胞类型，就设为6”，但Harmony内部的软簇数量不一定和生物学细胞类型完全对应。我建议大多数情况下保持 theta = NULL，让算法自己根据数据复杂度来选。
• 内存问题：当细胞数量达到几十万甚至上百万时，Harmony 的软聚类矩阵 W（细胞数 × 聚类数）会变得非常大。我处理过一个50万细胞的数据集，R内存占用直冲40GB。我的建议是：先用 Seurat 的标准流程做一次粗分（比如分到50个簇），再在分簇后的层面上跑 Harmony，而不是直接在全部细胞上跑。

3.6 一个小例子：手写一笔 Harmony 的核心计算

为了帮你理解，我手动模拟一下简化版的矫正逻辑。假设只有3个细胞、2个批次、2个簇：

# 伪代码演示，不是真正可运行的 R 代码
# 初始 PCA 嵌入 Z（3个细胞 x 2维）
Z <- matrix(c(
  1.0, 2.0,   # 细胞1，批次A
  1.2, 1.8,   # 细胞2，批次A
  3.0, 1.0    # 细胞3，批次B
), nrow=3, byrow=TRUE)

# 第1次迭代：软聚类得到隶属度 W
# 假设两个簇的质心：μ1=(1.1,1.9), μ2=(2.8,1.0)
# 细胞1的隶属度：w11=0.95, w12=0.05（很像属于簇1）
# 细胞3的隶属度：w31=0.10, w32=0.90（很像属于簇2）

# 批次A的簇1质心偏移 Δ=(1.1,1.9) - 全局质心
# 批次B的簇2质心偏移 等
# 细胞1的矫正向量 = 0.95 * Δ + 0.05 * 其他偏移
# 实际上，每个细胞的矫正量是其隶属度加权的线性组合

# 更新后的 Z_new 会让同类型细胞靠得更近

我建议你可以在自己的数据上尝试调低 lambda 值（默认为1），比如设为0.5，看看UMAP图的变化。你可能会发现某些细胞类型的分离度变高了——但也要小心，这也许意味着矫正不足。