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

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

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

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

???

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

??

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

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

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

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

???

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

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

Gao等（2017）的蛋白質表達矩陣表格可以在原文獻的補充材料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

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

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

?

示例數據集概要

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

?

使用Mfuzz包執行時間序列的聚類分析

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

#使用 Bioconductor 安裝 Mfuzz 包 #install.packages('BiocManager') #BiocManager::install('Mfuzz')#加載 Mfuzz 包 library(Mfuzz)#讀取矩陣表格，在我網盤中，示例數據為“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)#標準化數據 mfuzz_class <- standardise(mfuzz_class)#Mfuzz 基于 fuzzy c-means 的算法進行聚類，詳情 ?mfuzz #需手動定義目標聚類群的個數，例如這里我們為了重現原作者的結果，設定為 10，即期望獲得 10 組聚類群 #需要設定隨機數種子，以避免再次運行時獲得不同的結果 set.seed(123) cluster_num <- 10 mfuzz_cluster <- mfuzz(mfuzz_class, c = cluster_num, m = mestimate(mfuzz_class))#作圖，詳情 ?mfuzz.plot2 #time.labels 參數設置時間軸，需要和原基因表達數據集中的列對應 #顏色、線寬、坐標軸、字體等細節也可以添加其他參數調整，此處略，詳見函數幫助 mfuzz.plot2(mfuzz_class, cl = mfuzz_cluster, mfrow = c(2, 5), time.labels = colnames(protein))

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

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

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

?

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

接下來繼續在上述已獲得的聚類結果中，提取10個聚類群中包含的蛋白質集合。

#查看每個聚類群中各自包含的蛋白數量 cluster_size <- mfuzz_cluster$size names(cluster_size) <- 1:cluster_num cluster_size#查看每個蛋白所屬的聚類群 head(mfuzz_cluster$cluster)#Mfuzz 通過計算一個叫 membership 的統計量判斷蛋白質所屬的聚類群，以最大的 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)#如果您想提取數據分析過程中，標準化后的表達值（繪制曲線圖用的那個值，而非原始蛋白表達值） 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)

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

???

其它一些常見問題

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

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

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

存在組內生物學重復時怎么處理？

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

其它類型的數據能用Mfuzz包分析嗎？

當然可以。以上主要是以蛋白質組數據為例的展示，如果您換成轉錄組、微生物組、代謝組，或者其它非生物學數據，也是可行的，只不過可能在數據預處理或標準化方式上需要變更。并且，如果不是時間序列，而是其它類型的“梯度”的數據，如不同藥物處理濃度下基因表達數據、不同環境梯度下的物種豐度數據，這些情況下也存在一種“梯度序列”，理論上也都可以嘗試用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.

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

機器學習

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。