sanitize.italics <- function(str) {
  str_replace_all(str, c("_" = " ", 
                         # "ital\\{" = "\\\\textit{", 
                         "°" = "\\\\degree",
                         "±"="$\\pm$"))
}


# See: https://www.rdocumentation.org/packages/xtable/versions/1.8-4/topics/print.xtable

latex_table <- function(table, caption, label, file, align, 
                        scalebox = getOption("xtable.scalebox", NULL), 
                        digits = 0, 
                        add.to.row = getOption("xtable.add.to.row", NULL),
                        include.colnames = TRUE, 
                        tabular.environment = getOption("xtable.tabular.environment", "tabular") ) {
  
    table <- xtable::xtable(table,
                            label=label,
                            caption=caption,
                            align = align,
                            digits=digits)
    print(table,  scalebox = scalebox,
          caption.placement = "top",
          include.rownames = FALSE,
          add.to.row = add.to.row,
          include.colnames = include.colnames, 
          tabular.environment = tabular.environment,
          file=file,
          sanitize.text.function = sanitize.italics)
    }

full_file_name <- function(file_name) {
    str_c(table_path, file_name)
}

latex_table(table = table_primers, caption = "Type of primers listed in the pr2-primers database. General primers target all eukaryotes and specific primers only certain taxonomic groups.", 
    label = "tab:primers", file = full_file_name("table_primers.tex"), align = c("l", "c", "c", "c"), scalebox = 1, digits = 0)

latex_table(table = table_primer_sets, caption = "Regions of the 18S rRNA gene targeted by the primer sets from the pr2-primers database.", label = "tab:primer_sets", 
    file = full_file_name("table_primer_sets.tex"), align = c("l", "l", "c", "c"), scalebox = 1, digits = 0)

latex_table(table = table_primer_sets_ampli, caption = "Overall characteristics of primer sets listed in the pr2-primers database.", label = "tab:primer_sets_ampli", 
    file = full_file_name("table_primer_sets_ampli.tex"), align = c("l", "l", "c", "c"), scalebox = 1, digits = 1)

ggplot(filter(pr2_match_final, kingdom == "Eukaryota"), aes(x = primer_set_label_short, y = ampli_size, group = primer_set_id, 
    color = as.factor(mismatch_number))) + geom_point(size = 3) + theme(axis.text.x = element_text(angle = 45, 
    hjust = 1)) + labs(x = "Primer set") + scale_color_viridis_d() + facet_wrap(vars(specific), nrow = 2, 
    scales = "free_x")

ggplot(filter(pr2_match_final, kingdom == "Eukaryota"), aes(x = primer_set_label_short, y = ampli_size, group = primer_set_id, 
    fill = as.factor(mismatch_number))) + geom_boxplot(outlier.alpha = 0.3) + theme(axis.text.x = element_text(angle = 45, 
    hjust = 1)) + labs(x = "Primer set") + scale_fill_viridis_d() + facet_wrap(vars(specific), nrow = 2, 
    scales = "free_x")

