currently this code 


colvec <- c("White", "#0073C2","#A730C2")

df3$Treatment <- factor(df3$Treatment, levels = c("Unt","CD122-IL2cx", "CD25-IL2cx"))
# 


comparisons <- list(
  c("Unt", "CD25-IL2cx"),
  c("Unt", "CD122-IL2cx"),
  c("CD122-IL2cx", "CD25-IL2cx")
)



df3$measurement <- as.numeric(df3$measurement)

if (!exists(paste0(path.out, "Singleplots_gMFI"))) {
  dir.create(paste0(path.out, "Singleplots_gMFI"))
}

y_axis_settings <- scale_y_continuous(expand = expansion(mult = c(0, 0.05)), limits = c(0, NA))
plot_margin <- margin(15, 10, 20, 10)

unique_short <- unique(df3$short)
unique_marker <- unique(df3$marker)
unique_type <- unique(df3$type)

set.seed(21)
for (t in unique_type) {
  for (x in unique_marker) {
    int <- df3 %>% filter(type == t, marker == x)

    # Add needed cols for plot
    int$shape_group <- ifelse(int$Treatment == "Unt", "Unt", "Other")
    int$stroke_size <- ifelse(int$Treatment == "Unt", 0.75, 0.5)

    significant_valuesANOVA <- int %>%
      tukey_hsd(measurement ~ Treatment) %>%
      filter(paste(group1, group2, sep = " - ") %in% unlist(lapply(comparisons, function(x) paste(x, collapse = " - ")))) %>%
      add_xy_position()

    significant_valuesANOVA <- significant_valuesANOVA %>%
      filter(p.adj.signif != "ns")

    significant_valuesANOVA$p.adj <- as.character(significant_valuesANOVA$p.adj)
    significant_valuesANOVA$p.adj <- ifelse(as.numeric(significant_valuesANOVA$p.adj) < 0.0001, "<0.0001",
                                            format(as.numeric(significant_valuesANOVA$p.adj), digits = 5))

    significant_valuesANOVA$p.adj <- as.numeric(significant_valuesANOVA$p.adj)
    significant_valuesANOVA$p.adj <- format(significant_valuesANOVA$p.adj, scientific = FALSE)
    significant_valuesANOVA$p.adj[is.na(significant_valuesANOVA$p.adj)] <- "<0.0001"
    significant_valuesANOVA$p.adj <- as.character(significant_valuesANOVA$p.adj)
    significant_valuesANOVA$p.adj <- ifelse(as.numeric(significant_valuesANOVA$p.adj) < 0.0001, "<0.0001",
                                            format(as.numeric(significant_valuesANOVA$p.adj), digits = 4))
    significant_valuesANOVA$p.adj[is.na(significant_valuesANOVA$p.adj)] <- "<0.0001"
    significant_valuesANOVA$p.adj <- sub("\\.?0+$", "", significant_valuesANOVA$p.adj)

    # Set fixed y.position if you want all p-values to appear at the same vertical height
    max_y <- max(int$measurement, na.rm = TRUE)
    significant_valuesANOVA$y.position <- max_y * 1.05  # 5% above the max
    
    if (nrow(significant_valuesANOVA) > 0) {
  max_y <- max(int$measurement, na.rm = TRUE)
  spacing <- (max_y * 0.05)  # adjust vertical spacing here
  significant_valuesANOVA <- significant_valuesANOVA %>%
    arrange(p.adj) %>%
    mutate(y.position = max_y + (row_number() * spacing))
}
    
    
    

    p2 <- ggplot(int %>% filter(organ == "Spleen"), aes(x = Treatment, y = measurement)) +
      stat_summary(fun = mean, geom = "col", aes(fill = Treatment), na.rm = TRUE, colour = "black", size = 0.5, alpha = 0.4) +
      geom_point(aes(fill = Treatment, shape = shape_group), 
                 color = "black", size = 2,
                 stroke = int$stroke_size,
                 position = position_jitterdodge(jitter.width = 1, jitter.height = 0, dodge.width = 0)) +
      stat_summary(fun = mean, geom = "crossbar", width = 0.5, linewidth = 0.4) +
      theme_bw() +
      theme(axis.title.x = element_blank(),
            strip.background = element_blank()) +
      scale_color_manual(values = colvec) +
      scale_shape_manual(values = c("Unt" = 1, "Other" = 21)) +
      labs(y = "Frequency") +
      stat_pvalue_manual(significant_valuesANOVA, tip.length = 0, label = "p.adj") +
      theme(text = element_text(size = 4),
            plot.margin = plot_margin) +
      theme_prism() +
      theme(legend.position = "none",
            legend.title = element_text(colour = "black", size = 10, face = "bold")) +
      labs(color = 'Treatment') +
      facet_wrap("marker", scales = "free_y") +
      expand_limits(y = 0) +
      annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf) +
      annotate("segment", x = -Inf, xend = -Inf, y = -Inf, yend = Inf) +
      theme(axis.line = element_line()) +
      guides(color = guide_legend(override.aes = list(size = 5))) +
      theme(legend.title = element_text(size = 14)) +
      guides(fill = "none") +
      scale_x_discrete(name = "") +
      labs(x = " ") +
      theme(axis.title.x = element_blank(),
            axis.text.x = element_blank(),
            axis.ticks.x = element_blank()) +
      stat_summary(geom = "errorbar",
                   fun.data = mean_se,
                   fun.min = function(x) base::mean(x) - std.error(x),
                   fun.max = function(x) base::mean(x) + std.error(x),
                   width = 0.3) +
      y_axis_settings +
        # coord_cartesian(ylim = c(0, 99))+
      scale_fill_manual(values = colvec)

    tiff(paste0(path.out, "Singleplots_gMFI/Singleplots_gMFI", t, "_", x, ".tif"), height = 2, width = 2, res = 300, units = "in")
    print(p2)
    dev.off()

    pdf(paste0(path.out, "Singleplots_gMFI/Singleplots_gMFI_", t, "_", x, ".pdf"), onefile = FALSE, height = 2, width = 2)
    print(p2)
    dev.off()
  }
}

p2


makes a lot of plots for each unique type that I have in my dataframe. But I want to generate one big plot, where all of those unique plots are merged in one graph, so they share x and y axis so that I can make an overview on how those unique types compare to each other