Mfuzz包最初是為了研究具有時間序列特征的轉(zhuǎn)錄組和蛋白組數(shù)據(jù)中基因或蛋白表達的時間趨勢的一個工具包，它可以將基因按時間序列表達模式進行聚類，可以直觀地看到相同表達模式的基因在不同時間點的表達變化趨勢。
安裝
BiocManager::install("Mfuzz")
library("Mfuzz")
data(yeast)
yeast@assayData$exprs
酵母基因表達矩陣存儲在yeast@assayData$exprs，行為基因，列為不同時間點。
缺失值的處理
yeast.r <- filter.NA(yeast, thres=0.25)
thres參數(shù)設(shè)定閾值，如果某個基因的缺失值（NA）的百分比大于該閾值，則排除該基因。
填補缺失值
yeast.f <- fill.NA(yeast.r,mode="mean")
上一步驟還遺留了一部分缺失值，用該基因在所有樣本中的平均值替代缺失值NA，還可以是median（中位數(shù)），knn和wknn。如果沒有缺失值可以跳過該步驟。
數(shù)據(jù)過濾
根據(jù)其標準差（standard deviation）過濾表達水平低或者波動小的數(shù)據(jù)
yeast.s <- filter.std(yeast.f,min.std=0)
標準化數(shù)據(jù)
yeast.s <- standardise(yeast.f)
此處標準化實際為歸一化，使每個基因/蛋白的平均表達值為零，標準差為1。 此外，還有standardise2函數(shù)可以選擇，此時的到的結(jié)果為 standardise函數(shù)的結(jié)果沿著y軸平移，使其起點為零。
m值評估
Mfuzz包聚類需要兩個參數(shù)，一個是m，另一個是聚類數(shù)c。先用函數(shù)計算m值。
m1 <- mestimate(yeast.s)
m1
#[1] 1.149465
評估C聚類簇數(shù)
然后評估c值。
tmp <- Dmin(yeast.s,m=m1,crange=seq(2,40,1),repeats=3,visu=TRUE)
tmp
#[1] 2.889000 2.828634 2.837718 2.683517 2.680458 2.571425 2.648475 2.573232
#[9] 2.487565 2.284549 2.234671 2.216442 2.208910 2.150243 2.191863 2.047843
#[17] 2.055428 2.042039 1.925507 1.904427 1.833664 1.885157 1.594121 #1.840121
#[25] 1.740238 1.770424 1.719631 1.670364 1.664545 1.760566 1.736006 1.780919
#[33] 1.687090 1.751894 1.752642 1.620106 1.718774 1.769282 1.706513
此時的x軸對應(yīng)的為c值（Cluster number），可以看到隨著聚類數(shù)的增加，centroid distance越來越小。
可以通過兼顧Cluster number和centroid distance都最小的原則肉眼選擇c值，也可以通過函數(shù)來看。
Cluster= seq(2,40,1)
s=which(tmp==min(tmp))
Cluster[s]
mfuzz_cluster <- mfuzz(mfuzz_class, c = cluster_num, m = mestimate(mfuzz_class))
mfuzz.plot(yeast.s,cl = mfuzz_cluster,mfrow=c(4,4),time.labels=seq(0,160,10))
參考資料：
使用Mfuzz包進行基因表達的時間趨勢分析并劃分聚類群_生信寶典的博客-CSDN博客
https://mp.weixin.qq.com/s/J5hNJlNj2-fSIOkFLlfhMg
