R/plot_functions.R

Defines functions plot_var_part_pq treemap_pq ridges_sam_pq ridges_pq tax_bar_pq diff_fct_diff_class upset_test_pq upset_pq iNEXT_pq SRS_curve_pq plot_tsne_pq tsne_pq multitax_bar_pq plot_tax_pq multi_biplot_pq biplot_pq rotl_pq summary_plot_pq ggbetween_pq hill_pq multiplot ggvenn_pq venn_pq sankey_pq circle_pq accu_samp_threshold accu_plot_balanced_modality accu_plot plot_mt

Documented in accu_plot accu_plot_balanced_modality accu_samp_threshold biplot_pq circle_pq diff_fct_diff_class ggbetween_pq ggvenn_pq hill_pq iNEXT_pq multi_biplot_pq multiplot multitax_bar_pq plot_mt plot_tax_pq plot_tsne_pq plot_var_part_pq ridges_pq ridges_sam_pq rotl_pq sankey_pq SRS_curve_pq summary_plot_pq tax_bar_pq treemap_pq tsne_pq upset_pq upset_test_pq venn_pq

################################################################################
#' Plot the result of a mt test [phyloseq::mt()]
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-maturing-blue" alt="lifecycle-maturing"></a>
#'
#' Graphical representation of mt test.
#'
#' @param mt (required) Result of a mt test from the function [phyloseq::mt()].
#' @param pval (default: 0.05) Choose the cut off p-value to plot taxa.
#' @param color_tax (default: "Class") A taxonomic level to color the points.
#' @param taxa (default: "Species") The taxonomic level you choose for x-positioning.
#' @author Adrien Taudière
#' @examples
#' \donttest{
# #'  Filter samples that don't have Time
#' data_fungi_mini2 <- subset_samples(data_fungi_mini, !is.na(Time))
#' res <- mt(data_fungi_mini2, "Time", method = "fdr", test = "f", B = 300)
#' plot_mt(res)
#' plot_mt(res, taxa = "Genus", color_tax = "Order")
#' }
#' @return a \code{\link[ggplot2]{ggplot}}2 plot of result of a mt test
#' @export
#' @seealso [phyloseq::mt()]

plot_mt <-
  function(mt = NULL, pval = 0.05, color_tax = "Class", taxa = "Species") {
    d <- mt[mt$plower < pval, ]
    d$tax_col <- factor(as.character(d[, color_tax]))
    d$tax_col[is.na(d$tax_col)] <- "unidentified"
    d$tax <- as.character(d[, taxa])
    d$tax[is.na(d$tax)] <- "unidentified"
    d$tax <-
      factor(
        d$tax,
        levels = unique(factor(as.character(d[, taxa]))[rev(order(d$teststat))])
      )

    p <-
      ggplot(d, aes(x = tax, y = teststat, color = tax_col)) +
      geom_point(size = 6) +
      theme(
        axis.text.x = element_text(
          angle = -90,
          hjust = 0,
          vjust = 0.5
        )
      )
    p
  }
################################################################################

################################################################################
#' Plot accumulation curves for \code{\link[phyloseq]{phyloseq-class}} object
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-maturing-blue" alt="lifecycle-maturing"></a>
#'
#' Note that as most bioinformatic pipeline discard singleton, accumulation curves from metabarcoding
#' cannot be interpreted in the same way as with conventional biodiversity sampling techniques.
#'
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor in `physeq@sam_data` used to plot
#'    different lines
#' @param add_nb_seq (default: TRUE, logical)
#' Either plot accumulation curves using sequences or using samples
#' @param step (Integer) distance among points calculated to plot lines. A
#'  low value give better plot but is more time consuming.
#'  Only used if `add_nb_seq` = TRUE.
#' @param by.fact (default: FALSE, logical)
#' First merge the OTU table by factor to plot only one line by factor
#' @param ci_col Color vector for confidence interval.
#'   Only use if `add_nb_seq` = FALSE.
#'   If `add_nb_seq` = TRUE, you can use ggplot to modify the plot.
#' @param col Color vector for lines. Only use if `add_nb_seq` = FALSE.
#'   If `add_nb_seq` = TRUE, you can use ggplot to modify the plot.
#' @param lwd  (default: 3) thickness for lines. Only use if `add_nb_seq` = FALSE.
#' @param leg (default: TRUE, logical) Plot legend or not. Only use if `add_nb_seq` = FALSE.
#' @param print_sam_names (default: FALSE, logical) Print samples names or not?
#'    Only use if `add_nb_seq` = TRUE.
#' @param ci (default: 2, integer) Confidence interval value used to multiply the
#'   standard error to plot confidence interval
#' @param ... Additional arguments passed on to \code{\link[ggplot2]{ggplot}}
#' if `add_nb_seq` = TRUE or to \code{\link{plot}} if `add_nb_seq` = FALSE
#'
#' @examples
#' \donttest{
#' data("GlobalPatterns", package = "phyloseq")
#' GP <- subset_taxa(GlobalPatterns, GlobalPatterns@tax_table[, 1] == "Archaea")
#' GP <- rarefy_pq(subset_samples_pq(GP, sample_sums(GP) > 3000), replace = TRUE)
#' p <- accu_plot(GP, "SampleType", add_nb_seq = TRUE, by.fact = TRUE, step = 10)
#' p <- accu_plot(GP, "SampleType", add_nb_seq = TRUE, step = 10)
#'
#' p + theme(legend.position = "none")
#'
#' p + xlim(c(0, 400))
#' }
#' @return A \code{\link[ggplot2]{ggplot}}2 plot representing the richness
#' accumulation plot if add_nb_seq = TRUE, else, if add_nb_seq = FALSE
#' return a base plot.
#'
#' @export
#' @author Adrien Taudière
#' @seealso \code{\link[vegan]{specaccum}} [accu_samp_threshold()]
accu_plot <-
  function(
    physeq,
    fact = NULL,
    add_nb_seq = TRUE,
    step = NULL,
    by.fact = FALSE,
    ci_col = NULL,
    col = NULL,
    lwd = 3,
    leg = TRUE,
    print_sam_names = FALSE,
    ci = 2,
    ...
  ) {
    if (!inherits(physeq, "phyloseq")) {
      stop("physeq must be a phyloseq object")
    }

    if (!is.null(fact) && nlevels(as.factor(physeq@sam_data[[fact]])) < 2) {
      stop(
        "The factor '",
        fact,
        "' must have at least two levels for accu_plot ",
        "(species accumulation curves require at least 2 groups)."
      )
    }

    if (!taxa_are_rows(physeq)) {
      physeq@otu_table <-
        otu_table(t(physeq@otu_table), taxa_are_rows = TRUE)
    }

    if (!add_nb_seq) {
      factor_interm <-
        eval(parse(text = paste("physeq@sam_data$", fact, sep = "")))
      factor_interm <- as.factor(factor_interm)

      physeq_accu <- as(t(physeq@otu_table), "matrix")
      physeq_accu[physeq_accu > 0] <- 1
      accu_all <- vegan::specaccum(physeq_accu)

      accu <- vector("list", nlevels(factor_interm))
      for (i in seq_along(levels(factor_interm))) {
        accu[[i]] <-
          vegan::specaccum(physeq_accu[
            factor_interm == levels(factor_interm)[i],
          ])
      }

      if (is.null(col)) {
        col <- funky_color(nlevels(factor_interm) + 1)
      }
      if (is.null(ci_col)) {
        transp <- function(col, alpha = 0.5) {
          res <-
            apply(grDevices::col2rgb(col), 2, function(c) {
              grDevices::rgb(c[1] / 255, c[2] / 255, c[3] / 255, alpha)
            })
          return(res)
        }
        ci_col <-
          transp(funky_color(nlevels(factor_interm) + 1), 0.3)
      }

      plot(
        accu_all,
        # ci_type = "poly",
        # ci_col = ci_col[1],
        col = col[1],
        lwd = lwd,
        # ci_lty = 0,
        xlab = "Sample",
        ...
      )

      for (i in seq_along(levels(factor_interm))) {
        graphics::lines(accu[[i]], col = col[i + 1], lwd = lwd)
      }
      if (leg) {
        graphics::legend(
          "bottomright",
          c("all", levels(factor_interm)),
          col = col,
          lty = 1,
          lwd = 3
        )
      }
    }

    if (add_nb_seq) {
      fact_interm <-
        as.factor(unlist(unclass(physeq@sam_data[, fact])[fact]))

      if (by.fact) {
        x <- apply(physeq@otu_table, 1, function(x) {
          tapply(x, fact_interm, sum)
        })
      } else {
        x <- t(physeq@otu_table)
      }

      tot <- rowSums(x)
      nr <- nrow(x)

      if (is.null(step)) {
        step <- round(max(tot) / 30, 0)
      }

      n_max <- seq(1, max(tot), by = step)
      out <- lapply(seq_len(nr), function(i) {
        n <- seq(1, tot[i], by = step)
        if (n[length(n)] != tot[i]) {
          n <- c(n, tot[i])
        }
        res_interm <-
          vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE)
        res <-
          cbind(as.matrix(res_interm)[1, ], as.matrix(res_interm)[2, ])
        return(res)
      })

      names(out) <- names(tot)

      df <- plyr::ldply(out, data.frame)

      cond <- vector(mode = "logical")
      for (i in seq_along(levels(as.factor(df$.id)))) {
        cond <- c(cond, 1:table(df$.id)[i])
      }

      df$x <- n_max[cond]

      if (by.fact) {
        df$fact <- df$.id
      } else {
        df$fact <-
          as.factor(unlist(unclass(physeq@sam_data[
            match(df$.id, sample_names(physeq)),
            fact
          ])[fact]))
      }

      df$ymin <- df$X1 - df$X2 * ci
      df$ymin[is.na(df$ymin)] <- df$X1[is.na(df$ymin)]
      df$ymax <- df$X1 + df$X2 * ci
      df$ymax[is.na(df$ymax)] <- df$X1[is.na(df$ymax)]
      dff <- data.frame(matrix(nrow = length(tot)))
      dff$xlab <- tapply(df$x, df$.id, max)
      dff$xlab <- dff$xlab + max(dff$xlab, na.rm = TRUE) / 20
      dff$ylab <- tapply(df$X1, df$.id, max)
      dff$.id <- names(dff$ylab)
      p <-
        ggplot(
          data = df,
          aes(
            x = x,
            y = X1,
            group = .id,
            col = fact
          )
        ) +
        geom_ribbon(
          aes(
            ymin = ymin,
            ymax = ymax,
            col = NULL,
            fill = fact
          ),
          alpha = 0.2
        ) +
        geom_line() +
        xlab("Number of sequences") +
        ylab("Number of OTUs (with standard error)")

      if (print_sam_names) {
        p +
          geom_text(
            data = dff,
            aes(
              x = xlab,
              y = ylab,
              label = .id,
              col = NULL
            )
          )
      } else {
        p
      }
      return(p)
    }
  }
################################################################################

################################################################################
#' Plot accumulation curves with balanced modality and depth rarefaction
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#'   This function (i) rarefy (equalize) the number of samples per modality of a
#'   factor and (ii) rarefy the number of sequences per sample (depth). The
#'   seed is set to 1:nperm. Thus, with exacly the same parameter, including
#'   nperm values, results must be identical.
#'
#' @inheritParams clean_pq
#' @param fact (required) The variable to rarefy. Must be present in
#'   the `sam_data` slot of the physeq object.
#' @param nperm (int) The number of permutations to perform.
#' @param step 	(int) distance among points calculated to plot lines.
#' A low value give better plot but is more time consuming.
#' @param by.fact (logical, default TRUE)
#' First merge the OTU table by factor to plot only one line by factor
#' @param progress_bar (logical, default TRUE) Do we print progress during
#'   the calculation?
#' @param quantile_prob (float, `[0:1]`) the value to compute the quantile.
#'   Minimum quantile is compute using 1-quantile_prob.
#' @param rarefy_by_sample_before_merging (logical, default TRUE):
#'    rarefy_by_sample_before_merging = FALSE is buggy for the moment.Please
#'    only use rarefy_by_sample_before_merging = TRUE
#' @param sample.size (int) A single integer value equal to the number of
#'   reads being simulated, also known as the depth. See
#'   [phyloseq::rarefy_even_depth()] and [rarefy_even_depth_pq()].
#' @param verbose (logical). If TRUE, print additional information.
#' @param ... Other params for be passed on to [accu_plot()] function
#'
#' @export
#' @author Adrien Taudière
#' @seealso [accu_plot()], [rarefy_sample_count_by_modality()], [phyloseq::rarefy_even_depth()]
#'
#' @return A ggplot2 plot representing the richness accumulation plot
#' @examples
#' \donttest{
#' data_fungi_woNA4Time <-
#'   subset_samples(data_fungi_mini, !is.na(Time))
#' data_fungi_woNA4Time@sam_data$Time <-
#'   paste0("time-", data_fungi_woNA4Time@sam_data$Time)
#' accu_plot_balanced_modality(data_fungi_woNA4Time, "Time", nperm = 3)
#'
#' data_fungi_woNA4Height <-
#'   subset_samples(data_fungi_mini, !is.na(Height))
#' accu_plot_balanced_modality(data_fungi_woNA4Height, "Height", nperm = 3)
#' }
accu_plot_balanced_modality <- function(
  physeq,
  fact,
  nperm = 99,
  step = 2000,
  by.fact = TRUE,
  progress_bar = TRUE,
  quantile_prob = 0.975,
  rarefy_by_sample_before_merging = TRUE,
  sample.size = 1000,
  verbose = FALSE,
  ...
) {
  if (nlevels(as.factor(physeq@sam_data[[fact]])) < 2) {
    stop(
      "The factor '",
      fact,
      "' must have at least two levels for ",
      "accu_plot_balanced_modality (balanced accumulation curves require ",
      "at least 2 groups)."
    )
  }

  if (rarefy_by_sample_before_merging) {
    p_for_dim <- accu_plot(
      rarefy_sample_count_by_modality(
        rarefy_even_depth_pq(
          physeq,
          rngseed = 1,
          sample_size = sample.size
        ),
        fact,
        rngseed = 1
      ),
      fact = fact,
      step = step,
      by.fact = by.fact
    )$data
  } else {
    p_for_dim <- accu_plot(
      physeq,
      fact = fact,
      step = step,
      by.fact = by.fact
    )$data
    dim_for_plist <-
      max(tapply(sample_sums(physeq), physeq@sam_data[[fact]], sum))
  }

  dim_for_plist <- dim(p_for_dim)
  plist <- array(dim = c(dim_for_plist[1], 5, nperm))

  if (progress_bar) {
    pb <- txtProgressBar(
      min = 0,
      max = nperm,
      style = 3,
      width = 50,
      char = "="
    )
  }
  for (i in 1:nperm) {
    if (rarefy_by_sample_before_merging) {
      plist[,, i] <-
        as.matrix(suppressWarnings(suppressMessages(
          accu_plot(
            rarefy_sample_count_by_modality(
              rarefy_even_depth_pq(
                physeq,
                rngseed = i,
                sample_size = sample.size
              ),
              fact,
              rngseed = i
            ),
            fact = fact,
            step = step,
            by.fact = by.fact,
            ...
          )
        ))$data[, c(2:4, 6, 7)])
    } else {
      res_interm <-
        as.matrix(suppressWarnings(suppressMessages(
          accu_plot(
            rarefy_sample_count_by_modality(
              physeq,
              fact,
              rngseed = i,
              verbose = verbose
            ),
            fact = fact,
            step = step,
            by.fact = by.fact,
            ...
          )
        ))$data[, c(2:4, 6, 7)])
      plist[seq_along(nrow(res_interm)), , i] <- res_interm
    }
    if (progress_bar) {
      setTxtProgressBar(pb, i)
    }
  }

  res_mean <- data.frame(apply(plist, 1:2, mean, na.rm = TRUE))
  colnames(res_mean) <- colnames(p_for_dim[, c(2:4, 6, 7)])
  res_mean$fact <- p_for_dim$.id

  res_mean$X1_lim1 <- apply(plist, 1:2, function(x) {
    quantile(x, probs = quantile_prob, na.rm = TRUE)
  })[, 1]
  res_mean$X1_lim2 <- apply(plist, 1:2, function(x) {
    quantile(x, probs = 1 - quantile_prob, na.rm = TRUE)
  })[, 1]

  res_mean$factor <- p_for_dim$.id

  res_mean <- res_mean |>
    dplyr::filter(!is.na(X1)) |>
    arrange(X1)

  p <- ggplot(res_mean, aes(x = x, y = X1, color = factor)) +
    geom_line(linewidth = 1.5) +
    geom_ribbon(
      aes(
        ymin = X1_lim1,
        ymax = X1_lim2,
        fill = factor
      ),
      alpha = 0.2,
      linetype = 2,
      linewidth = 0.2
    )
  return(p)
}
################################################################################

################################################################################
#' Compute the number of sequence to obtain a given proportion of ASV in
#'  accumulation curves
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Note that as most bioinformatic pipeline discard singleton, accumulation curves from metabarcoding
#' cannot be interpreted in the same way as with conventional biodiversity sampling techniques.
#'
#' @param res_accuplot the result of the function accu_plot()
#' @param threshold the proportion of ASV to obtain in each samples
#'
#' @return a value for each sample of the number of sequences needed
#'   to obtain `threshold` proportion of the ASV
#'
#' @examples
#' \donttest{
#' data("GlobalPatterns", package = "phyloseq")
#' GP <- subset_taxa(GlobalPatterns, GlobalPatterns@tax_table[, 1] == "Archaea")
#' GP <- rarefy_pq(subset_samples_pq(GP, sample_sums(GP) > 3000), replace = TRUE)
#' p <- accu_plot(GP, "SampleType", add_nb_seq = TRUE, by.fact = TRUE, step = 10)
#'
#' val_threshold <- accu_samp_threshold(p)
#'
#' summary(val_threshold)
#'
#' ##'  Plot the number of sequences needed to accumulate 0.95% of ASV in 50%, 75%
#' ##'  and 100% of samples
#' p + geom_vline(xintercept = quantile(val_threshold, probs = c(0.50, 0.75, 1)))
#' }
#' @export
#' @author Adrien Taudière
#' @seealso [accu_plot()]
accu_samp_threshold <- function(res_accuplot, threshold = 0.95) {
  res <- vector("list", length(unique(res_accuplot$data$.id)))
  names(res) <- unique(res_accuplot$data$.id)
  for (id in unique(res_accuplot$data$.id)) {
    data <- res_accuplot$data |> dplyr::filter(.id == id)
    proportion <- data$X1 / max(data$X1)
    res[[id]] <- data$x[proportion > threshold][1]
  }
  return(unlist(res))
}


################################################################################

################################################################################
#' Plot OTU circle for \code{\link[phyloseq]{phyloseq-class}} object
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-maturing-blue" alt="lifecycle-maturing"></a>
#'
#' Graphical representation of distribution of taxa across a factor.
#'
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor to cluster samples by modalities.
#'        Need to be in \code{physeq@sam_data}.
#' @param taxa (default: 'Order') Name of the taxonomic rank of interest
#' @param nproc (default 1)
#'   Set to number of cpus/processors to use for parallelization
#' @param add_nb_seq (logical, default TRUE) Represent the number of sequences or the
#'    number of OTUs (add_nb_seq = FALSE)
#' @param rarefy (logical) Does each samples modalities need to be rarefy in
#'               order to compare them with the same amount of sequences?
#' @param min_prop_tax (default: 0.01) The minimum proportion for taxa to be
#'                     plotted
#' @param min_prop_mod (default: 0.1) The minimum proportion for modalities
#'                     to be plotted
#' @param gap_degree Gap between two neighbour sectors.
#'                    It can be a single value or a vector. If it is a vector,
#'                     the first value corresponds to the gap after the first
#'                     sector.
#' @param start_degree The starting degree from which the circle begins to
#'   draw. Note this degree is measured in the standard polar coordinate
#'   which means it is always reverse-clockwise.
#' @param row_col Color vector for row
#' @param grid_col Grid colors which correspond to sectors. The length of the
#'    vector should be either 1 or the number of sectors.
#'    It's preferred that grid_col is a named vector of which names
#'    correspond to sectors. If it is not a named vector, the
#'    order of grid_col corresponds to order of sectors.
#' @param log10trans (logical) Should sequence be log10 transformed
#'                   (more precisely by log10(1+x))?
#' @param ... Additional arguments passed on to
#'   \code{\link[circlize]{chordDiagram}} or \code{\link[circlize]{circos.par}}
#'
#' @examples
#' \donttest{
#' data("GlobalPatterns", package = "phyloseq")
#' GP <- subset_taxa(GlobalPatterns, GlobalPatterns@tax_table[, 1] == "Archaea")
#' circle_pq(GP, "SampleType")
#' }
#' \dontrun{
#' circle_pq(GP, "SampleType", add_nb_seq = FALSE)
#' circle_pq(GP, "SampleType", taxa = "Class")
#' }
#' @author Adrien Taudière
#'
#' @return A \code{\link[circlize]{chordDiagram}} plot representing the
#'   distribution of OTUs or sequences in the different modalities of the factor
#'   fact
#'
#' @export
#' @seealso \code{\link[circlize]{chordDiagram}}
#' @seealso \code{\link[circlize]{circos.par}}

circle_pq <-
  function(
    physeq = NULL,
    fact = NULL,
    taxa = "Order",
    nproc = 1,
    add_nb_seq = TRUE,
    rarefy = FALSE,
    min_prop_tax = 0.01,
    min_prop_mod = 0.1,
    gap_degree = NULL,
    start_degree = NULL,
    row_col = NULL,
    grid_col = NULL,
    log10trans = FALSE,
    ...
  ) {
    if (!inherits(physeq, "phyloseq")) {
      stop("physeq must be an object of class 'phyloseq'")
    }

    if (physeq@otu_table@taxa_are_rows) {
      otu_tab <- physeq@otu_table
    } else {
      otu_tab <- t(physeq@otu_table)
    }

    if (!add_nb_seq) {
      otu_tab[otu_tab > 0] <- 1
    }

    taxcol <- match(taxa, colnames(physeq@tax_table))
    if (is.na(taxcol)) {
      stop("The taxa argument do not match any taxa rank in physeq@tax_table")
    }

    taxsamp <- match(fact, colnames(physeq@sam_data))
    if (is.na(taxsamp)) {
      stop(
        "The samples argument do not match any sample attributes
           in physeq@sam_data"
      )
    }

    tax_groups <- as.character(physeq@tax_table[, taxcol])
    tax_groups[is.na(tax_groups)] <- "NA"
    otu_table_tax <- rowsum(
      as.matrix(otu_tab),
      group = tax_groups,
      na.rm = TRUE,
      reorder = TRUE
    )

    sam_groups <- as.character(unlist(physeq@sam_data[, taxsamp]))
    sam_groups[is.na(sam_groups)] <- "NA"
    otu_table_ech <- rowsum(
      t(otu_table_tax),
      group = sam_groups,
      na.rm = TRUE,
      reorder = TRUE
    )
    if (!is.matrix(otu_table_ech)) {
      otu_table_ech <- matrix(
        otu_table_ech,
        nrow = 1,
        dimnames = list(
          levels(as.factor(
            unlist(unclass(physeq@sam_data[, fact]))
          )),
          names(otu_table_ech)
        )
      )
    }
    if (rarefy) {
      otu_table_ech_interm <-
        vegan::rrarefy(otu_table_ech, min(rowSums(otu_table_ech)))
      message(
        paste(
          "Rarefaction by modalities deletes ",
          sum(otu_table_ech) - sum(otu_table_ech_interm),
          " (",
          round(
            100 *
              (sum(otu_table_ech) - sum(otu_table_ech_interm)) /
              sum(otu_table_ech),
            2
          ),
          "%) sequences.",
          sep = ""
        )
      )
      otu_table_ech <- otu_table_ech_interm
    }

    otu_table_ech <- otu_table_ech[, colSums(otu_table_ech) > 0, drop = FALSE]

    # Keep only taxa and modalities with a sufficient proportion (min_prop_tax,
    # min_prop_mod) to plot
    o_t_e_interm <-
      otu_table_ech[
        (rowSums(otu_table_ech) / sum(otu_table_ech)) > min_prop_mod,
        (colSums(otu_table_ech) / sum(otu_table_ech)) > min_prop_tax,
        drop = FALSE
      ]
    if (nrow(o_t_e_interm) != nrow(otu_table_ech)) {
      message(
        paste(
          "Only ",
          nrow(o_t_e_interm),
          " modalities are plot (",
          round(
            100 *
              nrow(o_t_e_interm) /
              nrow(otu_table_ech),
            2
          ),
          "%). Use 'min_prop_mod' to plot more samples.",
          sep = ""
        )
      )
    }

    if (ncol(o_t_e_interm) != ncol(otu_table_ech)) {
      message(
        paste(
          "Only ",
          ncol(o_t_e_interm),
          " taxa are plot (",
          round(
            100 * ncol(o_t_e_interm) / ncol(otu_table_ech),
            2
          ),
          "%). Use 'min_prop_tax' to plot more taxa",
          sep = ""
        )
      )
    }
    otu_table_ech <- o_t_e_interm

    if (log10trans) {
      otu_table_ech <- apply(otu_table_ech, 2, function(x) {
        log10(1 + x)
      })
    }

    if (is.null(gap_degree)) {
      col2keep <- rep(1, ncol(otu_table_ech) - 1)
      row2keep <- rep(1, nrow(otu_table_ech) - 1)
      gap_degree <- c(row2keep, 10, col2keep, 10)
    }
    if (is.null(start_degree)) {
      start_degree <- 170
    }

    if (is.null(grid_col)) {
      grid_col <-
        c(funky_color(nrow(otu_table_ech)), rep("grey", ncol(otu_table_ech)))
    }

    if (is.null(row_col)) {
      row_col <-
        c(funky_color(nrow(otu_table_ech)), rep("grey", ncol(otu_table_ech)))
    }

    circlize::circos.par(
      gap.degree = gap_degree,
      start.degree = start_degree,
      ...
    )
    circlize::chordDiagram(
      otu_table_ech,
      row.col = row_col,
      grid.col = grid_col,
      ...
    )
    circlize::circos.clear()
  }
################################################################################

################################################################################
#' Sankey plot of \code{\link[phyloseq]{phyloseq-class}} object
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-maturing-blue" alt="lifecycle-maturing"></a>
#'
#' Graphical representation of distribution of taxa across Taxonomy and (optionnaly a factor).
#'
#' @inheritParams clean_pq
#' @param fact Name of the factor to cluster samples by modalities.
#' Need to be in \code{physeq@sam_data}.
#' @param taxa a vector of taxonomic rank to plot
#' @param add_nb_seq Represent the number of sequences or the
#'   number of OTUs (add_nb_seq = FALSE). Note that plotting the number of
#'   sequences is slower.
#' @param min_prop_tax (default: 0) The minimum proportion for taxa to be
#'  plotted. EXPERIMENTAL. For the moment each links below the min.prop.
#'  tax is discard from the sankey network resulting in sometimes weird plot.
#' @param tax2remove  a vector of taxonomic groups to remove from the analysis
#'   (e.g. \code{c('Incertae sedis', 'unidentified')})
#' @param units  character string describing physical units (if any) for Value
#' @param symbol2sub (default: c('\\.', '-')) vector of symbol to delete in
#'   the taxonomy
#' @param ... Additional arguments passed on to
#'   \code{\link[networkD3]{sankeyNetwork}}
#'
#' @examples
#' data("GlobalPatterns", package = "phyloseq")
#' GP <- subset_taxa(GlobalPatterns, GlobalPatterns@tax_table[, 1] == "Archaea")
#' if (requireNamespace("networkD3")) {
#'   sankey_pq(GP, fact = "SampleType")
#' }
#' \donttest{
#' if (requireNamespace("networkD3")) {
#'   sankey_pq(GP, taxa = 1:4, min_prop_tax = 0.01)
#'   sankey_pq(GP, taxa = 1:4, min_prop_tax = 0.01, add_nb_seq = TRUE)
#' }
#' }
#' @author Adrien Taudière
#'
#' @return A \code{\link[networkD3]{sankeyNetwork}} plot representing the
#'  taxonomic distribution of OTUs or sequences. If \code{fact} is set,
#'  represent the distribution of the last taxonomic level in the modalities
#'  of \code{fact}
#'
#' @export
#' @seealso \code{\link[networkD3]{sankeyNetwork}}, [ggaluv_pq()]

