简介生物信息学概论TCGA-BRCA数据聚类分析任务的完整解决方案包面向需要掌握R语言聚类分析、PCA降维方法的高校学生与生信初学者。压缩包内含R源码、基因表达矩阵与临床注释数据以及层次聚类热图、碎石图、PCA聚类结果等全套输出图表可帮助读者快速复现「按基因表达水平对乳腺癌病人聚类并结合ER_Status_nature2012评估聚类效果」的完整流程。资源共22个文件以png、pdf图表、md说明文档、txt数据文件及一个r脚本为主包体10.91MB目录结构清晰README中对数据字段、代码逻辑与聚类参数做了简要说明便于对照学习。已有687人学习下载适合作为课程作业参考或自学实操案例尤其能直观理解不同距离/聚类参数对结果的影响以及主成分数目的选择依据。1. 聚类分析TCGA-BRCA数据一个压缩包背后是完整的生信分析链路你拿到的这份「生物信息学概论——聚类分析TCGA-BRCA数据.zip」大概率不是一份可以直接运行的成品而是一套教学性质的数据包里面装着TCGA乳腺癌BRCA的表达谱数据、临床注释以及一段聚类分析的示例脚本。它的价值不在“解压即用”而在于帮你走通一条从原始数据到分子分型的完整链路——这也是生物信息学入门最典型的实战场景。很多人以为聚类分析就是调个KMeans跑两行代码真正动手才发现数据清洗、特征筛选、k值确定、结果验证每一步都能让结果看起来合理却完全不可复现。这篇文章不会假装我读过你那个zip里的源码我只按照TCGA-BRCA聚类分析最常规的落地路径把每个环节的命令、参数和坑位讲清楚。适合谁看刚起步的生信研究生、湿实验想转干实验的医生、以及所有被“聚类分析”四个字劝退过的人。2. 拿到TCGA-BRCA数据之后先做这三步数据清洗2.1 表达谱矩阵的读取与样本名转换TCGA的BRCA项目里有超过一千个肿瘤样本和数量不等的癌旁样本表达谱文件通常以基因ID为行名、以样本条码barcode为列名。你从GDC Data Portal或GDC客户端下载到的表达谱常见的是HTSeq-Counts和HTSeq-FPKM两种格式。压缩包里的数据大概率已经做了初步整理但第一步永远是确认行名和列名的准确含义。我在本地一般先把表达矩阵读成DataFrame前几列看一眼再说。这里的关键动作是把ENSEMBL ID转成基因名同时把样本条码裁剪成12位的TCGA患者ID。import pandas as pd import re # 读取表达谱第一列为基因ID expr pd.read_csv(BRCA_expression.tsv, sep\t, index_col0) print(expr.iloc[:5, :5]) print(expr.shape) # ENSEMBL ID转基因名常见做法是用 mygene 在线注释 # 如果压缩包里已经给了基因名这一步可以跳过 from mygene import MyGeneInfo mg MyGeneInfo() ids expr.index.tolist() annot mg.querymany(ids, scopesensembl.gene, fieldssymbol, specieshuman) symbol_map {item[query]: item.get(symbol, item[query]) for item in annot} expr.index [symbol_map.get(i, i) for i in expr.index] # 条件样本名TCGA-XX-XXXX-01 → 取前12位去掉样本类型编号 expr.columns [c[:12] for c in expr.columns] # 保留肿瘤样本第14、15位为01癌旁样本11可根据需要剔除或保留 tumor_mask [c[13:15] 01 for c in expr.columns] expr_tumor expr.loc[:, tumor_mask] print(肿瘤样本数:, expr_tumor.shape[1])这段代码里querymany返回的列表顺序不一定和输入一致所以用字典按query字段回填。样本条码的规则是TCGA-XX-XXXX-XX第14、15位01是肿瘤组织11是癌旁聚类分析通常只在肿瘤样本上做否则样本分组的第一个维度永远是“癌与癌旁”根本看不到肿瘤内部的异质性。很多时候压缩包里给的表达矩阵列名已经被处理成了纯字母数字那expr.columns提取出来的片段长度和第13、14位的取法可能对不上。先用print(expr.columns[:10])肉眼确认一遍再去套正则。这个步骤花不了两分钟能省掉后面大半的坑。2.2 基因过滤与缺失值处理表达谱矩阵里大量基因在所有样本中的表达量都是零这些基因不携带任何分组信息却会稀释后续特征筛选的信号。我习惯先把全零基因删掉再过滤掉缺失比例过高的基因。import numpy as np # 删除全零基因 expr_filtered expr_tumor.loc[(expr_tumor ! 0).sum(axis1) 0, :] # 删除缺失比例超过20%的基因TCGA表达谱一般不做缺失但以防万一 miss_ratio expr_filtered.isna().sum(axis1) / expr_filtered.shape[1] expr_filtered expr_filtered.loc[miss_ratio 0.2, :] # 对剩余缺失值做中位数填充 expr_filtered expr_filtered.apply(lambda row: row.fillna(row.median()), axis1) print(过滤后基因数:, expr_filtered.shape[0])全零基因过滤是硬性条件缺失率过滤则看你拿到的数据类型——TCGA的官方表达矩阵基本没有缺失但一些二次处理过的文件里可能存在异常列。中位数填充是一种保守策略但它只适用于少量缺失如果某基因缺失超过一半填充多少都没意义不如直接删。这里特别注意不要用均值填充表达谱数据存在离群样本均值会被拉偏中位数稳得多。2.3 归一化为什么不能拿原始Counts直接聚类这是整个数据清洗环节里最“玄学”也最关键的一步。HTSeq-Counts是原始read数数值范围从几千到几千万跟基因长度、测序深度都强相关。你拿原始Counts做聚类得到的第一个主成分几乎永远是测序深度而不是生物学差异。常见做法是取log2(x1)把数据压缩到相对平稳的尺度。如果压缩包里给的是FPKM或TPM建议也做一次log2变换再进入下游。# log2(x1)变换 expr_log np.log2(expr_filtered 1) print(expr_log.iloc[:5, :5]) # 也可以做z-score标准化但这一步一般放到特征筛选之后 # 如果一定要先看全局分布可以用StandardScaler from sklearn.preprocessing import StandardScaler scaler StandardScaler() expr_z scaler.fit_transform(expr_log.T).T # 按基因维度标准化log2(x1)里的1是为了避免零的对数无定义这个常数不能省。至于为什么不是log2(x0.5)或log2(x2)实际影响很小业界通用做法就是1。先log2再z-score是标准流水线log2解决偏态分布z-score让每个基因在聚类时有相同的权重避免高表达基因主导距离计算。这两步的顺序不能反先标准化再取log会把零值残留的离散度放大效果差很多。3. 选特征而不是选算法聚类前的高变基因筛选3.1 为什么全基因组两万个基因直接聚类会翻车TCGA-BRCA表达谱过滤后通常还剩一万多个基因。把这一万多个维度直接丢进KMeans计算量是一方面更致命的是信号被稀释只有少数基因与乳腺癌分子分型相关其余大部分基因的表达波动来自技术噪声和个体差异。你得到的结果可能每次都稳定地分成三组但分组边界跟生物学亚型毫无关系。这就是高变基因筛选的意义。它的逻辑很朴素在所有样本间波动最大的基因最可能携带分组信息——因为如果只是个体间随机差异不会让几百个基因同时出现大幅波动。我一般取变异系数CV排名前2000到5000的基因这个区间既保留了足够信号又把噪声控制在一定范围内。3.2 用Python计算变异系数并筛选Top基因# 变异系数 标准差 / 均值只对log2之后的数据计算 mean_val expr_log.mean(axis1) std_val expr_log.std(axis1) cv std_val / mean_val # 按CV降序取前3000 n_top 3000 top_genes cv.sort_values(ascendingFalse).head(n_top).index expr_top expr_log.loc[top_genes, :] print(筛选后维度:, expr_top.shape) # 保存精选矩阵供聚类使用 expr_top.to_csv(BRCA_top3000_log2.tsv, sep\t)变异系数在这里比单纯方差更合适因为高表达基因的绝对方差天然就大但它们的倍性变化往往不如中低表达基因有区分度。CV做了归一化处理把表达水平的影响消掉了。2000到5000这个区间没有硬性规定TCGA泛癌聚类文章常用前25%或前50%的高变基因换算到两万基因维度就是5000到10000。但特征是越少越好解读我通常先用3000试跑如果聚类轮廓系数不理想再调整。有的教程会建议用方差排名那也不是不行只是结果会更偏向高表达基因。两种方式都可以跑关键是要明白自己选的是什么。如果你拿到的zip里给的特征基因列表只有几百个那大概率是作者提前做过一轮差异表达或相关性筛选了直接沿用就好。3.3 样本标准化让每个基因在同一尺度上说话特征筛选之后还必须做一次按基因的标准化。上一个环节我们只做了log2没做跨基因归一化因为高变基因筛选需要保留基因间的表达水平差异来计算CV。但聚类时用的是距离矩阵如果一个基因的均值是10、标准差是2另一个基因的均值是3、标准差是0.1前者会天然主导欧氏距离。这显然不是你想要的。# 按基因行做z-score标准化 from sklearn.preprocessing import StandardScaler scaler StandardScaler() expr_std pd.DataFrame( scaler.fit_transform(expr_top.T).T, indexexpr_top.index, columnsexpr_top.columns ) print(expr_std.iloc[:5, :5]) # 保存最终输入矩阵 expr_std.to_csv(BRCA_cluster_input.tsv, sep\t)这里fit_transform(expr_top.T).T两次转置是因为StandardScaler默认按列标准化而我们的基因为行、样本为列转了两次正好实现逐基因标准化。标准化后每个基因的均值被拉到0、方差被拉到1聚类距离计算里每个基因的贡献趋于平等。这一步做完你的输入矩阵才真正达到了“可以直接丢给聚类算法”的标准。很多人前面的预处理都对偏偏忘了标准化结果层次聚类的树状图乍一看合理换两个样本又完全变形就是这个原因。4. 三种聚类跑出分子亚型KMeans、层次聚类、一致性聚类4.1 KMeansk值怎么定轮廓系数和肘部法则KMeans是聚类分析最常用的入门算法但它的k值不会自己告诉你。回答“分几类”这个问题统计指标和生物学背景要结合看。TCGA乳腺癌的经典分型PAM50有五个亚型Luminal A、Luminal B、HER2富集、Basal样、Normal样所以很多做BRCA聚类的人会直接设k5这没问题但你最好先跑一遍肘部法则和轮廓系数看看数据支不支持。from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score import matplotlib.pyplot as plt X expr_std.T # 样本为行基因为列 # 肘部法则 轮廓系数 inertia_list [] sil_list [] k_range range(2, 9) for k in k_range: km KMeans(n_clustersk, random_state42, n_init10) labels km.fit_predict(X) inertia_list.append(km.inertia_) sil_list.append(silhouette_score(X, labels)) fig, axes plt.subplots(1, 2, figsize(12, 4)) axes[0].plot(k_range, inertia_list, markero) axes[0].set_xlabel(k); axes[0].set_ylabel(Inertia) axes[1].plot(k_range, sil_list, markero) axes[1].set_xlabel(k); axes[1].set_ylabel(Silhouette) plt.tight_layout() plt.savefig(k_selection.png, dpi150)random_state42是为了让结果可复现n_init10表示KMeans会跑10次随机初始化选最优这两个参数必须固定否则每次运行的聚类结果都可能不同。肘部法则看inertia下降的拐点——曲线由陡变缓的“肘部”对应合适的k轮廓系数则是越大越好取值范围[-1,1]0.3以上可以接受0.5以上说明聚类结构比较清晰。实际跑BRCA数据时肘部法则经常只能看出两到三类的拐点轮廓系数随k增大单调下降。这不是代码写错了而是高维数据的常态。我的处理策略是指标给出一个合理区间比如3到6再用生物学知识做约束——乳腺癌至少能分ER阳性和ER阴性两大组所以k不能小于2PAM50是五类所以k5可以作为一个候选。两个维度一交叉就能锁定一两个k值做后续验证。4.2 层次聚类Ward法与热图KMeans给的是硬分区但乳腺癌分子分型更像是连续谱系层次聚类能把样本间的亲疏关系完整保留下来。生信文章里最常见的层次聚类配套图是带注释的热图行是基因、列是样本树状图分别展示基因和样本的聚类关系。scipy可以直接完成。from scipy.cluster.hierarchy import linkage, dendrogram from scipy.spatial.distance import pdist import seaborn as sns # 样本间距离 Ward连接 dist_mat pdist(X, metriceuclidean) link linkage(dist_mat, methodward) # 画树状图 plt.figure(figsize(12, 6)) dendrogram(link, labelsX.index, leaf_rotation90, leaf_fontsize4) plt.savefig(hierarchical_dendrogram.png, dpi150) # 用seaborn出带聚类的热图 sns.clustermap(expr_std, methodward, metriceuclidean, cmapvlag, figsize(14, 10)) plt.savefig(clustermap.png, dpi150)Ward法的特点是最小化合并后的类内方差对团状数据表现好是基因表达聚类最常用的连接方法。metriceuclidean和Ward是官方推荐组合——Ward法要求输入是欧氏距离换成别的距离会得出无明显意义的树状图。整条链路是pdist算距离矩阵linkage跑层次聚类dendrogram画树。clustermap则是把热图和两侧树状图一起画出来节省不少排版功夫。层次聚类优于KMeans的地方在于你能从树状图上看到聚类的层级结构。乳腺癌数据里经常出现的一个现象是前两大分支把样本分成Basal样和非Basal样两组再往下细分才拆出Luminal A、Luminal B和HER2富集。这种“先粗分再细分”的结构是KMeans给不出来的。但层次聚类的短板也很明显——一万个样本时画树状图就是一团黑线所以它更适合样本量在几百到两千范围内的BRCA分析。4.3 一致性聚类在R里一遍跑出稳定性单独跑一次KMeans或层次聚类你无法回答一个关键问题这个分组是数据里真实存在的结构还是算法对这个特定数据集的偶然结果TCGA泛癌文章里常用的一致性聚类ConsensusClusterPlus就是专门回答这个问题的。它能重复抽样、多次聚类计算样本对的共现频率最终给出每个k值下的稳定性指标。这种方法在R生态里最成熟Python没有完全等价的包。你不需要二选一——我的做法是Python做预处理和初筛R做一致性聚类验证。# R脚本需要安装ConsensusClusterPlus library(ConsensusClusterPlus) # 读取Python导出的标准化表达矩阵 expr_std - read.table(BRCA_cluster_input.tsv, headerTRUE, row.names1, sep\t) # 转置为样本×基因 expr_t - t(expr_std) set.seed(42) results - ConsensusClusterPlus( d expr_t, maxK 6, # 尝试k2到6 reps 50, # 重复抽样次数正式分析建议100 pItem 0.8, # 每次抽样80%的样本 clusterAlg km, # 使用KMeans distance euclidean, seed 42, plot png ) # 输出每个k值的样本聚类标签 for (k in 2:6) { write.csv(results[[k]]$consensusClass, file paste0(consensus_k, k, .csv)) }reps50是演示用的参数正式分析里我至少跑到100次。pItem0.8意味着每次抽样只用80%的样本这是为了评估聚类对样本波动的鲁棒性。ConsensusClusterPlus会算出一个PAC分数共识矩阵的累积分布函数下的面积分数越低说明聚类越稳定。选k的标准是PAC最小且曲线开始出现平台的位置跟你Python里算的轮廓系数放在一起看趋势一致基本可以放心。一致性聚类跑出来的样本标签可以直接传入后续的差异表达分析和生存分析。注意R的read.table读入时默认行名从第一列取列名从第一行取千万别在导出时多写了一个索引列导致错位。5. 聚类结果避坑检查五个容易翻车的地方5.1 样本ID匹配错位聚类结果漂亮但标签全错现象聚类热图看起来很规则分组的轮廓系数也不低但是所有样本的临床亚型注释对不上号生存分析结果完全无意义。原因样本表型文件比如ER状态、PR状态、PAM50注释和表达谱的样本顺序不一致直接按列合并后整体错位。TCGA样本条码本身是唯一的锚点但许多教程数据会把样本重命名成Sample_1、Sample_2这类序号一旦两头排序方式不同匹配就错位。解决任何时候都以12位TCGA患者ID为key做左连接做完后随机抽5个样本打印核对一遍确认ER状态和已知标记基因比如ESR1的表达相关后再做分析。这个核对成本极低却是我见过最多人栽的地方。5.2 拿原始Counts或FPKM直接聚类现象聚类第一主成分把所有高总表达量的样本聚在一起低表达样本另成一簇跟任何已知亚型都对不上。原因前面第2.3节说过原始Counts受测序深度影响不代表真实生物学差异。解决统一按log2(x1)变换后再进下游。这里有个补充如果数据是FPKM或TPM本身已经做过长度校正可以不用再除以总reads数直接log2即可。千万别在log2之后又做一次CPM归一化那相当于二次标准化会把真实差异抹掉一部分。5.3 特征基因误用差异表达基因现象聚类结果能分出来但分出来的组跟肿瘤纯度比如免疫细胞浸润比例强相关而不是分子亚型。原因做聚类之前先用所有样本做了差异表达分析、选那些p值最小的基因作为特征。这种做法有个逻辑漏洞差异表达是分组后算的而你还没有分组——用全体样本找差异基因找到的往往是肿瘤和正常组织差异或者高肿瘤纯度样本和低纯度样本的差异一进聚类就分出个“纯度轴”。解决用无监督的特征筛选比如第3章的高变基因或者方差排名。差异表达应该在聚类分组完成之后做用来描述每个组的分子特征而不是反过来指导分组。5.4 k值只信算法指标不看生物学意义现象轮廓系数说k2最好但也只能分出个大样本组和一个小样本组硬切到k5轮廓系数掉到0.2以下以为做错了。原因这在生物数据里太常见了。BRCA的分子亚型不是五个离散球而是连续谱系Luminal A和Luminal B之间有许多中间态样本。算法指标理论上要找“紧密且分离”的簇但生物数据里很多簇是“紧密但连续”。解决把轮廓系数、肘部法则、PAC分数和PAM50亚型注释放在一起判读。分区能大致对应已知亚型、有清晰的临床特征差异比单看统计指标更重要。我的习惯是先跑k3到k6把每个k下分组的临床注释分布画出来人眼看一眼就知道哪个k讲得出故事。5.5 聚类结果不可复现每次运行分组都略有不同现象同一个输入数据、同一段脚本上午跑和下午跑分组不完全一致尤其是在组的边界上。原因KMeans的初始质心是随机选的层次聚类本身确定但样本顺序改动会影响中间过程的距离计算。到处都有随机性。解决全局固定random_state和n_init。KMeans里设random_state42, n_init100n_init越大结果越稳定但越慢100的计算量对BRCA数据完全可接受。R里一律set.seed(42)。另外每次跑完聚类把标签输出为CSV存档后面对照分析直接用存档标签不要每次都重新聚类。5.6 内存爆炸两万基因乘一千样本的矩阵放不进内存现象脚本跑到距离矩阵计算时直接卡死或者系统内存占用飙升到百分之百。原因pdist计算距离矩阵的复杂度是O(n²)级别样本数一千左右还勉强接受如果你把两万基因全保留KMeans每次迭代都要计算一万维的距离速度极慢。解决高变基因筛选后降到3000个基因再做聚类。如果做了特征筛选还不够可以考虑用PCA先降到50个主成分再聚类但主成分的生物学可解释性比基因列表差我一般不建议做主成分聚类除非只是做可视化。6. 验证聚类结果热图、PCA和生存曲线的三板斧聚类跑完不等于分析结束。你给出一组分组标签审稿人第一个问题一定是这个分组有生物学意义吗验证工作有三板斧热图看标记基因是否按组分布、PCA看样本在低维空间是否成团、生存曲线看各组预后是否有差异。三者都过聚类才算真正落地。# 标记基因验证把PAM50相关基因的表达分布画成热图 pam50_markers [ESR1, ERBB2, MKI67, PGR, EGFR] marker_expr expr_log.loc[[g for g in pam50_markers if g in expr_log.index], :] plt.figure(figsize(10, 4)) sns.heatmap(marker_expr, cmapvlag, center0) plt.title(PAM50 marker genes across samples) plt.savefig(pam50_markers_heatmap.png, dpi150) # PCA降维可视化 from sklearn.decomposition import PCA pca PCA(n_components2) X_pca pca.fit_transform(X) pca_df pd.DataFrame(X_pca, columns[PC1, PC2], indexX.index) pca_df[cluster] labels plt.figure(figsize(8, 6)) sns.scatterplot(datapca_df, xPC1, yPC2, huecluster, palettetab10) plt.savefig(pca_clusters.png, dpi150)热图里ESR1高表达组对应Luminal亚型ERBB2高表达组对应HER2富集型MKI67高表达提示增殖活跃这些模式跟文献对得上分组就可信。PCA散点图则是最直观的“聚类有效性”佐证——如果不同簇在二维平面上混成一团聚类大概率只是高维空间里的数学产物换成别的数据同样能分出来。生存分析可以继续用Python的lifelines或者把标签导入R的survival包。这一步我不展开写完整代码但逻辑上要看各组Kaplan-Meier曲线的中位生存期和log-rank检验p值。Basal样亚型通常预后差这几乎是乳腺癌领域的共识如果你的聚类结果里Basal样组的生存曲线反而最好先别急着改结论回头排查一下样本ID匹配的问题——大概率是标签对错了患者。我现在拿到任何一个聚类分析项目第一件事就是先把样本标签和临床表的样本数量核对一遍然后再问聚类算得对不对。这个习惯救过我太多次了。聚类分析的价值不在那张树状图多漂亮而在你的分组能不能在独立数据集上复现、能不能解释患者的预后差异。顺着这条路走下去你会发现TCGA-BRCA数据只是起点同样的流程换到LUAD肺癌或KIRC肾癌上只需要改动下载的数据集ID和标记基因列表。希望帮到你。本文还有配套的精品资源点击获取