简介这份资源面向生物信息学入门学习者与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人学习下载。读者可获得完整分析流程、可运行脚本、可视化结果与评估思路适合课程作业与实战练手。1. 从一份 BRCA 表达矩阵说起聚类分析到底在解决什么问题拿到一份 TCGA-BRCA 的基因表达矩阵第一反应往往不是“我要跑什么算法”而是“这堆样本里到底藏了几种人”。乳腺癌在临床上从来不是一个单一疾病同样是 ER 阳性、同样做了内分泌治疗有人五年无事有人两年复发。这种异质性如果只用一张病理切片和几个免疫组化指标去描述信息量是不够的。聚类分析在这里的角色就是把“分子层面的相似性”变成“可以分组、可以命名、可以拿去和生存曲线对照”的东西。这份标题里的“生物信息学概论”加“聚类分析 TCGA-BRCA 数据”落到实操上是一条很具体的链路从公共数据库取表达谱和临床注释做标准化和降维选一种聚类算法把样本切开再用临床变量验证这些簇有没有生物学意义。它适合两类人一类是刚进生信方向、想找一个完整但不至于太重的项目练手的人另一类是做过差异表达、想再往前推一步看样本分型的人。下面我按自己实际跑这类分析的顺序把每一步的参数、坑和判断标准讲清楚。2. 数据获取与预处理把 TCGA-BRCA 变成能喂给聚类算法的矩阵2.1 为什么选 STAR-Counts 而不是 FPKMTCGA 的表达数据在公共镜像里有几种形态常见的是 HTSeq-Counts、STAR-Counts 和 FPKM/FPKM-UQ。聚类分析对“样本间可比性”的要求高于对“绝对表达量”的要求所以我的习惯是优先用 STAR-Counts 或 HTSeq-Counts 这类原始计数然后在 R 里自己做归一化。FPKM 已经做过一次长度和测序深度校正看起来省事但它把后续做 DESeq2 或 edgeR 的余地堵死了而且不同批次之间的 FPKM 可比性并不总是可靠。下载时通常用 TCGAbiolinks 这个 R 包它能把表达矩阵和临床注释一起拉下来。下面这段是我常用的取数骨架注意 project 和 workflow.type 两个参数。# 加载 TCGA 数据获取包 library(TCGAbiolinks) # 查询 BRCA 项目的 STAR-Counts 数据 query - GDCquery( project TCGA-BRCA, data.category Transcriptome Profiling, data.type Gene Expression Quantification, workflow.type STAR - Counts ) # 下载到本地缓存第一次会慢之后走缓存 GDCdownload(query, files.per.chunk 50) # 整理成表达矩阵行是基因列是样本 expr_data - GDCprepare(query) expr_matrix - assay(expr_data, unstranded)逻辑说明GDCquery 只是构造查询条件不产生下载GDCdownload 的 files.per.chunk 控制每次请求的文件数网络不稳时调小到 20 更稳GDCprepare 返回的是 SummarizedExperiment 对象用 assay 取“unstranded”这一层得到原始计数矩阵。参数上workflow.type 必须和 data.type 匹配写成 “STAR - Counts” 时中间有空格和短横写错会直接报错找不到数据。2.2 基因 ID 转换与低表达过滤原始矩阵的行名是 Ensembl ID带版本号直接拿去聚类没有可读性。常见做法是用 org.Hs.eg.db 把 Ensembl ID 映射成基因符号同时把没有对应符号的行去掉。低表达过滤这一步很多人跳过结果是聚类被一堆全零或接近全零的基因主导。library(org.Hs.eg.db) library(dplyr) # 去掉 Ensembl ID 的版本后缀 rownames(expr_matrix) - gsub(\\..*, , rownames(expr_matrix)) # Ensembl ID 转基因符号 symbols - mapIds( org.Hs.eg.db, keys rownames(expr_matrix), column SYMBOL, keytype ENSEMBL, multiVals first ) # 保留有符号且非 NA 的行 keep - !is.na(symbols) expr_matrix - expr_matrix[keep, ] rownames(expr_matrix) - symbols[keep] # 低表达过滤至少在 25% 样本中计数大于 10 min_samples - floor(0.25 * ncol(expr_matrix)) expr_filtered - expr_matrix[rowSums(expr_matrix 10) min_samples, ]逻辑说明mapIds 的 multiVals 设为 “first” 是为了避免一个 Ensembl ID 对应多个符号时返回列表导致后续报错低表达阈值 10 和 25% 不是金标准但对 BRCA 这种样本量在千级别的数据集能过滤掉大量噪声基因同时保留足够多的特征。如果样本量只有几十25% 这个比例要往上调否则过滤太狠。2.3 归一化log-CPM 还是 VST聚类之前必须让样本之间可比。原始计数直接算距离会被测序深度带偏所以要做归一化。常见两条路edgeR 的 log-CPM 和 DESeq2 的 VST。log-CPM 快适合先看整体分布VST 对高表达基因的压缩更温和聚类结果通常更稳。我一般两个都跑看样本相关性热图的差异。library(edgeR) # 构建 DGEList 并计算 log-CPM dge - DGEList(counts expr_filtered) cpm_log - cpm(dge, log TRUE, prior.count 2) # 检查归一化后的分布 boxplot(cpm_log[, 1:8], main log-CPM distribution, las 2)逻辑说明prior.count 是在取对数前给每个计数加的一个小伪计数防止 log(0) 出现设成 2 比默认的 0.25 更保守适合样本间差异较大的情况。箱线图只看前 8 个样本目的是确认中位数大致对齐如果中位数错位明显说明归一化没做好后面聚类基本白跑。3. 降维与聚类从 PCA 到一致性聚类怎么选、参数怎么定3.1 先做 PCA 看有没有明显的批次或离群样本在跑任何聚类之前我会先做一次 PCA。目的不是降维本身而是看样本在二维平面上有没有明显的离群点、有没有按批次聚成条带。如果 PCA 图里样本散成一团说明整体异质性不高强行聚类可能得到没有意义的簇。# 取方差最大的 2000 个基因做 PCA rv - apply(cpm_log, 1, var) top_genes - names(sort(rv, decreasing TRUE))[1:2000] pca_res - prcomp(t(cpm_log[top_genes, ]), scale. TRUE) # 画前两个主成分 plot(pca_res$x[, 1], pca_res$x[, 2], col steelblue, pch 19, cex 0.6, xlab PC1, ylab PC2, main PCA of BRCA samples)逻辑说明scale. TRUE 让每个基因在 PCA 里权重一致否则高表达基因会主导主成分。取 2000 个高变基因是经验值太少会丢信号太多会把噪声带进来。如果 PCA 图里出现一两个远离主群的样本先查它们的测序深度和临床注释确认是技术离群还是真实的特殊亚型再决定剔不剔。3.2 一致性聚类为什么它比 k-means 更适合样本分型k-means 需要预先指定簇数而且对初始中心敏感。样本分型场景下我们往往不知道到底该分几类所以更常用的是 ConsensusClusterPlus 这类一致性聚类。它的思路是反复重采样看哪些样本稳定地分在一起最后用一致性矩阵和 delta area 图帮你判断簇数。library(ConsensusClusterPlus) # 输入要求是矩阵行是基因列是样本 cc - ConsensusClusterPlus( as.matrix(cpm_log[top_genes, ]), maxK 6, reps 1000, pItem 0.8, pFeature 1, clusterAlg hc, distance pearson, seed 1234, plot png )逻辑说明maxK 设 6 是给一个上限实际选几类看后面的图reps 是重采样次数1000 次在千级样本上耗时可以接受样本少时可以降到 500pItem 0.8 表示每次抽 80% 的样本pFeature 1 表示用全部基因clusterAlg 用层次聚类 “hc”distance 用 “pearson”这两个组合在表达谱上比较常见。seed 固定是为了结果可复现不固定的话每次跑出来的簇边界会有细微差别。3.3 怎么判断该分几类delta area 和一致性矩阵一起看ConsensusClusterPlus 会输出每个 k 对应的一致性矩阵图。判断簇数时我一般看两个东西一是 delta area 曲线拐点出现在哪个 k二是一致性矩阵里对角线上的块是否清晰。如果 k3 和 k4 的 delta area 差不多但 k4 的矩阵里有一块明显模糊那就选 k3。# 查看 delta area 曲线 plot(cc[[3]]$consensusTree, main Consensus tree for k3) # 提取 k3 的簇标签 cluster_labels - cc[[3]]$consensusClass table(cluster_labels)逻辑说明cc[[k]] 里存的是对应 k 的结果consensusClass 是每个样本的簇归属。table 一下看每类样本数如果某一类只有个位数样本要么是真实的小亚型要么是离群点被单独切开需要结合临床注释判断。这一步没有绝对标准delta area 的拐点只是参考最终要回到“这些簇在生物学上说得通吗”。4. 聚类结果的临床验证把簇标签和生存、分型对上4.1 用生存曲线看簇间预后差异聚类做完最直接的验证是 KM 生存曲线。如果几个簇的生存曲线几乎重合那这个分型在临床上没有区分度要么换特征基因要么换算法。TCGA-BRCA 的临床数据里生存时间和事件状态是现成的整理成 data.frame 后接 survival 和 survminer 就行。library(survival) library(survminer) # 假设 clinical_df 里有 OS.time、OS 和样本 ID clinical_df$cluster - cluster_labels[match(clinical_df$sampleID, names(cluster_labels))] fit - survfit(Surv(OS.time, OS) ~ cluster, data clinical_df) ggsurvplot(fit, pval TRUE, risk.table TRUE, palette c(#2E9FDF, #E7B800, #FC4E07))逻辑说明match 这一步是把聚类标签按样本 ID 对齐到临床表顺序错了会导致标签错位这是血泪经验里最常见的一类翻车。pval TRUE 会显示 log-rank 检验的 p 值但 p 值小不代表分型一定有用还要看曲线分开的幅度和风险表里的样本数。如果某一类在后期样本数掉到个位数曲线尾部的差异要谨慎解读。4.2 和 PAM50 分型做交叉表TCGA-BRCA 自带 PAM50 分型注释这是现成的外部验证。把自己的簇标签和 PAM50 做交叉表如果某个簇明显富集 Luminal A 或 Basal-like说明聚类抓到了已知的生物学信号如果完全随机那就要回头检查特征基因和归一化。# 交叉表 table(cluster_labels, clinical_df$PAM50_subtype) # 卡方检验看关联是否显著 chisq.test(table(cluster_labels, clinical_df$PAM50_subtype))逻辑说明交叉表看的是富集模式不是要求一一对应。聚类得到的簇数和 PAM50 的类数不一定相同一个簇里混了两种 PAM50 也是正常的关键看有没有明显的富集方向。卡方检验的 p 值受样本量影响很大千级样本下几乎必然显著所以更该看的是标准化残差或 Cramer‘s V而不是只盯 p 值。4.3 簇特征基因的提取与富集有了稳定的簇下一步是找每个簇的特征基因然后做 GO 或 KEGG 富集看这些簇在通路层面有没有区别。这一步能把“数学上的分组”翻译成“生物学上的解释”。library(limma) # 以某一簇 vs 其余簇为例做差异分析 design - model.matrix(~ 0 factor(cluster_labels)) colnames(design) - levels(factor(cluster_labels)) fit - lmFit(cpm_log, design) contrast - makeContrasts(C1 - (C2 C3)/2, levels design) fit2 - eBayes(contrasts.fit(fit, contrast)) top_genes - topTable(fit2, number 100, adjust.method BH)逻辑说明design 用 0 截距写法是为了方便构造任意对比makeContrasts 里 C1、C2、C3 要和簇标签的 level 对应写错顺序结果就反了。topTable 取前 100 个基因做富集adjust.method 用 BH 控制 FDR。富集结果里如果出现大量和细胞周期、免疫应答相关的通路说明簇间差异有生物学基础如果全是核糖体或线粒体基因多半是技术噪声。5. 避坑与排查聚类分析 TCGA-BRCA 时最容易翻车的 5 个点5.1 样本 ID 对不上簇标签整体错位现象生存曲线看起来有差异但交叉表和 PAM50 完全对不上富集结果也乱七八糟。原因TCGA 的样本 ID 有“-01A”“-01B”这类后缀表达矩阵的列名和临床表的 sampleID 格式不一致match 时大量返回 NA 或错配。解决统一截取前 15 位作为样本 ID再 matchmatch 后检查 sum(is.na(cluster_labels)) 是否为 0不为 0 就先处理缺失。5.2 用 FPKM 直接聚类结果被少数高表达基因主导现象PCA 图里样本挤成一条线聚类出来的簇和测序批次高度相关。原因FPKM 没有做跨样本的分布校正少数极高表达基因如线粒体基因在距离计算里权重过大。解决回到原始计数用 log-CPM 或 VST 归一化如果只有 FPKM至少做一次 log 变换和分位数归一化并在 PCA 前取高变基因。5.3 低表达过滤太松聚类结果不稳定现象每次换 seed 或换重采样次数簇边界变化很大。原因大量低表达基因在全零和个位数之间波动给距离计算贡献的是噪声。解决提高过滤阈值比如要求至少在 50% 样本中计数大于 10同时把高变基因数量从 2000 降到 1000 到 1500看一致性矩阵是否变清晰。5.4 簇数选得太多每类只剩几个样本现象delta area 在 k5 或 k6 时还在下降但选出来的簇里有一类只有三五个样本。原因一致性聚类在样本异质性高时会倾向于切出小簇这些小簇可能是离群点而不是真实亚型。解决结合临床注释看小簇样本是否有共同的批次或临床特征如果没有优先选更小的 k或者先把离群样本剔除再跑。5.5 把聚类结果当成诊断结论现象簇间生存差异显著就直接说“发现了新的 BRCA 亚型”。原因聚类是无监督的它只保证样本在表达谱上相似不保证这种相似有临床因果。解决至少做三件事——和已知分型交叉验证、在独立数据集里复现、用多变量 Cox 校正临床变量后再看簇是否独立预后。缺了这三步结论只能停在“观察到分组”不能叫“发现亚型”。6. 进阶技巧用 NMF 做软聚类和特征基因模块提取一致性聚类给的是硬标签每个样本只属于一个簇。但肿瘤样本往往处在几种状态的过渡带上硬切会丢掉这种连续性。非负矩阵分解NMF是这类场景下我比较常用的补充手段它把表达矩阵分解成“特征基因模块”和“样本对模块的权重”权重是连续的可以看每个样本在多个模块上的占比。library(NMF) # 只取高变基因NMF 对矩阵维度敏感 nmf_input - as.matrix(cpm_log[top_genes, ]) # 指定 rank 范围跑多次看稳定性 nmf_res - nmf(nmf_input, rank 2:5, nrun 30, seed 1234) # 看不同 rank 的 cophenetic 系数 plot(nmf_res) # 取最优 rank 的模块基因 best_rank - 3 w - basis(nmf_res[[best_rank]]) h - coef(nmf_res[[best_rank]])逻辑说明rank 是模块数nrun 是每个 rank 跑的次数30 次是为了让结果稳定cophenetic 系数越高说明分解越稳定但也不是越高越好还要看模块基因有没有生物学意义。basis 矩阵是基因乘模块coef 矩阵是模块乘样本后者可以拿来做样本的连续分型。和硬聚类相比NMF 的优势是能处理混合状态代价是解释起来更复杂适合已经跑通硬聚类、想再挖一层的人。我自己的习惯是先用一致性聚类拿到一个能对上临床的硬分型再用 NMF 看这些簇内部的连续性最后把两种结果叠在一起看哪些样本在边界上。这样做的好处是既有一个能拿去做生存分析的离散标签又不会把样本的异质性一刀切掉。聚类分析这件事算法只是工具真正决定结果值不值得信的是你对数据的反复检查和对自己结论的克制。希望帮到你。本文还有配套的精品资源点击获取