sankey_pq <-
  function(
    physeq = NULL,
    fact = NULL,
    taxa = 1:4,
    add_nb_seq = FALSE,
    min_prop_tax = 0,
    tax2remove = NULL,
    units = NULL,
    symbol2sub = c("\\.", "-"),
    ...
  ) {
    if (!inherits(physeq, "phyloseq")) {
      stop("physeq must be an object of class 'phyloseq'")
    }

    if (physeq@otu_table@taxa_are_rows) {
      otu_tab <- physeq@otu_table
    } else {
      otu_tab <- t(physeq@otu_table)
    }

    if (!add_nb_seq) {
      otu_tab[otu_tab > 0] <- 1
      mat_list <- vector("list", length(taxa) - 1)
      for (i in seq_len(length(taxa) - 1)) {
        res_interm <-
          table(physeq@tax_table[, taxa[i]], physeq@tax_table[, taxa[i + 1]])
        mat_interm <- reshape2::melt(res_interm)
        mat_list[[i]] <- mat_interm[mat_interm[, 3] > 0, ]
      }
      mat <- do.call(rbind, mat_list)
      colnames(mat) <- c("Var1", "Var2", "value")
    } else if (add_nb_seq) {
      mat_list <- vector("list", length(taxa) - 1)
      tax_table_interm <-
        physeq@tax_table[rep(seq_len(ntaxa(physeq)), times = taxa_sums(physeq))]

      for (i in seq_len(length(taxa) - 1)) {
        res_interm <-
          table(tax_table_interm[, taxa[i]], tax_table_interm[, taxa[i + 1]])
        mat_interm <- reshape2::melt(res_interm)
        mat_list[[i]] <- mat_interm[mat_interm[, 3] > 0, ]
      }
      mat <- do.call(rbind, mat_list)
      colnames(mat) <- c("Var1", "Var2", "value")
    }

    if (!is.null(fact)) {
      net_matrix2links <- function(m = NULL) {
        # Pre-calculate dimensions and non-zero positions for efficiency
        dims <- dim(m)
        rows <- row(m)
        cols <- col(m)
        mask <- m > 0

        if (!any(mask)) {
          return(matrix(ncol = 3)[0, ]) # Return empty matrix with correct structure
        }

        # Vectorized approach - much more efficient than nested loops
        row_names <- rownames(m)[rows[mask]]
        col_names <- colnames(m)[cols[mask]]
        values <- m[mask]

        result <- cbind(row_names, col_names, values)
        colnames(result) <- c("Var1", "Var2", "value")
        return(result)
      }

      mat_interm <-
        apply(otu_tab, 1, function(x) {
          tapply(
            x,
            physeq@sam_data[, fact],
            sum
          )
        })
      if (!is.matrix(mat_interm)) {
        mat_interm <- matrix(
          mat_interm,
          nrow = 1,
          dimnames = list(
            levels(as.factor(
              unlist(unclass(physeq@sam_data[, fact]))
            )),
            names(mat_interm)
          )
        )
      }

      if (!add_nb_seq) {
        mat_interm <-
          apply(mat_interm, 1, function(x) {
            tapply(
              x,
              physeq@tax_table[,
                taxa[length(taxa)]
              ],
              function(x) {
                sum(x > 0)
              }
            )
          })
      } else if (add_nb_seq) {
        mat_interm <-
          apply(mat_interm, 1, function(x) {
            tapply(
              x,
              physeq@tax_table[,
                taxa[length(taxa)]
              ],
              sum
            )
          })
      }
      if (!is.matrix(mat_interm)) {
        mat_interm <- matrix(
          mat_interm,
          ncol = 1,
          dimnames = list(names(mat_interm), colnames(mat_interm))
        )
      }

      samp_links <- net_matrix2links(mat_interm)
      samp_links[, 2] <- toupper(samp_links[, 2])
      colnames(samp_links) <- colnames(mat)
      mat <- rbind(mat, samp_links)
    }

    mat <- as.data.frame(mat[rowSums(is.na(mat)) == 0, ])
    mat[, 3] <- as.numeric(as.vector(mat[, 3]))
    mat <- mat[rowSums(is.na(mat)) == 0, ]

    if (!is.null(tax2remove)) {
      mat <- mat[!mat[, 1] %in% tax2remove, ]
      mat <- mat[!mat[, 2] %in% tax2remove, ]
    }

    if (min_prop_tax != 0) {
      min_nb_tax <- min_prop_tax * sum(mat[, 3]) / length(taxa)
      mat <- mat[mat[, 3] >= min_nb_tax, ]
    }

    for (i in seq_along(symbol2sub)) {
      mat <- apply(mat, 2, function(x) {
        gsub(symbol2sub[i], "", x)
      })
    }

    tax_sank <- vector("list", 2)
    names_nodes <-
      unique(c(as.vector(mat[, 1]), as.vector(mat[, 2])))
    names_nodes <- names_nodes[!is.na(names_nodes)]
    tax_sank$nodes <-
      data.frame((seq_along(names_nodes)) - 1, names_nodes)
    names(tax_sank$nodes) <- c("code", "name")
    mat2 <- mat
    for (i in seq_len(nrow(tax_sank$nodes))) {
      mat2[, 1] <-
        gsub(
          paste("\\<", tax_sank$nodes[i, 2], "\\>", sep = ""),
          tax_sank$nodes[
            i,
            1
          ],
          mat2[, 1]
        )
      mat2[, 2] <-
        gsub(
          paste("\\<", tax_sank$nodes[i, 2], "\\>", sep = ""),
          tax_sank$nodes[
            i,
            1
          ],
          mat2[, 2]
        )
    }

    tax_sank$links <- apply(mat2, 2, as.numeric)
    tax_sank$links <-
      data.frame(tax_sank$links[rowSums(is.na(tax_sank$links)) == 0, ])
    tax_sank$nodes <-
      as.data.frame(as.character(tax_sank$nodes[, 2]))
    names(tax_sank$nodes) <- "name"
    names(tax_sank$links) <- c("source", "target", "value")
    if (is.null(units)) {
      if (!add_nb_seq) {
        units <- "OTUs"
      } else if (add_nb_seq) {
        units <- "Sequences"
      }
    }
    networkD3::sankeyNetwork(
      Links = tax_sank$links,
      Nodes = tax_sank$nodes,
      Source = "source",
      Target = "target",
      Value = "value",
      NodeID = "name",
      units = units,
      ...
    )
  }
################################################################################

################################################################################
#' Venn diagram of \code{\link[phyloseq]{phyloseq-class}} object
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-maturing-blue" alt="lifecycle-maturing"></a>
#'
#' Graphical representation of distribution of taxa across combined modality of a factor.
#'
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor to cluster samples by modalities.
#' Need to be in \code{physeq@sam_data}.
#' @param min_nb_seq (default: 0) minimum number of sequences by OTUs by
#'  samples to take into count this OTUs in this sample. For example,
#'  if min_nb_seq=2,each value of 2 or less in the OTU table
#'  will be change into 0 for the analysis
#' @param print_values (logical) Print (or not) the table of number of OTUs
#' for each combination.
#' If print_values is TRUE the object is not a ggplot object.
#' Please use print_values = FALSE if you want to add ggplot function
#' (cf example).
#'
#' @examplesIf tolower(Sys.info()[["sysname"]]) != "windows"
#' if (requireNamespace("venneuler")) {
#'   data("enterotype")
#'   venn_pq(enterotype, fact = "SeqTech")
#' }
#' \donttest{
#' if (requireNamespace("venneuler")) {
#'   venn_pq(enterotype, fact = "ClinicalStatus")
#'   venn_pq(enterotype, fact = "Nationality", print_values = FALSE)
#'   venn_pq(enterotype, fact = "ClinicalStatus", print_values = FALSE) +
#'     scale_fill_hue()
#'   venn_pq(enterotype, fact = "ClinicalStatus", print_values = FALSE) +
#'     scale_fill_hue()
#' }
#' }
#' @return A \code{\link[ggplot2]{ggplot}}2 plot representing Venn diagram of
#' modalities of the argument \code{factor}
#'
#' @export
#' @author Adrien Taudière
#' @seealso \code{\link[venneuler]{venneuler}}

venn_pq <-
  function(physeq, fact, min_nb_seq = 0, print_values = TRUE) {
    if (!inherits(physeq, "phyloseq")) {
      stop("physeq must be an object of class 'phyloseq'")
    }

    moda <-
      as.factor(unlist(unclass(physeq@sam_data[, fact])[fact]))

    if (nlevels(moda) < 2) {
      stop(
        "The factor '",
        fact,
        "' must have at least two levels for venn_pq ",
        "(Venn diagrams require at least 2 sets)."
      )
    }
    if (length(moda) != dim(physeq@otu_table)[1]) {
      data_venn <-
        t(apply(physeq@otu_table, 1, function(x) {
          by(x, moda, max)
        }))
    } else if (length(moda) != dim(physeq@otu_table)[2]) {
      data_venn <-
        t(apply(t(physeq@otu_table), 1, function(x) {
          by(x, moda, max)
        }))
    } else {
      stop("The factor length and the number of samples must be identical")
    }
    combinations <- data_venn > min_nb_seq

    e <- new.env(TRUE, emptyenv())
    cn <- colnames(combinations)
    for (i in seq.int(dim(combinations)[1])) {
      if (any(combinations[i, ])) {
        ec <- paste(cn[combinations[i, ]], collapse = "&")
        e[[ec]] <- if (is.null(e[[ec]])) {
          1L
        } else {
          (e[[ec]] + 1L)
        }
      }
    }

    en <- ls(e, all.names = TRUE)
    weights <- as.numeric(unlist(lapply(en, get, e)))
    combinations <- as.character(en)

    table_value <-
      data.frame(
        combinations = as.character(combinations),
        weights = as.double(weights),
        stringsAsFactors = FALSE
      )

    venn <- venneuler::venneuler(data_venn > min_nb_seq)
    venn_res <-
      data.frame(
        x = venn$centers[, 1],
        y = venn$centers[, 2],
        radius = venn$diameters / 2
      )

    nmod <- nrow(venn_res)
    x1 <- vector("list", nmod)
    for (i in seq(1, nmod)) {
      x1[[i]] <- grep(rownames(venn_res)[i], table_value$combinations)
    }

    for (i in seq_len(nrow(table_value))) {
      table_value$x[i] <-
        mean(venn$centers[, "x"][unlist(lapply(
          x1,
          function(x) {
            sum(x %in% i) > 0
          }
        ))])
      table_value$y[i] <-
        mean(venn$centers[, "y"][unlist(lapply(
          x1,
          function(x) {
            sum(x %in% i) > 0
          }
        ))])
    }

    df <- venn_res
    df$xlab <- df$x + (df$x - mean(df$x))
    df$ylab <- df$y + (df$y - mean(df$y))

    circularise <- function(d, n = 360) {
      angle <- seq(-pi, pi, length = n)
      make_circle <- function(x, y, r, modality) {
        data.frame(
          x = x + r * cos(angle),
          y = y + r * sin(angle),
          modality
        )
      }
      lmat <- mapply(
        make_circle,
        modality = rownames(d),
        x = d[, 1],
        y = d[, 2],
        r = d[, 3],
        SIMPLIFY = FALSE
      )
      do.call(rbind, lmat)
    }

    circles <- circularise(df)

    p <-
      ggplot() +
      geom_polygon(
        data = circles,
        aes(x, y, group = modality, fill = modality),
        alpha = 0.5
      ) +
      theme_void()

    if (print_values) {
      g_legend <- function(agplot) {
        tmp <- ggplot_gtable(ggplot_build(agplot))
        leg <-
          which(
            vapply(
              tmp$grobs,
              function(x) {
                x$name
              },
              character(1)
            ) ==
              "guide-box"
          )
        legend <- tmp$grobs[[leg]]
        return(legend)
      }
      legend <- g_legend(p)

      grid::grid.newpage()
      vp1 <- grid::viewport(
        width = 0.75,
        height = 1,
        x = 0.375,
        y = 0.5
      )
      vpleg <-
        grid::viewport(
          width = 0.25,
          height = 0.5,
          x = 0.85,
          y = 0.75
        )
      subvp <- grid::viewport(
        width = 0.3,
        height = 0.3,
        x = 0.85,
        y = 0.25
      )
      print(p + theme(legend.position = "none"), vp = vp1)
      grid::upViewport(0)
      grid::pushViewport(vpleg)
      grid::grid.draw(legend)
      grid::upViewport(0)
      grid::pushViewport(subvp)
      grid::grid.draw(gridExtra::tableGrob(table_value[, c(1, 2)], rows = NULL))
    } else {
      return(p)
    }
  }
################################################################################

################################################################################
#' Venn diagram of \code{\link[phyloseq]{phyloseq-class}} object using
#' `ggVennDiagram::ggVennDiagram` function
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-maturing-blue" alt="lifecycle-maturing"></a>
#'
#' Note that you can use ggplot2 function to customize the plot
#' for ex. `+ scale_fill_distiller(palette = "BuPu", direction = 1)`
#' and `+ scale_x_continuous(expand = expansion(mult = 0.5))`. See
#' examples.
#'
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor to cluster samples by modalities.
#'   Need to be in \code{physeq@sam_data}.
#' @param min_nb_seq minimum number of sequences by OTUs by
#'   samples to take into count this OTUs in this sample. For example,
#'   if min_nb_seq=2,each value of 2 or less in the OTU table
#'   will not count in the venn diagram
#' @param taxonomic_rank Name (or number) of a taxonomic rank
#'   to count. If set to Null (the default) the number of OTUs is counted.
#' @param split_by Split into multiple plot using variable split_by.
#'   The name of a variable must be present in `sam_data` slot
#'   of the physeq object.
#' @param add_nb_samples (logical, default TRUE) Add the number of samples to
#'    levels names
#' @param add_nb_seq (logical, default FALSE) Add the number of sequences to
#'    levels names
#' @param rarefy_before_merging Rarefy each sample before merging by the
#'   modalities of args `fact`. Use `phyloseq::rarefy_even_depth()` function
#' @param rarefy_after_merging Rarefy each sample after merging by the
#'   modalities of args `fact`.
#' @param rngseed (Optional). A single integer value passed to
#'   [phyloseq::rarefy_even_depth()], which is used to fix a seed for
#'   reproducibly random number generation (in this case, reproducibly
#'   random subsampling). If set to FALSE, then no fiddling with the RNG seed
#'   is performed, and it is up to the user to appropriately call set.seed
#'   beforehand to achieve reproducible results. Default is FALSE.
#' @param return_data_for_venn (logical, default FALSE) If TRUE, the plot is
#'   not returned, but the resulting dataframe to plot with ggVennDiagram package
#'   is returned.
#' @param verbose (logical, default TRUE) If TRUE, prompt some messages.
#' @param type If "nb_taxa" (default), the number of taxa (ASV, OTU or
#'   taxonomic_rank if `taxonomic_rank` is not NULL) is
#'   used in plot. If "nb_seq", the number of sequences is plotted.
#'   `taxonomic_rank` is never used if type = "nb_seq".
#' @param na_remove (logical, default TRUE) If set to TRUE, remove samples with
#'   NA in the variables set in `fact` param
#' @param ... Other arguments for the `ggVennDiagram::ggVennDiagram` function
#'   for ex. `category.names`.
#' @return A \code{\link[ggplot2]{ggplot}}2 plot representing Venn diagram of
#'   modalities of the argument \code{factor} or if split_by is set a list
#'   of plots.
#' @seealso [upset_pq()]
#' @examples
#' if (requireNamespace("ggVennDiagram")) {
#'   ggvenn_pq(data_fungi_mini, fact = "Height")
#' }
#' \donttest{
#' if (requireNamespace("ggVennDiagram")) {
#'   ggvenn_pq(data_fungi_mini, fact = "Height") +
#'     ggplot2::scale_fill_distiller(palette = "BuPu", direction = 1)
#'   pl <- ggvenn_pq(data_fungi_mini, fact = "Height", split_by = "Time")
#'   for (i in seq_along(pl)) {
#'     p <- pl[[i]] +
#'       scale_fill_distiller(palette = "BuPu", direction = 1) +
#'       theme(plot.title = element_text(hjust = 0.5, size = 22))
#'     print(p)
#'   }
#'
#'   data_fungi2 <- subset_samples(
#'     data_fungi_mini,
#'     data_fungi_mini@sam_data$Tree_name == "A10-005" |
#'       data_fungi_mini@sam_data$Height %in% c("Low", "High")
#'   )
#'   ggvenn_pq(data_fungi2, fact = "Height")
#'
#'   ggvenn_pq(data_fungi2, fact = "Height", type = "nb_seq")
#'
#'   ggvenn_pq(data_fungi_mini, fact = "Height", add_nb_seq = TRUE, set_size = 4)
#'   ggvenn_pq(data_fungi_mini, fact = "Height", rarefy_before_merging = TRUE)
#'   ggvenn_pq(data_fungi_mini, fact = "Height", rarefy_after_merging = TRUE) +
#'     scale_x_continuous(expand = expansion(mult = 0.5))
#'
#'   # For more flexibility, you can save the dataset for more precise construction
#'   # with ggplot2 and ggVennDiagramm
#'   # (https://gaospecial.github.io/ggVennDiagram/articles/fully-customed.html)
#'   res_venn <- ggvenn_pq(data_fungi_mini,
#'     fact = "Height",
#'     return_data_for_venn = TRUE
#'   )
#'
#'   ggplot() +
#'     # 1. region count layer
#'     geom_polygon(aes(X, Y, group = id, fill = name),
#'       data = ggVennDiagram::venn_regionedge(res_venn)
#'     ) +
#'     scale_fill_manual(values = funky_color(7)) +
#'     # 2. set edge layer
#'     geom_path(aes(X, Y, color = id, group = id),
#'       data = ggVennDiagram::venn_setedge(res_venn),
#'       show.legend = FALSE, linewidth = 2
#'     ) +
#'     scale_color_manual(values = c("red", "red", "blue")) +
#'     # 3. set label layer
#'     geom_text(aes(X, Y, label = name),
#'       data = ggVennDiagram::venn_setlabel(res_venn)
#'     ) +
#'     # 4. region label layer
#'     geom_label(
#'       aes(X, Y, label = paste0(
#'         count, " (",
#'         scales::percent(count / sum(count), accuracy = 2), ")"
#'       )),
#'       data = ggVennDiagram::venn_regionlabel(res_venn)
#'     ) +
#'     theme_void()
#' }
#' }
#' @export
#' @author Adrien Taudière

ggvenn_pq <- function(
  physeq = NULL,
  fact = NULL,
  min_nb_seq = 0,
  taxonomic_rank = NULL,
  split_by = NULL,
  add_nb_samples = TRUE,
  add_nb_seq = FALSE,
  rarefy_before_merging = FALSE,
  rarefy_after_merging = FALSE,
  rngseed = FALSE,
  return_data_for_venn = FALSE,
  verbose = TRUE,
  type = "nb_taxa",
  na_remove = TRUE,
  ...
) {
  if (!is.factor(physeq@sam_data[[fact]])) {
    physeq@sam_data[[fact]] <- as.factor(physeq@sam_data[[fact]])
  }

  if (nlevels(physeq@sam_data[[fact]]) < 2) {
    stop(
      "The factor '",
      fact,
      "' must have at least two levels for ggvenn_pq ",
      "(Venn diagrams require at least 2 sets)."
    )
  }

  if (na_remove) {
    new_physeq <- subset_samples_pq(physeq, !is.na(physeq@sam_data[[fact]]))
    if (nsamples(physeq) - nsamples(new_physeq) > 0 && verbose) {
      message(paste0(
        nsamples(physeq) - nsamples(new_physeq),
        " were discarded due to NA in variable fact"
      ))
    }
    physeq <- new_physeq
  }

  physeq <- taxa_as_columns(physeq)

  if (rarefy_before_merging) {
    if (as(rngseed, "logical")) {
      set.seed(rngseed)
      if (verbose) {
        message(
          "`set.seed(",
          rngseed,
          ")` was used to initialize repeatable random subsampling."
        )
        message("Please record this for your records so others can reproduce.")
        message(
          "Try `set.seed(",
          rngseed,
          "); .Random.seed` for the full vector",
          sep = ""
        )
        message("...")
      }
    } else if (verbose) {
      message(
        "You set `rngseed` to FALSE. Make sure you've set & recorded\n",
        " the random seed of your session for reproducibility.\n",
        "See `?set.seed`\n"
      )
      message("...")
    }
    physeq <- rarefy_even_depth_pq(physeq, rngseed = rngseed)
    physeq <- clean_pq(physeq)
  }

  nb_samples <- table(physeq@sam_data[[fact]])

  if (rarefy_after_merging) {
    physeq <- merge_samples2(physeq, fact)
    if (as(rngseed, "logical")) {
      set.seed(rngseed)
      if (verbose) {
        message(
          "`set.seed(",
          rngseed,
          ")` was used to initialize repeatable random subsampling."
        )
        message("Please record this for your records so others can reproduce.")
        message(
          "Try `set.seed(",
          rngseed,
          "); .Random.seed` for the full vector",
          sep = ""
        )
        message("...")
      }
    } else if (verbose) {
      message(
        "You set `rngseed` to FALSE. Make sure you've set & recorded\n",
        " the random seed of your session for reproducibility.\n",
        "See `?set.seed`\n"
      )
      message("...")
    }
    physeq <- rarefy_even_depth_pq(physeq, rngseed = rngseed)
    physeq <- clean_pq(physeq)
  }

  res <- vector("list", nlevels(physeq@sam_data[[fact]]))
  names(res) <- levels(physeq@sam_data[[fact]])
  nb_seq <- vector(mode = "integer")

  for (f in levels(physeq@sam_data[[fact]])) {
    newphyseq <- physeq
    new_DF <- newphyseq@sam_data[
      newphyseq@sam_data[[fact]] == f,
      ,
      drop = FALSE
    ]
    sample_data(newphyseq) <- sample_data(new_DF)
    newphyseq <- clean_pq(newphyseq)
    if (is.null(taxonomic_rank) || type == "nb_seq") {
      res[[f]] <- colnames(newphyseq@otu_table[,
        colSums(newphyseq@otu_table) > min_nb_seq
      ])
    } else {
      res[[f]] <-
        as.character(stats::na.exclude(unique(newphyseq@tax_table[
          colSums(newphyseq@otu_table) > min_nb_seq,
          taxonomic_rank
        ])))
    }
    nb_seq <-
      c(
        nb_seq,
        sum(physeq@otu_table[physeq@sam_data[[fact]] == f, ], na.rm = TRUE)
      )

    if (type == "nb_seq") {
      res[[f]] <- unlist(sapply(res[[f]], function(x) {
        paste0(x, "_", seq(1, taxa_sums(physeq)[[x]]))
      }))
    }
  }

  if (max(nb_seq) / min(nb_seq) > 2 && verbose) {
    message(
      paste0(
        "Two modalities differ greatly (more than x2) in their number of sequences (",
        max(nb_seq),
        " vs ",
        min(nb_seq),
        "). You may be interested by the parameter rarefy_after_merging"
      )
    )
  }

  if (add_nb_samples) {
    names(res) <- paste0(names(res), "\n (", nb_samples, " sam.)")
  }

  if (add_nb_seq) {
    names(res) <- paste0(names(res), "\n (", nb_seq, " seq.)")
  }

  if (is.null(split_by)) {
    p <- ggVennDiagram::ggVennDiagram(res, ...)
  } else {
    modalities <-
      as.factor(unlist(unclass(physeq@sam_data[[split_by]])))
    p <- vector("list", nlevels(modalities))
    names(p) <- levels(modalities)
    for (moda in levels(modalities)) {
      physeq_interm <-
        clean_pq(subset_samples_pq(physeq, modalities == moda), silent = TRUE)
      p[[moda]] <- ggvenn_pq(
        physeq_interm,
        fact = fact,
        min_nb_seq = 0,
        taxonomic_rank = NULL
      ) +
        ggtitle(moda)
    }
  }
  if (return_data_for_venn) {
    return(ggVennDiagram::process_data(ggVennDiagram::Venn(res)))
  } else {
    return(p)
  }
}
################################################################################

################################################################################
#' Multiple plot function
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-stable-green" alt="lifecycle-stable"></a>
#'
#' ggplot objects can be passed in ..., or to plotlist (as a list of ggplot
#' objects)
#'
#' If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
#' then plot 1 will go in the upper left, 2 will go in the upper right, and
#' 3 will go all the way across the bottom.
#'
#' @param ... list of ggplot objects
#' @param plotlist list of ggplot objects
#' @param cols number of columns
#' @param layout A matrix specifying the layout.
#'   If present, 'cols' is ignored.
#' @return Nothing. Print the list of ggplot objects
#' @export

multiplot <-
  function(..., plotlist = NULL, cols = 1, layout = NULL) {
    # Make a list from the ... arguments and plotlist
    plots <- c(list(...), plotlist)

    num_plots <- length(plots)

    # If layout is NULL, then use 'cols' to determine layout
    if (is.null(layout)) {
      # Make the panel
      # ncol: Number of columns of plots
      # nrow: Number of rows needed, calculated from # of cols
      layout <- matrix(
        seq(1, cols * ceiling(num_plots / cols)),
        ncol = cols,
        nrow = ceiling(num_plots / cols)
      )
    }

    if (num_plots == 1) {
      message(plots[[1]])
    } else {
      # Set up the page
      grid::grid.newpage()
      grid::pushViewport(grid::viewport(
        layout = grid::grid.layout(
          nrow(layout),
          ncol(layout)
        )
      ))

      # Make each plot, in the correct location
      for (i in seq(1, num_plots)) {
        # Get the i,j matrix positions of the regions that contain this subplot
        matchidx <-
          as.data.frame(which(layout == i, arr.ind = TRUE))

        print(
          plots[[i]],
          vp = grid::viewport(
            layout.pos.row = matchidx$row,
            layout.pos.col = matchidx$col
          )
        )
      }
    }
  }
################################################################################

################################################################################
#' Graphical representation of hill number 0, 1 and 2 across a factor
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Hill numbers are the number of equiprobable species giving the same
#'   diversity value as the observed distribution. The Hill number 0
#'   correspond to Species richness), the Hill number 1 to
#'   the exponential of Shannon Index and the Hill number 2 to the inverse
#'   of Simpson Index)
#'
#' Note that (if correction_for_sample_size is TRUE, default behavior)
#'   this function use a sqrt of the read numbers in the linear
#'   model in order to correct for uneven sampling depth. This correction
#'   is only done before tuckey HSD plot and do not change the hill number
#'   computed.
#'
#' @inheritParams clean_pq
#' @param fact (required) The variable to test. Must be present in
#'   the `sam_data` slot of the physeq object.
#' @param variable : Alias for factor. Kept only for backward compatibility.
#' @param q (vector) Hill diversity orders to compute. Default computes
#'   Hill number 0 (species richness), 1 (exponential of Shannon index) and
#'   2 (inverse of Simpson index). Hill numbers are more appropriate in DNA
#'   metabarcoding studies when `q > 0` (Alberdi & Gilbert, 2019;
#'   Calderón-Sanou et al., 2019).
#' @param hill_scales `r lifecycle::badge("deprecated")` Use `q` instead.
#' @param ... Additional arguments passed to [divent_hill_matrix_pq()] and
#'   hence to [divent::div_hill()] (e.g. `estimator = "naive"`).
#' @param color_fac (optional): The variable to color the barplot. For ex.
#'   same as fact. Not very useful because ggplot2 plot colors can be
#'   change using `scale_color_XXX()` function.
#' @param letters (optional, default FALSE): If set to TRUE, the plot
#'   show letters based on p-values for comparison. Use the
#'   \code{\link[multcompView]{multcompLetters}} function from the package
#'   multcompLetters. BROKEN for the moment. Note that na values in The
#'   variable param need to be removed (see examples) to use letters.
#' @param add_points (logical, default FALSE): add jitter point on boxplot
#' @param add_info (logical, default TRUE) Do we add a subtitle with
#'   information about the number of samples per modality ?
#' @param one_plot (logical, default FALSE) If TRUE, return a unique
#'   plot with the four plot inside using the patchwork package.
#'   Note that if letters and one_plot are both TRUE, tuckey HSD results
#'   are discarded from the unique plot. In that case, use one_plot = FALSE
#'   to see the tuckey HSD results in the fourth plot of the resulting list.
#' @param kruskal_test (logical, default TRUE) Do we test for global effect of
#'   our factor on each hill scales values? When kruskal_test is TRUE, the
#'   resulting test value are add in each plot in subtitle (unless add_info is
#'   FALSE). Moreover, if at
#'   least one hill scales is not significantly link to fact (pval>0.05),
#'   a message is prompt saying that Tuckey HSD plot is not informative for
#'   those Hill scales and letters are not printed.
#' @param plot_with_tuckey (logical, default TRUE). If one_plot is set to
#'   TRUE and letters to FALSE, allow to discard the tuckey plot part with
#'   plot_with_tuckey = FALSE
#' @param correction_for_sample_size (logical, default TRUE) This function
#'   use a sqrt of the read numbers in the linear model in order to
#'   correct for uneven sampling depth in the Tuckey TEST. This params
#'   do not change value of Hill number but only the test associated
#'   values (including the pvalues). To rarefy samples, you may use the
#'   function [phyloseq::rarefy_even_depth()].
#' @param na_remove (logical, default TRUE) Do we remove samples with NA in
#'   the factor fact ? Note that na_remove is always TRUE when using
#'   letters = TRUE
#' @param vioplot (logical, default FALSE) Do we plot violin plot instead of
#'   boxplot ?
#' @return Either an unique ggplot2 object (if one_plot is TRUE) or
#'  a list of n+1 ggplot2 plot (with n the number of hill scale value).
#'  For example, with the default scale value:
#' - plot_Hill_0 : the boxplot of Hill number 0 (= species richness)
#'     against the variable
#' - plot_Hill_1 : the boxplot of Hill number 1 (= Shannon index)
#'      against the variable
#' - plot_Hill_2 : the boxplot of Hill number 2 (= Simpson index)
#'     against the variable
#' - plot_tuckey : plot the result of the Tuckey HSD test
#'
#' @export
#' @author Adrien Taudière
#' @examples
#'
#' data_f <- prune_samples(
#'   sample_names(data_fungi_mini)[1:20],
#'   data_fungi_mini
#' )
#' p <- hill_pq(data_f, "Height", q = c(0, 1))
#' p[[1]] + theme(legend.position = "none")
#' \dontrun{
#' if (requireNamespace("multcompView")) {
#'   p2 <- hill_pq(data_fungi_mini, "Time",
#'     correction_for_sample_size = FALSE,
#'     letters = TRUE, add_points = TRUE,
#'     plot_with_tuckey = FALSE
#'   )
#'   if (requireNamespace("patchwork")) {
#'     patchwork::wrap_plots(p2, guides = "collect")
#'   }
#'   p3 <- hill_pq(data_fungi_mini, "Height",
#'     letters = TRUE, vioplot = TRUE,
#'     add_points = TRUE
#'   )
#' }
#' }
#' @seealso [psmelt_samples_pq()] and [ggbetween_pq()]
#' @references
#' Alberdi, A., & Gilbert, M. T. P. (2019). A guide to the application of
#'   Hill numbers to DNA-based diversity analyses. *Molecular Ecology Resources*.
#'   \doi{10.1111/1755-0998.13014}
#'
#' Calderón-Sanou, I., Münkemüller, T., Boyer, F., Zinger, L., & Thuiller, W.
#'   (2019). From environmental DNA sequences to ecological conclusions: How
#'   strong is the influence of methodological choices? *Journal of Biogeography*,
#'   47. \doi{10.1111/jbi.13681}
hill_pq <- function(
  physeq,
  fact = NULL,
  variable = NULL,
  q = c(0, 1, 2),
  hill_scales = lifecycle::deprecated(),
  color_fac = NA,
  letters = FALSE,
  add_points = FALSE,
  add_info = TRUE,
  kruskal_test = TRUE,
  one_plot = FALSE,
  plot_with_tuckey = TRUE,
  correction_for_sample_size = TRUE,
  na_remove = TRUE,
  vioplot = FALSE,
  ...
) {
  if (lifecycle::is_present(hill_scales)) {
    lifecycle::deprecate_warn(
      "0.15.1",
      "hill_pq(hill_scales=)",
      "hill_pq(q=)"
    )
    q <- hill_scales
  }
  if (!is.null(variable)) {
    if (!is.null(fact)) {
      stop(
        "You must set only one parameter of variable or fact. This 2
        parameters are strictly equivalent."
      )
    } else {
      variable_fac <- variable
    }
  } else {
    if (!is.null(fact)) {
      variable_fac <- fact
    } else {
      stop("You must set the parameter fact.")
    }
  }
  var <- sym(variable_fac)
  if (is.na(color_fac)) {
    color_fac <- sym(variable_fac)
  } else {
    color_fac <- sym(color_fac)
  }

  physeq <- taxa_as_rows(physeq)
  if (na_remove || letters) {
    physeq <- subset_samples_pq(physeq, !is.na(physeq@sam_data[[fact]]))
  }
  physeq@sam_data[[fact]] <- as.factor(physeq@sam_data[[fact]])

  if (nlevels(physeq@sam_data[[fact]]) < 2) {
    stop(
      "The factor '",
      fact,
      "' must have at least two levels for hill_pq ",
      "(Kruskal-Wallis and Tukey tests require at least 2 groups)."
    )
  }

  otu_hill <- divent_hill_matrix_pq(
    as.data.frame(t(as.matrix(physeq@otu_table))),
    q = q,
    ...
  )
  colnames(otu_hill) <- paste0("Hill_", q)

  df_hill <- data.frame(otu_hill, physeq@sam_data)
  df_hill[, seq_along(q)] <-
    apply(df_hill[, seq_along(q)], 2, as.numeric)

  p_var <-
    hill_tuckey_pq(
      physeq,
      modality = variable_fac,
      q = q,
      correction_for_sample_size = correction_for_sample_size,
      ...
    )
  p_list <- vector("list", length(q))

  if (kruskal_test) {
    kt_res <- vector("list", length(q))
    for (i in seq_along(q)) {
      kt_res[[i]] <- kruskal.test(
        df_hill[, paste0("Hill_", q[[i]])],
        df_hill[, fact]
      )
    }
    if (
      sum(sapply(kt_res, function(x) {
        x$p.value > 0.05
      })) >
        0
    ) {
      message(paste0(
        sum(sapply(kt_res, function(x) {
          x$p.value > 0.05
        })),
        " out of ",
        length(kt_res),
        " Hill scales do not show any global trends with you factor ",
        fact,
        ". Tuckey HSD plot is not informative for those Hill scales. Letters are not printed for those Hill scales"
      ))
    }
  }

  for (i in seq_along(q)) {
    if (vioplot) {
      p_list[[i]] <-
        ggplot(
          df_hill,
          aes(
            x = .data[[paste0("Hill_", q[[i]])]],
            y = !!var
          )
        ) +
        geom_violin(aes(colour = as.factor(!!color_fac))) +
        labs(x = paste0("Hill_", q[[i]]))
    } else {
      p_list[[i]] <-
        ggplot(
          df_hill,
          aes(group = !!var, x = .data[[paste0("Hill_", q[[i]])]])
        ) +
        geom_boxplot(
          outlier.size = 2,
          aes(colour = as.factor(!!color_fac), y = !!var)
        ) +
        labs(x = paste0("Hill_", q[[i]]))
    }

    if (add_points) {
      p_list[[i]] <-
        p_list[[i]] +
        geom_jitter(
          aes(y = !!var, colour = as.factor(!!color_fac)),
          alpha = 0.5
        )
    }
    if (add_info) {
      subtitle_plot <- paste0(
        "Nb of samples: '",
        paste0(
          names(table(physeq@sam_data[[variable_fac]])),
          sep = "' : ",
          table(physeq@sam_data[[variable_fac]]),
          collapse = " - '"
        )
      )
      if (kruskal_test) {
        subtitle_plot <- paste0(
          subtitle_plot,
          "\n",
          paste0(
            " Hill ",
            q[[i]],
            " -- Kruskal-Wallis chi-squared =",
            round(kt_res[[i]]$statistic, 2),
            "; df = ",
            kt_res[[i]]$parameter,
            "; p.value =",
            format.pval(kt_res[[i]]$p.value, 2)
          )
        )
      }
      p_list[[i]] <- p_list[[i]] + labs(subtitle = subtitle_plot)
    }

    if (letters) {
      data_h <-
        p_var$data[grep(paste0("Hill_", q[[i]]), p_var$data[, 5]), ]
      data_h_pval <- data_h$`p adj`
      names(data_h_pval) <- data_h$modality
      Letters <-
        multcompView::multcompLetters(data_h_pval, reversed = TRUE)$Letters

      dt <- data.frame(variab = names(Letters), Letters = Letters)
      names(dt) <- c(var, "Letters")
      data_letters <- p_list[[i]]$data |>
        group_by(!!var) |>
        summarize(
          pos_letters = max(.data[[paste0("Hill_", q[[i]])]]) + 1
        ) |>
        inner_join(dt, by = join_by(!!fact))

      if (!kruskal_test || kt_res[[i]]$p.value < 0.05) {
        p_list[[i]] <- p_list[[i]] +
          geom_label(
            data = data_letters,
            aes(
              x = pos_letters,
              label = Letters,
            ),
            y = unique(ggplot_build(p_list[[i]])$data[[1]]$y),
            size = 4,
            stat = "unique",
            parse = TRUE
          )
      }
    }
  }

  res <- p_list
  names(res) <- paste0("plot_Hill_", q)
  if (plot_with_tuckey) {
    res[["tuckey"]] <- p_var
  }

  if (one_plot) {
    requireNamespace("patchwork", quietly = TRUE)
    if (letters || !plot_with_tuckey) {
      res[["tuckey"]] <- NULL
    }
    res <- patchwork::wrap_plots(res)
  }
  return(res)
}
################################################################################

