Hello,
I have a question about how to add a horizontal bar graph to a ggtree object. I would like to merge a ggtree object with a geom_boxplot object, where the abundances in the boxplot correspond to taxa on the ggtree tips. My understanding is that the facet_plot function would be ideal for this, but I have tried pretty hard and can't seem to get it to work. A reproducible example of what I can (and can't do) is attached using the Phyloseq package in the R environment.
I am sure there are other ways to do this, so this i just a rough skelton
There are different groups of abundance values for the facet_plot 'barploth' function, but they will not plot using this function, even though the data is there. Even if the tip names are unique, I am unable to get this to work. Happy to answer any other questions, and thanks in advance for the help! --Matt
#!usr/bin/env RScript
rm (list = ls())
library(phyloseq);library(ggstance);library(ggthemes);library(plyr);library(dplyr)
library(ggtree);library(tibble)
#load sample Phyloseq dataset
data("GlobalPatterns")
GP <- GlobalPatterns
GP <- prune_taxa(taxa_sums(GP) > 500, GP)
humantypes <- c("Feces", "Skin") #make it a bit simpler
sample_data(GP)$human <- get_variable(GP, "SampleType") %in% humantypes
mergedGP <- merge_samples(GP, "SampleType")
mergedGP <- rarefy_even_depth(mergedGP,rngseed=394582)
mergedGP <- tax_glom(mergedGP,"Order") #make the tree smaller for the purposes of this example
print(mergedGP) #final phyloseq object
#make a tree
tr <- phy_tree(mergedGP)
spec <- tax_table(mergedGP) %>% data.frame(stringsAsFactors = FALSE) %>% tibble::rownames_to_column("otu")
gt <- ggtree(tr, branch.length = "y") %<+% spec
gd <- gt$data
melt_data<-psmelt(mergedGP)
tt <- melt_data%>%left_join(select(gd,OTU=label,x,y),by="OTU") %>%
arrange(Sample) %>% mutate(sample2=factor(Sample,levels=unique(Sample)),
col=as.numeric(sample2),x.col=scales::rescale(col,to=c(1.3,2)))
#this will demonstrate that one can make a ggtree object with abundance data (other examples show this tpp)
g <- gt + geom_tippoint(data=gd$istip,aes(color=Phylum),size=2) +
geom_text(data=gd$istip,aes(label=Genus,x=x+0.01,size='smaller[1]'),hjust=0.1,check_overlap = T)+
geom_tile(data=tt,aes(x=x.col,y=y,fill=Sample,alpha=Abundance)) +
theme(legend.position="right") #+ geom_tiplab2()
g
#make a simple dataframe to build the boxplots
melt_simple <- subset(melt_data,select=c(OTU,Sample,Abundance))
#this just makes a simple tree
g.left <- gt + geom_tippoint(data=gd$istip,aes(color=Class),size=2) +
geom_text(data=gd$istip,aes(label=Genus,x=x+0.01,size='smaller[1]'),hjust=0.1,check_overlap = T)
g.left
#This function will produce the boxplot that I want, with boxes for each sample on the right side corresponding to that tree tip
ggplot(melt_simple,aes(x=OTU,y=Abundance))+
geom_boxplot(position=position_dodge(),colour="black",size=.3,alpha=0.5) +
scale_color_manual() +
scale_fill_manual() + coord_flip() +
scale_y_continuous(limits = c(0,300))
###code works until here###
#This is the step where I have an error
facet_plot(g.left, panel="Boxplot", data=melt_simple, geom_boxploth,
mapping = aes(x=Abundance,group=label,color="377eb8"),
outlier.shape=NA)