## ----setup, include=FALSE------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(cache = TRUE)
isdevel <- Biobase::package.version("Biobase") > 2.37

## ---- eval=FALSE---------------------------------------------------------
#  library(BiocInstaller)
#  biocLite(c("waldronlab/MicrobiomeWorkshop", "epiviz/metavizr@bioc-workshop"),
#           dependencies=TRUE)

## ----NeededPackagespac---------------------------------------------------
if(!("metavizr" %in% installed.packages()) || "1.1.3" != installed.packages()["metavizr","Version"])
  devtools::install_github("epiviz/metavizr@bioc-workshop", build_vignettes = TRUE)

## ------------------------------------------------------------------------
cranpkgs=c("ggplot2","devtools","vegan","httr")
BioCpkgs=c( "curatedMetagenomicData", "phyloseq")
ghpkgs= "metavizr"
allpkg <- c(cranpkgs, BioCpkgs, ghpkgs)
suppressPackageStartupMessages(sapply(allpkg, require, character.only = TRUE))

## ---- message=FALSE------------------------------------------------------
# BiocInstaller::biocLite("curatedMetagenomicData")  #Bioconductor version
# BiocInstaller::biocLite("waldronlab/curatedMetagenomicData")  #bleeding edge version
suppressPackageStartupMessages(library(curatedMetagenomicData))

## ----eval=FALSE----------------------------------------------------------
#  ?combined_metadata
#  View(combined_metadata)

## ------------------------------------------------------------------------
table(combined_metadata$antibiotics_current_use)
table(combined_metadata$disease)

## ---- message=FALSE------------------------------------------------------
suppressPackageStartupMessages(library(curatedMetagenomicData))
zeller.eset = ZellerG_2014.metaphlan_bugs_list.stool()

## ---- message=FALSE, results='hide', eval = isdevel----------------------
oral <- c("BritoIL_2016.metaphlan_bugs_list.oralcavity",
          "Castro-NallarE_2015.metaphlan_bugs_list.oralcavity")
esl <- curatedMetagenomicData(oral, dryrun = FALSE)
esl
esl[[1]]
esl[[2]]

## ---- eval=TRUE----------------------------------------------------------
curatedMetagenomicData("*metaphlan_bugs_list.stool*", dryrun = TRUE)

## ---- eval=isdevel-------------------------------------------------------
eset <- mergeData(esl)
eset

## ------------------------------------------------------------------------
experimentData( zeller.eset )

## ------------------------------------------------------------------------
head( pData( zeller.eset ) )

## ------------------------------------------------------------------------
exprs( zeller.eset )[1:6, 1:5]  #first 6 rows and 5 columns

## ------------------------------------------------------------------------
ExpressionSet2phyloseq(zeller.eset)
ExpressionSet2MRexperiment(zeller.eset)

## ---- eval=isdevel-------------------------------------------------------
ExpressionSet2phyloseq(eset)
ExpressionSet2MRexperiment(eset)

## ---- eval=isdevel-------------------------------------------------------
zeller.tree <- ExpressionSet2phyloseq( zeller.eset, phylogenetictree = TRUE)

## ---- eval=isdevel-------------------------------------------------------
wt = UniFrac(zeller.tree, weighted=TRUE, normalized=FALSE, 
             parallel=FALSE, fast=TRUE)
plot(hclust(wt), main="Weighted UniFrac distances")

## ------------------------------------------------------------------------
grep("coli", rownames(zeller.eset), value=TRUE)

## ------------------------------------------------------------------------
x = exprs( zeller.eset )[grep("s__Escherichia_coli$", rownames( zeller.eset)), ]
summary( x )

## ------------------------------------------------------------------------
hist( x, xlab = "Relative Abundance", main="Prevalence of E. Coli",
      breaks="FD")

