github-actions[bot]
commented
2 years ago 代码分享| 使用Mfuzz聚类并进行基因表达趋势（时间序列）分析 by 生信益站
❝上篇文章提到了 Mfuzz 包，今天讲讲具体用法。
安装与加载
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("Mfuzz")
library(Mfuzz)
以差异表达基因聚类为例
首先读取差异基因表达 FPKM 矩阵
setwd("zzzzzzzzzz")
fpkm <- read.delim('degs.mean.fpkm.txt', row.names = 1, check.names = FALSE)
head(fpkm)
## 0d 10d 20d 30d
## gene_1 29.96 32.39 52.44 101.96
## gene_2 3.99 3.43 3.37 41.93
## gene_3 7.67 20.23 15.55 13.11
## gene_4 8.07 3.70 3.79 3.81
## gene_5 75.17 34.35 33.90 39.71
## gene_6 2.94 23.37 17.76 5.08
fpkm <- as.matrix(fpkm)
构建 mfuzz 对象并处理缺失及异常值
mfuzz_class <- new('ExpressionSet',exprs = fpkm)
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
set.seed(123)
cluster_num <- 10
mfuzz_cluster <- mfuzz(mfuzz_class, c = cluster_num, m = mestimate(mfuzz_class))
# 作图，详情 ?mfuzz.plot2
mfuzz.plot2(mfuzz_class, cl = mfuzz_cluster, mfrow = c(2, 5), time.labels = colnames(fpkm),xlab="")
查看聚类结果
#查看每个cl中的差异基因数量
cluster_size <- mfuzz_cluster$size
names(cluster_size) <- 1:cluster_num
cluster_size
## 1 2 3 4 5 6 7 8 9 10
## 498 801 604 646 448 676 1368 1184 743 571
#查看每个差异基因所属的cl
head(mfuzz_cluster$cluster)
## gene_1 gene_2 gene_3 gene_4 gene_5 gene_6
## 6 7 9 8 8 5
# Mfuzz 通过计算一个叫 membership 的统计量判断差异基因所属的聚类群，以最大的 membership 值为准
# 查看每个基因的 membership 值
head(mfuzz_cluster$membership)
## 1 2 3 4 5 6 7 8
## gene_1 0.016405195 0.01591740 0.025620318 0.023235058 0.022382779 0.57943002 0.18482359 0.017602650
## gene_2 0.003006479 0.00300517 0.003780348 0.004721877 0.003621172 0.03528505 0.92989849 0.003390942
## gene_3 0.087444202 0.04345911 0.075781878 0.036200686 0.261996328 0.05656512 0.04985371 0.037292965
## gene_4 0.023386340 0.07690291 0.016253066 0.099633345 0.015229423 0.01605005 0.01719577 0.704469536
## gene_5 0.037055966 0.10374028 0.026542572 0.247029213 0.024909633 0.02817521 0.03064016 0.450116898
## gene_6 0.099953270 0.04205534 0.073382183 0.030859526 0.364997418 0.03892693 0.03499994 0.034025999
整合表达信息与聚类结果并输出
# 最后，提取所有基因所属的聚类群，并与fpkm整合
gene_cluster <- mfuzz_cluster$cluster
gene_cluster <- cbind(fpkm[names(gene_cluster), ], gene_cluster)
head(gene_cluster)
## 0d 10d 20d 30d gene_cluster
## gene_1 29.96 32.39 52.44 101.96 6
## gene_2 3.99 3.43 3.37 41.93 7
## gene_3 7.67 20.23 15.55 13.11 9
## gene_4 8.07 3.70 3.79 3.81 8
## gene_5 75.17 34.35 33.90 39.71 8
write.table(gene_cluster, 'gene_cluster.xls', sep = '\t', col.names = NA, quote = FALSE)
Mfuzz 算法优势
Mfuzz 采用的模糊均值算法（Fuzzy C-Means Clustering，FCM）是普通 K-means 算法的改进。普通 K-means 算法对于数据的划分是硬性的，让它聚成 10 个类，结果就会出现 10 个类。Mfuzz 则是一种柔性的模糊划分，让它聚成 10 个类，结果可能是 6 个类，也可能是 15 个类。K-means 的硬聚类把每个待辨识的对象严格地划分到某类中，具有非此即彼的性质。Mfuzz 的模糊聚类建立了样本对类别的不确定性描述，会根据数据本身的特性来调节聚类结果，从而更有优势。
聚类后可以做什么
思路在往期文章中提到过：
绘图专题| Mfuzz表达模式聚类和networkD3可视化案例
组学服务广而告之
联系站长
❝欢迎来 QQ 群 837738558 进行交流，群内可直接进行提问、唠嗑、学习、资料下载。
更多优质内容请扫码或点击下方名片，关注“生信益站”公众号。
https://mp.weixin.qq.com/s/nI_lKqwHUTeTvpm5cHFAsA