I am writing on this forum because the phyloseq github issues page is not being answered and I need help with a very specific formating question. Ultimately, I want to have a grid of graphs like this where the column headers (where asia, europe, and north america are) is body site, sex, and pregnancy status, while my 3 row headers are the alpha diversity metrics (shannon, observed, and Inverse Simpson). I want to use the plot_richness function with phyloseq and create diversity metrics for each of my categories.

I achieved getting three plots for body site with the code:

alpha_body_site <- plot_richness(ps_rarefied, x="body.site", measures = c("Shannon", "Observed", "InvSimpson"), color = "body.site")

#remove phyloseq dots
alpha_body_site$layers <- alpha_body_site$layers[-1]

alpha_body_site <- alpha_body_site +
  geom_boxplot(aes(fill = body.site), outlier.size = -1, color = "black") +
  scale_fill_manual(values = ghibli_palette("PonyoMedium")) +
  theme(strip.background = element_rect(color = "black", fill = "grey")) +
  theme(strip.text.x = element_text(color = "black", face = "bold", size = 12)) +
  theme(panel.background = element_rect(fill = "NA"), panel.grid.major = element_line(color = "grey90")) +
  theme(panel.border = element_rect(color = "black", fill = NA), legend.position = "none") +
  theme(axis.title.y.left = element_text(size = 14), axis.text.x.bottom = element_text(size = 12), axis.title.x.bottom = element_text(size = 14), axis.text.y.left = element_text(size = 12))

which creates a row of 3 alpha diversity plots by body site. I now want to add significance bars above the box plots that are significantly different, so I thought I need to expand the y limits to create more space at the top of the plots to fit the bars. When I try to expand the y limits with ylim() it will change all 3 plots to be the same metric which throws off the box plot sizes because shannon for example only goes up to 6 while observed is from 0-1000. I believe I will need to make a separate box plot for every diversity metric to achieve the look im going for, but I want to be able to format the final figure to look like the grid of graphs I mentioned above. I was able to create the plots I want by doing the following code for each diversity metric:

shannon_alpha_body_site <- plot_richness(ps_rarefied, x="body.site", measures = c("Shannon"), color = "body.site") + ylim(1,6.6)

shannon_alpha_body_site$layers <- shannon_alpha_body_site$layers[-1]

shannon_alpha_body_site <- shannon_alpha_body_site +
  geom_boxplot(aes(fill = body.site), outlier.size = -1, color = "black") +
  scale_fill_manual(values = ghibli_palette("PonyoMedium")) +
  theme(strip.background = element_rect(color = "black", fill = "grey")) +
  theme(strip.text.x = element_text(color = "black", face = "bold", size = 12)) +
  theme(panel.background = element_rect(fill = "NA"), panel.grid.major = element_line(color = "grey90")) +
  theme(panel.border = element_rect(color = "black", fill = NA), legend.position = "none") +
  theme(axis.title.y.left = element_text(size = 14), axis.text.x.bottom = element_text(size = 12, face = "bold")) +
  theme(axis.title.x.bottom = element_text(size = 14), axis.text.y.left = element_text(size = 10, face = "bold")) +
  theme(aspect.ratio = 2) + stat_pvalue_manual(shannon_df_filtered, y.position = c(5.5,6,6.5), hide.ns = TRUE, label = "p.signif", position = position_dodge(0.8)) +
  labs(y = NULL, x = NULL)

I see people use facet_grid() to create the grid look, but I can’t figure out how to make it work with separate plots, or re-format my data in a way that would make this grid work. Is there a way I can achieve this grid look, while also having the significance bars for each plot? Thanks for your help.