使用Mfuzz包做时间序列分析
下面是《张娟》的分享
既然是讲解时间序列分析,那么就不得不提一下Mfuzz包了,恰好生信技能树创始人jimmy的200篇生物信息学文献阅读活动分享过的一篇文章就有这个,作者主要使用了第一个结果中差异表达分析得到的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)
最后复现出来的图如下,可以看的出名字虽然不同,每个类别的基因个数也差了一点点,但是趋势基本是一样的。
- 洛谷P2197 nim游戏(Nim游戏)
- SQL Server 索引和表体系结构(聚集索引+非聚集索引)
- 3384/1750: [Usaco2004 Nov]Apple Catching 接苹果
- 1702: [Usaco2007 Mar]Gold Balanced Lineup 平衡的队列
- 1455: 罗马游戏
- SQL Server 高性能写入的一些总结
- 3389: [Usaco2004 Dec]Cleaning Shifts安排值班
- 1754: [Usaco2005 qua]Bull Math
- 3377: [Usaco2004 Open]The Cow Lineup 奶牛序列
- 3301: [USACO2011 Feb] Cow Line
- SQL Server 索引和表体系结构(包含列索引)
- TiDB 源码阅读系列文章(七)基于规则的优化
- 博弈论入门之nim游戏
- 1477: 青蛙的约会
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- Android——RecyclerView自定义OnScrollListener实现下拉刷新监听,上拉加载更多功能
- Java——类集框架:Map集合的详解及应用举例
- Java——类集框架:Set集合接口的详解及应用举例
- 如何有效地进行代码 Review?
- Java——对象序列化
- Android——MPAndroidChart折线图/柱状图/饼形图的使用
- Java——对String类型的时间进行加减操作
- Java——枚举基础应用总结(多例设计模式、Enum类、枚举的实际应用)
- Java——Annotation注解基本总结(简介、覆写、过期声明、压制警告)
- JavaWeb——一文快速入门BootStrap(栅格系统、全局CSS样式、组件、插件、基于BootStrap的官网案例实战)
- JavaWeb——XML入门详解(概述、语法、约束、Jsoup解析、Xpath解析)
- JavaWeb——CSS应用实例详解(概述、语法、选择器、属性、用户登录界面实例)
- JavaWeb——JavaScript精讲之事件监听机制与表单校验案例实战
- Java——扩展概念(匿名内部类、包装类、装箱与拆箱、数据类型的转换)
- Java——接口的基本总结(基本定义、使用接口定义标准、工厂设计模式、代理设计模式、抽象类与接口的区别)