cell_6.jpg
跟着 Cell 学作图 | 6.时间序列分析(Mfuzz包)
“实践是检验真理的唯一标准。”
“复现是学习R语言的最好办法。”
2021.4.12_1
这篇2020年发表在cell上关于新冠的组学文章里面有大量的生信内容。今天带大家复现其中的一个Supplemental Figure：时间序列分析图。
20210425_1
22
本文代码及示例数据领取：后台回复“20210426”
时间序列分析
在研究基因表达谱或者蛋白表达谱时，经常会涉及到对时间序列的分析。例如，不同的基因或蛋白表达水平随时间表现出怎样的动力学特征，怎样挖掘潜在的时间特征？Mfuzz是用来进行不同时间点转录组数据表达模式聚类分析的R包。它能够识别表达谱的潜在时间序列模式，并将相似模式的基因聚类，以帮助我们了解基因的动态模式和它们功能的联系。本图中1,2,3,4,分别表示健康、非新冠肺炎、非重度新冠肺炎、重度新冠肺炎。
数据格式
210421_2
注：示例数据仅作展示用，无实际意义！
绘制
#------
title: "Mfuzz"
author: "MZBJ"
date: "2020/4/25"
#-----
setwd("F:/~/mzbj/cell/20210426")
df <- read.csv(file = "df_umap.csv")
df1 <- df[,-1]
#分组求均值
#aggregate用法上一篇原创推文已经介绍过了，不熟悉的可以去回顾一下
df2<-aggregate(df1[,colnames(df1)[2:ncol(df1)]],by=list(df1$label),mean,na.rm= TRUE)
#把第一列设为行名
row.names(df2)<-df6[,1]
df3<-data.frame(t(df2[,-1]))
#第一次使用要下载
BiocManager::install("Mfuzz")
library("Mfuzz")
#构建对象
df3a<-as.matrix(df3)
df3Ex<- ExpressionSet(assayData = df3a)
#排除了超过25%的测量缺失的基因
df3F <- filter.NA(df3Ex,thres = 0.25)
#用相应基因的平均值表达值替换剩余的缺失值
df3F <- fill.NA(df3F,mode = 'mean')
#标准化
df3F <- standardise(df3F)
#聚类
set.seed(2021)
#手动定义聚类个数 c
cl <- mfuzz(df3F,c=8,m=1.25)
#作图
pdf("mfuzz.pdf")
mfuzz.plot2(df3F, cl=cl,mfrow=c(4,4),centre=TRUE,x11=F,centre.lwd=0.2)
dev.off()
#批量导出每个聚类所包含的基因
dir.create(path="mfuzz",recursive = TRUE)
for(i in 1:8){
potname<-names(cl$cluster[unname(cl$cluster)==i])
write.csv(cl[[4]][potname,i],paste0("mfuzz","/mfuzz_",i,".csv"))
}
出图：
20210425_2
一模一样4个字我说累了哈哈哈~
写在后面：
本系列重在复现，所以有些细节可能讲的不是很详细。大家有问题可以后台私信，或者在我的B站：
木舟笔记进行互动！制作不易，谢谢大家多多支持！
往期内容：
跟着Cell学作图 | 2.柱状图+误差棒+散点+差异显著性检验
