专栏首页生信技能树使用Mfuzz包做时间序列分析

使用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)

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

本文分享自微信公众号 - 生信技能树(biotrainee),作者:生信技能树

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-09-26

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 惊!3个同样的数据挖掘策略居然同时发表

    这个问题怎么说呢,生命科学领域的数据挖掘课题的发表主要是靠工作量,很少有新颖或者前沿,无非就是替换癌症替换分子替换生物学功能基因集,我整理过大家耳熟能详的策略,...

    生信技能树
  • 有趣的基因命名

    还有,如果你看到HS.开头的基因,它是unigene的ID了,已经不再是symbol啦。

    生信技能树
  • 关键问题答疑:WGCNA的输入矩阵到底是什么格式

    这样的问题我其实被问过好多次了,因为这次是学员提问,虽然已经过了一个月的答疑期,但是情谊还在,所以就系统性的回复一下。

    生信技能树
  • Mysql order by 优化

    本节描述MySQL何时可以使用索引来满足ORDER BY子句,当不能使用索引时使用filesort,以及优化器中有关ORDER BY的执行计划信息。

    XING辋
  • 【论文推荐】最新5篇视频分类相关论文—上下文门限、深度学习、时态特征、结构化标签推理、自动机器学习调优

    【导读】专知内容组整理了最近五篇视频分类(Video Classification)相关文章,为大家进行介绍,欢迎查看! 1.Learnable pooling...

    WZEARW
  • [日常] 小白来装机基本概念BIOS与硬盘分区

    这两天因为在linux进行测试,先是搞坏了linux的系统,然后在重装linux系统后搞坏了引导。在修复引导的过程中,搞坏了本机的win8系统,再次修复引导与重...

    陶士涵
  • tf.compat.v1.get_default_graph

    版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。 ...

    于小勇
  • Java自动化测试(接口操作优化 15)

    由于发起post请求的时候,它可能为json格式,也可能为form表单格式。所以对他进行提取

    zx钟
  • [抄录]Unicode & UTF

    Fundamentally, computers just deal with numbers. They store letters and other ch...

    雪飞鸿
  • Science发布最疯癫视频!AI守门员以最搞笑的方式让对手抓狂

    人工智能一直被视为学习能力极强、学习速度极快的“超级物种”,秒杀人类不解释、无商量;无论是在国际围棋比赛,还是在星际争霸游戏对决,均多次战胜人类顶级选手,所向披...

    新智元

扫码关注云+社区

领取腾讯云代金券