ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

TCGA-BRCA聚类分析R源码包:从表达矩阵到ER评估全流程

TCGA-BRCA聚类分析R源码包:从表达矩阵到ER评估全流程 简介这份资源面向生物信息学入门学习者与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先把文件按用途分三类后面排错会快很多。文件/目录类型用途GeneMatrix.txt输入数据基因表达矩阵行是基因列是病人clinical_data.txt输入数据病人临床信息含ER_Status_nature2012cluster.R源码主分析脚本聚类、PCA、绘图都在这里README.md说明运行顺序和依赖提示Figures-PNG/输出8 张 PNG 图快速预览用Figures-PDF/输出8 张 PDF 图写报告/论文用LICENSE授权使用前确认许可范围图目录里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.namesFALSE 防止列名里的连字符被改掉 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 的选择可以结合树状图高度和临床预期。切完之后用热图看分群是否在表达层面清晰。# 按 k2 切分对应 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 里最常用的实现默认对变量做中心化。# PCAscale.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()交集不为空第二距离度量先保持euclideanaverage跑通后再试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 图。从那以后我每次换数据都强制走一遍这个检查顺序省掉了很多“图好看但结论错”的后悔药。希望帮到你。本文还有配套的精品资源点击获取
返回列表