################################################################################
#' Box/Violin plots for between-subjects comparisons of Hill Number
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Note that contrary to [hill_pq()], this function does not take into
#' account for difference in the number of sequences per samples/modalities.
#' You may use rarefy_by_sample = TRUE if the mean number of sequences per
#' samples differs among modalities.
#'
#' Basically a wrapper of function [ggstatsplot::ggbetweenstats()] for
#' object of class phyloseq
#' @inheritParams clean_pq
#' @param fact (required) The variable to test. Must be present in
#'   the `sam_data` slot of the physeq object.
#' @param one_plot (logical, default FALSE) If TRUE, return a unique
#'   plot with the three plot inside using the patchwork package.
#' @param rarefy_by_sample (logical, default FALSE) If TRUE, rarefy
#'   samples using [phyloseq::rarefy_even_depth()] function
#' @param rngseed (Optional). A single integer value passed to
#'   [phyloseq::rarefy_even_depth()], which is used to fix a seed for
#'   reproducibly random number generation (in this case, reproducibly
#'   random subsampling). If set to FALSE, then no fiddling with the RNG seed
#'   is performed, and it is up to the user to appropriately call set.seed
#'   beforehand to achieve reproducible results. Default is FALSE.
#' @param verbose (logical). If TRUE, print additional information.
#' @param q (numeric vector, default `c(0, 1, 2)`) Hill diversity orders to
#'   compute. One plot is produced per value. Hill numbers are more appropriate
#'   in DNA metabarcoding studies when `q > 0` (Alberdi & Gilbert, 2019;
#'   Calderón-Sanou et al., 2019).
#' @param ... Additional arguments passed on to [ggstatsplot::ggbetweenstats()] function.

#' @return Either an unique ggplot2 object (if one_plot is TRUE) or
#'  a list of ggplot2 plots, one per Hill order in `q`. With default `q`:
#' - plot_Hill_0 : the ggbetweenstats of Hill number 0 (= species richness)
#'     against the variable fact
#' - plot_Hill_1 : the ggbetweenstats of Hill number 1 (= Shannon index)
#'      against the variable fact
#' - plot_Hill_2 : the ggbetweenstats of Hill number 2 (= Simpson index)
#'     against the variable fact
#'
#' @export
#' @examples
#' \donttest{
#' library("divent")
#' if (requireNamespace("ggstatsplot")) {
#'   data_f <- clean_pq(prune_samples(
#'   sample_names(data_fungi_sp_known)[1:10],
#'   data_fungi_sp_known
#' ))
#'   p <- ggbetween_pq(data_f, fact = "Time", p.adjust.method = "BH")
#'   p[[1]]
#' }
#' }
#' \dontrun{
#' if (requireNamespace("ggstatsplot")) {
#'   ggbetween_pq(data_fungi, fact = "Height", one_plot = TRUE)
#'   ggbetween_pq(data_fungi, fact = "Height", one_plot = TRUE, rarefy_by_sample = TRUE)
#' }
#' }
#' @author Adrien Taudière
#' @details This function is mainly a wrapper of the work of others.
#'   Please make a reference to `ggstatsplot::ggbetweenstats()` if you
#'   use this function.
#' @references
#' Alberdi, A., & Gilbert, M. T. P. (2019). A guide to the application of
#'   Hill numbers to DNA-based diversity analyses. *Molecular Ecology Resources*.
#'   \doi{10.1111/1755-0998.13014}
#'
#' Calderón-Sanou, I., Münkemüller, T., Boyer, F., Zinger, L., & Thuiller, W.
#'   (2019). From environmental DNA sequences to ecological conclusions: How
#'   strong is the influence of methodological choices? *Journal of Biogeography*,
#'   47. \doi{10.1111/jbi.13681}

ggbetween_pq <-
  function(
    physeq,
    fact,
    one_plot = FALSE,
    rarefy_by_sample = FALSE,
    rngseed = FALSE,
    verbose = TRUE,
    q = c(0, 1, 2),
    ...
  ) {
    verify_pq(physeq)

    if (nlevels(as.factor(physeq@sam_data[[fact]])) < 2) {
      stop(
        "The factor '",
        fact,
        "' must have at least two levels for ",
        "ggbetween_pq (between-group comparison requires at least 2 groups)."
      )
    }

    physeq <- taxa_as_columns(physeq)

    if (rarefy_by_sample) {
      if (as(rngseed, "logical")) {
        set.seed(rngseed)
        if (verbose) {
          message(
            "`set.seed(",
            rngseed,
            ")` was used to initialize repeatable random subsampling."
          )
          message(
            "Please record this for your records so others can reproduce."
          )
          message(
            "Try `set.seed(",
            rngseed,
            "); .Random.seed` for the full vector",
            sep = ""
          )
          message("...")
        }
      } else if (verbose) {
        message(
          "You set `rngseed` to FALSE. Make sure you've set & recorded\n",
          " the random seed of your session for reproducibility.\n",
          "See `?set.seed`\n"
        )
        message("...")
      }
      physeq <- clean_pq(rarefy_even_depth_pq(physeq, rngseed = rngseed))
    }

    if (are_modality_even_depth(physeq, fact)$p.value < 0.05) {
      warning(
        paste0(
          "The mean number of sequences per samples vary across modalities of the variable '",
          fact,
          "' You should use rarefy_by_sample = TRUE or try hill_pq() with correction_for_sample_size = TRUE"
        )
      )
    }

    hill_mat <- divent_hill_matrix_pq(
      as.data.frame(physeq@otu_table),
      q = q
    )
    colnames(hill_mat) <- paste0("hill_", q)
    df <- cbind(
      "nb_taxa" = sample_sums(physeq@otu_table),
      physeq@sam_data,
      hill_mat
    )
    fact_sym <- sym(fact)
    res <- lapply(q, function(qi) {
      col_name <- paste0("hill_", qi)
      ggstatsplot::ggbetweenstats(df, !!fact_sym, !!sym(col_name), ...)
    })
    names(res) <- paste0("plot_Hill_", q)

    if (one_plot) {
      requireNamespace("patchwork", quietly = TRUE)
      res <- patchwork::wrap_plots(res)
    }
    return(res)
  }


################################################################################
#' Summarize a \code{\link[phyloseq]{phyloseq-class}} object using a plot.
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-maturing-blue" alt="lifecycle-maturing"></a>
#'
#' Graphical representation of a phyloseq object.
#'
#' @inheritParams clean_pq
#' @param add_info Does the bottom down corner contain
#'   extra informations?
#' @param min_seq_samples (int): Used only when add_info is set
#'   to true to print the number of samples with less sequences than
#'   this number.
#' @param clean_pq (logical): Does the phyloseq
#'   object is cleaned using the [clean_pq()] function?
#' @param text_size (Num, default 1) A size factor to expand or minimize
#'   text size.
#' @param text_size_info (Num, default 1) A size factor to expand or minimize
#'   text size for extra informations.
#' @examples
#'
#' summary_plot_pq(data_fungi_mini)
#' summary_plot_pq(data_fungi_mini, add_info = FALSE) + scale_fill_viridis_d()
#' \donttest{
#' if (requireNamespace("patchwork")) {
#'   (summary_plot_pq(data_fungi, text_size = 0.5, text_size_info = 0.6) +
#'     summary_plot_pq(data_fungi_mini, text_size = 0.5, text_size_info = 0.6)) /
#'     (summary_plot_pq(data_fungi_sp_known, text_size = 0.5, text_size_info = 0.6) +
#'       summary_plot_pq(subset_taxa(data_fungi_sp_known, Phylum == "Ascomycota"),
#'         text_size = 0.5, text_size_info = 0.6
#'       ))
#' }
#' }
#' @return A ggplot2 object
#' @export
summary_plot_pq <- function(
  physeq,
  add_info = TRUE,
  min_seq_samples = 500,
  clean_pq = TRUE,
  text_size = 1,
  text_size_info = 1
) {
  if (clean_pq) {
    physeq <- clean_pq(physeq)
  }
  if (physeq@otu_table@taxa_are_rows) {
    otu_tab <- t(physeq@otu_table)
  } else {
    otu_tab <- physeq@otu_table
  }
  d <- data.frame(
    y1 = c(1, 1, 5.3, 1),
    y2 = c(5, 5, 7.5, 5),
    x1 = c(1, 3.15, 1, 4.3),
    x2 = c(3, 4.15, 3, 5.5),
    names = c("@otu_table", "@tax_table", "@sam_data", "@refseq"),
    nb_values = c(
      paste(
        format(ncol(otu_tab), big.mark = " "),
        "taxa\n",
        format(sum(otu_tab), big.mark = " "),
        "sequences\n",
        format(sum(otu_tab > 0), big.mark = " "),
        "occurrences"
      ),
      paste(ncol(physeq@tax_table), "taxonomic \n levels"),
      paste(
        ncol(physeq@sam_data),
        "variables\n",
        nsamples(physeq),
        "samples"
      ),
      paste(
        "Sequences length:\n",
        ifelse(
          is.null(physeq@refseq),
          "No refseq slot",
          paste(
            round(
              mean(
                Biostrings::width(physeq@refseq)
              ),
              2
            ),
            "+/-",
            round(
              stats::sd(
                Biostrings::width(physeq@refseq)
              ),
              2
            )
          )
        )
      )
    ),
    stringsAsFactors = FALSE
  )

  p <- ggplot() +
    scale_x_continuous(name = "x") +
    scale_y_reverse(name = "y") +
    theme_void() +
    theme(legend.position = "none") +
    scale_fill_manual(values = c("#aa4c26", "#003f5f", "khaki4", "#c8a734")) +
    geom_rect(
      data = d,
      mapping = aes(
        xmin = x1,
        xmax = x2,
        ymin = y1,
        ymax = y2,
        fill = names
      ),
      color = "black",
      alpha = 0.5
    ) +
    geom_text(
      data = d,
      aes(
        x = x1 + (x2 - x1) / 2,
        y = y1 + (y2 - y1) / 1.7,
        label = nb_values
      ),
      size = 4.5 * text_size
    ) +
    geom_text(
      data = d,
      aes(
        x = x1 + (x2 - x1) / 2,
        y = y1 + (y2 - y1) / 5,
        label = names
      ),
      size = 6 * text_size
    ) +
    annotate(
      geom = "text",
      x = 0.65,
      y = 3,
      label = "Taxa",
      size = 6 * text_size,
      color = "#aa4c26",
      angle = 90,
      fontface = 2
    ) +
    annotate(
      geom = "text",
      x = 0.85,
      y = 3,
      label = "(OTUs, ASVs, ...)",
      size = 5 * text_size,
      color = "#aa4c26",
      angle = 90
    ) +
    annotate(
      geom = "text",
      x = 2,
      y = 0.65,
      label = "Samples",
      size = 6 * text_size,
      fontface = 2,
      color = "khaki4"
    )

  if (add_info) {
    supplementary_info <-
      data.frame(
        y1 = 5.3,
        y2 = 7.5,
        x1 = 3.15,
        nb_values = paste0(
          "Min nb seq per sample (",
          stringr::str_trunc(
            names(sort(
              sample_sums(otu_tab)
            ))[1],
            15,
            "right"
          ),
          min(sample_sums(otu_tab)),
          "\n",
          "Nb samples with less than ",
          min_seq_samples,
          " seq : ",
          sum(sample_sums(otu_tab) < min_seq_samples),
          "\n",
          "Min nb seq per taxa: ",
          min(taxa_sums(otu_tab)),
          " (",
          sum(taxa_sums(otu_tab) == min(taxa_sums(otu_tab))),
          " Taxons)",
          "\n",
          "Min seq length: ",
          ifelse(
            is.null(physeq@refseq),
            "No refseq slot",
            min(Biostrings::width(physeq@refseq))
          ),
          "\n",
          "Max nb seq 1 taxa in 1 sample: ",
          max(otu_tab),
          "\n",
          "Max nb of sample for one taxon (",
          names(sort(taxa_sums(otu_tab > 0), decreasing = TRUE))[1],
          "): ",
          max(taxa_sums(otu_tab > 0)),
          "\n",
          "Nb of taxa present in 1 sample only: ",
          sum(taxa_sums(otu_tab > 0) == 1)
        )
      )

    p <- p +
      geom_text(
        data = supplementary_info,
        aes(
          x = x1,
          y = y1 + (y2 - y1) / 2.1,
          label = nb_values
        ),
        size = 3.5 * text_size_info,
        hjust = 0
      )
  }

  return(p)
}
################################################################################

################################################################################
#' rotl wrapper for phyloseq data
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#'   Make a taxonomic tree using the ASV names of a physeq object and the
#'   Open Tree of Life tree.
#'
#' @inheritParams clean_pq
#' @param taxonomic_rank (Character)
#'   The column(s) present in the @tax_table slot of the phyloseq object. Can
#'   be a vector of two columns (e.g. the default c("Genus", "Species")). If only
#'   one column is set it need to be format in this way ("Genus species" for ex.
#'   "Quercus robur") with a space.
#' @param context_name : can bue used to select only a part of the Open Tree
#'   of Life. See `?rotl::tnrs_contexts()` for available values
#' @param discard_genus_alone (logical) If TRUE (default), genus without
#'   information at the species level are discarded.
#' @param pattern_to_remove_tip (character regex string) A regex to remove
#'   unwanted part of tip names. If set to null, tip names are left intact.
#' @param pattern_to_remove_node (character regex string) A regex to remove
#'   unwanted part of node names. If set to null, node names are left intact.
#' @return A plot
#' @export
#' @author Adrien Taudière
#' @details
#' This function is mainly a wrapper of the work of others.
#'   Please make a reference to `rotl` package if you
#'   use this function.
#' @examples
#' \dontrun{
#' if (requireNamespace("rotl")) {
#'   tr <- rotl_pq(data_fungi_mini, pattern_to_remove_tip = NULL)
#'   plot(tr)
#'
#'   tr_Asco <- rotl_pq(data_fungi,
#'     taxonomic_rank = c("Genus", "Species"),
#'     context_name = "Ascomycetes"
#'   )
#'   plot(tr_Asco)
#' }
#' }
rotl_pq <- function(
  physeq,
  taxonomic_rank = c("Genus", "Species"),
  context_name = "All life",
  discard_genus_alone = TRUE,
  pattern_to_remove_tip = c("ott\\d+|_ott\\d+"),
  pattern_to_remove_node = c("_ott.*|mrca*")
) {
  if (sum(!taxonomic_rank %in% colnames(physeq@tax_table)) != 0) {
    stop(
      "The taxonomic_rank parameter do not fit with the @tax_table column of your phyloseq object."
    )
  }

  taxnames <- apply(
    physeq@tax_table[, taxonomic_rank],
    1,
    paste,
    collapse = " "
  )
  taxnames <- taxnames[!is.na(taxnames)]
  if (discard_genus_alone) {
    taxnames <- taxnames[grepl(pattern = " ", taxnames)]
    taxnames <- taxnames[!grepl(pattern = "NA", taxnames)]
  } else {
    taxnames <- taxnames[!grepl(pattern = "NA NA", taxnames)]
    taxnames <- gsub(" NA", "", taxnames)
  }

  taxa_names_rotl <- as.vector(unique(taxnames))

  resolved_names <- httr::with_config(
    httr::config(ssl_verifypeer = FALSE),
    rotl::tnrs_match_names(taxa_names_rotl)
  )
  resolved_names <- resolved_names[resolved_names$flags == "", ]
  clean_taxa_names_rotl <-
    taxa_names_rotl[taxa_names_rotl %in% resolved_names$unique_name]

  resolved_names2 <- httr::with_config(
    httr::config(ssl_verifypeer = FALSE),
    rotl::tnrs_match_names(clean_taxa_names_rotl, context_name = context_name)
  )

  tr <- httr::with_config(
    httr::config(ssl_verifypeer = FALSE),
    rotl::tol_induced_subtree(ott_ids = rotl::ott_id(resolved_names2))
  )

  if (!is.null(pattern_to_remove_tip)) {
    tr$tip.label <- stringr::str_remove(tr$tip.label, pattern_to_remove_tip)
  }

  if (!is.null(pattern_to_remove_node)) {
    tr$node.label <- stringr::str_remove(tr$node.label, pattern_to_remove_node)
  }

  return(tr)
}
################################################################################

################################################################################
# #'  Heat tree from `metacoder` package using `tax_table` slot
# #'  @description
# #'
# #'  <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
# #'  <img src="https://img.shields.io/badge/lifecycle-maturing-blue" alt="lifecycle-maturing"></a>
# #'
# #'  Note that the number of ASV is store under the name `n_obs`
# #'  and the number of sequences under the name `nb_sequences`
# #'
# #'  @inheritParams clean_pq
# #'  @param taxonomic_level (default: NULL) a vector of selected
# #'  taxonomic level using
# #'    their column numbers (e.g. taxonomic_level = 1:7)
# #'  @param ... Arguments passed on to \code{\link[metacoder]{heat_tree}}
# #'
# #'  @return A plot
# #'  @export
# #'  @author Adrien Taudière
# #'
# #'  @examples
# #'  \donttest{
# #'  if (requireNamespace("metacoder")) {
# #'    library("metacoder")
# #'    data("GlobalPatterns", package = "phyloseq")
# #'
# #'    GPsubset <- subset_taxa(
# #'      GlobalPatterns,
# #'      GlobalPatterns@tax_table[, 1] == "Bacteria"
# #'    )
# #'
# #'    GPsubset <- subset_taxa(
# #'      GPsubset,
# #'      rowSums(GPsubset@otu_table) > 5000
# #'    )
# #'
# #'    GPsubset <- subset_taxa(
# #'      GPsubset,
# #'      rowSums(is.na(GPsubset@tax_table)) == 0
# #'    )
# #'
# #'    heat_tree_pq(GPsubset,
# #'      node_size = n_obs,
# #'      node_color = n_obs,
# #'      node_label = taxon_names,
# #'      tree_label = taxon_names,
# #'      node_size_trans = "log10 area"
# #'    )
# #'
# #'    heat_tree_pq(GPsubset,
# #'      node_size = nb_sequences,
# #'      node_color = n_obs,
# #'      node_label = taxon_names,
# #'      tree_label = taxon_names,
# #'      node_size_trans = "log10 area"
# #'    )
# #'  }
# #'  }
# heat_tree_pq <- function(physeq, taxonomic_level = NULL, ...) {
#   requireNamespace("metacoder", quietly = TRUE)
#   if (!is.null(taxonomic_level)) {
#     physeq@tax_table <- physeq@tax_table[, taxonomic_level]
#   }
#
#   data_metacoder <- metacoder::parse_phyloseq(physeq)
#   data_metacoder$data$taxon_counts <-
#     metacoder::calc_taxon_abund(data_metacoder, data = "otu_table")
#   data_metacoder$data$taxon_counts$nb_sequences <-
#     rowSums(data_metacoder$data$taxon_counts[, -1])
#
#   p <- heat_tree(data_metacoder, ...)
#
#   return(p)
# }

################################################################################

