Skip to content

Commit

Permalink
BiocCheck suggestions
Browse files Browse the repository at this point in the history
  • Loading branch information
@xec-cm
xec-cm committed Nov 14, 2023
1 parent 0e0fa96 commit 9b99bb8
Show file tree
Hide file tree
Showing 10 changed files with 225 additions and 201 deletions.
63 changes: 38 additions & 25 deletions R/ancom.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -269,19 +269,7 @@ run_ancom <- function(rec,
eval()

## Remove OTUs with zero variance
to_remove <-
phyloseq::sample_data(s_phy) %>%
tibble::as_tibble(rownames = "sample_id") %>%
dplyr::group_split(!!dplyr::sym(vars)) %>%
purrr::map_dfr(~ {
.x %>%
dplyr::pull(sample_id) %>%
phyloseq::prune_samples(s_phy) %>%
phyloseq::filter_taxa(function(x) var(x) > 0) %>%
tibble::as_tibble(rownames = "otu_id")
}) %>%
dplyr::filter(value == FALSE) %>%
dplyr::pull(otu_id)
to_remove <- rm_zero_variance(s_phy, vars)

s_phy <-
rownames(phyloseq::otu_table(s_phy)) %>%
Expand Down Expand Up @@ -309,18 +297,43 @@ run_ancom <- function(rec,
trend = trend
)

res$res %>%
dplyr::select(taxa = taxon, dplyr::contains(!!var)) %>%
dplyr::right_join(tax_table(rec), ., by = "taxa") %>%
stats::setNames(
stringr::str_remove_all(names(.), stringr::str_c("_", var, ".*"))
) %>%
dplyr::mutate(
comparison = stringr::str_c(comparison, collapse = "_"),
var = var,
effect = lfc,
signif = diff
)
ancom_stats_tbl(res$res, var, rec, comparison)
})
})
}

#' @noRd
#' @keywords internal
#' @autoglobal
ancom_stats_tbl <- function(ancom_res, var, rec, comparison) {
ancom_res %>%
dplyr::select(taxa = taxon, dplyr::contains(!!var)) %>%
dplyr::right_join(tax_table(rec), ., by = "taxa") %>%
stats::setNames(
stringr::str_remove_all(names(.), stringr::str_c("_", var, ".*"))
) %>%
dplyr::mutate(
comparison = stringr::str_c(comparison, collapse = "_"),
var = var,
effect = lfc,
signif = diff
)
}

#' @noRd
#' @keywords internal
#' @autoglobal
rm_zero_variance <- function(phy, vars) {
phyloseq::sample_data(phy) %>%
tibble::as_tibble(rownames = "sample_id") %>%
dplyr::group_split(!!dplyr::sym(vars)) %>%
purrr::map_dfr(~ {
.x %>%
dplyr::pull(sample_id) %>%
phyloseq::prune_samples(phy) %>%
phyloseq::filter_taxa(function(x) var(x) > 0) %>%
tibble::as_tibble(rownames = "otu_id")
}) %>%
dplyr::filter(value == FALSE) %>%
dplyr::pull(otu_id)
}
155 changes: 86 additions & 69 deletions R/corncob.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -266,82 +266,99 @@ run_corncob <- function(rec,
)
}, error = function(e) { conditionMessage(e) })

## Skip error for no convergenc
if (!methods::is(corncob_res, "differentialTest")) {
if (stringr::str_detect(corncob_res, "failed to converge")) {
rlang::abort(c(
glue::glue(
"{crayon::bgMagenta('corncob')}: All models failed to ",
"converge!"
),
glue::glue(
"{crayon::bgMagenta('corncob')}: If you are seeing this, it",
" is likely that your model is overspecified. This occurs ",
"when your sample size is not large enough to estimate all ",
"the parameters of your model. This is most commonly due to",
" categorical variables that include many categories."
),
glue::glue(
"Please remove or edit the ",
"{crayon::bgMagenta('step_corncob()')} and rerun."
)
))
} else {
rlang::abort(
c(
glue::glue("{crayon::bgMagenta('corncob')}: Internal error!"),
glue::glue(
"Please remove or edit the ",
"{crayon::bgMagenta('step_corncob()')} and rerun."
),
glue::glue(
"Please report this bug on GitHub: ",
"https://github.com/MicrobialGenomics-IrsicaixaOrg/dar/",
"issues"
)
)
)
}
}
## Skip error for no convergence
check_non_convergence(corncob_res)

