时间序列基因表达数据的聚类 （Mfuzz包）

今天学习这篇文章

1.png
作者使用定量质谱，从每个阶段的 8,000 个小鼠胚胎（合子、2 细胞、4 细胞、8 细胞、桑椹胚和胚泡）中鉴定了近 5,000 种蛋白质。 发现受精卵、桑椹胚和胚泡中的蛋白质表达不同于 2 到 8 细胞胚胎。 蛋白质磷酸化分析确定了关键激酶和信号转导途径。 作者强调关键因素及其在胚胎发育中的重要作用。 转录组和蛋白质组数据的综合分析揭示了 RNA 降解、转录和翻译的协调控制，并确定了以前未定义的外显子连接衍生肽。
今天就来复现文章中的这个图
2.png
数据格式如下：
图片.png
数据为文章中的补充信息中table2。
链接：https://www.cell.com/cell-reports/fulltext/S2211-1247(17)31795-3?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS2211124717317953%3Fshowall%3Dtrue#supplementaryMaterial
图片.png
我们首先看看Mfuzz包的数据格式，找到包自带的数据看一下格式：
图片.png
能够看到每个时间点都只有一列表达值，所以本文中的数据胚胎的每个发育时期只需要一列表达值，我们需要用均值进行Mfuzz分析。
首先导入数据

#我们只需要一列均值在35列到40列，所以将首列蛋白名称+35-40列的数据另存为mmc2.txt
protein <- read.table("mmc2.txt",sep = "\t",header = T)
##需要转换成矩阵然后再进行数据的预处理和画图
protein <- as.matrix(protein)

#安装Mfuzz包
#使用 Bioconductor 安装 Mfuzz 包
install.packages('BiocManager')
BiocManager::install('Mfuzz')
#加载 Mfuzz 包
library(Mfuzz)

#构建对象
mfuzz <- new('ExpressionSet',exprs = protein)

数据处理以及数据的标准化

#预处理缺失值或者异常值
mfuzz_class <- filter.NA(mfuzz)
mfuzz_class <- fill.NA(mfuzz_class, mode = 'mean')
mfuzz_class <- filter.std(mfuzz_class, min.std = 0)

#标准化数据
mfuzz_class <- standardise(mfuzz_class)
# display of cluster cores with alpha = 0.5
#mfuzz.plot(mfuzz_class ,cl=cl,mfrow=c(2,2),min.mem=0.5)

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
#mfuzz.plot2(eset,cl,mfrow=c(1,1),colo,min.mem=0,time.labels,time.points,
ylim.set=c(0,0), xlab="Time",ylab="Expression changes",x11=TRUE,
                        ax.col="black",bg = "white",col.axis="black",col.lab="black",
                        col.main="black",col.sub="black",col="black",centre=FALSE,
                        centre.col="black",centre.lwd=2,
                        Xwidth=5,Xheight=5,single=FALSE,...)
#time.labels 参数设置时间轴，需要和原基因表达数据集中的列对应
#颜色、线宽、坐标轴、字体等细节也可以添加其他参数调整，此处略，详见函数帮助
mfuzz.plot2(mfuzz_class, cl = mfuzz_cluster, mfrow = c(2, 5), time.labels = colnames(protein))

查看每个聚类群中各自包含的蛋白数量

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)

原始值raw_data

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)

将数据导出每个cluster的数据，方便看每一个cluster中有相同随着时间变化趋势相同的基因或者蛋白

write.table(protein_standard_cluster, 'protein_standard_cluster.txt', sep = '\t', col.names = NA, quote = FALSE)