### R code from vignette source 'fitTimeSeries.Rnw'

###################################################
### code chunk number 1: fitTimeSeries.Rnw:15-17
###################################################
require(knitr)
opts_chunk$set(concordance=TRUE,tidy=TRUE)


###################################################
### code chunk number 2: config
###################################################
options(width = 65)
options(continue=" ")
options(warn=-1)
set.seed(42)


###################################################
### code chunk number 3: requireMetagenomeSeq
###################################################
library(metagenomeSeq)
library(gss)


###################################################
### code chunk number 4: dataset2
###################################################
data(mouseData)
mouseData


###################################################
### code chunk number 5: createMRexperiment1
###################################################
# Creating mock sample replicates
sampleID = rep(paste("sample",1:10,sep=":"),times=20)
# Creating mock class membership
class = rep(c(rep(0,5),rep(1,5)),times=20)
# Creating mock time
time = rep(1:20,each=10)

phenotypeData = AnnotatedDataFrame(data.frame(sampleID,class,time))
# Creating mock abundances
set.seed(1)
# No difference
measurement1 = rnorm(200,mean=100,sd=1)
# Some difference
measurement2 = rnorm(200,mean=100,sd=1)
measurement2[1:5]=measurement2[1:5] + 100
measurement2[11:15]=measurement2[11:15] + 100
measurement2[21:25]=measurement2[21:25] + 50
mat = rbind(measurement1,measurement2)
colnames(mat) = 1:200
mat[1:2,1:10]


###################################################
### code chunk number 6: createMRexperiment2
###################################################
# This is an example of potential lvl's to aggregate by.
data(mouseData)
colnames(fData(mouseData))


###################################################
### code chunk number 7: createMRexperiment3
###################################################
obj = newMRexperiment(counts=mat,phenoData=phenotypeData)
obj
res1 = fitTimeSeries(obj,feature=1,
              class='class',time='time',id='sampleID',
              B=10,norm=FALSE,log=FALSE)
res2 = fitTimeSeries(obj,feature=2,
              class='class',time='time',id='sampleID',
              B=10,norm=FALSE,log=FALSE)


classInfo = factor(res1$data$class)
par(mfrow=c(3,1))
plotClassTimeSeries(res1,pch=21,bg=classInfo)
plotTimeSeries(res2)
plotClassTimeSeries(res2,pch=21,bg=classInfo)


###################################################
### code chunk number 8: timeSeries
###################################################
res = fitTimeSeries(obj=mouseData,lvl="class",feature="Actinobacteria",class="status",id="mouseID",time="relativeTime",B=10)

# We observe a time period of differential abundance for "Actinobacteria"
res$timeIntervals

str(res)


###################################################
### code chunk number 9: timeSeriesAllClasses
###################################################
classes = unique(fData(mouseData)[,"class"])

timeSeriesFits = lapply(classes,function(i){
        fitTimeSeries(obj=mouseData,
            feature=i,
            class="status",
            id="mouseID",
            time="relativeTime",
            lvl='class',
            C=.3,# a cutoff for 'interesting' 
            B=1) # B is the number of permutations and should clearly not be 1
    })
names(timeSeriesFits) = classes

# Removing classes of bacteria without a potentially
# interesting time interval difference.
timeSeriesFits = sapply(timeSeriesFits,function(i){i[[1]]})[-grep("No",timeSeriesFits)]

# Naming the various interesting time intervals.
for(i in 1:length(timeSeriesFits)){
    rownames(timeSeriesFits[[i]]) = 
      paste(
        paste(names(timeSeriesFits)[i]," interval",sep=""),
        1:nrow(timeSeriesFits[[i]]),sep=":"
      )
}

# Merging into a table.
timeSeriesFits = do.call(rbind,timeSeriesFits)

# Correcting for multiple testing.
pvalues = timeSeriesFits[,"p.value"]
adjPvalues = p.adjust(pvalues,"bonferroni")
timeSeriesFits = cbind(timeSeriesFits,adjPvalues)

head(timeSeriesFits)


###################################################
### code chunk number 10: timeSeriesPlotting
###################################################
par(mfrow=c(2,1))
plotClassTimeSeries(res,pch=21,
                    bg=res$data$class,ylim=c(0,8))
plotTimeSeries(res)


###################################################
### code chunk number 11: cite
###################################################
citation("metagenomeSeq")


###################################################
### code chunk number 12: sessionInfo
###################################################
sessionInfo()