inst/doc/philr-intro.R

In philr: Phylogenetic partitioning based ILR transform for metagenomics data

## ---- echo=TRUE, message=FALSE------------------------------------------------
library(philr); packageVersion("philr")
library(phyloseq); packageVersion("phyloseq")
library(ape); packageVersion("ape")
library(ggplot2); packageVersion("ggplot2")

data(GlobalPatterns)

## ---- message=FALSE, warning=FALSE--------------------------------------------
GP <-  filter_taxa(GlobalPatterns, function(x) sum(x > 3) > (0.2*length(x)), TRUE)
GP <-  filter_taxa(GP, function(x) sd(x)/mean(x) > 3.0, TRUE)
GP <- transform_sample_counts(GP, function(x) x+1)

## ---- echo=FALSE--------------------------------------------------------------
GP

## ---- message=FALSE, warning=FALSE--------------------------------------------
is.rooted(phy_tree(GP)) # Is the tree Rooted?
is.binary.tree(phy_tree(GP)) # All multichotomies resolved?

## ---- message=FALSE, warning=FALSE--------------------------------------------
phy_tree(GP) <- makeNodeLabel(phy_tree(GP), method="number", prefix='n')

## -----------------------------------------------------------------------------
name.balance(phy_tree(GP), tax_table(GP), 'n1')

## -----------------------------------------------------------------------------
otu.table <- t(otu_table(GP))
tree <- phy_tree(GP)
metadata <- sample_data(GP)
tax <- tax_table(GP)

otu.table[1:2,1:2] # OTU Table
tree # Phylogenetic Tree
head(metadata,2) # Metadata
head(tax,2) # taxonomy table

## -----------------------------------------------------------------------------
gp.philr <- philr(otu.table, tree, 
                  part.weights='enorm.x.gm.counts', 
                  ilr.weights='blw.sqrt')
gp.philr[1:5,1:5]

## -----------------------------------------------------------------------------
gp.dist <- dist(gp.philr, method="euclidean")
gp.pcoa <- ordinate(GP, 'PCoA', distance=gp.dist)
plot_ordination(GP, gp.pcoa, color='SampleType') + geom_point(size=4)

## ---- message=FALSE, warning=FALSE--------------------------------------------
sample_data(GP)$human <- factor(get_variable(GP, "SampleType") %in% c("Feces", "Mock", "Skin", "Tongue"))

## ---- message=FALSE, warning=FALSE--------------------------------------------
library(glmnet); packageVersion('glmnet')
glmmod <- glmnet(gp.philr, sample_data(GP)$human, alpha=1, family="binomial")

## -----------------------------------------------------------------------------
top.coords <- as.matrix(coefficients(glmmod, s=0.2526))
top.coords <- rownames(top.coords)[which(top.coords != 0)]
(top.coords <- top.coords[2:length(top.coords)]) # remove the intercept as a coordinate

## -----------------------------------------------------------------------------
tc.names <- sapply(top.coords, function(x) name.balance(tree, tax, x))
tc.names

## -----------------------------------------------------------------------------
votes <- name.balance(tree, tax, 'n730', return.votes = c('up', 'down'))

votes[[c('up.votes', 'Family')]] # Numerator at Family Level
votes[[c('down.votes', 'Family')]] # Denominator at Family Level

## ---- message=FALSE, warning=FALSE--------------------------------------------
library(ggtree); packageVersion("ggtree")
library(dplyr); packageVersion('dplyr')

## ---- message=FALSE, warning=FALSE--------------------------------------------
tc.nn <- name.to.nn(tree, top.coords)
tc.colors <- c('#a6cee3', '#1f78b4', '#b2df8a', '#33a02c', '#fb9a99')
p <- ggtree(tree, layout='fan') +
  geom_balance(node=tc.nn[1], fill=tc.colors[1], alpha=0.6) +
  geom_balance(node=tc.nn[2], fill=tc.colors[2], alpha=0.6) +
  geom_balance(node=tc.nn[3], fill=tc.colors[3], alpha=0.6) +
  geom_balance(node=tc.nn[4], fill=tc.colors[4], alpha=0.6) +
  geom_balance(node=tc.nn[5], fill=tc.colors[5], alpha=0.6)
p <- annotate_balance(tree, 'n16', p=p, labels = c('n16+', 'n16-'),
                 offset.text=0.15, bar=FALSE)
annotate_balance(tree, 'n730', p=p, labels = c('n730+', 'n730-'),
                 offset.text=0.15, bar=FALSE)

## -----------------------------------------------------------------------------
gp.philr.long <- convert_to_long(gp.philr, get_variable(GP, 'human')) %>%
  filter(coord %in% top.coords)

ggplot(gp.philr.long, aes(x=labels, y=value)) +
  geom_boxplot(fill='lightgrey') +
  facet_grid(.~coord, scales='free_x') +
  xlab('Human') + ylab('Balance Value') +
  theme_bw()

## ---- message=FALSE, warning=FALSE--------------------------------------------
library(tidyr); packageVersion('tidyr')

gp.philr.long %>%
  rename(Human=labels) %>%
  filter(coord %in% c('n16', 'n730')) %>%
  spread(coord, value) %>%
  ggplot(aes(x=n16, y=n730, color=Human)) +
  geom_point(size=4) +
  xlab(tc.names['n16']) + ylab(tc.names['n730']) +
  theme_bw()