使用Mfuzz包做时间序列分析

生信入门学员群有一个提问，关于做时间序列分析的， 一般的提问，我们都会让大家去搜索我们的公众号历史教程，搜索是有技巧的，见：

尤其是《生信技能树》公众号的号内搜索，居家旅行必备神器，好的关键词肯定能解决大家的问题。

但是时间序列分析我确实没有写过，我们的转录组授课讲师张娟看到了学员提问，里面写了个教程投稿！

下面是《张娟》的分享

既然是讲解时间序列分析，那么就不得不提一下Mfuzz包了，恰好生信技能树创始人jimmy的200篇生物信息学文献阅读活动分享过的一篇文章就有这个（2018年的文献列表在：https://zhuanlan.zhihu.com/c_1024966446748618752）

该文献题目：Dynamic network biomarker indicates pulmonary metastasis at the tipping point of hepatocellular carcinoma

杂志：NATURE COMMUNICATIONS | (2018) 9:678

文章中的图如下：

 

作者主要使用了第一个结果中差异表达分析得到的13,247 个差异基因列表（使用的是传统的T检验，对任意两组的组合找差异，最后合并）。

1 数据下载

表达谱数据：文章中的测序数据上传到了GEO：GSE94016

差异基因列表：文章的附件41467_2018_3024_MOESM4_ESM.xlsx，网页版文章中可以直接下载。

2 数据预处理

我们在GEO下载下来的数据是探针水平的，那么首先要将探针水平的表达谱处理成基因水平的，代码如下：

# 清空当前环境变量
rm(list=ls())
options(stringsAsFactors = F)

# 读取表达谱
probe_exp <- read.table("../data/GSE94016_series_matrix.txt",comment.char = "!",strip.white=T,sep="\t",header = T)
probe_exp[1:4,1:4]

# 读取平台注释文件
anno <- read.table("../data/GPL15207-17536.txt",header = T, comment.char = "#",sep="\t",quote = "\"",strip.white = T,fill = T)
colnames(anno)

anno1 <- anno[,c("ID","Gene.Symbol")]
head(anno1)

expdata <- merge(anno1,probe_exp,by.x = "ID",by.y = "ID_REF")

# 去掉一对多的探针
expdata <- expdata[-(grep("///",expdata$Gene.Symbol)),]

# 去掉一对空的探针
expdata <- expdata[-which(expdata$Gene.Symbol==""),]

# 对多个探针注释到一个基因上的取均值
# 最后剩下18836个基因
library(limma)
expdata1 <- limma::avereps(expdata[,-c(1:2)],ID=expdata[,2])
expdata1[1:4,1:4]
dim(expdata1)

# 保存数据
save(expdata1,file = "../data/expdata1.RData")

然后根据差异基因列表得到差异表达基因，在这里还发现，有一些基因在我前面的探针水平处理到基因水平的时候，丢失了一些差异表达基因，可能是由于预处理方式跟作者不太一样，不过对结果的影响不会很大。代码如下：

## 读取作者的差异表达基因列表
DEGs <- read.table("../data/DEG.txt",header = T,sep="\t") 
head(DEGs)

# 提取差异表达基因的表达
loc <- match(DEGs$Symbols,rownames(expdata1))
loc <- loc[!is.na(loc)]
DEGs_exp <- expdata1[loc,]

看文章中的图，我们发现横坐标是时间节点，那么我们根据样本的时间节点信息，需要将差异基因表达谱处理一下，变成时间节点的表达，时间节点信息来自GEO的matrix文件的表头注释。

# 读入样本时间节点
time <- read.table("../data/sample_type.txt",header = T,sep="\t")
head(time)
  geo_accession type
1    GSM2467146   W2
2    GSM2467147   W2
3    GSM2467148   W2
4    GSM2467149   W2
5    GSM2467150   W2
6    GSM2467151   W3

tmp <- data.frame(colnames(DEGs_exp),t(DEGs_exp))
temp <- data.frame(time[match(tmp[,1],time[,1]),],tmp)
temp[1:5,1:4]
DEGs_exp_averp <- t(limma::avereps(temp[,-c(1:3)],ID=temp[,2]))

# 最后的时间节点表达谱如下
head(DEGs_exp_averp)
                W2       W3       W4       W5
