简介：这份资源面向生物信息学入门学习者与R语言数据分析实践者，围绕TCGA-BRCA乳腺癌基因表达数据展开聚类分析练习。内容涵盖利用层次聚类按基因表达水平对病人分型、选择average距离度量、绘制heatmap等图示，并实现PCA降维后重新聚类，与原始结果对比评估，同时借助临床数据中的ER_Status_nature2012标签检验聚类是否符合预期。压缩包共22个文件，约10.91MB，包含8个png与8个pdf结果图、2个md说明、2个txt数据文件、1个license及1个R源码脚本，图表与代码配套便于复现。已有687人学习下载。读者可获得完整的数据分析流程、可运行的R脚本、聚类与降维的可视化结果，以及基于临床标签的评估思路，适合作为课程作业或生信聚类实战参考。

1. 拿到 TCGA-BRCA 表达矩阵之后：这份 R 聚类源码包到底能跑出什么

如果你手上正好有一份 TCGA-BRCA 的基因表达矩阵，却卡在“怎么把病人按表达谱分群”这一步，这个压缩包值得先拆开看看。它把生物信息学概论里最经典的聚类分析流程做成了可复现的 R 工程：输入是GeneMatrix.txt和clinical_data.txt，输出是一整套聚类图、PCA 图、热图，代码集中在cluster.R，图同时给了 PNG 和 PDF 两种格式。换句话说，它不是只给你一段演示代码，而是把“层次聚类 + PCA 降维 + 再聚类 + 用 ER 状态评估”这条链路完整落到了文件和图上。适合正在做课程设计、想复现 TCGA 聚类流程、或者需要一份能直接改参数跑自己数据的从业者。下面我按实际拆包顺序，把这份资源怎么用、参数怎么设、哪里容易翻车讲清楚。

2. 拆开压缩包先看什么：文件结构与数据格式核对

2.1 目录里每个文件对应哪一步分析

拿到生物信息学概论——聚类分析TCGA-BRCA数据.zip，解压后不要急着跑cluster.R，先把文件按用途分三类，后面排错会快很多。

图目录里cluster.PNG、pca_cluster.PNG、heatmap.PNG、ER_heatmap.PNG、ER_cluster.PNG这几张是评估聚类效果的核心，pca_screeplot.PNG和pca_cumulative.PNG用来定主成分数目，pca_heatmap.PNG是降维后的热图。先看这些图，能大致判断脚本跑出来的结果是否符合预期。

2.2 读入前必须核对的两个数据格式

GeneMatrix.txt和clinical_data.txt的列名必须能对上，否则后面按 ER 状态评估时会直接错位。常见做法是先单独读一遍，确认行名、列名和维度。

# 读入表达矩阵，check.names=FALSE 防止列名里的连字符被改掉 gene_matrix <- read.table("GeneMatrix.txt", header = TRUE, row.names = 1, sep = "\t", check.names = FALSE) # 读入临床信息，第一列通常是病人编号 clinical <- read.table("clinical_data.txt", header = TRUE, row.names = 1, sep = "\t", check.names = FALSE) # 核对维度：行是基因，列是病人 dim(gene_matrix) dim(clinical) # 核对病人编号交集，这一步不做后面 ER 评估必翻车 common_samples <- intersect(colnames(gene_matrix), rownames(clinical)) length(common_samples)

逻辑说明：row.names = 1把第一列当作行名，check.names = FALSE保留原始列名，避免 R 自动把TCGA-XX-XXXX这类编号改得对不上。intersect用来确认表达矩阵和临床信息里共有的病人数量，如果交集明显小于表达矩阵列数，说明两份数据的编号体系不一致，需要先统一。参数上，sep = "\t"对应制表符分隔，如果你的文件是逗号分隔，改成sep = ","。

提示：先跑dim()和head()看数据，不要直接进聚类。数据没对齐，后面所有图都是错的。

3. 层次聚类怎么落地：距离选 average 的完整代码链路

3.1 为什么先做样本聚类而不是基因聚类

这份数据的分析目标是“按基因表达水平把病人分类”，所以聚类对象是样本（列），不是基因（行）。常见做法是先转置矩阵，让行变成病人、列变成基因，再算样本间距离。距离用average（即 UPGMA），这是题目明确要求的，原因是它在样本量不大、表达谱噪声较高时比complete更稳健，不会因为个别极端值把整棵树拉偏。

# 转置：行变成病人，列变成基因 expr_t <- t(gene_matrix) # 算样本间距离，method 可选 "euclidean"、"manhattan"、"correlation" dist_mat <- dist(expr_t, method = "euclidean") # 层次聚类，average 即 UPGMA hc <- hclust(dist_mat, method = "average") # 画树状图，main 里写清楚用的是 average plot(hc, main = "Hierarchical Clustering (average)", xlab = "", sub = "")