for (one_primer_set in c(8, 27)) {
    
    if (one_primer_set == 8) {
        xmin = 400
        xmax = 450
        xmax2 = 2000
    } else {
        xmin = 140
        xmax = 190
        xmax2 = 1000
        
    }
    
    pr2_match_final_one <- pr2_match_final %>% filter(primer_set_id == one_primer_set) %>% filter(kingdom == 
        "Eukaryota") %>% filter(!(supergroup %in% c("Apusozoa", "Eukaryota_X", "Protoalveoalata")))
    
    primer_set_label_long = pr2_match_final_one$primer_set_label_long[1]
    
    
    g <- ggplot(pr2_match_final_one, aes(x = ampli_size)) + geom_density(fill = "blue", alpha = 0.9) + xlab("Amplicon size") + 
        ggtitle(str_c("Primer set - ", primer_set_label_long)) + xlim(xmin, xmax)
    
    print(g)
    
    g <- ggplot(pr2_match_final_one, aes(x = ampli_size, fill = supergroup)) + geom_density(alpha = 0.9) + 
        theme_bw() + # theme(legend.position = 'none') +
    guides(fill = guide_legend(nrow = 3)) + theme(legend.position = "top", legend.box = "horizontal") + scale_fill_viridis_d(option = "inferno") + 
        labs(x = "Amplicon size (bp)", y = "Density", fill = "Supergroup") + # ggtitle(str_c('Primer set - ', primer_set_label_long)) +
    xlim(xmin, xmax)
    
    print(g)
    fig3[[str_c(one_primer_set, " size_distri")]] <- g
    
    g <- ggplot(pr2_match_final_one, aes(x = supergroup, y = ampli_size)) + geom_boxplot(outlier.alpha = 0.3) + 
        theme_bw() + coord_flip() + # theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    xlab("Supergroup") + ylab("Amplicon size (bp)") + ylim(0, xmax2) + geom_hline(yintercept = c(450, 550), 
        linetype = c(2, 3))
    
    print(g)
    fig3[[str_c(one_primer_set, " size")]] <- g
    
    g <- ggplot(pr2_match_final_one, aes(x = supergroup, y = ampli_size)) + geom_violin() + coord_flip() + 
        # theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    xlab("Supergroup") + ylab("Amplicon size (bp)")
    print(g)
    
}

for (one_primer_set in c(8, 27)) {
    
    pr2_match_summary_filtered <- pr2_match_summary_primer_set_sg %>% filter((n_seq > 20) & (primer_set_id == one_primer_set) & (kingdom == "Eukaryota") & 
        !(supergroup %in% c("Apusozoa", "Eukaryota_X")))
    
    
    
    g <- ggplot(pr2_match_summary_filtered) + geom_col(data = pr2_match_summary_filtered, aes(x = str_c(supergroup, " - n = ", n_seq), y = ampli_pct, 
        fill = supergroup), position = "dodge") + theme_bw() + # theme(legend.position = 'none') +
    guides(fill = guide_legend(nrow = 2)) + theme(legend.position = "top", legend.box = "horizontal") + scale_fill_viridis_d(option = "inferno") + 
        coord_flip() + labs(x = "% of sequences amplified", y = "", legend = "Supergroup") + scale_y_continuous(minor_breaks = seq(0, 100, by = 10), 
        breaks = seq(0, 100, by = 20), limits = c(0, 100))
    
    # ggtitle (str_c('Set - ', one_primer_set, ' - % amplified per Supergroup') )
    
    print(g)
    fig3[[str_c(one_primer_set, " pct")]] <- g
    
    g <- ggplot(filter(pr2_match_summary_filtered, !is.nan(ampli_size_mean))) + geom_point(aes(x = str_c(supergroup, " - n = ", n_seq), y = ampli_size_mean), 
        colour = "black") + coord_flip() + geom_errorbar(aes(x = str_c(supergroup, " - n = ", n_seq), ymax = ampli_size_mean + ampli_size_sd, ymin = ampli_size_mean - 
        ampli_size_sd)) + xlab("Supergroup") + ylab("Amplicon size (bp)")  # +
    # scale_y_continuous(breaks = (1:8)*200, limits = c(0,1500)) + ggtitle (str_c('Set - ', one_primer_set, ' - Amplicon sizes') ) +
    # geom_hline(yintercept = c(450,550) , linetype= 2)
    
    print(g)
    
    
    
    
}

