简介：这份资源面向生物信息学入门学习者与R语言数据分析实践者，围绕TCGA-BRCA乳腺癌基因表达数据展开聚类分析教学案例。内容涵盖利用层次聚类按基因表达水平对病人样本分型，采用average距离度量，并通过热图等图示呈现聚类效果；同时实现PCA降维，讨论主成分数目选择依据，再基于降维特征重新聚类并与原始结果对比，最后借助临床数据中的ER_Status_nature2012字段评估聚类是否符合预期。压缩包共22个文件，约10.91MB，包含8个PNG与8个PDF结果图、2个markdown说明、2个txt数据文件、1个license及1个R源码脚本，图表与代码配套便于复现。已有687人学习下载，适合希望掌握聚类与降维完整流程、对照图示理解分析思路的读者参考。

1. 生物信息学概论里的聚类分析：为什么 TCGA-BRCA 是练手的第一站

如果你刚接触生物信息学概论，又恰好拿到一份 TCGA-BRCA 的表达矩阵，最想干的事大概就是「把样本分个组看看」。聚类分析就是干这个的：不给你标签，让数据自己说话，把表达模式相近的样本或基因归到一起。TCGA-BRCA 是乳腺癌队列，样本量大、临床注释相对完整，表达谱维度高，天然适合拿来练无监督聚类。它解决的核心问题是：在一堆看不出规律的表达值里，找出分子层面的亚型结构，或者筛出一批协同表达的基因模块。适合谁？适合已经会读表达矩阵、会跑 Python 或 R、但对「怎么选特征、怎么定簇数、怎么判断结果靠不靠谱」还没底的人。这一篇不讲空泛概念，直接按能复现的路径走一遍，把参数和坑都摊开。

2. 拿到表达矩阵先别急着聚类：数据清洗与特征筛选的四个动作

聚类对输入极其敏感。同一份 TCGA-BRCA 数据，你直接丢进 KMeans 和先做一轮清洗再聚类，结果可能完全是两回事。这一步的目标是把「噪声大、方差低、缺失多」的维度砍掉，让距离度量真正反映生物学差异。

2.1 表达矩阵的读入与基因 ID 对齐

TCGA-BRCA 的表达数据常见形式是每行一个基因、每列一个样本，基因 ID 可能是 Ensembl 也可能是基因符号。第一步是确认索引唯一、样本列名可追溯。下面用 pandas 读入并做基本检查。

import pandas as pd import numpy as np # 读入表达矩阵，第一列是基因 ID expr = pd.read_csv("tcga_brca_expr.csv", index_col=0) # 基本体检：维度、缺失、重复索引 print("矩阵维度:", expr.shape) print("缺失值总数:", expr.isna().sum().sum()) print("重复基因 ID 数:", expr.index.duplicated().sum()) # 基因 ID 去重：同一基因多个探针时取表达均值 expr = expr.groupby(expr.index).mean() # 样本列名去掉多余空白，避免后续匹配临床数据时对不上 expr.columns = expr.columns.str.strip()

逻辑说明：groupby(index).mean()处理的是同一基因对应多行的情况，取均值是最稳妥的默认策略。参数上，index_col=0要求你的 CSV 第一列确实是基因标识；如果读进来发现列名是Unnamed: 0，说明导出时索引没处理好，需要回退检查。缺失值这里只做统计，不急着填，因为下一步的筛选会把高缺失基因直接剔除。

2.2 低表达基因过滤与对数变换

测序数据里大量基因在所有样本中表达接近零，这些基因不携带区分信息，还会拉高维度。常见做法是保留在至少一定比例样本中表达量超过阈值的基因。同时，原始 count 或 FPKM 的分布右偏严重，聚类前一般做 log2 变换压缩动态范围。

# 过滤：至少在 50% 样本中表达量 > 1 min_samples = int(expr.shape[1] * 0.5) mask = (expr > 1).sum(axis=1) >= min_samples expr_filt = expr.loc[mask] # log2 变换，加 1 防止 log(0) expr_log = np.log2(expr_filt + 1) print("过滤后基因数:", expr_log.shape[0])