TNFAIP8L1 6.372295 6.129076 5.944875 6.085894
OTOP2     6.152929 5.683012 4.943879 5.070979
SAMD7     4.215502 4.090231 3.668044 3.706270
ARRDC5    6.612363 6.129062 5.665847 5.974294
FAM86C1   7.558140 7.155209 6.769227 6.862116
C2CD4B    3.326672 3.163326 2.835682 2.858878

这样就挑选到了作者分析好的13,247 个差异基因列表的表达矩阵啦！

3 时序分析

文章中使用的是R包Mfuzz，这个包是时序分析最常用的，如下所示：

Mfuzz采用了一种新的聚类算法fuzzy c-means algorithm。在文献中称这种聚类算法为soft clustering算法，相比K-means等hard clustering算法，一定程度上降低了噪声对聚类结果的干扰，而且这种算法有效的定义了基因和cluster之间的关系，即基因是否属于某个cluster, 对应的值为memebership。

对于分析而言，我们只需要提供基因表达量的数据就可以了，需要注意的是，Mfuzz默认你提供的数据是归一化之后的表达量，这意味着表达量必须可以直接在样本间进行比较，对于FPKM, TPM这两种定量方式而言，是可以直接在样本间进行比较的；但是对于count的定量结果，我们必须先进行归一化，可以使用edgeR或者DESeq先得到归一化之后的数据在进行后续分析。

我们得到的GEO中的表达谱是经过了MAS5.0处理的affy的芯片数据，正好可以直接使用。

通过以下几个步骤就可以得到聚类的结果。

# 首先安装R包并加载
BiocManager::install("Mfuzz")
library(Mfuzz)

## 1. 预处理：去除表达量太低或者在不同时间点间变化太小的基因等步骤
# Mfuzz聚类时要求是一个ExpressionSet类型的对象，所以需要先用表达量构建这样一个对象。
eset <- new("ExpressionSet",exprs = DEGs_exp_averp)

# 根据标准差去除样本间差异太小的基因
eset <- filter.std(eset,min.std=0)


## 2. 标准化：聚类时需要用一个数值来表征不同基因间的距离，Mfuzz中采用的是欧式距离，
# 由于普通欧式距离的定义没有考虑不同维度间量纲的不同，所以需要先进行标准化
eset <- standardise(eset)


## 3. 聚类：Mfuzz中的聚类算法需要提供两个参数，第一个参数为希望最终得到的聚类的个数，这个参数由我们直接指定；
# 第二个参数称之为fuzzifier值，用小写字母m表示，可以通过函数评估一个最佳取值
c <- 6 # 聚类个数，文章中用的6个聚类，我们也用6个
m <- mestimate(eset) #  评估出最佳的m值
cl <- mfuzz(eset, c = c, m = m) # 聚类

# 在cl这个对象中就保存了聚类的完整结果，对于这个对象的常见操作如下
cl$size # 查看每个cluster中的基因个数，看的出来与文章每个类别基因个数差了一些
[1] 2269 1982 2289 1653 1393 1729

cl$cluster[cl$cluster == 1] # 提取某个cluster下的基因
cl$membership # 查看基因和cluster之间的membership

## 4. 可视化
library(RColorBrewer)
color.2 <- colorRampPalette(rev(c("#ff0000", "Yellow", "OliveDrab1")))(1000)
mfuzz.plot(eset,cl,mfrow=c(2,3),new.window= FALSE,time.labels=colnames(DEGs_exp_averp),colo = color.2)

最后复现出来的图如下，可以看的出名字虽然不同，每个类别的基因个数也差了一点点，但是趋势基本是一样的。

 

如果你完全看不懂我们教程里面的代码，时间和精力不够，自学起来困难重重（虽然有那么多资料），而且又确实有经费支持你参加学习班，不妨考虑我们生信技能树官方举办的学习班：

• 数据挖掘学习班第7期（线上直播3周，马拉松式陪伴，带你入门），原价4800的数据挖掘全套课程， 疫情期间半价即可抢购。
• 生信爆款入门-第9期（线上直播4周，马拉松式陪伴，带你入门），原价9600的生信入门全套课程，疫情期间3.3折即可抢购。