用MATLAB进行K-means聚类分析时，如何选择合适的K值？

冤大头编程之路
117 阅读9分钟

在MATLAB中为K-means聚类选择合适的K值，核心是平衡“数学最优（簇内紧凑、簇间分离）”和“生物学意义（贴合基因表达的功能模块）”（如上调、下调、无差异基因簇）。针对基因表达数据的高维、有生物学模块特征，以下从「核心原则、量化方法、实战案例、注意事项」展开，每个方法附MATLAB可运行代码和适配解读。

一、核心原则（适配基因表达数据）

选择K值需满足两个条件：

  1. 数学层面：簇内平方和（SSE）尽可能小（簇内紧凑），簇间距离尽可能大（簇间分离）；
  2. 生物学层面：K值贴合基因的真实功能模块（如上调/下调/无差异模块、通路相关模块），避免过度追求数学最优而脱离生物学意义。

二、量化选择方法（按优先级+适配性排序）

方法1：肘部法则（Elbow Method）—— 最常用，快速初筛

原理：计算不同K值对应的类内平方和（SSE）（所有样本到其簇中心的平方距离之和），SSE随K增大递减，当K超过“肘部点”后，SSE下降速率骤减（边际收益递减），该点即为最优K。 适配场景：无生物学先验、快速确定K值范围（基因表达数据初筛）。

% ===================== 肘部法则实现（基因表达数据） =====================
% 1. 数据准备：预处理后的差异基因表达矩阵（基因×样本，行=基因，列=样本）
rng(1); % 固定随机种子
n_genes = 100;  % 100个差异基因
n_samples = 60; % 60个样本
expr_norm = zscore(randn(n_samples, n_genes)); % 标准化表达矩阵（样本×基因）
gene_expr = expr_norm'; % 转置为基因×样本（聚类对象：基因）

% 2. 定义K值范围（通常测试1~10，基因表达多为2~5）
k_range = 1:10;
sse = zeros(length(k_range), 1); % 存储不同K的SSE

% 3. 循环计算每个K的SSE
for i = 1:length(k_range)
    k = k_range(i);
    % K-means聚类（用K-means++初始，避免局部最优）
    [~, ~, sumd] = kmeans(gene_expr, k, ...
        'Start','plus', ...          % 最优初始中心
        'Distance','euclidean', ...  % 基因表达用欧氏距离
        'Replicates',5);             % 5次重复，取最优结果
    sse(i) = sum(sumd); % 累计类内平方和
end

% 4. 可视化肘部曲线，标注最优K
figure('Color','w');
plot(k_range, sse, 'o-', 'LineWidth',1.5, 'MarkerSize',8);
xlabel('聚类数K'); ylabel('类内平方和（SSE）');
title('肘部法则：基因表达数据K值选择');
grid on;
% 手动标注肘部点（示例：K=3，对应基因的3个功能模块）
hold on; 
plot(3, sse(3), 'rs', 'MarkerSize',10); % 红色标注肘部点
text(3.2, sse(3), '最优K=3（肘部点）', 'FontSize',12);
hold off;

结果解读：示例中K=3时，SSE从陡峭下降变为平缓，是典型的“肘部点”，对应基因的上调、下调、无差异3个模块。

方法2：轮廓系数法（Silhouette Coefficient）—— 更严谨，兼顾簇内/簇间

原理：每个样本的轮廓系数= (簇间平均距离 - 簇内平均距离) / max(簇间平均距离, 簇内平均距离)，取值[-1,1]：

  • 接近1：样本分配合理（簇内紧、簇间分）；
  • 接近-1：样本分配错误；
  • 接近0：簇重叠。 最优K：所有K值中，轮廓系数均值最大的K（基因表达数据优先选此方法，避免肘部法则的主观判断）。
% ===================== 轮廓系数法实现 =====================
% 1. 复用上述基因表达数据（gene_expr）和K范围（k_range）
sil_mean = zeros(length(k_range), 1); % 存储不同K的轮廓系数均值

% 2. 循环计算每个K的轮廓系数
for i = 1:length(k_range)
    k = k_range(i);
    if k == 1 % 轮廓系数仅适用于K≥2
        sil_mean(i) = 0;
        continue;
    end
    [cluster_idx, ~] = kmeans(gene_expr, k, 'Start','plus', 'Replicates',5);
    % 计算轮廓系数
    sil = silhouette(gene_expr, cluster_idx, 'euclidean');
    sil_mean(i) = mean(sil); % 取均值（避免单个样本干扰）
end

