簡介這份資源面向生物信息學入門學習者與R語言數(shù)據(jù)分析實踐者圍繞TCGA-BRCA乳腺癌基因表達數(shù)據(jù)展開聚類分析教學案例。內(nèi)容涵蓋利用層次聚類按基因表達水平對病人分型、選擇average距離度量、繪制heatmap可視化聚類效果并實現(xiàn)PCA降維、依據(jù)主成分數(shù)目重新聚類再結(jié)合臨床數(shù)據(jù)中的ER_Status_nature2012標簽評估聚類是否符合預期幫助讀者理解無監(jiān)督學習在腫瘤分型中的應用邏輯。資源包共22個文件約10.91MB包含8個png與8個pdf圖示、2個md說明、2個txt數(shù)據(jù)文件、1個license及1個R源碼腳本圖示與文檔便于對照結(jié)果與復現(xiàn)流程。目前已有687人學習下載適合需要完整賽題式分析方案、可運行代碼與可視化輸出參考的讀者。1. 從一份 TCGA-BRCA 聚類作業(yè)說起R 軟件怎么把基因表達矩陣變成分型圖如果你手頭有一份基因表達矩陣行是基因、列是病人導師讓你“按表達水平把病人分個類”你大概率會先愣一下從哪下手這份「生物信息學概論——聚類分析TCGA-BRCA數(shù)據(jù).zip」就是干這個的。它把一套完整的 R 分析流程、示例數(shù)據(jù)、出圖腳本和結(jié)果圖打包在一起核心任務是用層次聚類和 PCA 對 TCGA-BRCA 的基因表達數(shù)據(jù)做無監(jiān)督分型再用臨床里的 ER 狀態(tài)做外部驗證。適合正在做生物信息課程設計、想跑通“表達矩陣 → 聚類 → 熱圖 → PCA 降維 → 再聚類”這條鏈的人。我拆完的第一感覺是代碼不長但坑全在數(shù)據(jù)預處理和參數(shù)選擇上下面按能復現(xiàn)的順序講。2. 數(shù)據(jù)拆包與矩陣對齊GeneMatrix 和 clinical_data 怎么接上拿到壓縮包先別急著跑cluster.R把文件結(jié)構看清楚。包里有GeneMatrix.txt、clinical_data.txt、cluster.R、README.md以及Figures-PNG和Figures-PDF兩個結(jié)果目錄。圖里包含cluster、pca_cluster、pca_screeplot、heatmap、ER_heatmap、pca_cumulative、pca_heatmap、ER_cluster這幾類基本覆蓋了作業(yè)要求的全部輸出。2.1 兩個核心文件的格式與讀入方式GeneMatrix.txt是表達矩陣一行一個基因一列一個病人帶行名和列名。clinical_data.txt是臨床信息一行一個病人病人編號與表達矩陣的列名對應但只覆蓋其中一部分。這意味著你不能直接merge完事得先取交集。# 讀入表達矩陣check.namesFALSE 防止列名里的連字符被改掉 expr - read.table(GeneMatrix.txt, header TRUE, row.names 1, sep \t, check.names FALSE) # 讀入臨床信息 clin - read.table(clinical_data.txt, header TRUE, row.names 1, sep \t, check.names FALSE) # 取病人交集順序以表達矩陣列為準 common_samples - intersect(colnames(expr), rownames(clin)) expr - expr[, common_samples] clin - clin[common_samples, , drop FALSE] # 確認沒有全零或全 NA 的基因 expr - expr[rowSums(is.na(expr)) 0, ] expr - expr[rowSums(expr) 0, ]邏輯說明check.names FALSE很關鍵TCGA 的樣本編號常帶-R 默認會把連字符轉(zhuǎn)成點號后面按編號取子集就會對不上。取交集后必須讓clin的行順序跟expr的列順序一致否則后面 ER 狀態(tài)和聚類標簽會錯位這種錯位不會報錯但熱圖上的注釋條會全亂屬于典型的“玄學翻車”。參數(shù)上sep \t對應制表符分隔如果你的文件是逗號分隔就改sep ,。row.names 1表示第一列是行名若第一列是數(shù)字索引而非基因名需要改成row.names NULL再手動賦值。2.2 表達值分布檢查與標準化選擇聚類對量綱敏感。如果矩陣里既有原始 counts 又有 TPM或者不同樣本測序深度差異大直接算距離會被高表達基因主導。常見做法是先做log2(x 1)轉(zhuǎn)換再按行做 z-score 標準化。# 判斷是否需要 log 轉(zhuǎn)換若最大值遠大于 50 且分布右偏先 log summary(as.numeric(expr[1, ])) # log2 轉(zhuǎn)換加 1 避免 log(0) expr_log - log2(expr 1) # 按基因行做 z-score讓每個基因在樣本間均值為 0、標準差為 1 expr_scaled - t(scale(t(expr_log))) # 檢查是否有 NaN標準差為 0 的基因會產(chǎn)生 NaN sum(is.nan(expr_scaled)) expr_scaled - expr_scaled[!is.nan(rowSums(expr_scaled)), ]這里有個選擇理由作業(yè)要求“按基因表達水平聚類病人”距離用 average即 UPGMA。如果不對基因做標準化高表達基因會主導距離矩陣聚類結(jié)果基本只反映那幾個高變基因。z-score 之后每個基因權重接近聚類更能反映整體表達模式。注意scale默認按列所以要用t()轉(zhuǎn)置兩次這是 R 里最容易寫反的一步。提示如果sum(is.nan(expr_scaled))不為 0說明有基因在所有樣本中表達值完全相同標準化后標準差為 0必須剔除否則dist()會報 NA 錯誤。3. 層次聚類實操average 距離、樹切分與熱圖輸出數(shù)據(jù)對齊后進入核心分析。作業(yè)明確要求用層次聚類、average 距離并輸出 heatmap。這一章把hclust、cutree和pheatmap的配合講清楚。3.1 距離矩陣與 hclust 參數(shù)# 計算樣本間距離用歐氏距離 dist_mat - dist(t(expr_scaled), method euclidean) # 層次聚類average 即 UPGMA hc - hclust(dist_mat, method average) # 查看聚類樹 plot(hc, labels FALSE, main Hierarchical Clustering (average))邏輯說明dist作用于t(expr_scaled)因為我們要聚類的是病人列不是基因行。method average對應作業(yè)要求它的合并準則是對兩類間所有點對距離取平均比 single 更穩(wěn)健比 complete 更不容易被離群點拉偏。hclust的method還可選ward.D2如果你想讓類內(nèi)方差最小化ward 通常更緊湊但作業(yè)指定 average 就按 average 來。參數(shù)上dist的method可選euclidean、manhattan、correlation。對表達數(shù)據(jù)有人偏好correlation因為它關注表達模式形狀而非絕對幅度。但作業(yè)沒指定用歐氏距離配合 z-score 是常見組合。3.2 切樹數(shù)目與熱圖注釋切幾類作業(yè)沒給死數(shù)但臨床里有 ER 狀態(tài)可以拿它做參照。常見做法是看樹高突變點或者直接切 2 到 4 類對比。# 切成 2 類 clusters - cutree(hc, k 2) # 與 ER 狀態(tài)做交叉表 table(clusters, clin$ER_Status_nature2012) # 用 pheatmap 出圖帶 ER 注釋條 library(pheatmap) annotation_col - data.frame(ER_Status clin$ER_Status_nature2012, Cluster factor(clusters)) rownames(annotation_col) - colnames(expr_scaled) pheatmap(expr_scaled, annotation_col annotation_col, show_rownames FALSE, show_colnames FALSE, clustering_distance_cols euclidean, clustering_method average, filename Figures-PNG/heatmap.png)邏輯說明cutree的k是切分簇數(shù)也可以給h按樹高切。交叉表能直觀看出聚類標簽和 ER 狀態(tài)是否吻合。pheatmap的annotation_col必須是有行名的數(shù)據(jù)框行名與矩陣列名一致否則注釋條對不上。clustering_distance_cols和clustering_method要與前面hclust保持一致不然熱圖上的樹和你的聚類結(jié)果不是一回事。參數(shù)上show_rownames FALSE是因為基因數(shù)通常上萬顯示會糊成一片。filename直接輸出 PNG省去手動保存。包里的ER_heatmap.PNG應該就是加了 ER 注釋的版本heatmap.PNG是不帶注釋的版本。注意pheatmap默認會對行做聚類如果你只想展示樣本聚類把cluster_rows設為FALSE否則行樹會占用大量繪圖空間。4. PCA 降維與再聚類主成分數(shù)目怎么定、和層次聚類比什么作業(yè)的第二問是 PCA 降維后再聚類并與第一問比較。這一步的坑在于主成分數(shù)目選擇選少了信息不夠選多了等于沒降維。4.1 PCA 實現(xiàn)與碎石圖、累積方差圖# 轉(zhuǎn)置后做 PCAprcomp 默認中心化 pca_res - prcomp(t(expr_scaled), center TRUE, scale. FALSE) # 碎石圖 png(Figures-PNG/pca_screeplot.png, width 800, height 600) plot(pca_res, type l, main Scree Plot) dev.off() # 累積方差貢獻 var_explained - summary(pca_res)$importance[3, ] png(Figures-PNG/pca_cumulative.png, width 800, height 600) plot(var_explained, type b, xlab PC, ylab Cumulative Variance, main Cumulative Variance Explained) abline(h 0.8, col red, lty 2) dev.off()邏輯說明prcomp的scale. FALSE是因為前面已經(jīng)做過 z-score再 scale 會重復標準化。碎石圖看拐點累積方差圖看達到 80% 需要幾個主成分。常見做法是選累積方差到 80% 或碎石圖拐點處的主成分數(shù)。包里的pca_screeplot.PNG和pca_cumulative.PNG就是這兩張圖。參數(shù)上center TRUE是必須的PCA 對中心化敏感。scale.如果設為TRUE會對每個基因再除標準差若前面沒做 z-score這里設TRUE也能達到類似效果但兩者疊加會過度標準化。4.2 用前幾個主成分再聚類并對比# 假設選前 3 個主成分根據(jù)累積方差圖調(diào)整 pc_use - 3 pca_data - pca_res$x[, 1:pc_use] # 對主成分得分做層次聚類 dist_pca - dist(pca_data, method euclidean) hc_pca - hclust(dist_pca, method average) clusters_pca - cutree(hc_pca, k 2) # 與原始聚類比較 table(Original clusters, PCA clusters_pca) # PCA 聚類熱圖 pheatmap(t(pca_data), annotation_col data.frame(Cluster factor(clusters_pca)), show_rownames TRUE, show_colnames FALSE, filename Figures-PNG/pca_heatmap.png)邏輯說明pca_res$x是樣本在主成分上的得分取前pc_use列作為新特征。再聚類后與原始clusters做交叉表看一致性。如果 PCA 聚類結(jié)果與原始層次聚類高度一致說明前幾個主成分已經(jīng)抓住了主要變異如果差異大可能是主成分選少了或者數(shù)據(jù)本身沒有明顯分群。參數(shù)上pc_use的選擇要有依據(jù)不能拍腦袋。我一般會看累積方差圖選達到 80% 的最小主成分數(shù)同時參考碎石圖拐點。包里的pca_cluster.PNG應該是 PCA 得分的散點圖按聚類標簽著色。提示PCA 聚類和原始聚類比較時簇的編號是任意的交叉表要允許標簽置換不能只看對角線。5. 避坑與排查從數(shù)據(jù)錯位到聚類不收斂的 5 個血淚經(jīng)驗這一章記錄我復現(xiàn)這類作業(yè)時踩過的坑每條按現(xiàn)象、原因、解決寫?，F(xiàn)象一熱圖注釋條顏色和聚類樹完全對不上。原因clin取子集后行順序沒跟expr列順序?qū)R或者annotation_col的行名與矩陣列名不一致。解決取交集后強制clin - clin[colnames(expr), , drop FALSE]并在pheatmap前用all(rownames(annotation_col) colnames(expr_scaled))檢查返回FALSE就停下來查。現(xiàn)象二hclust報錯“NA/NaN/Inf in foreign function call”。原因表達矩陣里有 NA、Inf或者 z-score 后出現(xiàn) NaN。解決讀入后先expr - expr[complete.cases(expr), ]再檢查sum(is.nan(expr_scaled))把標準差為 0 的基因剔除。別用na.omit直接刪樣本會丟病人?，F(xiàn)象三聚類結(jié)果每次跑都不一樣。原因?qū)哟尉垲惐旧硎谴_定性的但如果前面用了隨機初始化比如 k-means或者數(shù)據(jù)里有重復列導致距離矩陣不穩(wěn)定。解決確認用的是hclust而非kmeans檢查duplicated(t(expr_scaled))有重復樣本就合并或剔除。現(xiàn)象四PCA 累積方差圖顯示前兩個主成分就到 95%但聚類效果很差。原因數(shù)據(jù)沒做標準化少數(shù)高表達基因主導了方差前幾個主成分只反映那幾個基因。解決回到 2.2 節(jié)做 z-score再重新 PCA。標準化前后各跑一次對比累積方差曲線你會看到明顯差別?，F(xiàn)象五ER 狀態(tài)交叉表里 NA 一大堆。原因clinical_data.txt里 ER 狀態(tài)列有缺失值或者病人編號匹配時因為連字符/點號問題沒對上。解決讀入時check.names FALSE匹配前用gsub(\\., -, colnames(expr))統(tǒng)一格式對 ER 缺失的樣本在交叉表里單獨列為 NA不要直接刪否則樣本量對不上。6. 進階技巧用輪廓系數(shù)和 ER 一致性給聚類結(jié)果打分跑完圖不算完作業(yè)要求“評估聚類是否符合預期”。除了看 ER 交叉表還可以算輪廓系數(shù)和調(diào)整蘭德指數(shù)。輪廓系數(shù)衡量簇內(nèi)緊密度和簇間分離度值越接近 1 越好調(diào)整蘭德指數(shù)衡量聚類標簽與 ER 狀態(tài)的一致性越接近 1 越吻合。library(cluster) library(mclust) # 輪廓系數(shù)需要距離矩陣 sil - silhouette(clusters, dist_mat) png(Figures-PNG/silhouette.png, width 800, height 600) plot(sil, main Silhouette Plot) dev.off() mean(sil[, 3]) # 調(diào)整蘭德指數(shù)比較聚類與 ER 狀態(tài) eri - adjustedRandIndex(clusters, clin$ER_Status_nature2012) print(paste(Adjusted Rand Index:, round(eri, 3)))邏輯說明silhouette的第一個參數(shù)是簇標簽第二個是距離矩陣必須與聚類時用的距離一致。mean(sil[, 3])是平均輪廓系數(shù)大于 0.5 說明結(jié)構清晰0.2 到 0.5 算一般低于 0.2 就要懷疑是否真的存在分群。adjustedRandIndex對標簽置換不敏感比直接看交叉表對角線更客觀。參數(shù)上silhouette對大規(guī)模數(shù)據(jù)計算慢如果病人數(shù)上千可以先抽樣。mclust包提供adjustedRandIndex如果不想裝額外包也可以手寫交叉表算 Rand 指數(shù)但調(diào)整版更嚴謹。我一般會把這幾個指標和 ER 交叉表放一起看如果輪廓系數(shù)低但 ARI 高說明聚類結(jié)構弱但恰好和 ER 吻合如果兩者都低說明這批數(shù)據(jù)在表達層面沒有明顯分型或者預處理有問題。從那以后我每次跑完聚類都強制走一遍“輪廓系數(shù) ARI 交叉表”三件套不再只看熱圖顏色好不好看。希望幫到你。本文還有配套的精品資源點擊獲取