#' Plot amplicon sequence variants #' #' This function plots the amplicon variants per sample using a provided distance clustering table. The plot can be customized by specifying the sequence type and whether hierarchical clustering should be performed. The function also creates a custom color scheme based on the viridis color palette. #' #' @param distclust_table A tibble containing the distance clustering table (default: mytibble) #' @param sequence_type The type of sequence to consider (default: 'bp') #' @param hierarchical_clustering Whether to perform hierarchical clustering (default: FALSE) #' #' @return A ggplot object representing the variants per sample plot #' #' @import dplyr ggplot2 tidyr #' @importFrom cowplot get_legend plot_grid #' @importFrom vegan vegdist #' @importFrom viridisLite viridis #' #' @examples #' plot_variants_per_sample() #' #' @export plot_variants_per_sample <- function(distclust_table = mytibble, sequence_type = 'bp', hierarchical_clustering = FALSE){ # Make a custom color scheme that assigns viridis 0 (purple) to ref/major, # viridis 1 (yellow) to major/minor # and a gradient between 0.7 and 0.3 to the rest. custom_colors <- c(viridisLite::viridis(1, begin = 0), viridisLite::viridis(1, begin = 1), viridisLite::viridis( length(unique(distclust_table$seq_var))-2, option = 'D', begin = 0.7 , end = 0.3)) names(custom_colors) <- distclust_table$seq_var %>% unique() %>% sort() if(hierarchical_clustering){ for_vegan <- distclust_table %>% select(var_type, seqnames, rel_var_abundance) %>% pivot_wider(id_cols = var_type, names_from = seqnames, values_from = rel_var_abundance) # rearrange x_axis variable according to order in hierarchical cluster dat <- dendrogram_hclust(veganify_generic_wide_tbl(for_vegan), seed = 1) distclust_table <- dplyr::mutate( distclust_table, var_type = factor(var_type, levels = dat$labels$label)) } length_measure <- ifelse(sequence_type != 'aa', 'seqs_lengths', 'longest_aa_seq') legend_text_master <- function(){ element_text(size = 10, margin = margin(t = 0, r = 15, b = 0, l = 0, unit = 'pt')) } # plot pl <- ggplot(distclust_table, aes(x = .data[[length_measure]], y = var_type, fill = seq_var, alpha = rel_var_abundance)) + geom_col(position = position_dodge2(padding = 0)) + theme_classic() + facet_grid(cols = vars(clus_name), scales = 'free') + scale_fill_manual(values = custom_colors) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), axis.text.y = element_text(hjust = 0), strip.text.x.top = element_text(angle = 90, hjust = 0, vjust = 0.5), strip.background = element_blank(), axis.title.y = element_blank()) + guides(fill = guide_legend(title = 'Sequence variants', title.theme = legend_text_master(), label.theme = legend_text_master(), keywidth = 0.5, keyheight = 0.5), alpha = guide_legend(title = 'Prevalence', title.theme = legend_text_master(), label.theme = legend_text_master(), keywidth = 0.5, keyheight = 0.5)) if(hierarchical_clustering){ d_pl <- plot_dendrogram(distclust_table = distclust_table) # extract a legend that is laid out horizontally legend <- get_legend(pl + theme(legend.position = "bottom")# + #guides(fill = guide_legend(nrow = 6)) ) pl_p <- cowplot::plot_grid(d_pl, pl + theme(legend.position = "none"), align = "h", axis = "bt", rel_widths = c(1, 5), nrow = 1) pl <- cowplot::plot_grid(pl_p, legend, ncol = 1, rel_heights = c(1, .15)) + theme(plot.background = element_rect(fill = "white")) } return(pl) }