maaslin2_plot

maaslin2_plot

make_results_plot <- function(path, DEBUG = FALSE, reorder=FALSE, qval_thresh=.1){
  alpha_levels <- c(".0001 or less", ".001", ".01", ".05", ".05 or above")
  alpha_breaks <- rev(c(.09, .5, .62, .7, .85))
  nonscaled <- read.csv(path, sep="\t") 
  if(nrow(nonscaled) == 0){
    warning("Error running maaslin, results file empty")
    return(NULL)
  }
  any_sig <- nonscaled %>% group_by(feature) %>% filter(any(qval < qval_thresh)) %>% pull(feature) %>% unique()  #changed any qval <.5 from .05
  if (length(any_sig) == 0){
    warning("No significant results to plot")
    return (NULL)
  } else if (length(any_sig) == 1){
    warning("not reording; single significant result")
    reorder <- FALSE
  }
  if(reorder){
    nonscaled$feature <- factor(nonscaled$feature, levels = labels(plot_dat_clustering))
    plot_dat_clustering <- nonscaled %>% filter(feature %in% any_sig) %>% 
      select(feature, metadata, value, coef ) %>%
      pivot_wider(names_from=c(metadata, value), values_from=coef) %>% 
      column_to_rownames("feature") %>%
      as.matrix() %>%
      dist(.) %>% 
      hclust() %>% 
      as.dendrogram()
  }
  pdata <- nonscaled %>%  
    filter(feature %in% any_sig) %>% 
    mutate(
      fac = metadata, 
      qval_bin = as.factor(
        case_when(
          qval <.0001 ~ alpha_levels[1],
          qval <.001 ~ alpha_levels[2],
          qval <.01 ~ alpha_levels[3],
          qval <.05 ~ alpha_levels[4],
          qval >=.05 ~ alpha_levels[5]
        )),
      x=value)
  
  thisxlab="Parameter"
  p <- ggplot(pdata,
              aes(x=x, 
                  y=feature, 
                  alpha = qval_bin, 
                  color = ifelse(coef >0, "positive", "negative"), 
                  #fill = sign(coef), 
                  size=abs(coef) )) +
    geom_point(stroke=0) + #, guide = guide_legend(override.aes = list(size = 3) )) +
    scale_color_manual(values=c("blue", "red")) + 
    theme(plot.caption = element_text(hjust=0) ,
          axis.text.x = element_text(angle=45, hjust=1),
          plot.background = element_rect(fill="white"),
          panel.background = element_rect(fill="white"),
          plot.title = element_text(size=4),
          plot.subtitle = element_text(size=6)
    ) + 
    scale_alpha_manual(values=alpha_breaks, breaks = alpha_levels) + 
    facet_grid(~fac, scales="free_x",space = "free_x") + 
    labs(x=thisxlab, y="Species", color="Relative Change", size="Magnitude", alpha="P-value\nBH-corrected",  
         #subtitle=thissubtitle,
         title = dirname(path))
  if (DEBUG) p <- p + scale_size(limits = c(0, 1)) 
  p
}