使用Mfuzz包分析基因表達的時間趨勢并劃分聚類群

本篇簡介一個R包，Mfuzz（http://mfuzz.sysbiolab.eu）。Mfuzz包最初是為處理基因表達或蛋白表達譜數(shù)據(jù)而開發(fā)的一種聚類方法，核心算法基于模糊c均值聚類（Fuzzy C-Means Clustering，FCM），用于在具有時間序列特征的轉(zhuǎn)錄組、蛋白質(zhì)組數(shù)據(jù)中分析基因或蛋白表達的時間趨勢，并將具有相似表達模式的基因或蛋白劃分聚類，幫助了解這些生物學分子的動態(tài)模式以及與功能的聯(lián)系。

盡管Mfuzz包在一開始是為處理基因表達或蛋白表達譜數(shù)據(jù)而開發(fā)的，但實際應用中也可以對其它類型的生物學或非生物學數(shù)據(jù)進行聚類分析，或者“其它非時間梯度”的情形，這些在本篇的最后也有簡單提及。

本篇不涉及Mfuzz的詳細計算細節(jié)，主要簡介如何在R語言中使用Mfuzz包執(zhí)行聚類分析。

???

一篇使用到Mfuzz包聚類的相關文獻案例

??

首先來看一篇文獻的部分內(nèi)容，我當初也是在這篇文獻中第一次看到了使用Mfuzz包對時間序列劃分聚類群。

Gao等（2017）基于蛋白質(zhì)譜的方法，研究了小鼠胚胎著床前發(fā)育過程中的蛋白質(zhì)組。共涉及了6個發(fā)育階段，受精卵（Zygote）、二分胚（2-cell）、四分胚（4-cell）、八分胚（8-cell）、桑葚胚（Morula）和囊胚（Blastocyst）。為了將蛋白質(zhì)功能與胚胎發(fā)育相結(jié)合，作者首先表征了蛋白質(zhì)豐度與胚胎發(fā)育階段的時間關系，根據(jù)所有蛋白質(zhì)在每個階段的豐度信息，通過Mfuzz包對這些蛋白質(zhì)執(zhí)行了時間序列的聚類。最終獲得了10組聚類群（即10組蛋白群），代表了胚胎蛋白質(zhì)的10種動力學模式，并在隨后對這10組蛋白群的豐度變化與其功能展開了更細致的討論（如基因集富集分析，蛋白網(wǎng)絡分析等）。

圖片節(jié)選自Gao等（2017）原文。左側(cè)來自概要圖，展示了小鼠胚胎著床前發(fā)育的6個階段（受精卵、二分胚、四分胚、八分胚、桑葚胚和囊胚）的蛋白質(zhì)組試驗流程；右上側(cè)來自原文圖1C，為使用Mfuzz包對蛋白質(zhì)組的聚類分析，根據(jù)蛋白質(zhì)豐度隨胚胎發(fā)育階段的時間關系，獲得了10組不同動力學模式的蛋白群；右下側(cè)來自原文圖5C，聯(lián)合蛋白質(zhì)表達的時間模式和蛋白質(zhì)功能，對小鼠胚胎發(fā)育階段中蛋白質(zhì)組功能的概括。

當然，本篇不介紹那么多東西，只關注Mfuzz包的聚類分析那部分（上圖紅色框區(qū)域）。

???

使用Mfuzz包分析基因表達的時間趨勢并劃分聚類群的簡單演示

接下來，我們不妨就以上述Gao等（2017）的蛋白質(zhì)組數(shù)據(jù)為例，展示使用Mfuzz包對時間序列類型數(shù)據(jù)的聚類過程。

Gao等（2017）的蛋白質(zhì)表達矩陣表格可以在原文獻的補充材料Table S1中自行下載（Excel表格中，Sheet名稱為“union_all_protein_exp_cluster”）：

https://www.cell.com/cms/10.1016/j.celrep.2017.11.111/attachment/bc1eedf3-89d1-473f-9b66-b44a17994029/mmc2.xlsx

我在網(wǎng)盤中也保存了一份（已重命名為“mmc2.union_all_protein_exp.txt”，提取碼4mtu）：

https://pan.baidu.com/s/1wXYop3wql8SZqcM2CzpJIg

?

示例數(shù)據(jù)集概要

所示的蛋白質(zhì)表達矩陣文件內(nèi)容如下，第一列為蛋白質(zhì)名稱，隨后幾列依次為這些蛋白質(zhì)在小鼠受精卵（Zygote）、二分胚（2-cell）、四分胚（4-cell）、八分胚（8-cell）、桑葚胚（Morula）和囊胚（Blastocyst）時期的相對豐度數(shù)值。

