## ---- 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()