2. Harmony算法原理（一）：迭代聚类与校正的数学直觉

这一章我打算跟大家聊聊Harmony的核心数学直觉。别看它名字叫“和谐”，其实背后的逻辑特别直白——说白了就是：聚类、看偏差、挪一下、再聚类。

我第一次看到这个思路的时候，心里想的是：就这么简单？后来用多了才发现，简单的循环往往最有效。嗯，咱们慢慢拆开看。

2.1 从问题出发：为什么要迭代？

单细胞数据整合最大的麻烦是批次效应。不同实验、不同测序仪、不同操作员，都会给数据带上“烙印”。

你想想，如果我们直接拿原始数据去降维聚类，相同的细胞类型在不同的批次里可能会被分到不同的簇。这就是批次效应的“罪恶”。

Harmony的做法很聪明：它不试图一次性抹掉所有差异，而是通过迭代逐步修正。每轮迭代只做两件事：

• 软聚类：给每个细胞分配隶属度（属于每个簇的概率）
• 校正：根据批次隶属度差异，把细胞往“中心”拉一拉

这样循环下去，批次效应就被一点点削掉了。我之前在项目里试过一次跑100轮，发现其实20轮左右就够了，再跑梯度基本不动。

2.2 核心思想：软聚类与校正

经典的K-means做的是硬聚类——一个细胞只能属于一个簇。但Harmony用的是软聚类，学名叫模糊C均值聚类（Fuzzy C-Means）的变体。

每个细胞对每个簇都有一个隶属度 R_ik，数值在0到1之间，加起来总和为1。比如某个T细胞可能70%属于簇1，20%属于簇2，10%属于簇3。这其实更符合真实情况——细胞类型之间的边界很多时候是模糊的。

校正怎么做呢？当Harmony发现某个批次里的细胞在某个簇里扎堆太多，就会把这个批次的所有细胞往全局中心推一推。具体的数学公式我先不说太细，你只需记住一个直觉：

为什么会这样？因为批次效应本质上是一个加性偏差，用向量减法就能处理。当然实际还要考虑数据在低维空间（比如PCA）的表现形式。

2.3 数学直觉：期望最大化视角

懂EM算法（期望最大化）的朋友应该已经看出来了，Harmony的迭代过程本质上是个EM框架：

• E步（期望步）：根据当前的校正参数，估计每个细胞的簇隶属度
• M步（最大化步）：固定隶属度，更新批次校正参数，使得“混合数据”中的批次混叠最小化

这里我习惯用一个比喻：好比你在调一个多声道音响，每个批次相当于一个声道，你需要在每个频段（簇）上找到那个声道的声音有多大，然后把过强的那一轨音量降下来。来回试几次，各个声道的音量就均衡了。

迭代步骤	对应EM	数学直觉
软聚类	E步	当前参数下，细胞属于各簇的概率
计算批次偏移	M步的一部分	每个批次在该簇中的平均偏差
校正数据	M步另一部分	把偏差从细胞坐标中减去

2.4 算法流程简述

把上面的直觉串起来，就是下面这张图的核心逻辑。我个人非常建议大家先看懂这张图，再去啃代码。很多细节你一看流程就明白了。

2.5 关键参数怎么理解

Harmony最常用的参数就三个。我列个表方便你对照着看：

参数名	缺省值	我的理解
lambda	1	正则化强度，越大校正越保守。我用0.5处理过高度异质的批次，效果比默认好
sigma	0.1	软聚类时隶属度分布的宽窄。越小聚类越“硬”，越大越“软”。一般不需要改
nclust	自动	聚类中心数量，如果数据量小可以手动设少一点（比如10~20）

举个实际例子。我之前一个项目有3个批次，每个批次包含8万个细胞，原始数据用PCA降到50维。我直接用HarmonyMatrix跑，nclust用默认值（等于min(50, n_batches*7)），迭代了大概28轮就收敛了。整合后UMAP图上的批次混叠基本消失。

2.6 一个小坑：初始化的敏感性

Harmony用K-means++做初始化。理论上不同初始种子会得到不同结果，但后续迭代会收敛到相近的解。我个人的习惯是设置random.seed为固定值（比如42L），保证可复现。如果你不关心可复现性，或者想验证稳定性，可以多跑几次看看结果是否一致。

总的来说，Harmony的迭代思路非常稳健。它不像一些方法需要先找“锚点”或做批次间配对，而是通过全局的软聚类和校正，让所有批次在迭代中自己“协商”出一个一致的结构。这种“民主协商”式的整合，在实际数据上表现相当靠谱。