### R code from vignette source 'vignettes/metagenomeSeq/inst/doc/metagenomeSeq.Rnw'

###################################################
### code chunk number 1: metagenomeSeq.Rnw:16-18
###################################################
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)


###################################################
### code chunk number 4: dataset1
###################################################
data(lungData)
lungData


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


###################################################
### code chunk number 6: loadData
###################################################
dataDirectory <- system.file("extdata", package="metagenomeSeq")
lung = load_meta(file.path(dataDirectory,"CHK_NAME.otus.count.csv")) 
dim(lung$counts)


###################################################
### code chunk number 7: loadTaxa
###################################################
taxa = read.delim(file.path(dataDirectory,"CHK_otus.taxonomy.csv"),stringsAsFactors=F)[,2]
otu  = read.delim(file.path(dataDirectory,"CHK_otus.taxonomy.csv"),stringsAsFactors=F)[,1]


###################################################
### code chunk number 8: loadClin
###################################################
clin = load_phenoData(file.path(dataDirectory,"CHK_clinical.csv"),tran=TRUE)
ord = match(colnames(lung$counts),rownames(clin)) 
clin = clin[ord,]
head(clin[1:2,])


###################################################
### code chunk number 9: createMRexperiment1
###################################################
phenotypeData = as(clin,"AnnotatedDataFrame")
phenotypeData


###################################################
### code chunk number 10: createMRexperiment2
###################################################
OTUdata = as(lung$taxa,"AnnotatedDataFrame")
varLabels(OTUdata) = "taxa"
OTUdata


###################################################
### code chunk number 11: createMRexperiment3
###################################################
obj = newMRexperiment(lung$counts,phenoData=phenotypeData,featureData=OTUdata)
# Links to a paper providing further details can be included optionally.
# experimentData(obj) = annotate::pmid2MIAME("21680950")
obj


###################################################
### code chunk number 12: calculateNormFactors
###################################################
data(lungData)
p=cumNormStat(lungData)


###################################################
### code chunk number 13: normalizeData
###################################################
lungData = cumNorm(lungData,p=p)


###################################################
### code chunk number 14: saveData
###################################################
mat = MRcounts(lungData,norm=TRUE,log=TRUE)[1:5,1:5]
exportMat(mat,output=file.path(dataDirectory,"tmp.tsv"))


###################################################
### code chunk number 15: exportStats
###################################################
exportStats(lungData[,1:5],output=file.path(dataDirectory,"tmp.tsv"))
head(read.csv(file=file.path(dataDirectory,"tmp.tsv"),sep="\t"))


###################################################
### code chunk number 16: removeData
###################################################
system(paste("rm",file.path(dataDirectory,"temp.tsv")))


###################################################
### code chunk number 17: preprocess
###################################################
controls = grep("Extraction.Control",pData(lungData)$SampleType)
lungTrim = lungData[,-controls]
sparseFeatures = which(rowSums(MRcounts(lungTrim)>0)<10)
lungTrim = lungTrim[-sparseFeatures,]
lungp = cumNormStat(lungTrim,pFlag=TRUE,main="Trimmed lung data")
lungTrim = cumNorm(lungTrim,p=lungp)


###################################################
### code chunk number 18: zigTesting
###################################################
smokingStatus = pData(lungTrim)$SmokingStatus
bodySite = pData(lungTrim)$SampleType
mod = model.matrix(~smokingStatus+bodySite)
settings = zigControl(maxit=10,verbose=TRUE)
fit = fitZig(obj = lungTrim,mod=mod,control=settings)


###################################################
### code chunk number 19: fittedResult
###################################################
taxa = 
  sapply(strsplit(as.character(fData(lungTrim)$taxa),split=";"),
         function(i){i[length(i)]})
head(MRcoefs(fit,taxa=taxa,coef=2))


###################################################
### code chunk number 20: presenceAbsence
###################################################
data(mouseData)
classes = pData(mouseData)$diet
res = MRfisher(mouseData[1:5,],cl=classes)
# Warning - the p-value is calculating 1 despite a high odd's ratio.
head(res)


###################################################
### code chunk number 21: heatmapData
###################################################
data(mouseData)
trials = pData(mouseData)$diet
heatmapColColors=brewer.pal(12,"Set3")[as.integer(factor(trials))];
heatmapCols = colorRampPalette(brewer.pal(9, "RdBu"))(50)

# plotMRheatmap
plotMRheatmap(obj=mouseData,n=200,cexRow = 0.4,cexCol = 0.4,trace="none",
                col = heatmapCols,ColSideColors = heatmapColColors)

# plotCorr
plotCorr(obj=mouseData,n=200,cexRow = 0.25,cexCol = 0.25,
         trace="none",dendrogram="none",col=heatmapCols)


###################################################
### code chunk number 22: MDSandRareplots
###################################################
data(mouseData)
cl = factor(pData(mouseData)$diet)

# plotOrd
plotOrd(mouseData,tran=TRUE,usePCA=FALSE,useDist=TRUE,bg=cl,pch=21,xlab="1st coordinate",ylab="2nd coordinate")

# plotRare
res = plotRare(mouseData,cl=cl,ret=TRUE,pch=21,bg=cl)

# Linear fits for plotRare / legend
tmp=lapply(levels(cl), function(lv) 
  lm(res[,"ident"]~res[,"libSize"]-1, subset=cl==lv))
for(i in 1:length(levels(cl))){
   abline(tmp[[i]], col=i)
}
legend("topleft", c("Diet 1","Diet 2"), text.col=c(1,2),box.col=NA)


###################################################
### code chunk number 23: plotOTUData
###################################################
head(MRtable(fit,coef=2,taxa=1:length(fData(lungTrim)$taxa)))
patients=sapply(strsplit(rownames(pData(lungTrim)),split="_"),
          function(i){
            i[3]
          })
pData(lungTrim)$patients=patients
classIndex=list(smoker=which(pData(lungTrim)$SmokingStatus=="Smoker"))
classIndex$nonsmoker=which(pData(lungTrim)$SmokingStatus=="NonSmoker")
otu = 779

# plotOTU
plotOTU(lungTrim,otu=otu,classIndex,main="Neisseria meningitidis")

# Now multiple OTUs annotated similarly
x = fData(lungTrim)$taxa[otu]
otulist = grep(x,fData(lungTrim)$taxa)

# plotGenus
plotGenus(lungTrim,otulist,classIndex,labs=FALSE,
            main="Neisseria meningitidis")

lablist<- c("S","NS")
axis(1, at=seq(1,6,by=1), labels = rep(lablist,times=3))


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


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