Mfuzz包最初是为了研究具有时间序列特征的转录组和蛋白组数据中基因或蛋白表达的时间趋势的一个工具包，它可以将基因按时间序列表达模式进行聚类，可以直观地看到相同表达模式的基因在不同时间点的表达变化趋势。
安装
BiocManager::install("Mfuzz")
library("Mfuzz")
data(yeast)
yeast@assayData$exprs
酵母基因表达矩阵存储在yeast@assayData$exprs，行为基因，列为不同时间点。
缺失值的处理
yeast.r <- filter.NA(yeast, thres=0.25)
thres参数设定阈值，如果某个基因的缺失值（NA）的百分比大于该阈值，则排除该基因。
填补缺失值
yeast.f <- fill.NA(yeast.r,mode="mean")
上一步骤还遗留了一部分缺失值，用该基因在所有样本中的平均值替代缺失值NA，还可以是median（中位数），knn和wknn。如果没有缺失值可以跳过该步骤。
数据过滤
根据其标准差（standard deviation）过滤表达水平低或者波动小的数据
yeast.s <- filter.std(yeast.f,min.std=0)
标准化数据
yeast.s <- standardise(yeast.f)
此处标准化实际为归一化，使每个基因/蛋白的平均表达值为零，标准差为1。 此外，还有standardise2函数可以选择，此时的到的结果为 standardise函数的结果沿着y轴平移，使其起点为零。
m值评估
Mfuzz包聚类需要两个参数，一个是m，另一个是聚类数c。先用函数计算m值。
m1 <- mestimate(yeast.s)
m1
#[1] 1.149465
评估C聚类簇数
然后评估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
Cluster number
此时的x轴对应的为c值（Cluster number），可以看到随着聚类数的增加，centroid distance越来越小。
可以通过兼顾Cluster number和centroid distance都最小的原则肉眼选择c值，也可以通过函数来看。
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))
Cluster1.png
参考资料：
使用Mfuzz包进行基因表达的时间趋势分析并划分聚类群_生信宝典的博客-CSDN博客
https://mp.weixin.qq.com/s/J5hNJlNj2-fSIOkFLlfhMg