for (one_primer_set in c(8, 27)) {
    
    pr2_mismatches <- pr2_match_summary_primer_set_sg %>% tidyr::pivot_longer(names_to = "mismatch_number", values_to = "mismatch_pct", cols = contains("ampli_mismatch"), 
        names_prefix = "ampli_mismatch_") %>% select(-contains("ampli_mismatch")) %>% mutate(mismatch_number = str_sub(mismatch_number, 1, 1)) %>% 
        mutate(mismatch_number = str_replace(mismatch_number, "5", "5+")) %>% filter((n_seq > 20) & (primer_set_id == one_primer_set) & (kingdom == 
        "Eukaryota") & !(supergroup %in% c("Apusozoa", "Eukaryota_X")))
    
    g <- ggplot(pr2_mismatches, aes(x = str_c(supergroup, " - n = ", n_seq), y = mismatch_pct, group = primer_set_id, fill = as.factor(mismatch_number))) + 
        geom_col() + theme_bw() + theme(axis.text.y = element_text(angle = 0, hjust = 0, vjust = 0)) + labs(x = "Supergroup", y = "% of sequences with mismatches", 
        fill = "Mismatches") + scale_fill_viridis_d(direction = -1, alpha = 0.85) + guides(fill = guide_legend(nrow = 1)) + theme(legend.position = "top", 
        legend.box = "horizontal") + coord_flip()
    
    print(g)
    
    fig3[[str_c(one_primer_set, "mismatches", sep = " ")]] <- g
    
    
    print(g)
    
    
}

for (one_primer_set in c(8, 27)) {
    
    df <- pr2_match_final %>% filter(primer_set_id == one_primer_set) %>% filter(kingdom == "Eukaryota") %>% filter(!(supergroup %in% c("Apusozoa", 
        "Eukaryota_X")))
    
    g <- ggplot() + geom_histogram(data = filter(df, !is.na(fwd_mismatch_primer_position_5prime)), aes(x = fwd_mismatch_primer_position_5prime, fill = supergroup), 
        binwidth = 1, stat = "bin", alpha = 1) + scale_fill_viridis_d(option = "inferno") + theme_bw() + scale_x_continuous(minor_breaks = 0:30) + 
        labs(y = "Density", x = "Position of mismatches from 5' end", title = "fwd")
    
    print(g)
    fig3[[str_c(one_primer_set, "mismatches positions fwd", sep = " ")]] <- g
    
    g <- ggplot() + geom_histogram(data = filter(df, !is.na(rev_mismatch_primer_position_5prime)), aes(x = rev_mismatch_primer_position_5prime, fill = supergroup), 
        binwidth = 1, stat = "bin", alpha = 1) + scale_fill_viridis_d(option = "inferno") + theme_bw() + scale_x_continuous(minor_breaks = 0:30) + 
        labs(y = "Density", x = "Position of mismatches from 5' end", title = "rev")
    
    print(g)
    fig3[[str_c(one_primer_set, "mismatches positions rev", sep = " ")]] <- g
    
}

