Mfuzz能够识别表达谱的潜在时间序列模式，并将相似模式的基因聚类，以帮助我们了解基因的动态模式和它们功能的联系。Mfuzz的核心算法是模糊c均值聚类分析，用于识别相似的基因表达谱。此外，Mfuzz提供了绘图功能，除了实现基因表达谱的聚类外，还能绘制时间序列，清晰地为我们呈现基因表达的动力学特征。
图片出自Immunity. 2020 Sep 15;53(3):685-696.e3.
Mfuzz包的使用
1. 安装R包，加载演示数据集
#Bioconductor 安装 Mfuzz 包
BiocManager::install('Mfuzz')
library(Mfuzz) #加载
data(yeast) #Mfuzz包内部数据集，酵母基因表达矩阵，行为基因，列为时间样本（按时间顺序来）
dat <- as.matrix(yeast@assayData$exprs)
View(dat) #查看表达矩阵
在示例的基因表达矩阵中，每一行是一种基因，列为样本。其中，第一列为第一个时间点的样本，第二列为第二个时间点的样本，以此类推。
2. 对基因表达矩阵进行标准化处理，并处理缺失值后，执行聚类分析，将具有相似的时间表达模式的基因聚在一类并绘图。
#构建对象
dat <- new('ExpressionSet',exprs = dat)
#处理 NA 值
dat <- filter.NA(dat, thres = 0.25)
# 49 genes excluded.
dat <- fill.NA(dat, mode = 'mean')
#根据标准差去除样本间差异太小的基因
dat <- filter.std(dat, min.std = 0)
#0 genes excluded.
#标准化
dat <- standardise(dat)
Standard deviation of gene expression vectors before standardisation.
#fuzzy c-means 聚类，需手动定义聚类个数，比方说设置 16 个簇
n <- 16
#评估出最佳的 m 值，防止随机数据聚类
m <- mestimate(dat)
m
# 1.15
#聚类
set.seed(2021)
cl <- mfuzz(dat, c = n, m = m)
#作图，time.labels 参数设置时间轴，和原基因表达数据集中的列对应
mfuzz.plot(dat, cl = cl, mfrow = c(4, 4), time.labels = seq(0, 160, 10))
如上过程基于基因表达值进行了聚类，对于每个簇中的基因，具有相似的时间表达特征。随后，即可从图中识别一些重要的聚类簇，比方说簇中基因随时间表达趋势增加或减少，以及在特定时间出现了更高或更低的表达等，以建立和观察的表型的联系。
3. 获取各簇中包含的基因集。
#每个簇下基因数量
cl$size
#每个基因所属簇
head(cl$cluster)
#基因和 cluster 之间的 membership，用于判断基因所属簇，对应最大值的那个簇
head(cl$membership)
#整合关系输出
gene_cluster <- cbind(cl$cluster, cl$membership)
colnames(gene_cluster)[1] <- 'cluster'
write.table(gene_cluster, 'gene_cluster.txt', sep = '\t', col.names = NA, quote = FALSE)
View(gene_cluster)
前两列就是基因名称和聚类簇的对应关系
这样，就将基因名称和其所属的聚类簇对应起来了。根据上文的折线图判断重要的时间表达模式的基因集，并在该表中进一步筛选出更具体的基因名称就可以了。
得到时间序列模式后，就可以对不同时间动力学特征基因的功能分析了。
网友评论