TCGA-BRCA聚类分析：R语言层次聚类与PCA实战指南-北京尧图网络

简介面向生物信息学初学者与R数据分析实践者该资源围绕TCGA-BRCA乳腺癌基因表达数据设计了一套完整的聚类分析方案涵盖层次聚类距离默认取average与PCA降维两大核心任务并利用临床信息中的ER_Status_nature2012分类标准对聚类结果进行验证帮助读者理解如何根据基因表达水平对病人分群。压缩包共22个文件大小约10.91MB主要包括R源码、基因表达矩阵与临床数据文件、8张PNG结果图及对应PDF版图以及README说明文档便于对照代码与图表复盘每一步分析逻辑。资源已有687人学习下载适合正在完成生物信息学课程作业或希望掌握R聚类与PCA实操的读者。通过研读源码与输出图可快速复现层次聚类热图、PCA碎石图与主成分累积贡献图理解主成分数目的选择依据并掌握降维后重新聚类及与原始分类比较的方法。1. TCGA-BRCA 聚类分析一份能直接复现的 R 代码包拿到「生物信息学概论——聚类分析TCGA-BRCA数据.zip」这份资源的时候我第一反应是这又是一个课程作业的参考答案。但解压之后发现里面不只是代码还有完整的 Figures-PNG 和 Figures-PDF 两套输出图片、基因表达矩阵和临床注释文件结构上已经接近一个微型生信项目的标准布局。这份资源解决的核心问题是给定一个基因表达矩阵如何用层次聚类和 PCA 两种方法对乳腺癌病人分组并用 ER 状态这个临床金标准去验证聚类结果是否合理。适合正在学生信概论、被作业卡住的同学也适合想看看 TCGA 数据标准分析流程的 R 入门者。我要说的是这份代码能跑通但里面有几个决定结果好坏的细节——距离矩阵选错、数据没标准化、病人 ID 没对齐——如果不处理聚类图就是一团浆糊。2. 数据进场先看格式GeneMatrix 与 clinical_data 的对齐关系2.1 两个文件的字段含义与读取方式GeneMatrix.txt 是表达值矩阵行是基因列是病人。注意它和 clinical_data.txt 的病人 ID 不是一一对应关系GeneMatrix 里的病人只是 clinical 数据的一个子集。这意味着第一步不是读文件而是做病人交集筛选。我在处理 TCGA 数据时习惯先看维度再看内容避免读完发现几十个病人的样本名对不上。# 读取表达矩阵行为基因列为病人第一列是基因名 expr_raw - read.delim(GeneMatrix.txt, row.names 1, check.names FALSE) # 读取临床信息每行一个病人 clin_raw - read.delim(clinical_data.txt, row.names 1, check.names FALSE) # 取交集病人这一步是后续所有分析正确的前提 common_patients - intersect(colnames(expr_raw), rownames(clin_raw)) expr - expr_raw[, common_patients] clin - clin_raw[common_patients, , drop FALSE]这里check.names FALSE是关键。R 默认会把列名里的-转成.如果你不关掉这个行为后面取交集时永远匹配不上。read.delim默认分隔符是制表符TCGA 的矩阵文件一般也是这个格式如果你打开发现是空格分隔的改成read.table(sep )就行。交集做完之后建议顺手检查一下length(common_patients)确认剩下多少个病人。我见过不少人忽略这一步到画图时才发现病人数比预期多了一倍。2.2 ER_Status_nature2012 字段评估聚类结果的依据临床数据里有一列叫ER_Status_nature2012这是 2012 年 Nature 论文里定义的 ER 状态注释。ER 即雌激素受体乳腺癌病人按照 ER 阳性或阴性在治疗方案和预后上有显著差异。如果聚类能把 ER 阳性和阴性的病人大致分开说明你用的基因表达特征确实携带了生物学信号。# 查看 ER 状态分布 table(clin$ER_Status_nature2012) # 如果存在空值或 Unknown决定是否剔除 # 有些作业数据里会有 NA 或 Unknown建议直接剔除 valid_idx - clin$ER_Status_nature2012 %in% c(Positive, Negative) expr - expr[, valid_idx] clin - clin[valid_idx, , drop FALSE] expr - expr[apply(expr, 1, function(x) all(is.finite(x))), ]这一步做两件额外的事第一过滤掉 ER 状态不明确的样本避免后期评估时出现无法归类的点第二过滤表达矩阵中非有限值的基因行。TCGA 的表达矩阵偶尔会出现空值或极端值如果不处理后面的距离计算直接报错。我一般还会顺手做个基因表达量的分布直方图确认数据是不是已经做过 log2 归一化——这决定了后面聚类时是否需要额外的标准化处理。3. 层次聚类实战从距离矩阵到热图3.1 为什么选择 average 距离题目明确要求距离选择 average这其实是层次聚类里最稳的选择。ward 倾向生成紧凑的球状簇在基因表达数据上容易过度分割single 容易产生链状效应一组样本被逐个拽进一个大簇里聚类树看起来像一条长链。average 即 UPGMA计算两个簇之间所有样本对距离的均值对离群点不那么敏感在转录组数据上表现相对稳健。表达矩阵聚类的另一个关键决策是对行做标准化还是对列做标准化或者都不做。基因表达数据的量纲差异很大高表达的基因会主导距离计算。常见做法是先对每个基因做 z-score 标准化让所有基因在同一尺度上比较。# 转置后按行基因做 z-score 标准化 expr_z - t(scale(t(expr))) # 计算样本间距离矩阵method average 是题目指定 dist_mat - dist(t(expr_z), method euclidean) # 层次聚类 hc - hclust(dist_mat, method average) # 聚类树可视化的最基本形态 plot(hc, labels FALSE, main Hierarchical Clustering of TCGA-BRCA Samples)scale(t(expr))做的事情是对列做标准化转置之后就成了对每个基因做标准化。这一步非常容易反了我见过有人直接scale(expr)那样是对病人做标准化聚类结果会变得毫无意义。对基因做 z-score 的含义是每个基因的表达值分布被拉到均值为 0、方差为 1 的尺度上这样高表达基因和低表达基因在距离计算中是平等的。dist函数默认计算欧氏距离这也是表达谱聚类最常用的选择。3.2 热图判断聚类效果最直观的武器聚类树只是给出了样本的分组关系真正能看出生物学意义的还是热图。R 里做热图有很多选择基础版的heatmap函数简单直接pheatmap更灵活。这份资源里的 heatmap.PNG 就是典型的样本-基因双聚类图。# 选择变异最大的前 1000 个基因降低噪声 gene_var - apply(expr_z, 1, var) top_genes - names(sort(gene_var, decreasing TRUE))[1:1000] # 生成热图列按层次聚类结果排序 pheatmap::pheatmap( expr_z[top_genes, ], cluster_cols hc, cluster_rows TRUE, scale none, show_rownames FALSE, show_colnames FALSE, annotation_col data.frame( ER_Status factor(clin$ER_Status_nature2012), row.names colnames(expr_z) ), main Heatmap of Top 1000 Variable Genes )选前 1000 个高变异基因是一个常见做法理由是在所有基因上用聚类不表达或表达量极低的基因会稀释真正的信号。高变异基因更可能携带区分样本的生物学差异。annotation_col参数把 ER 状态作为热图顶部的注释条展示出来这能直观看到聚类结果与 ER 状态的对应关系。如果聚类效果好你应该能看到热图左侧的样本树分成两到三个大簇注释条的颜色块与簇的划分有清晰的对应。3.3 为什么你的聚类树总是乱糟糟的聚类树乱通常不是代码问题而是数据预处理问题。第一没有做基因过滤全部基因参与聚类噪声淹没信号第二没有做标准化个别高表达基因主导了距离矩阵导致所有样本之间的距离都差不多大。如果你发现聚类树最后才分成两簇即大部分样本在树顶部还没有分离大概率是标准化没做对。第三样本量太少。TCGA-BRCA 的完整数据集有上千例但这份资源里只给了部分病人如果只有几十例聚类结构不稳定是正常的。# 检查聚类树的分簇质量计算共表型相关系数 cor_coph - cor(cophenetic(hc), dist_mat) cat(Cophenetic correlation:, cor_coph)共表型相关系数是一个简单有效的聚类质量指标数值越接近 1说明聚类树保留的距离信息越完整。如果这个值低于 0.75说明你的数据预处理方式可能不太适合层次聚类可以尝试换标准化方式或过滤策略。4. PCA 降维与二次聚类从高维空间找主结构4.1 主成分数目的选择碎石图与累计方差解释率PCA 在基因表达数据上的作用是把成千上万个基因的表达模式压缩成几个综合变量——主成分。每个主成分是基因的线性组合第一个主成分捕获数据中最大的方差方向。选多少个主成分合适没有绝对标准但两个常用判据是Kaiser 准则特征值大于 1和累计方差解释率达到 80% 左右。# 对转置后的表达矩阵做 PCA样本为行基因为列 pca_result - prcomp(t(expr_z), center TRUE, scale. FALSE) # 查看每个主成分解释的方差比例 var_explained - summary(pca_result)$importance[2, ] cum_var - summary(pca_result)$importance[3, ] # 碎石图主成分编号与方差解释率 plot(var_explained[1:20], type b, xlab Principal Component, ylab Proportion of Variance Explained, main Scree Plot for TCGA-BRCA Expression Data) # 累计方差解释率 plot(cum_var[1:20], type b, xlab Number of Principal Components, ylab Cumulative Variance Explained, main Cumulative Variance Explained)这份资源里的 pca_screeplot.PNG 和 pca_cumulative.PNG 就是这两张图的标准输出。我的习惯是看累计方差曲线在哪个位置开始趋于平缓选取拐点处的主成分数目。一般取前 5 到 10 个主成分能解释 60% 到 70% 的方差就够了不用强求 80%因为基因表达数据的噪声很大强行追求高解释率会把噪声也纳入模型。4.2 基于 PCA 特征的聚类与基因空间的对比PCA 降维之后每个样本被表示成一个低维向量。这个向量也可以作为聚类分析的基础和直接用全部基因聚类形成对比。两者的差异在于直接基因聚类等于给每个基因一样的权重PCA 聚类则是按方差贡献加权低信息量的基因被自动降权。# 提取前 5 个主成分作为新特征 pca_scores - pca_result$x[, 1:5] # 基于 PCA 特征重新计算距离并聚类 dist_pca - dist(pca_scores, method euclidean) hc_pca - hclust(dist_pca, method average) # 对比两次聚类的分组一致性用 k3 的分组 group_gene - cutree(hc, k 3) group_pca - cutree(hc_pca, k 3) table(group_gene, group_pca)cutree把层次聚类树切成指定数量的簇。用交叉表对比两次聚类结果可以看到 PCA 降维之后的分组和基因空间的分组重合度如何。如果重合度很高说明前 5 个主成分已经抓住了数据的主结构如果重合度低说明你选的主成分数目可能太少丢失了分组相关的信息。这个对比逻辑非常重要因为它连接了两种聚类方案并给出了量化评估。4.3 PCA 散点图与 ER 状态的叠加PCA 最直观的展示方式是散点图把前两个主成分作为坐标轴每个样本一个点用 ER 状态给点上色。如果 ER 阳性和阴性的病人在散点图上形成两个分离的区域说明数据中的主要变异方向确实和 ER 状态相关。# 前两个主成分的散点图颜色按 ER 状态 plot(pca_result$x[, 1], pca_result$x[, 2], col ifelse(clin$ER_Status_nature2012 Positive, dodgerblue, firebrick), pch 16, cex 0.8, xlab paste0(PC1 (, round(var_explained[1] * 100, 1), %)), ylab paste0(PC2 (, round(var_explained[2] * 100, 1), %)), main PCA of TCGA-BRCA Samples Colored by ER Status) legend(topright, legend c(ER, ER-), col c(dodgerblue, firebrick), pch 16)x 轴标签里带上 PC1 解释的方差百分比这是让图更专业的小技巧。如果第一主成分就能解释 20% 以上的方差且 ER 和 ER- 在 PC1 方向上有明显分离说明这个数据集的聚类基础主要是 ER 相关的转录程序差异。反之如果两个颜色完全混在一起要么是数据质量有问题要么是 ER 状态本身不是这个数据集中最强的分组信号。5. 避坑指南TCGA 表达数据聚类的四个高频翻车点5.1 病人 ID 对不上导致聚类结果全错现象聚类热图画出来了但 ER 注释条的颜色与簇的划分完全随机看不出任何规律。原因是最初读取数据时没有取 GeneMatrix 和 clinical_data 的交集导致表达矩阵和临床信息中病人的顺序不一致ER 状态标注是错位的。解决在画图前强制用intersect取交集然后用rownames(clin)对样本顺序重新排序确保expr的列名和clin的行名顺序完全一致。# 强制顺序对齐 clin - clin[colnames(expr), , drop FALSE]这行代码看着简单但它是热图注释正确的前提。我自己的习惯是每次画图前都跑一下这个对齐它比任何复杂的分析技巧都重要。5.2 基因没做标准化导致距离矩阵失效现象聚类树看起来所有样本都挤在一起分支很短cophenetic 相关系数低于 0.6。原因是原始表达值中高表达基因的绝对数值较大在欧氏距离计算中完全压倒了低表达基因的贡献。解决先对每个基因做 z-score 标准化再用标准化后的矩阵计算距离。注意scale默认是对列操作转置顺序不能错。5.3 表达值没做 log 转换导致数据偏态现象热图颜色一片深红或一片深蓝看不出梯度差异。原因是TCGA 的 RNA-Seq 表达值一般是 Count 或 TPM数值范围从 0 到几十万直接画热图时颜色标尺会被少数高表达基因拉爆。解决检查数据分布如果不是 log 化的先做log2(expr 1)再标准化。1是为了避免 log(0) 产生负无穷。# 检查表达值范围 range(expr_raw) # 如果最大值超过 1000考虑做 log 转换 expr_log - log2(expr_raw 1)这一步看着简单但它是热图可读性的决定性因素。我看过不少人跳过 log 直接聚类聚类树也凑合能看但那是因为高表达基因恰好是区分样本的关键基因——这种运气不常有。5.4 ER 状态里有坑不是只有 Positive 和 Negative现象ER 状态评估时出现第三类或缺失值导致注释条颜色对不上分组评估的表格里多出一行。原因是临床数据里存在 NA、Unknown、Indeterminate 等额外取值直接把整个列转成因子时这些值也变成了一个新的类别。解决在分析前就过滤掉非正负两类的样本或者单独处理成 NA 后在热图注释里也标成灰色。我一般倾向直接过滤因为这本来就是验证聚类的辅助信息不需要为此增加分析复杂度。6. 进阶验证轮廓系数定量评估聚类的合理性层次聚类和 PCA 聚类都做完之后需要回答一个量化问题到底分几类是合理的以及两次聚类哪个分组质量更高这时候轮廓系数Silhouette Coefficient是最好用的工具。它同时衡量了样本与自身簇的紧密度和与最近邻簇的分离度取值在 -1 到 1 之间越接近 1 说明分组越合理。# 安装并加载 cluster 包 # 评估不同 k 值簇数下的平均轮廓宽度 library(cluster) # 基于基因空间的聚类评估 2 到 6 簇 sil_scores - sapply(2:6, function(k) { groups - cutree(hc, k k) sil - silhouette(groups, dist_mat) mean(sil[, 3]) }) names(sil_scores) - 2:6 print(sil_scores)对这份 TCGA-BRCA 数据我跑下来的典型结果是 k2 时轮廓系数最高k3 时略有下降。这与 ER 状态的生物学背景是吻合的ER 和 ER- 是乳腺癌中最基本的分型。如果 k3 或 k4 得分也很高说明数据里可能存在更细的分组结构那就可以去看每个簇中 ER 状态的组成比例找出额外的生物学解释。# 按 k3 分组看每个簇中 ER 状态的分布 groups_final - cutree(hc, k 3) table(groups_final, clin$ER_Status_nature2012)还有一个用交叉验证思路做分组的技巧就是把数据随机分成训练集和验证集只在训练集上做聚类然后用验证集去检测分组是否稳定。常见做法是用 bootstrap 重采样每次抽样后重新聚类看哪些样本对始终被分到同一簇。R 的pvclust包可以做这件事它会给出每个分支的 bootstrap 支持率支持率低于 95% 的分支不要当作真实分组来解释。不过要注意pvclust在几十个样本的小数据上计算量不大但如果全量基因参与速度会很慢建议还是先过滤到前 1000 个高变异基因再跑。# pvclust 的快速示例基于标准化表达矩阵的转置 library(pvclust) set.seed(123) # 注意pvclust 需要行为样本、列为变量这里转置为病人行、基因列 pv - pvclust(t(expr_z[top_genes, ]), method.dist euclidean, method.hclust average, nboot 100) plot(pv)bootstrap 数设 100 就够了如果设 1000小数据集可能要跑十几分钟。看输出的聚类树时重点关注au值Approximately Unbiased它比bp值更准确。从那以后我每次做聚类分析都会强制走一遍这个验证流程先看共表型相关系数确认聚类树质量再用轮廓系数对比不同 k 值最后用 bootstrap 检查关键分支的稳定性。这套流程帮我挡掉了很多看着漂亮其实站不住脚的聚类结果希望也能帮到你。这份资源里的 R 代码和输出图已经覆盖了上述所有核心步骤——数据读入、层次聚类、PCA 降维、热图可视化、ER 状态验证——直接照着跑一遍再按你选定的主成分数目定制参数就能得到一套完整可用的分析结果。本文还有配套的精品资源点击获取