for (one_kingdom in kingdoms) {
    
    for (one_region in gene_regions) {
        
        pr2_match_summary_primer_set_region_long <- pr2_match_summary_primer_set_long %>% filter(kingdom == one_kingdom) %>% filter(gene_region == 
            one_region) %>% # Remove the group specific primers
        filter(specific == "general") %>% filter(pct_category %in% c("ampli_pct", "fwd_pct", "rev_pct"))
        
        pr2_match_summary_primer_set_region <- pr2_match_summary_primer_set %>% filter(kingdom == one_kingdom) %>% filter(gene_region == one_region) %>% 
            # Remove the group specific primers)
        filter(specific == "general")
        
        pr2_match_region <- pr2_match_final %>% filter(kingdom == one_kingdom) %>% filter(gene_region == one_region) %>% # Remove the group specific primers
        filter(specific == "general")
        
        # % Ampli
        
        g <- ggplot(pr2_match_summary_primer_set_region_long) + geom_col(aes(x = str_replace(str_c(sprintf("%03d", primer_set_id), primer_set_label_long, 
            sep = " - "), " - general", ""), y = pct_seq, fill = fct_reorder(pct_category, pct_category_order)), width = 0.7, position = "dodge") + 
            theme_bw() + xlab("Primer set") + ylab("% of sequences amplified") + scale_fill_manual(name = "", values = c(ampli_pct = "black", fwd_pct = "blue", 
            rev_pct = "red"), labels = c("Amplicons", "rev", "fwd")) + ggtitle(str_c(one_kingdom, one_region, "Percentage of sequences recovered", 
            sep = " - ")) + theme(axis.text.x = element_text(angle = 0, hjust = 1)) + theme(axis.text.y = element_text(angle = 0, hjust = 0, vjust = 0)) + 
            ylim(0, 100) + coord_flip() + guides(fill = guide_legend(nrow = 1)) + theme(legend.position = "top", legend.box = "horizontal")
        
        print(g)
        
        # Size
        
        fig1[[str_c(one_kingdom, one_region, "pct", sep = " ")]] <- g
        
        
        g <- ggplot(filter(pr2_match_summary_primer_set_region, !is.nan(ampli_size_mean))) + geom_point(aes(x = str_c(sprintf("%03d", primer_set_id), 
            primer_set_label_long, sep = " - "), y = ampli_size_mean), colour = "black") + geom_errorbar(aes(x = str_c(sprintf("%03d", primer_set_id), 
            primer_set_label_long, sep = " - "), ymax = ampli_size_mean + ampli_size_sd, ymin = ampli_size_mean - ampli_size_sd)) + theme_bw() + xlab("Primer set") + 
            ylab("Amplicon size (bp)") + # scale_y_continuous(breaks = (1:8)*200, limits = c(0,1500)) +
        ggtitle(str_c(one_kingdom, one_region, "Amplicon size", sep = " - ")) + geom_hline(yintercept = c(450, 550), linetype = c(2, 3)) + ylim(0, 
            850) + theme(axis.text.x = element_text(angle = 0, hjust = 1)) + theme(axis.text.y = element_text(angle = 0, hjust = 0, vjust = 0)) + coord_flip()
        
        print(g)
        
        fig1[[str_c(one_kingdom, one_region, "size", sep = " ")]] <- g
        
        # Mismatches
        
        pr2_mismatches <- pr2_match_summary_primer_set %>% filter(kingdom == one_kingdom) %>% tidyr::pivot_longer(names_to = "mismatch_number", values_to = "mismatch_pct", 
            cols = contains("ampli_mismatch"), names_prefix = "ampli_mismatch_") %>% select(-contains("ampli_mismatch")) %>% mutate(mismatch_number = str_sub(mismatch_number, 
            1, 1)) %>% mutate(mismatch_number = str_replace(mismatch_number, "5", "5+")) %>% filter(gene_region == one_region) %>% # Remove the group specific primers
        filter(specific == "general")
        
        g <- ggplot(pr2_mismatches, aes(x = str_c(sprintf("%03d", primer_set_id), primer_set_label_long, sep = " - "), y = mismatch_pct, group = primer_set_id, 
            fill = as.factor(mismatch_number))) + geom_col() + theme_bw() + theme(axis.text.y = element_text(angle = 0, hjust = 0, vjust = 0)) + labs(x = "Primer set", 
            y = "% of sequences with mismatches", fill = "Mismatches", title = str_c(one_kingdom, one_region, sep = " - ")) + scale_fill_viridis_d(direction = -1, 
            alpha = 0.85) + guides(fill = guide_legend(nrow = 1)) + theme(legend.position = "top", legend.box = "horizontal") + coord_flip()
        
        print(g)
        
        fig1[[str_c(one_kingdom, one_region, "mismatches", sep = " ")]] <- g
        
        
    }
}