逻辑说明：min_samples取样本数的一半是经验值，样本量小的时候可以放宽到 30%。阈值> 1针对的是 FPKM/TPM 类数据；如果你手上是原始 count，这个阈值要相应调大，或者改用 CPM 归一化后再判断。log2 变换是聚类前的标准动作，但注意：如果数据已经是对数尺度，再取一次 log 会把差异压平，动手前先看数值范围。

2.3 高变基因筛选：聚类的信号主要来自这里

全基因集聚类计算量大且噪声多，通常先选高变基因（HVG）。思路是计算每个基因的方差或变异系数，取排名靠前的一部分。这一步直接决定聚类结果有没有生物学解释力。

from sklearn.feature_selection import VarianceThreshold # 按方差选 top 2000 高变基因 variances = expr_log.var(axis=1) top_genes = variances.sort_values(ascending=False).head(2000).index expr_hvg = expr_log.loc[top_genes] print("高变基因矩阵:", expr_hvg.shape)

逻辑说明：head(2000)是常用起点，样本量几百到上千时，1000 到 3000 之间都合理。选太多会把噪声带进来，选太少可能丢掉关键亚型信号。判断标准可以看后续聚类轮廓系数，或者直接看高变基因里有没有已知的乳腺癌标志物。注意方差计算前要确认数据已经过对数变换，否则高表达基因会主导排名。

2.4 样本与基因的双向标准化取舍

聚类前要不要对样本做 z-score，取决于你想让什么主导距离。对基因维度标准化（每个基因减均值除标准差）能消除表达量绝对高低的影响，让所有基因等权参与。对样本维度标准化则会改变样本间距离的含义，一般不做。

from sklearn.preprocessing import StandardScaler # 转置成样本 x 基因，对基因维度标准化 X = expr_hvg.T scaler = StandardScaler() X_scaled = scaler.fit_transform(X) print("标准化后矩阵:", X_scaled.shape)

逻辑说明：expr_hvg.T把矩阵转成样本为行、基因为列，这是 sklearn 的输入格式要求。StandardScaler默认按列（即每个基因）做 z-score，正好符合需求。参数上不需要额外调整，但要注意：标准化后数据不再有生物学上的绝对表达含义，解释结果时看的是相对模式。如果后续要用相关性距离，标准化这一步可以省略，因为相关系数本身对尺度不敏感。

3. 聚类算法怎么选：KMeans、层次聚类与共识聚类的实战对比

数据准备好之后，下一个问题是算法。TCGA-BRCA 这类表达谱聚类，常见选择是 KMeans、层次聚类和共识聚类。它们不是随便挑一个就行，样本量、簇形状假设、结果稳定性都不一样。

3.1 KMeans 的最小可跑通流程与簇数选择

KMeans 快、易解释，适合先跑一个基线。它的假设是簇近似球形、大小相近，表达谱数据不一定满足，但作为起点足够。

from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score inertias = [] sil_scores = [] K_range = range(2, 9) for k in K_range: km = KMeans(n_clusters=k, n_init=10, random_state=42) labels = km.fit_predict(X_scaled) inertias.append(km.inertia_) sil_scores.append(silhouette_score(X_scaled, labels)) # 输出每个 k 的惯性和轮廓系数 for k, inertia, sil in zip(K_range, inertias, sil_scores): print(f"k={k}, inertia={inertia:.1f}, silhouette={sil:.3f}")

逻辑说明：n_init=10表示用不同初始化跑 10 次取最优，避免陷入局部最优，这个参数在新版 sklearn 里默认值有变化，显式写出来更稳。random_state固定后结果可复现。簇数选择看两个指标：惯性下降的拐点和轮廓系数峰值。实际中两者不一定一致，轮廓系数更贴近「分得开不开」，优先参考它。注意轮廓系数在样本量很大时计算开销不小，可以抽样估算。

3.2 层次聚类：距离度量和连接方式比算法本身更重要

层次聚类不预设簇数，输出树状图，适合观察样本的层级关系。它的关键参数是距离度量和连接方式，这两个选错，树形结构会完全失真。