################################################################################
#' Visualization of two samples for comparison
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-maturing-blue" alt="lifecycle-maturing"></a>
#'
#' Graphical representation of distribution of taxa across two samples.
#'
#' @inheritParams clean_pq
#' @param fact (default: NULL) Name of the factor in `physeq@sam_data`.
#'   If left to NULL use the `left_name` and `right_name` parameter as modality.
#' @param merge_sample_by (default: NULL) if not `NULL` samples of
#'   physeq are merged using the vector set by `merge_sample_by`. This
#'   merging used the [merge_samples2()]. In the case of
#'   [biplot_pq()] this must be a factor with two levels only.
#' @param rarefy_after_merging Rarefy each sample after merging by the
#'   modalities merge_sample_by
#' @param rngseed (Optional). A single integer value passed to
#'   [phyloseq::rarefy_even_depth()], which is used to fix a seed for
#'   reproducibly random number generation (in this case, reproducibly
#'   random subsampling). If set to FALSE, then no fiddling with the RNG seed
#'   is performed, and it is up to the user to appropriately call set.seed
#'   beforehand to achieve reproducible results. Default is FALSE.
#' @param verbose (logical). If TRUE, print additional information.
#' @param inverse_side Inverse the side (put the right modality in the left side).
#' @param left_name Name fo the left sample.
#' @param left_name_col Color for the left name
#' @param left_fill Fill fo the left sample.
#' @param left_col Color fo the left sample.
#' @param right_name Name fo the right sample.
#' @param right_name_col Color for the right name
#' @param right_fill Fill fo the right sample.
#' @param right_col Color fo the right sample.
#' @param log10trans (logical) Does abundancy is log10 transformed ?
#' @param nudge_y A parameter to control the y position of abundancy values.
#'   If a vector of two values are set. The first value is for the left side.
#'   and the second value for the right one. If one value is set,
#'   this value is used for both side.
#' @param geom_label (default: FALSE, logical) if TRUE use the [ggplot2::geom_label()] function
#'   instead of [ggplot2::geom_text()] to indicate the numbers of sequences.
#' @param text_size size for the number of sequences
#' @param size_names size for the names of the 2 samples
#' @param y_names y position for the names of the 2 samples. If NA (default),
#'   computed using the maximum abundances values.
#' @param ylim_modif vector of two values. Modificator (by a multiplication)
#'   of ylim. If one value is set, this value is used for both limits.
#' @param nb_samples_info (default: TRUE, logical) if TRUE and merge_sample_by is set,
#'   add the number of samples merged for both levels.
#' @param split_by_sample (default: FALSE, logical) if TRUE and
#'   `merge_sample_by` is set, the bars are not merged but stacked by sample,
#'   with borders between segments so that the distribution of sequences
#'   across samples is visible. The border color and width are controlled by
#'   `sample_border_col` and `sample_border_width`.
#' @param sample_border_col (default: "white") Color of the border between
#'   sample segments when `split_by_sample = TRUE`.
#' @param sample_border_width (default: 0.3) Width of the border between
#'   sample segments when `split_by_sample = TRUE`.
#' @param color_rank (default: NULL) Name of a taxonomic rank in `tax_table(physeq)`
#'   to use for coloring bars. When NULL (default), bars are colored by sample
#'   modality using `left_fill` and `right_fill`. When set (e.g. `"Class"`),
#'   each bar is colored according to its taxonomic assignment at that rank
#'   and the `left_fill`/`right_fill` color parameters are ignored.
#' @param taxa_names_rank (default: NULL) Name of a taxonomic rank in
#'   `tax_table(physeq)` to use as labels on the taxa axis instead of
#'   `taxa_names()`. When NULL (default), `taxa_names()` are used. When set
#'   (e.g. `"Genus"`), the genus name is displayed. OTUs sharing the same
#'   label at this rank will appear as a single merged bar.
#' @param plotly_version If TRUE, use [plotly::ggplotly()] to return
#'   a interactive ggplot.
#' @param ... Other arguments for the ggplot function
#' @return A plot
#'
#' @examples
#' data_fungi_2Height <- subset_samples(data_fungi_mini, Height %in% c("Low", "High"))
#' biplot_pq(data_fungi_2Height, "Height", merge_sample_by = "Height")
#' biplot_pq(data_fungi_2Height, "Height",
#'   merge_sample_by = "Height",
#'   split_by_sample = TRUE
#' )
#' biplot_pq(data_fungi_2Height, "Height",
#'   merge_sample_by = "Height",
#'   color_rank = "Order",
#'   taxa_names_rank = "Genus"
#' )
#' @export
#' @author Adrien Taudière
#'
biplot_pq <- function(
  physeq,
  fact = NULL,
  merge_sample_by = NULL,
  rarefy_after_merging = FALSE,
  rngseed = FALSE,
  verbose = TRUE,
  inverse_side = FALSE,
  left_name = NULL,
  left_name_col = "#4B3E1E",
  left_fill = "#4B3E1E",
  left_col = "#4B3E1E",
  right_name = NULL,
  right_name_col = "#1d2949",
  right_fill = "#1d2949",
  right_col = "#1d2949",
  log10trans = TRUE,
  nudge_y = c(0.3, 0.3),
  geom_label = FALSE,
  text_size = 3,
  size_names = 5,
  y_names = NA,
  ylim_modif = c(1, 1),
  nb_samples_info = TRUE,
  split_by_sample = FALSE,
  sample_border_col = "#d4d0acff",
  sample_border_width = 0.3,
  color_rank = NULL,
  taxa_names_rank = NULL,
  plotly_version = FALSE,
  ...
) {
  if (!is.null(merge_sample_by)) {
    if (nb_samples_info) {
      modality_1_nb <- table(physeq@sam_data[, merge_sample_by])[1]
      modality_2_nb <- table(physeq@sam_data[, merge_sample_by])[2]
    }
    if (!split_by_sample) {
      physeq <- merge_samples2(physeq, merge_sample_by)
      physeq <- clean_pq(physeq)
    }
  }

  if (!split_by_sample && nsamples(physeq) != 2) {
    stop(
      "biplot_pq needs only two samples in the
    physeq object or a valid merge_sample_by parameter"
    )
  }

  if (split_by_sample && is.null(merge_sample_by) && is.null(fact)) {
    stop(
      "split_by_sample requires either merge_sample_by or fact to be set"
    )
  }

  if (rarefy_after_merging) {
    if (as(rngseed, "logical")) {
      set.seed(rngseed)
      if (verbose) {
        message(
          "`set.seed(",
          rngseed,
          ")` was used to initialize repeatable random subsampling."
        )
        message("Please record this for your records so others can reproduce.")
        message(
          "Try `set.seed(",
          rngseed,
          "); .Random.seed` for the full vector",
          sep = ""
        )
        message("...")
      }
    } else if (verbose) {
      message(
        "You set `rngseed` to FALSE. Make sure you've set & recorded\n",
        " the random seed of your session for reproducibility.\n",
        "See `?set.seed`\n"
      )
      message("...")
    }
    physeq <- clean_pq(rarefy_even_depth_pq(physeq, rngseed = rngseed))
  }

  if (
    sample_sums(physeq)[1] / sample_sums(physeq)[2] > 2 ||
      sample_sums(physeq)[2] / sample_sums(physeq)[1] > 2
  ) {
    message(
      paste0(
        "The two modalities differ greatly (more than x2) in their number of sequences (",
        sample_sums(physeq)[1],
        " vs ",
        sample_sums(physeq)[2],
        "). You may be interested by the parameter rarefy_after_merging"
      )
    )
  }

  if (is.null(fact)) {
    if (split_by_sample && !is.null(merge_sample_by)) {
      fact <- merge_sample_by
      modality <-
        as.factor(eval(parse(
          text = paste("physeq@sam_data$", fact, sep = "")
        )))
    } else {
      if (is.null(left_name)) {
        left_name <- "A"
      }
      if (is.null(right_name)) {
        right_name <- "B"
      }
      modality <- factor(c(left_name, right_name))
    }
  } else {
    modality <-
      as.factor(eval(parse(
        text = paste("physeq@sam_data$", fact, sep = "")
      )))
  }

  if (inverse_side) {
    modality <- factor(modality, rev(levels(as.factor(modality))))
  }

  if (is.null(left_name)) {
    left_name <- levels(modality)[1]
  }
  if (is.null(right_name)) {
    right_name <- levels(modality)[2]
  }

  if (!is.null(merge_sample_by) && nb_samples_info) {
    left_name <- paste0(left_name, " (", modality_1_nb, " samples)")
    right_name <-
      paste0(right_name, " (", modality_2_nb, " samples)")
  }

  physeq@sam_data$modality <- modality

  mdf <- phyloseq::psmelt(physeq)
  mdf <- mdf[mdf$Abundance > 0, ]

  if (!is.null(taxa_names_rank)) {
    if (!taxa_names_rank %in% colnames(tax_table(physeq))) {
      stop(paste0(
        "'taxa_names_rank' must be a column of tax_table. ",
        "Valid ranks: ",
        paste(colnames(tax_table(physeq)), collapse = ", ")
      ))
    }
    rank_vals <- as.character(mdf[[taxa_names_rank]])
    taxa_label_map <- setNames(rank_vals, mdf$OTU)
    taxa_label_map <- taxa_label_map[!duplicated(names(taxa_label_map))]
  }
  mdf$taxa_label <- mdf$OTU

  if (!is.null(color_rank)) {
    if (!color_rank %in% colnames(tax_table(physeq))) {
      stop(paste0(
        "'color_rank' must be a column of tax_table. ",
        "Valid ranks: ",
        paste(colnames(tax_table(physeq)), collapse = ", ")
      ))
    }
    mdf$fill_var <- as.factor(mdf[[color_rank]])
  } else {
    mdf$fill_var <- mdf$modality
  }
  fill_legend_name <- if (!is.null(color_rank)) color_rank else ""

  if (length(ylim_modif) == 1) {
    ylim_modif <- c(ylim_modif, ylim_modif)
  }

  if (length(y_names) == 1) {
    y_names <- c(y_names, y_names)
  }

  if (length(nudge_y) == 1) {
    nudge_y <- c(nudge_y, nudge_y)
  }
  if (log10trans) {
    mdf$Ab <- log10(mdf$Abundance + 1)
  } else {
    mdf$Ab <- mdf$Abundance
  }

  if (split_by_sample) {
    mdf <- mdf[order(mdf$OTU, mdf$modality, -mdf$Abundance), ]
    mdf$.stack_order <- factor(seq_len(nrow(mdf)), levels = seq_len(nrow(mdf)))
    if (log10trans) {
      mdf <- do.call(
        rbind,
        lapply(
          split(mdf, list(mdf$OTU, mdf$modality), drop = TRUE),
          function(df) {
            total_ab <- sum(df$Abundance)
            total_height <- log10(total_ab + 1)
            df$Ab <- (df$Abundance / total_ab) * total_height
            df
          }
        )
      )
      rownames(mdf) <- NULL
    }
  }

  mdf$Ab[mdf$modality == levels(modality)[1]] <-
    -mdf$Ab[mdf$modality == levels(modality)[1]]
  mdf$Proportion <- paste0(
    round(
      100 *
        mdf$Abundance /
        sum(mdf$Abundance[mdf$modality == levels(modality)[2]]),
      2
    ),
    "%"
  )
  mdf$Proportion[mdf$modality == levels(modality)[1]] <-
    paste0(
      round(
        100 *
          mdf$Abundance[mdf$modality == levels(modality)[1]] /
          sum(mdf$Abundance[mdf$modality == levels(modality)[1]]),
        2
      ),
      "%"
    )

  if (split_by_sample) {
    ab_by_label_mod <- stats::aggregate(
      Ab ~ taxa_label + modality,
      data = mdf,
      FUN = sum
    )
    max_ab <- max(ab_by_label_mod$Ab)
    min_ab <- min(ab_by_label_mod$Ab)
  } else {
    agg_for_lim <- stats::aggregate(
      Ab ~ taxa_label + modality,
      data = mdf,
      FUN = sum
    )
    max_ab <- max(agg_for_lim$Ab)
    min_ab <- min(agg_for_lim$Ab)
  }

  if (split_by_sample) {
    p <- mdf |>
      ggplot(
        aes(
          x = stats::reorder(taxa_label, Abundance),
          y = Ab,
          fill = fill_var,
          group = .stack_order,
          names = taxa_label,
          Ab = Abundance,
          Proportion = Proportion
        ),
        ...
      )
  } else {
    p <- mdf |>
      ggplot(
        aes(
          x = stats::reorder(taxa_label, Abundance),
          y = Ab,
          fill = fill_var,
          names = taxa_label,
          Ab = Abundance,
          Proportion = Proportion
        ),
        ...
      )
  }

  p <- p +
    geom_bar(
      stat = "identity",
      width = 0.6,
      color = if (split_by_sample) sample_border_col else NA,
      linewidth = if (split_by_sample) sample_border_width else 0
    ) +
    annotate(
      "rect",
      xmin = "Samples",
      xmax = "Samples",
      ymin = -max_ab,
      ymax = max_ab
    ) +
    annotate(
      geom = "text",
      label = right_name,
      x = "Samples",
      y = ifelse(is.na(y_names), max_ab / 2, y_names[2]),
      hjust = 0.5,
      vjust = 0.5,
      size = size_names,
      fontface = "bold",
      col = right_name_col
    ) +
    annotate(
      geom = "text",
      label = left_name,
      x = "Samples",
      y = ifelse(is.na(y_names), (min_ab / 2), -y_names[1]),
      hjust = 0.5,
      vjust = 0.5,
      size = size_names,
      fontface = "bold",
      col = left_name_col
    ) +
    geom_hline(aes(yintercept = 0)) +
    scale_x_discrete(
      limits = c(
        names(sort(
          tapply(mdf$Abundance, mdf$taxa_label, sum)
        )),
        "Samples"
      ),
      labels = if (!is.null(taxa_names_rank)) {
        c(taxa_label_map, "Samples" = "")
      } else {
        c("Samples" = "")
      }
    ) +
    ylim(min_ab * 1.1, max_ab * 1.1)

  if (split_by_sample) {
    mdf_totals <- stats::aggregate(
      cbind(Ab, Abundance) ~ taxa_label + modality,
      data = mdf,
      FUN = sum
    )
    mdf_totals$Proportion <- ""
    mdf_totals$Proportion[mdf_totals$modality == levels(modality)[2]] <-
      paste0(
        round(
          100 *
            mdf_totals$Abundance[mdf_totals$modality == levels(modality)[2]] /
            sum(mdf_totals$Abundance[
              mdf_totals$modality == levels(modality)[2]
            ]),
          2
        ),
        "%"
      )
    mdf_totals$Proportion[mdf_totals$modality == levels(modality)[1]] <-
      paste0(
        round(
          100 *
            mdf_totals$Abundance[mdf_totals$modality == levels(modality)[1]] /
            sum(mdf_totals$Abundance[
              mdf_totals$modality == levels(modality)[1]
            ]),
          2
        ),
        "%"
      )
  }

  if (split_by_sample) {
    if (geom_label) {
      p <- p +
        geom_label(
          data = mdf_totals[mdf_totals$Ab > 0, ],
          aes(
            x = taxa_label,
            label = Abundance,
            fill = fill_var,
            alpha = 0.5,
            y = Ab
          ),
          inherit.aes = FALSE,
          color = right_col,
          hjust = -0.1,
          size = text_size
        ) +
        geom_label(
          data = mdf_totals[mdf_totals$Ab < 0, ],
          aes(
            x = taxa_label,
            label = Abundance,
            fill = fill_var,
            alpha = 0.5,
            y = Ab
          ),
          inherit.aes = FALSE,
          color = left_col,
          hjust = 1.1,
          size = text_size
        )
    } else {
      p <- p +
        geom_text(
          data = mdf_totals[mdf_totals$Ab > 0, ],
          aes(x = taxa_label, label = Abundance, y = Ab),
          inherit.aes = FALSE,
          color = right_col,
          hjust = -0.1,
          size = text_size
        ) +
        geom_text(
          data = mdf_totals[mdf_totals$Ab < 0, ],
          aes(x = taxa_label, label = Abundance, y = Ab),
          inherit.aes = FALSE,
          color = left_col,
          hjust = 1.1,
          size = text_size
        )
    }
  } else if (geom_label) {
    p <- p +
      geom_label(
        data = mdf[mdf$Ab > 0, ],
        aes(label = Abundance, fill = fill_var, alpha = 0.5, y = Ab),
        color = right_col,
        hjust = -0.1,
        size = text_size
      ) +
      geom_label(
        data = mdf[mdf$Ab < 0, ],
        aes(label = Abundance, fill = fill_var, alpha = 0.5, y = Ab),
        color = left_col,
        hjust = 1.1,
        size = text_size
      )
  } else {
    p <- p +
      geom_text(
        data = mdf[mdf$Ab > 0, ],
        aes(label = Abundance, y = Ab),
        color = right_col,
        hjust = -0.1,
        size = text_size
      ) +
      geom_text(
        data = mdf[mdf$Ab < 0, ],
        aes(label = Abundance, y = Ab),
        color = left_col,
        hjust = 1.1,
        size = text_size
      )
  }

  p <- p +
    coord_flip() +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5),
      axis.ticks = element_blank()
    ) +
    scale_color_manual(values = c(left_col, right_col), guide = "none") +
    ylim(
      c(
        layer_scales(p)$y$get_limits()[1] * ylim_modif[1],
        layer_scales(p)$y$get_limits()[2] * ylim_modif[2]
      )
    )

  if (is.null(color_rank)) {
    p <- p + scale_fill_manual(values = c(left_fill, right_fill))
  }
  p <- p + labs(fill = fill_legend_name)

  if (plotly_version) {
    if (split_by_sample) {
      p <- mdf |>
        ggplot(
          aes(
            x = stats::reorder(taxa_label, Abundance),
            y = Ab,
            fill = fill_var,
            group = .stack_order,
            names = taxa_label,
            Ab = Abundance,
            Proportion = Proportion,
            Family = Family,
            Genus = Genus,
            Species = Species
          ),
          ...
        )
    } else {
      p <- mdf |>
        ggplot(
          aes(
            x = stats::reorder(taxa_label, Abundance),
            y = Ab,
            fill = fill_var,
            names = taxa_label,
            Ab = Abundance,
            Proportion = Proportion,
            Family = Family,
            Genus = Genus,
            Species = Species
          ),
          ...
        )
    }

    p <- p +
      geom_bar(
        stat = "identity",
        width = 0.6,
        color = if (split_by_sample) sample_border_col else NA,
        linewidth = if (split_by_sample) sample_border_width else 0
      ) +
      annotate(
        "rect",
        xmin = "Samples",
        xmax = "Samples",
        ymin = -max_ab,
        ymax = max_ab
      ) +
      annotate(
        geom = "text",
        label = right_name,
        x = "Samples",
        y = ifelse(is.na(y_names), max_ab / 2, y_names[2]),
        hjust = 0.5,
        vjust = 0.5,
        size = size_names,
        fontface = "bold",
        col = right_name_col
      ) +
      annotate(
        geom = "text",
        label = left_name,
        x = "Samples",
        y = ifelse(is.na(y_names), (min_ab / 2), -y_names[1]),
        hjust = 0.5,
        vjust = 0.5,
        size = size_names,
        fontface = "bold",
        col = left_name_col
      ) +
      geom_hline(aes(yintercept = 0)) +
      scale_x_discrete(
        limits = c(
          names(sort(
            tapply(mdf$Abundance, mdf$taxa_label, sum)
          )),
          "Samples"
        ),
        labels = if (!is.null(taxa_names_rank)) {
          c(taxa_label_map, "Samples" = "")
        } else {
          c("Samples" = "")
        }
      ) +
      ylim(min_ab * 1.1, max_ab * 1.1)

    if (is.null(color_rank)) {
      p <- p + scale_fill_manual(values = c(left_fill, right_fill))
    }
    p <- p + labs(fill = fill_legend_name)

    p <- plotly::ggplotly(
      p,
      tooltip = c("names", "Ab", "Proportion", "Family", "Genus", "Species"),
      height = 1200,
      width = 800
    ) |>
      plotly::layout(
        xaxis = list(autorange = TRUE),
        yaxis = list(autorange = TRUE)
      ) |>
      plotly::config(locale = "fr") |>
      plotly::hide_legend()
  }
  return(p)
}
################################################################################

################################################################################
#' Visualization of a collection of couples of samples for comparison
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' This allow to plot all the possible [biplot_pq()] combination
#' using one factor.
#'
#' @inheritParams clean_pq
#' @param split_by (required if pairs is NULL) the name of the factor to make all combination
#'   of couples of values
#' @param pairs (required if split_by is NULL) the name of the factor in physeq@sam_data` slot
#'   to make plot by pairs of samples. Each level must be present only two times.
#'   Note that if you set pairs, you also must set fact arguments to passed on to [biplot_pq()].
#' @param na_remove (logical, default TRUE) if TRUE remove all the samples
#'   with NA in the `split_by` variable of the `physeq@sam_data` slot
#' @param ... Other parameters passed on to [biplot_pq()]
#'
#' @return a list of ggplot object
#' @export
#'
#' @examples
#' \donttest{
#' data_fungi_abun <- subset_taxa_pq(
#'   data_fungi_mini,
#'   taxa_sums(data_fungi_mini) > 1000
#' )
#' p <- multi_biplot_pq(data_fungi_abun, "Height")
#' lapply(p, print)
#' }
#' @author Adrien Taudière
multi_biplot_pq <- function(
  physeq,
  split_by = NULL,
  pairs = NULL,
  na_remove = TRUE,
  ...
) {
  if (is.null(pairs) && is.null(split_by)) {
    stop("You must set one of split_by or pairs.")
  } else if (!is.null(pairs) && !is.null(split_by)) {
    stop("You must set either split_by or pairs, not both.")
  } else if (
    !is.null(split_by) &&
      is.null(physeq@sam_data[[split_by]])
  ) {
    stop("split_by must be set and must be a variable in physeq@sam_data")
  } else if (!is.null(pairs) && is.null(physeq@sam_data[[pairs]])) {
    stop("pairs must be set and must be a variable in physeq@sam_data")
  }

  if (na_remove && !is.null(split_by)) {
    new_physeq <-
      subset_samples_pq(physeq, !is.na(physeq@sam_data[[split_by]]))
    if (nsamples(physeq) - nsamples(new_physeq) > 0) {
      message(
        paste0(
          nsamples(physeq) - nsamples(new_physeq),
          " were discarded due to NA in variables present in formula."
        )
      )
    }
    physeq <- new_physeq
  }

  if (!is.null(pairs)) {
    p <- vector("list", nlevels(as.factor(physeq@sam_data[[pairs]])))
    names(p) <- levels(as.factor(physeq@sam_data[[pairs]]))
    for (c in levels(as.factor(physeq@sam_data[[pairs]]))) {
      new_physeq <-
        subset_samples_pq(physeq, physeq@sam_data[[pairs]] %in% c)
      p[[c]] <- biplot_pq(new_physeq, ...) + ggtitle(c)
    }
  } else {
    names_split_by <- names(table(physeq@sam_data[[split_by]]))
    couples <- combn(names_split_by, 2)

    p <- vector("list", ncol(couples))
    names(p) <- apply(couples, 2, function(x) {
      paste0(x, collapse = "-")
    })
    for (c in seq_len(ncol(couples))) {
      names_p <- paste0(couples[1, c], "-", couples[2, c])
      new_physeq <-
        subset_samples_pq(
          physeq,
          physeq@sam_data[[split_by]] %in%
            c(couples[1, c], couples[2, c])
        )
      p[[names_p]] <- biplot_pq(
        new_physeq,
        fact = split_by,
        merge_sample_by = split_by,
        ...
      )
    }
  }
  return(p)
}
################################################################################

################################################################################
#' Plot taxonomic distribution in function of a factor with stacked bar in %
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#'   An alternative to `phyloseq::plot_bar()` function.
#'
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor to cluster samples by modalities.
#'   Need to be in \code{physeq@sam_data}.
#' @param merge_sample_by a vector to determine
#'   which samples to merge using the
#'   [merge_samples2()] function.
#'   Need to be in \code{physeq@sam_data}
#' @param type If "nb_seq" (default), the number of sequences is
#'   used in plot. If "nb_taxa", the number of ASV is plotted. If both,
#'   return a list of two plots, one for nbSeq and one for ASV.
#' @param taxa_fill (default: 'Order') Name of the taxonomic rank of interest
#' @param print_values (logical, default TRUE): Do we print some values on plot?
#' @param color_border color for the border
#' @param linewidth The line width of geom_bar
#' @param prop_print_value minimal proportion to print value (default 0.01)
#' @param nb_print_value number of higher values to print
#'    (replace prop_print_value if both are set).
#' @param add_info (logical, default TRUE) Do we add title and subtitle with
#'   information about the total number of sequences and the number of samples
#'   per modality.
#' @param na_remove (logical, default TRUE) if TRUE remove all the samples
#'   with NA in the `split_by` variable of the `physeq@sam_data` slot
#' @param clean_pq (logical)
#'   If set to TRUE, empty samples are discarded after subsetting ASV
#' @return A ggplot2 object
#' @export
#' @author Adrien Taudière
#' @seealso [tax_bar_pq()] and [multitax_bar_pq()]
#' @examples
#' data(data_fungi_sp_known)
#' plot_tax_pq(data_fungi_sp_known,
#'   "Time",
#'   merge_sample_by = "Time",
#'   taxa_fill = "Class"
#' )
#' \donttest{
#' plot_tax_pq(data_fungi_sp_known,
#'   "Height",
#'   merge_sample_by = "Height",
#'   taxa_fill = "Class",
#'   na_remove = TRUE,
#'   color_border = rgb(0, 0, 0, 0)
#' )
#'
#' plot_tax_pq(data_fungi_sp_known,
#'   "Height",
#'   merge_sample_by = "Height",
#'   taxa_fill = "Class",
#'   na_remove = FALSE,
#'   clean_pq = FALSE
#' )
#' }
plot_tax_pq <-
  function(
    physeq,
    fact = NULL,
    merge_sample_by = NULL,
    type = "nb_seq",
    taxa_fill = "Order",
    print_values = TRUE,
    color_border = "lightgrey",
    linewidth = 0.1,
    prop_print_value = 0.01,
    nb_print_value = NULL,
    add_info = TRUE,
    na_remove = TRUE,
    clean_pq = TRUE
  ) {
    if (na_remove) {
      new_physeq <-
        subset_samples_pq(physeq, !is.na(physeq@sam_data[[fact]]))
      if (nsamples(physeq) - nsamples(new_physeq) > 0) {
        message(
          paste0(
            nsamples(physeq) - nsamples(new_physeq),
            " were discarded due to NA in variables present in formula."
          )
        )
      }
      physeq <- new_physeq
    }

    if (clean_pq) {
      physeq <- clean_pq(physeq)
    }

    physeq_old <- physeq

    if (!is.null(merge_sample_by)) {
      physeq <- merge_samples2(physeq, merge_sample_by)
    }

    if (!is.null(nb_print_value)) {
      prop_print_value <-
        taxa_sums(physeq)[nb_print_value] / sum(physeq@otu_table)
    }

    if (type %in% c("nb_seq", "both")) {
      mdf <- psmelt(physeq)
      mdf <- mdf |> mutate(percent = Abundance / sum(Abundance))

      p_seq <-
        ggplot(
          mdf,
          aes(
            x = .data[[fact]],
            y = .data[["Abundance"]],
            fill = .data[[taxa_fill]]
          )
        ) +
        geom_bar(
          aes(fill = .data[[taxa_fill]]),
          stat = "identity",
          position = "fill",
          color = color_border,
          linewidth = linewidth
        ) +
        scale_y_continuous(labels = scales::percent) +
        ylab("Pseudo-abundance (nb of sequences)")

      if (print_values) {
        p_seq <- p_seq +
          geom_text(
            aes(label = round(Abundance)),
            color = ifelse(
              p_seq$data$percent > prop_print_value,
              "black",
              rgb(1, 1, 1, 0)
            ),
            position = position_fill(vjust = 0.5)
          )
      }
    }
    if (type %in% c("nb_taxa", "both")) {
      mdf <-
        psmelt(as_binary_otu_table(physeq))
      mdf <- mdf |> mutate(percent = Abundance / sum(Abundance))

      p_taxa <-
        ggplot(
          mdf,
          aes(
            x = .data[[fact]],
            y = .data[["Abundance"]],
            fill = .data[[taxa_fill]]
          )
        ) +
        geom_bar(
          aes(fill = .data[[taxa_fill]]),
          stat = "identity",
          position = "fill",
          color = color_border,
          linewidth = linewidth
        ) +
        scale_y_continuous(labels = scales::percent) +
        ylab("Nb_ASV")
    }

    if (add_info) {
      if (type %in% c("nb_seq", "both")) {
        p_seq <- p_seq +
          labs(
            title = paste("Total nb of sequences: ", sum(physeq_old@otu_table)),
            subtitle = paste0(
              "Nb of samples: '",
              paste0(
                names(table(physeq_old@sam_data[[fact]])),
                sep = "' : ",
                table(physeq_old@sam_data[[fact]]),
                collapse = " - '"
              )
            )
          )
      }
      if (type %in% c("nb_taxa", "both")) {
        p_taxa <- p_taxa +
          labs(
            title = paste("Total nb of sequences: ", sum(physeq_old@otu_table)),
            subtitle = paste0(
              "Nb of samples: '",
              paste0(
                names(table(physeq_old@sam_data[[fact]])),
                sep = "' : ",
                table(physeq_old@sam_data[[fact]]),
                collapse = " - '"
              )
            )
          )
      }
    }

    if (type == "nb_seq") {
      return(p_seq)
    } else if (type == "nb_taxa") {
      return(p_taxa)
    } else if (type == "both") {
      return(list(p_seq, p_taxa))
    }
  }
################################################################################

################################################################################
#' Plot taxonomic distribution across 3 taxonomic levels and optionally
#' one sample factor
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Note that lvl3 need to be nested in lvl2 which need to be nested
#' in lvl1
#'
#' @inheritParams clean_pq
#' @param lvl1 (required) Name of the first (higher) taxonomic rank of interest
#' @param lvl2 (required) Name of the second (middle) taxonomic rank of interest
#' @param lvl3 (required) Name of the first (lower) taxonomic rank of interest
#' @param fact Name of the factor to cluster samples by modalities.
#'   Need to be in \code{physeq@sam_data}. If not set, the taxonomic
#'   distribution is plot for all samples together.
#' @param nb_seq (logical; default TRUE) If set to FALSE, only the number of ASV
#'   is count. Concretely, physeq otu_table is transformed in a binary
#'   otu_table (each value different from zero is set to one)
#' @param log10trans (logical, default TRUE) If TRUE,
#'   the number of sequences (or ASV if nb_seq = FALSE) is log10
#'   transformed.
#' @return A ggplot2 object
#' @export
#'
#' @author Adrien Taudière
#' @examples
#' \donttest{
#' if (requireNamespace("ggh4x")) {
#'   multitax_bar_pq(data_fungi_sp_known, "Phylum", "Class", "Order", "Time")
#'   multitax_bar_pq(data_fungi_sp_known, "Phylum", "Class", "Order")
#'   multitax_bar_pq(data_fungi_sp_known, "Phylum", "Class", "Order",
#'     nb_seq = FALSE, log10trans = FALSE
#'   )
#' }
#' }
multitax_bar_pq <- function(
  physeq,
  lvl1,
  lvl2,
  lvl3,
  fact = NULL,
  nb_seq = TRUE,
  log10trans = TRUE
) {
  psm_1 <- psmelt(physeq) |>
    filter(Abundance > 0) |>
    filter(!is.na(.data[[lvl1]])) |>
    filter(!is.na(.data[[lvl3]])) |>
    filter(!is.na(.data[[lvl2]]))

  if (is.null(fact)) {
    psm_2 <- psm_1 |>
      group_by(OTU) |>
      summarise(Abundance = sum(Abundance))

    psm <- inner_join(
      psm_2,
      psm_1[, c("OTU", lvl1, lvl2, lvl3)],
      by = join_by("OTU" == "OTU"),
      multiple = "first"
    )

    if (!nb_seq) {
      psm$Abundance <- 1
    }

    data_gg <- tibble(
      "Abundance" = tapply(psm$Abundance, psm[[lvl3]], sum),
      "LVL1" = tapply(psm[[lvl1]], psm[[lvl3]], unique),
      "LVL2" = tapply(psm[[lvl2]], psm[[lvl3]], unique),
      "LVL3" = tapply(psm[[lvl3]], psm[[lvl3]], unique)
    )

    if (log10trans) {
      data_gg$Abundance <- log10(data_gg$Abundance)
    }

    p <- ggplot(
      data_gg,
      aes(
        x = Abundance,
        fill = LVL1,
        y = LVL3
      )
    ) +
      geom_bar(stat = "identity") +
      ggh4x::facet_nested(LVL1 + LVL2 ~ ., scales = "free", space = "free") +
      theme(strip.text.y.right = element_text(angle = 0)) +
      theme(legend.position = "none")
  } else {
    psm_2 <- psm_1 |>
      group_by(OTU, .data[[fact]]) |>
      summarise(Abundance = sum(Abundance)) |>
      filter(Abundance > 0)

    psm <- inner_join(
      psm_2,
      psm_1[, c("OTU", lvl1, lvl2, lvl3)],
      by = join_by("OTU" == "OTU"),
      multiple = "first"
    )

    if (!nb_seq) {
      psm$Abundance <- 1
    }

    data_gg <- tibble(
      "Abundance" = tapply(psm$Abundance, paste(psm[[fact]], psm[[lvl3]]), sum),
      "FACT" = tapply(psm[[fact]], paste(psm[[fact]], psm[[lvl3]]), unique),
      "LVL1" = tapply(psm[[lvl1]], paste(psm[[fact]], psm[[lvl3]]), unique),
      "LVL2" = tapply(psm[[lvl2]], paste(psm[[fact]], psm[[lvl3]]), unique),
      "LVL3" = tapply(psm[[lvl3]], paste(psm[[fact]], psm[[lvl3]]), unique)
    )

    if (log10trans) {
      data_gg$Abundance <- log10(data_gg$Abundance)
    }

    p <- ggplot(
      data_gg,
      aes(
        x = Abundance,
        fill = LVL1,
        y = LVL3
      )
    ) +
      geom_bar(stat = "identity") +
      ggh4x::facet_nested(LVL1 + LVL2 ~ FACT, scales = "free", space = "free") +
      theme(strip.text.y.right = element_text(angle = 0)) +
      theme(legend.position = "none")
  }
  return(p)
}
################################################################################

################################################################################
#' Compute tSNE position of samples from a phyloseq object
#'
#' @inheritParams clean_pq
#' @param method A method to calculate distance using `vegan::vegdist()` function
#' @param dims (Int) Output dimensionality (default: 2)
#' @param theta (Numeric) Speed/accuracy trade-off (increase for less accuracy), set to 0.0 for exact TSNE (default: 0.0 see details in the man page of `Rtsne::Rtsne`).
#' @param perplexity (Numeric) Perplexity parameter (should not be bigger than 3 * perplexity < nrow(X) - 1, see details in the man page of `Rtsne::Rtsne`)
#' @param ... Additional arguments passed on to `Rtsne::Rtsne()`
#'
#' @return A list of element including the matrix Y containing the new representations for the objects.
#'   See ?Rtsne::Rtsne() for more information
#' @export
#'
#' @examplesIf tolower(Sys.info()[["sysname"]]) != "windows"
#' if (requireNamespace("Rtsne")) {
#'   res_tsne <- tsne_pq(data_fungi_mini)
#' }
tsne_pq <-
  function(
    physeq,
    method = "bray",
    dims = 2,
    theta = 0.0,
    perplexity = 30,
    ...
  ) {
    physeq <- taxa_as_rows(physeq)

    res_tsne <-
      Rtsne::Rtsne(
        vegan::vegdist(as(t(physeq@otu_table), "matrix"), method = method),
        dims = dims,
        theta = theta,
        perplexity = perplexity,
        is_distance = TRUE,
        ...
      )

    return(res_tsne)
  }
################################################################################

################################################################################
#' Plot a tsne low dimensional representation of a phyloseq object
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Partially inspired by `phylosmith::tsne_phyloseq()` function developed by Schuyler D. Smith.
#'
#' @inheritParams clean_pq
#' @param method A method to calculate distance using `vegan::vegdist()` function (default: "bray")
#' @param dims (Int) Output dimensionality (default: 2)
#' @param theta (Numeric) Speed/accuracy trade-off (increase for less accuracy), set to 0.0 for exact TSNE (default: 0.0 see details in the man page of `Rtsne::Rtsne`).
#' @param perplexity (Numeric) Perplexity parameter (should not be bigger than 3 * perplexity < nrow(X) - 1, see details in the man page of `Rtsne::Rtsne`)
#' @param fact Name of the column in `physeq@sam_data` used to color points and compute ellipses.
#' @param ellipse_level The level used in stat_ellipse. Set to NULL to discard ellipse (default = 0.95)
#' @param plot_dims A vector of 2 values defining the rank of dimension to plot (default: c(1,2))
#' @param na_remove (logical, default TRUE) Does the samples with NA values in fact are removed? (default: true)
#' @param force_factor (logical, default TRUE) Force the fact column to be a factor.
#' @param ... Additional arguments passed on to `Rtsne::Rtsne()`
#'
#' @return
#' A ggplot object
#'
#' @export
#' @author Adrien Taudière
#'
#' @examplesIf tolower(Sys.info()[["sysname"]]) != "windows"
#' if (requireNamespace("Rtsne")) {
#'   plot_tsne_pq(data_fungi_mini, fact = "Height", perplexity = 15)
#' }
#' \donttest{
#' if (requireNamespace("Rtsne")) {
#'   plot_tsne_pq(data_fungi_mini, fact = "Time") +
#'     geom_label(aes(label = Sample_id, fill = Time))
#'   plot_tsne_pq(data_fungi_mini,
#'     fact = "Time", na_remove = FALSE,
#'     force_factor = FALSE
#'   )
#' }
#' }
#'
plot_tsne_pq <- function(
  physeq,
  method = "bray",
  dims = 2,
  theta = 0.0,
  perplexity = 30,
  fact = NA,
  ellipse_level = 0.95,
  plot_dims = c(1, 2),
  na_remove = TRUE,
  force_factor = TRUE,
  ...
) {
  if (
    !is.factor(physeq@sam_data[[fact]]) &&
      !is.na(fact) &&
      force_factor
  ) {
    physeq@sam_data[[fact]] <- as.factor(physeq@sam_data[[fact]])
  }

  if (na_remove && !is.na(fact)) {
    physeq <- subset_samples_pq(physeq, !is.na(physeq@sam_data[[fact]]))
  }

  if (!is.na(fact) && nlevels(as.factor(physeq@sam_data[[fact]])) < 2) {
    stop(
      "The factor '",
      fact,
      "' must have at least two levels for ",
      "plot_tsne_pq (t-SNE visualization requires at least 2 groups)."
    )
  }

  tsne <- tsne_pq(
    physeq = physeq,
    method = method,
    dims = dims,
    theta = theta,
    perplexity = perplexity,
    ...
  )

  res_tSNE_A <- tsne$Y[, plot_dims[1]] / 100
  res_tSNE_B <- tsne$Y[, plot_dims[2]] / 100

  df <- tibble(
    res_tSNE_A,
    res_tSNE_B,
    as(physeq@sam_data, "data.frame")
  )

  g <-
    ggplot(
      data = df,
      aes(.data[["res_tSNE_A"]], .data[["res_tSNE_B"]], group = .data[[fact]])
    ) +
    xlab(paste0("Dimension ", plot_dims[1], " of tSNE analysis")) +
    ylab(paste0("Dimension ", plot_dims[2], " of tSNE analysis"))

  g <- g +
    geom_point(
      aes(fill = .data[[fact]]),
      shape = 21,
      color = "black",
      size = 3,
      alpha = 1.0
    )

  if (!is.null(ellipse_level) && !is.na(fact)) {
    g <-
      g + stat_ellipse(aes(color = .data[[fact]]), level = ellipse_level)
  }

  return(g)
}
################################################################################

