基本打通RAIN+Mfuzz的流程

#导入数据-------------------------------------------------------------------------------------------
rm(list=ls())
setwd("")#getwd()
data <-read.csv("GSE201207_count_216.csv.gz")
kidney <- data[,c(1,26:37)]
#整理数据-----------------------------------------------------------------------------------------------------
library("biomaRt")
mart <- useMart("ENSEMBL_MART_ENSEMBL", dataset="mmusculus_gene_ensembl")
gene_symbols <- getBM(attributes=c("external_gene_name"),
                      filters="external_gene_name", 
                      values=kidney$X, 
                      mart=mart)
kidney <- merge(gene_symbols,kidney,by.x=1, by.y=1)
library(dplyr)
kidney <- na.omit(distinct(kidney, external_gene_name, .keep_all = T))
rownames(kidney) <- kidney$external_gene_name
kidney <- kidney[,-1]
write.table(kidney, file = "output.csv", sep = ",")
#绘制单一基因表达量的时间变化曲线----------------------------------------------------------
require('lattice')
#Pcnt <- kidney["Pcnt",(0:5 * 2 + rep(c(1, 2), each = 6))]
Pcnt <- kidney["Pcnt",]
xyplot(Pcnt~rep(0:11 * 4 + 2, each = 1), type = 'b', pch = 16, 
       xlab = 'time', ylab = 'expression value')
#rhythm analysis: RAIN--------------------------------------------------------------------- 
kidney.rain <- kidney[,c(1,7,2,8,3,9,4,10,5,11,6,12)]
library(rain)
results <- rain(t(kidney.rain),period = 24, deltat = 4, nr.series = 2,
                peak.border = c(0.3,0.7), verbose = FALSE)
best <- order(results$pVal)[1:10]
xyplot(as.matrix(kidney
                 [best, ])
       ~rep(0:11 * 4 + 2, each = 10) |rownames(kidney)[best],
       scales = list(y = list(relation = 'free')),
       layout = c(2, 5), type = 'b', pch = 16, xlab = 'time',
       ylab = 'expression value', cex.lab = 1)
select <- order(results$pVal)[1:50]
kidney.Mfuzz <- kidney[select, ]
#trend analysis: Mfuzz---------------------------------------------------------------------
require(Biobase)
bio.kidney.Mfuzz <- ExpressionSet(as.matrix(kidney.Mfuzz))
require("Mfuzz")
bio.kidney.Mfuzz <- standardise(
  filter.std(
     fill.NA(
      filter.NA(bio.kidney.Mfuzz, thres = 0.25), 
      mode = "mean"), 
      min.std = 0)
  )
m1 <- mestimate(bio.kidney.Mfuzz)
c1 <- mfuzz(bio.kidney.Mfuzz, c=3, m=m1)
mfuzz.plot(bio.kidney.Mfuzz, cl=c1, 
           mfrow = c(2,2),time.labels = seq(0,44,4))