### R code from vignette source 'metagenomeSeq.Rnw' ################################################### ### code chunk number 1: metagenomeSeq.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) ################################################### ### code chunk number 4: loadData ################################################### dataDirectory <- system.file("extdata", package="metagenomeSeq") lung = load_meta(file.path(dataDirectory,"CHK_NAME.otus.count.csv")) dim(lung$counts) ################################################### ### code chunk number 5: loadTaxa ################################################### taxa = read.delim(file.path(dataDirectory,"CHK_otus.taxonomy.csv"),stringsAsFactors=FALSE) ################################################### ### code chunk number 6: 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 7: createMRexperiment1 ################################################### phenotypeData = AnnotatedDataFrame(clin) phenotypeData ################################################### ### code chunk number 8: createMRexperiment2 ################################################### OTUdata = AnnotatedDataFrame(taxa) OTUdata ################################################### ### code chunk number 9: 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 10: dataset1 ################################################### data(lungData) lungData ################################################### ### code chunk number 11: dataset2 ################################################### data(mouseData) mouseData ################################################### ### code chunk number 12: pdata ################################################### phenoData(obj) head(pData(obj),3) ################################################### ### code chunk number 13: fdata ################################################### featureData(obj) head(fData(obj)[,-c(2,10)],3) ################################################### ### code chunk number 14: MRcounts ################################################### head(MRcounts(obj[,1:2])) ################################################### ### code chunk number 15: metagenomeSeq.Rnw:201-206 ################################################### featuresToKeep = which(rowSums(obj)>=100) samplesToKeep = which(pData(obj)$SmokingStatus=="Smoker") obj_smokers = obj[featuresToKeep,samplesToKeep] obj_smokers head(pData(obj_smokers),3) ################################################### ### code chunk number 16: filterData ################################################### data(mouseData) filterData(mouseData,present=10,depth=1000) ################################################### ### code chunk number 17: calculateNormFactors ################################################### data(lungData) p=cumNormStatFast(lungData) ################################################### ### code chunk number 18: normalizeData ################################################### lungData = cumNorm(lungData,p=p) ################################################### ### code chunk number 19: saveData ################################################### mat = MRcounts(lungData,norm=TRUE,log=TRUE)[1:5,1:5] exportMat(mat,file=file.path(dataDirectory,"tmp.tsv")) ################################################### ### code chunk number 20: exportStats ################################################### exportStats(lungData[,1:5],file=file.path(dataDirectory,"tmp.tsv")) head(read.csv(file=file.path(dataDirectory,"tmp.tsv"),sep="\t")) ################################################### ### code chunk number 21: removeData ################################################### system(paste("rm",file.path(dataDirectory,"tmp.tsv"))) ################################################### ### code chunk number 22: fitFeatureModel ################################################### data(lungData) lungData = lungData[,-which(is.na(pData(lungData)$SmokingStatus))] lungData=filterData(lungData,present=30,depth=1) lungData <- cumNorm(lungData, p=.5) s <- normFactors(lungData) pd <- pData(lungData) mod <- model.matrix(~1+SmokingStatus, data=pd) lungres1 = fitFeatureModel(lungData,mod) head(MRcoefs(lungres1)) ################################################### ### code chunk number 23: preprocess ################################################### data(lungData) controls = grep("Extraction.Control",pData(lungData)$SampleType) lungTrim = lungData[,-controls] rareFeatures = which(rowSums(MRcounts(lungTrim)>0)<10) lungTrim = lungTrim[-rareFeatures,] lungp = cumNormStat(lungTrim,pFlag=TRUE,main="Trimmed lung data") lungTrim = cumNorm(lungTrim,p=lungp) ################################################### ### code chunk number 24: zigTesting ################################################### smokingStatus = pData(lungTrim)$SmokingStatus bodySite = pData(lungTrim)$SampleType normFactor = normFactors(lungTrim) normFactor = log2(normFactor/median(normFactor) + 1) mod = model.matrix(~smokingStatus+bodySite + normFactor) settings = zigControl(maxit=10,verbose=TRUE) fit = fitZig(obj = lungTrim,mod=mod,useCSSoffset = FALSE, control=settings) # The default, useCSSoffset = TRUE, automatically includes the CSS scaling normalization factor. ################################################### ### code chunk number 25: contrasts ################################################### # maxit=1 is for demonstration purposes settings = zigControl(maxit=1,verbose=FALSE) mod = model.matrix(~bodySite) colnames(mod) = levels(bodySite) # fitting the ZIG model res = fitZig(obj = lungTrim,mod=mod,control=settings) # The output of fitZig contains a list of various useful items. hint: names(res). # # Probably the most useful is the limma 'MLArrayLM' object called fit. zigFit = res$fit finalMod = res$fit$design contrast.matrix = makeContrasts(BAL.A-BAL.B,OW-PSB,levels=finalMod) fit2 = contrasts.fit(zigFit, contrast.matrix) fit2 = eBayes(fit2) topTable(fit2) # See help pages on decideTests, topTable, topTableF, vennDiagram, etc. ################################################### ### code chunk number 26: 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 27: timeSeries ################################################### # vignette("fitTimeSeries") ################################################### ### code chunk number 28: perm ################################################### coeffOfInterest = 2 res = fitLogNormal(obj = lungTrim, mod = mod, useCSSoffset = FALSE, B = 10, coef = coeffOfInterest) # extract p.values and adjust for multiple testing # res$p are the p-values calculated through permutation adjustedPvalues = p.adjust(res$p,method="fdr") # extract the absolute fold-change estimates foldChange = abs(res$fit$coef[,coeffOfInterest]) # determine features still significant and order by the sigList = which(adjustedPvalues <= .05) sigList = sigList[order(foldChange[sigList])] # view the top taxa associated with the coefficient of interest. head(taxa[sigList]) ################################################### ### code chunk number 29: presenceAbsence ################################################### classes = pData(mouseData)$diet res = fitPA(mouseData[1:5,],cl=classes) # Warning - the p-value is calculating 1 despite a high odd's ratio. head(res) ################################################### ### code chunk number 30: discOdds ################################################### classes = pData(mouseData)$diet res = fitDO(mouseData[1:100,],cl=classes,norm=FALSE,log=FALSE) head(res) ################################################### ### code chunk number 31: corTest ################################################### cors = correlationTest(mouseData[55:60,],norm=FALSE,log=FALSE) head(cors) ################################################### ### code chunk number 32: uniqueFeatures ################################################### cl = pData(mouseData)[["diet"]] uniqueFeatures(mouseData,cl,nsamples = 10,nreads = 100) ################################################### ### code chunk number 33: aggTax ################################################### obj = aggTax(mouseData,lvl='phylum',out='matrix') head(obj[1:5,1:5]) ################################################### ### code chunk number 34: aggSamp ################################################### obj = aggSamp(mouseData,fct='mouseID',out='matrix') head(obj[1:5,1:5]) ################################################### ### code chunk number 35: interactiveDisplay ################################################### # Calling display on the MRexperiment object will start a browser session with interactive plots. # require(interactiveDisplay) # display(mouseData) ################################################### ### code chunk number 36: heatmapData ################################################### 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 37: MDSandRareplots ################################################### cl = factor(pData(mouseData)$diet) # plotOrd - can load vegan and set distfun = vegdist and use dist.method="bray" plotOrd(mouseData,tran=TRUE,usePCA=FALSE,useDist=TRUE,bg=cl,pch=21) # plotRare res = plotRare(mouseData,cl=cl,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 38: 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 39: plotFeatureData ################################################### classIndex=list(Western=which(pData(mouseData)$diet=="Western")) classIndex$BK=which(pData(mouseData)$diet=="BK") otuIndex = 8770 # par(mfrow=c(1,2)) dates = pData(mouseData)$date plotFeature(mouseData,norm=FALSE,log=FALSE,otuIndex,classIndex, col=dates,sortby=dates,ylab="Raw reads") ################################################### ### code chunk number 40: cite ################################################### citation("metagenomeSeq") ################################################### ### code chunk number 41: sessionInfo ################################################### sessionInfo()