使用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 type1 GSM2467146 W22 GSM2467147 W23 GSM2467148 W24 GSM2467149 W25 GSM2467150 W26 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 W5TNFAIP8L1 6.372295 6.129076 5.944875 6.085894OTOP2 6.152929 5.683012 4.943879 5.070979SAMD7 4.215502 4.090231 3.668044 3.706270ARRDC5 6.612363 6.129062 5.665847 5.974294FAM86C1 7.558140 7.155209 6.769227 6.862116C2CD4B 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)

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