for (one_kingdom in kingdoms) {
    
    for (specific_one in c("general", "specific")) {
        
        df <- pr2_match_summary_primer_set_long %>% filter(kingdom == one_kingdom) %>% filter(pct_category %in% c("ampli_pct", "fwd_pct", "rev_pct")) %>% 
            filter(specific == specific_one)
        
        g <- ggplot(df) + geom_col(aes(x = str_replace(str_c(primer_set_label_long, sprintf("%03d", primer_set_id), sep = " - "), " - general", ""), 
            y = pct_seq, fill = fct_reorder(pct_category, pct_category_order)), width = 0.7, position = "dodge") + theme_bw() + xlab("Primer set") + 
            ylab("% of sequences amplified") + scale_fill_manual(name = "", values = c(ampli_pct = "black", fwd_pct = "blue", rev_pct = "red"), labels = c("Amplicons", 
            "rev", "fwd")) + ggtitle(str_c(one_kingdom, specific_one, "Percentage of sequences recovered", sep = " - ")) + theme(axis.text.y = element_text(angle = 0, 
            hjust = 0, vjust = 0)) + ylim(0, 100) + coord_flip() + guides(fill = guide_legend(nrow = 1)) + theme(legend.position = "top", legend.box = "horizontal")
        
        
        print(g)
        
        fig1[[str_c(one_kingdom, specific_one, "pct", sep = " ")]] <- g
        
        df <- pr2_match_summary_primer_set %>% filter(kingdom == one_kingdom) %>% filter(!is.nan(ampli_size_mean)) %>% filter(specific == specific_one)
        
        g <- ggplot(df) + geom_point(aes(x = str_c(primer_set_label_long, sprintf("%03d", primer_set_id), sep = " - "), y = ampli_size_mean), colour = "black") + 
            geom_errorbar(aes(x = str_c(primer_set_label_long, sprintf("%03d", primer_set_id), sep = " - "), ymax = ampli_size_mean + ampli_size_sd, 
                ymin = ampli_size_mean - ampli_size_sd)) + theme_bw() + xlab("Primer set") + ylab("Amplicon size (bp)") + # scale_y_continuous(breaks = (1:8)*200, limits = c(0,1500)) +
        ggtitle(str_c(one_kingdom, specific_one, " - Amplicon size", sep = " - ")) + geom_hline(yintercept = c(450, 550), linetype = c(2, 3)) + ylim(0, 
            900) + theme(axis.text.y = element_text(angle = 0, hjust = 0)) + coord_flip()
        
        print(g)
        
        fig1[[str_c(one_kingdom, specific_one, "size", sep = " ")]] <- g
        
    }
}

for (one_kingdom in kingdoms) {
    
    for (specific_one in c("general", "specific")) {
        
        pr2_mismatches <- pr2_match_summary_primer_set %>% tidyr::pivot_longer(names_to = "mismatch_number", values_to = "mismatch_pct", cols = contains("ampli_mismatch"), 
            names_prefix = "ampli_mismatch_") %>% select(-contains("ampli_mismatch")) %>% mutate(mismatch_number = str_sub(mismatch_number, 1, 1)) %>% 
            mutate(mismatch_number = str_replace(mismatch_number, "5", "5+")) %>% filter(specific == specific_one) %>% filter(kingdom == one_kingdom)
        
        g <- ggplot(pr2_mismatches, aes(x = str_c(primer_set_label_long, sprintf("%03d", primer_set_id), sep = " - "), y = mismatch_pct, group = primer_set_id, 
            fill = as.factor(mismatch_number))) + geom_col() + theme_bw() + theme(axis.text.y = element_text(angle = 0, hjust = 0, vjust = 0)) + labs(x = "Primer set", 
            y = "% of sequences with mismatches", fill = "Mismatches", title = str_c(one_kingdom, specific_one, sep = " - ")) + scale_fill_viridis_d(direction = -1, 
            alpha = 0.85) + guides(fill = guide_legend(nrow = 1)) + theme(legend.position = "top", legend.box = "horizontal") + coord_flip()
        
        print(g)
        
        fig1[[str_c(one_kingdom, specific_one, "mismatches", sep = " ")]] <- g
        
        
    }
}

