Mfuzz 工具采用模糊 c 均值聚类（Fuzzy C-Means Clustering，FCM）算法，根据具有时间序列特征的转录组、蛋白质组数据中基因或蛋白表达的时间趋势，对具有相同表达模式的基因或蛋白划分 cluster。
相较于传统的统计学显著的上下调基因分组，Mfuzz 可以实现更为多元的生物学实验设计的基因聚类，如时间梯度或者浓度梯度处理的实验。
迄今为止的聚类大多为硬聚类（K-means 等算法），即将每个基因或蛋白质都完全分配给一个聚类，但在实际情况中，基因/蛋白质簇经?；岢鱿种氐?，且硬聚类算法通常对噪声十分敏感。为了克服硬聚类的局限性，Mfuzz 采用了软聚类的方法，不仅具有了更强的抗噪性，且可以避免对基因进行先验的预过滤，防止了生物学上相关的基因/蛋白质的丢失。
官网：http://193.136.227.155/sysbiolab/mfuzz/
安装
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("Mfuzz")
使用
1. 输入数据
Mfuzz 默认提供的数据是归一化之后的表达量，表达量必须可以直接在样本间进行比较，对于 FPKM 和 TPM 可直接使用；对于 count 需要进行归一化之后使用，可以使用 edgeR 或者 DESeq 得到归一化后的数据
library(Mfuzz)
## Mfuzz 聚类要求为 ExpressionSet 类型对象
# 载入测试数据，酵母基因表达矩阵，行为基因，列为时间样本（按时间顺序）
# 后续分析皆在测试数据上完成，真实数据同理
data(yeast) # yeast 为 ExpressionSet 对象
exprs_data <- as.matrix(yeast@assayData$exprs) # 酵母基因表达矩阵
# 载入真实数据（示例）
exprs_data <- read.delim('mmc2.union_all_protein_exp.txt', row.names = 1, check.names = FALSE)
exprs_data <- as.matrix(exprs_data) # 将数据框转换为数值矩阵
mfuzz_data <- new('ExpressionSet',exprs = exprs_data) # 构建 ExpressionSet 对象
补充：此处提供了一种归一化 count 的方法
# BiocManager::install("edgeR")
library(edgeR)
library(limma)
dat <- log2(edgeR::cpm(exprs_data)+1) # 注意 exprs_data 必须非 NA 且非负
# 此处 colnames(dat) 应为样本名，同一样本不同生物学重复列名相同
avereps_df <- t(limma::avereps(t(dat) , ID = colnames(dat))) # 取平均值
# 按时间顺序对列排序
avereps_df = avereps_df[,c( "time0", "10min","30min" ,
"1h" , "2h", "4h" , "8h" , "12h" )]
# 此处 avereps_df 即可用于构建 ExpressionSet 对象和后续分析
2. 数据预处理
2.1 处理NA值
# 排除丢失超过25％的测量值的基因，此处缺失值应在基因表达矩阵中用 NA 表示
yeast.r <- filter.NA(yeast, thres=0.25)
# 模糊c均值不允许缺失值，此处需要用相应基因的平均值表达值替换剩余缺失值
# 此处也可以使用（加权）k 最近邻法 （mode='knn'/'wknn'），性能通常较好，但计算量大
yeast.f <- fill.NA(yeast.r, mode="mean")
2.2 过滤低表达或变化小的基因
迄今大多数聚类需要根据标准差的阈值去除低水平表达或仅显示表达微小变化的基因，但由于过滤阈值的值设置没有严格的标准，Mfuzz 避免了对基因数据进行任何事先过滤，以防止可能具有生物学重要性的基因的丢失
# 根据标准差去除样本间差异太小的基因
tmp <- filter.std(yeast.f, min.std=0)
2.3 标准化
聚类时需要用一个数值来表征不同基因间的距离，Mfuzz中采用的是欧式距离；由于普通欧式距离的定义没有考虑不同维度间量纲的不同，所以需要先进行标准化
# 在欧几里得空间中进行聚类，因此基因的表达值为标准化后的平均值为0，标准差为1
yeast.s <- standardise(yeast.f)
注意：进行归一化是为了使不同的样本具有可比性，而标准化（在 Mfuzz 中）是为了使转录本（或基因或蛋白质）具有可比性
3. Mfuzz 聚类
3.1 FCM 聚类参数设置
对于 FCM 聚类，模糊参数 m 和聚类数 c 必须预先选择，因此 Mfuzz 需要提供两个参数：
- fuzzifier 值：用小写字母 m 表示，通过设置该值可以防止随机数据聚类，此处可通过函数评估最佳取值；可以使用函数
partcoef进行测试，是否针对m的特定设置对随机数据进行了聚类 - 希望最终得到的聚类的个数：用小写字母 c 表示，由用户直接指定；设置聚类 c 的最优数量通常很困难，尤其是对于短时间序列和重叠聚类的情况；此处可以通过
cselection测试 c 的范围，此过程会出现空聚类（merbership < 0.5）；此外聚类质心之间的最小距离Dmin可以用作聚类有效性指数，D.min 在达到最优 c 后下降得更慢；另外可以对一系列簇数进行聚类，然后对其生物学相关性进行评估选择最佳簇数，如GO分析
## m 值设置
m <- mestimate(yeast.s)
## c 值设置
# cselection 函数，repeats 为重复聚类，visu=TRUE 意味着输出可视化
tmp <- cselection(yeast.s, m=1.25, crange=seq(5,40,5), repeats=5, visu=TRUE) # 此处为出现空聚类的示例，由此可确定 c 的范围
# Dmin 函数确定最优 c, D.min 在达到最优 c 后下降得更慢
tmp <- Dmin(yeast.s, m=m, crange=seq(4,40,4), repeats=5, visu=TRUE)
c <- 12
3.2 聚类
cl <- mfuzz(yeast.s, c=c, m=m)
cl <- mfuzz(yeast.s, c=16, m=1.25)
## 可视化
# 默认可视化
mfuzz.plot(yeast.s, cl=cl, mfrow=c(4,4), time.labels=, time.labels=seq(0,160,10))
# 自定义颜色
library(RColorBrewer)
color <- colorRampPalette(rev(c("#ff0000", "Yellow", "OliveDrab1")))(1000)
mfuzz.plot(yeast.s, cl=cl, mfrow=c(4,4), time.labels=, time.labels=seq(0,160,10), colo = color)
黄色或绿色的线条对应于 membership 值较低的基因；红色和紫色线对应具有 membership 值较高的基因
3.3 查看聚类结果
# 在 cl 中保存了聚类的完整结果，对于这个对象的常见操作如下
cl$size # 查看每个cluster中的基因个数
cl$cluster[cl$cluster == 1] # 提取某个cluster下的基因
cl$membership # 查看基因和cluster之间的membership
4. 补充
4.1 membership 值
membership 值指示向量彼此之间的相似性，如果两个基因对于特定 cluster 具有较高的 membership 值，则它们通常彼此相似，此处定义 membership 值大于选定阈值的基因属于聚类的 α-core ，要提取属于聚类核心的基因列表，可以使用 acore 函数
acore.list <- acore(yeast.s, cl=cl, min.acore=0.7)
4.2 全局聚类
聚类之间的重叠为软聚类的特征之一，其表明两个簇之间具有共享的基因。重叠度低的 cluster 总体上表现出明显的差异模式；重叠度大的 cluster 则表现出更为相似的聚类模式，因此重叠度可以衡量成对 cluster 相似性，可以通过 overlap 和 overlap.plot 函数查看全局聚类的重叠情况
O <- overlap(cl)
Ptmp <- overlap.plot(cl,over=O,thres=0.05)
参考教程：
使用Mfuzz包做时间序列分析-腾讯云开发者社区-腾讯云 (tencent.com)
Mfuzz做转录变化的时间趋势分析后对每个趋势分组挑一个代表性基因-腾讯云开发者社区-腾讯云 (tencent.com)