?

使用Mfuzz包執(zhí)行時間序列的聚類分析

根據(jù)幫助文檔的操作過程，加載Mfuzz包后，將數(shù)據(jù)表讀取到R中，執(zhí)行數(shù)據(jù)轉(zhuǎn)換、標準化、聚類等一系列操作，將具有相似的時間表達特征的蛋白聚在一類。

#使用 Bioconductor 安裝 Mfuzz 包 #install.packages('BiocManager') #BiocManager::install('Mfuzz')#加載 Mfuzz 包 library(Mfuzz)#讀取矩陣表格，在我網(wǎng)盤中，示例數(shù)據(jù)為“mmc2.union_all_protein_exp.txt” #該示例中，行為基因或蛋白名稱，列為時間樣本（按時間順序提前排列好，若存在生物學重復需提前取均值） protein <- read.delim('mmc2.union_all_protein_exp.txt', row.names = 1, check.names = FALSE) protein <- as.matrix(protein)#構建對象 mfuzz_class <- new('ExpressionSet',exprs = protein)#預處理缺失值或者異常值 mfuzz_class <- filter.NA(mfuzz_class, thres = 0.25) mfuzz_class <- fill.NA(mfuzz_class, mode = 'mean') mfuzz_class <- filter.std(mfuzz_class, min.std = 0)#標準化數(shù)據(jù) mfuzz_class <- standardise(mfuzz_class)#Mfuzz 基于 fuzzy c-means 的算法進行聚類，詳情 ?mfuzz #需手動定義目標聚類群的個數(shù)，例如這里我們?yōu)榱酥噩F(xiàn)原作者的結(jié)果，設定為 10，即期望獲得 10 組聚類群 #需要設定隨機數(shù)種子，以避免再次運行時獲得不同的結(jié)果 set.seed(123) cluster_num <- 10 mfuzz_cluster <- mfuzz(mfuzz_class, c = cluster_num, m = mestimate(mfuzz_class))#作圖，詳情 ?mfuzz.plot2 #time.labels 參數(shù)設置時間軸，需要和原基因表達數(shù)據(jù)集中的列對應 #顏色、線寬、坐標軸、字體等細節(jié)也可以添加其他參數(shù)調(diào)整，此處略，詳見函數(shù)幫助 mfuzz.plot2(mfuzz_class, cl = mfuzz_cluster, mfrow = c(2, 5), time.labels = colnames(protein))

如上過程基于蛋白質(zhì)表達值的時間序列進行了聚類。根據(jù)預先指定的聚類數(shù)量，最終獲得了10組不同動力學模式的聚類群（蛋白群），如下圖所示。對于每個聚類群中的蛋白質(zhì)，它們具有相似的時間表達特征；而不同聚類群的蛋白質(zhì)之間的動力學模式則差異明顯。

由于直接使用的Gao等（2017）的蛋白質(zhì)組數(shù)據(jù)，我們期望重現(xiàn)原作者的分析，您可以將分析結(jié)果和原文獻進行比較，發(fā)現(xiàn)結(jié)果是基本一致的。聚類群的命名序號（cluster 1、2、3等）會存在區(qū)別，但這完全沒有影響，您可以重新匹配具有相似外觀的聚類群，肯定都可以找到完全相似的另一個，然后重新編號即可。極少數(shù)蛋白可能與原文獻所劃分的聚類群不完全一致，因為它們的時間特征比較模糊，而Mfuzz包實質(zhì)上基于模糊c均值聚類的算法，難以為它們鑒定準確的邊界，故極少數(shù)蛋白出現(xiàn)聚類不穩(wěn)定的情形。當然還可能是分析過程中，參數(shù)選擇和原作者不同等，導致的微小區(qū)別。

在獲得了聚類結(jié)果后，即可從圖中識別一些重要的或者感興趣的蛋白集合，比方說某些聚類群的蛋白質(zhì)出現(xiàn)了預期的隨時間增加而增加或減少的趨勢，在特定時間點出現(xiàn)了相對更高或更低的表達，或者觀察到明顯的拐點等。并繼續(xù)對這些感興趣的蛋白質(zhì)進行功能分析（如基因集富集分析，蛋白網(wǎng)絡分析等），以及建立和細胞或生物體的表型特征的聯(lián)系等，討論它們的生物學意義。

?

當然，討論蛋白質(zhì)的功能不是本篇的內(nèi)容，后續(xù)的分析需要做哪些，您自己根據(jù)實際情況來。在這之前，一個有待解決的問題是，如何獲得各聚類群中，都包含哪些蛋白呢？