from scipy.cluster.hierarchy import linkage, fcluster, dendrogram from scipy.spatial.distance import pdist # 计算样本间距离，用相关性距离更贴合表达谱 dist_matrix = pdist(X_scaled, metric="correlation") # Ward 连接 Z = linkage(dist_matrix, method="ward") # 按距离阈值切树，得到簇标签 labels_hc = fcluster(Z, t=4, criterion="maxclust") print("层次聚类簇分布:", np.bincount(labels_hc))

逻辑说明：metric="correlation"衡量的是表达模式形状的相似性，比欧氏距离更适合表达谱。method="ward"要求距离是欧氏距离才严格成立，用相关性距离时更稳妥的选择是average或complete。这里写 ward 是为了演示流程，实际用相关性距离时建议换成average。fcluster的criterion="maxclust"表示按最大簇数切，t=4就是切出 4 个簇。树状图可以用dendrogram(Z)画出来目视检查，但样本超过几百时叶子会糊成一片，需要截断显示。

3.3 共识聚类：为什么 TCGA 亚型分析更信它

共识聚类（consensus clustering）的思路是反复重采样，看哪些样本稳定地分在一起。它比单次 KMeans 更抗噪，在肿瘤亚型分析里被广泛采用。核心参数是重采样次数和每次抽样的特征比例。

from sklearn.cluster import KMeans from sklearn.utils import resample def consensus_cluster(X, n_iter=100, k=4, sample_ratio=0.8): n_samples = X.shape[0] co_matrix = np.zeros((n_samples, n_samples)) for i in range(n_iter): # 对样本重采样 idx = resample(range(n_samples), n_samples=int(n_samples * sample_ratio), random_state=i) km = KMeans(n_clusters=k, n_init=5, random_state=i) labels = km.fit_predict(X[idx]) # 同簇样本对计数 for a in range(len(idx)): for b in range(a + 1, len(idx)): if labels[a] == labels[b]: co_matrix[idx[a], idx[b]] += 1 co_matrix[idx[b], idx[a]] += 1 # 归一化 co_matrix = co_matrix / n_iter return co_matrix co_mat = consensus_cluster(X_scaled, n_iter=100, k=4) print("共识矩阵形状:", co_mat.shape)

逻辑说明：n_iter=100是常见起点，追求稳定可以加到 500 甚至 1000，代价是时间线性增长。sample_ratio=0.8表示每次抽 80% 样本，比例太低会导致共识矩阵稀疏，太高则重采样多样性不足。共识矩阵本身可以再做一次层次聚类得到最终标签，也可以直接看矩阵的热图判断簇边界清不清晰。这段代码是双层循环，样本量大时用向量化或并行会快很多，但作为理解原理的最小实现够用。

3.4 三种方法的结果怎么对照

跑完三种方法，不要只看某一个的标签就下结论。常见做法是算两两之间的调整兰德指数（ARI），看一致性有多高。

from sklearn.metrics import adjusted_rand_score # 假设 labels_km 是 KMeans 结果，labels_hc 是层次聚类结果 print("KMeans vs 层次聚类 ARI:", adjusted_rand_score(labels_km, labels_hc))

逻辑说明：ARI 取值在 -1 到 1 之间，越接近 1 表示两种划分越一致。如果三种方法结果差异很大，说明数据结构本身不清晰，或者特征筛选、标准化环节有问题，需要回头检查，而不是硬选一个结果往下走。

4. 聚类结果怎么验证与解释：从轮廓系数到临床注释

聚类是无监督的，没有标签告诉你对错。验证和解释是让结果站得住脚的关键，也是新手最容易糊弄过去的一步。

4.1 内部指标：轮廓系数、Calinski-Harabasz 与 Davies-Bouldin

内部指标只依赖数据和标签，用来评估簇的紧致度和分离度。三个常用指标各有侧重，建议一起看。

from sklearn.metrics import calinski_harabasz_score, davies_bouldin_score print("CH 指数:", calinski_harabasz_score(X_scaled, labels_km)) print("DB 指数:", davies_bouldin_score(X_scaled, labels_km))

逻辑说明：CH 指数和 DB 指数计算快，适合在 k 值扫描时一起输出。注意这些指标都基于距离，如果换了距离度量，数值不能直接跨方法比较。轮廓系数对簇形状敏感，KMeans 的球形假设下会偏高，换成层次聚类结果后可能下降，这不代表层次聚类更差。