################################################################################
#' Scaling with ranked subsampling (SRS) curve of phyloseq object
#'
#' @description
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' A wraper of [SRS::SRScurve()] function.
#' @inheritParams clean_pq
#' @param clean_pq (logical): Does the phyloseq
#'   object is cleaned using the [clean_pq()] function?
#' @param ... Additional arguments passed on to `SRS::SRScurve()`
#' @return A plot
#' @export
#'
#' @examples
#' if (requireNamespace("SRS")) {
#'   SRS_curve_pq(data_fungi_mini,
#'     max.sample.size = 200,
#'     rarefy.comparison = TRUE, rarefy.repeats = 3
#'   )
#' }
#' \donttest{
#' if (requireNamespace("SRS")) {
#'   SRS_curve_pq(data_fungi_mini, max.sample.size = 500, metric = "shannon")
#' }
#' }
SRS_curve_pq <- function(physeq, clean_pq = FALSE, ...) {
  if (clean_pq) {
    physeq <- clean_pq(physeq)
  }

  physeq <- taxa_as_rows(physeq)

  df <- data.frame(physeq@otu_table)

  SRS::SRScurve(df, ...)
}
################################################################################

################################################################################
#' iNterpolation and EXTrapolation of Hill numbers (with iNEXT)
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Note that this function is quite time-consuming due to high dimensionality in metabarcoding community matrix.
#'
#' @inheritParams clean_pq
#' @param merge_sample_by (default: NULL) if not `NULL` samples of
#'   physeq are merged using the vector set by `merge_sample_by`. This
#'   merging used the [merge_samples2()]. In the case of
#'   [biplot_pq()] this must be a factor with two levels only.
#' @param ... Other arguments for the [iNEXT::iNEXT()] function
#' @return see [iNEXT::iNEXT()] documentation
#' @export
#'
#' @examples
#' \dontrun{
#' if (requireNamespace("iNEXT")) {
#'   data("GlobalPatterns", package = "phyloseq")
#'   GPsubset <- subset_taxa(
#'     GlobalPatterns,
#'     GlobalPatterns@tax_table[, 1] == "Bacteria"
#'   )
#'   GPsubset <- subset_taxa(
#'     GPsubset,
#'     rowSums(GPsubset@otu_table) > 20000
#'   )
#'   GPsubset <- subset_taxa(
#'     GPsubset,
#'     rowSums(is.na(GPsubset@tax_table)) == 0
#'   )
#'   GPsubset@sam_data$human <- GPsubset@sam_data$SampleType %in%
#'     c("Skin", "Feces", "Tong")
#'   res_iNEXT <- iNEXT_pq(
#'     GPsubset,
#'     merge_sample_by = "human",
#'     q = 1,
#'     datatype = "abundance",
#'     nboot = 2
#'   )
#'   iNEXT::ggiNEXT(res_iNEXT)
#'   # iNEXT::ggiNEXT(res_iNEXT, type = 2)
#'   # iNEXT::ggiNEXT(res_iNEXT, type = 3)
#' }
#' }
#' @author Adrien Taudière
#' This function is mainly a wrapper of the work of others.
#'   Please make a reference to `iNEXT::iNEXT()` if you
#'   use this function.
#'
iNEXT_pq <- function(physeq, merge_sample_by = NULL, ...) {
  if (!is.null(merge_sample_by)) {
    physeq <- merge_samples2(physeq, merge_sample_by)
    physeq <- taxa_as_columns(physeq)
  }

  df <- data.frame(t(as.matrix(unclass(physeq@otu_table))))
  res_iNEXT <- iNEXT::iNEXT(df, ...)
  return(res_iNEXT)
}
################################################################################

################################################################################
#' Make upset plot for phyloseq object.
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Alternative to venn plot.
#'
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor to cluster samples by modalities.
#'   Need to be in \code{physeq@sam_data}.
#' @param min_nb_seq minimum number of sequences by OTUs by
#'   samples to take into count this OTUs in this sample. For example,
#'   if min_nb_seq=2,each value of 2 or less in the OTU table
#'   will not count in the venn diagram
#' @param taxa_fill (default NULL) fill the ASV upset using a column in
#'   `tax_table` slot.
#' @param na_remove : if TRUE (the default), NA values in fact are removed
#'   if FALSE, NA values are set to "NA"
#' @param numeric_fonction (default : sum) the function for numeric vector
#'   useful only for complex plot (see examples)
#' @param rarefy_after_merging Rarefy each sample after merging by the
#'   modalities of `fact` parameter
#' @param rngseed (Optional). A single integer value passed to
#'   [phyloseq::rarefy_even_depth()], which is used to fix a seed for
#'   reproducibly random number generation (in this case, reproducibly
#'   random subsampling). If set to FALSE, then no fiddling with the RNG seed
#'   is performed, and it is up to the user to appropriately call set.seed
#'   beforehand to achieve reproducible results. Default is FALSE.
#' @param verbose (logical). If TRUE, print additional information.
#' @param ... Additional arguments passed on to the [ComplexUpset::upset()]
#'
#' @return A \code{\link[ggplot2]{ggplot}}2 plot
#' @export
#' @author Adrien Taudière
#'
#' @seealso [ggvenn_pq()]
#' @examples
#' if (requireNamespace("ComplexUpset") && packageVersion("ggplot2") < "4.0.0") {
#'   upset_pq(data_fungi_mini,
#'     fact = "Height", width_ratio = 0.2,
#'     taxa_fill = "Class"
#'   )
#' }
#' \donttest{
#' if (requireNamespace("ComplexUpset") && packageVersion("ggplot2") < "4.0.0") {
#'   upset_pq(data_fungi_mini, fact = "Height", min_nb_seq = 1000)
#'   upset_pq(data_fungi_mini, fact = "Height", na_remove = FALSE)
#'
#'   upset_pq(data_fungi_mini, fact = "Time", width_ratio = 0.2, rarefy_after_merging = TRUE)
#'
#'   upset_pq(
#'     data_fungi_mini,
#'     fact = "Time",
#'     width_ratio = 0.2,
#'     annotations = list(
#'       "Sequences per ASV \n (log10)" = (
#'         ggplot(mapping = aes(y = log10(Abundance)))
#'         +
#'           geom_jitter(aes(
#'             color =
#'               Abundance
#'           ), na.rm = TRUE)
#'           +
#'           geom_violin(alpha = 0.5, na.rm = TRUE) +
#'           theme(legend.key.size = unit(0.2, "cm")) +
#'           theme(axis.text = element_text(size = 12))
#'       ),
#'       "ASV per phylum" = (
#'         ggplot(mapping = aes(fill = Phylum))
#'         +
#'           geom_bar() +
#'           ylab("ASV per phylum") +
#'           theme(legend.key.size = unit(0.2, "cm")) +
#'           theme(axis.text = element_text(size = 12))
#'       )
#'     )
#'   )
#'
#'   upset_pq(
#'     data_fungi_mini,
#'     fact = "Time",
#'     width_ratio = 0.2,
#'     numeric_fonction = mean,
#'     annotations = list(
#'       "Sequences per ASV \n (log10)" = (
#'         ggplot(mapping = aes(y = log10(Abundance)))
#'         +
#'           geom_jitter(aes(
#'             color =
#'               Abundance
#'           ), na.rm = TRUE)
#'           +
#'           geom_violin(alpha = 0.5, na.rm = TRUE) +
#'           theme(legend.key.size = unit(0.2, "cm")) +
#'           theme(axis.text = element_text(size = 12))
#'       ),
#'       "ASV per phylum" = (
#'         ggplot(mapping = aes(fill = Phylum))
#'         +
#'           geom_bar() +
#'           ylab("ASV per phylum") +
#'           theme(legend.key.size = unit(0.2, "cm")) +
#'           theme(axis.text = element_text(size = 12))
#'       )
#'     )
#'   )
#'
#'   upset_pq(
#'     subset_taxa(data_fungi_mini, Phylum == "Basidiomycota"),
#'     fact = "Time",
#'     width_ratio = 0.2,
#'     base_annotations = list(),
#'     annotations = list(
#'       "Sequences per ASV \n (log10)" = (
#'         ggplot(mapping = aes(y = log10(Abundance)))
#'         +
#'           geom_jitter(aes(
#'             color =
#'               Abundance
#'           ), na.rm = TRUE)
#'           +
#'           geom_violin(alpha = 0.5, na.rm = TRUE) +
#'           theme(legend.key.size = unit(0.2, "cm")) +
#'           theme(axis.text = element_text(size = 12))
#'       ),
#'       "ASV per phylum" = (
#'         ggplot(mapping = aes(fill = Class))
#'         +
#'           geom_bar() +
#'           ylab("ASV per Class") +
#'           theme(legend.key.size = unit(0.2, "cm")) +
#'           theme(axis.text = element_text(size = 12))
#'       )
#'     )
#'   )
#'
#'   data_fungi2 <- data_fungi_mini
#'   data_fungi2@sam_data[["Time_0"]] <- data_fungi2@sam_data$Time == 0
#'   data_fungi2@sam_data[["Height__Time_0"]] <-
#'     paste0(data_fungi2@sam_data[["Height"]], "__", data_fungi2@sam_data[["Time_0"]])
#'   data_fungi2@sam_data[["Height__Time_0"]][grepl("NA", data_fungi2@sam_data[["Height__Time_0"]])] <-
#'     NA
#'   upset_pq(data_fungi2, fact = "Height__Time_0", width_ratio = 0.2, min_size = 2)
#' }
#' }
upset_pq <- function(
  physeq,
  fact,
  taxa_fill = NULL,
  min_nb_seq = 0,
  na_remove = TRUE,
  numeric_fonction = sum,
  rarefy_after_merging = FALSE,
  rngseed = FALSE,
  verbose = TRUE,
  ...
) {
  if (nlevels(as.factor(physeq@sam_data[[fact]])) < 2) {
    stop(
      "The factor '",
      fact,
      "' must have at least two levels for upset_pq ",
      "(UpSet plots require at least 2 sets)."
    )
  }

  if (!is.null(min_nb_seq)) {
    physeq@otu_table[physeq@otu_table < min_nb_seq] <- 0
  }

  if (na_remove) {
    physeq <-
      subset_samples_pq(physeq, !is.na(physeq@sam_data[[fact]]))
  } else {
    physeq@sam_data[[fact]][is.na(physeq@sam_data[[fact]])] <-
      "NA"
  }

  physeq <- merge_samples2(physeq, fact)

  if (rarefy_after_merging) {
    if (as(rngseed, "logical")) {
      set.seed(rngseed)
      if (verbose) {
        message(
          "`set.seed(",
          rngseed,
          ")` was used to initialize repeatable random subsampling."
        )
        message("Please record this for your records so others can reproduce.")
        message(
          "Try `set.seed(",
          rngseed,
          "); .Random.seed` for the full vector",
          sep = ""
        )
        message("...")
      }
    } else if (verbose) {
      message(
        "You set `rngseed` to FALSE. Make sure you've set & recorded\n",
        " the random seed of your session for reproducibility.\n",
        "See `?set.seed`\n"
      )
      message("...")
    }
    physeq <- clean_pq(rarefy_even_depth_pq(physeq, rngseed = rngseed))
  }

  psm <- psmelt(physeq)
  samp_names <- unique(psm$Sample)
  psm <-
    psm |>
    mutate(val = TRUE) |>
    tidyr::pivot_wider(names_from = Sample, values_from = val)
  psm[samp_names][is.na(psm[samp_names])] <- FALSE

  psm <- psm |> filter(Abundance != 0)
  psm[[fact]] <- as.character(psm[[fact]])

  psm2 <- data.frame(lapply(psm, function(col) {
    tapply(col, paste0(psm$OTU), function(vec) {
      diff_fct_diff_class(
        vec,
        numeric_fonction = numeric_fonction,
        na.rm = TRUE
      )
    })
  })) |>
    arrange(desc(Abundance))

  colnames(psm2) <- colnames(psm)

  if (is.null(taxa_fill)) {
    p <-
      ComplexUpset::upset(psm2, intersect = samp_names, ...) + xlab(fact)
  } else {
    p <- ComplexUpset::upset(
      psm2,
      intersect = samp_names,
      base_annotations = list(),
      annotations = list(
        "Taxa" = (ggplot(mapping = aes(fill = .data[[taxa_fill]])) +
          geom_bar() +
          ylab("Taxa per Class") +
          theme(legend.key.size = unit(0.2, "cm")) +
          theme(axis.text = element_text(size = 12)))
      ),
      ...
    ) +
      xlab(fact)
  }

  return(p)
}
################################################################################

################################################################################
#' Test for differences between intersections
#'
#' @description
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' See [upset_pq()] to plot upset. There is a bug with ggplot2 >= 4.0.0. See issue
#'  <https://github.com/krassowski/complex-upset/issues/213> for more details.
#'
#' @inheritParams upset_pq
#' @param var_to_test (default c("OTU")) : a vector of column present in
#'   the tax_table slot from the physeq object
#' @param ... Additional arguments passed on to the [ComplexUpset::upset_test()]
#'
#' @return A \code{\link[ggplot2]{ggplot}}2 plot
#' @export
#' @author Adrien Taudière
#'
#' @seealso [upset_pq()]
#' @examples
#' if (requireNamespace("ComplexUpset")) {
#'   upset_test_pq(data_fungi_mini, "Height",
#'     var_to_test = c("OTU", "Class", "Guild")
#'   )
#'   upset_test_pq(data_fungi_mini, "Time")
#' }
upset_test_pq <-
  function(
    physeq,
    fact,
    var_to_test = "OTU",
    min_nb_seq = 0,
    na_remove = TRUE,
    numeric_fonction = sum,
    ...
  ) {
    if (!is.null(min_nb_seq)) {
      physeq <- subset_taxa_pq(physeq, taxa_sums(physeq) >= min_nb_seq)
    }

    if (na_remove) {
      physeq <-
        subset_samples_pq(physeq, !is.na(physeq@sam_data[[fact]]))
    } else {
      physeq@sam_data[[fact]][is.na(physeq@sam_data[[fact]])] <-
        "NA"
    }

    physeq <- merge_samples2(physeq, fact)

    psm <- psmelt(physeq)
    samp_names <- unique(psm$Sample)
    psm <-
      psm |>
      mutate(val = TRUE) |>
      tidyr::pivot_wider(names_from = Sample, values_from = val)
    psm[samp_names][is.na(psm[samp_names])] <- FALSE

    psm <- psm |> filter(Abundance != 0)
    psm[[fact]] <- as.character(psm[[fact]])

    psm2 <- data.frame(lapply(psm, function(col) {
      tapply(col, paste0(psm$OTU), function(vec) {
        diff_fct_diff_class(
          vec,
          numeric_fonction = numeric_fonction,
          na.rm = TRUE
        )
      })
    })) |>
      arrange(desc(Abundance))

    colnames(psm2) <- colnames(psm)

    res_test <-
      ComplexUpset::upset_test(
        psm2[, c(var_to_test, samp_names)],
        intersect = samp_names,
        ...
      )

    return(res_test)
  }
################################################################################

################################################################################
#' Compute different functions for different class of vector.
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Mainly an internal function useful in "sapply(..., tapply)" methods
#'
#' @param x : a vector
#' @param numeric_fonction : a function for numeric vector. For ex. `sum` or `mean`
#' @param logical_method : A method for logical vector. One of :
#'   - TRUE_if_one (default)
#'   - NA_if_not_all_TRUE
#'   - FALSE_if_not_all_TRUE
#' @param character_method : A method for character vector (and factor). One of :
#'   - unique_or_na (default)
#'   - more_frequent
#'   - more_frequent_without_equality
#' @param ... Additional arguments passed on to the numeric function (ex. na.rm=TRUE)
#' @return a single value
#' @export
#'
#' @examples
#'
#' diff_fct_diff_class(
#'   data_fungi@sam_data$Sample_id,
#'   numeric_fonction = sum,
#'   na.rm = TRUE
#' )
#' diff_fct_diff_class(
#'   data_fungi@sam_data$Time,
#'   numeric_fonction = mean,
#'   na.rm = TRUE
#' )
#' diff_fct_diff_class(
#'   data_fungi@sam_data$Height == "Low",
#'   logical_method = "TRUE_if_one"
#' )
#' diff_fct_diff_class(
#'   data_fungi@sam_data$Height == "Low",
#'   logical_method = "NA_if_not_all_TRUE"
#' )
#' diff_fct_diff_class(
#'   data_fungi@sam_data$Height == "Low",
#'   logical_method = "FALSE_if_not_all_TRUE"
#' )
#' diff_fct_diff_class(
#'   data_fungi@sam_data$Height,
#'   character_method = "unique_or_na"
#' )
#' diff_fct_diff_class(
#'   c("IE", "IE"),
#'   character_method = "unique_or_na"
#' )
#' diff_fct_diff_class(
#'   c("IE", "IE", "TE", "TE"),
#'   character_method = "more_frequent"
#' )
#' diff_fct_diff_class(
#'   c("IE", "IE", "TE", "TE"),
#'   character_method = "more_frequent_without_equality"
#' )
#' @author Adrien Taudière
diff_fct_diff_class <-
  function(
    x,
    numeric_fonction = mean,
    logical_method = "TRUE_if_one",
    character_method = "unique_or_na",
    ...
  ) {
    if (is.character(x) || is.factor(x)) {
      if (length(unique(x)) == 1) {
        return(unique(x))
      } else if (character_method == "unique_or_na") {
        return(NA_character_)
      } else if (character_method == "more_frequent") {
        return(names(sort(table(x), decreasing = TRUE)[1]))
      } else if (character_method == "more_frequent_without_equality") {
        if (
          sort(table(x), decreasing = TRUE)[1] ==
            sort(table(x), decreasing = TRUE)[2]
        ) {
          return(NA_character_)
        } else {
          return(names(sort(table(x), decreasing = TRUE)[1]))
        }
      } else {
        stop(paste0(
          character_method,
          " is not a valid method for character_method params."
        ))
      }
    } else if (is.numeric(x)) {
      return(numeric_fonction(x, ...))
    } else if (is.logical(x)) {
      if (logical_method == "TRUE_if_one") {
        if (sum(x, na.rm = TRUE) > 0) {
          return(TRUE)
        } else {
          return(FALSE)
        }
      }
      if (logical_method == "NA_if_not_all_TRUE") {
        if (sum(x, na.rm = TRUE) > 0 && sum(!x, na.rm = TRUE) == 0) {
          return(TRUE)
        } else if (
          sum(!x, na.rm = TRUE) > 0 &&
            sum(x, na.rm = TRUE) > 0
        ) {
          return(NA)
        } else if (
          sum(!x, na.rm = TRUE) > 0 &&
            sum(x, na.rm = TRUE) == 0
        ) {
          return(FALSE)
        }
      }
      if (logical_method == "FALSE_if_not_all_TRUE") {
        if (sum(x, na.rm = TRUE) > 0 && sum(!x, na.rm = TRUE) == 0) {
          return(TRUE)
        } else {
          return(FALSE)
        }
      } else {
        stop(paste0(
          logical_method,
          " is not a valid method for character_method params."
        ))
      }
    } else {
      stop("At least one column is neither numeric nor character or logical")
    }
  }
################################################################################

################################################################################
#' Plot the distribution of sequences or ASV in one taxonomic levels
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Graphical representation of distribution of taxonomy, optionnaly across a factor.
#'
#' @inheritParams clean_pq
#' @param fact Name of the factor to cluster samples by modalities.
#'   Need to be in \code{physeq@sam_data}.
#' @param taxa (default: 'Order') Name of the taxonomic rank of interest
#' @param percent_bar (default FALSE) If TRUE, the stacked bar fill all
#'   the space between 0 and 1. It just set position = "fill" in the
#'   `ggplot2::geom_bar()` function
#' @param nb_seq (logical; default TRUE) If set to FALSE, only the number of ASV
#'   is count. Concretely, physeq otu_table is transformed in a binary
#'   otu_table (each value different from zero is set to one)
#' @param add_ribbon (logical; default FALSE) If TRUE and `fact` is not
#'   "Sample", add curved ribbons connecting matching taxa between
#'   adjacent bars. Only meaningful when `fact` has more than one level.
#' @param ribbon_alpha (numeric; default 0.3) Transparency of the ribbons.
#' @param label_taxa (logical; default FALSE) If TRUE, replace the legend
#'   with direct labels on the right side of the last bar. Taxa that appear
#'   in the first bar but are absent from the last bar are additionally
#'   labelled on the left side of the first bar. Segments are drawn to
#'   resolve overlapping labels.
#' @param void_theme (logical; default TRUE) If TRUE, use
#'   [ggplot2::theme_void()] when `label_taxa` is TRUE.
#' @param show_values (logical; default FALSE) If TRUE, display
#'   abundance values (or percentages when `percent_bar = TRUE`) inside
#'   bar segments that exceed `minimum_value_to_show`.
#' @param minimum_value_to_show (numeric; default 0) When
#'   `show_values = TRUE`, only segments with a value strictly above
#'   this threshold get a label.
#' @param label_size (numeric; default 3.2) Font size (in ggplot2 mm
#'   units) for taxa labels when `label_taxa = TRUE`.
#' @param value_size (numeric; default 3) Font size (in ggplot2 mm
#'   units) for value labels when `show_values = TRUE`.
#' @param top_label_size (numeric; default 3.2) Font size (in ggplot2 mm
#'   units) for the top group labels when `fact` is not "Sample".
#' @param bar_width (numeric; default NULL set 0.9 if `add_ribbon = FALSE`, 0.5 if
#'   `add_ribbon = TRUE` and `fact != "Sample"`, and 0.6 if fact is only a one-level
#'   factor). Width of the bars. Set to 0 to have no visible bars
#'   and only ribbons.
#' @param bar_internal_color (default NA) Color of bar borders. Use `NA` (default)
#'   to remove borders, which avoids thin white lines in PDF output.
#'   Set to e.g. `"black"` or `"grey30"` for visible borders.
#' @param linewidth_bar_internal (default 0 if `bar_internal_color` is `NA`, otherwise 0.5)
#'  Line width of bar borders.
#' @param show_n_samples (logical; default `TRUE`) If `TRUE`, the number of
#'   samples per group is displayed below each bar as `"(n=X)"`.
#' @param n_sample_text_size (numeric; default `3`) Font size (in ggplot2 mm
#'   units) for the `(n=X)` label displayed below each bar when
#'   `show_n_samples = TRUE`.
#'
#' @return A \code{\link[ggplot2]{ggplot}}2 plot  with bar representing the
#'   number of sequence en each taxonomic groups
#' @export
#'
#' @examples
#'
#' data_fungi_ab <- subset_taxa_pq(
#'   data_fungi_mini,
#'   taxa_sums(data_fungi_mini) > 1000
#' )
#' tax_bar_pq(data_fungi_ab) + theme(legend.position = "none")
#' tax_bar_pq(data_fungi_ab,
#'   taxa = "Class", fact = "Height",
#'   show_n_samples = TRUE
#' )
#' \donttest{
#' tax_bar_pq(data_fungi_ab, taxa = "Class")
#' tax_bar_pq(data_fungi_ab, taxa = "Class", percent_bar = TRUE)
#' tax_bar_pq(data_fungi_ab, taxa = "Class", fact = "Time")
#' tax_bar_pq(data_fungi_ab,
#'   taxa = "Class", fact = "Time",
#'   percent_bar = TRUE, add_ribbon = TRUE
#' )
#' tax_bar_pq(data_fungi_ab,
#'   taxa = "Class", fact = "Time",
#'   percent_bar = TRUE, add_ribbon = TRUE, label_taxa = TRUE
#' )
#' tax_bar_pq(data_fungi_ab,
#'   taxa = "Class", fact = "Time",
#'   show_values = TRUE, minimum_value_to_show = 10000
#' )
#' tax_bar_pq(data_fungi_ab,
#'   fact = "Height", taxa = "Class",
#'   nb_seq = FALSE, percent_bar = TRUE, label_taxa = TRUE,
#'   add_ribbon = TRUE, value_size = 7, ribbon_alpha = .6,
#'   show_values = TRUE, label_size = 4, top_label_size = 6,
#'   minimum_value_to_show = 0.05
#' ) |>
#'   reorder_distinct_colors(alternate_lightness = TRUE)
#'
#' tax_bar_pq(data_fungi_mini,
#'   fact = "Height", taxa = "Order",
#'   nb_seq = TRUE, percent_bar = TRUE, label_taxa = TRUE,
#'   add_ribbon = TRUE, value_size = 5,
#'   ribbon_alpha = .6, show_values = TRUE,
#'   label_size = 4, top_label_size = 8,
#'   minimum_value_to_show = 0.05, bar_width = NULL,
#'   linewidth_bar_internal = 0.1, bar_internal_color = "black"
#' ) |>
#'   reorder_distinct_colors(alternate_lightness = TRUE)
#' }
#' @author Adrien Taudière
#' @seealso [plot_tax_pq()] and [multitax_bar_pq()]
#'
tax_bar_pq <-
  function(
    physeq,
    fact = "Sample",
    taxa = "Order",
    percent_bar = FALSE,
    nb_seq = TRUE,
    add_ribbon = FALSE,
    ribbon_alpha = 0.3,
    label_taxa = FALSE,
    void_theme = TRUE,
    show_values = FALSE,
    minimum_value_to_show = 0,
    label_size = 3.2,
    value_size = 3,
    top_label_size = 3.2,
    bar_width = NULL,
    bar_internal_color = NA,
    linewidth_bar_internal = ifelse(is.na(bar_internal_color), 0, 0.5),
    show_n_samples = TRUE,
    n_sample_text_size = 3
  ) {
    if (!nb_seq) {
      physeq <- as_binary_otu_table(physeq)
    }
    psm <- psmelt(physeq)

    psm[[fact]] <- factor(psm[[fact]])

    if (show_n_samples && fact != "Sample") {
      n_per_group <- psm |>
        dplyr::distinct(.data[[fact]], Sample) |>
        dplyr::count(.data[[fact]])
      n_lookup <- stats::setNames(
        n_per_group$n,
        as.character(n_per_group[[fact]])
      )
    }

    # When nb_seq = FALSE and grouping by a factor other than Sample, each OTU
    # present in multiple samples of the same modality would be counted once per
    # sample (summing binary 1s). The intended value is the number of distinct
    # OTUs in each taxa rank per group: first collapse to presence per OTU per
    # group (max across samples), then sum OTUs per (group x taxa rank).
    if (!nb_seq && fact != "Sample") {
      psm <- psm |>
        dplyr::group_by(.data[[fact]], OTU, .data[[taxa]]) |>
        dplyr::summarise(Abundance = max(Abundance), .groups = "drop") |>
        dplyr::group_by(.data[[fact]], .data[[taxa]]) |>
        dplyr::summarise(Abundance = sum(Abundance), .groups = "drop")
    }

    if (nlevels(psm[[fact]]) < 2) {
      add_ribbon <- FALSE
    }

    if (is.null(bar_width)) {
      bar_width <- if (add_ribbon && fact != "Sample") {
        0.5
      } else if (nlevels(psm[[fact]]) == 1) {
        0.6
      } else {
        0.9
      }
    }

    bar_pos <- if (percent_bar) "fill" else "stack"

    p <- ggplot(psm) +
      geom_bar(
        aes(x = .data[[fact]], fill = .data[[taxa]], y = Abundance),
        stat = "identity",
        position = bar_pos,
        width = bar_width,
        color = bar_internal_color,
        linewidth = linewidth_bar_internal
      )

    if (add_ribbon && fact != "Sample") {
      hw <- (bar_width %||% 0.9) / 2
      sigmoid <- \(x) 1 / (1 + exp(-12 * (x - 0.5)))

      pb <- ggplot_build(p)
      bar_df <- pb$data[[1]]
      taxa_chr <- as.character(psm[[taxa]])
      taxa_chr[is.na(taxa_chr)] <- "NA"
      bar_df$taxa_name <- taxa_chr

      bar_agg <- bar_df |>
        dplyr::group_by(x, taxa_name) |>
        dplyr::summarise(
          ymin = min(ymin),
          ymax = max(ymax),
          .groups = "drop"
        )

      x_vals <- sort(unique(bar_agg$x))

      ribbon_data <- do.call(
        rbind,
        lapply(
          seq_len(length(x_vals) - 1),
          \(i) {
            left <- bar_agg[bar_agg$x == x_vals[i], ]
            right <- bar_agg[bar_agg$x == x_vals[i + 1], ]
            common <- intersect(left$taxa_name, right$taxa_name)
            do.call(
              rbind,
              lapply(common, \(g) {
                l <- left[left$taxa_name == g, ]
                r <- right[right$taxa_name == g, ]
                t_seq <- seq(0, 1, length.out = 50)
                s <- sigmoid(t_seq)
                x_left <- x_vals[i] + hw
                x_right <- x_vals[i + 1] - hw
                x_seq <- x_left + t_seq * (x_right - x_left)
                data.frame(
                  x = c(x_seq, rev(x_seq)),
                  y = c(
                    l$ymin + s * (r$ymin - l$ymin),
                    rev(l$ymax + s * (r$ymax - l$ymax))
                  ),
                  taxa_fill = g,
                  group_id = paste(i, g)
                )
              })
            )
          }
        )
      )

      ribbon_data$taxa_fill[ribbon_data$taxa_fill == "NA"] <- NA

      bar_tops <- bar_agg |>
        dplyr::group_by(x) |>
        dplyr::summarise(ymax = max(ymax), .groups = "drop")
      x_labels <- pb$layout$panel_params[[1]]$x$get_labels()
      x_labels[is.na(x_labels)] <- "NA"
      bar_tops$label <- x_labels[bar_tops$x]

      p <- p +
        geom_polygon(
          data = ribbon_data,
          aes(x = x, y = y, fill = .data[["taxa_fill"]], group = group_id),
          alpha = ribbon_alpha,
          inherit.aes = FALSE
        ) +
        geom_text(
          data = bar_tops,
          aes(x = x, y = ymax, label = label),
          vjust = -0.5,
          inherit.aes = FALSE,
          size = top_label_size
        )
      if (show_n_samples) {
        n_vals_ribbon <- n_lookup[bar_tops$label]
        n_label_ribbon <- data.frame(
          x = bar_tops$x[!is.na(n_vals_ribbon)],
          label = paste0("(n=", n_vals_ribbon[!is.na(n_vals_ribbon)], ")")
        )
        if (nrow(n_label_ribbon) > 0) {
          p <- p +
            geom_text(
              data = n_label_ribbon,
              aes(x = x, y = 0, label = label),
              vjust = 1.5,
              inherit.aes = FALSE,
              size = n_sample_text_size
            )
        }
      }
    }

    if (label_taxa) {
      if (!exists("pb")) {
        pb <- ggplot_build(p)
      }
      if (!exists("bar_df")) {
        bar_df <- pb$data[[1]]
        taxa_chr <- as.character(psm[[taxa]])
        taxa_chr[is.na(taxa_chr)] <- "NA"
        bar_df$taxa_name <- taxa_chr
      }
      if (!exists("bar_agg")) {
        bar_agg <- bar_df |>
          dplyr::group_by(x, taxa_name) |>
          dplyr::summarise(
            ymin = min(ymin),
            ymax = max(ymax),
            .groups = "drop"
          )
      }

      x_max <- max(bar_agg$x)
      last_bar <- bar_agg[bar_agg$x == x_max, ]
      last_bar <- last_bar[last_bar$ymax - last_bar$ymin > 1e-4, ]
      last_bar$ymid <- (last_bar$ymin + last_bar$ymax) / 2

      last_bar <- last_bar[order(last_bar$ymid), ]

      min_gap <- diff(range(c(last_bar$ymin, last_bar$ymax))) /
        (nrow(last_bar) * 1.8)
      label_y <- last_bar$ymid
      for (j in seq_along(label_y)[-1]) {
        if (label_y[j] - label_y[j - 1] < min_gap) {
          label_y[j] <- label_y[j - 1] + min_gap
        }
      }
      last_bar$label_y <- label_y

      hw_bar <- (bar_width %||% 0.9) / 10
      label_df <- data.frame(
        x_bar = x_max + hw_bar,
        x_label = x_max + hw_bar + 0.3,
        y_bar = last_bar$ymid,
        y_label = last_bar$label_y,
        taxa_name = last_bar$taxa_name
      )
      label_df$taxa_fill <- label_df$taxa_name
      label_df$taxa_fill[label_df$taxa_fill == "NA"] <- NA

      needs_segment <- abs(label_df$y_bar - label_df$y_label) > 1e-4

      p <- p +
        geom_text(
          data = label_df,
          aes(
            x = x_label,
            y = y_label,
            label = taxa_name,
            color = taxa_fill
          ),
          hjust = 0,
          size = label_size,
          inherit.aes = FALSE,
          show.legend = FALSE
        ) +
        scale_color_discrete(na.value = "grey50")

      if (void_theme) {
        p <- p + theme_void()
      }
      p <- p +
        theme(
          legend.position = "none",
          plot.margin = margin(5.5, 80, 5.5, 5.5)
        ) +
        coord_cartesian(clip = "off")

      if (any(needs_segment)) {
        seg_df <- label_df[needs_segment, ]
        p <- p +
          geom_segment(
            data = seg_df,
            aes(
              x = x_bar,
              xend = x_label - 0.05,
              y = y_bar,
              yend = y_label,
              color = taxa_fill
            ),
            linewidth = 0.3,
            inherit.aes = FALSE,
            show.legend = FALSE
          )
      }

      # Left-side labels for taxa in first bar but absent from last bar
      x_min <- min(bar_agg$x)
      if (x_min < x_max) {
        first_bar <- bar_agg[bar_agg$x == x_min, ]
        first_bar <- first_bar[first_bar$ymax - first_bar$ymin > 1e-4, ]
        first_bar <- first_bar[
          !first_bar$taxa_name %in% last_bar$taxa_name,
        ]
        if (nrow(first_bar) > 0) {
          first_bar$ymid <- (first_bar$ymin + first_bar$ymax) / 2
          first_bar <- first_bar[order(first_bar$ymid), ]

          min_gap_l <- diff(range(c(first_bar$ymin, first_bar$ymax))) /
            (nrow(first_bar) * 1.8)
          label_y_l <- first_bar$ymid
          for (j in seq_along(label_y_l)[-1]) {
            if (label_y_l[j] - label_y_l[j - 1] < min_gap_l) {
              label_y_l[j] <- label_y_l[j - 1] + min_gap_l
            }
          }
          first_bar$label_y <- label_y_l

          label_df_left <- data.frame(
            x_bar = x_min - hw_bar,
            x_label = x_min - hw_bar - 0.3,
            y_bar = first_bar$ymid,
            y_label = first_bar$label_y,
            taxa_name = first_bar$taxa_name
          )
          label_df_left$taxa_fill <- label_df_left$taxa_name
          label_df_left$taxa_fill[label_df_left$taxa_fill == "NA"] <- NA

          needs_segment_left <-
            abs(label_df_left$y_bar - label_df_left$y_label) > 1e-4

          p <- p +
            geom_text(
              data = label_df_left,
              aes(
                x = x_label,
                y = y_label,
                label = taxa_name,
                color = taxa_fill
              ),
              hjust = 1,
              size = label_size,
              inherit.aes = FALSE,
              show.legend = FALSE
            ) +
            theme(plot.margin = margin(5.5, 80, 5.5, 80))

          if (any(needs_segment_left)) {
            seg_df_left <- label_df_left[needs_segment_left, ]
            p <- p +
              geom_segment(
                data = seg_df_left,
                aes(
                  x = x_bar,
                  xend = x_label + 0.05,
                  y = y_bar,
                  yend = y_label,
                  color = taxa_fill
                ),
                linewidth = 0.3,
                inherit.aes = FALSE,
                show.legend = FALSE
              )
          }
        }

        # Warn about taxa only in intermediate bars (neither first nor last)
        all_taxa_with_data <- unique(
          bar_agg$taxa_name[bar_agg$ymax - bar_agg$ymin > 1e-4]
        )
        first_taxa <- unique(
          bar_agg$taxa_name[
            bar_agg$x == x_min & bar_agg$ymax - bar_agg$ymin > 1e-4
          ]
        )
        unlabeled <- setdiff(
          all_taxa_with_data,
          union(last_bar$taxa_name, first_taxa)
        )
        if (length(unlabeled) > 0) {
          warning(
            length(unlabeled),
            " taxon/taxa only appear in intermediate levels and will not ",
            "be labelled: ",
            paste(unlabeled, collapse = ", "),
            ". Consider using label_taxa = FALSE.",
            call. = FALSE
          )
        }
      }
    } else {
      if (void_theme) {
        p <- p + theme_void()
      }
    }

    if (show_values) {
      pb_val <- ggplot_build(p)
      val_df <- pb_val$data[[1]]

      # Aggregate per bar x-position and fill group
      val_agg <- val_df |>
        dplyr::group_by(x, group) |>
        dplyr::summarise(
          ymin = min(ymin),
          ymax = max(ymax),
          fill = fill[1],
          .groups = "drop"
        )
      val_agg$seg_value <- val_agg$ymax - val_agg$ymin
      val_agg$ymid <- (val_agg$ymin + val_agg$ymax) / 2

      if (percent_bar) {
        val_agg$label_text <- paste0(
          round(val_agg$seg_value * 100, 1),
          "%"
        )
      } else {
        val_agg$label_text <- as.character(round(val_agg$seg_value))
      }

      val_agg <- val_agg[val_agg$seg_value >= minimum_value_to_show, ]

      if (nrow(val_agg) > 0) {
        p <- p +
          geom_text(
            data = val_agg,
            aes(x = x, y = ymid, label = label_text),
            inherit.aes = FALSE,
            size = value_size,
            color = "white"
          )
      }
    }

    if (fact != "Sample" && !add_ribbon) {
      pb_n <- ggplot_build(p)
      bar_tops_n <- pb_n$data[[1]] |>
        dplyr::group_by(x) |>
        dplyr::summarise(ymax = max(ymax), .groups = "drop")
      x_labs_n <- pb_n$layout$panel_params[[1]]$x$get_labels()
      x_labs_n[is.na(x_labs_n)] <- "NA"
      bar_tops_n$group_label <- x_labs_n[bar_tops_n$x]
      bar_tops_n$label <- bar_tops_n$group_label
      p <- p +
        geom_text(
          data = bar_tops_n,
          aes(x = x, y = ymax, label = label),
          vjust = -0.5,
          inherit.aes = FALSE,
          size = top_label_size
        )
      if (show_n_samples) {
        n_vals_n <- n_lookup[bar_tops_n$group_label]
        n_label_n <- data.frame(
          x = bar_tops_n$x[!is.na(n_vals_n)],
          label = paste0("(n=", n_vals_n[!is.na(n_vals_n)], ")")
        )
        if (nrow(n_label_n) > 0) {
          p <- p +
            geom_text(
              data = n_label_n,
              aes(x = x, y = 0, label = label),
              vjust = 1.5,
              inherit.aes = FALSE,
              size = n_sample_text_size
            )
        }
      }
    }

    p
  }