逻辑说明：dist()默认是欧氏距离，如果基因表达值量纲差异大，可以先做标准化再算距离。hclust()的method参数支持"ward.D2"、"complete"、"single"等，这里按题目要求用"average"。plot()出来的树状图就是cluster.PNG对应的内容。参数上，dist()的method和hclust()的method是两个独立选择，不要混为一谈。

3.2 切树得到病人分群并输出热图

树画出来只是第一步，真正要的是每个病人属于哪一类。用cutree()按 k 切分，k 的选择可以结合树状图高度和临床预期。切完之后用热图看分群是否在表达层面清晰。

# 按 k=2 切分，对应 ER 阳性/阴性两类预期 clusters <- cutree(hc, k = 2) table(clusters) # 热图：需要矩阵形式，先转回基因 x 病人 heatmap(as.matrix(gene_matrix), ColSideColors = ifelse(clusters == 1, "blue", "red"), scale = "row", main = "Heatmap with cluster assignment")

逻辑说明：cutree()的k是切分簇数，也可以换成h按高度切。table(clusters)先看每类有多少样本，如果一类只有一两个样本，说明 k 选大了或者距离度量不合适。heatmap()里scale = "row"表示按基因标准化，这样热图颜色反映的是相对表达高低，而不是绝对数值。ColSideColors把聚类结果标在列上方，方便和 ER 状态对比。

注意：heatmap()是 R 基础包函数，样本量超过几百时渲染会很慢，可以考虑pheatmap或ComplexHeatmap，但这份源码用的是基础函数，先按原脚本跑通再换。

4. PCA 降维与再聚类：主成分数目怎么定、和第一次聚类差在哪

4.1 PCA 实现与碎石图、累计方差图

PCA 的目的是把高维基因表达压到少数几个主成分，再用这些主成分重新聚类，看结果是否和第一次一致。prcomp()是 R 里最常用的实现，默认对变量做中心化。

# PCA，scale.=TRUE 表示同时做标准化 pca_res <- prcomp(expr_t, scale. = TRUE) # 碎石图：看每个主成分解释的方差 plot(pca_res, type = "l", main = "Scree Plot") # 累计方差图：辅助决定保留几个主成分 cum_var <- cumsum(pca_res$sdev^2 / sum(pca_res$sdev^2)) plot(cum_var, type = "b", xlab = "Principal Component", ylab = "Cumulative Variance Explained", main = "Cumulative Variance") abline(h = 0.8, col = "red", lty = 2)

逻辑说明：prcomp()的scale. = TRUE会在 PCA 前对每个基因做标准化，避免高表达基因主导主成分。pca_res$sdev是每个主成分的标准差，平方后除以总和就是方差解释比例。碎石图对应pca_screeplot.PNG，累计方差图对应pca_cumulative.PNG。abline(h = 0.8)是常见的 80% 累计方差参考线，但具体保留几个主成分要看拐点和实际聚类效果。

4.2 选几个主成分：拐点法加累计方差双条件

主成分数目没有唯一正确答案，但可以用两个条件交叉判断：累计方差达到 80% 左右，且碎石图出现明显拐点。常见做法是取前 2 到 5 个主成分，因为这份数据的样本量不大，主成分太多会引入噪声，太少又可能丢掉分群信息。

# 取前 3 个主成分 pc_use <- pca_res$x[, 1:3] # 用主成分重新算距离并聚类 dist_pca <- dist(pc_use, method = "euclidean") hc_pca <- hclust(dist_pca, method = "average") plot(hc_pca, main = "Hierarchical Clustering on PCA (average)") # 再切分，和第一次聚类对比 clusters_pca <- cutree(hc_pca, k = 2) table(clusters, clusters_pca)

逻辑说明：pca_res$x是样本在主成分上的得分矩阵，取前几列就是降维后的特征。table(clusters, clusters_pca)是两次聚类的交叉表，如果大部分样本落在对角线上，说明 PCA 降维后保留了主要分群结构；如果交叉表很乱，可能是主成分数目不合适，或者第一次聚类本身就不稳定。参数上，1:3可以改成1:2或1:5做敏感性测试。

4.3 用 ER 状态评估聚类：交叉表与热图双重验证

clinical_data.txt里的ER_Status_nature2012是评估聚类是否合理的标签。把聚类结果和 ER 状态做交叉表，再看ER_heatmap.PNG和ER_cluster.PNG对应的图，就能判断聚类是否“符合预期”。