signif_taxa <- corncob::otu_to_taxonomy(
OTU = corncob_res$significant_taxa,
data = corncob_res$data,
level = tax_level %>% dplyr::pull()
) %>% stringr::str_c(" (", corncob_res$significant_taxa, ")")

qval <- stats::qnorm(.975)
corncob_res$significant_models %>%
purrr::map2_dfr(stats::na.omit(corncob_res$p_fdr), function(m, p) {
coefs_mu <- m$coefficients[seq_len(m$np.mu), , drop = FALSE]
coefs_mu <- coefs_mu[2, , drop = FALSE]

tibble::tibble(
log2FC = coefs_mu[1, 1],
log2FC_min = coefs_mu[1, 1] - qval * coefs_mu[1, 2],
log2FC_xmax = coefs_mu[1, 1] + qval * coefs_mu[1, 2],
padj = p
)
}) %>%
dplyr::mutate(
taxa = signif_taxa,
comparison = stringr::str_c(comparison, collapse = "_"),
var = var
) %>%
tidyr::separate(
taxa, c("taxa", "taxa_id"), sep = " ", remove = TRUE
) %>%
dplyr::mutate(
taxa_id = stringr::str_remove_all(taxa_id, "[(]|[)]")
) %>%
dplyr::mutate(
effect = log2FC,
signif = ifelse(
padj < fdr_cutoff & abs(log2FC) >= log2FC,
TRUE,
FALSE
)
)
corncob_stats_tbl(
corncob_res,
signif_taxa,
comparison,
var,
log2FC,
fdr_cutoff
)
})
})
}

#' @noRd
#' @keywords internal
#' @autoglobal
corncob_stats_tbl <- function(corncob_res,
signif_taxa,
comparison,
var,
log2FC,
fdr_cutoff) {

qval <- stats::qnorm(0.975)
corncob_res$significant_models %>%
purrr::map2_dfr(stats::na.omit(corncob_res$p_fdr), function(m, p) {
coefs_mu <- m$coefficients[seq_len(m$np.mu), , drop = FALSE]
coefs_mu <- coefs_mu[2, , drop = FALSE]
tibble::tibble(
log2FC = coefs_mu[1, 1],
log2FC_min = coefs_mu[1, 1] - qval * coefs_mu[1, 2],
log2FC_xmax = coefs_mu[1, 1] + qval * coefs_mu[1, 2],
padj = p
)
}) %>%
dplyr::mutate(
taxa = signif_taxa,
comparison = stringr::str_c(comparison, collapse = "_"),
var = var
) %>%
tidyr::separate(taxa, c("taxa", "taxa_id"), sep = " ", remove = TRUE) %>%
dplyr::mutate(taxa_id = stringr::str_remove_all(taxa_id, "[(]|[)]")) %>%
dplyr::mutate(
effect = log2FC,
signif = ifelse(padj < fdr_cutoff & abs(log2FC) >= log2FC, TRUE, FALSE)
)
}