################################################################################

################################################################################
################################################################################
#' Ridge plot of a phyloseq object
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Graphical representation of distribution of taxa across a factor using ridges.
#'
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor in `physeq@sam_data` used to plot
#'    different lines
#' @param nb_seq (logical; default TRUE) If set to FALSE, only the number of ASV
#'   is count. Concretely, physeq `otu_table` is transformed in a binary
#'   `otu_table` (each value different from zero is set to one)
#' @param log10trans (logical, default TRUE) If TRUE,
#'   the number of sequences (or ASV if nb_seq = FALSE) is log10
#'   transformed.
#' @param tax_level The taxonomic level to fill ridges
#' @param type Either "density" (the default) or "ecdf" to plot a
#'   plot a cumulative version using [ggplot2::stat_ecdf()]
#' @param ... Other params passed on to [ggridges::geom_density_ridges()]
#'
#' @return A \code{\link[ggplot2]{ggplot}}2 plot  with bar representing the number of sequence en each
#'   taxonomic groups
#' @export
#' @author Adrien Taudière
#' @examples
#' if (requireNamespace("ggridges")) {
#'   ridges_pq(data_fungi_mini, "Time", alpha = 0.5, log10trans = FALSE) + xlim(c(0, 1000))
#' }
#' \donttest{
#' if (requireNamespace("ggridges")) {
#'   ridges_pq(data_fungi_mini, "Time", alpha = 0.5, scale = 0.9)
#'   ridges_pq(data_fungi_mini, "Time", alpha = 0.5, scale = 0.9, type = "ecdf")
#'   ridges_pq(data_fungi_mini, "Sample_names", log10trans = TRUE) + facet_wrap("~Height")
#'
#'   ridges_pq(data_fungi_mini,
#'     "Time",
#'     jittered_points = TRUE,
#'     position = ggridges::position_points_jitter(width = 0.05, height = 0),
#'     point_shape = "|", point_size = 3, point_alpha = 1, alpha = 0.7,
#'     scale = 0.8
#'   )
#' }
#' }
ridges_pq <- function(
  physeq,
  fact,
  nb_seq = TRUE,
  log10trans = TRUE,
  tax_level = "Class",
  type = "density",
  ...
) {
  psm <- psmelt(physeq)
  psm <- psm |> dplyr::filter(Abundance > 0)

  if (log10trans) {
    psm$Abundance <- log10(psm$Abundance)
  }
  if (nb_seq) {
    p <- ggplot(
      psm,
      aes(
        y = factor(.data[[fact]]),
        x = Abundance,
        fill = .data[[tax_level]],
        color = .data[[tax_level]]
      )
    )
  } else {
    psm_asv <-
      psm |>
      group_by(.data[[fact]], OTU, .data[[tax_level]]) |>
      summarise("count" = n())

    p <- ggplot(
      psm_asv,
      aes(
        y = factor(.data[[fact]]),
        x = count,
        fill = .data[[tax_level]],
        color = .data[[tax_level]]
      )
    )
  }

  if (type == "density") {
    p <- p +
      ggridges::geom_density_ridges(aes(), ...) +
      xlim(c(0, NA))
  } else if (type == "ecdf") {
    p <- p +
      stat_ecdf(aes(y = NULL)) +
      facet_wrap(fact, ncol = 2) +
      theme_minimal() +
      ylab("Probability")
  }
  return(p)
}
################################################################################

################################################################################
#' Ridges plot of sample distribution across taxa
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Graphical representation of distribution of samples across taxa using ridges.
#' This is the sample-centric counterpart of [ridges_pq()]: each ridge
#' represents a taxon (at `tax_level`) and the x-axis shows the abundance
#' distribution across samples, optionally colored by a sample factor.
#'
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor in `physeq@sam_data` used to color
#'   the ridges
#' @param nb_seq (logical; default TRUE) If set to FALSE, only the number of
#'   samples is counted. Concretely, physeq `otu_table` is transformed in a
#'   binary `otu_table` (each value different from zero is set to one)
#' @param log10trans (logical, default TRUE) If TRUE,
#'   the abundance is log10 transformed.
#' @param tax_level The taxonomic level used for grouping taxa on the y-axis
#' @param type Either "density" (the default) or "ecdf" to plot a
#'   cumulative version using [ggplot2::stat_ecdf()]
#' @param ... Other params passed on to [ggridges::geom_density_ridges()]
#'
#' @return A \code{\link[ggplot2]{ggplot}}2 plot with ridges representing the
#'   distribution of samples for each taxon
#' @export
#' @author Adrien Taudière
#' @examples
#' if (requireNamespace("ggridges")) {
#'   ridges_sam_pq(data_fungi_mini, "Height",
#'     alpha = 0.5,
#'     log10trans = FALSE, tax_level = "Genus"
#'   ) +
#'     xlim(c(0, 1000))
#' }
#' \donttest{
#' if (requireNamespace("ggridges")) {
#'   ridges_sam_pq(data_fungi_mini, "Height", alpha = 0.5, scale = 0.9)
#'   ridges_sam_pq(data_fungi_mini, "Height",
#'     alpha = 0.5, scale = 0.9,
#'     type = "ecdf"
#'   )
#' }
#' }
ridges_sam_pq <- function(
  physeq,
  fact,
  nb_seq = TRUE,
  log10trans = TRUE,
  tax_level = "Class",
  type = "density",
  ...
) {
  psm <- psmelt(physeq)
  psm <- psm |> dplyr::filter(Abundance > 0)

  if (log10trans) {
    psm$Abundance <- log10(psm$Abundance)
  }
  if (nb_seq) {
    p <- ggplot(
      psm,
      aes(
        y = factor(.data[[tax_level]]),
        x = Abundance,
        fill = .data[[fact]],
        color = .data[[fact]]
      )
    )
  } else {
    psm_sam <-
      psm |>
      group_by(.data[[tax_level]], Sample, .data[[fact]]) |>
      summarise("count" = n())

    p <- ggplot(
      psm_sam,
      aes(
        y = factor(.data[[tax_level]]),
        x = count,
        fill = .data[[fact]],
        color = .data[[fact]]
      )
    )
  }

  if (type == "density") {
    p <- p +
      ggridges::geom_density_ridges(aes(), ...) +
      xlim(c(0, NA))
  } else if (type == "ecdf") {
    p <- p +
      stat_ecdf(aes(y = NULL)) +
      facet_wrap(tax_level, ncol = 2) +
      theme_minimal() +
      ylab("Probability")
  }
  return(p)
}
################################################################################

################################################################################
#' Plot treemap of 2 taxonomic levels
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Note that lvl2need to be nested in lvl1
#'
#' @inheritParams clean_pq
#' @param lvl1 (required) Name of the first (higher) taxonomic rank of interest
#' @param lvl2 (required) Name of the second (lower) taxonomic rank of interest
#' @param nb_seq (logical; default TRUE) If set to FALSE, only the number of ASV
#'   is count. Concretely, physeq otu_table is transformed in a binary
#'   otu_table (each value different from zero is set to one)
#' @param log10trans (logical, default TRUE) If TRUE,
#'   the number of sequences (or ASV if nb_seq = FALSE) is
#'   log10(x + 1) transformed. The +1 ensures that taxa with a
#'   count of 1 still have a visible tile area.
#' @param plot_legend (logical, default FALSE) If TRUE, plot che
#'   legend of color for lvl 1
#' @param show_count (logical, default FALSE) If TRUE, appends the raw
#'   count in parentheses after each `lvl2` label, e.g. `"Agaricus (42)"`.
#' @param facet_by (character, default NULL) Name of a column in
#'   `sample_data(physeq)` to facet by. Each level produces its own
#'   treemap panel via [ggplot2::facet_wrap()].
#' @param growing_text (logical, default TRUE) If FALSE, all tile labels are
#'   drawn at the same font size (disables per-tile text growing), which
#'   corresponds to the smallest size that would otherwise be computed.
#' @param text_size (numeric, default 15) Base font size for tile labels.
#'   Mostly useful when `growing_text = FALSE`, as it sets the size of all
#'   labels.
#' @param show_na (logical, default TRUE) If TRUE, taxa with NA values for
#'   `lvl1` or `lvl2` are kept and displayed as a grey "NA" area. If FALSE,
#'   they are removed (previous default behavior).
#' @param na_label (character, default "NA") Label used to replace NA values
#'   in `lvl1` and `lvl2` when `show_na = TRUE`.
#' @param min_text_size (numeric, default 0) Minimum font size in points
#'   for tile labels. Labels that would be smaller than this are hidden.
#'   Set to 0 to always show all labels.
#' @param ... Additional arguments passed on to
#'   [treemapify::geom_treemap()] function.
#'
#' @return A ggplot2 object
#' @export
#'
#' @author Adrien Taudière
#' @examples
#' data(data_fungi_sp_known)
#' if (requireNamespace("treemapify")) {
#'   treemap_pq(
#'     clean_pq(subset_taxa(
#'       data_fungi_sp_known,
#'       Phylum == "Basidiomycota"
#'     )),
#'     "Order", "Class",
#'     plot_legend = TRUE
#'   )
#' }
#' \donttest{
#' if (requireNamespace("treemapify")) {
#'   treemap_pq(
#'     clean_pq(subset_taxa(
#'       data_fungi_sp_known,
#'       Phylum == "Basidiomycota"
#'     )),
#'     "Order", "Class",
#'     log10trans = FALSE
#'   )
#'   treemap_pq(
#'     clean_pq(subset_taxa(
#'       data_fungi_sp_known,
#'       Phylum == "Basidiomycota"
#'     )),
#'     "Order", "Class",
#'     nb_seq = FALSE, log10trans = FALSE
#'   )
#'   treemap_pq(
#'     clean_pq(subset_taxa(
#'       data_fungi_sp_known,
#'       Phylum == "Basidiomycota"
#'     )),
#'     "Order", "Class",
#'     show_count = TRUE, log10trans = FALSE
#'   )
#' }
#' }
treemap_pq <- function(
  physeq,
  lvl1,
  lvl2,
  nb_seq = TRUE,
  log10trans = TRUE,
  plot_legend = FALSE,
  show_count = FALSE,
  facet_by = NULL,
  growing_text = TRUE,
  text_size = 15,
  show_na = TRUE,
  na_label = "NA",
  min_text_size = 0,
  ...
) {
  if (!nb_seq) {
    physeq <- as_binary_otu_table(physeq)
  }

  if (!is.null(facet_by)) {
    sam <- as.data.frame(sample_data(physeq))
    if (!facet_by %in% colnames(sam)) {
      stop(
        "Column '",
        facet_by,
        "' not found in sample_data. ",
        "Available: ",
        paste(colnames(sam), collapse = ", ")
      )
    }
  }

  psm <- psmelt(physeq)

  if (show_na) {
    psm[[lvl1]] <- ifelse(is.na(psm[[lvl1]]), na_label, psm[[lvl1]])
    psm[[lvl2]] <- ifelse(is.na(psm[[lvl2]]), na_label, psm[[lvl2]])
  } else {
    psm <- psm |>
      filter(!is.na(.data[[lvl2]])) |>
      filter(!is.na(.data[[lvl1]]))
  }

  if (!is.null(facet_by)) {
    psm2 <- psm |>
      group_by(.data[[lvl2]], .data[[facet_by]]) |>
      reframe(Abundance = sum(Abundance), LVL1 = unique(.data[[lvl1]]))
  } else {
    psm2 <- psm |>
      group_by(.data[[lvl2]]) |>
      reframe(Abundance = sum(Abundance), LVL1 = unique(.data[[lvl1]]))
  }

  psm2$raw_count <- psm2$Abundance

  if (log10trans) {
    psm2$Abundance <- log10(psm2$Abundance + 1)
  }

  if (show_count) {
    psm2$label <- paste0(psm2[[lvl2]], "\n(", psm2$raw_count, ")")
  } else {
    psm2$label <- psm2[[lvl2]]
  }

  p <-
    ggplot(
      psm2,
      aes(
        area = Abundance,
        fill = LVL1,
        label = label,
        subgroup = LVL1
      )
    ) +
    treemapify::geom_treemap(...) +
    treemapify::geom_treemap_subgroup_border(colour = "white", size = 4) +
    treemapify::geom_treemap_text(
      colour = "white",
      place = "centre",
      size = text_size,
      grow = growing_text,
      min.size = min_text_size
    )

  if (show_na && na_label %in% psm2$LVL1) {
    lvl1_values <- unique(psm2$LVL1)
    non_na <- setdiff(lvl1_values, na_label)
    n_colors <- length(non_na)
    default_colors <- scales::hue_pal()(n_colors)
    fill_colors <- stats::setNames(default_colors, non_na)
    fill_colors[[na_label]] <- "grey70"
    p <- p + scale_fill_manual(values = fill_colors)
  }

  if (!is.null(facet_by)) {
    p <- p + facet_wrap(vars(.data[[facet_by]]))
  }

  if (!plot_legend) {
    p <- p + theme(legend.position = "none")
  }

  if (nb_seq) {
    if (log10trans) {
      p <-
        p +
        ggtitle(paste0(
          "Nb of sequences (log10 transformed) by ",
          lvl1,
          " and ",
          lvl2
        ))
    } else {
      p <- p + ggtitle(paste0("Nb of sequences by ", lvl1, " and ", lvl2))
    }
  } else {
    if (log10trans) {
      p <- p +
        ggtitle(paste0(
          "Nb of ASV (log10 transformed) by ",
          lvl1,
          " and ",
          lvl2
        ))
    } else {
      p <- p + ggtitle(paste0("Nb of ASV by ", lvl1, " and ", lvl2))
    }
  }

  return(p)
}
################################################################################

################################################################################
#' Plot the partition the variation of a phyloseq object
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Graphical representation of the partition of variation obtain with [var_par_pq()].
#' @param res_varpart (required) the result of the functions [var_par_pq()]
#'   or [var_par_rarperm_pq()]
#' @param cutoff The values below cutoff will not be displayed.
#' @param digits The number of significant digits.
#' @param digits_quantile The number of significant digits for quantile.
#' @param fill_bg Fill colours of ellipses.
#' @param show_quantiles Do quantiles are printed ?
#' @param filter_quantile_zero Do we filter out value with quantile encompassing
#'   the zero value?
#' @param show_dbrda_signif Do dbrda significance for each component is printed
#'   using *?
#' @param show_dbrda_signif_pval (float, `[0:1]`) The value under which the
#'  dbrda is considered significant.
#' @param alpha (int, `[0:255]`) Transparency of the fill colour.
#' @param id.size A numerical value giving the character expansion factor for the names of circles or ellipses.
#' @param min_prop_pval_signif_dbrda (float, `[0:1]`) Only used if using the
#'   result of [var_par_rarperm_pq()] function. The * for dbrda_signif is only add if
#'   at least `min_prop_pval_signif_dbrda` of permutations show significance.
#'
#' @return A plot
#' @export
#' @author Adrien Taudière
#' @seealso [var_par_rarperm_pq()], [var_par_pq()]
#' @examples
#' \donttest{
#' if (requireNamespace("vegan")) {
#'   data_fungi_woNA <- subset_samples(
#'     data_fungi_mini,
#'     !is.na(Time) & !is.na(Height)
#'   )
#'   res_var0 <- var_par_pq(data_fungi_woNA,
#'     list_component = list(
#'       "Time" = c("Time"),
#'       "Size" = c("Height", "Diameter")
#'     )
#'   )
#'   plot_var_part_pq(res_var0)
#' }
#' }
#' \dontrun{
#' if (requireNamespace("vegan")) {
#'   res_var_2 <- var_par_rarperm_pq(
#'     data_fungi_woNA,
#'     list_component = list(
#'       "Time" = c("Time"),
#'       "Size" = c("Height", "Diameter")
#'     ),
#'     nperm = 2,
#'     dbrda_computation = TRUE
#'   )
#'   plot_var_part_pq(res_var0, digits_quantile = 2, show_dbrda_signif = TRUE)
#'   plot_var_part_pq(
#'     res_var_2,
#'     digits = 5,
#'     digits_quantile = 2,
#'     cutoff = 0,
#'     show_quantiles = TRUE
#'   )
#' }
#' }
#' @importFrom stats anova as.formula quantile
#' @details
#' This function is mainly a wrapper of the work of others.
#'   Please make a reference to `vegan::varpart()` if you
#'   use this function.
plot_var_part_pq <-
  function(
    res_varpart,
    cutoff = 0,
    digits = 1,
    digits_quantile = 2,
    fill_bg = c("seagreen3", "mediumpurple", "blue", "orange"),
    show_quantiles = FALSE,
    filter_quantile_zero = TRUE,
    show_dbrda_signif = FALSE,
    show_dbrda_signif_pval = 0.05,
    alpha = 63,
    id.size = 1.2,
    min_prop_pval_signif_dbrda = 0.95
  ) {
    if (show_dbrda_signif_pval > 1 || show_dbrda_signif_pval < 0) {
      stop("show_dbrda_signif_pval value must be within the range [0-1]")
    }
    if (
      min_prop_pval_signif_dbrda > 1 ||
        min_prop_pval_signif_dbrda < 0
    ) {
      stop("show_dbrda_signif_pval value must be within the range [0-1]")
    }
    x <- res_varpart$part
    vals <- x$indfract$Adj.R.square
    is.na(vals) <- vals < cutoff
    vals <- round(vals, digits + 1)
    labs_text <- format(vals, digits = digits, nsmall = digits + 1)
    labs_text <- gsub("NA", "", labs_text)
    if (show_quantiles) {
      labs_text <- paste0(
        labs_text,
        "\n (",
        round(x$indfract$Adj.R.squared_quantil_min, digits_quantile + 1),
        "...",
        round(x$indfract$Adj.R.squared_quantil_max, digits_quantile + 1),
        ")"
      )
      labs_text[is.na(vals)] <- ""
    }

    if (filter_quantile_zero) {
      labs_text[x$indfract$Adj.R.squared_quantil_min < 0] <- ""
    }

    if (show_dbrda_signif) {
      if (is.null(res_varpart$dbrda_result_prop_pval_signif)) {
        cond <-
          seq_along(res_varpart$dbrda_result)[sapply(
            res_varpart$dbrda_result,
            function(x) {
              x$`Pr(>F)`[[1]] < show_dbrda_signif_pval
            }
          )]
        res_varpart$Xnames[cond] <-
          paste0(res_varpart$Xnames[cond], "*")
      } else {
        cond <-
          seq_along(res_varpart$dbrda_result)[
            res_varpart$dbrda_result_prop_pval_signif >=
              min_prop_pval_signif_dbrda
          ]
        res_varpart$Xnames[cond] <-
          paste0(res_varpart$Xnames[cond], "*")
      }
    }

    vegan::showvarparts(
      x$nsets,
      labs_text,
      bg = fill_bg,
      alpha = alpha,
      id.size = id.size,
      Xnames = res_varpart$Xnames
    )
    if (anyNA(vals)) {
      graphics::mtext(paste("Values <", cutoff, " not shown", sep = ""), 1)
    }
    if (
      sum(x$indfract$Adj.R.squared_quantil_min) > 0 &&
        filter_quantile_zero
    ) {
      graphics::mtext(
        paste("Values with min quantile <0 not shown", sep = ""),
        side = 1,
        line = 1
      )
    }
    if (show_dbrda_signif) {
      if (is.null(res_varpart$dbrda_result_prop_pval_signif)) {
        graphics::mtext(
          paste(
            "* indicate significant anova of dbRDA for each component at p=",
            show_dbrda_signif_pval,
            sep = ""
          )
        )
      } else {
        graphics::mtext(
          paste(
            "* indicate significant anova of dbRDA, for each component, in at least ",
            round(min_prop_pval_signif_dbrda * 100, 2),
            "% of rarefaction permutations",
            sep = ""
          )
        )
      }
    }
    return(invisible())
  }
################################################################################

################################################################################
#' Scatterplot with marginal distributions and statistical results against
#' Hill diversity of phyloseq object
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Basically a wrapper of function [ggstatsplot::ggscatterstats()] for
#' object of class phyloseq and Hill number.
#'
#' @inheritParams clean_pq
#' @param num_modality (required) Name of the numeric column in
#'   `physeq@sam_data` to plot and test against hill number
#' @param q (a vector of integer) The list of q values to compute
#'   the hill number H^q. If Null, no hill number are computed. Default value
#'   compute the Hill number 0 (Species richness), the Hill number 1
#'   (exponential of Shannon Index) and the Hill number 2 (inverse of Simpson
#'   Index). Hill numbers are more appropriate in DNA metabarcoding studies
#'   when `q > 0` (Alberdi & Gilbert, 2019; Calderón-Sanou et al., 2019).
#' @param rarefy_by_sample (logical, default FALSE) If TRUE, rarefy
#'   samples using [phyloseq::rarefy_even_depth()] function.
#' @param rngseed (Optional). A single integer value passed to
#'   [phyloseq::rarefy_even_depth()], which is used to fix a seed for
#'   reproducibly random number generation (in this case, reproducibly
#'   random subsampling). If set to FALSE, then no fiddling with the RNG seed
#'   is performed, and it is up to the user to appropriately call set.seed
#'   beforehand to achieve reproducible results. Default is FALSE.
#' @param verbose (logical). If TRUE, print additional information.
#' @param one_plot (logical, default FALSE) If TRUE, return a unique
#'   plot with the three plot inside using the patchwork package.
#' @param ... Additional arguments passed on to [ggstatsplot::ggscatterstats()]
#'   function.
#'
#' @return Either an unique ggplot2 (when `one_plot` is TRUE) or
#'  a list of ggplot2 plot for each q.
#' @export
#' @author Adrien Taudière
#'
#' @examples
#' if (requireNamespace("ggstatsplot")) {
#'   library("divent")
#'   ggscatt_pq(data_fungi_mini, "Time", q = 0, type = "non-parametric")
#' }
#' \donttest{
#' if (requireNamespace("ggstatsplot")) {
#'   ggscatt_pq(data_fungi_mini, "Sample_id",
#'     q = 0,
#'     one_plot = FALSE
#'   )
#' }
#' }
#' @details
#' This function is mainly a wrapper of the work of others.
#'   Please make a reference to `ggstatsplot::ggscatterstats()` if you
#'   use this function.
#' @seealso [ggbetween_pq()]
ggscatt_pq <- function(
  physeq,
  num_modality,
  q = c(0, 1, 2),
  rarefy_by_sample = FALSE,
  rngseed = FALSE,
  verbose = TRUE,
  one_plot = TRUE,
  ...
) {
  verify_pq(physeq)
  physeq <- clean_pq(physeq, force_taxa_as_columns = TRUE)

  if (rarefy_by_sample) {
    if (as(rngseed, "logical")) {
      set.seed(rngseed)
      if (verbose) {
        message(
          "`set.seed(",
          rngseed,
          ")` was used to initialize repeatable random subsampling."
        )
        message("Please record this for your records so others can reproduce.")
        message(
          "Try `set.seed(",
          rngseed,
          "); .Random.seed` for the full vector",
          sep = ""
        )
        message("...")
      }
    } else if (verbose) {
      message(
        "You set `rngseed` to FALSE. Make sure you've set & recorded\n",
        " the random seed of your session for reproducibility.\n",
        "See `?set.seed`\n"
      )
      message("...")
    }
    physeq <- clean_pq(rarefy_even_depth_pq(physeq, rngseed = rngseed))
  }

  p_list <- vector("list", length(q))
  psm_res <- psmelt_samples_pq(physeq, q = q)
  for (i in seq_along(q)) {
    p_list[[i]] <-
      ggstatsplot::ggscatterstats(
        psm_res,
        !!paste0("Hill_", q[[i]]),
        !!num_modality,
        ...
      )
  }

  if (one_plot) {
    return(patchwork::wrap_plots(p_list))
  } else {
    return(p_list)
  }
}
################################################################################

