inst/doc/fitTimeSeries.R

In metagenomeSeq: Statistical analysis for sparse high-throughput sequencing

## ----include=FALSE------------------------------------------------------------
require(knitr)
opts_chunk$set(concordance=TRUE,tidy=TRUE)

## ----config,echo=FALSE-----------------------------------------
options(width = 65)
options(continue=" ")
options(warn=-1)
set.seed(42)

## ----requireMetagenomeSeq,warning=FALSE,message=FALSE----------
library(metagenomeSeq)
library(gss)

## ----dataset2,tidy=FALSE---------------------------------------
data(mouseData)
mouseData

## ----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]

## ----createMRexperiment2---------------------------------------
# This is an example of potential lvl's to aggregate by.
data(mouseData)
colnames(fData(mouseData))

## ----createMRexperiment3,tidy=FALSE----------------------------
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)

## ----plotMRexperiment3,tidy=FALSE------------------------------
par(mfrow=c(3,1))
plotClassTimeSeries(res1,pch=21,bg=classInfo)
plotTimeSeries(res2)
plotClassTimeSeries(res2,pch=21,bg=classInfo)

## ----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)

## ----timeSeriesAllClasses, tidy=FALSE--------------------------
set.seed(123)
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 = lapply(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)

## ----timeSeriesPlotting----------------------------------------
par(mfrow=c(2,1))
plotClassTimeSeries(res,pch=21,
                    bg=res$data$class,ylim=c(0,8))
plotTimeSeries(res)

## ----cite------------------------------------------------------
citation("metagenomeSeq")

## ----sessionInfo-----------------------------------------------
sessionInfo()