ggplot(filter(pr2_match_final, kingdom == "Eukaryota"), aes(x = primer_set_label_long)) + geom_boxplot(aes(y = fwd_mismatch_primer_position_5prime), 
    outlier.alpha = 0.3, color = "blue") + coord_flip() + theme_bw() + theme(axis.text.y = element_text(angle = 0, vjust = 0, hjust = 0)) + scale_y_continuous(minor_breaks = 0:25, 
    breaks = seq(0, 25, by = 5), limits = c(0, 25)) + labs(x = "Primer set", y = "Position of mismatches", title = "Primer fwd") + facet_wrap(vars(specific), 
    nrow = 2, scales = "free_y")

ggplot(filter(pr2_match_final, kingdom == "Eukaryota"), aes(x = primer_set_label_long)) + geom_boxplot(aes(y = rev_mismatch_primer_position_5prime), 
    outlier.alpha = 0.3, color = "red") + coord_flip() + theme_bw() + theme(axis.text.y = element_text(angle = 0, vjust = 0, hjust = 0)) + scale_y_continuous(minor_breaks = 0:25, 
    breaks = seq(0, 25, by = 5), limits = c(0, 25)) + labs(x = "Primer set", y = "Position of mismatches", title = "Primer rev") + facet_wrap(vars(specific), 
    nrow = 2, scales = "free_y")

ggplot(filter(pr2_match_final, !is.na(rev_mismatch_primer_position_5prime), kingdom == "Eukaryota")) + geom_histogram(aes(x = rev_mismatch_primer_position_5prime), 
    color = "red", binwidth = 1, stat = "density") + theme_bw() + scale_x_continuous(minor_breaks = 0:30, breaks = seq(0, 30, by = 5), limits = c(0, 
    30)) + labs(y = "Primer set", x = "Position of mismatches", title = "Primer rev") + facet_wrap(vars(primer_set_label_long), ncol = ncol)

fig_supergroup <- list()
ncol = 8

for (specific_one in c("general", "specific")) {
    
    pr2_match_summary_primer_set_sg_region <- pr2_match_summary_primer_set_sg %>% filter(kingdom == "Eukaryota") %>% 
        filter(specific == specific_one) %>% filter((n_seq > 20)) %>% filter(!(supergroup %in% c("Apusozoa", 
        "Eukaryota_X")))
    
    
    
    g <- ggplot(pr2_match_summary_primer_set_sg_region) + geom_col(aes(x = supergroup, y = ampli_pct, fill = supergroup), 
        position = "dodge") + theme_bw() + coord_flip() + ylab("% of sequences amplified") + xlab("Supergroup") + 
        ggtitle(str_c(specific_one, " primers")) + scale_fill_viridis_d(option = "inferno") + ylim(0, 100) + 
        facet_wrap(~primer_set_label_short, scales = "fixed", ncol = ncol) + guides(fill = guide_legend(nrow = 1)) + 
        theme(legend.position = "top", legend.box = "horizontal") + theme(legend.position = "none")
    
    print(g)
    
    list_label <- str_c("pct", specific_one, sep = " ")
    fig_supergroup[[list_label]] <- g
    
    
    
    g <- ggplot(filter(pr2_match_summary_primer_set_sg_region, !is.nan(ampli_size_mean))) + geom_point(aes(x = supergroup, 
        y = ampli_size_mean), colour = "black") + theme_bw() + coord_flip() + geom_errorbar(aes(x = supergroup, 
        ymax = ampli_size_mean + ampli_size_sd, ymin = ampli_size_mean - ampli_size_sd)) + xlab("Supergroup") + 
        ylab("Amplicon size (bp)") + ggtitle(str_c(specific_one, " primers")) + geom_hline(yintercept = c(450, 
        550), linetype = 2) + facet_wrap(~primer_set_label_short, scales = "fixed", ncol = ncol)
    
    print(g)
    
    list_label <- str_c("size", specific_one, sep = " ")
    
    fig_supergroup[[list_label]] <- g
    
}

ncol = 8