################################################################################
#' Alluvial plot for taxonomy and samples factor vizualisation
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' Basically a wrapper of [ggalluvial](https://corybrunson.github.io/ggalluvial/index.html)
#' package
#'
#' @inheritParams clean_pq
#' @param taxa_ranks A vector of taxonomic ranks. For examples c("Family","Genus").
#'   If taxa ranks is not set
#'   (default value = c("Phylum", "Class", "Order", "Family")).
#' @param wrap_factor A name to determine
#'   which samples to merge using [merge_samples2()] function.
#'   Need to be in \code{physeq@sam_data}.
#'   Need to be use when you want to wrap by factor the final plot
#'   with the number of taxa (type="nb_taxa")
#' @param by_sample (logical) If FALSE (default), sample information is not taking
#'   into account, so the taxonomy is studied globally. If fact is not NULL, by_sample
#'   is automatically set to TRUE.
#' @param rarefy_by_sample (logical, default FALSE) If TRUE, rarefy
#'   samples using [phyloseq::rarefy_even_depth()] function.
#' @param rngseed (Optional). A single integer value passed to
#'   [phyloseq::rarefy_even_depth()], which is used to fix a seed for
#'   reproducibly random number generation (in this case, reproducibly
#'   random subsampling). If set to FALSE, then no fiddling with the RNG seed
#'   is performed, and it is up to the user to appropriately call set.seed
#'   beforehand to achieve reproducible results. Default is FALSE.
#' @param verbose (logical). If TRUE, print additional information.
#' @param fact (required) Name of the factor in `physeq@sam_data` used to plot  the last column
#' @param type If "nb_seq" (default), the number of sequences is
#'   used in plot. If "nb_taxa", the number of ASV is plotted.
#' @param width (passed on to [ggalluvial::geom_flow()]) the width of each stratum,
#'   as a proportion of the distance between axes. Defaults to 1/3.
#' @param min.size (passed on to [ggfittext::geom_fit_text()]) Minimum font size,
#'   in points. Text that would need to be shrunk below this size to fit the box will
#'   be hidden. Defaults to 4 pt.
#' @param na_remove (logical, default FALSE) If set to TRUE, remove samples with
#'   NA in the variables set in formula.
#' @param use_ggfittext (logical, default FALSE) Do we use ggfittext to plot labels?
#' @param use_geom_label (logical, default FALSE) Do we use geom_label to plot labels?
#' @param size_lab Size for label if use_ggfittext is FALSE
#' @param ... Additional arguments passed on to [ggalluvial::geom_flow()] function.
#'
#' @return A ggplot object
#' @export
#' @author Adrien Taudière
#' @examples
#' if (requireNamespace("ggalluvial")) {
#'   ggaluv_pq(data_fungi_mini)
#' }
#' \donttest{
#' if (requireNamespace("ggalluvial")) {
#'   library(ggalluvial)
#'   ggaluv_pq(data_fungi_mini)
#'
#'   ggaluv_pq(data_fungi_mini, type = "nb_taxa") +
#'     geom_text(stat = "stratum", size = 1.8)
#'
#'   ggaluv_pq(data_fungi_mini,
#'     wrap_factor = "Height",
#'     by_sample = TRUE,
#'     type = "nb_taxa"
#'   ) +
#'     facet_wrap("Height")
#'
#'   ggaluv_pq(data_fungi_mini,
#'     width = 0.9, min.size = 10,
#'     type = "nb_taxa", taxa_ranks = c("Phylum", "Class", "Order", "Family", "Genus")
#'   ) + coord_flip() +
#'     scale_x_discrete(limits = rev)
#' }
#' }
#' @details
#' This function is mainly a wrapper of the work of others.
#'   Please make a reference to `ggalluvial` package if you
#'   use this function.
#'
#' When you want to add text to the plot, this function requires
#'  ggalluvial to be loaded with before use (`library(ggalluvial)`).
#' @seealso [sankey_pq()]
ggaluv_pq <- function(
  physeq,
  taxa_ranks = c("Phylum", "Class", "Order", "Family"),
  wrap_factor = NULL,
  by_sample = FALSE,
  rarefy_by_sample = FALSE,
  rngseed = FALSE,
  verbose = TRUE,
  fact = NULL,
  type = "nb_seq",
  width = 1.2,
  min.size = 3,
  na_remove = FALSE,
  use_ggfittext = FALSE,
  use_geom_label = FALSE,
  size_lab = 2,
  ...
) {
  verify_pq(physeq)
  if (rarefy_by_sample) {
    if (as(rngseed, "logical")) {
      set.seed(rngseed)
      if (verbose) {
        message(
          "`set.seed(",
          rngseed,
          ")` was used to initialize repeatable random subsampling."
        )
        message("Please record this for your records so others can reproduce.")
        message(
          "Try `set.seed(",
          rngseed,
          "); .Random.seed` for the full vector",
          sep = ""
        )
        message("...")
      }
    } else if (verbose) {
      message(
        "You set `rngseed` to FALSE. Make sure you've set & recorded\n",
        " the random seed of your session for reproducibility.\n",
        "See `?set.seed`\n"
      )
      message("...")
    }
    physeq <- rarefy_even_depth_pq(physeq, rngseed = rngseed)
  }

  if (na_remove && !is.null(fact)) {
    physeq <- subset_samples_pq(physeq, !is.na(physeq@sam_data[[fact]]))
  }

  if (!is.null(wrap_factor)) {
    physeq <-
      merge_samples2(physeq, physeq@sam_data[[wrap_factor]])
  } else if (!by_sample || !is.null(fact)) {
    physeq <-
      merge_samples2(
        physeq,
        group = rep("all_samples_together", nsamples(physeq))
      )
  }

  if (type == "nb_taxa") {
    physeq <- as_binary_otu_table(physeq)
  } else if (type != "nb_seq") {
    stop("Type must be eiter nb_seq or nb_taxa")
  }

  psm_samp <-
    psmelt_samples_pq(
      physeq,
      taxa_ranks = taxa_ranks,
      q = NULL,
      rarefy_by_sample = FALSE
    )

  if (is.null(fact)) {
    psm_samp <- ggalluvial::to_lodes_form(psm_samp, axes = taxa_ranks)
  } else {
    psm_samp <- ggalluvial::to_lodes_form(psm_samp, axes = c(taxa_ranks, fact))
  }

  p <- ggplot(
    data = psm_samp,
    aes(
      alluvium = alluvium,
      x = x,
      stratum = stratum,
      y = Abundance,
      fill = after_stat(stratum),
      label = after_stat(stratum)
    )
  ) +
    ggalluvial::geom_flow(...) +
    ggalluvial::geom_stratum() +
    theme_minimal() +
    theme(legend.position = "none")

  if (use_ggfittext) {
    if (!"package:ggalluvial" %in% search()) {
      message("Please load ggalluvial with: library(ggalluvial)")
      stop("ggalluvial must be attached when use_ggfittext is TRUE")
    }
    p <- p +
      ggalluvial::geom_stratum() +
      ggfittext::geom_fit_text(
        stat = "stratum",
        width = width,
        min.size = min.size
      )
  } else if (use_geom_label) {
    if (!"package:ggalluvial" %in% search()) {
      message("Please load ggalluvial with: library(ggalluvial)")
      stop("ggalluvial must be attached when use_geom_label is TRUE")
    }
    p <- p +
      ggalluvial::geom_stratum() +
      geom_label(,
        stat = "stratum",
        size = size_lab
      )
  }

  if (!is.null(wrap_factor)) {
    p <- p + facet_wrap(wrap_factor)
  }
  return(p)
}
################################################################################

################################################################################
#' Plot the nucleotide proportion at both extremity of the sequences
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#'   It is a useful function to check for the absence of unwanted patterns caused
#'   for example by Illumina adaptator or bad removal of primers.
#'
#'   If `hill_scale` is not null, Hill diversity number are used to represent the distribution
#'   of the diversity (equitability) along the sequences.
#'
#' @inheritParams clean_pq
#' @param first_n (int, default 10) The number of nucleotides to plot the 5' extremity.
#' @param last_n (int, default 10) The number of nucleotides to plot the 3' extremity.
#' @param q (vector) A vector defining the Hill number wanted. Set to NULL if
#'   you don't want to plot Hill diversity metrics. Hill numbers are more
#'   appropriate in DNA metabarcoding studies when `q > 0` (Alberdi & Gilbert,
#'   2019; Calderón-Sanou et al., 2019).
#' @param min_width (int, default 0) Select only the sequences from physeq@refseq with using a
#'   minimum length threshold. If `first_n` is superior to the minimum length of the
#'   references sequences, you must use min_width to filter out the narrower sequences
#' @return A list of 4 objects
#'  - p_start and p_last are the ggplot object representing respectively the start and
#'   the end of the sequences.
#'  - df_start and df_last are the data.frame corresponding to the ggplot object.
#' @export
#' @author Adrien Taudière
#' @examples
#' data_f <- prune_samples(
#'   sample_names(data_fungi_mini)[1:20],
#'   data_fungi_mini
#' )
#' library("divent")
#' res1 <- plot_refseq_extremity_pq(data_f, q = 1)
#' names(res1)
#' \donttest{
#' res1$plot_start
#' res1$plot_last
#'
#' res2 <- plot_refseq_extremity_pq(data_f, first_n = 200, last_n = 100)
#' res2$plot_start
#' res2$plot_last
#'
#' plot_refseq_extremity_pq(data_f,
#'   first_n = NULL,
#'   last_n = 200,
#'   min_width = 200,
#'   q = c(3)
#' )$plot_last
#' }
plot_refseq_extremity_pq <- function(
  physeq,
  first_n = 10,
  last_n = 10,
  q = c(1, 2),
  min_width = 0
) {
  if (min_width > 0) {
    cond <- Biostrings::width(physeq@refseq) > min_width
    names(cond) <- taxa_names(physeq)
    physeq <- clean_pq(subset_taxa_pq(physeq, cond))
  }

  if (!is.null(first_n)) {
    end_n <- Biostrings::width(physeq@refseq)
    end_n[end_n < first_n] <- first_n

    letters_sequences <-
      c(strsplit(
        as.character(IRanges::narrow(physeq@refseq, end = end_n)),
        split = ""
      ))
    letters_sequences <- lapply(
      letters_sequences,
      `length<-`,
      max(lengths(letters_sequences))
    )

    tib_interm <- t(data.frame(letters_sequences))
    nucleotide_first_interm <- data.frame(
      "nb_A" = colSums(tib_interm == "A", na.rm = TRUE) / nrow(tib_interm),
      "nb_C" = colSums(tib_interm == "C", na.rm = TRUE) / nrow(tib_interm),
      "nb_G" = colSums(tib_interm == "G", na.rm = TRUE) / nrow(tib_interm),
      "nb_T" = colSums(tib_interm == "T", na.rm = TRUE) / nrow(tib_interm),
      "seq_id" = seq_len(ncol(tib_interm))
    ) |>
      rowwise() |>
      mutate("max_letter_prob" = max(across(starts_with("nb_"))))

    nucleotide_first <- nucleotide_first_interm |>
      tidyr::pivot_longer(cols = starts_with("nb_"))

    p_start <- ggplot(nucleotide_first) +
      geom_point(aes(x = seq_id, y = value, color = name)) +
      xlim(c(0, first_n)) +
      labs(
        subtitle = paste0(
          "Proportion of nucleotide along the ",
          first_n,
          " first nucleotides \n in sequences of ",
          ntaxa(physeq),
          " taxa representing ",
          sum(physeq@otu_table),
          " sequences in ",
          nsamples(physeq),
          " samples."
        )
      )

    if (!is.null(q)) {
      suppressMessages(
        hill_nucleotide <- divent_hill_matrix_pq(
          nucleotide_first_interm[, c("nb_A", "nb_C", "nb_G", "nb_T")],
          q = q
        )
      )
      if (sum(q == 0) > 0) {
        hill_nucleotide <- hill_nucleotide |>
          rename("Hill 0 (Richness)" = "0")
      }
      if (sum(q == 1) > 0) {
        hill_nucleotide <- hill_nucleotide |>
          rename("Hill 1 (Shannon)" = "1")
      }
      if (sum(q == 2) > 0) {
        hill_nucleotide <- hill_nucleotide |>
          rename("Hill 2 (Simpson)" = "2")
      }

      hill_nucleotide <- hill_nucleotide |>
        tidyr::pivot_longer(cols = everything())

      hill_nucleotide$seq_id <-
        sort(rep(
          seq_len(ncol(tib_interm)),
          times = nrow(hill_nucleotide) / ncol(tib_interm)
        ))

      p_start <- p_start +
        geom_line(
          data = hill_nucleotide,
          aes(x = seq_id, y = value, color = name)
        )
    }
  } else {
    p_start <- NULL
    nucleotide_first_interm <- NULL
  }

  if (!is.null(last_n)) {
    tib_interm_last <- t(data.frame(
      letters = c(
        strsplit(
          as.vector(
            IRanges::narrow(
              physeq@refseq,
              start = Biostrings::width(physeq@refseq) - last_n
            )
          ),
          split = ""
        )
      )
    ))

    nucleotide_last_interm <- data.frame(
      "nb_A" = colSums(tib_interm_last == "A") / nrow(tib_interm_last),
      "nb_C" = colSums(tib_interm_last == "C") / nrow(tib_interm_last),
      "nb_G" = colSums(tib_interm_last == "G") / nrow(tib_interm_last),
      "nb_T" = colSums(tib_interm_last == "T") / nrow(tib_interm_last),
      "seq_id" = seq_len(ncol(tib_interm_last))
    )

    nucleotide_last <- nucleotide_last_interm |>
      tidyr::pivot_longer(cols = starts_with("nb_"))

    p_last <- ggplot(nucleotide_last) +
      geom_point(aes(x = seq_id, y = value, color = name)) +
      labs(
        subtitle = paste0(
          "Proportion of nucleotide along the ",
          last_n,
          " last nucleotides \n in sequences of ",
          ntaxa(physeq),
          " taxa representing ",
          sum(physeq@otu_table),
          " sequences in ",
          nsamples(physeq),
          " samples."
        )
      )

    if (!is.null(q)) {
      suppressMessages(
        hill_nucleotide <- divent_hill_matrix_pq(
          nucleotide_last_interm[, c("nb_A", "nb_C", "nb_G", "nb_T")],
          q = q
        )
      )
      if (sum(q == 0) > 0) {
        hill_nucleotide <- hill_nucleotide |>
          rename("Hill 0 (Richness)" = "0")
      }
      if (sum(q == 1) > 0) {
        hill_nucleotide <- hill_nucleotide |>
          rename("Hill 1 (Shannon)" = "1")
      }
      if (sum(q == 2) > 0) {
        hill_nucleotide <- hill_nucleotide |>
          rename("Hill 2 (Simpson)" = "2")
      }

      hill_nucleotide <- hill_nucleotide |>
        tidyr::pivot_longer(cols = everything())

      hill_nucleotide$seq_id <-
        sort(rep(
          seq_len(ncol(tib_interm_last)),
          times = nrow(hill_nucleotide) / ncol(tib_interm_last)
        ))

      p_last <- p_last +
        geom_line(
          data = hill_nucleotide,
          aes(
            x = seq_id,
            y = value,
            color = name
          )
        )
    }
  } else {
    p_last <- NULL
    nucleotide_last_interm <- NULL
  }

  return(list(
    "plot_start" = p_start,
    "plot_last" = p_last,
    "df_start" = nucleotide_first_interm,
    "df_end" = nucleotide_last_interm
  ))
}
################################################################################

################################################################################
#' Plot the nucleotide proportion of references sequences
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#'  It is a wrapper of the function `plot_refseq_extremity_pq()`. See
#'  ?plot_refseq_extremity_pq for more examples.
#'
#'   If `hill_scale` is not null, Hill diversity number are used to represent the distribution
#'   of the diversity (equitability) along the sequences.
#'
#' @inheritParams clean_pq
#' @param first_n (int, default 10) The number of nucleotides to plot the 5' extremity.
#' @param last_n (int, default 10) The number of nucleotides to plot the 3' extremity.
#' @param q (vector) A vector defining the Hill number wanted. Set to NULL if
#'   you don't want to plot Hill diversity metrics. Hill numbers are more
#'   appropriate in DNA metabarcoding studies when `q > 0` (Alberdi & Gilbert,
#'   2019; Calderón-Sanou et al., 2019).
#' @param min_width (int, default 0) Select only the sequences from physeq@refseq with using a
#'   minimum length threshold. If `first_n` is superior to the minimum length of the
#'   references sequences, you must use min_width to filter out the narrower sequences
#' @return A ggplot2 object
#' @export
#' @author Adrien Taudière
#' @examples
#' plot_refseq_pq(data_fungi_mini)
#' \dontrun{
#' plot_refseq_pq(data_fungi_mini, q = c(2), first_n = 300)
#' }
#'
plot_refseq_pq <- function(
  physeq,
  q = NULL,
  first_n = min(Biostrings::width(physeq@refseq)),
  last_n = NULL,
  min_width = first_n
) {
  res <- plot_refseq_extremity_pq(
    physeq,
    q = q,
    first_n = first_n,
    last_n = last_n,
    min_width = min_width
  )
  return(res$plot_start)
}


################################################################################
#' Discard legend in ggplot2
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-stable-green" alt="lifecycle-stable"></a>
#'
#'  A more memorable shortcut for theme(legend.position = "none").
#'
#' @export
#' @return A ggplot2 object
#' @author Adrien Taudière
#' @examples
#' plot_refseq_pq(data_fungi_mini)
#' plot_refseq_pq(data_fungi_mini) + no_legend()
no_legend <- function() {
  list(theme(legend.position = "none"))
}
################################################################################

################################################################################
#' Hill Diversities and Corresponding Accumulation Curves for phyloseq
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#'   Basically a wrapper of [vegan::renyi()] and
#'   [vegan::renyiaccum()] functions
#'
#' @inheritParams clean_pq
#' @param merge_sample_by a vector to determine which samples to merge using
#'   the [merge_samples2()] function.  Need to be in `physeq@sam_data`
#' @param color_fac (optional): The variable to color the barplot. For ex.
#'   same as fact. If merge_sample_by is set, color_fac must be nested in
#'   the mq_by factor. See examples.
#' @param q Scales of Rényi diversity.
#' @param nperm (int Default NULL) If a integer is set to nperm, nperm
#'   permutation are computed to draw confidence interval for each curves.
#'   The function use [vegan::renyi()] if nperm is NULL and
#'   [vegan::renyiaccum()] else.
#' @param na_remove (logical, default FALSE) If set to TRUE, remove samples with
#'   NA in the variables set in merge_sample_by. Not used if merge_sample_by is
#'   NULL.
#' @param wrap_factor (logical, default TRUE) Do the plot is wrap by the factor
#' @param plot_legend (logical, default TRUE) If set to FALSE,
#'   no legend are plotted.
#' @param linewidth (int, default 2) The linewidth of lines.
#' @param size_point (int, default 1) The size of the point.
#'
#' @param ... Additional arguments passed on to [vegan::renyi()] function or
#'   [vegan::renyiaccum()] if nperm is not NULL.
#'
#' @export
#' @author Adrien Taudière
#' @return A ggplot2 object
#' @examples
#' \donttest{
#' if (requireNamespace("vegan")) {
#'   hill_curves_pq(data_fungi_mini, merge_sample_by = "Time")
#'   hill_curves_pq(data_fungi_mini, color_fac = "Time", plot_legend = FALSE)
#'   hill_curves_pq(data_fungi_mini,
#'     color_fac = "Time", plot_legend = FALSE,
#'     nperm = 9, size_point = 1, linewidth = 0.5
#'   )
#'
#'   hill_curves_pq(data_fungi_mini,
#'     nperm = 9, plot_legend = FALSE, size_point = 1,
#'     linewidth = 0.5
#'   )
#'   hill_curves_pq(data_fungi_mini, "Height",
#'     q = c(0, 1, 2, 8), plot_legend = FALSE
#'   )
#'   hill_curves_pq(data_fungi_mini, "Height",
#'     q = c(0, 0.5, 1, 2, 4, 8),
#'     nperm = 9
#'   )
#'   hill_curves_pq(data_fungi_mini, "Height", nperm = 9, wrap_factor = FALSE)
#'
#'   data_fungi_mini@sam_data$H_T <- paste0(
#'     data_fungi_mini@sam_data$Height,
#'     "_", data_fungi_mini@sam_data$Time
#'   )
#'   merge_samples2(data_fungi_mini, "H_T")
#'   hill_curves_pq(data_fungi_mini, "H_T", color_fac = "Time", nperm = 9)
#' }
#' }
#' @details
#' This function is mainly a wrapper of the work of others.
#'   Please make a reference to [vegan::renyi()] or
#'   [vegan::renyiaccum()] functions
#'
hill_curves_pq <- function(
  physeq,
  merge_sample_by = NULL,
  color_fac = NULL,
  q = c(0, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, Inf),
  nperm = NULL,
  na_remove = TRUE,
  wrap_factor = TRUE,
  plot_legend = TRUE,
  linewidth = 2,
  size_point = 2,
  ...
) {
  verify_pq(physeq)

  if (na_remove && !is.null(merge_sample_by)) {
    new_physeq <-
      subset_samples_pq(physeq, !is.na(physeq@sam_data[[merge_sample_by]]))
    if (nsamples(physeq) - nsamples(new_physeq) > 0) {
      message(
        paste0(
          nsamples(physeq) - nsamples(new_physeq),
          " were discarded due to NA in variables present in formula."
        )
      )
    }
    physeq <- new_physeq
  }
  if (!is.null(merge_sample_by)) {
    physeq <- merge_samples2(physeq, merge_sample_by)
  }

  physeq <- clean_pq(
    physeq,
    force_taxa_as_rows = TRUE,
    remove_empty_samples = FALSE,
    remove_empty_taxa = FALSE,
    clean_samples_names = FALSE
  )

  otu_mat <- as(t(physeq)@otu_table, "matrix")

  if (!is.null(nperm)) {
    df_hill <-
      vegan::renyiaccum(
        otu_mat,
        scales = q,
        permutation = nperm,
        hill = TRUE,
        ...
      )
    what <- c("Collector", "mean", "Qnt 0.025", "Qnt 0.975")
    what <- what[what %in% dimnames(df_hill)[[3]]]
    if (any(what %in% dimnames(df_hill)[[3]])) {
      df_hill <- df_hill[,, what, drop = FALSE]
    }
    dm <- dim(df_hill)
    dnam <- dimnames(df_hill)
    lin <- rep(dnam[[3]], each = dm[1] * dm[2])
    alp <- factor(dnam[[2]], levels = dnam[[2]])
    alpha <- rep(rep(alp, each = dm[1]), len = prod(dm))
    diversity <- as.vector(df_hill)

    if (is.null(color_fac)) {
      modality <- rep(sample_names(physeq), len = prod(dm))
    } else {
      modality <- rep(
        levels(as.factor(physeq@sam_data[, color_fac][[1]])),
        len = prod(dm)
      )
    }

    samp_names <- rep(sample_names(physeq), len = prod(dm))

    df_plot <- data.frame(
      diversity = diversity,
      Type = lin,
      Modality = modality,
      alpha_hill = alpha,
      samp_names = samp_names
    )

    p <- ggplot(
      df_plot,
      aes(
        y = diversity,
        x = alpha_hill,
        color = Modality,
        group = samp_names
      )
    ) +
      geom_point(data = subset(df_plot, Type == "mean"), size = size_point) +
      geom_line(
        data = subset(df_plot, Type == "mean"),
        linetype = 1,
        linewidth = linewidth
      ) +
      geom_line(
        data = subset(df_plot, Type == "Qnt 0.025"),
        linetype = 2,
        linewidth = linewidth / 3
      ) +
      geom_line(
        data = subset(df_plot, Type == "Qnt 0.975"),
        linetype = 2,
        linewidth = linewidth / 3
      )
    if (wrap_factor) {
      p <- p +
        facet_wrap("Modality")
    }
  } else {
    df_hill <- vegan::renyi(otu_mat, scales = q, hill = TRUE)
    if (inherits(df_hill, "data.frame")) {
      if (is.null(color_fac)) {
        modality <- sample_names(physeq)
      } else {
        modality <- factor(
          rep(rownames(df_hill), ncol(df_hill)),
          levels = rownames(df_hill)
        )
      }

      alp <- factor(
        rep(colnames(df_hill), each = nrow(df_hill)),
        levels = colnames(df_hill)
      )
      div <- as.vector(as.matrix(df_hill))
      df_plot <- data.frame(
        diversity = div,
        Modality = modality,
        alpha_hill = alp
      )
    } else {
      df_plot <- data.frame(
        diversity = x,
        alpha = factor(names(x), levels = names(x)),
        plot = "plot"
      )
      lo <- hi <- med <- NA
    }
    p <- ggplot(
      df_plot,
      aes(
        y = diversity,
        x = alpha_hill,
        color = Modality,
        group = Modality
      )
    ) +
      geom_point(size = 2) +
      geom_line()
  }

  if (!plot_legend) {
    p <- p + no_legend()
  }
  return(p)
}
################################################################################

################################################################################
#' Computes a manifold approximation and projection (UMAP) for
#' phyloseq object
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' https://journals.asm.org/doi/full/10.1128/msystems.00691-21
#'
#' @inheritParams clean_pq
#' @param pkg Which R packages to use, either "umap" or "uwot".
#' @param ... Additional arguments passed on to [umap::umap()] or
#'   [uwot::umap2()] function.
#'   For example `n_neighbors` set the number of nearest neighbors (Default 15).
#'   See [umap::umap.defaults()] or [uwot::umap2()] for the list of
#'   parameters and default values.
#'
#' @return A dataframe with samples informations and the x_umap and y_umap position
#' @author Adrien Taudière
#' @export
#' @seealso [umap::umap()], [tsne_pq()], [phyloseq::plot_ordination()]
#' @examples
#' library("umap")
#' data_f <- prune_samples(
#'   sample_names(data_fungi_mini)[1:20],
#'   data_fungi_mini
#' )
#' df_umap <- umap_pq(data_f, n_neighbors = 3)
#' ggplot(df_umap, aes(x = x_umap, y = y_umap, col = Height)) +
#'   geom_point(size = 2)
#'
#' \dontrun{
#' df_uwot <- umap_pq(data_fungi_mini, pkg = "uwot")
#' library(patchwork)
#' physeq <- data_fungi_mini
#' df_umap <- umap_pq(physeq, n_neighbors = 3)
#' res_tsne <- tsne_pq(data_fungi_mini)
#' df_umap_tsne <- df_umap
#' df_umap_tsne$x_tsne <- res_tsne$Y[, 1]
#' df_umap_tsne$y_tsne <- res_tsne$Y[, 2]
#' ((ggplot(df_umap, aes(x = x_umap, y = y_umap, col = Height)) +
#'   geom_point(size = 2) +
#'   ggtitle("UMAP")) +
#'   (plot_ordination(physeq,
#'     ordination = ordinate(physeq, method = "PCoA", distance = "bray"),
#'     color = "Height"
#'   ) + ggtitle("PCoA"))) /
#'   ((ggplot(df_umap_tsne, aes(x = x_tsne, y = y_tsne, col = Height)) +
#'     geom_point(size = 2) +
#'     ggtitle("tsne")) +
#'     (plot_ordination(physeq,
#'       ordination = ordinate(physeq, method = "NMDS", distance = "bray"),
#'       color = "Height"
#'     ) + ggtitle("NMDS"))) +
#'   patchwork::plot_layout(guides = "collect")
#'
#' (ggplot(df_umap, aes(x = x_umap, y = y_umap, col = Height)) +
#'   geom_point(size = 2) +
#'   ggtitle("umap::umap")) /
#'   (ggplot(df_uwot, aes(x = x_umap, y = y_umap, col = Height)) +
#'     geom_point(size = 2) +
#'     ggtitle("uwot::umap2"))
#' }
#' @details
#' This function is mainly a wrapper of the work of others.
#'   Please make a reference to `umap::umap()` if you
#'   use this function.

umap_pq <- function(physeq, pkg = "umap", ...) {
  verify_pq(physeq)
  physeq <- MiscMetabar::taxa_as_columns(physeq)

  psm_samp <- psmelt_samples_pq(physeq)
  if (pkg == "umap") {
    res_umap <- umap::umap(as.matrix(unclass(physeq@otu_table)), ...)
    umap_layout <- as_tibble(res_umap$layout, .name_repair = "minimal")
    umap_layout$Sample <- rownames(res_umap$layout)
    names(umap_layout) <- c("x_umap", "y_umap", "Sample")
  } else if (pkg == "uwot") {
    res_umap <- uwot::umap2(as.matrix(unclass(physeq@otu_table)), ...)
    umap_layout <- as_tibble(res_umap, .name_repair = c("minimal"))
    umap_layout$Sample <- rownames(res_umap)
    names(umap_layout) <- c("x_umap", "y_umap", "Sample")
  } else {
    stop("Param pkg must be set to 'umap' or 'uwot'.")
  }

  df_umap <- left_join(umap_layout, psm_samp)

  return(df_umap)
}
################################################################################

################################################################################
#' Plot kmer complexity of references sequences of a phyloseq object
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#'   Basically a wrapper of [dada2::seqComplexity()]
#'
#' @inheritParams clean_pq
#' @param kmer_size int (default 2) The size of the kmers
#'   (or "oligonucleotides" or "words") to use.
#' @param window (int, default NULL) The width in nucleotides of the moving
#'    window. If NULL the whole sequence is used.
#' @param by (int, default 5) The step size in nucleotides between each moving
#'    window tested.
#' @param bins (int, default 100). The number of bins to use for the histogram.
#' @param aggregate (logical, default FALSE) If TRUE, compute an aggregate quality profile
#'    for all samples
#' @param vline_random_kmer (logical, default TRUE) If TRUE, add a vertical line
#'   at the value for random kmer (equal to 4^kmerSize))
#' @param ... Arguments passed on to geom_histogram.
#'
#' @return A ggplot2 object
#' @export
#' @author Adrien Taudière
#' @seealso [dada2::seqComplexity()], [dada2::plotComplexity()]
#' @examples
#' plot_complexity_pq(subset_samples(data_fungi_mini, Height == "High"),
#'   vline_random_kmer = FALSE
#' )
#' # plot_complexity_pq(subset_samples(data_fungi_mini, Height == "Low"),
#' #  aggregate = FALSE, kmer_size = 4
#' # )
#' # plot_complexity_pq(subset_samples(data_fungi, Height == "Low"),
#' #  kmer_size = 4)
#'
#' @details
#' This function is mainly a wrapper of the work of others.
#'   Please make a reference to [dada2::seqComplexity()]