#' @noRd
#' @keywords internal
#' @autoglobal
check_non_convergence <- function(corncob_res) {
if (!methods::is(corncob_res, "differentialTest")) {
if (stringr::str_detect(corncob_res, "failed to converge")) {
rlang::abort(c(
glue::glue(
"{crayon::bgMagenta('corncob')}: All models failed to converge!"
),
glue::glue(
"{crayon::bgMagenta('corncob')}: If you are seeing this, it",
" is likely that your model is overspecified. This occurs ",
"when your sample size is not large enough to estimate all ",
"the parameters of your model. This is most commonly due to",
" categorical variables that include many categories."
),
glue::glue(
"Please remove or edit the ",
"{crayon::bgMagenta('step_corncob()')} and rerun."
)
))
} else {
rlang::abort(
c(
glue::glue("{crayon::bgMagenta('corncob')}: Internal error!"),
glue::glue(
"Please remove or edit the ",
"{crayon::bgMagenta('step_corncob()')} and rerun."
),
glue::glue(
"Please report this bug on GitHub: ",
"https://github.com/MicrobialGenomics-IrsicaixaOrg/dar/",
"issues"
)
)
)
}
}
}
51 changes: 28 additions & 23 deletions R/deseq2.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ methods::setGeneric(
name = "step_deseq",
def = function(rec,
test = "Wald",
fitType = "parametric",
fitType = "local",
betaPrior = FALSE,
type = "ashr",
max_significance = 0.05,
Expand Down Expand Up @@ -194,32 +194,37 @@ run_deseq <- function(rec,
phy <- phyloseq::rarefy_even_depth(phy, rngseed = 1234, verbose = FALSE)
}

phy <- phyloseq::tax_glom(phy, taxrank = tax_level, NArm = FALSE)

phy <- phyloseq::tax_glom(phy, taxrank = tax_level, NArm = FALSE)
phyloseq::sample_data(phy) <-
to_tibble(phyloseq::sample_data(phy), "sample_id") %>%
dplyr::mutate(dplyr::across(where(is.character), as.factor)) %>%
data.frame(row.names = 1) %>%
phyloseq::sample_data()

vars %>%
purrr::set_names() %>%
purrr::map(function(var) {
dds <-
suppressWarnings(
suppressMessages(
phy %>%
phyloseq::phyloseq_to_deseq2(
design = stats::as.formula(stringr::str_c("~", var))
) %>%
DESeq2::estimateSizeFactors(geoMeans = apply(
X = DESeq2::counts(.),
MARGIN = 1,
FUN = function(x) {
exp(sum(log(x[x > 0]), na.rm = TRUE) / length(x))
}
)) %>%
DESeq2::DESeq(
fitType = fitType,
test = test,
betaPrior = betaPrior,
quiet = TRUE
)
dds <- suppressMessages(
phy %>%
phyloseq::phyloseq_to_deseq2(
design = stats::as.formula(stringr::str_c("~", var))
)
)

dds <-
dds %>%
DESeq2::estimateSizeFactors(geoMeans = apply(
X = DESeq2::counts(.),
MARGIN = 1,
FUN = function(x) {
exp(sum(log(x[x > 0]), na.rm = TRUE) / length(x))
}
)) %>%
DESeq2::DESeq(
fitType = fitType,
test = test,
betaPrior = betaPrior,
quiet = TRUE
)

contrasts_df <- get_comparisons(var, phy, as_list = FALSE, n_cut = 1)
Expand Down
15 changes: 8 additions & 7 deletions R/globals.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ utils::globalVariables(c(
# <intersection_df>
# <cool>
# <run_ancom>
# <ancom_stats_tbl>
# <run_bake>
# <run_deseq>
# <prepro_lefse>
Expand Down Expand Up @@ -39,7 +40,7 @@ utils::globalVariables(c(
"feature",
# <prepro_lefse>
"Kingdom",
# <run_ancom>
# <ancom_stats_tbl>
"lfc",
# <run_deseq>
"log2FoldChange",
Expand All @@ -66,9 +67,9 @@ utils::globalVariables(c(
"no_zero",
# <run_lefse>
"otu",
# <run_ancom>
# <rm_zero_variance>
"otu_id",
# <run_corncob>
# <corncob_stats_tbl>
# <run_deseq>
# <run_metagenomeseq>
# <run_wilcox>
Expand Down Expand Up @@ -103,7 +104,7 @@ utils::globalVariables(c(
# <recipe>
"tax_lev",
# <run_bake>
# <run_corncob>
# <corncob_stats_tbl>
# <.abundance_boxplot>
# <.abundance_heatmap>
# <mutual_plt>
Expand All @@ -112,7 +113,7 @@ utils::globalVariables(c(
# <intersection_df>
# <run_aldex>
# <run_bake>
# <run_corncob>
# <corncob_stats_tbl>
# <zero_otu>
# <zero_otu>
# <run_maaslin>
Expand All @@ -124,13 +125,13 @@ utils::globalVariables(c(
# <mutual_plt>
# <.otu_method_count>
"taxa_id",
# <run_ancom>
# <ancom_stats_tbl>
"taxon",
# <zero_otu>
# <zero_otu>
# <exclusion_plt>
"total",
# <run_ancom>
# <rm_zero_variance>
# <zero_otu>
# <zero_otu>
# <kruskal_test>
Expand Down
9 changes: 3 additions & 6 deletions R/metagenomeseq.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,12 +199,9 @@ run_metagenomeseq <- function(rec,
unique() %>%
phyloseq::prune_taxa(phy2)

suppressMessages(
mr_obj <-
phy2 %>%
phyloseq_to_MRexperiment() %>%
metagenomeSeq::cumNorm(p = metagenomeSeq::cumNormStat(.))
)
mr_exp <- phyloseq_to_MRexperiment(phy2)
p <- suppressMessages(metagenomeSeq::cumNormStat(mr_exp))
mr_obj <- metagenomeSeq::cumNorm(mr_exp, p = p)

vct_var <- phyloseq::sample_data(phy2)[[var]]
norm_factor <- metagenomeSeq::normFactors(mr_obj)
Expand Down
Loading

0 comments on commit 9b99bb8

Please sign in to comment.