3. K-means聚类原理与实现归纳总结

K-means 应该是空间转录组分析里最常用的无监督方法了。为什么？因为它简单、快，结果也还算直观。我当年刚开始接触 Visium 数据时，第一个试的就是 K-means，跑完一看聚类图——嗯，组织学结构基本能对上，当时就觉得这门算法有点东西。

但别急着上手。你得先搞清楚它到底在算什么，以及为什么有些坑踩进去就爬不出来。

3.1 算法核心——一句话讲明白

K-means 的目标很简单：把高维空间里的点分成 K 堆，每堆内部尽量紧凑，堆与堆之间尽量疏远。说白了就是“物以类聚，人以群分”，只不过是用数学方式算出来的。

它的迭代逻辑只有四步：

  1. 初始化：随机选 K 个点作为初始聚类中心。
  2. 分配：每个样本点找离自己最近的中心，归入那个簇。
  3. 更新：计算每个簇的均值，作为新的中心。
  4. 重复：跳到第 2 步，直到中心不再变化或达到最大迭代次数。

你想想看，整个过程就是一个“猜中心 → 划分子 → 算平均 → 再猜”的循环。数学上，它是在最小化每个样本到其所在簇中心点距离的平方和——也就是惯性（inertia）。

K-means 的损失函数：
J = Σ (x - μ_k)²
其中 μ_k 是第 k 个簇的中心，x 属于该簇。最小化这个值，就是让簇内点尽可能密集。

3.2 距离度量与数据标准化

默认都用欧氏距离。但我在 Visium 项目中遇到过一个问题：基因表达矩阵里不同基因的方差差了几个数量级，直接用原始值跑 K-means，高方差基因会主导聚类结果。你懂吧，那结果基本就是“谁数值大谁说了算”。

所以建议先做标准化：要么 Z-score 按基因缩放，要么用对数归一化 + 标准化。我个人习惯在降维到前 50 个 PCA 主成分后再跑 K-means，这样既能去噪声又能降维，速度也快不少。

场景推荐处理原因
原始表达量LogNormalize + Scale消除文库大小差异，均衡基因方差
PCA 降维后不需要再标准化PCA 本身已经对主成分做了方差归一
空间坐标保持原始坐标或 min-max 归一避免组织形状被缩放扭曲

3.3 手动实现就是最好的理解方式

我建议每个学聚类的人都手写一遍 K-means。别怕麻烦，写一次比看十遍原理都管用。下面这个实现去掉了花哨的加速，保留了最核心的逻辑：

import numpy as np

def k_means(X, K, max_iter=100, tol=1e-4):
    # 随机初始化中心
    idx = np.random.choice(len(X), K, replace=False)
    centers = X[idx]
    
    for i in range(max_iter):
        # 分配：计算每个点到所有中心的欧氏距离
        distances = np.linalg.norm(X[:, None] - centers, axis=2)
        labels = np.argmin(distances, axis=1)
        
        # 更新中心
        new_centers = np.array([X[labels == k].mean(axis=0) for k in range(K)])
        
        # 检查收敛
        shift = np.linalg.norm(new_centers - centers)
        if shift < tol:
            break
        centers = new_centers
        
    return labels, centers

这段代码你拿去跑小样本数据完全能用。但在 Visium 数据上——通常有几千个 spot、两万个基因——别这么干。建议直接用 sklearn.cluster.KMeans，它优化了初始中心选择（k-means++）和迭代加速，还支持并行。

sklearn 用法：
from sklearn.cluster import KMeans
km = KMeans(n_clusters=8, init='k-means++', n_init=10, random_state=42)
labels = km.fit_predict(X_pca)

3.4 在 Visium 数据上的实战经验

Visium 每个 spot 是一个 55μm 的圈，每个圈里混合了几到十几个细胞。所以 spot 的表达其实是个混合信号。这时候做 K-means，你会发现斑点之间的边界通常对应组织层的分界。嗯，很好理解——不同层域的细胞类型不同，表达谱自然有差异。