# 提取 ER 状态，确保病人顺序和聚类结果一致 er_status <- clinical[names(clusters), "ER_Status_nature2012"] table(clusters, er_status) # 用 ER 状态给热图列上色 heatmap(as.matrix(gene_matrix), ColSideColors = ifelse(er_status == "Positive", "blue", "red"), scale = "row", main = "Heatmap colored by ER status")

逻辑说明：clinical[names(clusters), ]按聚类结果的病人顺序取临床信息，避免顺序错位。table(clusters, er_status)看聚类簇和 ER 阳性/阴性的对应关系，如果某一簇里 ER 阳性占绝大多数，说明聚类捕捉到了生物学信号。ColSideColors换成 ER 状态后，热图能直观看出 ER 相关基因是否在两类病人间差异表达。

提示：ER 状态只是评估标签之一，不是绝对标准。聚类结果和 ER 不完全一致时，先检查数据标准化和主成分数目，不要直接否定聚类。

5. 避坑与排查：这份源码跑不通时先查这五处

5.1 现象：读入后列名变成TCGA.XX.XXXX，和临床信息对不上

原因：read.table()默认check.names = TRUE，会把连字符、空格等替换成点号。解决：读入时显式设置check.names = FALSE，并在读入后立刻用intersect()核对病人编号交集。

5.2 现象：热图报错“cannot allocate vector of size”，或画出来一片糊

原因：heatmap()对矩阵维度敏感，基因数太多或样本太多时会内存不足，且不做筛选时热图没有可读性。解决：先按方差或表达量筛掉低变异基因，再画热图；样本量超过 200 时改用pheatmap，并设置cluster_rows = TRUE、cluster_cols = TRUE。

5.3 现象：PCA 碎石图和累计方差图对不上，不知道选几个主成分

原因：只看累计方差 80% 可能保留过多主成分，只看拐点又可能太少。解决：两个条件一起用，先看拐点，再确认累计方差是否接近 80%，最后用table(clusters, clusters_pca)验证降维后聚类是否稳定。主成分数目在 2 到 5 之间做敏感性测试。

5.4 现象：两次聚类结果差异很大，交叉表几乎不对角

原因：第一次聚类用的是全部基因，第二次用的是前几个主成分，如果主成分数目太少，降维后丢掉了分群信息；如果太多，又引入了噪声。解决：先检查 PCA 前是否做了标准化，再调整主成分数目，同时确认两次聚类都用了average距离。必要时对基因做方差筛选后再跑 PCA。

5.5 现象：ER 评估交叉表里某一类样本数为 0

原因：cutree()的 k 选得太大，或者聚类树本身没有分出对应分支。解决：先看table(clusters)每类样本数，再结合树状图调整 k 或改用h切树。如果 ER 阳性/阴性本身不平衡，交叉表要按行或列算比例，不要只看绝对数。

6. 把这份源码改成自己的数据：参数替换与结果验证的固定习惯

这份资源最大的价值不是那几张图，而是cluster.R里那条可替换参数的链路。换成自己的表达矩阵时，我一般会按固定顺序改四处：第一，读入部分改文件名和分隔符，确认intersect()交集不为空；第二，距离度量先保持euclidean+average，跑通后再试correlation；第三，PCA 的scale.保持TRUE，主成分数目从 3 开始，用累计方差图和交叉表各验证一次；第四，热图先筛低变异基因，再画ER_heatmap对应的版本。

# 换成自己的数据时，按这个顺序改参数 gene_matrix <- read.table("YourMatrix.txt", header = TRUE, row.names = 1, sep = "\t", check.names = FALSE) clinical <- read.table("YourClinical.txt", header = TRUE, row.names = 1, sep = "\t", check.names = FALSE) # 筛低变异基因：保留方差前 2000 个 gene_var <- apply(gene_matrix, 1, var) gene_matrix <- gene_matrix[order(gene_var, decreasing = TRUE)[1:2000], ] # 后续聚类、PCA、热图代码不变，只改输入和 k

逻辑说明：apply(gene_matrix, 1, var)按行算方差，order(..., decreasing = TRUE)[1:2000]取方差最大的 2000 个基因。这一步不是必须，但能显著提升热图可读性和聚类稳定性。参数上，2000 可以按数据规模调整，样本少时取 1000 到 3000 都常见。

验证方法上，我习惯每次改完参数都跑三件事：table(clusters)看每类样本数是否合理，table(clusters, er_status)看聚类和 ER 的对应关系，plot(hc)看树状图有没有明显异常分支。三件事都过了，再去看热图和 PCA 图。从那以后我每次换数据都强制走一遍这个检查顺序，省掉了很多“图好看但结论错”的后悔药。希望帮到你。

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