跟著 Cell 學(xué)作圖 | 6.時(shí)間序列分析(Mfuzz包)
“實(shí)踐是檢驗(yàn)真理的唯一標(biāo)準(zhǔn)堤器。”
“復(fù)現(xiàn)是學(xué)習(xí)R語(yǔ)言的最好辦法÷骞茫”
這篇2020年發(fā)表在cell上關(guān)于新冠的組學(xué)文章里面有大量的生信內(nèi)容。今天帶大家復(fù)現(xiàn)其中的一個(gè)Supplemental Figure：時(shí)間序列分析圖皮服。
22
本文代碼及示例數(shù)據(jù)領(lǐng)壤惆：后臺(tái)回復(fù)“20210426”
時(shí)間序列分析
在研究基因表達(dá)譜或者蛋白表達(dá)譜時(shí)，經(jīng)常會(huì)涉及到對(duì)時(shí)間序列的分析冰更。例如产徊，不同的基因或蛋白表達(dá)水平隨時(shí)間表現(xiàn)出怎樣的動(dòng)力學(xué)特征，怎樣挖掘潛在的時(shí)間特征蜀细？Mfuzz是用來(lái)進(jìn)行不同時(shí)間點(diǎn)轉(zhuǎn)錄組數(shù)據(jù)表達(dá)模式聚類(lèi)分析的R包舟铜。它能夠識(shí)別表達(dá)譜的潛在時(shí)間序列模式，并將相似模式的基因聚類(lèi)奠衔，以幫助我們了解基因的動(dòng)態(tài)模式和它們功能的聯(lián)系谆刨。本圖中1,2,3,4,分別表示健康、非新冠肺炎归斤、非重度新冠肺炎痊夭、重度新冠肺炎。
數(shù)據(jù)格式
注：示例數(shù)據(jù)僅作展示用脏里，無(wú)實(shí)際意義她我！
繪制
#------
title: "Mfuzz"
author: "MZBJ"
date: "2020/4/25"
#-----
setwd("F:/~/mzbj/cell/20210426")
df <- read.csv(file = "df_umap.csv")
df1 <- df[,-1]
#分組求均值
#aggregate用法上一篇原創(chuàng)推文已經(jīng)介紹過(guò)了，不熟悉的可以去回顧一下
df2<-aggregate(df1[,colnames(df1)[2:ncol(df1)]],by=list(df1$label),mean,na.rm= TRUE)
#把第一列設(shè)為行名
row.names(df2)<-df6[,1]
df3<-data.frame(t(df2[,-1]))
#第一次使用要下載
BiocManager::install("Mfuzz")
library("Mfuzz")
#構(gòu)建對(duì)象
df3a<-as.matrix(df3)
df3Ex<- ExpressionSet(assayData = df3a)
#排除了超過(guò)25%的測(cè)量缺失的基因
df3F <- filter.NA(df3Ex,thres = 0.25)
#用相應(yīng)基因的平均值表達(dá)值替換剩余的缺失值
df3F <- fill.NA(df3F,mode = 'mean')
#標(biāo)準(zhǔn)化
df3F <- standardise(df3F)
#聚類(lèi)
set.seed(2021)
#手動(dòng)定義聚類(lèi)個(gè)數(shù) 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()
#批量導(dǎo)出每個(gè)聚類(lèi)所包含的基因
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"))
}
出圖：
一模一樣4個(gè)字我說(shuō)累了哈哈哈~
寫(xiě)在后面：
本系列重在復(fù)現(xiàn)，所以有些細(xì)節(jié)可能講的不是很詳細(xì)番舆。大家有問(wèn)題可以后臺(tái)私信酝碳，或者在我的B站：
木舟筆記進(jìn)行互動(dòng)！制作不易恨狈，謝謝大家多多支持疏哗！
往期內(nèi)容：
跟著Cell學(xué)作圖 | 2.柱狀圖+誤差棒+散點(diǎn)+差異顯著性檢驗(yàn)
