Nothing
################################################################################
#' Plot the result of a mt test [phyloseq::mt()]
#'
#' @description
#' `r lifecycle::badge("maturing")`
#'
#' @param mt (required) Result of a mt test from the function [phyloseq::mt()].
#' @param alpha (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{ggplot}}2 plot of result of a mt test
#' @export
#' @seealso [phyloseq::mt()]
plot_mt <-
function(mt = NULL,
alpha = 0.05,
color_tax = "Class",
taxa = "Species") {
d <- mt[mt$plower < alpha, ]
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-class}} object
#'
#' @description
#' `r lifecycle::badge("maturing")`
#' @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{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_even_depth(subset_samples_pq(GP, sample_sums(GP) > 3000))
#' 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{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 (!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.matrix(t(physeq@otu_table))
physeq_accu[physeq_accu > 0] <- 1
accu_all <- vegan::specaccum(physeq_accu)
accu <- list()
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 <- t(physeq@otu_table)
} else {
x <- apply(physeq@otu_table, 1, function(x) {
tapply(x, fact_interm, sum)
})
}
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) {
res <- list()
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 <-
as.factor(unlist(unclass(physeq@sam_data
[match(df$.id, sample_names(physeq)), fact])
[fact]))
} else {
df$fact <- df$.id
}
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
#' `r lifecycle::badge("experimental")`
#'
#' 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()].
#' @param verbose (logical). If TRUE, print additional informations.
#' @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, !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, !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 (rarefy_by_sample_before_merging) {
p_for_dim <- accu_plot(
rarefy_sample_count_by_modality(
rarefy_even_depth(
physeq,
rngseed = 1,
sample.size = sample.size,
verbose = verbose
),
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(
physeq,
rngseed = i,
sample.size =
sample.size,
verbose = verbose
),
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[1: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
#' `r lifecycle::badge("experimental")`
#' @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_even_depth(subset_samples_pq(GP, sample_sums(GP) > 3000))
#' 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 <- list()
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-class}} object
#' @description
#' `r lifecycle::badge("maturing")`
#' @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 (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 taxon 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{
#' if (requireNamespace("pbapply")) {
#' data("GlobalPatterns", package = "phyloseq")
#' GP <- subset_taxa(GlobalPatterns, GlobalPatterns@tax_table[, 1] == "Archaea")
#' circle_pq(GP, "SampleType")
#' 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 <- t(physeq@otu_table)
} else {
otu_tab <- 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")
}
otu_table_tax <-
pbapply::pbapply(otu_tab, 2, function(x) {
pbapply::pbtapply(
x, physeq@tax_table[, taxcol],
function(xx) {
sum(xx, na.rm = TRUE)
}
)
}, cl = nproc)
otu_table_ech <-
pbapply::pbapply(otu_table_tax, 1, function(x) {
pbapply::pbtapply(
x, physeq@sam_data[, taxsamp],
function(xx) {
sum(xx, na.rm = TRUE)
}
)
}, cl = nproc)
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]
# 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
]
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-class}} object
#' @description
#' `r lifecycle::badge("maturing")`
#' @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 taxon 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 <- t(physeq@otu_table)
} else {
otu_tab <- physeq@otu_table
}
if (!add_nb_seq) {
otu_tab[otu_tab > 0] <- 1
mat_interm <- matrix()
mat <- matrix(ncol = 3)
colnames(mat) <- c("Var1", "Var2", "value")
for (i in 1:(length(taxa) - 1)) {
res_interm <-
table(physeq@tax_table[, taxa[i]], physeq@tax_table[, taxa[i + 1]])
mat_interm <- reshape2::melt(res_interm)
mat_interm <- mat_interm[mat_interm[, 3] > 0, ]
mat <- rbind(mat, mat_interm)
}
} else if (add_nb_seq) {
mat_interm <- matrix()
mat <- matrix(ncol = 3)
colnames(mat) <- c("Var1", "Var2", "value")
tax_table_interm <-
physeq@tax_table[rep(seq(1, ntaxa(physeq)),
times = taxa_sums(physeq)
)]
for (i in 1:(length(taxa) - 1)) {
res_interm <-
table(tax_table_interm[, taxa[i]], tax_table_interm[, taxa[i + 1]])
mat_interm <- reshape2::melt(res_interm)
mat_interm <- mat_interm[mat_interm[, 3] > 0, ]
mat <- rbind(mat, mat_interm)
}
}
if (!is.null(fact)) {
net_matrix2links <- function(m = NULL) {
res <- matrix(ncol = 3)
for (i in seq_len(dim(m)[1])) {
for (j in seq_len(dim(m)[2])) {
if (m[i, j] > 0) {
res <- rbind(res, c(rownames(m)[i], colnames(m)[j], m[i, j]))
}
}
}
return(res)
}
mat_interm <-
apply(otu_tab, 1, function(x) {
tapply(
x, physeq@sam_data[, fact],
sum
)
})
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)
})
}
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_len(length(symbol2sub))) {
mat <- apply(mat, 2, function(x) {
gsub(symbol2sub[i], "", x)
})
}
tax_sank <- list()
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_len(length(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-class}} object
#' @description
#' `r lifecycle::badge("maturing")`
#' @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{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 (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)
)
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 <- list()
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 = .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-class}} object using
#' `ggVennDiagram::ggVennDiagram` function
#' @description
#' `r lifecycle::badge("maturing")`
#'
#' 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 ... Other arguments for the `ggVennDiagram::ggVennDiagram` function
#' for ex. `category.names`.
#' @return A \code{\link{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, fact = "Height")
#' }
#' \donttest{
#' if (requireNamespace("ggVennDiagram")) {
#' ggvenn_pq(data_fungi, fact = "Height") +
#' ggplot2::scale_fill_distiller(palette = "BuPu", direction = 1)
#' pl <- ggvenn_pq(data_fungi, fact = "Height", split_by = "Time")
#' for (i in 1:length(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, data_fungi@sam_data$Tree_name == "A10-005" |
#' data_fungi@sam_data$Height %in% c("Low", "High"))
#' ggvenn_pq(data_fungi2, fact = "Height")
#'
#' ggvenn_pq(data_fungi, fact = "Height", add_nb_seq = TRUE, set_size = 4)
#' ggvenn_pq(data_fungi, fact = "Height", rarefy_before_merging = TRUE)
#' ggvenn_pq(data_fungi, fact = "Height", rarefy_after_merging = TRUE) +
#' scale_x_continuous(expand = expansion(mult = 0.5))
#' }
#' }
#' @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,
...) {
if (!is.factor(physeq@sam_data[[fact]])) {
physeq@sam_data[[fact]] <- as.factor(physeq@sam_data[[fact]])
}
physeq <- taxa_as_columns(physeq)
if (rarefy_before_merging) {
physeq <- rarefy_even_depth(physeq)
physeq <- clean_pq(physeq)
}
nb_samples <- table(physeq@sam_data[[fact]])
if (rarefy_after_merging) {
physeq <- merge_samples2(physeq, fact)
physeq <- rarefy_even_depth(physeq)
physeq <- clean_pq(physeq)
}
res <- list()
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, ]
sample_data(newphyseq) <- sample_data(new_DF)
if (is.null(taxonomic_rank)) {
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 (max(nb_seq) / min(nb_seq) > 2) {
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 {
p <- list()
modalities <-
as.factor(unlist(unclass(physeq@sam_data[[split_by]])))
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)
}
}
return(p)
}
################################################################################
################################################################################
#' Multiple plot function
#' @description
#' `r lifecycle::badge("stable")`
#'
# 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
#' `r lifecycle::badge("experimental")`
#' 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 hill_scales (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).
#' @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): 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 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
#' @return Either an unique ggplot2 object (if one_plot is TRUE) or
#' a list of 4 ggplot2 plot:
#' - 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
#'
#' p <- hill_pq(data_fungi_mini, "Height", hill_scales = 1:2)
#' p_h1 <- p[[1]] + theme(legend.position = "none")
#' p_h2 <- p[[2]] + theme(legend.position = "none")
#' multiplot(plotlist = list(p_h1, p_h2, p[[3]]), cols = 4)
#' \donttest{
#' if (requireNamespace("multcompView")) {
#' p2 <- hill_pq(data_fungi, "Time",
#' correction_for_sample_size = FALSE,
#' letters = TRUE, add_points = TRUE, plot_with_tuckey = FALSE
#' )
#' if (requireNamespace("patchwork")) {
#' patchwork::wrap_plots(p2, guides = "collect")
#' }
#' # Artificially modify data_fungi to force alpha-diversity effect
#' data_fungi_modif <- clean_pq(subset_samples_pq(data_fungi, !is.na(data_fungi@sam_data$Height)))
#' data_fungi_modif@otu_table[data_fungi_modif@sam_data$Height == "High", ] <-
#' data_fungi_modif@otu_table[data_fungi_modif@sam_data$Height == "High", ] +
#' sample(c(rep(0, ntaxa(data_fungi_modif) / 2), rep(100, ntaxa(data_fungi_modif) / 2)))
#' p3 <- hill_pq(data_fungi_modif, "Height", letters = TRUE)
#' p3[[1]]
#' }
#' }
#' @seealso [psmelt_samples_pq()] and [ggbetween_pq()]
hill_pq <- function(physeq,
fact = NULL,
variable = NULL,
hill_scales = c(0, 1, 2),
color_fac = NA,
letters = FALSE,
add_points = FALSE,
add_info = TRUE,
one_plot = FALSE,
plot_with_tuckey = TRUE,
correction_for_sample_size = TRUE,
na_remove = TRUE) {
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 {
if (!is.null(fact)) {
variable <- fact
} else {
stop("You must set the parameter fact.")
}
}
var <- sym(variable)
if (is.na(color_fac)) {
color_fac <- sym(variable)
} 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]])
otu_hill <-
vegan::renyi(t(physeq)@otu_table,
scales = hill_scales,
hill = TRUE
)
colnames(otu_hill) <- paste0("Hill_", hill_scales)
df_hill <- data.frame(otu_hill, physeq@sam_data)
df_hill[, seq_along(hill_scales)] <-
apply(df_hill[, seq_along(hill_scales)], 2, as.numeric)
p_var <-
hill_tuckey_pq(
physeq,
modality = variable,
hill_scales = hill_scales,
correction_for_sample_size = correction_for_sample_size
)
p_list <- list()
for (i in seq_along(hill_scales)) {
p_list[[i]] <-
ggplot(df_hill, aes(group = !!var, .data[[paste0("Hill_", hill_scales[[i]])]])) +
geom_boxplot(outlier.size = 2, aes(colour = as.factor(!!color_fac), y = !!var)) +
labs(x = paste0("Hill_", hill_scales[[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]])),
sep = "' : ",
table(physeq@sam_data[[variable]]),
collapse = " - '"
)
)
p_list[[i]] <- p_list[[i]] + labs(subtitle = subtitle_plot)
}
if (letters) {
data_h <-
p_var$data[grep(paste0("Hill_", hill_scales[[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_", hill_scales[[i]])]]) + 1) %>%
inner_join(dt)
p_list[[i]] <- p_list[[i]] +
geom_label(
data = data_letters,
aes(
x = pos_letters,
label = Letters
),
y = ggplot_build(p_list[[i]])$data[[1]]$y,
size = 4,
stat = "unique",
parse = TRUE
)
}
}
res <- p_list
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
#' `r lifecycle::badge("experimental")`
#'
#' 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 ... Other arguments passed on to [ggstatsplot::ggbetweenstats()] function.
#' @return Either an unique ggplot2 object (if one_plot is TRUE) or
#' a list of 3 ggplot2 plot:
#' - 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{
#' if (requireNamespace("ggstatsplot")) {
#' p <- ggbetween_pq(data_fungi, fact = "Time", p.adjust.method = "BH")
#' p[[1]]
#' 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.
ggbetween_pq <-
function(physeq,
fact,
one_plot = FALSE,
rarefy_by_sample = FALSE,
...) {
verify_pq(physeq)
physeq <- taxa_as_columns(physeq)
if (rarefy_by_sample) {
physeq <- clean_pq(rarefy_even_depth(physeq))
}
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"
)
)
}
df <- cbind(
"nb_asv" = sample_sums(physeq@otu_table),
physeq@sam_data,
"hill_0" = vegan::renyi(physeq@otu_table, scales = 0, hill = TRUE),
"hill_1" = vegan::renyi(physeq@otu_table, scales = 1, hill = TRUE),
"hill_2" = vegan::renyi(physeq@otu_table, scales = 2, hill = TRUE)
)
fact <- sym(fact)
p0 <- ggstatsplot::ggbetweenstats(df, !!fact, hill_0, ...)
p1 <- ggstatsplot::ggbetweenstats(df, !!fact, hill_1, ...)
p2 <- ggstatsplot::ggbetweenstats(df, !!fact, hill_2, ...)
res <- list(
"plot_Hill_0" = p0,
"plot_Hill_1" = p1,
"plot_Hill_2" = p2
)
if (one_plot) {
requireNamespace("patchwork", quietly = TRUE)
res <- res[[1]] + res[[2]] + res[[3]]
}
return(res)
}
################################################################################
#' Summarize a \code{\link{phyloseq-class}} object using a plot.
#' @description
#' `r lifecycle::badge("maturing")`
#' @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?
#' @examples
#'
#' summary_plot_pq(data_fungi)
#' summary_plot_pq(data_fungi, add_info = FALSE) + scale_fill_viridis_d()
#' @return A ggplot2 object
#' @export
summary_plot_pq <- function(physeq,
add_info = TRUE,
min_seq_samples = 500,
clean_pq = TRUE) {
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)
)
)
)
)
)
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
) +
geom_text(
data = d,
aes(
x = x1 + (x2 - x1) / 2,
y = y1 + (y2 - y1) / 5,
label = names
),
size = 6
) +
annotate(
geom = "text",
x = 0.65,
y = 3,
label = "Taxa",
size = 6,
color = "#aa4c26",
angle = 90,
fontface = 2
) +
annotate(
geom = "text",
x = 0.85,
y = 3,
label = "(OTUs, ASVs, ...)",
size = 5,
color = "#aa4c26",
angle = 90
) +
annotate(
geom = "text",
x = 2,
y = 0.65,
label = "Samples",
size = 6,
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))),
" ASV)",
"\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,
hjust = 0
)
}
return(p)
}
################################################################################
################################################################################
#' rotl wrapper for phyloseq data
#' @description
#' `r lifecycle::badge("experimental")`
#' Make a phylogenetic tree using the ASV names of a physeq object and the
#' Open Tree of Life tree.
#'
#' @inheritParams clean_pq
#' @param species_colnames (default: "Genus_species"): the name of the column
#' where the species binominal name is stored in `@tax_table` slot.
#' Can also be a vector of two columns names e.g. `c("Genus", "Species")`
#' @param context_name : can bue used to select only a part of the Open Tree
#' of Life. See `?rotl::tnrs_contexts()` for available values
#' @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.
#' @examplesIf tolower(Sys.info()[["sysname"]]) != "windows"
#' \donttest{
#' if (requireNamespace("rotl")) {
#' tr <- rotl_pq(data_fungi_mini, species_colnames = "Genus_species")
#' plot(tr)
#'
#' tr_Asco <- rotl_pq(data_fungi, species_colnames = "Genus_species", context_name = "Ascomycetes")
#' plot(tr_Asco)
#' }
#' }
rotl_pq <- function(physeq,
species_colnames = "Genus_species",
context_name = "All life") {
if (length(species_colnames) == 2) {
physeq@tax_table <- tax_table(cbind(
physeq@tax_table,
"Genus_species" = paste(physeq@tax_table[, species_colnames[1]],
physeq@tax_table[, species_colnames[2]],
sep = "_"
)
))
species_colnames <- "Genus_species"
}
taxa_names_rotl <- physeq@tax_table[, species_colnames]
taxa_names_rotl <- taxa_names_rotl[!grepl("NA", taxa_names_rotl)]
taxa_names_rotl <- c(unclass(gsub("_", " ", taxa_names_rotl)))
resolved_names <- 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 <-
rotl::tnrs_match_names(clean_taxa_names_rotl, context_name = context_name)
tr <-
rotl::tol_induced_subtree(ott_ids = rotl::ott_id(resolved_names2))
return(tr)
}
################################################################################
################################################################################
#' Heat tree from `metacoder` package using `tax_table` slot
#' @description
#' `r lifecycle::badge("maturing")`
#'
#' 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
#' `r lifecycle::badge("maturing")`
#' @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 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 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")
#' @export
#' @author Adrien Taudière
#'
biplot_pq <- function(physeq,
fact = NULL,
merge_sample_by = NULL,
rarefy_after_merging = FALSE,
inverse_side = FALSE,
left_name = NULL,
left_name_col = "#4B3E1E",
left_fill = "#4B3E1E",
left_col = "#f3f2d9",
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,
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]
}
physeq <- merge_samples2(physeq, merge_sample_by)
physeq <- clean_pq(physeq)
}
if (nsamples(physeq) != 2) {
stop(
"biplot_pq needs only two samples in the
physeq object or a valid merge_sample_by parameter"
)
}
if (rarefy_after_merging) {
physeq <- clean_pq(rarefy_even_depth(physeq))
}
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 (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, ]
# mdf <- dplyr::rename(mdf, Abundance = Abundance)
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
nudge_y <- mean(mdf$Abundance) * nudge_y
}
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), "%")
p <- mdf %>%
ggplot(
aes(
x = stats::reorder(OTU, Abundance),
y = Ab,
fill = modality,
names = OTU,
Ab = Abundance,
Proportion = Proportion,
Family = Family,
Genus = Genus,
Species = Species
),
...
) +
geom_bar(stat = "identity", width = .6) +
annotate(
"rect",
xmin = "Samples",
xmax = "Samples",
ymin = -max(mdf$Ab),
ymax = max(mdf$Ab)
) +
annotate(
geom = "text",
label = right_name,
x = "Samples",
y = ifelse(is.na(y_names), max(mdf$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(mdf$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$OTU, sum)
)), "Samples")) +
ylim(min(mdf$Ab), max(mdf$Ab) * 1.1)
if (geom_label) {
p <- p +
geom_label(
aes(
label = Abundance,
color = modality,
fill = modality,
alpha = 0.5,
y = ifelse(Ab > 0, Ab + nudge_y[2], Ab + nudge_y[1])
),
size = text_size
)
} else {
p <- p +
geom_text(
aes(
label = Abundance,
color = modality,
y = ifelse(Ab > 0, Ab + nudge_y[2], Ab + nudge_y[1])
),
size = text_size
)
}
p <- p + coord_flip() +
theme_minimal() +
theme(plot.title = element_text(hjust = .5), axis.ticks = element_blank()) +
scale_fill_manual(values = c(left_fill, right_fill)) +
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 (plotly_version) {
p <- plotly::ggplotly(
p,
tooltip = c("OTU", "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
#' `r lifecycle::badge("experimental")`
#'
#' 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 pairs 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 pass 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, taxa_sums(data_fungi) > 10000)
#' 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 <- list()
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 <- list()
for (c in seq_along(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
#' `r lifecycle::badge("experimental")`
#'
#' 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_asv", 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 graphic
#' @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_asv", "both")) {
mdf <-
psmelt(as_binary_otu_table(physeq))
mdf <- mdf %>% mutate(percent = Abundance / sum(Abundance))
p_asv <-
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_asv", "both")) {
p_asv <- p_asv +
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_asv") {
return(p_asv)
} else if (type == "both") {
return(list(p_seq, p_asv))
}
}
################################################################################
################################################################################
#' Plot taxonomic distribution across 3 taxonomic levels and optionally
#' one sample factor
#'
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' 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 graphic
#' @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 ... Other 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)
#' }
tsne_pq <-
function(physeq,
method = "bray",
dims = 2,
theta = 0.0,
perplexity = 30,
...) {
physeq <- taxa_as_rows(physeq)
res_tsne <-
Rtsne::Rtsne(
vegan::vegdist(t(physeq@otu_table), method = method),
dims = dims,
theta = theta,
perplexity = perplexity,
is_distance = TRUE,
...
)
return(res_tsne)
}
################################################################################
################################################################################
#' Plot a tsne low dimensional representation of a phyloseq object
#'
#' `r lifecycle::badge("experimental")`
#'
#' 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 ... Other arguments passed on to `Rtsne::Rtsne()`
#'
#' @return
#' A ggplot object
#'
#' @export
#' @author Adrien Taudière
#'
#' @examplesIf tolower(Sys.info()[["sysname"]]) != "windows"
#' data(data_fungi)
#' if (requireNamespace("Rtsne")) {
#' plot_tsne_pq(data_fungi, fact = "Height", perplexity = 15)
#' }
#' \donttest{
#' if (requireNamespace("Rtsne")) {
#' plot_tsne_pq(data_fungi, fact = "Time") + geom_label(aes(label = Sample_id, fill = Time))
#' plot_tsne_pq(data_fungi, 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]]))
}
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
#' `r lifecycle::badge("experimental")`
#' @inheritParams clean_pq
#' @param clean_pq (logical): Does the phyloseq
#' object is cleaned using the [clean_pq()] function?
#' @param ... Other 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
#' )
#' 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
#' `r lifecycle::badge("experimental")`
#'
#' @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
#' \donttest{
#' 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
#' `r lifecycle::badge("experimental")`
#'
#' 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 ... Other arguments passed on to the [ComplexUpset::upset()]
#'
#' @return A \code{\link{ggplot}}2 plot
#' @export
#' @author Adrien Taudière
#'
#' @seealso [ggvenn_pq()]
#' @examples
#' if (requireNamespace("ComplexUpset")) {
#' upset_pq(data_fungi_mini,
#' fact = "Height", width_ratio = 0.2,
#' taxa_fill = "Class"
#' )
#' }
#' \donttest{
#' if (requireNamespace("ComplexUpset")) {
#' 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,
...) {
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)
if (rarefy_after_merging) {
physeq <- clean_pq(rarefy_even_depth(physeq))
}
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("ASV" = (
ggplot(mapping = aes(fill = .data[[taxa_fill]])) +
geom_bar() +
ylab("ASV 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
#' `r lifecycle::badge("experimental")`
#'
#' @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 ... Other arguments passed on to the [ComplexUpset::upset_test()]
#'
#' @return A \code{\link{ggplot}}2 plot
#' @export
#' @author Adrien Taudière
#'
#' @seealso [upset_pq()]
#' @examples
#' data(data_fungi)
#' if (requireNamespace("ComplexUpset")) {
#' upset_test_pq(data_fungi, "Height", var_to_test = c("OTU", "Class", "Guild"))
#' upset_test_pq(data_fungi, "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
#' `r lifecycle::badge("experimental")`
#' 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 ... Other 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)
} 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)
} 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
#' `r lifecycle::badge("experimental")`
#'
#' @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)
#' @return A \code{\link{ggplot}}2 plot with bar representing the number of sequence en each
#' taxonomic groups
#' @export
#'
#' @examples
#'
#' data_fungi_ab <- subset_taxa_pq(data_fungi, taxa_sums(data_fungi) > 10000)
#' tax_bar_pq(data_fungi_ab) + theme(legend.position = "none")
#' \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")
#' }
#' @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) {
if (!nb_seq) {
physeq <- as_binary_otu_table(physeq)
}
psm <- psmelt(physeq)
if (percent_bar) {
ggplot(psm) +
geom_bar(aes(x = .data[[fact]], fill = .data[[taxa]], y = Abundance),
stat = "identity",
position = "fill"
)
} else {
ggplot(psm) +
geom_bar(aes(x = .data[[fact]], fill = .data[[taxa]], y = Abundance),
stat = "identity"
)
}
}
################################################################################
################################################################################
################################################################################
#' Ridge plot of a phyloseq object
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' @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 ... Other params passed on to [ggridges::geom_density_ridges()]
#'
#' @return A \code{\link{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, "Sample_names", log10trans = TRUE)
#'
#' 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,
...) {
psm <- psmelt(physeq)
psm <- psm %>% filter(Abundance > 0)
if (log10trans) {
psm$Abundance <- log10(psm$Abundance)
}
if (nb_seq) {
p <- ggplot(psm, aes(y = factor(.data[[fact]]), x = Abundance)) +
ggridges::geom_density_ridges(aes(fill = Class), ...) +
xlim(c(0, NA))
} else {
psm_asv <-
psm %>%
group_by(.data[[fact]], OTU, Class) %>%
summarise("count" = n())
p <-
ggplot(psm_asv, aes(y = factor(.data[[fact]]), x = count)) +
ggridges::geom_density_ridges(
aes(fill = Class),
...
) +
xlim(c(0, NA))
}
return(p)
}
################################################################################
################################################################################
################################################################################
#' Plot treemap of 2 taxonomic levels
#'
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' 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
#' transformed.
#' @param plot_legend (logical, default FALSE) If TRUE, plot che
#' legend of color for lvl 1
#' @param ... Other arguments passed on to [treemapify::geom_treemap()] function.
#'
#' @return A ggplot2 graphic
#' @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 <- function(physeq,
lvl1,
lvl2,
nb_seq = TRUE,
log10trans = TRUE,
plot_legend = FALSE,
...) {
if (!nb_seq) {
physeq <- as_binary_otu_table(physeq)
}
psm <- psmelt(physeq) %>%
filter(!is.na(.data[[lvl2]])) %>%
filter(!is.na(.data[[lvl1]]))
psm2 <- psm %>%
group_by(.data[[lvl2]]) %>%
reframe(Abundance = sum(Abundance), LVL1 = unique(.data[[lvl1]]))
if (log10trans) {
psm2$Abundance <- log10(psm2$Abundance)
}
p <-
ggplot(
psm2,
aes(
area = Abundance,
fill = LVL1,
label = .data[[lvl2]],
subgroup = LVL1
)
) +
treemapify::geom_treemap(...) +
treemapify::geom_treemap_subgroup_border(colour = "white", size = 4) +
treemapify::geom_treemap_text(
colour = "white",
place = "centre",
size = 15,
grow = TRUE
)
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
#' `r lifecycle::badge("experimental")`
#'
#' @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, !is.na(Time) & !is.na(Height))
#' res_var_9 <- var_par_rarperm_pq(
#' data_fungi_woNA,
#' list_component = list(
#' "Time" = c("Time"),
#' "Size" = c("Height", "Diameter")
#' ),
#' nperm = 9,
#' dbrda_computation = TRUE
#' )
#' 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
#' )
#' res_var0 <- var_par_pq(data_fungi_woNA,
#' list_component = list(
#' "Time" = c("Time"),
#' "Size" = c("Height", "Diameter")
#' ),
#' dbrda_computation = TRUE
#' )
#' plot_var_part_pq(res_var0, digits_quantile = 2, show_dbrda_signif = TRUE)
#' plot_var_part_pq(res_var_9,
#' digits_quantile = 2, show_quantiles = TRUE,
#' 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 <-
c(1:length(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 <-
c(1:length(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 (any(is.na(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
#' `r lifecycle::badge("experimental")`
#'
#' 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 numberk
#' @param hill_scales (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).
#' @param rarefy_by_sample (logical, default FALSE) If TRUE, rarefy
#' samples using [phyloseq::rarefy_even_depth()] function.
#' @param one_plot (logical, default FALSE) If TRUE, return a unique
#' plot with the three plot inside using the patchwork package.
#' @param ... Other arguments passed on to [ggstatsplot::ggscatterstats()]
#' function.
#'
#' @return Either an unique ggplot2 object (if one_plot is TRUE) or
#' a list of ggplot2 plot for each hill_scales.
#' @export
#' @author Adrien Taudière
#'
#' @examples
#' if (requireNamespace("ggstatsplot")) {
#' ggscatt_pq(data_fungi_mini, "Time", type = "non-parametric")
#' ggscatt_pq(data_fungi_mini, "Time", hill_scales = 1:4, type = "parametric")
#' ggscatt_pq(data_fungi_mini, "Sample_id",
#' hill_scales = c(0, 0.5),
#' 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,
hill_scales = c(0, 1, 2),
rarefy_by_sample = FALSE,
one_plot = TRUE,
...) {
verify_pq(physeq)
physeq <- clean_pq(physeq, force_taxa_as_columns = TRUE)
if (rarefy_by_sample) {
physeq <- clean_pq(rarefy_even_depth(physeq))
}
p_list <- list()
psm_res <- psmelt_samples_pq(physeq, hill_scales = hill_scales)
for (i in seq_along(hill_scales)) {
p_list[[i]] <-
ggstatsplot::ggscatterstats(
psm_res, !!paste0("Hill_", hill_scales[[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
#' `r lifecycle::badge("experimental")`
#'
#' 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_asv")
#' @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 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_asv", 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 ... Other 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")) {
#' ggaluv_pq(data_fungi_mini, type = "nb_asv")
#'
#' ggaluv_pq(data_fungi_mini, wrap_factor = "Height", by_sample = TRUE, type = "nb_asv") +
#' facet_wrap("Height")
#'
#' ggaluv_pq(data_fungi_mini,
#' width = 0.9, min.size = 10,
#' type = "nb_asv", 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.
#' @seealso [sankey_pq()]
ggaluv_pq <- function(physeq,
taxa_ranks = c("Phylum", "Class", "Order", "Family"),
wrap_factor = NULL,
by_sample = FALSE,
rarefy_by_sample = FALSE,
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) {
physeq <- rarefy_even_depth(physeq)
}
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_asv") {
physeq <- as_binary_otu_table(physeq)
} else if (type != "nb_seq") {
stop("Type must be eiter nb_seq or nb_asv")
}
psm_samp <-
psmelt_samples_pq(
physeq,
taxa_ranks = taxa_ranks,
hill_scales = 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)
)
) +
ggalluvial::geom_flow(...) +
ggalluvial::geom_stratum() +
theme_minimal() +
theme(legend.position = "none")
if (use_ggfittext) {
p <- p +
ggfittext::geom_fit_text(
aes(label = stratum),
stat = "stratum",
width = width,
min.size = min.size
)
} else if (use_geom_label) {
p <- p +
geom_label(
aes(label = stratum),
stat = "stratum",
size = size_lab
)
}
if (!is.null(wrap_factor)) {
p <- p + facet_wrap(wrap_factor)
}
return(p)
}
################################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.