rnaseq_基因表达的时间序列分析（mfuzz）

mfuzz是用来进行不同时间点转录组数据表达模式聚类分析的R包。它能够识别表达谱的潜在时间序列模式，并将相似模式的基因聚类，以帮助我们了解基因的动态模式和它们功能的联系。

设置工作路径以及加载包

setwd('C:\\Users\\DI\\Desktop\\mfuzz')
library(Mfuzz)
library(RColorBrewer)
library(openxlsx)
library(Deseq2)

读取count矩阵并对其进行标准化

count <- read.csv("C:\\Users\\DI\\Desktop\\PRJNA516151_lnc\\lncRNA_Gastrocnemius_transcript_matrix.csv",header=T,row.names=1)
###### VST标准化##########
## 过滤在所有重复样本中小于1的基因，表达量太低也没研究意义
count <- count[rowMeans(count)>1,]
##载入样本信息
data <- read.table("C:\\Users\\DI\\Desktop\\PRJNA516151_lnc\\sample_data.txt",header = T,row.names = 1)
#一定要变为因子数据，否者用DESeq2包分析时候会出错
data[,1] <- as.factor(data$Type)
all(rownames(data) %in% colnames(count))
all(rownames(data) == colnames(count))
dds <-  DESeqDataSetFromMatrix(countData = count,colData = data,design = ~ Type)
dim(dds)
#过滤
dds <- dds[rowSums(counts(dds)) > 1,]
nrow(dds) 
#方差稳定变换，The variance stabilizing transformation
vsd <- vst(object=dds,blind=T) 
head(assay(vsd),3)
rs <- assay(vsd)
colnames(rs) <- gsub("X", "", colnames(rs))
colnames(rs) <- gsub(".{4}$", "", colnames(rs))
rs <- as.data.frame(rs)
rs <- as.data.frame( # sapply returns a list here, so we convert it to a data.frame
  sapply(unique(names(rs)), # for each unique column name
         function(col) rowMeans(rs[names(rs) == col]) # calculate row means
  )
)
d_1 <- rs[c(1:5)]
d_2 <- rs[c(6,7)]
rs <-  cbind(d_2,d_1)
rs <- as.matrix(rs)

mfuzz时序分析

#####mfuzz######
eset <- new("ExpressionSet",exprs = rs)

boxplot(rs)
eset <- filter.NA(eset, thres = 0.25)
eset <- fill.NA(eset, mode = 'mean')
eset <- filter.std(eset, min.std = 0)
eset <- standardise(eset)
c <- 8 #fuzzy c-means 聚类，需手动定义聚类个数，比方说设置 8 个簇
m <- mestimate(eset) #  评估出最佳的m值
set.seed(1234)
cl <- mfuzz(eset, c = c, m = m) # 聚类
#Color <- colorRampPalette(rev(c("#ff0000", "Yellow", "OliveDrab1")))(1000) ##修改颜色
#mfuzz.plot(eset,cl,mfrow = c(4,4),colo = Color)
mfuzz.plot(eset,cl,mfrow = c(4,4),time.labels = c(6,9,12,18,21,24,27))

得到时序分析图，分为八个cluster

 可以看出，cluster7和cluster8的变化较为显著，于是提取cluster7和cluster8中的基因进行后续分析

cl$size # 查看每个cluster中的基因个数
gene_7 <- cl$cluster[cl$cluster == 7] # 提取某个cluster下的基因
gene_8 <- cl$cluster[cl$cluster == 8] 
cl$membership # 查看基因和cluster之间的membership，用于判断基因所属簇，对应最大值的那个簇
output <- data.frame(DEGs_exp_averp)
output$cluster <- cl$cluster
#write.xlsx(as.data.frame(output),file = "./Mfuzz.xlsx",colNames = TRUE,rowNames=FALSE)