## ------------------------------------------------------------------------
zeller.counts = sweep(exprs( zeller.eset ), 2, zeller.eset$number_reads / 100, "*")
zeller.counts = round(zeller.counts)
zeller.counts[1:6, 1:5]

## ------------------------------------------------------------------------
zeller.eset2 = curatedMetagenomicData("ZellerG_2014.metaphlan_bugs_list.stool",
                                     counts = TRUE, 
                                     dryrun = FALSE)[[1]]
all.equal(exprs(zeller.eset2), zeller.counts)

## ---- warning=FALSE------------------------------------------------------
suppressPackageStartupMessages(library(phyloseq))
zeller.pseq = ExpressionSet2phyloseq( zeller.eset )

## ---- eval=FALSE---------------------------------------------------------
#  zeller_MR_expr = ExpressionSet2MRexperiment(zeller.eset)
#  zeller_MR_expr <- zeller_MR_expr[, which(pData(zeller_MR_expr)$age < 60)]
#  zeller_MR_expr <- zeller_MR_expr[,-which(pData(zeller_MR_expr)$study_condition == "adenoma")]
#  
#  public_ipv4 <- try(httr::content(httr::GET("http://169.254.169.254/latest/meta-data/public-ipv4")), silent = TRUE)
#  
#  if(grepl("Timeout was reached", simpleMessage(public_ipv4))){
#    app <- startMetaviz(host="http://metaviz-dev.cbcb.umd.edu")
#  } else{
#    app <- startMetaviz(host="http://metaviz-dev.cbcb.umd.edu", ws_host = public_ipv4, daemonized = FALSE, try_ports = TRUE)
#  }

## ---- eval = FALSE-------------------------------------------------------
#  facetZoom <- app$plot(zeller_MR_expr, type = "InnerNodeCounts", datasource_name = "zeller", feature_order = colnames(fData(zeller_MR_expr)))
#  app$service()

## ---- eval = FALSE-------------------------------------------------------
#  heatmap_plot <- app$chart_mgr$revisualize(chart_type = "HeatmapPlot", chart = facetZoom)
#  app$service()

## ---- eval=FALSE---------------------------------------------------------
#  app$chart_mgr$rm_all_charts()

## ---- eval=FALSE---------------------------------------------------------
#  # Removes counts and features all levels of the taxonomy except for the lowest level
#  zeller_MR_expr <- zeller_MR_expr[-unname(which(is.na(fData(zeller_MR_expr)$Strain))),]
#  
#  # Remove the t__ from each leaf level feature in the MRexperiment counts data
#  rownames(zeller_MR_expr) <- sapply(strsplit(rownames(zeller_MR_expr), "__"), function(i){unname(i)[2]})

## ---- eval=FALSE---------------------------------------------------------
#  # Normalize counts using cummulative sum scaling method.
#  zeller_MR_expr <- filterData(zeller_MR_expr, present = 5)
#  zeller_MR_expr <- cumNorm(zeller_MR_expr, p = 0.75)
#  
#  # Set model to study_condition and use the fitFeatureModel to test
#  zeller_sample_data <-  pData(zeller_MR_expr)
#  mod <-  model.matrix(~1+study_condition, data = zeller_sample_data)
#  results_zeller <-  fitFeatureModel(zeller_MR_expr, mod)
#  logFC_zeller <- MRcoefs(results_zeller, number = nrow(results_zeller))

