Nothing
################################################################################
#' @title Performs graph-based permutation tests on phyloseq object
#' @description
#' `r lifecycle::badge("maturing")`
#'
#' A wrapper of [phyloseqGraphTest::graph_perm_test()] for quick plot with
#' important statistics
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor to cluster samples by modalities.
#' Need to be in \code{physeq@sam_data}. This should be a factor
#' with two or more levels.
#' @param merge_sample_by a vector to determine
#' which samples to merge using [merge_samples2()] function.
#' Need to be in \code{physeq@sam_data}
#' @param nperm (int) The number of permutations to perform.
#' @param return_plot (logical) Do we return only the result
#' of the test or do we plot the result?
#' @param title The title of the Graph.
#' @param na_remove (logical, default FALSE) If set to TRUE, remove samples with
#' NA in the variables set in formula.
#' @param ... Other params for be passed on to
#' [phyloseqGraphTest::graph_perm_test()] function
#'
#' @examples
#' \donttest{
#' if (requireNamespace("phyloseqGraphTest")) {
#' data(enterotype)
#' graph_test_pq(enterotype, fact = "SeqTech")
#' graph_test_pq(enterotype, fact = "Enterotype", na_remove = TRUE)
#' }
#' }
#' @author Adrien Taudière
#'
#' @return A \code{\link{ggplot}}2 plot with a subtitle indicating the pvalue
#' and the number of permutations
#' @details
#' This function is mainly a wrapper of the work of others.
#' Please cite `phyloseqGraphTest` package.
#' @export
graph_test_pq <- function(physeq,
fact,
merge_sample_by = NULL,
nperm = 999,
return_plot = TRUE,
title = "Graph Test",
na_remove = FALSE,
...) {
verify_pq(physeq)
if (!is.null(merge_sample_by)) {
physeq <- merge_samples2(physeq, merge_sample_by)
physeq <- clean_pq(physeq)
}
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
}
res_graph_test <- phyloseqGraphTest::graph_perm_test(physeq,
sampletype = fact,
nperm = nperm,
...
)
if (!return_plot) {
return(res_graph_test)
} else {
p <- phyloseqGraphTest::plot_test_network(res_graph_test) +
labs(
title = title,
subtitle = paste(
"pvalue = ",
res_graph_test$pval,
"(",
length(res_graph_test$perm),
" permutations)"
)
)
return(p)
}
}
################################################################################
################################################################################
#' @title Permanova on a phyloseq object
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' A wrapper for the [vegan::adonis2()] function in the case of `physeq` object.
#' @inheritParams clean_pq
#' @param formula (required) the right part of a formula for [vegan::adonis2()].
#' Variables must be present in the `physeq@sam_data` slot.
#' @param dist_method (default "bray") the distance used. See
#' [phyloseq::distance()] for all available distances or run
#' [phyloseq::distanceMethodList()].
#' For aitchison and robust.aitchison distance, [vegan::vegdist()]
#' function is directly used.
#' @param merge_sample_by a vector to determine
#' which samples to merge using the [merge_samples2()]
#' function. Need to be in `physeq@sam_data`
#' @param na_remove (logical, default FALSE) If set to TRUE, remove samples with
#' NA in the variables set in formula.
#' @param correction_for_sample_size (logical, default FALSE) If set to TRUE,
#' the sample size (number of sequences by samples) is add to formula in
#' the form `y~Library_Size + Biological_Effect` following recommendation of
#' [Weiss et al. 2017](https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-017-0237-y).
#' `correction_for_sample_size` overcome `rarefy_nb_seqs` if both are TRUE.
#' @param rarefy_nb_seqs (logical, default FALSE) Rarefy each sample
#' (before merging if merge_sample_by is set) using
#' `phyloseq::rarefy_even_depth()`.
#' if `correction_for_sample_size` is TRUE, rarefy_nb_seqs will have no
#' effect.
#' @param verbose (logical, default TRUE) If TRUE, prompt some messages.
#' @param ... Other arguments passed on to [vegan::adonis2()] function.
#' @return The function returns an anova.cca result object with a
#' new column for partial R^2. See help of [vegan::adonis2()] for
#' more information.
#' @examples
#' data(enterotype)
#' \donttest{
#' adonis_pq(enterotype, "SeqTech*Enterotype", na_remove = TRUE)
#' adonis_pq(enterotype, "SeqTech", dist_method = "jaccard")
#' adonis_pq(enterotype, "SeqTech", dist_method = "robust.aitchison")
#' }
#' @export
#' @author Adrien Taudière
#' @details
#' This function is mainly a wrapper of the work of others.
#' Please make a reference to `vegan::adonis2()` if you
#' use this function.
adonis_pq <- function(physeq,
formula,
dist_method = "bray",
merge_sample_by = NULL,
na_remove = FALSE,
correction_for_sample_size = FALSE,
rarefy_nb_seqs = FALSE,
verbose = TRUE,
...) {
physeq <- taxa_as_columns(physeq)
if (dist_method %in% c("aitchison", "robust.aitchison")) {
phy_dist <-
paste0(
'vegan::vegdist(as.matrix(physeq@otu_table), method="',
dist_method,
'")'
)
} else {
phy_dist <-
paste0('phyloseq::distance(physeq, method="', dist_method, '")')
}
.formula <- stats::reformulate(formula, response = phy_dist)
termf <- stats::terms(.formula)
term_lab <- attr(termf, "term.labels")[attr(termf, "order") == 1]
verify_pq(physeq)
if (na_remove) {
new_physeq <- physeq
for (tl in term_lab) {
if (verbose) {
message(tl)
}
new_physeq <-
subset_samples_pq(new_physeq, !is.na(physeq@sam_data[[tl]]))
}
if (nsamples(physeq) - nsamples(new_physeq) > 0) {
message(
paste0(
nsamples(physeq) - nsamples(new_physeq),
" were discarded due to NA in variables present in formula."
)
)
}
physeq <- new_physeq
}
if (!is.null(merge_sample_by)) {
physeq <- merge_samples2(physeq, merge_sample_by)
physeq <- clean_pq(physeq)
}
if (correction_for_sample_size) {
formula <- paste0("sample_size+", formula)
.formula <- stats::reformulate(formula, response = phy_dist)
} else if (rarefy_nb_seqs) {
physeq <- rarefy_even_depth(physeq)
physeq <- clean_pq(physeq)
}
verify_pq(physeq)
metadata <- as(sample_data(physeq), "data.frame")
if (correction_for_sample_size) {
metadata$sample_size <- sample_sums(physeq)
}
res_ado <- vegan::adonis2(.formula, data = metadata, ...)
return(res_ado)
}
################################################################################
################################################################################
#' Permanova (adonis) on permutations of rarefaction even depth
#'
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' @inheritParams adonis_pq
#' @param nperm (int) The number of permutations to perform.
#' @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 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 ... Other params for be passed on to [adonis_pq()] function
#'
#' @return A list of three dataframe representing the mean, the minimum quantile
#' and the maximum quantile value for adonis results. See [adonis_pq()].
#' @export
#' @author Adrien Taudière
#' @seealso [adonis_pq()]
#' @examples
#' if (requireNamespace("vegan")) {
#' data_fungi_woNA <-
#' subset_samples(data_fungi, !is.na(Time) & !is.na(Height))
#' adonis_rarperm_pq(data_fungi_woNA, "Time*Height", na_remove = TRUE, nperm = 3)
#' }
adonis_rarperm_pq <- function(physeq,
formula,
dist_method = "bray",
merge_sample_by = NULL,
na_remove = FALSE,
rarefy_nb_seqs = FALSE,
verbose = TRUE,
nperm = 99,
progress_bar = TRUE,
quantile_prob = 0.975,
sample.size = min(sample_sums(physeq)),
...) {
res_perm <- list()
if (progress_bar) {
pb <- txtProgressBar(
min = 0,
max = nperm,
style = 3,
width = 50,
char = "="
)
}
for (i in 1:nperm) {
res_perm[[i]] <-
adonis_pq(
rarefy_even_depth(
physeq,
rngseed = i,
sample.size = sample.size,
verbose = verbose
),
formula,
dist_method = dist_method,
merge_sample_by = merge_sample_by,
na_remove = na_remove,
correction_for_sample_size = FALSE,
rarefy_nb_seqs = rarefy_nb_seqs,
verbose = verbose,
sample.size = sample.size,
...
)
if (progress_bar) {
setTxtProgressBar(pb, i)
}
}
res_adonis <- list()
res_adonis[["mean"]] <-
apply(array(unlist(res_perm), c(dim(
as.data.frame(res_perm[[1]])
), nperm)), c(1, 2), mean)
colnames(res_adonis[["mean"]]) <- colnames(res_perm[[1]])
rownames(res_adonis[["mean"]]) <- rownames(res_perm[[1]])
res_adonis[["quantile_min"]] <-
apply(
array(unlist(res_perm), c(dim(
as.data.frame(res_perm[[1]])
), nperm)),
c(1, 2),
quantile,
na.rm = TRUE,
probs = 1 - quantile_prob
)
colnames(res_adonis[["quantile_min"]]) <- colnames(res_perm[[1]])
rownames(res_adonis[["quantile_min"]]) <- rownames(res_perm[[1]])
res_adonis[["quantile_max"]] <-
apply(
array(unlist(res_perm), c(dim(
as.data.frame(res_perm[[1]])
), nperm)),
c(1, 2),
quantile,
na.rm = TRUE,
probs = quantile_prob
)
colnames(res_adonis[["quantile_max"]]) <- colnames(res_perm[[1]])
rownames(res_adonis[["quantile_max"]]) <- rownames(res_perm[[1]])
return(res_adonis)
}
################################################################################
################################################################################
#' @title Compute and test local contributions to beta diversity (LCBD) of
#' samples
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' A wrapper for the [adespatial::beta.div()] function in the case of `physeq`
#' object.
#' @inheritParams clean_pq
#'
#' @param p_adjust_method (chr, default "BH"): the method used to adjust p-value
#' @param ... Others arguments passed on to [adespatial::beta.div()] function
#'
#' @return An object of class `beta.div` see [adespatial::beta.div()] function
#' for more information
#' @export
#' @seealso [plot_LCBD_pq], [adespatial::beta.div()]
#' @examples
#' if (requireNamespace("adespatial")) {
#' res <- LCBD_pq(data_fungi_sp_known, nperm = 5)
#' str(res)
#' length(res$LCBD)
#' length(res$SCBD)
#' }
#' \donttest{
#' if (requireNamespace("adespatial")) {
#' LCBD_pq(data_fungi_sp_known, nperm = 5, method = "jaccard")
#' }
#' }
#'
#' @author Adrien Taudière
#' This function is mainly a wrapper of the work of others.
#' Please make a reference to `adespatial::beta.div()` if you
#' use this function.
LCBD_pq <- function(physeq,
p_adjust_method = "BH",
...) {
physeq <- taxa_as_columns(physeq)
mat <- as.matrix(unclass(physeq@otu_table))
resBeta <- adespatial::beta.div(mat, adj = FALSE, ...)
resBeta$p.adj <-
p.adjust(resBeta$p.LCBD, method = p_adjust_method)
return(resBeta)
}
################################################################################
################################################################################
#' @title Plot and test local contributions to beta diversity (LCBD) of samples
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' A wrapper for the [adespatial::beta.div()] function in the case of `physeq`
#' object.
#' @inheritParams clean_pq
#'
#' @param p_adjust_method (chr, default "BH"): the method used to adjust p-value
#' @param pval (int, default 0.05): the value to determine the significance of
#' LCBD
#' @param sam_variables A vector of variables names present in the `sam_data`
#' slot to plot alongside the LCBD value
#' @param only_plot_significant (logical, default TRUE) Do we plot all LCBD
#' values or only the significant ones
#' @param ... Others arguments passed on to [adespatial::beta.div()] function
#'
#' @return A ggplot2 object build with the package patchwork
#' @export
#' @seealso [LCBD_pq], [adespatial::beta.div()]
#'
#' @examples
#' data(data_fungi)
#' if (requireNamespace("adespatial")) {
#' plot_LCBD_pq(data_fungi_mini,
#' nperm = 100, only_plot_significant = FALSE,
#' pval = 0.2
#' )
#' }
#' \donttest{
#' if (requireNamespace("adespatial")) {
#' plot_LCBD_pq(data_fungi_mini,
#' nperm = 100, only_plot_significant = TRUE,
#' pval = 0.2
#' )
#' if (requireNamespace("patchwork")) {
#' plot_LCBD_pq(data_fungi_mini,
#' nperm = 100, only_plot_significant = FALSE,
#' sam_variables = c("Time", "Height")
#' )
#' plot_LCBD_pq(data_fungi_mini,
#' nperm = 100, only_plot_significant = TRUE, pval = 0.2,
#' sam_variables = c("Time", "Height", "Tree_name")
#' ) &
#' theme(
#' legend.key.size = unit(0.4, "cm"),
#' legend.text = element_text(size = 10),
#' axis.title.x = element_text(size = 6)
#' )
#' }
#' }
#' }
#' @author Adrien Taudière
#' @details
#' This function is mainly a wrapper of the work of others.
#' Please make a reference to `vegan::beta.div()` if you
#' use this function.
plot_LCBD_pq <- function(physeq,
p_adjust_method = "BH",
pval = 0.05,
sam_variables = NULL,
only_plot_significant = TRUE,
...) {
resBeta <- LCBD_pq(physeq,
p_adjust_method = p_adjust_method,
...
)
sam_data <- data.frame(physeq@sam_data)
if (sum(is.na(match(
rownames(sam_data), names(resBeta$LCBD)
))) > 0) {
warning("At least one sample was removed by the beta.div function")
sam_data <-
sam_data[rownames(sam_data) %in% names(resBeta$LCBD), ]
}
resLCBD <- tibble(
"LCBD" = resBeta$LCBD,
"p.LCBD" = resBeta$p.LCBD,
"p.adj" = resBeta$p.adj,
sam_data
)
if (only_plot_significant) {
p_LCBD <-
ggplot(
filter(resLCBD, p.adj < pval),
aes(
x = LCBD,
y = reorder(Sample_names, -LCBD, sum)
)
) +
geom_point(size = 4)
if (is.null(sam_variables)) {
return(p_LCBD)
} else {
p_heatmap <- list()
for (i in seq_len(length(sam_variables))) {
p_heatmap[[i]] <- ggplot(filter(resLCBD, p.adj < pval)) +
geom_tile(inherit.aes = FALSE, aes(
y = reorder(Sample_names, -LCBD, sum),
x = i,
fill = .data[[sam_variables[[i]]]]
)) +
theme_void() +
theme(axis.title.x = element_text()) +
xlab(sam_variables[[i]])
}
p <-
p_LCBD + patchwork::wrap_plots(p_heatmap) + patchwork::plot_layout(
widths = c(3, 1), guides =
"collect"
)
return(p)
}
} else {
p_LCBD <-
ggplot(resLCBD, aes(x = LCBD, y = reorder(Sample_names, -LCBD, sum))) +
geom_point(aes(color = p.adj < pval))
if (is.null(sam_variables)) {
return(p_LCBD)
} else {
p_heatmap <- list()
for (i in seq_len(length(sam_variables))) {
p_heatmap[[i]] <- ggplot(resLCBD) +
geom_tile(inherit.aes = FALSE, aes(
y = reorder(Sample_names, -LCBD, sum),
x = i,
fill = .data[[sam_variables[[i]]]]
)) +
theme_void() +
theme(axis.title.x = element_text()) +
xlab(sam_variables[[i]])
}
p <-
p_LCBD + patchwork::wrap_plots(p_heatmap) + patchwork::plot_layout(
widths = c(3, 1), guides =
"collect"
)
return(p)
}
}
}
################################################################################
################################################################################
#' @title Plot species contributions to beta diversity (SCBD) of samples
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' A wrapper for the [adespatial::beta.div()] function in the case of `physeq`
#' object.
#' @inheritParams clean_pq
#' @param tax_level Taxonomic level to used in y axis
#' @param tax_col Taxonomic level to colored points
#' @param min_SCBD (default 0.01) the minimum SCBD value
#' to plot the taxa
#' @param ... Others arguments passed on to [adespatial::beta.div()] function
#'
#' @return A ggplot2 object build with the package patchwork
#' @export
#' @seealso [LCBD_pq], [adespatial::beta.div()]
#'
#' @examples
#' data(data_fungi)
#' if (requireNamespace("adespatial")) {
#' plot_SCBD_pq(data_fungi) +
#' geom_text(aes(label = paste(Genus, Species)), hjust = 1, vjust = 2) +
#' xlim(c(0, NA))
#' }
#' \donttest{
#' if (requireNamespace("adespatial")) {
#' plot_SCBD_pq(data_fungi, tax_level = "Class", tax_col = "Phylum", min_SCBD = 0) +
#' geom_jitter()
#' }
#' }
#'
#' @author Adrien Taudière
#' @details
#' This function is mainly a wrapper of the work of others.
#' Please make a reference to `vegan::beta.div()` if you
#' use this function.
plot_SCBD_pq <- function(physeq,
tax_level = "ASV",
tax_col = "Order",
min_SCBD = 0.01,
...) {
resBeta <- LCBD_pq(physeq, nperm = 0, ...)
tax_tab <- data.frame(physeq@tax_table)
resSCBD <- tibble(
"ASV" = taxa_names(physeq),
"SCBD" = resBeta$SCBD,
tax_tab
)
p_SCBD <- ggplot(
filter(resSCBD, SCBD > min_SCBD),
aes(
x = SCBD, y = factor(.data[[tax_level]]),
color = .data[[tax_col]]
)
) +
geom_point()
return(p_SCBD)
}
################################################################################
################################################################################
#' @title Test and plot multipatt result
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' A wrapper for the [indicspecies::multipatt()] function in the case of
#' `physeq` object.
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor in `physeq@sam_data` used to plot
#' different lines
#' @param p_adjust_method (chr, default "BH"): the method used to adjust p-value
#' @param pval (int, default 0.05): the value to determine the significance of
#' LCBD
#' @param control see `?indicspecies::multipatt()`
#' @param ... Others arguments passed on to [indicspecies::multipatt()] function
#'
#' @return A ggplot2 object
#' @export
#' @examplesIf tolower(Sys.info()[["sysname"]]) != "windows"
#' if (requireNamespace("indicspecies")) {
#' data(data_fungi)
#' data_fungi_ab <- subset_taxa_pq(data_fungi, taxa_sums(data_fungi) > 10000)
#' multipatt_pq(subset_samples(data_fungi_ab, !is.na(Time)), fact = "Time")
#' }
#' \donttest{
#' if (requireNamespace("indicspecies")) {
#' multipatt_pq(subset_samples(data_fungi_ab, !is.na(Time)),
#' fact = "Time",
#' max.order = 1, control = permute::how(nperm = 99)
#' )
#' }
#' }
#'
#' @author Adrien Taudière
#' @details
#' This function is mainly a wrapper of the work of others.
#' Please make a reference to `indicspecies::multipatt()` if you
#' use this function.
multipatt_pq <- function(physeq,
fact,
p_adjust_method = "BH",
pval = 0.05,
control = permute::how(nperm = 999),
...) {
physeq <- taxa_as_columns(physeq)
res <-
indicspecies::multipatt(as.matrix(physeq@otu_table),
physeq@sam_data[[fact]],
control = control,
...
)
res_df <- res$sign
res_df$p.adj <- p.adjust(res_df$p.value, method = p_adjust_method)
res_df$ASV_names <- rownames(res_df)
res_df_signif <-
res_df %>%
filter(p.adj < pval) %>%
tidyr::pivot_longer(cols = starts_with("s."))
p <- ggplot(
res_df_signif,
aes(
x = ASV_names,
y = name,
size = 2 * value,
color = stat
)
) +
geom_point() +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
))
return(p)
}
################################################################################
################################################################################
#' Run ANCOMBC2 on phyloseq object
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' A wrapper for the [ANCOMBC::ancombc2()] function
#'
#' @inheritParams clean_pq
#' @param fact (required) Name of the factor in `physeq@sam_data` used to plot
#' different lines
#' @param levels_fact (default NULL) The order of the level in the factor.
#' Used for reorder levels and select levels (filter out levels not present
#' en levels_fact)
#' @param tax_level The taxonomic level passed on to [ANCOMBC::ancombc2()]
#' @param ... Other arguments passed on to [ANCOMBC::ancombc2()] function.
#' @return The result of [ANCOMBC::ancombc2()] function
#' @export
#'
#' @examples
#' \donttest{
#' if (requireNamespace("mia")) {
#' data_fungi_mini@tax_table <- phyloseq::tax_table(cbind(
#' data_fungi_mini@tax_table,
#' "taxon" = taxa_names(data_fungi_mini)
#' ))
#' res_height <- ancombc_pq(
#' data_fungi_mini,
#' fact = "Height",
#' levels_fact = c("Low", "High"),
#' verbose = TRUE
#' )
#'
#' ggplot(
#' res_height$res,
#' aes(
#' y = reorder(taxon, lfc_HeightHigh),
#' x = lfc_HeightHigh,
#' color = diff_HeightHigh
#' )
#' ) +
#' geom_vline(xintercept = 0) +
#' geom_segment(aes(
#' xend = 0, y = reorder(taxon, lfc_HeightHigh),
#' yend = reorder(taxon, lfc_HeightHigh)
#' ), color = "darkgrey") +
#' geom_point()
#'
#' res_time <- ancombc_pq(
#' data_fungi_mini,
#' fact = "Time",
#' levels_fact = c("0", "15"),
#' tax_level = "Family",
#' verbose = TRUE
#' )
#' }
#' }
#' @author Adrien Taudière
#' @details
#' This function is mainly a wrapper of the work of others.
#' Please make a reference to `ANCOMBC::ancombc2()` if you
#' use this function.
ancombc_pq <- function(physeq, fact, levels_fact = NULL, tax_level = "Class", ...) {
if (!is.null(levels_fact)) {
physeq <- subset_samples_pq(physeq, as.vector(physeq@sam_data[, fact])[[1]] %in% levels_fact)
}
tse <- mia::makeTreeSummarizedExperimentFromPhyloseq(physeq)
if (!is.null(levels_fact)) {
SummarizedExperiment::colData(tse)[[fact]] <- factor(tse[[fact]], levels = levels_fact)
}
res_ancomb <- ANCOMBC::ancombc2(
data = tse,
assay_name = "counts",
fix_formula = fact,
tax_level = tax_level,
group = fact,
...
)
return(res_ancomb)
}
################################################################################
################################################################################
#' Filter ancombc_pq results
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' @param ancombc_res (required) the result of the ancombc_pq function
#' For the moment only bimodal factors are possible.
#' @param filter_passed (logical, default TRUE) Do we filter using the column
#' passed_ss? The passed_ss value is TRUE if the taxon passed the sensitivity
#' analysis, i.e., adding different pseudo-counts to 0s would not change the results.
#' @param filter_diff (logical, default TRUE) Do we filter using the column
#' diff? The diff value is TRUE if the taxon is significant
#' (has q less than alpha)
#' @param min_abs_lfc (integer, default0) Minimum absolute value to filter
#' results based on Log Fold Change. For ex. a value of 1 filter out taxa
#' for which the abundance in a given level of the modality is not at least
#' the double of the abundance in the other level.
#'
#' @return A data.frame with the same number of columns than the `ancombc_res`
#' param but with less (or equal) numbers of rows
#' @export
#'
#' @seealso [ancombc_pq()], [plot_ancombc_pq()]
#'
#' @examples
#' \donttest{
#' if (requireNamespace("mia")) {
#' data_fungi_mini@tax_table <- phyloseq::tax_table(cbind(
#' data_fungi_mini@tax_table,
#' "taxon" = taxa_names(data_fungi_mini)
#' ))
#'
#' res_time <- ancombc_pq(
#' data_fungi_mini,
#' fact = "Time",
#' levels_fact = c("0", "15"),
#' tax_level = "taxon",
#' verbose = TRUE
#' )
#'
#' signif_ancombc(res_time)
#' }
#' }
#' @details
#' This function is mainly a wrapper of the work of others.
#' Please make a reference to `ANCOMBC::ancombc2()` if you
#' use this function.
signif_ancombc <- function(ancombc_res,
filter_passed = TRUE,
filter_diff = TRUE,
min_abs_lfc = 0) {
signif_ancombc_res <- ancombc_res$res
clnames <- colnames(signif_ancombc_res)
name_modality <-
gsub(
"passed_ss", "",
clnames[grepl("passed_ss", clnames) &
!grepl("Intercept", clnames)]
)
if (filter_passed) {
signif_ancombc_res <- signif_ancombc_res %>%
filter(.data[[paste0("passed_ss", name_modality)]])
}
if (filter_diff) {
signif_ancombc_res <- signif_ancombc_res %>%
filter(.data[[paste0("diff", name_modality)]])
}
signif_ancombc_res <- signif_ancombc_res %>%
filter(abs(.data[[paste0("lfc", name_modality)]]) > min_abs_lfc)
return(signif_ancombc_res)
}
################################################################################
################################################################################
#' Plot ANCOMBC2 result for phyloseq object
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' @inheritParams clean_pq
#' @param ancombc_res (required) the result of the ancombc_pq function
#' For the moment only bimodal factors are possible.
#' @param filter_passed (logical, default TRUE) Do we filter using the column
#' passed_ss? The passed_ss value is TRUE if the taxon passed the sensitivity
#' analysis, i.e., adding different pseudo-counts to 0s would not change the results.
#' @param filter_diff (logical, default TRUE) Do we filter using the column
#' diff? The diff value is TRUE if the taxon is significant
#' (has q less than alpha)
#' @param min_abs_lfc (integer, default 0) Minimum absolute value to filter
#' results based on Log Fold Change. For ex. a value of 1 filter out taxa
#' for which the abundance in a given level of the modality is not at least
#' the double of the abundance in the other level.
#' @param tax_col The taxonomic level (must be present in `tax_table` slot)
#' to color the points
#' @param tax_label The taxonomic level (must be present in `tax_table` slot)
#' to add label
#' @param add_marginal_vioplot (logical, default TRUE) Do we add a marginal
#' vioplot representing all the taxa lfc from ancombc_res.
#' @param add_label (logical, default TRUE) Do we add a label?
#' @param add_hline_cut_lfc (logical, default NULL) Do we add two horizontal
#' lines when min_abs_lfc is set (different from zero)?
#'
#' @return A ggplot2 object. If add_marginal_vioplot is TRUE, this is a
#' patchworks of plot made using `patchwork::plot_layout()`.
#' @export
#'
#' @examples
#' \donttest{
#' if (requireNamespace("mia")) {
#' data_fungi_mini@tax_table <- phyloseq::tax_table(cbind(
#' data_fungi_mini@tax_table,
#' "taxon" = taxa_names(data_fungi_mini)
#' ))
#'
#' res_time <- ancombc_pq(
#' data_fungi_mini,
#' fact = "Time",
#' levels_fact = c("0", "15"),
#' tax_level = "taxon",
#' verbose = TRUE
#' )
#'
#' plot_ancombc_pq(data_fungi_mini, res_time,
#' filter_passed = FALSE,
#' tax_label = "Genus", tax_col = "Order"
#' )
#' plot_ancombc_pq(data_fungi_mini, res_time, tax_col = "Genus")
#' plot_ancombc_pq(data_fungi_mini, res_time,
#' filter_passed = FALSE,
#' filter_diff = FALSE, tax_col = "Family", add_label = FALSE
#' )
#' }
#' }
#' @details
#' This function is mainly a wrapper of the work of others.
#' Please make a reference to `ANCOMBC::ancombc2()` if you
#' use this function.
#' @author Adrien Taudière
plot_ancombc_pq <-
function(physeq,
ancombc_res,
filter_passed = TRUE,
filter_diff = TRUE,
min_abs_lfc = 0,
tax_col = "Genus",
tax_label = "Species",
add_marginal_vioplot = TRUE,
add_label = TRUE,
add_hline_cut_lfc = NULL) {
verify_pq(physeq)
if (is.null(add_hline_cut_lfc)) {
if (min_abs_lfc != 0) {
add_hline_cut_lfc <- TRUE
} else {
add_hline_cut_lfc <- FALSE
}
}
signif_ancombc_res <- signif_ancombc(
ancombc_res,
filter_passed = filter_passed,
filter_diff = filter_diff,
min_abs_lfc = min_abs_lfc
)
clnames <- colnames(ancombc_res$res)
name_modality <-
gsub(
"passed_ss", "",
clnames[grepl("passed_ss", clnames) &
!grepl("Intercept", clnames)]
)
taxtable <- data.frame(physeq@tax_table)
taxtable$taxon <- taxa_names(physeq)
df <-
left_join(signif_ancombc_res, taxtable, by = join_by("taxon" == "taxon"))
df[[tax_label]] <- gsub("unidentified", NA, df[[tax_label]])
p <- ggplot(
df,
aes(
y = reorder(taxon, .data[[paste0("lfc", name_modality)]]),
x = .data[[paste0("lfc", name_modality)]],
color = .data[[tax_col]]
)
) +
geom_vline(xintercept = 0) +
geom_segment(aes(
xend = 0,
y = reorder(taxon, .data[[paste0("lfc", name_modality)]]),
yend = reorder(taxon, .data[[paste0("lfc", name_modality)]])
), color = "darkgrey") +
geom_point(size = 3)
if (add_label) {
p <- p + geom_label(aes(
label = .data[[tax_label]], x =
0
))
}
if (add_hline_cut_lfc) {
p <- p +
geom_vline(
xintercept = min_abs_lfc,
color = "grey20",
lty = 2
) +
geom_vline(
xintercept = -min_abs_lfc,
color = "grey20",
lty = 2
)
}
if (add_marginal_vioplot) {
marg_vio <-
(
ggplot(data = ancombc_res$res, aes(x = .data[[paste0("lfc", name_modality)]], y = 1)) +
geom_violin(fill = transp("grey", 0.5)) +
geom_point(aes(color = .data[[paste0("diff", name_modality)]] &
.data[[paste0("passed_ss", name_modality)]]), position = position_jitter()) +
scale_color_manual(values = c("grey70", "grey20")) +
theme(legend.position = "none")
)
p <-
marg_vio / (p + xlim(layer_scales(marg_vio)$x$get_limits())) +
patchwork::plot_layout(heights = c(1, 7), axes = "collect")
}
return(p)
}
################################################################################
################################################################################
#' Show taxa which are present in only one given level of a modality
#'
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' @inheritParams clean_pq
#' @param modality (required) The name of a column present in the `@sam_data` slot
#' of the physeq object. Must be a character vector or a factor.
#' @param level (required) The level (must be present in modality) of interest
#' @param min_nb_seq_taxa (default 0 = no filter) The minimum number of sequences per taxa
#' @param min_nb_samples_taxa (default 0 = no filter) The minimum number of samples per taxa
#'
#' @return A vector of taxa names
#' @export
#'
#' @author Adrien Taudière
#' @examples
#' data_fungi_mini_woNA4height <- subset_samples(
#' data_fungi_mini,
#' !is.na(data_fungi_mini@sam_data$Height)
#' )
#' taxa_only_in_one_level(data_fungi_mini_woNA4height, "Height", "High")
taxa_only_in_one_level <- function(physeq,
modality,
level,
min_nb_seq_taxa = 0,
min_nb_samples_taxa = 0) {
if (min_nb_seq_taxa > 0) {
physeq <-
subset_taxa_pq(physeq, taxa_sums(physeq) >= min_nb_seq_taxa)
}
if (min_nb_samples_taxa > 0) {
physeq <-
subset_taxa_pq(
physeq,
taxa_sums(as_binary_otu_table(physeq)) >= min_nb_samples_taxa
)
}
physeq_merged <- clean_pq(merge_samples2(physeq, modality))
physeq_merged_only_one_level <-
subset_taxa_pq(physeq_merged, taxa_sums(as_binary_otu_table(physeq_merged)) ==
1)
physeq_merged_only_level_given <-
clean_pq(subset_samples_pq(
physeq_merged_only_one_level,
rownames(physeq_merged_only_one_level@sam_data) == level
))
return(taxa_names(physeq_merged_only_level_given))
}
################################################################################
################################################################################
#' Distribution of sequences across a factor for one taxon
#'
#' @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 taxa_name (required): the name of the taxa
#' @param digits (default = 2) integer indicating the number of decimal places
#' to be used (see `?round` for more information)
#'
#' @return a dataframe with levels as rows and information as column :
#' - the number of sequences of the taxa (nb_seq)
#' - the number of samples of the taxa (nb_samp)
#' - the mean (mean_nb_seq) and standard deviation (sd_nb_seq) of the *nb_seq*
#' - the mean (mean_nb_seq_when_present) *nb_seq* excluding samples with zero
#' - the total number of samples (nb_total_samp)
#' - the proportion of samples with the taxa
#'
#' @export
#' @author Adrien Taudière
#' @examples
#' distri_1_taxa(data_fungi, "Height", "ASV2")
#' distri_1_taxa(data_fungi, "Time", "ASV81", digits = 1)
#' @importFrom stats sd
distri_1_taxa <- function(physeq, fact, taxa_name, digits = 2) {
physeq <- taxa_as_rows(physeq)
df <-
data.frame(
"nb_seq" = tapply(
as.vector(physeq@otu_table[taxa_name, ]),
physeq@sam_data[[fact]], sum
),
"nb_samp" = tapply(
as.vector(as_binary_otu_table(physeq)@otu_table[taxa_name, ]),
physeq@sam_data[[fact]],
sum
),
"mean_nb_seq" = round(tapply(
as.vector(physeq@otu_table[taxa_name, ]),
physeq@sam_data[[fact]], mean
), digits = digits),
"sd_nb_seq" = round(tapply(
as.vector(physeq@otu_table[taxa_name, ]),
physeq@sam_data[[fact]], sd
), digits = digits)
) %>%
mutate("mean_nb_seq_when_present" = round(nb_seq / nb_samp,
digits =
digits
))
df$nb_total_samp <- table(physeq@sam_data[[fact]])
df$prop_samp <-
round(df$nb_samp / df$nb_total_samp, digits = digits)
return(df)
}
################################################################################
################################################################################
#' Partition the Variation of a phyloseq object by 2, 3, or 4 Explanatory Matrices
#' @description
#' `r lifecycle::badge("experimental")`
#' The function partitions the variation in otu_table using
#' distance (Bray per default) with respect to two, three, or four explanatory
#' tables, using
#' adjusted R² in redundancy analysis ordination (RDA) or distance-based
#' redundancy analysis. If response is a single vector, partitioning is by
#' partial regression. Collinear variables in the explanatory tables do NOT
#' have to be removed prior to partitioning. See [vegan::varpart()] for more
#' information.
#'
#' @inheritParams clean_pq
#' @param list_component (required) A named list of 2, 3 or four vectors with
#' names from the `@sam_data` slot.
#' @param dist_method (default "bray") the distance used. See
#' [phyloseq::distance()] for all available distances or run
#' [phyloseq::distanceMethodList()].
#' For "aitchison" and "robust.aitchison" distance, [vegan::vegdist()]
#' function is directly used.
#' @param dbrda_computation (logical) Do dbrda computations are runned for each
#' individual component (each name of the list component) ?
#'
#' @return an object of class "varpart", see [vegan::varpart()]
#' @export
#' @author Adrien Taudière
#' @seealso [var_par_rarperm_pq()], [vegan::varpart()], [plot_var_part_pq()]
#' @examples
#' \donttest{
#' if (requireNamespace("vegan")) {
#' data_fungi_woNA <-
#' subset_samples(data_fungi, !is.na(Time) & !is.na(Height))
#' res_var <- var_par_pq(data_fungi_woNA,
#' list_component = list(
#' "Time" = c("Time"),
#' "Size" = c("Height", "Diameter")
#' ),
#' dbrda_computation = TRUE
#' )
#' }
#' }
#' @details
#' This function is mainly a wrapper of the work of others.
#' Please make a reference to `vegan::varpart()` if you
#' use this function.
var_par_pq <-
function(physeq,
list_component,
dist_method = "bray",
dbrda_computation = TRUE) {
physeq <- taxa_as_columns(physeq)
verify_pq(physeq)
if (dist_method %in% c("robust.aitchison", "aitchison")) {
dist_physeq <-
vegan::vegdist(as(otu_table(physeq, taxa_are_rows = FALSE), "matrix"),
method = dist_method
)
} else {
dist_physeq <- phyloseq::distance(physeq, method = dist_method)
}
for (i in 1:length(list_component)) {
assign(
names(list_component)[i],
as.data.frame(unclass(physeq@sam_data[, list_component[[i]]]))
)
}
if (length(list_component) == 2) {
res_varpart <-
vegan::varpart(
dist_physeq,
eval(sym(names(list_component)[1])),
eval(sym(names(list_component)[2]))
)
} else if (length(list_component) == 3) {
res_varpart <-
vegan::varpart(
dist_physeq,
eval(sym(names(list_component)[1])),
eval(sym(names(list_component)[2])),
eval(sym(names(list_component)[3]))
)
} else if (length(list_component) == 4) {
res_varpart <-
vegan::varpart(
dist_physeq,
eval(sym(names(list_component)[1])),
eval(sym(names(list_component)[2])),
eval(sym(names(list_component)[3])),
eval(sym(names(list_component)[4]))
)
} else {
stop("The list_component must be of length 2, 3 or 4")
}
if (dbrda_computation) {
res_varpart$dbrda_result <- list()
for (i in 1:length(list_component)) {
res_varpart$dbrda_result[[i]] <-
anova(vegan::dbrda(
as.formula(paste0(
"dist_physeq ~ ", paste(c(list_component[[i]]), collapse = " + ")
)),
data = as.data.frame(unclass(
physeq@sam_data
))
))
}
}
res_varpart$Xnames <- names(list_component)
return(res_varpart)
}
################################################################################
################################################################################
#' Partition the Variation of a phyloseq object with rarefaction permutations
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' This is an extension of the function [var_par_pq()]. The main addition is
#' the computation of nperm permutations with rarefaction even depth by
#' sample. The return object
#'
#'
#' @inheritParams clean_pq
#' @param list_component (required) A named list of 2, 3 or four vectors with
#' names from the `@sam_data` slot.
#' @param dist_method (default "bray") the distance used. See
#' [phyloseq::distance()] for all available distances or run
#' [phyloseq::distanceMethodList()].
#' For aitchison and robust.aitchison distance, [vegan::vegdist()]
#' function is directly used.#' @param fill_bg
#' @param nperm (int) The number of permutations to perform.
#' @param quantile_prob (float, `[0:1]`) the value to compute the quantile.
#' Minimum quantile is compute using 1-quantile_prob.
#' @param dbrda_computation (logical) Do dbrda computations are runned for each
#' individual component (each name of the list component) ?
#' @param dbrda_signif_pval (float, `[0:1]`) The value under which the dbrda is
#' considered significant.
#' @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 progress_bar (logical, default TRUE) Do we print progress during
#' the calculation?
#'
#' @return A list of class varpart with additional information in the
#' `$part$indfract` part. Adj.R.square is the mean across permutation.
#' Adj.R.squared_quantil_min and Adj.R.squared_quantil_max represent
#' the quantile values of adjuste R squared
#' @export
#' @seealso [var_par_pq()], [vegan::varpart()], [plot_var_part_pq()]
#' @author Adrien Taudière
#' @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
#' )
#' }
#' }
#' @details
#' This function is mainly a wrapper of the work of others.
#' Please make a reference to `vegan::varpart()` if you
#' use this function.
var_par_rarperm_pq <-
function(physeq,
list_component,
dist_method = "bray",
nperm = 99,
quantile_prob = 0.975,
dbrda_computation = FALSE,
dbrda_signif_pval = 0.05,
sample.size = min(sample_sums(physeq)),
verbose = FALSE,
progress_bar = TRUE) {
physeq <- taxa_as_columns(physeq)
verify_pq(physeq)
if (progress_bar) {
pb <- txtProgressBar(
min = 0,
max = nperm,
style = 3,
width = 50,
char = "="
)
}
if (dist_method %in% c("robust.aitchison", "aitchison")) {
dist_physeq <-
vegdist(as(otu_table(physeq, taxa_are_rows = FALSE), "matrix"), method = dist_method)
} else {
dist_physeq <- phyloseq::distance(physeq, method = dist_method)
}
res_perm <- list()
for (i in 1:nperm) {
res_perm[[i]] <-
var_par_pq(
physeq =
rarefy_even_depth(
physeq,
rngseed = i,
sample.size = sample.size,
verbose = verbose
),
list_component = list_component,
dist_method = dist_method,
dbrda_computation = dbrda_computation
)
if (progress_bar) {
setTxtProgressBar(pb, i)
}
}
res_varpart <- var_par_pq(
physeq = physeq,
list_component = list_component,
dist_method = dist_method,
dbrda_computation = dbrda_computation
)
if (dbrda_computation) {
res_varpart$dbrda_result_prop_pval_signif <-
rowSums(sapply(res_perm, function(x) {
sapply(x$dbrda_result, function(xx) {
xx$`Pr(>F)`[[1]]
})
}) < dbrda_signif_pval) / nperm
}
res_varpart$part$indfract$R.square <-
rowMeans(sapply(res_perm, function(x) {
(x$part$indfract$R.square)
}))
res_varpart$part$indfract$R.square_quantil_max <-
apply(sapply(res_perm, function(x) {
(x$part$indfract$R.square)
}), 1, function(xx) {
quantile(xx, probs = quantile_prob, na.rm = TRUE)
})
res_varpart$part$indfract$R.square_quantil_min <-
apply(sapply(res_perm, function(x) {
(x$part$indfract$R.square)
}), 1, function(xx) {
quantile(xx, probs = 1 - quantile_prob, na.rm = TRUE)
})
res_varpart$part$indfract$Adj.R.square <-
rowMeans(sapply(res_perm, function(x) {
(x$part$indfract$Adj.R.square)
}))
res_varpart$part$indfract$Adj.R.squared_quantil_max <-
apply(sapply(res_perm, function(x) {
(x$part$indfract$Adj.R.square)
}), 1, function(xx) {
quantile(xx, probs = quantile_prob, na.rm = TRUE)
})
res_varpart$part$indfract$Adj.R.squared_quantil_min <-
apply(sapply(res_perm, function(x) {
(x$part$indfract$Adj.R.square)
}), 1, function(xx) {
quantile(xx, probs = 1 - quantile_prob, na.rm = TRUE)
})
return(res_varpart)
res_varpart$Xnames <- names(list_component)
return(res_varpart)
}
################################################################################
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.