本篇简介一个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中，执行数据转换、标准化、聚类等一系列操作，将具有相似的时间表达特征的蛋白聚在一类。
代码语言：javascript复制#使用 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个聚类群中包含的蛋白质集合。
代码语言：javascript复制#查看每个聚类群中各自包含的蛋白数量
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.