TCGA-BRCA聚类分析：R语言层次聚类与PCA实战指南

相关推荐

OpenClaw 搭建智能运维巡检工作流实践：用 Skill 编排 Agent 自动巡检

编译原理课程设计：词法分析器、LL(1)与LR(1) Java项目全解

Pandoc 老将 vs MarkItDown 新王：AI 数据流水线到底该选谁？

更多资讯

工业连接器选型详解：EDAC矩形与D-Sub接口如何避坑

BLE广播、扫描与连接：协议详解与工程组合实践

神经网络信号流解析：从输入层到输出层的可调试模块拆解

论文降重工具真实测评：查重机制与组合避坑指南

MATLAB电力系统故障仿真：从序网建模到暂态波形全链路实现

滑动窗口算法全解析：三种形态、单调队列与工程应用

今日更新

再见，SSE！你好，Streamable HTTP！轻松开发 Streamable HTTP MCP Server 并接入 TaoToken

工业软件标准化路线图：国产替代的落地施工图

VCMI安卓版实操指南：原生运行英雄无敌3的3步技术落地

本周精选

CMOS设计三大物理效应：体效应、沟道长度调制与亚阈值导电

微信小程序在线阅读毕设全解析：Java后端+MySQL+避坑实战

Android音频问题五层分析法：从App到Codec的系统化诊断

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？