pr2_mismatches_sg <- pr2_match_summary_primer_set_sg %>% filter(kingdom == "Eukaryota") %>% tidyr::pivot_longer(names_to = "mismatch_number", values_to = "mismatch_pct", 
    cols = contains("ampli_mismatch"), names_prefix = "ampli_mismatch_") %>% select(-contains("ampli_mismatch")) %>% mutate(mismatch_number = str_sub(mismatch_number, 
    1, 1)) %>% mutate(mismatch_number = str_replace(mismatch_number, "5", "5+")) %>% filter((n_seq > 20)) %>% filter(!(supergroup %in% c("Apusozoa", 
    "Eukaryota_X")))

for (specific_one in c("general", "specific")) {
    
    pr2_mismatches_sg_one <- pr2_mismatches_sg %>% filter(specific == specific_one)
    
    g <- ggplot(pr2_mismatches_sg_one, aes(x = supergroup, y = mismatch_pct, fill = fct_rev(mismatch_number))) + geom_col() + theme_bw() + theme(axis.text.x = element_text(angle = 0, 
        hjust = 1)) + labs(x = "Group", y = "% of sequences with mismatches", title = str_c(specific_one, " primers"), fill = "Mismatches") + scale_fill_viridis_d(direction = 1, 
        alpha = 0.85) + guides(fill = guide_legend(nrow = 1)) + theme(legend.position = "top", legend.box = "horizontal") + coord_flip() + facet_wrap(vars(primer_set_label_short), 
        ncol = ncol)
    
    print(g)
    
    fig_supergroup[[str_c("mismatches", specific_one, sep = " ")]] <- g
}

fig_1 <- (fig1[["Eukaryota V4 pct"]] + ggtitle("") + ylab("") + fig1[["Eukaryota V4 mismatches"]] + ggtitle("") + xlab("") + ylab("") + theme(axis.text.y = element_blank()) + 
    fig1[["Eukaryota V4 size"]] + ggtitle("") + xlab("") + ylab("") + theme(axis.text.y = element_blank()))/(fig1[["Eukaryota V9 pct"]] + ggtitle("") + 
    theme(legend.position = "none") + fig1[["Eukaryota V9 mismatches"]] + ggtitle("") + theme(legend.position = "none") + xlab("") + theme(axis.text.y = element_blank()) + 
    fig1[["Eukaryota V9 size"]] + ggtitle("") + xlab("") + theme(axis.text.y = element_blank())) + plot_layout(heights = c(20, 4))
fig_1

ggsave(plot = fig_1, filename = "figs/fig_pct_sizes_mismatches_V4_V9_euk.pdf", width = 17, height = 12, scale = 2, units = "cm", useDingbats = FALSE)

fig_1_bact <- (fig1[["Bacteria V4 pct"]] + ggtitle("") + ylab("") + fig1[["Bacteria V4 mismatches"]] + ggtitle("") + xlab("") + ylab("") + theme(axis.text.y = element_blank()) + 
    fig1[["Bacteria V4 size"]] + ggtitle("") + xlab("") + ylab("") + theme(axis.text.y = element_blank()))
fig_1_bact

ggsave(plot = fig_1_bact, filename = "figs/fig_pct_sizes_mismatches_V4_bact.pdf", width = 17, height = 10, scale = 2, units = "cm", useDingbats = FALSE)

fig_1_arch <- (fig1[["Archaea V4 pct"]] + ggtitle("") + ylab("") + fig1[["Archaea V4 mismatches"]] + ggtitle("") + xlab("") + ylab("") + theme(axis.text.y = element_blank()) + 
    fig1[["Archaea V4 size"]] + ggtitle("") + xlab("") + ylab("") + theme(axis.text.y = element_blank()))
fig_1_arch

ggsave(plot = fig_1_arch, filename = "figs/fig_pct_sizes_mismatches_V4_arch.pdf", width = 17, height = 10, scale = 2, units = "cm", useDingbats = FALSE)