## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----eval = FALSE-------------------------------------------------------------
# # if (!require("remotes", quietly = TRUE))
# #     install.packages("remotes")
# #
# # remotes::install_github("statdivlab/radEmu")

## ----setup, message = FALSE---------------------------------------------------
library(magrittr)
library(dplyr)
library(ggplot2)
library(stringr)
library(radEmu)

## ----results = 'hide'---------------------------------------------------------
phy <- requireNamespace("phyloseq", quietly = TRUE) == TRUE

## ----echo = FALSE-------------------------------------------------------------
print(paste0("phyloseq is installed: ", phy))

## ----eval = phy---------------------------------------------------------------
data(wirbel_sample)
data(wirbel_otu)
data(wirbel_taxonomy)
wirbel_phylo <- phyloseq::phyloseq(phyloseq::sample_data(wirbel_sample),
                                   phyloseq::otu_table(wirbel_otu, taxa_are_rows = FALSE),
                                   phyloseq::tax_table(wirbel_taxonomy))
wirbel_phylo

## ----eval = phy---------------------------------------------------------------
dim(phyloseq::sample_data(wirbel_phylo))
head(phyloseq::sample_data(wirbel_phylo))

## ----eval = phy---------------------------------------------------------------
dim(phyloseq::otu_table(wirbel_phylo))
# let's check out a subset
phyloseq::otu_table(wirbel_phylo)[1:5, 1:3]

## ----eval = phy---------------------------------------------------------------
mOTU_names <- colnames(phyloseq::otu_table(wirbel_phylo))

## ----eval = phy---------------------------------------------------------------
head(phyloseq::tax_table(wirbel_phylo))

## ----eval = phy---------------------------------------------------------------
phyloseq::sample_data(wirbel_phylo)$Group <- factor(phyloseq::sample_data(wirbel_phylo)$Group, levels = c("CTR","CRC"))

## ----eval = phy---------------------------------------------------------------
chosen_genera <- c("Eubacterium", "Faecalibacterium", "Fusobacterium", "Porphyromonas")
wirbel_restrict <- phyloseq::subset_taxa(wirbel_phylo, genus %in% chosen_genera)

## ----eval = phy---------------------------------------------------------------
wirbel_china <- phyloseq::subset_samples(wirbel_restrict, Country == "CHI")

## -----------------------------------------------------------------------------
sum(rowSums(phyloseq::otu_table(wirbel_china)) == 0) # no samples have a count sum of 0 
sum(colSums(phyloseq::otu_table(wirbel_china)) == 0) # one category has a count sum of 0 
category_to_rm <- names(which(colSums(phyloseq::otu_table(wirbel_china)) == 0))
wirbel_china <- phyloseq::subset_taxa(wirbel_china, species != category_to_rm)
sum(colSums(phyloseq::otu_table(wirbel_china)) == 0) # now no categories have a count sum of 0 

## ----eval = phy---------------------------------------------------------------
ch_fit <- emuFit(formula = ~ Group, 
                 Y = wirbel_china, 
                 run_score_tests = FALSE) 

## ----fig.height = 6, fig.width = 6, eval = phy--------------------------------
plot(ch_fit)$plots