接下來繼續(xù)在上述已獲得的聚類結(jié)果中，提取10個聚類群中包含的蛋白質(zhì)集合。

#查看每個聚類群中各自包含的蛋白數(shù)量 cluster_size <- mfuzz_cluster$size names(cluster_size) <- 1:cluster_num cluster_size#查看每個蛋白所屬的聚類群 head(mfuzz_cluster$cluster)#Mfuzz 通過計算一個叫 membership 的統(tǒng)計量判斷蛋白質(zhì)所屬的聚類群，以最大的 membership 值為準 #查看各蛋白的 membership 值 head(mfuzz_cluster$membership)

#最后，提取所有蛋白所屬的聚類群，并和它們的原始表達值整合在一起 protein_cluster <- mfuzz_cluster$cluster protein_cluster <- cbind(protein[names(protein_cluster), ], protein_cluster) head(protein_cluster) write.table(protein_cluster, 'protein_cluster.txt', sep = '\t', col.names = NA, quote = FALSE)#如果您想提取數(shù)據(jù)分析過程中，標準化后的表達值（繪制曲線圖用的那個值，而非原始蛋白表達值） protein_cluster <- mfuzz_cluster$cluster protein_standard <- mfuzz_class@assayData$exprs protein_standard_cluster <- cbind(protein_standard[names(protein_cluster), ], protein_cluster) head(protein_standard_cluster) #write.table(protein_standard_cluster, 'protein_standard_cluster.txt', sep = '\t', col.names = NA, quote = FALSE)

這樣，就將蛋白名稱、蛋白表達值以及其所屬的聚類群對應起來了。如果根據(jù)上文的折線圖挑選出了感興趣的時間表達特征的聚類群，就可以在該表中進一步將這些聚類群中的蛋白質(zhì)信息提取出來。再往后，分析這些蛋白的功能等，不再多說。

???

其它一些常見問題

我應該劃分多少個聚類群？

可以手動逐一嘗試去評估。如果您看到有些聚類群的折線圖在外觀上比較相似，它們大致反映了同一種時間趨勢時，表明存在冗余的聚類群，可能意味著需減少聚類群數(shù)量。如果您看到有些聚類群的折線圖內(nèi)的折線并不整齊，波動幅度過大，或者出現(xiàn)了多個密度中心，可能意味著需增加聚類群數(shù)量。多次嘗試后，再結(jié)合實際情況，挑選出感覺最合適的那個出來。

有一些機器學習方法，可以幫助自動評估最優(yōu)的聚類群數(shù)量。例如在前文“k均值劃分聚類”中，曾簡單提到過一些，如NbClust包的NbClust()、vegan包的cascadeKM()等。Mfuzz的內(nèi)部原理是模糊c均值聚類，它實際上和k均值劃分聚類等比較相像，因此這些方法也可以直接借鑒。此外還有很多方法可以實現(xiàn)這一點，大家有興趣可自行查閱相關資料了解。

存在組內(nèi)生物學重復時怎么處理？

以上示例數(shù)據(jù)中，每個時間點都只有一列數(shù)據(jù)。如果您的數(shù)據(jù)中包含生物學重復樣本，也就是一個時間點對應多列數(shù)據(jù)時，需要提前將生物學重復樣本進行合并，例如取均值等。聚類函數(shù)mfuzz()的幫助文檔里也是這樣建議的。

其它類型的數(shù)據(jù)能用Mfuzz包分析嗎？

當然可以。以上主要是以蛋白質(zhì)組數(shù)據(jù)為例的展示，如果您換成轉(zhuǎn)錄組、微生物組、代謝組，或者其它非生物學數(shù)據(jù)，也是可行的，只不過可能在數(shù)據(jù)預處理或標準化方式上需要變更。并且，如果不是時間序列，而是其它類型的“梯度”的數(shù)據(jù)，如不同藥物處理濃度下基因表達數(shù)據(jù)、不同環(huán)境梯度下的物種豐度數(shù)據(jù)，這些情況下也存在一種“梯度序列”，理論上也都可以嘗試用Mfuzz包進行聚類。

參考文獻

Gao Yawei, Liu Xiaoyu, Tang Bin et al. Protein Expression Landscape of Mouse Embryos during Pre-implantation Development. Cell Rep, 2017, 21: 3957-3969.

往期精品(點擊圖片直達文字對應教程)

機器學習

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結(jié)

以上是生活随笔為你收集整理的使用Mfuzz包进行基因表达的时间趋势分析并划分聚类群的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。