I made a plot that shows median expression of some markers across different cell populations. I made it using a Bioconductor package that is based on ggplot and the final result is a ggplot object. Here is the plot and I want to be able to show an analysis of variance between the conditions (so between "Ref" and "BCRXL") for every antigen within every cell type. I was thinking I could use ggpubr since my plot is ggplot based, but I'm having difficulty.

p <- plotMedExprs(sce, k = "merging1", 
    facet_by = "cluster_id", shape_by = "patient_id")
p$facet$params$ncol <- 2
p

I know this is a complex plot, but any help would be appreciated! :)