if (interactive()){
set.seed(1)
data(yeast) # loading the reduced CDC28 yeast set (from the Mfuzz package)
# Data preprocessing
yeast <- filter.NA(yeast) # filters genes with more than 25% of the expression values missing
yeast <- fill.NA(yeast) # for illustration only; rather use knn method for
yeast <- standardise(yeast)
#
T.yeast <- 85 # cell cycle period (t=85min)
times.yeast <- pData(yeast)$time # time of measurements
#
yeast.test <- yeast[1:600,] # To speed up the example
#
NN <- 50 # number of generated background models
# Here, a small number was chosen for demonstration purpose.
# For the actual analysis, rather set N = 1000
# Calculation of FDRs
# i) based on random permutation as background model
fdr.rr <- fdrfourier(eset=yeast.test,T=T.yeast,
times=times.yeast,background.model="rr",N=NN,progress=TRUE)
# ii) based on Gaussian distribution
fdr.g <- fdrfourier(eset=yeast.test,T=T.yeast,
times=times.yeast,background.model="gauss",N=NN,progress=TRUE)
# iii) based on AR(1) models as background
fdr.ar1 <- fdrfourier(eset=yeast.test,T=T.yeast,
times=times.yeast,background.model="ar1",N=NN,progress=TRUE)
# Number of significant genes based on diff. background models
sum(fdr.rr$fdr < 0.1)
sum(fdr.g$fdr < 0.1)
sum(fdr.ar1$fdr < 0.1)
# Plot top scoring gene
plot(times.yeast,exprs(yeast.test)[order(fdr.ar1$fdr)[1],],type="o",
xlab="Time",ylab="Expression",
main=paste(featureNames(yeast.test)[order(fdr.ar1$fdr)[1]],"-- FDR:",
fdr.ar1$fdr[order(fdr.ar1$fdr)[1]]))
# List significant genes
fdr.ar1$fdr[which(fdr.ar1$fdr < 0.1)]
}
Run the code above in your browser using DataLab