4.2 外部验证：用临床注释做交叉表

TCGA-BRCA 有临床数据，比如受体状态、分期、生存信息。把聚类标签和这些变量做交叉表，看簇是否对应已知的生物学分组，是最有说服力的验证。

import pandas as pd from scipy.stats import chi2_contingency # 假设 clinical 是临床数据框，含 sample_id 和 receptor_status clinical = pd.read_csv("tcga_brca_clinical.csv", index_col=0) merged = clinical.loc[X.index].copy() merged["cluster"] = labels_km # 交叉表与卡方检验 ct = pd.crosstab(merged["cluster"], merged["receptor_status"]) print(ct) chi2, p, dof, expected = chi2_contingency(ct) print(f"卡方 p 值: {p:.4f}")

逻辑说明：clinical.loc[X.index]保证临床数据和聚类样本顺序对齐，这一步不对齐后面全错。卡方检验的 p 值小于 0.05 说明聚类标签和受体状态不是独立分布的，即簇有生物学意义。但要注意样本量小时卡方检验的期望频数可能不满足假设，需要看expected里有没有小于 5 的格子。

4.3 生存分析：聚类标签有没有预后区分度

如果临床数据里有生存时间和事件状态，可以用 Kaplan-Meier 曲线看不同簇的生存差异。这是把聚类结果往临床价值上靠的关键一步。

from lifelines import KaplanMeierFitter from lifelines.statistics import logrank_test import matplotlib.pyplot as plt kmf = KaplanMeierFitter() for c in sorted(merged["cluster"].unique()): mask = merged["cluster"] == c kmf.fit(merged.loc[mask, "survival_time"], merged.loc[mask, "event"], label=f"Cluster {c}") kmf.plot_survival_function() plt.savefig("km_curve.png")

逻辑说明：survival_time和event是生存分析的标准输入，前者是随访时间，后者是事件是否发生。logrank_test可以两两比较或整体比较，p 值显著说明簇间生存有差异。注意生存分析对样本量和随访完整性敏感，如果某个簇样本很少，曲线会很不稳定，解释时要谨慎。

4.4 基因模块解释：簇间差异表达基因

聚类分的是样本，但解释往往要落到基因上。对每个簇做差异表达分析，看哪些基因在簇间变化最大，再做过表达富集，能给出机制层面的线索。

from scipy.stats import ttest_ind # 以簇 0 vs 其余为例 group0 = X_scaled[labels_km == 0] group_rest = X_scaled[labels_km != 0] t_stat, p_val = ttest_ind(group0, group_rest, axis=0) # 取 p 值最小的前 50 个基因 top_idx = np.argsort(p_val)[:50] top_gene_names = expr_hvg.index[top_idx] print("簇 0 差异基因:", list(top_gene_names[:10]))

逻辑说明：这里用的是 t 检验，对表达谱数据来说是最简单的差异分析，实际中更常用 limma 或 DESeq2 的框架。axis=0表示沿样本方向比较。p 值没有做多重检验校正，只是演示筛选逻辑，正式分析要用 Benjamini-Hochberg 等方法控制 FDR。取出的基因列表可以丢进富集分析工具，看有没有富集到已知通路。

5. 避坑与排查：TCGA-BRCA 聚类里最容易翻车的五个地方

这一章是血泪经验集中区。下面五条都是实际跑数据时反复遇到的，按「现象 → 原因 → 解决」写清楚。

5.1 聚类结果每次跑都不一样

现象：同样的数据和代码，换一次运行，样本分组变了。原因：KMeans 的初始化是随机的，没有固定random_state；或者共识聚类的重采样种子没控制。解决：所有涉及随机的步骤都显式设random_state，KMeans 的n_init设大一点（10 到 50），共识聚类的每次迭代用固定种子序列。如果固定种子后结果仍然不稳定，说明数据结构本身模糊，需要回到特征筛选环节。

5.2 所有样本被分到一个簇

现象：不管 k 设多少，某个簇吞掉了 90% 以上的样本。原因：数据没有标准化，少数高表达基因主导了距离计算；或者高变基因选得不对，选进来的都是低方差基因。解决：确认标准化步骤在特征筛选之后执行，检查高变基因的方差排名是否合理。可以画一个样本距离分布图，如果大多数样本对的距离都很小，说明特征没有区分度。