但我踩过一个坑：直接用全部基因（~20k）跑 K-means，结果每次运行聚类都不一样。为什么？高维空间下距离度量失效（“维度灾难”），加上初始中心随机，结果极不稳定。

解决办法：

  • 先做特征选择：选高变基因（top 2000~5000）。
  • 再 PCA 降到 30~50 维。
  • 最后跑 K-means，设置 n_init=20 以上，选最优解。

有一回我忘了对坐标做对齐，把两个不同切片拼一起跑了 K-means——结果聚出来的簇完全按照切片ID分开，毫无生物学意义。你想想，这多误人子弟。

3.5 K 值怎么选？肘部法和轮廓系数

K-means 最头疼的就是要自己定 K。没有标准答案。我习惯先用肘部法画惯性曲线：

import matplotlib.pyplot as plt
from sklearn.cluster import KMeans

inertias = []
K_range = range(2, 15)
for k in K_range:
    km = KMeans(n_clusters=k, init='k-means++', n_init=10, random_state=42)
    km.fit(X_pca)
    inertias.append(km.inertia_)

plt.plot(K_range, inertias, 'bo-')
plt.xlabel('K')
plt.ylabel('Inertia')
plt.show()

找那个"拐点"——惯性下降速度突然变缓的地方。但说实话，很多 Visium 数据并不存在明显拐点，这时候我会配合轮廓系数（silhouette score）一起看。

指标适用场景注意
惯性（inertia）初筛，看曲线下降趋势越大不一定越差，得看拐点
轮廓系数评估簇内紧密度和簇间分离度值在[-1,1]，越大越好
生物学先验如果知道组织有几层手动预设 K，再微调

千万别迷信某个单一指标选 K。我记得有次肘部法给出的 K=4，轮廓系数却给了 K=7。最后结合病理图像一看，K=7 把白质里几个亚结构分开了，更有意义。所以，多用几个视角，包括空间可视化检查。

3.6 缺点和避坑清单

  • 对初始中心敏感：每次运行结果可能不同。用 k-means++ 初始化 + 多次重复取最优。
  • 假设簇是球形且大小相近：Visium 中有些细胞类型很少（比如免疫细胞浸润），这些簇会被其他大簇吞并。遇到这种情况，可以考虑用 DBSCAN 或层次聚类。
  • 无法处理离群点：个别离群的 spot 会被强行归入某个簇，干扰边界。我一般先做一轮 DBSCAN 过滤掉异常点后再跑 K-means。
  • 尺度敏感：一定要对特征进行标准化，否则高数值特征主导结果。
  • K 值难以确定：上面已经说了，结合多种方法，并做最终的空间可视化确认。

3.7 本章知识结构

下面这张图帮你把 K-means 在 Visium 上的完整流程串起来，从数据准备到聚类后验证。我画的时候特别强调了“标准化”和“降维”这两步——很多人跳过它们直接跑 K-means，然后来问我为什么结果不好。

K-means 在 Visium 空间转录组中的完整流程 原始计数矩阵 (Spots × Genes) 标准化 + 高变基因筛选 (LogNormalize + Scale) PCA 降维 (30~50 主成分) K 值选择 肘部法 / 轮廓系数 生物学先验 K-means 聚类 k-means++ 初始化 n_init=20 取最优 结果验证与可视化 空间映射 + 差异基因分析 关键：标准化和降维是 K-means 在 Visium 数据上稳定工作的前提，跳过这两步结果必然漂移

K-means 是空间转录组分析工具箱里的“瑞士军刀”——好用、通用，但别指望它处理所有问题。下一把刀是层次聚类，它不需要预设 K，而且能画出好看的树状图。不过那是另一个话题了，先把 K-means 练熟再说。

← 上一章 📖 返回目录 下一章 →