适用条件和背景
Mfuzz包最初是为了研究具有时间序列特征的转录组和蛋白组数据中基因或蛋白表达的时间趋势的一个工具包，能将具有相似表达模式的基因或蛋白进行软聚类，从而可以了解不同组别的基因集或蛋白集在生物学上的动态模式与功能的联系。简而言之，Mfuzz包可以对一个表达矩阵的基因进行聚类，然后可以看到不同组别的基因集或蛋白集在具有时间维度的样本上的表达动态变化。事实上，这个包可以不限于仅具有时间特征的样本，只要具有一定顺序特征的样本，理论上都可以使用，例如，具有空间顺序特征的样本。
安装Mfuzz
使用BiocManager::install安装
BiocManager::install('Mfuzz')
library(Mfuzz)
主函数plot_mfuzz
既然这个包是一开始是用于转录组数据的，那么也能用于单细胞转录组数据。分析思路，这里用的平均表达量矩阵作为输入，然后进行标准化，使用函数mfuzz进行软聚类，mfuzz.plot和mfuzz.plot2是自带的可视化函数。我写了一个函数，以下是函数参数介绍：
- obj, Seurat 对象；
-assay, 计算平均表达量的assay，默认为 'RNA'； - mode='knn', 填补缺失值的模式，默认为‘knn’，还可以选‘mean’和‘wknn’等
- nclusters, 可指定的聚类数目，默认为NULL，若为NULL则根据有多少个样本组别就聚成多少个clusters
- prefix, 输出文件前缀，默认为'out'
- gene.use, 指定聚类的基因集，默认为NULL，将使用高变基因
- nfeatures, FindVariableFeatures函数指定的基因数，默认为3000
- group.by, 细胞组别，默认为'ident'
plot_mfuzz <- function(obj,assay='RNA',mode='knn',nclusters=NULL,prefix='out',gene.use=NULL,nfeatures=3000,group.by='ident'){
library(Mfuzz)
if (is.null(gene.use)) {
obj <- FindVariableFeatures(obj,nfeatures=nfeatures)
gene.use <- VariableFeatures(obj)
}
ave <- AverageExpression(obj,group.by = group.by,assays=c(assay),features=gene.use)
ingene <- ave$RNA
gene_tpm <- data.matrix(ingene)
eset <- new("ExpressionSet",exprs = gene_tpm)
gene.r <- filter.NA(eset, thres=0.25)
#gene.f <- fill.NA(gene.r,mode="mean")
gene.f <- fill.NA(gene.r,mode=mode)
#gene.f <- fill.NA(gene.r,mode="wknn")
tmp <- filter.std(gene.f,min.std=0)
gene.s <- standardise(tmp)
if (is.null(nclusters)) {
nclusters=length(colnames(gene_tpm))
}
print(paste0('nclusters is:',nclusters))
c <- nclusters
m <- mestimate(gene.s)
cl <- mfuzz(gene.s, c = c, m = m)
saveRDS(cl,paste0(prefix,'.rds'))
write.table(cl$cluster,"output.txt",quote=F,row.names=T,col.names=F,sep="\t")
nrow <- ceiling(nclusters/3)
pdf(paste0(prefix,'.pdf'),18,ceiling(nclusters/6)*10)
mfuzz.plot(gene.s,cl,mfrow=c(nrow,3),new.window= FALSE)
#mfuzz.plot(gene.s,cl,mfrow=c(2,3),new.window= FALSE,time.labels=as.vector(colnames(gene_tpm)))
mfuzz.plot2(gene.s,cl,mfrow=c(nrow,3),time.labels=colnames(gene_tpm),x11 = FALSE)
dev.off()
}
使用示例
这里使用pbmc_3k数据集
library(Seurat)
library(SeuratData)
data("pbmc3k")
obj <- pbmc3k.final
运行上述函数：
plot_mfuzz(obj,assay='RNA',mode='knn',nclusters=NULL,prefix='out',gene.use=NULL,nfeatures=3000,group.by='ident')
#简化版
plot_mfuzz(obj)
结果说明
一共会生成3个文件：
1.output.txt
左边是基因名，右边是这个基因所属的类别，因为使用的组别把细胞分成了9组，所以默认把基因分成9类：
PPBP 8
LYZ 1
S100A9 1
IGLL5 6
GNLY 5
FTL 4
PF4 8
FTH1 4
GNG11 8
S100A8 1
FCER1A 7
HLA-DRA 7
CD74 7
CLU 8
NKG7 5
GZMB 5
CST3 7
CCL4 5
C1QA 4
HLA-DPB1 7
2.out.pdf
里面有两张图，一张是使用默认的x轴标签：
图1
另一张是自定义的标签：
图2
3.out.rds
mfuzz函数运行得到的对象，可以使用readRDS读入：
cl <- readRDS('out.rds')
#查看每类基因数目
cl$size
# 查看每类基因ID
cl$cluster[cl$cluster == 1]
总结与讨论
以上就是Mfuzz在单细胞转录组数据的应用，以此为模板，应该可以用到很多其它场景。在使用的过程中发现一个问题，自定义x周坐标时，标签显示不完整，目前找不到办法，可能只能手动调了，如果有大神知道解决办法欢迎评论区留言。
########参考博客
https://www.yisu.com/zixun/658715.html
https://cloud.tencent.com/developer/article/1845577