5.3 轮廓系数很高但临床交叉表毫无关联

现象：内部指标看起来不错，但聚类标签和受体状态、分期完全不相关。原因：聚类捕捉到的是技术噪声，比如批次效应、测序深度差异，而不是生物学信号。解决：检查数据有没有做批次校正，样本来源是否混杂。可以先用 PCA 看前几个主成分和批次变量的相关性，如果 PC1 和批次高度相关，聚类前必须做校正。

5.4 层次聚类树状图完全糊在一起

现象：树状图底部密密麻麻，看不出任何簇结构。原因：距离度量或连接方式选错，比如用欧氏距离处理未标准化的表达谱；或者样本量太大，树状图本身就不适合目视。解决：换相关性距离，连接方式用 average；样本超过 200 时不要硬看全图，用truncate_mode截断显示，或者直接看共识矩阵的热图。

5.5 差异表达基因列表里全是同一类基因

现象：簇间差异基因 top50 里大部分是同一家族或同一通路。原因：这些基因本身高度共表达，聚类时作为一个模块被选进来，差异分析时又一起显著。解决：这不是错误，但解释时要意识到基因间的共线性。可以做基因聚类先把共表达模块合并，或者用富集分析看整体通路，而不是盯着单个基因。

6. 把聚类分析做成可复用的流程：参数固化与结果归档

跑通一次不难，难的是下次换一份数据还能稳定复现。这一章讲怎么把上面的步骤收成一个可复用的流程，以及一个我常用的结果检查习惯。

6.1 用配置文件固化关键参数

把样本过滤比例、高变基因数、簇数范围、随机种子这些写进一个 YAML 或 JSON，代码只读配置，不硬编码。这样换数据集时只改配置，不动逻辑。

import yaml config = { "min_sample_ratio": 0.5, "min_expr_threshold": 1, "n_hvg": 2000, "k_range": [2, 3, 4, 5, 6, 7, 8], "random_state": 42, "consensus_n_iter": 100, "consensus_sample_ratio": 0.8 } with open("cluster_config.yaml", "w") as f: yaml.dump(config, f)

逻辑说明：min_sample_ratio和min_expr_threshold控制基因过滤，n_hvg控制特征数，k_range控制簇数扫描范围。把这些参数外置的好处是，当你换一个队列时，只需要调整配置值，不用去代码里翻哪一行改了阈值。random_state统一管理，保证所有随机步骤可复现。

6.2 结果归档：标签、指标、图一个都不能少

每次运行输出一个带时间戳的目录，里面存样本标签、各 k 的指标表、共识矩阵、树状图、KM 曲线。这样过一个月回头看，还能知道当时用的什么参数、结果长什么样。

import os import datetime run_id = datetime.datetime.now().strftime("%Y%m%d_%H%M%S") out_dir = f"results/{run_id}" os.makedirs(out_dir, exist_ok=True) # 保存标签 pd.DataFrame({"sample": X.index, "cluster": labels_km}).to_csv(f"{out_dir}/labels.csv", index=False) # 保存指标 pd.DataFrame({"k": list(K_range), "inertia": inertias, "silhouette": sil_scores}).to_csv(f"{out_dir}/metrics.csv", index=False) print("结果已归档到:", out_dir)

逻辑说明：run_id用时间戳保证每次运行不覆盖。标签文件里保留样本 ID 和簇编号，方便后续和临床数据合并。指标文件记录了每个 k 的表现，下次想换簇数时可以直接对比。图文件用savefig存到同一目录，不要只留在内存里。

6.3 一个我常用的结果检查习惯

每次跑完聚类，我会先做一件事：把共识矩阵的热图、PCA 散点图（按簇着色）、KM 曲线三张图并排看一遍。如果共识矩阵的块状结构清晰、PCA 里簇有分离趋势、KM 曲线有区分度，这三个信号同时出现，我才认为这次聚类值得往下做。如果只有其中一个成立，我会回头检查特征筛选和标准化，而不是急着解释结果。这个习惯帮我省了很多在噪声结果上浪费时间的情况。希望帮到你。

本文还有配套的精品资源，点击获取