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