% 3. 可视化轮廓系数曲线，标注最优K
figure('Color','w');
plot(k_range, sil_mean, 'o-', 'LineWidth',1.5, 'MarkerSize',8);
xlabel('聚类数K'); ylabel('轮廓系数均值');
title('轮廓系数法：基因表达数据K值选择');
grid on;
% 找最大值对应的K
[max_sil, max_idx] = max(sil_mean);
best_k_sil = k_range(max_idx);
hold on;
plot(best_k_sil, max_sil, 'rs', 'MarkerSize',10);
text(best_k_sil+0.2, max_sil, sprintf('最优K=%d（均值=%.2f）', best_k_sil, max_sil), 'FontSize',12);
hold off;

结果解读：示例中K=3时轮廓系数均值最大（≈0.7），验证了肘部法则的结论，且数学上更严谨。

方法3：间隙统计量法（Gap Statistic）—— 高维数据首选，消除主观判断

原理：对比“真实数据的SSE”和“随机均匀分布数据的SSE”，计算“间隙值（Gap）”，Gap最大的K为最优（考虑数据的绝对分布，而非相对变化）。 适配场景：高维基因表达数据（维度灾难导致肘部法则不明显），MATLAB无内置函数，需自定义实现。

% ===================== 间隙统计量法实现（自定义函数） =====================
function [best_k_gap, gap, sdk] = gap_statistic(X, k_max, B)
    % X：数据矩阵（行=样本，列=特征）
    % k_max：最大测试K值
    % B：随机数据生成次数（通常取10~20）
    [n, p] = size(X);
    gap = zeros(k_max, 1);
    sdk = zeros(k_max, 1);
    % 计算真实数据的SSE
    sse_real = zeros(k_max, 1);
    for k = 1:k_max
        [~, ~, sumd] = kmeans(X, k, 'Start','plus', 'Replicates',5);
        sse_real(k) = sum(sumd);
    end
    % 生成B组随机数据，计算SSE
    sse_rand = zeros(k_max, B);
    for b = 1:B
        % 生成均匀随机数据（值域与原数据一致）
        X_rand = rand(n,p) .* (max(X,[],1)-min(X,[],1)) + min(X,[],1);
        for k = 1:k_max
            [~, ~, sumd] = kmeans(X_rand, k, 'Start','plus', 'Replicates',5);
            sse_rand(k,b) = sum(sumd);
        end
    end
    % 计算Gap和标准差
    for k = 1:k_max
        log_sse_rand = log(sse_rand(k,:));
        gap(k) = mean(log_sse_rand) - log(sse_real(k));
        sdk(k) = std(log_sse_rand) * sqrt(1 + 1/B);
    end
    % 找最优K（第一个Gap(k) ≥ Gap(k+1) - sdk(k+1)）
    best_k_gap = 1;
    for k = 1:k_max-1
        if gap(k) >= gap(k+1) - sdk(k+1)
            best_k_gap = k;
            break;
        end
    end
    if best_k_gap == 1 && gap(k_max) > gap(1)
        best_k_gap = k_max;
    end
end

% ===================== 调用间隙统计量函数 =====================
k_max = 10;
B = 10; % 随机数据生成次数（越大越准，耗时越长）
[best_k_gap, gap, sdk] = gap_statistic(gene_expr, k_max, B);

% 可视化间隙统计量
figure('Color','w');
plot(1:k_max, gap, 'o-', 'LineWidth',1.5, 'MarkerSize',8);
errorbar(1:k_max, gap, sdk, 'Color','k', 'Alpha',0.5);
xlabel('聚类数K'); ylabel('间隙值（Gap）');
title('间隙统计量法：基因表达数据K值选择');
grid on;
hold on;
plot(best_k_gap, gap(best_k_gap), 'rs', 'MarkerSize',10);
text(best_k_gap+0.2, gap(best_k_gap), sprintf('最优K=%d', best_k_gap), 'FontSize',12);
hold off;

结果解读：间隙统计量法对高维基因数据更严谨，示例中最优K仍为3，验证了前两种方法的结论。

方法4：基于生物学先验/ARI（有标签时）—— 最贴合实验逻辑

原理：

  • 有生物学先验：直接指定K（如已知基因有3个功能模块→K=3；上调+下调→K=2）；
  • 有实验分组标签（如Control/Treatment）：计算不同K的调整兰德指数（ARI）（衡量聚类结果与真实标签的一致性，取值[0,1]），ARI最大值对应最优K。
% ===================== 基于ARI的K值选择 =====================
% 1. 构造真实生物学标签（模拟：3个模块）
true_label = [ones(30,1); 2*ones(40,1); 3*ones(30,1)]; % 真实模块标签
k_range = 1:10;
ari = zeros(length(k_range), 1);

