Phyloseq Barplot, customizable ggplot2

ggphylobarplot.R

ggphylobarplot <- function(physobj, 
      barlevel = "lowest", 
      colorlevel = "phylum", 
      splitvar = "Time", 
      splitlevelorder = c("1w", "1m", "3m"), 
      splitlevellabels = c("1 week", "1 month", "3 months"), 
      sortlevel = "1w", 
      mracut = .009, 
      agglomerate = FALSE, 
      agglomeratelevel = "genus", 
      decr = FALSE, 
      ycut = NULL, 
      colors = c("orange3","red4","chartreuse4","darkorchid2","dodgerblue3"),
      justDF = FALSE,
      breakinterval = 0.1) {
  physobj <- physobj %>% transform_sample_counts(function(x) x/sum(x))

  if(agglomerate)
    physobj <- tax_glom(physobj, barlevel)

  if(is.null(splitlevelorder))
      splitlevelorder <- get_variable(physobj, splitvar) %>% as.character %>% unique

  subphys <- lapply(splitlevelorder, function(level) {
    index <- get_variable(physobj, splitvar) %>% as.character %in% level
    if(sum(index) == 0)
      stop("Invalid split level: ", level)
    prune_samples(index, physobj)
  }) %>% setNames(splitlevelorder)

  means <- lapply(subphys, function(x) taxa_sums(x)/nsamples(x)) %>% do.call("cbind", .) %>% data.frame(bug = rownames(.), ., stringsAsFactors = FALSE) %>% setNames(c("bug", splitlevelorder))
  means$keep <- rowSums(select(means, -bug) > mracut) > 0

  others <- means %>% filter(!keep) %>% select(-bug, -keep) %>% colSums %>% data.frame %>% t
  means <- means %>% filter(keep) %>% select(-keep) %>% mutate(name = tax_table(physobj)[.[,"bug"],barlevel] %>% as.character) %>% rbind(cbind(bug = "Other", others, name = "Other"))
  means <- means[order(means[,sortlevel], decreasing = decr),] %>% arrange(xor(decr,name != "Other")) %>% mutate(name = factor(name, levels = name))
  means[means$bug != "Other","colvar"] <- tax_table(physobj)[means %>% filter(bug != "Other") %>% .[,"bug"],colorlevel] %>% as.character

  real_levels <- sort(unique(means$colvar %>% na.omit))
  if(length(real_levels) > length(colors) & !is.null(colors))
    stop(paste0("Not enough colors. Need ", length(real_levels), " colors!"))

  plotdat <- gather(means, group, value, -bug, -name, -colvar) %>% 
    mutate( value = as.numeric(value), 
            group = factor(group, levels = splitlevelorder, labels = splitlevellabels),
            percent = paste0(sprintf("%.2f",value*100), "%"))

  if(is.null(ycut))
    ycut <- max(plotdat$value, na.rm = T)
  if(is.null(colors))
    colors <- rainbow(length(real_levels))
  colors <- c(colors, "grey")
  plotdat <- plotdat %>% mutate(colvar = factor(ifelse(is.na(colvar), "Other", colvar), levels = c(real_levels, "Other")))

  if(justDF)
    return(plotdat)

  ggplot(plotdat, aes(y=as.numeric(name), x=value)) + facet_grid(.~group, scales = "free", space = "free") + 
    geom_rect(aes(fill = colvar, xmin = 0, xmax = value, ymin = as.numeric(name) + 0.4, ymax = as.numeric(name) - 0.4)) + 
    geom_text(aes(label = percent), nudge_x = 0.1, hjust = 1) +
    ylab(NULL) + 
    scale_x_continuous(name = "Mean relative abundance", labels = scales::percent, breaks = seq(0,1,breakinterval)) + 
    scale_y_continuous(breaks = 1:length(unique(plotdat$name)), labels = levels(plotdat$name), limits = c(1-.41,length(unique(plotdat$name))+.41)) +
    scale_fill_manual(values = colors, name = NULL, breaks = real_levels) +
    theme_bw() + theme(strip.background = element_blank(), panel.grid.minor.y = element_blank())
}