plot_complexity_pq <- function(
  physeq,
  kmer_size = 2,
  window = NULL,
  by = 5,
  bins = 100,
  aggregate = FALSE,
  vline_random_kmer = TRUE,
  ...
) {
  if (aggregate) {
    refseq_complex <- dada2::seqComplexity(
      physeq@refseq,
      kmerSize = kmer_size,
      window = window,
      by = by
    )
  } else {
    refseq_complex <- vector("list", length = nsamples(physeq))
    names(refseq_complex) <- sample_names(physeq)
    for (sam in sample_names(physeq)) {
      physeq_interm <- subset_samples_pq(physeq, sample_names(physeq) == sam)
      refseq_complex[[sam]] <- dada2::seqComplexity(
        physeq_interm@refseq,
        kmerSize = kmer_size,
        window = window,
        by = by
      )
    }
    df <- data.frame(
      complexity = unlist(refseq_complex),
      file = rep(sample_names(physeq), times = sapply(refseq_complex, length))
    )
  }

  p <- ggplot(data = df, aes(x = complexity)) +
    geom_histogram(
      bins = bins,
      na.rm = TRUE,
      ...
    ) +
    ylab("Count") +
    xlab("Effective Oligonucleotide Number") +
    theme_bw() +
    facet_wrap(~file) +
    scale_x_continuous(
      limits = c(0, 4^kmer_size),
      breaks = seq(0, 4^kmer_size, (4^kmer_size) / 4)
    )
  if (vline_random_kmer) {
    p <- p +
      geom_vline(xintercept = 4^kmer_size, color = "red")
  }
  return(p)
}
################################################################################

################################################################################
#' A diagnostic plot of the number of sequences per samples
#'
#' @inheritParams clean_pq
#' @param min_nb_seq (int) The minimum number of sequences per samples to compare
#'   the ratio.
#' @param annotations (logical, default TRUE). If FALSE, no annotations are
#'   plotted
#'
#' @returns A ggplot2 object
#' @export
#' @author Adrien Taudière
#' @details The x axis depict the number of sequences per samples and the y
#'   axis depicted the ratio of the number of sequences for a given sample
#'   divide by the number of sequences of the previous sample when ordered by
#'   the number of sequences. A high ratio indicate an important and quick
#'   increase of the number of sequence which may indicate that below this
#'   ratio, samples are suspicious.
#'
#'   The general idea is to first removed all samples with definitively not
#'   enough sequences and then, among the kept samples, find the higher
#'   augmentation (ratio) to possibly detect suspicious samples.
#'
#' @examples
#' plot_seq_ratio_pq(data_fungi_mini, min_nb_seq = 10, annotations = FALSE)
#' \donttest{
#' plot_seq_ratio_pq(data_fungi, min_nb_seq = 200)
#' data(GlobalPatterns)
#' plot_seq_ratio_pq(GlobalPatterns, min_nb_seq = 100000)
#' }
plot_seq_ratio_pq <- function(physeq, min_nb_seq = 1000, annotations = TRUE) {
  if (min_nb_seq < min(sample_sums(physeq))) {
    stop(
      "You must specify a min_nb_seq below the minimum value of sample_sums in your phyloseq object."
    )
  }

  cutof_index <- sort(sample_sums(physeq)) |>
    as_tibble(.name_repair = c("minimal")) |>
    mutate(diff = c(0, diff(value))) |>
    filter(value < min_nb_seq) |>
    pull(diff) |>
    which.max()

  n_cutoff <- cutof_index

  df <- tibble(
    "value" = sort(sample_sums(physeq)),
    "name" = names(sort(sample_sums(physeq)))
  ) |>
    mutate(diff = c(0, diff(value))) |>
    mutate(ratio = value / dplyr::lag(value, default = 0))

  cutof_value <- df$value[-c(1:n_cutoff)][which.max(df$ratio[-c(1:n_cutoff)])]
  cutof_ratio <- max(df$ratio[-c(1:n_cutoff)])

  df <- df |>
    mutate(
      color_group = case_when(
        value >= cutof_value ~ "keep",
        value >= min_nb_seq ~ "suspicious",
        .default = "discard"
      )
    )

  p <- ggplot(df) +
    geom_point(
      aes(x = value, y = ratio, color = color_group),
      size = 2,
      alpha = 0.8
    ) +
    scale_x_log10() +
    geom_vline(xintercept = cutof_value, alpha = 0.8, color = "darkgreen") +
    geom_vline(xintercept = min_nb_seq, alpha = 0.8, color = "grey") +
    geom_hline(yintercept = cutof_ratio, alpha = 0.8, color = "darkgreen") +
    scale_color_manual(values = c("orange", "darkgreen", "grey20")) +
    guides(color = "none") +
    coord_cartesian(clip = "off")

  if (annotations) {
    p <- p +
      annotate(
        geom = "segment",
        x = min_nb_seq,
        y = 1.04 * cutof_ratio,
        xend = min(df$value),
        yend = 1.04 * cutof_ratio,
        arrow = arrow(length = unit(2, "mm"))
      ) +
      annotate(
        geom = "text",
        x = min_nb_seq,
        y = 1.08 * cutof_ratio,
        label = paste0("Samples with less than \n", min_nb_seq, " sequences"),
        hjust = "right",
        size = 3
      ) +
      annotate(
        geom = "curve",
        x = 4 * cutof_value,
        y = 1.15 * cutof_ratio,
        xend = 1.05 * cutof_value,
        yend = 1.01 * cutof_ratio,
        curvature = 0.3,
        arrow = arrow(length = unit(2, "mm"))
      ) +
      annotate(
        geom = "text",
        x = 4.1 * cutof_value,
        y = 1.16 * cutof_ratio,
        label = df$name[-c(1:n_cutoff)][which.max(df$ratio[-c(1:n_cutoff)])],
        hjust = "left",
        size = 3
      ) +
      annotate(
        geom = "segment",
        x = 0.98 * cutof_value,
        y = 0.98 * cutof_ratio,
        xend = min_nb_seq,
        yend = 0.98 * cutof_ratio,
        arrow = arrow(length = unit(2, "mm"))
      ) +
      annotate(
        geom = "text",
        x = cutof_value * 0.98,
        y = 0.94 * cutof_ratio,
        label = "Samples with \n suspicious ratio",
        hjust = "right",
        size = 3
      ) +
      xlab("Number of sequences per samples (log10)") +
      ylab("Ratio of the number of sequences with the previous sample") +
      labs(
        caption = paste0(
          "For the ratio (y-axis), a value of 2 indicate that sample i contains twice the number of sequences compared to the sample i-1 \n (samples ordered by their number of sequences). Run `subset_samples_pq(physeq, sample_sums(physeq)>=",
          cutof_value,
          ")` to keep only green point \n or `subset_samples_pq(physeq, sample_sums(physeq)>=",
          min_nb_seq,
          ")` to discarded only orange samples."
        )
      )
  }

  return(p)
}
################################################################################

################################################################################
#' Reorder fill and color scales to maximize perceptual contrast between
#' adjacent segments
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' In stacked bar plots, ggplot2's default discrete palette assigns colors
#' using level ordered (sometimes alphabetically), which often places perceptually
#' similar colors next to
#' each other. This function reassigns the **same set of colors** to factor
#' levels so that visually adjacent segments receive maximally different
#' colors. Both the fill and color scales are updated so that direct
#' labels (e.g. from `label_taxa = TRUE`) stay in sync with the bars.
#'
#' @param p A ggplot object that uses a discrete fill aesthetic. Can be
#'   omitted when using the `+` operator (e.g.
#'   `p + reorder_distinct_colors()`).
#' @param alternate_lightness (logical, default FALSE) If TRUE, darken every
#'   other level to add a luminance alternation cue on top of hue
#'   differences.
#' @param lightness_amount (numeric, default 0.15) Intensity of the
#'   lightness alternation (proportion to darken). Only used when
#'   `alternate_lightness = TRUE`.
#' @param colorblind (logical, default FALSE) If TRUE, compute perceptual
#'   distances under simulated deuteranopia so that the reordering
#'   optimizes contrast for colorblind viewers.
#'
#' @return A new ggplot object with [ggplot2::scale_fill_manual()] and
#'   (if a color scale is present) [ggplot2::scale_color_manual()]
#'   replacing the original scales. When `p` is omitted, returns an
#'   object that can be added to a ggplot with `+`.
#' @export
#' @author Adrien Taudière
#' @importFrom grDevices convertColor
#' @importFrom stats dist
#' @examples
#' p <- tax_bar_pq(data_fungi_mini, taxa = "Class", fact = "Time")
#' reorder_distinct_colors(p)
#' reorder_distinct_colors(p, colorblind = TRUE)
#' p + reorder_distinct_colors(alternate_lightness = TRUE)
#'
#' tax_bar_pq(data_fungi_mini,
#'   fact = "Height", taxa = "Order",
#'   nb_seq = FALSE, percent_bar = TRUE, label_taxa = TRUE,
#'   add_ribbon = TRUE, value_size = 7, ribbon_alpha = .6,
#'   show_values = TRUE, label_size = 4, top_label_size = 8,
#'   minimum_value_to_show = 0.05
#' ) |>
#'   reorder_distinct_colors(alternate_lightness = TRUE)
reorder_distinct_colors <- function(
  p = NULL,
  alternate_lightness = FALSE,
  lightness_amount = 0.15,
  colorblind = FALSE
) {
  spec <- structure(
    list(
      alternate_lightness = alternate_lightness,
      lightness_amount = lightness_amount,
      colorblind = colorblind
    ),
    class = "reorder_distinct_colors_spec"
  )
  if (is.null(p)) {
    return(spec)
  }
  if (!inherits(p, "gg")) {
    stop("p must be a ggplot object")
  }

  pb <- ggplot_build(p)
  fill_scale <- pb$plot$scales$get_scales("fill")
  if (is.null(fill_scale) || !fill_scale$is_discrete()) {
    stop("The plot must have a discrete fill scale")
  }

  levels <- fill_scale$get_limits()
  n <- length(levels)
  if (n <= 1) {
    return(p)
  }

  na_val <- fill_scale$na.value %||% "grey50"
  colors <- fill_scale$palette(n)
  if (is.null(names(colors))) {
    names(colors) <- levels
  }

  # Convert hex to sRGB matrix (rows = colors, cols = R/G/B in [0,1])
  rgb_mat <- t(col2rgb(colors)) / 255

  # Optionally simulate deuteranopia before computing distances
  if (colorblind) {
    # Brettel 1997 deuteranopia simulation matrix for sRGB
    deutan_mat <- matrix(
      c(
        0.625,
        0.375,
        0.0,
        0.7,
        0.3,
        0.0,
        0.0,
        0.3,
        0.7
      ),
      nrow = 3,
      byrow = TRUE
    )
    rgb_for_dist <- rgb_mat %*% t(deutan_mat)
  } else {
    rgb_for_dist <- rgb_mat
  }

  # Convert to CIE Lab for perceptual distance
  lab_mat <- convertColor(rgb_for_dist, from = "sRGB", to = "Lab")

  # Pairwise Euclidean distances in Lab space
  dist_mat <- as.matrix(dist(lab_mat))

  # Greedy reordering: start with the color having the largest mean distance
  avg_dist <- rowMeans(dist_mat)
  order_idx <- integer(n)
  order_idx[1] <- which.max(avg_dist)
  remaining <- setdiff(seq_len(n), order_idx[1])

  for (i in 2:n) {
    prev <- order_idx[i - 1]
    dists_to_prev <- dist_mat[prev, remaining]
    best <- which.max(dists_to_prev)
    order_idx[i] <- remaining[best]
    remaining <- setdiff(remaining, order_idx[i])
  }

  reordered_colors <- colors[order_idx]

  # Optional: alternate lightness (darken even, lighten odd)
  if (alternate_lightness) {
    rgb_reordered <- t(col2rgb(reordered_colors)) / 255
    for (i in seq_along(reordered_colors)) {
      if (i %% 2 == 0) {
        # Darken
        rgb_reordered[i, ] <- pmax(
          rgb_reordered[i, ] * (1 - lightness_amount),
          0
        )
      } else {
        # Lighten
        rgb_reordered[i, ] <- pmin(
          rgb_reordered[i, ] + (1 - rgb_reordered[i, ]) * lightness_amount,
          1
        )
      }
    }
    reordered_colors <- rgb(
      rgb_reordered[, 1],
      rgb_reordered[, 2],
      rgb_reordered[, 3]
    )
  }

  # Build named vector: level -> reordered color
  new_colors <- stats::setNames(reordered_colors, levels)

  # Remove existing fill scale and add the new one
  p$scales$scales <- p$scales$scales[
    !vapply(p$scales$scales, \(s) "fill" %in% s$aesthetics, logical(1))
  ]
  p <- p + scale_fill_manual(values = new_colors, na.value = na_val)

  # Also update the color scale if one exists
  color_scale <- pb$plot$scales$get_scales("colour")
  if (!is.null(color_scale) && color_scale$is_discrete()) {
    color_na_val <- color_scale$na.value %||% "grey50"
    p$scales$scales <- p$scales$scales[
      !vapply(
        p$scales$scales,
        \(s) "colour" %in% s$aesthetics,
        logical(1)
      )
    ]
    p <- p + scale_color_manual(values = new_colors, na.value = color_na_val)
  }

  p
}

#' @exportS3Method ggplot2::ggplot_add
ggplot_add.reorder_distinct_colors_spec <- function(object, plot, ...) {
  reorder_distinct_colors(
    p = plot,
    alternate_lightness = object$alternate_lightness,
    lightness_amount = object$lightness_amount,
    colorblind = object$colorblind
  )
}
################################################################################

################################################################################
#' A wrapper of plot_ordination with vegan distance matrix
#'
#' @details
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#'   Basically a wrapper of [phyloseq::plot_ordination()] to use aitchison and
#'  robust.aitchison distances from vegan package.
#'
#' @inheritParams clean_pq
#' @param method (string, default "robust.aitchison") The distance method to use
#'   from vegan::vegdist(). See ?vegan::vegdist for more details.
#' @param ordination_method (string, default "NMDS") The ordination method to use
#'   in phyloseq::ordinate(). See ?phyloseq::ordinate for more details.
#' @param ... Additional arguments passed on to phyloseq::plot_ordination()
#' @returns A ggplot2 object
#' @export
#' @author Adrien Taudière
#'
#' @examples
#' library(patchwork)
#' plot_ordination_pq(data_fungi_mini, method = "robust.aitchison", color = "Height") +
#'   plot_ordination_pq(data_fungi_mini, method = "bray", color = "Height")
plot_ordination_pq <- function(
  physeq,
  method = "robust.aitchison",
  ordination_method = "NMDS",
  ...
) {
  verify_pq(physeq)
  physeq <- taxa_as_columns(physeq)
  dist_mat <- vegan::vegdist(unclass(physeq@otu_table), method = method)
  attr(dist_mat, "Labels") <- sample_names(physeq)
  p <- plot_ordination(
    physeq,
    ordination = ordinate(
      physeq,
      method = ordination_method,
      distance = dist_mat
    ),
    ...
  ) +
    labs(
      title = paste(ordination_method, "ordination"),
      subtitle = paste("Using ", method, "distance")
    )
  return(p)
}
################################################################################

################################################################################
# Default y-axis labels for common Hill orders (internal)
.hill_y_lab <- function(q) {
  labs <- c(
    "0" = "Richness (Hill q=0)",
    "1" = "Shannon diversity (Hill q=1)",
    "2" = "Simpson diversity (Hill q=2)"
  )
  lab <- labs[as.character(q)]
  if (is.na(lab)) paste0("Hill index (q=", q, ")") else unname(lab)
}

# Single-panel bar-plot engine used by hill_bar_pq() (internal)
.hill_bar_single <- function(
  data,
  x_name,
  y_name,
  fill_name,
  x_lab,
  y_lab,
  alpha,
  point_size,
  base_size,
  jitter_width,
  bar_width,
  add_letters,
  p_threshold,
  letter_size,
  letters_top_offset,
  y_lab_size,
  x_lab_size,
  show_n_samples,
  palette,
  error_fun,
  error_fun_lab,
  error_bar_alpha,
  point_alpha,
  letters_below_bar
) {
  # --- Kruskal-Wallis test ---
  data[[x_name]] <- as.factor(data[[x_name]])
  kw <- kruskal.test(reformulate(x_name, response = y_name), data = data)
  kw_subtitle <- sprintf(
    "Kruskal-Wallis: X-squared(%d) = %.2f, p = %s",
    kw$parameter,
    kw$statistic,
    format.pval(kw$p.value, digits = 3, eps = 0.001)
  )

  # --- Summary stats ---
  summary_data <- data |>
    dplyr::reframe(
      mean = mean(.data[[y_name]], na.rm = TRUE),
      lower = error_fun(.data[[y_name]])[[1L]],
      upper = error_fun(.data[[y_name]])[[2L]],
      .by = dplyr::all_of(x_name)
    )

  # --- Compact letter display (Tukey HSD after Kruskal-Wallis) ---
  # NA groups are excluded from statistical comparisons; they receive "n.d."
  .grp_chr <- function(x) ifelse(is.na(x), "<NA>", as.character(x))

  tukey_run <- FALSE
  if (add_letters) {
    if (kw$p.value < p_threshold) {
      tukey <- TukeyHSD(aov(
        reformulate(x_name, response = y_name),
        data = data
      ))
      tuk_mat <- tukey[[x_name]]
      pvals <- stats::setNames(tuk_mat[, "p adj"], rownames(tuk_mat))
      if (!is.null(names(pvals)) && length(pvals) > 0L) {
        letters_vec <- multcompView::multcompLetters(pvals)$Letters
        tukey_run <- TRUE
      } else {
        groups <- unique(.grp_chr(data[[x_name]]))
        letters_vec <- stats::setNames(rep("a", length(groups)), groups)
      }
    } else {
      groups <- unique(.grp_chr(data[[x_name]]))
      letters_vec <- stats::setNames(rep("a", length(groups)), groups)
    }

    # Assign "n.d." to any group absent from the letters (e.g. NA groups)
    all_groups <- unique(.grp_chr(data[[x_name]]))
    missing <- setdiff(all_groups, names(letters_vec))
    if (length(missing) > 0L) {
      letters_vec[missing] <- "n.d."
    }

    point_max <- data |>
      dplyr::summarise(
        max_y = max(.data[[y_name]], na.rm = TRUE),
        .by = dplyr::all_of(x_name)
      )
    y_offset <- diff(range(data[[y_name]], na.rm = TRUE)) * letters_top_offset

    summary_data <- summary_data |>
      dplyr::left_join(point_max, by = x_name) |>
      dplyr::mutate(
        letter = letters_vec[.grp_chr(.data[[x_name]])],
        letter_y = if (letters_below_bar) {
          -y_offset
        } else {
          pmax(upper, max_y) + y_offset
        }
      )
  }

  # --- Base plot ---
  p <- ggplot2::ggplot(
    summary_data,
    ggplot2::aes(x = .data[[x_name]], y = mean, fill = .data[[fill_name]])
  ) +
    ggplot2::geom_col(alpha = alpha, width = bar_width) +
    ggplot2::geom_errorbar(
      ggplot2::aes(ymin = lower, ymax = upper),
      width = 0.2,
      linewidth = 0.8
    ) +
    ggplot2::geom_jitter(
      data = data,
      ggplot2::aes(
        x = .data[[x_name]],
        y = .data[[y_name]],
        fill = .data[[fill_name]]
      ),
      shape = 21,
      size = point_size,
      alpha = point_alpha,
      width = jitter_width,
      height = 0,
      inherit.aes = FALSE
    ) +
    ggplot2::geom_errorbar(
      ggplot2::aes(ymin = mean, ymax = upper),
      width = 0.2,
      linewidth = 0.8,
      alpha = error_bar_alpha
    )

  if (add_letters) {
    p <- p +
      ggplot2::geom_text(
        data = summary_data,
        ggplot2::aes(x = .data[[x_name]], y = letter_y, label = letter),
        inherit.aes = FALSE,
        size = letter_size,
        fontface = "bold"
      )
  }

  if (show_n_samples) {
    n_per_group <- data |>
      dplyr::summarise(n = dplyr::n(), .by = dplyr::all_of(x_name))
    x_labels <- stats::setNames(
      paste0(n_per_group[[x_name]], "\n(n=", n_per_group$n, ")"),
      n_per_group[[x_name]]
    )
    p <- p + ggplot2::scale_x_discrete(labels = x_labels)
  }

  p +
    ggplot2::scale_fill_manual(values = palette) +
    ggplot2::labs(
      x = x_lab,
      y = y_lab,
      subtitle = kw_subtitle,
      caption = if (add_letters && tukey_run) {
        paste0(
          "Error bars: ",
          error_fun_lab,
          "\nletters from Tukey HSD pairwise comparisons"
        )
      } else if (add_letters && !tukey_run) {
        paste0(
          "Error bars: ",
          error_fun_lab,
          "; Kruskal-Wallis p \u2265 ",
          p_threshold,
          "\nTukey HSD pairwise comparisons not run (no global significance)"
        )
      } else {
        paste0("Error bars: ", error_fun_lab)
      }
    ) +
    ggplot2::theme_bw(base_size = base_size) +
    ggplot2::theme(
      panel.grid.major = ggplot2::element_blank(),
      panel.grid.minor = ggplot2::element_blank(),
      panel.border = ggplot2::element_blank(),
      axis.line.x = ggplot2::element_line(linewidth = 0.4),
      axis.line.y = ggplot2::element_line(linewidth = 0.4),
      axis.ticks = ggplot2::element_line(linewidth = 0.3),
      strip.background = ggplot2::element_blank(),
      strip.text = ggplot2::element_text(face = "bold"),
      legend.key = ggplot2::element_blank(),
      legend.background = ggplot2::element_blank(),
      legend.position = "none",
      axis.text.y = ggplot2::element_text(
        size = if (is.null(y_lab_size)) base_size else y_lab_size
      ),
      axis.text.x = ggplot2::element_text(
        size = if (is.null(x_lab_size)) base_size else x_lab_size
      ),
      plot.subtitle = ggplot2::element_text(
        size = base_size * 0.8,
        colour = "grey40"
      ),
      plot.caption = ggplot2::element_text(
        size = base_size * 0.7,
        colour = "grey50"
      ),
      plot.margin = ggplot2::margin(5, 5, 5, 5, "pt")
    )
}

################################################################################
#' Bar plot of Hill diversity with SE, jittered points, and Kruskal-Wallis test
#'
#' @description
#'
#' <a href="https://adrientaudiere.github.io/MiscMetabar/articles/Rules.html#lifecycle">
#' <img src="https://img.shields.io/badge/lifecycle-experimental-orange" alt="lifecycle-experimental"></a>
#'
#' For each Hill diversity order in `q`, draws a bar at the group mean (±1 SE)
#' with jittered individual points. A Kruskal-Wallis test is reported in the
#' subtitle; when the global effect is significant, Tukey HSD pairwise
#' comparisons produce compact letter displays above the bars. Multiple values
#' of `q` are assembled into a [patchwork] layout automatically.
#'
#' @inheritParams clean_pq
#' @param x Name (unquoted) of the grouping variable in `sam_data` (x-axis).
#' @param q Numeric vector of Hill diversity orders to plot. The corresponding
#'   `Hill_<q>` columns are computed by [psmelt_samples_pq()].
#'   Default `c(0, 2)`.
#' @param fill Name (unquoted) of the fill aesthetic column. Defaults to `x`.
#' @param x_lab Label for the x-axis. Defaults to the column name of `x`.
#' @param y_labs Named character vector of y-axis labels keyed by `Hill_<q>`
#'   column name (e.g. `c(Hill_0 = "Richness")`). Unspecified orders receive
#'   a default label.
#' @param ncol Number of columns in the patchwork layout when `length(q) > 1`.
#'   Default `NULL` (automatic).
#' @param alpha Transparency of bars. Default `0.6`.
#' @param point_size Size of jittered points. Default `3`.
#' @param base_size Base font size in pts. Default `13`.
#' @param jitter_width Horizontal jitter width. Default `0.15`.
#' @param bar_width Width of bars. Default `0.7`.
#' @param add_letters Logical. Add compact letter display above bars.
#'   Requires the \pkg{multcompView} package. Default `TRUE`.
#' @param p_threshold Significance threshold for the Kruskal-Wallis test.
#'   Below this value, Tukey HSD pairwise comparisons are run and letters
#'   assigned; above it all groups receive `"a"`. Default `0.05`.
#' @param letter_size Size of letter labels in ggplot2 units. Default `5`.
#' @param letters_top_offset Fraction of the y-range added above the highest
#'   point / error-bar to position letters. Default `0.05`.
#' @param y_lab_size Size of y-axis tick labels in pts. Defaults to `base_size`.
#' @param x_lab_size Size of x-axis tick labels in pts. Defaults to `base_size`.
#' @param show_n_samples Logical. If `TRUE`, the number of samples per group is
#'   appended below each x-axis tick label as `(n=X)`. Default `TRUE`.
#' @param palette Character vector of fill colours. Defaults to the Okabe-Ito
#'   palette.
#' @param error_fun Function taking a numeric vector and returning a 2-element
#'   numeric vector `c(lower, upper)` with the actual y-axis bounds of the
#'   error bar (not offsets from the mean). The first element is the lower
#'   bound, the second is the upper bound. This allows asymmetric intervals
#'   such as quantile ranges. Default computes mean ± SE. Example for a 95%
#'   quantile interval: `function(x) quantile(x, c(0.025, 0.975),
#'   na.rm = TRUE)`.
#' @param error_fun_lab Label for the error bar used in the plot caption.
#'   Default `"mean ± SE"`.
#' @param error_bar_alpha Transparency of the secondary top-half error bar
#'   drawn over the jittered points to hint at the upper extent without
#'   obscuring data. Default `0.35`.
#' @param point_alpha Transparency of the jittered data points. Default `0.5`.
#' @param letters_below_bar Logical. When `TRUE`, compact letters are placed
#'   below the x-axis (at `y = -letters_top_offset * y_range`), giving a clean
#'   fixed position independent of data spread. When `FALSE` (default), letters
#'   are placed above whichever is higher: the error bar top or the highest
#'   data point.
#' @param ... Additional arguments passed to [psmelt_samples_pq()] and hence
#'   to [divent::div_hill()] (e.g. `estimator = "naive"`).
#'
#' @return A `ggplot` object when `length(q) == 1`, or a `patchwork` object
#'   when `length(q) > 1`.
#'
#' @export
#' @author Adrien Taudière
#'
#' @examples
#' hill_bar_pq(data_fungi_mini, Height, q = 1)
#' \dontrun{
#' hill_bar_pq(data_fungi_mini, Height, q = 0)
#' hill_bar_pq(data_fungi_mini, Height, q = c(0, 1, 2), ncol = 1)
#' hill_bar_pq(data_fungi_mini, Height,
#'   q = c(0, 2),
#'   y_labs = c(Hill_0 = "Richness", Hill_2 = "Simpson diversity")
#' )
#' hill_bar_pq(data_fungi_mini, Height, add_letters = FALSE)
#' }
#'
#' @seealso [hill_pq()], [psmelt_samples_pq()], [ggbetween_pq()]
hill_bar_pq <- function(
  physeq,
  x,
  q = c(0, 2),
  fill,
  x_lab = NULL,
  y_labs = NULL,
  ncol = NULL,
  alpha = 0.6,
  point_size = 3,
  base_size = 13,
  jitter_width = 0.15,
  bar_width = 0.7,
  add_letters = TRUE,
  p_threshold = 0.05,
  letter_size = 5,
  letters_top_offset = 0.05,
  y_lab_size = NULL,
  x_lab_size = NULL,
  show_n_samples = TRUE,
  palette = c(
    "#E69F00",
    "#56B4E9",
    "#009E73",
    "#F0E442",
    "#0072B2",
    "#D55E00",
    "#CC79A7",
    "#000000"
  ),
  error_fun = function(x) {
    m <- mean(x, na.rm = TRUE)
    se <- sd(x, na.rm = TRUE) / sqrt(sum(!is.na(x)))
    c(lower = m - se, upper = m + se)
  },
  error_fun_lab = "mean \u00b1 SE",
  error_bar_alpha = 0.35,
  point_alpha = 0.5,
  letters_below_bar = FALSE,
  ...
) {
  verify_pq(physeq)

  data <- psmelt_samples_pq(physeq, q = q, ...)

  x_var <- rlang::ensym(x)
  fill_var <- if (missing(fill)) x_var else rlang::ensym(fill)
  x_name <- rlang::as_string(x_var)
  fill_name <- rlang::as_string(fill_var)
  x_lab <- if (is.null(x_lab)) x_name else x_lab

  ys <- paste0("Hill_", q)

  plot_args <- list(
    data = data,
    x_name = x_name,
    fill_name = fill_name,
    x_lab = x_lab,
    alpha = alpha,
    point_size = point_size,
    base_size = base_size,
    jitter_width = jitter_width,
    bar_width = bar_width,
    add_letters = add_letters,
    p_threshold = p_threshold,
    letter_size = letter_size,
    letters_top_offset = letters_top_offset,
    y_lab_size = y_lab_size,
    x_lab_size = x_lab_size,
    show_n_samples = show_n_samples,
    palette = palette,
    error_fun = error_fun,
    error_fun_lab = error_fun_lab,
    error_bar_alpha = error_bar_alpha,
    point_alpha = point_alpha,
    letters_below_bar = letters_below_bar
  )

  if (length(ys) == 1) {
    y_lab <- if (!is.null(y_labs) && ys %in% names(y_labs)) {
      y_labs[[ys]]
    } else {
      .hill_y_lab(q)
    }
    do.call(.hill_bar_single, c(plot_args, list(y_name = ys, y_lab = y_lab)))
  } else {
    plots <- lapply(seq_along(ys), function(i) {
      y_name <- ys[[i]]
      y_lab <- if (!is.null(y_labs) && y_name %in% names(y_labs)) {
        y_labs[[y_name]]
      } else {
        .hill_y_lab(q[[i]])
      }
      do.call(
        .hill_bar_single,
        c(plot_args, list(y_name = y_name, y_lab = y_lab))
      )
    })
    patchwork::wrap_plots(plots, ncol = ncol)
  }
}
################################################################################

Try the MiscMetabar package in your browser

Any scripts or data that you put into this service are public.

MiscMetabar documentation built on June 8, 2026, 5:07 p.m.