% 2. 循环计算每个K的ARI
for i = 1:length(k_range)
    k = k_range(i);
    if k == 1
        ari(i) = 0;
        continue;
    end
    [cluster_idx, ~] = kmeans(gene_expr, k, 'Start','plus', 'Replicates',5);
    ari(i) = adjustedRandIndex(cluster_idx, true_label); % MATLAB内置函数
end

% 3. 可视化ARI曲线
figure('Color','w');
plot(k_range, ari, 'o-', 'LineWidth',1.5, 'MarkerSize',8);
xlabel('聚类数K'); ylabel('调整兰德指数（ARI）');
title('基于ARI的K值选择（贴合生物学标签）');
grid on;
% 找ARI最大值对应的K
[max_ari, max_idx] = max(ari);
best_k_ari = k_range(max_idx);
hold on;
plot(best_k_ari, max_ari, 'rs', 'MarkerSize',10);
text(best_k_ari+0.2, max_ari, sprintf('最优K=%d（ARI=%.2f）', best_k_ari, max_ari), 'FontSize',12);
hold off;

结果解读：示例中K=3时ARI=0.92（接近1），说明聚类结果与真实生物学模块高度匹配，是最具解释性的最优K。

三、实战案例（基因表达数据完整流程）

% ===================== 完整流程：综合确定最优K =====================
% 1. 数据预处理（筛选差异基因，避免高维噪声）
expr_raw = randn(n_samples, 500); % 原始基因表达数据（样本×500基因）
% 过滤低表达基因
expr_filtered = expr_raw(:, mean(expr_raw>0.1,1)>=0.8);
% 标准化
expr_norm = zscore(expr_filtered);
% 筛选差异基因（FC>1.5）
group = categorical([repmat('Control',20,1); repmat('Treatment',40,1)]);
fc = mean(expr_norm(group=='Treatment',:)) - mean(expr_norm(group=='Control',:));
diff_genes = abs(fc)>1.5;
gene_expr = expr_norm(:, diff_genes)'; % 差异基因×样本

% 2. 多方法验证K值
% 肘部法则：K=3
% 轮廓系数：K=3
% 间隙统计量：K=3
% ARI（真实标签）：K=3

% 3. 最终聚类（K=3）
[cluster_idx, centroid] = kmeans(gene_expr, 3, 'Start','plus', 'Replicates',5);

% 4. 验证聚类效果
sil = silhouette(gene_expr, cluster_idx);
ari = adjustedRandIndex(cluster_idx, true_label);
fprintf('最终K=3，轮廓系数均值=%.2f，ARI=%.2f\n', mean(sil), ari);

% 5. 可视化聚类结果（PCA降维后）
[coeff, score] = pca(gene_expr);
score_2d = score(:,1:2);
figure('Color','w');
gscatter(score_2d(:,1), score_2d(:,2), cluster_idx, viridis(3), 'o', 8);
xlabel('PCA维度1'); ylabel('PCA维度2');
title('K=3时基因聚类结果（PCA降维可视化）');
grid on;

四、关键注意事项（适配基因表达数据）

  1. 高维预处理优先： 聚类前必须筛选差异基因（FC>1.5+FDR<0.05）或PCA降维（保留前2~5主成分），避免“维度灾难”导致K值选择失真。

  2. 方法优先级：

    • 有生物学先验（如已知模块数）→ 直接指定K；
    • 有实验标签→ ARI法；
    • 无先验→ 轮廓系数法（主）+ 间隙统计量法（验证）+ 肘部法则（初筛）。

  3. 避免过聚类： 若数学方法推荐K=5，但生物学仅存在3个功能模块，优先选K=3（聚类的核心是解释生物学意义，而非纯数学最优）。

  4. 耗时平衡： 间隙统计量的B值（随机数据次数）建议取10（基因数据），过大（如50）会显著增加耗时，且增益有限。

五、方法选型总结

方法优势劣势适配场景
肘部法则简单、快速、易实现主观、高维数据不明显初筛K值范围、小维度基因数据
轮廓系数法客观、兼顾簇内/簇间对异常值敏感无先验、高维基因数据（核心方法）
间隙统计量法最严谨、适配高维数据耗时、需自定义函数高维基因数据（验证最优K）
生物学先验/ARI贴合实验逻辑、结果易解释依赖先验/标签有功能模块知识/实验分组的基因数据

通过以上方法，可确保K值既满足数学聚类效果，又贴合基因表达的生物学模块特征，让聚类结果具备统计意义和实验解释性。

avatar
文章
阅读
粉丝