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 }