## ---- eval=FALSE---------------------------------------------------------
#  # Creating a separate MRexperiment to pass to fitZig
#  zig_zeller_MR_expr = ExpressionSet2MRexperiment(zeller.eset)
#  zig_zeller_MR_expr <- zig_zeller_MR_expr[, which(pData(zig_zeller_MR_expr)$age < 60)]
#  zig_zeller_MR_expr <- zig_zeller_MR_expr[,-which(pData(zig_zeller_MR_expr)$study_condition == "adenoma")]
#  zig_zeller_MR_expr <- zig_zeller_MR_expr[-unname(which(is.na(fData(zig_zeller_MR_expr)$Strain))),]
#  rownames(zig_zeller_MR_expr) <- sapply(strsplit(rownames(zig_zeller_MR_expr), "__"), function(i){unname(i)[2]})
#  zig_zeller_MR_expr <- filterData(zig_zeller_MR_expr, present = 5)
#  zig_zeller_MR_expr <- cumNorm(zig_zeller_MR_expr, p = 0.75)
#  
#  zig_zeller_sample_data <-  pData(zig_zeller_MR_expr)
#  
#  mod_zig <-  model.matrix(~1+study_condition+gender, data = zig_zeller_sample_data)
#  results_zeller_zig <-  fitZig(zig_zeller_MR_expr, mod_zig)
#  coefs_zeller_zig <- MRtable(results_zeller_zig, number = nrow(results_zeller_zig))

## ---- eval=FALSE---------------------------------------------------------
#  features <- rownames(logFC_zeller)
#  
#  featuresToKeep_names <- features[which(logFC_zeller[which(abs(logFC_zeller$logFC) > 2),]$adjPvalues <.1)]
#  featuresToKeep <- rep(2, length(featuresToKeep_names))
#  names(featuresToKeep) <-  featuresToKeep_names
#  
#  featuresToRemove_names <- features[!(features %in% featuresToKeep_names)]
#  featuresToRemove <- rep(0, length(featuresToRemove_names))
#  names(featuresToRemove) <- featuresToRemove_names
#  featureSelection = c(featuresToKeep, featuresToRemove)

## ---- eval=FALSE---------------------------------------------------------
#  facetZoom <- app$plot(zeller_MR_expr, type = "LeafCounts", datasource_name="zeller_leaf_level", feature_order = colnames(fData(zeller_MR_expr)))
#  app$service()
#  heatmap <- app$chart_mgr$revisualize(chart_type = "HeatmapPlot", chart = facetZoom)
#  app$service()
#  stackedPlot <- app$chart_mgr$revisualize(chart_type ="StackedLinePlot", chart = facetZoom)
#  app$service()

## ---- eval = FALSE-------------------------------------------------------
#  app$get_ms_object(facetZoom)$propagateHierarchyChanges(featureSelection, request_with_labels = TRUE)

## ---- eval=FALSE---------------------------------------------------------
#  #Aggregate to the Class level
#  aggregation_level = "Class"
#  agg_zeller_MR_expr <- aggregateByTaxonomy(zeller_MR_expr, lvl=aggregation_level)
#  
#  agg_zeller_MR_expr <-  cumNorm(agg_zeller_MR_expr, p = 0.75)
#  agg_mod <-  model.matrix(~1+study_condition, data = pData(agg_zeller_MR_expr))
#  agg_results_zeller <-  fitFeatureModel(agg_zeller_MR_expr, agg_mod)
#  agg_logFC_zeller <- MRcoefs(agg_results_zeller, number = nrow(agg_results_zeller))
#  
#  agg_features <- rownames(agg_logFC_zeller)
#  agg_featuresToKeep_names <- agg_features[which(agg_logFC_zeller[which(abs(agg_logFC_zeller$logFC) > 1),]$adjPvalues <.1)]
#  agg_featuresToKeep <- rep(2, length(agg_featuresToKeep_names))
#  names(agg_featuresToKeep) <- agg_featuresToKeep_names
#  
#  agg_featuresToRemove_names <- agg_features[!(agg_features %in% agg_featuresToKeep_names)]
#  agg_featuresToRemove <- rep(0, length(agg_featuresToRemove_names))
#  names(agg_featuresToRemove) <- agg_featuresToRemove_names
#  
#  agg_selectionUpdate <- c(agg_featuresToKeep, agg_featuresToRemove)
#  app$get_ms_object(facetZoom)$propagateHierarchyChanges(agg_selectionUpdate, request_with_labels = TRUE)