From 138ff2b54f6b93357f17bdccdd5cd6267666e088 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bol=C3=ADvar=20Aponte=20Rol=C3=B3n?= <64099468+jibarozzo@users.noreply.github.com> Date: Wed, 6 May 2026 08:56:28 -0500 Subject: [PATCH 01/10] Release/v2.0.4 to `main` (#104) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * fix: Update vignette and docs for BRCore 2.0.2 CRAN submission * chore: remove parallel dependencies and clean up docs for CRAN v2.0.2 - Removal of "multiple" rarefaction concept. Replace by "multiple iterations" * docs: update cran-comments.md for v2.0.2 resubmission * Increment version number to 2.0.2 * docs: expand package description with methodology references and citations * docs: improve documentation formatting and examples across neutral model functions * docs: improve examples, add `seealso` links, and remove \donttest wrappers * chore: Remove parallel computing utilities from pkgdown reference and fix OTU table formatting * docs: Add `switchgrass_core` precomputed dataset and refactor examples to use it. Updated citation style to APA 7th edition. * docs: update `seealso` links and remove global warning suppression in `sncm.fit()` * fix: update DOI links to use canonical doi: format in docs and DESCRIPTION * docs: add switchgrass_core to pkgdown reference * chore: update CRAN comments for v2.0.3 resubmission and add Quarto to .gitignore * Increment version number to 2.0.3 * fix: Suppress warnings from `dnorm` calls in `sncm.fit` likelihood functions * docs: Add `scnm.fit()` to README function list * Update release notes * Increment version number to 2.0.4 * fix: typo `scnm.fit` → `sncm.fit` and use HTTPS for pkgdown URL * fix: Quote 'BRCore' in DESCRIPTION and fix README badge links * Increment version number to 2.0.4.9000 --- .gitignore | 2 + DESCRIPTION | 9 +- NAMESPACE | 12 - NEWS.md | 24 +- R/add_rarefaction_metrics.R | 25 +- R/brcore_theme.R | 2 - R/data.R | 54 ++- R/fit_neutral_model.R | 41 +- R/identify_core.R | 55 ++- R/multi_rarefy.R | 41 +- R/parallel_helpers.R | 189 --------- R/plot_abundance_occupancy.R | 19 +- R/plot_core_distribution.R | 23 +- R/plot_identified_core.R | 38 +- R/plot_neutral_model.R | 18 +- R/plot_rarefaction_metrics.R | 5 +- R/plot_variance_propagation.R | 28 +- R/sncm.fit.R | 103 ++--- R/update_otu_table.R | 19 +- README.Rmd | 13 +- README.md | 15 +- _pkgdown.yml | 14 +- cran-comments.md | 48 ++- data/switchgrass_core.rda | Bin 0 -> 35964 bytes man/add_rarefaction_metrics.Rd | 26 +- man/bcse.Rd | 8 +- man/bean.Rd | 6 +- man/detect_hpc_type.Rd | 19 - man/fit_neutral_model.Rd | 40 +- man/get_available_cores.Rd | 23 -- man/identify_core.Rd | 36 +- man/is_hpc_environment.Rd | 18 - man/mimulus.Rd | 6 +- man/multi_rarefy.Rd | 21 +- man/plot_abundance_occupancy.Rd | 19 +- man/plot_core_distribution.Rd | 23 +- man/plot_identified_core.Rd | 26 +- man/plot_neutral_model.Rd | 18 +- man/plot_rarefaction_metrics.Rd | 5 +- man/plot_variance_propagation.Rd | 12 +- man/setup_parallel_backend.Rd | 22 - man/sncm.fit.Rd | 96 +++-- man/switchgrass.Rd | 6 +- man/switchgrass_core.Rd | 35 ++ man/update_otu_table.Rd | 9 +- tests/testthat/test-vignette-workflow.R | 2 +- .../testthat/test_sets/regenerate_test_data.R | 2 +- vignettes/BRCore-vignette.R | 52 ++- vignettes/BRCore-vignette.Rmd | 383 +++++++++++++----- vignettes/BRCore-vignette.Rmd.orig | 139 +++---- vignettes/images/fig1_plot_metrics-1.png | Bin 41669 -> 43915 bytes .../images/fig2_rarefaction_variance-1.png | Bin 34055 -> 35129 bytes vignettes/images/fig3_identified_core-1.png | Bin 20995 -> 21092 bytes vignettes/images/fig4_plot_increase-1.png | Bin 37955 -> 39960 bytes vignettes/images/fig5_plot_elbow-1.png | Bin 37517 -> 39349 bytes vignettes/images/fig6_plot_type_bar-1.png | Bin 23005 -> 38060 bytes vignettes/images/fig7_plot_type_line-1.png | Bin 46795 -> 62802 bytes vignettes/images/fig8_plot_type_heatmap-1.png | Bin 33767 -> 46924 bytes vignettes/images/fig9_plot_neutral_fit-1.png | Bin 48493 -> 51171 bytes 59 files changed, 908 insertions(+), 941 deletions(-) delete mode 100644 R/parallel_helpers.R create mode 100644 data/switchgrass_core.rda delete mode 100644 man/detect_hpc_type.Rd delete mode 100644 man/get_available_cores.Rd delete mode 100644 man/is_hpc_environment.Rd delete mode 100644 man/setup_parallel_backend.Rd create mode 100644 man/switchgrass_core.Rd diff --git a/.gitignore b/.gitignore index f1a61b7..b48f8b3 100644 --- a/.gitignore +++ b/.gitignore @@ -41,3 +41,5 @@ renv.lock /doc/ /Meta/ docs + +/.quarto/ diff --git a/DESCRIPTION b/DESCRIPTION index 7f9349d..a54712a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BRCore Title: Unified Framework for Identification and Ecological Interpretation of Microbial Data from Bioenergy Research Centers -Version: 2.0.1 +Version: 2.0.4.9000 Authors@R: c( person("Bolívar", "Aponte Rolón", email = "bolaponte@pm.me", @@ -27,10 +27,10 @@ Authors@R: c( person("Great Lakes Bioenergy Research Center", role = c("cph", "fnd")) ) -Description: A unified framework for identification and ecological interpretation of core microbiomes across time and space, enhancing robustness and reproducibility in microbiome data analysis. +Description: A unified framework for identification and ecological interpretation of core microbiomes across time and space, enhancing robustness and reproducibility in microbiome data analysis. 'BRCore' implements the workflow proposed by Shade and Stopnisek (2019) and incorporates additional rarefaction steps. The proposed workflow aims to identify persistent microbiomes using abundance-occupancy distributions and neutral community model fitting. For more details on abundance-occupancy distributions see Shade A, Stopnisek N (2019) , for neutral models, see Sloan et al. (2006) and Burns et al. (2015) . License: MIT + file LICENSE Encoding: UTF-8 -Language: en +Language: en-US Roxygen: list(markdown = TRUE) RoxygenNote: 7.3.3 Imports: @@ -39,12 +39,9 @@ Imports: ggplot2 (>= 4.0.1), ggpubr (>= 0.6.2), ggrepel (>= 0.9.6), - ggsci (>= 4.1.0), Hmisc (>= 5.2.4), magrittr (>= 2.0.4), minpack.lm (>= 1.2.4), - parallel (>= 4.5.2), - parallelly (>= 1.45.1), phyloseq (>= 1.54.0), rlang (>= 1.1.6), scales (>= 1.4.0), diff --git a/NAMESPACE b/NAMESPACE index 433a4af..7de4704 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,11 +2,8 @@ export("%>%") export(add_rarefaction_metrics) -export(detect_hpc_type) export(fit_neutral_model) -export(get_available_cores) export(identify_core) -export(is_hpc_environment) export(multi_rarefy) export(plot_abundance_occupancy) export(plot_core_distribution) @@ -14,7 +11,6 @@ export(plot_identified_core) export(plot_neutral_model) export(plot_rarefaction_metrics) export(plot_variance_propagation) -export(setup_parallel_backend) export(sncm.fit) export(update_otu_table) importFrom(Hmisc,binconf) @@ -83,17 +79,9 @@ importFrom(ggplot2,theme_bw) importFrom(ggplot2,theme_classic) importFrom(ggpubr,ggarrange) importFrom(ggrepel,geom_text_repel) -importFrom(ggsci,scale_fill_npg) importFrom(grid,unit) importFrom(magrittr,"%>%") importFrom(minpack.lm,nlsLM) -importFrom(parallel,clusterEvalQ) -importFrom(parallel,clusterExport) -importFrom(parallel,makeCluster) -importFrom(parallel,parLapply) -importFrom(parallel,stopCluster) -importFrom(parallelly,availableCores) -importFrom(parallelly,makeClusterPSOCK) importFrom(phyloseq,"sample_data<-") importFrom(phyloseq,otu_table) importFrom(phyloseq,phy_tree) diff --git a/NEWS.md b/NEWS.md index f24139a..bf6468e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,9 +1,31 @@ +# BRCore (development version) +* Fixes to quote 'BRCore' in DESCRIPTION and typo fix of `sncm.fit` in _pkgdown.yaml + +# BRCore 2.0.4 +Date: 2026-04-30 + +* Fixes to documentation: + * DOI links updated to correct format + * New example for `scnm.fit()` + * New `@seealso` cross-refences in documentation + * Unwrapped examples from `\dontrun` and `\donttest` + * Gained new `swithchgrass_core` object, identified microbial core, for examples + * New APA 7th edition style for references + +# BRCore 2.0.2 +Date: 2026-04-28 + +* Fixes to vignette and function documentation prior to CRAN submission. +* Removed "multiple rarefaction" concept from documentation. +* Removed `parallel_helpers.R` since we have moved away from parallel implementations. +* Cleaned NAMESPACE and DESCRIPTION. + # BRCore 2.0.1 Date: 2026-04-25 * Resubmission to CRAN to address auto-check issues on Debian and Windows. -BRCore 2.0.0 +# BRCore 2.0.0 Date: 2026-04-24 ## Breaking Changes diff --git a/R/add_rarefaction_metrics.R b/R/add_rarefaction_metrics.R index b1eaff0..b9c6cd9 100644 --- a/R/add_rarefaction_metrics.R +++ b/R/add_rarefaction_metrics.R @@ -28,29 +28,34 @@ #' 0.95, means that 5% of the reads in that sample are from OTUs #' that appear only once. #' +#' @seealso [plot_rarefaction_metrics()] for visualizing these metrics, and +#' [multi_rarefy()] for performing rarefaction on a `phyloseq` object. +#' #' @examples -#' \donttest{ #' library(phyloseq) #' library(BRCore) -#' # From an object class "phyloseq" with added alpha metrics #' #' data("bcse", package = "BRCore") +#' +#' # Adding metrics to a "phyloseq" object #' bcse_metrics <- add_rarefaction_metrics(data = bcse) -#' sample_data(bcse_metrics) +#' sample_data(bcse_metrics)|> +#' head(10) #' -#' # From a class "data.frame" count table object #' +#' # Adding metrics to a "data.frame" count table object #' bcse_otutable <- as.data.frame( -#' as(otu_table(bcse), "matrix") +#' as.matrix(otu_table(bcse)) #' ) -#' test_otutable_metrics <- add_rarefaction_metrics( +#' +#' bcse_otutable_metrics <- add_rarefaction_metrics( #' data = bcse_otutable #' ) -#' test_otutable_metrics[ -#' utils::tail(seq_len(nrow(test_otutable_metrics)), 10), -#' utils::tail(seq_len(ncol(test_otutable_metrics)), 20) +#' bcse_otutable_metrics[ +#' head(seq_len(nrow(bcse_otutable_metrics)), 10), +#' tail(seq_len(ncol(bcse_otutable_metrics)), 20) #' ] -#' } +#' #' #' @importFrom phyloseq otu_table sample_data sample_data<- taxa_are_rows #' @importFrom tibble rownames_to_column column_to_rownames diff --git a/R/brcore_theme.R b/R/brcore_theme.R index 9bf0f36..a9c32f4 100644 --- a/R/brcore_theme.R +++ b/R/brcore_theme.R @@ -34,7 +34,6 @@ #' to any \code{ggplot} with \code{+}. #' #' @examples -#' \donttest{ #' library(ggplot2) #' ggplot(mtcars, aes(wt, mpg)) + #' geom_point() + @@ -47,7 +46,6 @@ #' dashed_grid = TRUE, #' extra_themes = list(no_legend = theme(legend.position = "none")) #' ) -#' } #' #' @importFrom ggplot2 theme_classic theme element_text element_line #' element_blank element_rect diff --git a/R/data.R b/R/data.R index d0e6f5e..1f63fcf 100644 --- a/R/data.R +++ b/R/data.R @@ -7,13 +7,39 @@ #' @docType data #' @format An object class "phyloseq". #' -#' @references Shade A, Stopnisek N. Abundance-occupancy distributions to -#' prioritize plant core microbiome membership. Curr Opin Microbiol. 2019 Jun -#' 49:50-58. doi:10.1016/j.mib.2019.09.008. Epub 2019 Nov 10. PMID: 31715441. +#' @references Shade A, Stopnisek N (2019) Abundance-occupancy distributions to +#' prioritize plant core microbiome membership. Current Opinion in +#' Microbiology, 49:50-58 #' #' @source Internal test dataset. "switchgrass" +#' Identified core microbiome members for the switchgrass dataset +#' +#' @description A phyloseq containing the identified core microbiome members +#' for the switchgrass dataset. This object is used for BRCore function +#' examples. +#' It was generated by running the following on BRCore 2.0.2 (2026-04-30): +#' ```R +#' switchgrass_core <- identify_core( +#' physeq_obj = switchgrass, +#' priority_var = "sampling_date", +#' increase_value = 0.02, +#' seed = 092825 +#' ) +#' ``` +#' +#' @usage data("switchgrass_core") +#' @docType data +#' @format An object class "phyloseq". +#' +#' @references Shade A, Stopnisek N (2019) Abundance-occupancy distributions to +#' prioritize plant core microbiome membership. Current Opinion in +#' Microbiology, 49:50-58 +#' +#' @source Internal test dataset. +"switchgrass_core" + #' 16S amplicon dataset from yellow monkeyflower (Mimulus guttatus) #' #' @description A phyloseq containing 11085 ASVs (Amplicon Sequence Variants) @@ -23,9 +49,9 @@ #' @docType data #' @format An object class "phyloseq". #' -#' @references Shade A, Stopnisek N. Abundance-occupancy distributions to -#' prioritize plant core microbiome membership. Curr Opin Microbiol. 2019 Jun -#' 49:50-58. doi:10.1016/j.mib.2019.09.008. Epub 2019 Nov 10. PMID: 31715441 +#' @references Shade A, Stopnisek N (2019) Abundance-occupancy distributions to +#' prioritize plant core microbiome membership. Current Opinion in +#' Microbiology, 49:50-58 #' #' @source Internal test dataset. "mimulus" @@ -39,9 +65,10 @@ #' @docType data #' @format An object class "phyloseq". #' -#' @references Shade A, Stopnisek N. Abundance-occupancy distributions to -#' prioritize plant core microbiome membership. Curr Opin Microbiol. 2019 Jun -#' 49:50-58. doi:10.1016/j.mib.2019.09.008. Epub 2019 Nov 10. PMID: 31715441. +#' @references Shade A, Stopnisek N (2019) Abundance-occupancy distributions to +#' prioritize plant core microbiome membership. Current Opinion in +#' Microbiology, 49:50-58 +#' #' #' @source Internal test dataset. "bean" @@ -57,10 +84,11 @@ #' @docType data #' @format An object class "phyloseq". #' -#' @references Haan NL, Benucci GNM, Fiser CM, Bonito G, Landis DA. -#' Contrasting effects of bioenergy crops on biodiversity. Sci Adv. 2023 Sep -#' 22;9(38):eadh7960. doi: 10.1126/sciadv.adh7960. Epub 2023 Sep 22. PMID: -#' 37738354; PMCID: PMC10516493. +#' @references Haan, N. L., Benucci, G. N. M., Fiser, C. M., Bonito, G., & +#' Landis, D. A. (2023). Contrasting effects of bioenergy crops on +#' biodiversity. Science Advances, 9(38), eadh7960. +#' + #' #' @source Internal test dataset. "bcse" diff --git a/R/fit_neutral_model.R b/R/fit_neutral_model.R index 32aac9f..bc5c8dd 100644 --- a/R/fit_neutral_model.R +++ b/R/fit_neutral_model.R @@ -1,7 +1,7 @@ #' Fit a Neutral Model to Microbial Community Data #' #' The function fits the neutral distribution model developed by Sloan et al. -#' 2006, and implemented in R by Burns et al. 2015, to an OTU/ASV table and +#' (2006), and implemented in R by Burns et al. (2015), to an OTU/ASV table and #' returns several goodness of fit statistics alongside a data.frame #' with predicted occurrence frequencies for each OTU/ASV based on their #' abundance in the metacommunity for plotting the abundance-occupancy @@ -25,36 +25,26 @@ #' joined to per-taxon predictions #' #' @references -#' Sloan WT, Lunn M, Woodcock S, Head IM, Nee S, Curtis TP. (2006) Quantifying -#' the roles of immigration and chance in shaping prokaryote community -#' structure. Environ Microbiol. 8(4):732-40. -#' doi: https://doi.org/10.1111/j.1462-2920.2005.00956.x +#' Sloan, W. T., Lunn, M., Woodcock, S., Head, I. M., Nee, S., & Curtis, T. P. +#' (2006). Quantifying the roles of immigration and chance in shaping +#' prokaryote community structure. Environmental Microbiology, 8(4), 732–740. +#' #' -#' Burns AR, Stephens WZ, Stagaman K, Wong S, Rawls JF, Guillemin K, Bohannan -#' BJ. (2015) Contribution of neutral processes to the assembly of gut microbial -#' communities in the zebrafish over host development. ISME J. -#' 10(3):655-64. doi: https://doi.org/10.1038/ismej.2015.142 +#' Burns, A. R., Stephens, W. Z., Stagaman, K., Wong, S., Rawls, J. F., +#' Guillemin, K., & Bohannan, B. J. M. (2016). Contribution of neutral +#' processes to the assembly of gut microbial communities in the zebrafish over +#' host development. The ISME Journal, 10(3), 655–664. +#' #' #' Shade A, Stopnisek N (2019) Abundance-occupancy distributions to prioritize #' plant core microbiome membership. Current Opinion in Microbiology, 49:50-58 -#' doi: https://doi.org/10.1016/j.mib.2019.09.008 +#' #' -#' @seealso \code{\link{plot_neutral_model}} +#' @seealso [plot_neutral_model()] #' #' @examples -#' \donttest{ -#' library(phyloseq) #' library(BRCore) -#' # Example using your switchgrass phyloseq object and grouping variable -#' # 'sampling_date' -#' data("switchgrass", package = "BRCore") -#' -#' switchgrass_core <- identify_core( -#' physeq_obj = switchgrass, -#' priority_var = "sampling_date", -#' increase_value = 0.02, -#' seed = 092825 -#' ) +#' data("switchgrass_core", package = "BRCore") #' #' switchgrass_core_fit <- fit_neutral_model( #' otu_table = switchgrass_core$otu_table, @@ -64,13 +54,14 @@ #' #' str(switchgrass_core_fit) #' switchgrass_core_fit$goodness_of_fit -#' switchgrass_core_fit$model_prediction |> head() -#' } +#' switchgrass_core_fit$model_prediction |> +#' head() #' #' @importFrom stats confint pbeta pbinom ppois dnorm AIC BIC coef #' @importFrom stats4 mle #' @importFrom minpack.lm nlsLM #' @importFrom Hmisc binconf +#' #' @export fit_neutral_model <- function( otu_table, diff --git a/R/identify_core.R b/R/identify_core.R index c04e1fb..59fa9b6 100644 --- a/R/identify_core.R +++ b/R/identify_core.R @@ -19,7 +19,8 @@ #' @param priority_var The column name in the `sample_data` (e.g. #' sampling_date", "site") that is used for prioritizing the core microbiome. #' @param increase_value Increase value (numeric, scalar) used in the -#' calculation (default 0.02) for "increase". The "elbow" is always calculated and returned as \code{elbow_core} (see below for details). +#' calculation (default 0.02) for "increase". The "elbow" is always calculated +#' and returned as \code{elbow_core} (see below for details). #' @param abundance_weight Numeric in `[0,1]`; how much to weight mean relative #' abundance in the ranking score. `0` (default) uses occupancy/composite only. #' `1` ranks purely by abundance. Values in between blend the two (e.g., @@ -87,29 +88,23 @@ #' @references Shade A, Stopnisek N (2019) Abundance-occupancy #' distributions to prioritize plant core microbiome membership. Current #' Opinion in Microbiology, 49:50-58 -#' doi:https://doi.org/10.1016/j.mib.2019.09.008 +#' #' #' @section Dependencies: #' Requires \pkg{phyloseq}, \pkg{dplyr}, \pkg{tidyr}, \pkg{tibble}, \pkg{rlang}, #' and \pkg{vegan}. #' -#' @importFrom phyloseq sample_sums taxa_are_rows otu_table sample_data tax_table -#' @importFrom dplyr left_join group_by summarise transmute arrange desc mutate n last select slice_head -#' @importFrom tidyr pivot_longer -#' @importFrom tibble rownames_to_column column_to_rownames -#' @importFrom rlang ensym as_name .data -#' @importFrom vegan decostand -#' @importFrom cli cli_text cli_warn cli_abort cli_alert_success cli_alert_info -#' @importFrom utils combn tail +#' @seealso [multi_rarefy()], [plot_identified_core()], and +#' [plot_core_distribution()] #' #' @examples #' \donttest{ -#' library(phyloseq) #' library(BRCore) -#' # Example using your switchgrass phyloseq object and grouping variable -#' # 'sampling_date' +#' #' data("switchgrass", package = "BRCore") +#' data("bcse", package = "BRCore") #' +#' # With rarefied data #' res <- identify_core( #' physeq_obj = switchgrass, #' priority_var = "sampling_date", @@ -117,10 +112,39 @@ #' seed = 091825 #' ) #' -#' # Inspect results #' str(res) +#' +#' # With unrarefied data (requires multi_rarefy step) +#' rarefied_list <- multi_rarefy( +#' physeq_obj = bcse, +#' depth_level = 1000, +#' num_iter = 3, +#' .as = "list", +#' set_seed = 7642 +#' ) +#' +#' +#' res_rare <- identify_core( +#' physeq_obj = bcse, +#' rarefied_list = rarefied_list, +#' priority_var = "Crop", +#' increase_value = 0.02, +#' seed = 091825 +#' ) +#' +#' str(res_rare) #' } #' +#' @importFrom phyloseq sample_sums taxa_are_rows otu_table sample_data +#' @importFrom phyloseq tax_table +#' @importFrom dplyr left_join group_by summarise transmute arrange desc mutate +#' @importFrom dplyr n last select slice_head +#' @importFrom tidyr pivot_longer +#' @importFrom tibble rownames_to_column column_to_rownames +#' @importFrom rlang ensym as_name .data +#' @importFrom vegan decostand +#' @importFrom cli cli_text cli_warn cli_abort cli_alert_success cli_alert_info +#' @importFrom utils combn tail #' @export identify_core <- function( physeq_obj, @@ -144,7 +168,8 @@ identify_core <- function( .phyloseq_class_check(physeq_obj) # Define arguments ---- - # Check if samples are rarefied (all have same depth, accounting for floating-point precision) + # Check if samples are rarefied (all have same depth, accounting for + # floating-point precision) ## Validate / generate rarefied_list ---- min_sum <- min(sample_sums(physeq_obj)) max_sum <- max(sample_sums(physeq_obj)) diff --git a/R/multi_rarefy.R b/R/multi_rarefy.R index b3bbcf4..eef9e4a 100644 --- a/R/multi_rarefy.R +++ b/R/multi_rarefy.R @@ -1,9 +1,9 @@ -#' Run multiple rarefaction for microbiome count tables +#' Run rarefaction for microbiome count tables #' -#' This function performs multiple rarefaction on a `phyloseq` object by randomly -#' sub-sampling OTUs/ASVs within samples without replacement. The process is -#' repeated for a specified number of iterations, and the results are averaged. -#' Samples with fewer OTUs/ASVs than the specified `depth_level` are discarded. +#' This function performs rarefaction on a `phyloseq` object by randomly +#' sub-sampling OTUs/ASVs within samples without replacement for a number of +#' iterations specified by the user. Samples with fewer OTUs/ASVs than the +#' specified `depth_level` are discarded. #' #' @param physeq_obj A `phyloseq` object containing an OTU/ASV table. #' @param depth_level An integer specifying the sequencing depth (number of @@ -22,21 +22,11 @@ #' represent the average sequence counts calculated across all iterations. #' Samples with less than `depth_level` sequences are discarded. #' -#' -#' @importFrom parallelly availableCores -#' @importFrom parallel makeCluster stopCluster clusterExport parLapply -#' @importFrom parallel clusterEvalQ -#' @importFrom dplyr group_by summarise across everything filter near -#' @importFrom dplyr where -#' @importFrom tibble rownames_to_column column_to_rownames -#' @importFrom phyloseq otu_table -#' @importFrom cli cli_h1 cli_h2 cli_alert_info cli_alert_warning cli_alert_success cli_alert_danger -#' @importFrom utils head -#' @importFrom vegan rrarefy +#' @seealso [update_otu_table()] for updating the OTU table in a `phyloseq` +#' object and [vegan::rrarefy()] for the underlying rarefaction method used in +#' this function. #' #' @examples -#' \donttest{ -#' library(phyloseq) #' library(BRCore) #' data("bcse", package = "BRCore") #' @@ -44,13 +34,22 @@ #' otu_table_rare <- multi_rarefy( #' physeq_obj = bcse, #' depth_level = 1000, -#' num_iter = 100, +#' num_iter = 10, #' .as = "list", #' set_seed = 7642 #' ) #' #' rowSums(otu_table_rare[[1]]) -#' } +#' +#' +#' @importFrom dplyr group_by summarise across everything filter near +#' @importFrom dplyr where +#' @importFrom tibble rownames_to_column column_to_rownames +#' @importFrom phyloseq otu_table +#' @importFrom cli cli_h1 cli_h2 cli_alert_info cli_alert_warning +#' @importFrom cli cli_alert_success cli_alert_danger +#' @importFrom utils head +#' @importFrom vegan rrarefy #' #' @export multi_rarefy <- function( @@ -61,7 +60,7 @@ multi_rarefy <- function( set_seed = NULL ) { # Input validation ---- - cli::cli_h1("Multiple Rarefaction") + cli::cli_h1("Rarefaction iterations starting...") cli::cli_h2("Input Validation") if (!requireNamespace("phyloseq", quietly = TRUE)) { diff --git a/R/parallel_helpers.R b/R/parallel_helpers.R deleted file mode 100644 index 38997ea..0000000 --- a/R/parallel_helpers.R +++ /dev/null @@ -1,189 +0,0 @@ -#' Detect HPC environment type -#' -#' Detects if code is running in an HPC environment and identifies which scheduler. -#' -#' @return Character string identifying the HPC environment type: -#' "slurm", "pbs", "sge", "lsf", "condor", "generic_hpc", or "local" -#' -#' @examples -#' detect_hpc_type() -#' -#' @export -detect_hpc_type <- function() { - # Map of environment variables to HPC types - hpc_type_vars <- list( - slurm = c("SLURM_JOB_ID", "SLURM_CPUS_PER_TASK", "SLURM_JOB_NAME"), - pbs = c("PBS_JOBID", "PBS_NUM_PPN", "PBS_O_WORKDIR"), - sge = c("SGE_TASK_ID", "NSLOTS", "SGE_ROOT"), - lsf = c("LSF_JOBID", "LSB_DJOB_NUMPROC", "LSB_JOBID"), - condor = c("CONDOR_JOB_ID", "_CONDOR_JOB_AD", "CONDOR_SLOT") - ) - - # Generic HPC variables (not specific to a scheduler) - generic_vars <- c("OMP_NUM_THREADS", "JOB_ID", "HOSTNAME") - - # Check each HPC type - for (type in names(hpc_type_vars)) { - if ( - any(vapply( - hpc_type_vars[[type]], - function(var) { - val <- Sys.getenv(var, "") - !identical(val, "") - }, - FUN.VALUE = logical(1) - )) - ) { - return(type) - } - } - - # Check generic HPC indicators - if ( - any(vapply( - generic_vars, - function(var) { - val <- Sys.getenv(var, "") - !identical(val, "") - }, - FUN.VALUE = logical(1) - )) - ) { - return("generic_hpc") - } - - # Not in an HPC environment - return("local") -} - -#' Check if running in an HPC environment -#' -#' @return Logical indicating whether we're in any HPC environment -#' @examples -#' is_hpc_environment() -#' -#' @export -is_hpc_environment <- function() { - !identical(detect_hpc_type(), "local") -} - -#' Get the number of available CPU cores based on environment -#' -#' Determines the appropriate number of cores to use based on whether -#' we're in an HPC environment or a local machine. -#' -#' @param default Default number of cores to use if detection fails -#' @return Integer number of cores available -#' @examples -#' get_available_cores() -#' get_available_cores(default = 2) -#' -#' @importFrom parallelly availableCores makeClusterPSOCK -#' @importFrom parallel makeCluster stopCluster clusterExport parLapply -#' -#' @export -get_available_cores <- function(default = 1) { - if (is_hpc_environment()) { - # Check HPC-specific environment variables in a prioritized order - hpc_core_vars <- c( - "SLURM_CPUS_PER_TASK", # SLURM - "PBS_NUM_PPN", # PBS/TORQUE - "NSLOTS", # SGE - "LSB_DJOB_NUMPROC", # LSF - "OMP_NUM_THREADS" # OpenMP (generic) - ) - - # Try each variable in order - for (var in hpc_core_vars) { - core_str <- Sys.getenv(var, "") - if (core_str != "") { - cores <- suppressWarnings(as.integer(core_str) - 1L) # Keeping 1 core available - if (!is.na(cores) && cores > 0) { - cli::cli_alert_info( - "HPC environment detected. Using {cores} core(s)/worker(s) from {var}. - Keeping 1 core/worker available for system." - ) - return(cores) - } - } - } - - # If no specific HPC variable found, try parallelly if available - if (requireNamespace("parallelly", quietly = TRUE)) { - cores <- parallelly::availableCores() - 1L - cli::cli_alert_info( - "HPC environment detected. Using {cores} core(s)/worker(s) via parallelly::availableCores(). - Keeping 1 core/worker available for system." - ) - return(cores) - } - - # Last resort for HPC: use a safe default - cli::cli_alert_warning( - "HPC environment detected but could not determine core count. Using {default} core(s)." - ) - return(default) - } else { - # Local machine: use parallelly if available, otherwise safer detectCores - if (requireNamespace("parallelly", quietly = TRUE)) { - cores <- parallelly::availableCores() - 1L - cli::cli_alert_info( - "Local environment. Using {cores} core(s)/worker(s) via parallelly::availableCores(). - Keeping 1 core/worker available for system." - ) - return(cores) - } else { - # Safe local default: leave one core free for system - cores <- max(1, parallel::detectCores() - 1L) - if (is.na(cores)) { - cores <- default - } - cli::cli_alert_info( - "Local environment. Using {cores} core(s)/worker(s) via parallel::detectCores(). - Keeping 1 core/worker available for system." - ) - return(cores) - } - } -} - -#' Set up an appropriate parallel backend -#' -#' Creates a parallel cluster if in an HPC environment or returns -#' the number of cores to use if on a local machine. -#' -#' @param default Default number of cores to use if detection fails -#' @return Either a cluster object (HPC) or integer number of cores (local) -#' @examples -#' setup_parallel_backend() -#' -#' @importFrom parallelly availableCores makeClusterPSOCK -#' @importFrom parallel makeCluster stopCluster clusterExport parLapply -#' @export -setup_parallel_backend <- function(default = 1) { - if (is_hpc_environment()) { - if (requireNamespace("parallelly", quietly = TRUE)) { - # Create a PSOCK cluster with available workers - workers <- parallelly::availableWorkers() - cli::cli_alert_info( - "Creating HPC cluster with {length(workers) - 1L} core(s)/worker(s) via parallelly::makeClusterPSOCK(). - Keeping 1 core/worker available for system." - ) - cl <- parallelly::makeClusterPSOCK(length(workers) - 1L, autoStop = TRUE) # Keep one for other processes - return(cl) - } else { - # Fallback if parallelly is not available - cores <- get_available_cores(default) - cli::cli_alert_info( - "Creating HPC cluster with {cores - 1L} core(s)/worker(s) using parallel::makeCluster(). - Keeping 1 core/worker available for system." - ) - cl <- parallel::makeCluster(cores - 1L) - return(cl) - } - } else { - # Just return the number of cores for local processing - cores <- get_available_cores(default) - return(cores) - } -} diff --git a/R/plot_abundance_occupancy.R b/R/plot_abundance_occupancy.R index e4ba022..cdd8ef2 100644 --- a/R/plot_abundance_occupancy.R +++ b/R/plot_abundance_occupancy.R @@ -26,29 +26,18 @@ #' highlighted in red and non-core taxa in grey. The x-axis shows log10-transformed #' mean abundance and the y-axis shows occupancy (0-1). #' -#' @seealso \code{\link{plot_core_distribution}} and \code{\link{identify_core}} +#' @seealso [plot_core_distribution()] and [identify_core()] #' #' @examples -#' \donttest{ -#' library(phyloseq) #' library(BRCore) -#' # Generate an object from `identify_core()` and then plot #' -#' data("switchgrass", package = "BRCore") +#' data("switchgrass_core", package = "BRCore") #' -#' switchgrass_core <- identify_core( -#' physeq_obj = switchgrass, -#' priority_var = "sampling_date", -#' increase_value = 0.02, -#' abundance_weight = 0, -#' seed = 1234 -#' ) -#' -#' plot_abundance_occupancy( +#' p <- plot_abundance_occupancy( #' core_result = switchgrass_core, #' core_set = "increase" #' ) -#' } +#' print(p) #' #' @importFrom dplyr mutate #' @importFrom ggplot2 ggplot aes geom_point scale_fill_manual theme_classic diff --git a/R/plot_core_distribution.R b/R/plot_core_distribution.R index 4baabf8..ed7e9cb 100644 --- a/R/plot_core_distribution.R +++ b/R/plot_core_distribution.R @@ -1,6 +1,6 @@ #' Plot Core Taxa Occupancy Across Metadata Groups #' -#' Creates a bar plot showing core taxa (i.e. OTUs/ASVs) occupancy patterns across +#' Creates a plot showing core taxa (i.e. OTUs/ASVs) occupancy patterns across #' a grouping variable. #' #' @param core_result A list object returned by \code{\link{identify_core}}, @@ -27,28 +27,20 @@ #' #' @return A ggplot2 object that can be further customized. #' +#' @seealso [identify_core()], [plot_abundance_occupancy()], and +#' [plot_identified_core()] +#' #' @examples -#' \donttest{ -#' library(phyloseq) #' library(BRCore) -#' # Generate an object from identify_core and then plot -#' data("switchgrass", package = "BRCore") -#' -#' switchgrass_core <- identify_core( -#' physeq_obj = switchgrass, -#' priority_var = "sampling_date", -#' increase_value = 0.02, -#' abundance_weight = 0, -#' seed = 1234 -#' ) +#' data("switchgrass_core", package = "BRCore") #' -#' plot_core_distribution( +#' p <- plot_core_distribution( #' core_result = switchgrass_core, #' core_set = "increase", #' group_var = "sampling_date", #' plot_type = "bar" #' ) -#' } +#' print(p) #' #' @importFrom tibble rownames_to_column #' @importFrom tidyr pivot_longer @@ -57,7 +49,6 @@ #' @importFrom ggplot2 theme_classic theme element_text element_line #' @importFrom ggplot2 element_blank labs geom_line geom_point facet_wrap #' @importFrom ggplot2 scale_y_continuous geom_tile -#' @importFrom ggsci scale_fill_npg #' @importFrom vegan decostand #' @importFrom grid unit #' diff --git a/R/plot_identified_core.R b/R/plot_identified_core.R index 9d56fc4..eb7990b 100644 --- a/R/plot_identified_core.R +++ b/R/plot_identified_core.R @@ -25,34 +25,25 @@ #' @seealso \link[=identify_core]{identify_core()} #' #' @examples -#' \donttest{ -#' library(phyloseq) #' library(BRCore) -#' # Example with the package switchgrass dataset -#' data("switchgrass", package = "BRCore") +#' data("switchgrass_core", package = "BRCore") #' -#' # Identify core taxa -#' res <- identify_core( -#' physeq_obj = switchgrass, -#' priority_var = "sampling_date", -#' increase_value = 0.02, -#' seed = 48821 +#' p <- plot_identified_core( +#' bray_curtis_ranked = switchgrass_core$bray_curtis_ranked, +#' elbow = switchgrass_core$elbow, +#' lastCall = switchgrass_core$bc_increase, +#' increase_value = switchgrass_core$increase_value #' ) -#' -#' # Plot using the returned curve and cut indices; label from increase_value -#' plot_identified_core( -#' bray_curtis_ranked = res$bray_curtis_ranked, -#' elbow = res$elbow, -#' lastCall = res$bc_increase, -#' increase_value = res$increase_value -#' ) -#' } +#' print(p) #' #' @seealso \code{\link{identify_core}} #' #' @importFrom grid unit #' @importFrom dplyr mutate filter -#' @importFrom ggplot2 ggplot aes geom_point geom_vline annotate labs theme_classic theme element_text scale_x_continuous scale_y_continuous expansion +#' @importFrom ggplot2 ggplot aes geom_point geom_vline annotate labs +#' @importFrom ggplot2 theme_classic theme element_text scale_x_continuous +#' @importFrom ggplot2 scale_y_continuous expansion +#' #' @export plot_identified_core <- function( bray_curtis_ranked, @@ -176,13 +167,6 @@ plot_identified_core <- function( .brcore_theme() + scale_x_continuous(breaks = seq(0, max(new_data$rank_num), by = 5)) - # To add some more fancy graphics later... - # axis.title = element_text(size = 11), - # axis.text = element_text(size = 9), - # panel.grid.minor = element_blank(), - # panel.grid.major.x = element_line(color = "gray90", size = 0.3), - # panel.grid.major.y = element_line(color = "gray90", size = 0.3)) - return(list( df_for_plot = new_data, plot_identified_core = p diff --git a/R/plot_neutral_model.R b/R/plot_neutral_model.R index 3405d61..0c01c00 100644 --- a/R/plot_neutral_model.R +++ b/R/plot_neutral_model.R @@ -1,8 +1,9 @@ #' Plot a fitted Neutral Model to Microbial Community Data #' #' This function plots ASV/OTUs log abundances into a fitted neutral model of -#' microbial abundance-occupancy distribution. ASV/OTUs that are *Core*, *As predicted* (i.e. Neutral), -#' *Below* and *Above* model predictions are drawn with distinct point colors, see details. +#' microbial abundance-occupancy distribution. ASV/OTUs that are *Core*, *As +#' predicted* (i.e. Neutral), *Below* and *Above* model predictions are drawn +#' with distinct point colors, see details. #' #' @param fit_result A list-like object returned by \link[=fit_neutral_model]{fit_neutral_model()}. #' @@ -32,14 +33,17 @@ #' predictions and metrics. #' #' @examples -#' \dontrun{ -#' # Assuming you have run a neutral model fit on the switchgrass data.frames -#' is.data.frame(switchgrass_core_fit$model_prediction) -#' is.data.frame(switchgrass_core_fit$goodness_of_fit) +#' library(BRCore) +#' data("switchgrass_core", package = "BRCore") +#' +#' switchgrass_core_fit <- fit_neutral_model( +#' otu_table = switchgrass_core$otu_table, +#' core_set = switchgrass_core$increase_core, +#' abundance_occupancy = switchgrass_core$abundance_occupancy +#' ) #' #' p <- plot_neutral_model(switchgrass_core_fit) #' print(p) -#' } #' #' @seealso \code{\link{fit_neutral_model}} #' diff --git a/R/plot_rarefaction_metrics.R b/R/plot_rarefaction_metrics.R index 317a646..bb1f95a 100644 --- a/R/plot_rarefaction_metrics.R +++ b/R/plot_rarefaction_metrics.R @@ -9,8 +9,9 @@ #' #' @return A `ggarrange` object with six plots. #' +#' @seealso [add_rarefaction_metrics()] +#' #' @examples -#' \donttest{ #' library(phyloseq) #' library(BRCore) #' @@ -26,7 +27,6 @@ #' # pre-computed read_num, goods_cov, and outlier columns #' sample_data_df <- data.frame(sample_data(bcse_metrics)) #' plot_rarefaction_metrics(sample_data_df) -#' } #' #' @importFrom ggplot2 ggplot aes geom_histogram geom_point geom_jitter #' @importFrom ggplot2 geom_boxplot geom_bar theme_bw theme labs @@ -39,7 +39,6 @@ #' @importFrom scales comma #' @importFrom cli cli_alert_info cli_progress_step cli_alert_success cli_abort #' @importFrom phyloseq sample_data - #' #' @export #' diff --git a/R/plot_variance_propagation.R b/R/plot_variance_propagation.R index 0a47d66..f55b23b 100644 --- a/R/plot_variance_propagation.R +++ b/R/plot_variance_propagation.R @@ -1,6 +1,6 @@ #' Variance propagation diagnostic for rarefaction #' -#' This function evaluate the variance generated during multiple rarefaction by plotting comparing raw vs. rarefied alpha diversity metrics calculated at each iterations. It is possible to plot observed richness (q=0), Shannon diversity (q=1), or Simpson diversity (q=2) by setting the `q` parameter to "richness" or `q` = 0, "shannon" or `q` = 1, or "shannon" or `q` = 2. The plot is faceted by method (raw vs rarefied) and colored by a specified grouping variable from the sample data. +#' This function evaluate the variance between rarefaction iterations from \code{multi_rarefy()} by visually comparing raw vs. rarefied alpha diversity metrics calculated at each iterations. It is possible to plot observed richness (q=0), Shannon diversity (q=1), or Simpson diversity (q=2) by setting the `q` parameter to "richness" or `q` = 0, "shannon" or `q` = 1, or "shannon" or `q` = 2. The plot is faceted by method (raw vs rarefied) and colored by a specified grouping variable from the sample data. #' #' @param physeq_obj Raw phyloseq object #' @param rarefied Output from multi_rarefy(). Either a list of dataframes or and array. @@ -10,31 +10,28 @@ #' @param convert_to_factor Logical. If \code{TRUE}, both \code{group_var} and #' \code{group_color} are coerced to \code{factor} before plotting, which is #' useful when those columns are numeric/continuous (e.g. dates, counts) but -#' should be treated as discrete groups. When \code{TRUE} a discrete colour +#' should be treated as discrete groups. When \code{TRUE} a discrete color #' scale (\code{scale_color_viridis_d}) is used; otherwise the continuous #' scale (\code{scale_color_viridis_c}) is used. Default \code{FALSE}. #' -#' @importFrom tibble rownames_to_column column_to_rownames -#' @importFrom phyloseq otu_table sample_data taxa_are_rows -#' @importFrom cli cli_h1 cli_h2 cli_alert_info cli_alert_warning cli_alert_success cli_alert_danger -#' @importFrom utils head -#' @importFrom vegan rrarefy -#' @importFrom dplyr left_join all_of + #' #' @return ggplot object comparing raw vs rarefied diversity distributions across iterations. #' +#' #' @examples -#' \donttest{ #' library(phyloseq) #' library(BRCore) -#' # Example using the bcse dataset, comparing hill q=1 between Poplar and Switchgrass plots +#' # Example comparing hill q=1 between Poplar and Switchgrass plots +#' data("bcse", package = "BRCore") #' bcse_filt <- bcse |> #' subset_samples(Crop %in% c("Poplar", "Switchgrass")) +#' #' bcse_rarefied_otutable_filt <- #' multi_rarefy( #' physeq_obj = bcse_filt, #' depth_level = 1000, -#' num_iter = 100, +#' num_iter = 10, #' .as = "list", #' set_seed = 7643 #' ) @@ -46,7 +43,14 @@ #' group_var = "Crop", #' group_color = "Plot" #') -#' } +#' +#' @importFrom tibble rownames_to_column column_to_rownames +#' @importFrom phyloseq otu_table sample_data taxa_are_rows +#' @importFrom cli cli_h1 cli_h2 cli_alert_info cli_alert_warning +#' @importFrom cli cli_alert_success cli_alert_danger +#' @importFrom utils head +#' @importFrom vegan rrarefy +#' @importFrom dplyr left_join all_of #' #' @export plot_variance_propagation <- function( diff --git a/R/sncm.fit.R b/R/sncm.fit.R index 778ed1f..5c83810 100644 --- a/R/sncm.fit.R +++ b/R/sncm.fit.R @@ -1,60 +1,71 @@ #' Fit Sloan Neutral Community Model (SNCM) #' -#' This function fits the Sloan Neutral Community Model to community data to assess -#' neutral processes in community assembly. It compares the neutral model against -#' binomial and Poisson null models using various statistical measures. +#' The function implements the Sloan Neutral Community Model which +#' predicts the occurrence frequency of taxa based on their relative abundance +#' in the source pool and a migration parameter (\code{m}). It compares the +#' neutral model against binomial and Poisson null models using various +#' statistical measures. #' -#' @param spp A matrix or data frame where rows represent communities (samples) and -#' columns represent species/taxa. Values should be abundances or counts. +#' @param spp A matrix or data frame where rows represent communities (samples) +#' and columns represent species/taxa. Values should be abundances or counts. #' @param pool Optional. A matrix or data frame representing the source pool for -#' calculating relative abundances. If NULL, the source pool is calculated from -#' the spp matrix. Default is NULL. -#' @param stats Logical. If TRUE, returns fit statistics including AIC, BIC, R-squared, -#' and RMSE for model comparison. If FALSE, returns predicted vs observed frequencies -#' with confidence intervals. Default is TRUE. -#' @param taxon Optional. A data frame containing taxonomic information to merge with -#' results when stats=FALSE. Should have row names matching species names. Default is NULL. +#' calculating relative abundances. If NULL, the source pool is calculated +#' from the spp matrix. Default is \code{NULL}. +#' @param stats Logical. If \code{TRUE}, returns fit statistics including AIC, BIC, +#' R-squared, and RMSE for model comparison. If \code{FALSE}, returns predicted vs. +#' observed frequencies with confidence intervals. Default is \code{TRUE}. +#' @param taxon Optional. A data frame containing taxonomic information to +#' merge with results when \code{stats = FALSE}. Should have row names matching species +#' names. Default is \code{NULL}. #' -#' @return If stats=TRUE, returns a data frame with model fit statistics including: +#' @return If \code{stats = TRUE}, returns a data frame with model fit statistics including: #' \itemize{ -#' \item m: Migration rate parameter from NLS fit -#' \item m.ci: Confidence interval for m parameter -#' \item m.mle: Migration rate from maximum likelihood estimation -#' \item maxLL, binoLL, poisLL: Log-likelihood values for SNCM, binomial, and Poisson models -#' \item Rsqr, Rsqr.bino, Rsqr.pois: R-squared values for each model -#' \item RMSE, RMSE.bino, RMSE.pois: Root mean squared error for each model -#' \item AIC, BIC: Information criteria for model selection -#' \item N: Average number of individuals per community -#' \item Samples: Number of samples/communities -#' \item Richness: Number of taxa analyzed -#' \item Detect: Detection limit (1/N) +#' \item \code{m}: Migration rate parameter from NLS fit +#' \item \code{m.ci}: Confidence interval for m parameter +#' \item \code{m.mle}: Migration rate from maximum likelihood estimation +#' \item \code{maxLL}, \code{binoLL}, \code{poisLL}: Log-likelihood values for SNCM, binomial, and Poisson models +#' \item \code{Rsqr}, \code{Rsqr.bino}, \code{Rsqr.pois}: R-squared values for each model +#' \item \code{RMSE}, \code{RMSE.bino}, \code{RMSE.pois}: Root mean squared error for each model +#' \item \code{AIC}, \code{BIC}: Information criteria for model selection +#' \item \code{N}: Average number of individuals per community +#' \item \code{Samples}: Number of samples/communities +#' \item \code{Richness}: Number of taxa analyzed +#' \item \code{Detect}: Detection limit (1/N) #' } -#' If stats=FALSE, returns a data frame with observed and predicted frequencies -#' along with confidence intervals for visualization. +#' If \code{stats = FALSE}, returns a data frame with observed and predicted +#' frequencies along with confidence intervals for visualization. #' -#' @details The function implements the Sloan Neutral Community Model which predicts -#' the occurrence frequency of taxa based on their relative abundance in the source -#' pool and a migration parameter (m). The model assumes neutral processes govern -#' community assembly. Three models are compared: SNCM (beta distribution), -#' binomial null model, and Poisson null model. +#' @details The model assumes neutral processes govern community assembly. Three models are compared: SNCM (beta distribution), binomial null model, and Poisson null model. #' #' @references -#' Sloan, W.T. et al. (2006) Quantifying the roles of immigration and chance in -#' shaping prokaryote community structure. Environmental Microbiology, 8, 732-740. +#' Sloan WT, Lunn M, Woodcock S, Head IM, Nee S, Curtis TP. (2006) Quantifying +#' the roles of immigration and chance in shaping prokaryote community +#' structure. Environ Microbiol. 8(4):732-40. +#' #' #' @examples -#' \dontrun{ -#' # Generate example community data -#' spp_data <- matrix(rpois(100 * 50, 5), nrow = 100, ncol = 50) -#' colnames(spp_data) <- paste0("Species_", 1:50) +#' # Generate community data +#' set.seed(42) +#' n_samples <- 100 +#' n_species <- 80 #' -#' # Fit SNCM and get statistics -#' fit_stats <- sncm.fit(spp_data, stats = TRUE) -#' print(fit_stats) +#' # Log-normal abundances +#' mean_abund <- rlnorm(n_species, meanlog = 2, sdlog = 1.5) +#' # Simulate community matrix with multinomial sampling +#' spp_data <- t(sapply(seq_len(n_samples), function(i) { +#' rmultinom( +#' 1, +#' size = sample(500:2000, 1), +#' prob = mean_abund / sum(mean_abund) +#' )[, 1] +#'})) +#' colnames(spp_data) <- paste0("Species_", seq_len(n_species)) #' -#' # Get predictions for plotting +#' fit_stats <- sncm.fit(spp_data, stats = TRUE) #' predictions <- sncm.fit(spp_data, stats = FALSE) -#' } +#' +#' +#' @seealso [plot_neutral_model()], [fit_neutral_model()] #' #' @importFrom minpack.lm nlsLM #' @importFrom Hmisc binconf @@ -63,8 +74,6 @@ #' #' @export sncm.fit <- function(spp, pool = NULL, stats = TRUE, taxon = NULL) { - base::options(warn = -1) - # Calculate the number of individuals per community N <- base::mean(base::apply(spp, 1, sum)) N.int <- base::round(N) @@ -109,7 +118,7 @@ sncm.fit <- function(spp, pool = NULL, stats = TRUE, taxon = NULL) { sncm.LL <- function(m, sigma) { R <- freq - stats::pbeta(d, N * m * p, N * m * (1 - p), lower.tail = FALSE) - R <- stats::dnorm(R, 0, sigma) + R <- suppressWarnings(stats::dnorm(R, 0, sigma)) -base::sum(base::log(R)) } m.mle <- stats4::mle( @@ -147,7 +156,7 @@ sncm.fit <- function(spp, pool = NULL, stats = TRUE, taxon = NULL) { ## Calculate AIC for binomial model bino.LL <- function(mu, sigma) { R <- freq - stats::pbinom(d, N.int, p, lower.tail = FALSE) - R <- stats::dnorm(R, mu, sigma) + R <- suppressWarnings(stats::dnorm(R, mu, sigma)) -base::sum(base::log(R)) } bino.mle <- stats4::mle( @@ -179,7 +188,7 @@ sncm.fit <- function(spp, pool = NULL, stats = TRUE, taxon = NULL) { ## Calculate AIC for Poisson model pois.LL <- function(mu, sigma) { R <- freq - stats::ppois(d, N * p, lower.tail = FALSE) - R <- stats::dnorm(R, mu, sigma) + R <- suppressWarnings(stats::dnorm(R, mu, sigma)) -base::sum(base::log(R)) } pois.mle <- stats4::mle( diff --git a/R/update_otu_table.R b/R/update_otu_table.R index a667bf8..b87d681 100644 --- a/R/update_otu_table.R +++ b/R/update_otu_table.R @@ -16,20 +16,14 @@ #' #' @return A `phyloseq` object. #' -#' @importFrom phyloseq otu_table sample_data tax_table phy_tree refseq -#' @importFrom phyloseq phyloseq sample_names prune_taxa prune_samples -#' @importFrom phyloseq sample_sums taxa_sums -#' @importFrom cli cli_inform cli_alert_success cli_alert_warning cli_alert_info -#' @importFrom cli cli_abort -#' @importFrom magrittr %>% +#' @seealso [multi_rarefy()] #' #' @examples -#' \donttest{ #' library(phyloseq) #' library(BRCore) #' data(GlobalPatterns, package = "phyloseq") #' -#' # List output (.as = "list") +#' # List output #' otu_list <- #' multi_rarefy( #' physeq_obj = GlobalPatterns, @@ -42,7 +36,7 @@ #' # Extract iteration 2 #' rarefied_gp <- update_otu_table(GlobalPatterns, otu_list, iteration = 2) #' -#' # Array output (.as = "array") +#' # Array output #' otu_array <- #' multi_rarefy( #' physeq_obj = GlobalPatterns, @@ -54,8 +48,13 @@ #' #' # Extract iteration 1 #' rarefied_gp2 <- update_otu_table(GlobalPatterns, otu_array, iteration = 1) -#' } #' +#' @importFrom phyloseq otu_table sample_data tax_table phy_tree refseq +#' @importFrom phyloseq phyloseq sample_names prune_taxa prune_samples +#' @importFrom phyloseq sample_sums taxa_sums +#' @importFrom cli cli_inform cli_alert_success cli_alert_warning cli_alert_info +#' @importFrom cli cli_abort +#' @importFrom magrittr %>% #' @export update_otu_table <- function(physeq_obj, rarefied_otus, iteration = NULL) { suppressMessages(.phyloseq_class_check(physeq_obj)) diff --git a/README.Rmd b/README.Rmd index 48ce392..c127331 100644 --- a/README.Rmd +++ b/README.Rmd @@ -16,7 +16,7 @@ knitr::opts_chunk$set( -[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental) [![](https://www.r-pkg.org/badges/version-last-release/BRCore)](https://www.r-pkg.org/badges/version-last-release/BRCore) [![](https://cranlogs.r-pkg.org/badges/grand-total/BRCore)](https://cran.r-project.org/package=BRCore) [![Codecov test coverage](https://codecov.io/gh/germs-lab/BRCore/graph/badge.svg)](https://app.codecov.io/gh/germs-lab/BRCore) [![R-CMD-check](https://github.com/germs-lab/BRCore/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/germs-lab/BRCore/actions/workflows/R-CMD-check.yaml) +[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental) [![](https://www.r-pkg.org/badges/version-last-release/BRCore)](https://cran.r-project.org/package=BRCore) [![](https://cranlogs.r-pkg.org/badges/grand-total/BRCore)](https://www.r-pkg.org/pkg/BRCore) [![Codecov test coverage](https://codecov.io/gh/germs-lab/BRCore/graph/badge.svg)](https://app.codecov.io/gh/germs-lab/BRCore) [![R-CMD-check](https://github.com/germs-lab/BRCore/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/germs-lab/BRCore/actions/workflows/R-CMD-check.yaml) @@ -26,10 +26,11 @@ BRCore provides a unified framework for identification and ecological interpreta BRCore provides tools for: -**Rarefaction analysis**: Calculate pre-rarefaction metrics and perform multiple rarefactions +**Rarefaction analysis**: Calculate pre-rarefaction metrics and perform rarefactions - `add_rarefaction_metrics()` - `multi_rarefy()` +- `plot_variance_propagation()` **Core microbiome identification**: Identify core microbial taxa using abundance-occupancy distributions @@ -39,6 +40,7 @@ BRCore provides tools for: - `fit_neutral_model()` - `plot_neutral_model()` +- `sncm.fit()` **Visualization**: Plot rarefaction diagnostics, abundance-occupancy curves, and core distributions @@ -77,7 +79,7 @@ data("bcse", package = "BRCore") # Add rarefaction metrics bcse_metrics <- add_rarefaction_metrics(data = bcse) -# Perform multiple rarefaction +# Perform rarefaction bcse_rarefied_list <- multi_rarefy( physeq_obj = bcse, depth_level = 1000, @@ -86,7 +88,10 @@ bcse_rarefied_list <- multi_rarefy( ) # Update phyloseq object with rarefied data -bcse_rare_single <- update_otu_table(physeq_obj = bcse, rarefied_otus = bcse_rarefied_list, iteration = 2) # Your preffered iteration can be used here +bcse_rare_single <- update_otu_table( + physeq_obj = bcse, + rarefied_otus = bcse_rarefied_list, + iteration = 2) # Your preffered iteration can be used here # Identify core microbiome diff --git a/README.md b/README.md index 86a9a9f..ec24332 100644 --- a/README.md +++ b/README.md @@ -7,8 +7,8 @@ [![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental) -[![](https://www.r-pkg.org/badges/version-last-release/BRCore)](https://www.r-pkg.org/badges/version-last-release/BRCore) -[![](https://cranlogs.r-pkg.org/badges/grand-total/BRCore)](https://cran.r-project.org/package=BRCore) +[![](https://www.r-pkg.org/badges/version-last-release/BRCore)](https://cran.r-project.org/package=BRCore) +[![](https://cranlogs.r-pkg.org/badges/grand-total/BRCore)](https://www.r-pkg.org/pkg/BRCore) [![Codecov test coverage](https://codecov.io/gh/germs-lab/BRCore/graph/badge.svg)](https://app.codecov.io/gh/germs-lab/BRCore) [![R-CMD-check](https://github.com/germs-lab/BRCore/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/germs-lab/BRCore/actions/workflows/R-CMD-check.yaml) @@ -24,10 +24,11 @@ robustness and reproducibility in microbiome data analysis. BRCore provides tools for: **Rarefaction analysis**: Calculate pre-rarefaction metrics and perform -multiple rarefactions +rarefactions - `add_rarefaction_metrics()` - `multi_rarefy()` +- `plot_variance_propagation()` **Core microbiome identification**: Identify core microbial taxa using abundance-occupancy distributions @@ -38,6 +39,7 @@ abundance-occupancy distributions - `fit_neutral_model()` - `plot_neutral_model()` +- `sncm.fit()` **Visualization**: Plot rarefaction diagnostics, abundance-occupancy curves, and core distributions @@ -78,7 +80,7 @@ data("bcse", package = "BRCore") # Add rarefaction metrics bcse_metrics <- add_rarefaction_metrics(data = bcse) -# Perform multiple rarefaction +# Perform rarefaction bcse_rarefied_list <- multi_rarefy( physeq_obj = bcse, depth_level = 1000, @@ -87,7 +89,10 @@ bcse_rarefied_list <- multi_rarefy( ) # Update phyloseq object with rarefied data -bcse_rare_single <- update_otu_table(physeq_obj = bcse, rarefied_otus = bcse_rarefied_list, iteration = 2) # Your preffered iteration can be used here +bcse_rare_single <- update_otu_table( + physeq_obj = bcse, + rarefied_otus = bcse_rarefied_list, + iteration = 2) # Your preffered iteration can be used here # Identify core microbiome diff --git a/_pkgdown.yml b/_pkgdown.yml index 912024e..710e0c9 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,4 +1,4 @@ -url: http://www.germslab.org/BRCore/ +url: https://www.germslab.org/BRCore/ template: bootstrap: 5 bootswatch: minty @@ -28,24 +28,22 @@ reference: - add_rarefaction_metrics - multi_rarefy - plot_rarefaction_metrics - - title: Parallel Computing Utilities - contents: - - detect_hpc_type - - get_available_cores - - is_hpc_environment - - setup_parallel_backend - title: OTU Table Utilities - contents: update_otu_table + contents: + - update_otu_table - title: Datasets contents: - bcse - bean - mimulus - switchgrass + - switchgrass_core authors: Bolívar Aponte Rolón: href: "https://www.bolaponte.com/" + Gian Maria Niccolò Benucci: + href: "https://github.com/Gian77" Adina Howe: href: "https://www.germslab.org/about/" Center for Advanced Bioenergy and Bioproducts Innovation: diff --git a/cran-comments.md b/cran-comments.md index 1c67c05..83ac0a7 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,6 +1,25 @@ -## R CMD check results (v2.0.1 resubmission) +## Resubmission notes (v2.0.4) -0 errors | 0 warnings | 2 notes +This is a resubmission. Previous submissions: + +- v2.0.3: Passed auto-checks: authors found error in `scnm.fit()` and retracted by maintainer. +- v2.0.2: Passed auto-checks; requested resubmission by manual inspection +- v2.0.1: Passed auto-checks; retracted by maintainer due to typos +- v2.0.0: Rejected by CRAN auto-check service + +**The following issues from v2.0.3 were resolved:** + +1. Adding `suppressWarnings()` around `stats::dnorm()` calls in `scnm.fit()` to handle warnings previously suppressed by `options(warn=-1)` used a debugging phase. + +2. **DOI** links in the documentation were updated to use the correct format (``). +3. Minor typos and formatting issues in the documentation were corrected. +4. Examples were unwrapped from `\dontrun` and `\donttest`. Only example for `identify_core()` remains wrapped due to long execution time (> 5s). + 6. An appropriate example for `scnm.fit()` was included. + + +### R CMD check results + +0 errors | 0 warnings | 3 notes ### Notes @@ -10,26 +29,14 @@ the local check environment (`/usr/lib/aspell-0.60/en_US`). This is a local environment issue and does not affect the package. -2. **HTML version of manual**: `V8` package unavailable in the local check +2. **Checking for future file timestamps**: Unable to verify current time. This is a + local environment issue and does not affect the package. + +3. **HTML version of manual**: `V8` package unavailable in the local check environment, skipping math rendering check. This is a local environment issue. -### Resubmission notes (v2.0.1) - -This is a resubmission of v2.0.0, which was rejected by the CRAN auto-check -service (win-builder, r-devel-windows-x86_64). The following issues were -identified and resolved: - -1. **URL in DESCRIPTION**: Changed `http://` to `https://` for the package - website URL. - -2. **Description field**: Reworded to not start with "This package" per CRAN - policy. - -3. **Possibly misspelled words** (`BRCore`, `Bioenergy`, `microbiome`, and `microbiomes`): These are intentional. `BRCore` is the package name and `Bioenergy` , `microbiome`, and `microbiomes` are domain-specific terms. - -4. **Invalid file URIs in README.md**: Relative links to `CONTRIBUTING.md` and - `CODE_OF_CONDUCT.md` replaced with full GitHub URLs. +--- ### win-builder notes @@ -39,5 +46,6 @@ declared imports (`phyloseq`, `Hmisc`, `knitr`) that were absent in the win-builder guest environment. The CRAN incoming checks on both Windows and Debian passed with only NOTEs, confirming these are not direct imports in package code. + + -* This is a new submission to CRAN. diff --git a/data/switchgrass_core.rda b/data/switchgrass_core.rda new file mode 100644 index 0000000000000000000000000000000000000000..a966935f28071c75db4270df0c0dcf9e054f6705 GIT binary patch literal 35964 zcmV(tKvQ&2UKVgRpfklJ z)Hd78#Oenx$cBnP&=P&%8luk|GfWgvE;$P|kp};b!HNZ{g?j>I=g4AIWbnoC)WFmR z=j^Up6P|ieX-DY`s8c>86?LY|v8TH*e3%1Sq^@tJj(4*4(SZv6b0$R#JA>}ALx|U9 zZzdfSOaz;Cx2^fLflT17<$3O}l}@b)b^NiOI%1656`$DW&xEExR((vw$nM`Bs8($d zOS5B!o=tXjpg@Z7wv>()v`=4Xl+J@AIWz_P1i8mIwD(&~L8VA?EFqXT`Y1N%UqgJE zU%Mp!5VHViG#Yh3@hgr1)y$#B?GPXAp!(FPwAU_p!#0E{?rpl9aG(lxI8H!?nq=Ov z2lDn^`?*K{XO2V1ZM&penAQKONEq2wPS%`b(k}nkatx39jZt&F|FSy^>1_f|bP_8h zL!VpO$QAYf@D73Rqfo>rG7@F%$|1hYlC(n|)uFqG*DTd6cvCWSwpySCn=(pRZ>B^m}r4GBQx1r0wrETd5s-pp@k)}G8^yOzOZWf*ZS{LrfLr&V; zf1${SSp&4EtFV4FjJ*6HYBBsy8`naJ$oINNo4(n62nU&r% z!N_lnu8o4z_erM}Wk*wpF>lB~LC2-V7_`h?y3im0GFV^s-rw-ItIxlZSg%%rV(BuJ zbM*FLoOCW!^xu#OY6$zbx;v#l_3wmVd~guiKzatsk_Q>ZQ3~bSL{JQ(QtLdul_g)# zL&(Yh{vZ{yG$H`+KZR#Q6QF;-;+GuBMd;RgxYUiUe`m^!jCto4tuL_SJX2dc(gWPN zLhKe65|$P88BUD5bZEt4`2kkgD;qYSduW#WtFK{WhWD zqGB)*YJN~aX2s*?2070KwLXsAe{fUUs9+$RWO}KcH~EH@aB)v#S`C#Tvg9<+jSvOZ z2v)_OpEHj>6yRi5aWwbjWWO2#SQAPpT?@}jGpC{Jqbcg!>R4{?az$c51TaIMNaH+D zS}hNwnD_LBv|9P2IslvM@YA;C1~T3`)B}r?OnJ&LpMRU4Ia&unU;aV}y39Qk0F^H_ zg{;y*A2fOGUO-p`GY}^v*}|yi8Od+~WdU2dC$Xq1lFW&@{1P`o6Q);vM#Z>4|L}}u z!PX(EdV?RA_%sT&Q|*4Z+Ipr!hD;0cUUx0D@xGA}sz`viN)ELKai6<-LTs$4=JJm0 zFjAUA(^scagE4QeQIJQ#3LVs$u--XcSX5F_;ST)gq^a zhM<`T&P7NRHrOBiH>-q)w`0ELcW*cNfad1($xsL za#6OVQjdWz=>_W#Az5Z!eT_Jl_W6dnn5%A$K;k=d0R2O6v+Q5PT^HF%FQ6^sceLKv z;*(qTb3c{X4{-F0FSbOe4Oe1E9U*we+zUL1HLadnVf&Sdis0JdNCqq+-UMy5udKJ$DMA=cA?o5BIEl1kEfjEFzmvCZyX_ zEeC4~@jsxColT5xo74scR{Kl6k^!aiiGNZk(rfyrV`1svd};y3fQHZ0Dd z+S^Al&{oU8Rb)Qr@oh#Y^bC)Axk{N*$hZnidu(;MUXMx6&cqIPRvoqPrRnFDLLG;K zF#|;%98pzkd3rtlg_EYzK5qRTaF9Y=cA?Fo1J>!r1}nzIWX;R7d{3aYJ>!Tq1>yWG~~0NaD|e&Ho4|rA$e^17?J2E{oGE9EArE&&gN%d2U=5ux6I1TtlAgV&|!i9 z$#2JC`e}!O-1ykK)U@_h7-OGV37g(ya1GwMAb1n%*Szj-I)QMYZ81F|qcub?v)kCT z+;J5Nf#Kw=lr_zQ#)w(<&sP$}=u*k(cCL2AJ|KQ#-xf{!y?eQ|yTwDTJ5c-(qg*SP zW}gKF&{!a=Ar?NpbY<${gc*S)$J~SE=XLm>DpL>^Edj<~+}+=6M*GK6L+{_zPcM;? zHtw80IZYgs5?mqlhyv#}R0%jqFPRD!vmk8W1xS!ty5~3IyqAG1letja{LkBhR(^9Q zGX3jvXAAa#^2uoSJNJC_>D|TkxikYyxqX%$ zf!L~-qf-pukRW|i&W<={7y}qN?H3+_g_+-x@<^`=kuYX8N*pAAVoDdO6uo$-m&nD{kRJcL$P9j?!nw@hT2P?>oZA51JRRW)I(PUUT%MU_p$`uMZSqcl#2ow|t74l?y1B=C!OY(ug zoY%NS4ta85!d=)x8M{d6MLfQwz4#gni=l?rQ!BbXJrEW%3ooen!Bn!bbRf<_;}?Mn zUF%D`x;8EZ4nvpf2f>vI99!7E%ywvcP3ZMtf>tMdOy+@Qm`%C8f}`RRQ6)96c*-ss zg6?|WwhaX;5{fPI5_v~2k)km-!^B}WF8+Q1w1+X!<6VhoC;N9;|6W?cqzC?tC;CsW zf2MPNR>8CROx|ss~7I z{UV{BVKD(-XEWQ|bswk_zAp!Bh9_EQL5(@m-x%;LWesa(3QvR?>3DKGP8LND?(exj z0qg#QlSQKaE!JYsiRG)JWY@;dAe}K|Kol-E>G2_}2M2Z38LIY@?|l!#M8l{}jrGp+ zGfOf@m#^sz-Wg1IIK!bp=G(Rc*hx4#jnMjYQ~9eq0S54Jc|{TmaiUSAWI_;z3qbI| zWUoaGdm*)@7;jM%ft(gunN&;9WUR-9ROMki=Bp5P-d9Cm+C2=ceu1&lqYc3n0_FQ{ z1qYIAqV}&>`w)F+jkFQ$zcb8L=--L41w@}hv-DvG2Db6(P-?Z(Y?3DUR$SZZ`KPug z_T7Kg!*w@SW{1ddkmtvX?f3F;R8QB=JM6H~C6Rn{!i;CNc>hq2y#lhn1a1r3fq~CDE8$;6rdJ?|_#s?6nNyFtJFXlc`zcZ{K#=wxGKDJOE11^G&-!y(SN z?SHl^%a!%x5t!jFmAdtVzfFKDnrOer!OPSdaEyPN`)~zxp1V_>eAryQX4gdZWxUAI zy%Gb=pxITYkJQo*gW9yLZz`M1TqgvX^F-MsQ*Rd&7Oc$$vlGcnk7b zRgP~3Bl`oXuz;DD^h7ywBbAJ?RnC)c3j7h^Tzf7nTsx9<<5XoRq>6vzErej@4TWQ_ zXpS8*a+kCXCk z^y2lU-FMQ`DAl04tV^|t06SG|NXh}^8)9iN6T<_+v9>Y4tgeI3Lc6hJy1#D%Y}}9# z!f|JRs!g`vF@V@aN!E$qH$N<2b<`xi$=v~tMt3WNYY4vBuXEVquY|Q@joiloS=feN zfyf(cStV6v@AZ&f^E=1BU7Wp&gQ2D(c)X#+Bsr{jNM!b%4v=Ff^gypV@x52kN)R}F z{n(ip`E>^YbWu%z!2Y>-)Qz1?eX5R#W6(>$KReSZM%%g+cMV=xxg`%qJ{KAtDAVZko;HW#hECKSgs@NaJ>kL zrc1erHG71gZ`H?osLYvjK3MRPzsx6R$0z_FwG6nuznq3g0;1jT_M-|*!4^6`ZjzSM ze=VBpS(6z9#s0Z1VHF72`9bN?e>!EXhC2!c$z3j~ zSE(4O%B;}+$xvD^VKb~E>D}OUg6L^h`J4LI>f z%oiA_H9yM^6wQe+gB7rBWfZIy&R!%6Xee?agOMeo^Hu@39W-@%EZA1m7Pd>UJG2+} z?r#f5Jo>h69?DzcC{hZ;UX=$``O00K$RKp^?a$yv8*!RG*!u1Q&K_63T+#|$s z0`gwf^-t^^aHmd6?&aSXC>@43N@q_Y%I|$F?UI+Vn<}s~kA=Xb_ zEbsyZk3esd`u$~fr%d3muV~ zfGfd8ver1iDTg~d+kZSIpXd)WR>GPyc#*DVLTA=2RJMRjf9AAk?Hi}eW2ilF>0G7G zc12v$xpFT)8~+{lcF@FByw{++0e)lPNmVdL5><&FGpIetK1enytYA@l^j_YfwFLDQHt^`1xQ@sjaX@3cXz7<<9Q=;bS+2o9NBlW;cYvp(M!7xWd0L$Ng)3;50+TH++ zup)OJxBdl+LuZNI6x}v?BwO%%7nF5p5tE3v(*&E0ghr|yYxc_r{8AhGo1JXmdj$ce zxz2XT?K==>%=P@7@;+!$-#_Xl7I#rbiU>_7ZXOP+Axr@1k7sCdu`~Vtt=4*fctc6~ zH+DZE+E-=UrZBq}0=Qi?RL?Y?E`lsqd=Jp&-xPt1bO2`R0fXxvfnK8S2SJqYiD$hE zVM|`nKfC`!ciSAwQ?k%cIvJV{B$t0Oy8!PVY_9)6_=7^NN%ejVPMDUiAZym2o0$YA zkqLM$tf9wYk&ZVA_qGq3Gckmo-AP+Ol8;{>{~G7al^!b)~x=fh@3^4`>O&{6@#=4mk{u+!xPvyF+T=$cDX zOzx+-V49;7`pj=psR9`|{A#>|Imaz0#6jeX+nM|E(|_n(J(^34L}L-1N(?QAvBQ(J z72sw43_D+dcr-JtCQC!-$U#Tf((zVAHPL=<&!hGg;BL_&{tUE+r`bJ85uL~DhH?}R z%&rB$Owu*R;)Mx2LRom9BZ$YS2Q!omXL?5F?l>4tuS!0ORXm~0-t6?K;P}$H_ttp$9S5g|Yl6Nn?Dsn8LLqIH962;;Cr!BHRw(I_q#M8gQ_=pOqa1iS~I-?3&gN4o~G(^2=jb@)dhq zv!xhE(cx}`Ko*UPS7hJd8Iz#x9bVUJ<`~4O1}PaW#^qnKMGIh8_MJ|@ejrud9IqU& zZ+*U(ZD9w~O3IvhQ{$1c+2)GQk10Im{b0)+j-;VOR6R9oZDJ60f9-Km+i?{*(6hm% zm+;`^+t9x+p`Yk(VGH#5%1eUscztUx^a^R?kzW~NGA=t$P0?+hv9%0kk1Zdwb;Fr3 z7P*xM+v|kH8_}jI8q4f#9<#gBlun@-4Ku=Vgf=e?F1?xW*R7W`im<8HOxW#l-hp>bw>Gr;S4&;SGrkwaRP6N=kV$_(&7LG20~jdJsX&)C_l^ zq)*Mm4j?L8jT669=&c&xjmmp=$0U7*h%^Dk1jc7&YLo5)?!Z}?$_)$BGa}PAQOtuj zo=C#>F=(70J7~M@Cgb?`E=LPX|Rr$PH{2t@D` zPRA4V(K&9SIzeFS3Rpg79~R$yxEO|=!jRH+9+nLak2Qm)vU?rBBNKec95|`;u-WBR zWs|{>Rm78!J@Id2tMjc=VemPj2*tB;h{y?Xj(SDgfZ@o7oAUStD7aP(WuN#yHk1zg zYx`=+25j24j^63%^X(aS4!?%{Wwe`WIZCzTv547@iSAbB-dwR%>I2-SC9*&uVd2^q zmQ9uuVdngA){OhQ68EW02`n-h(D2sDc3avIJV!h%LbSB+Cn2pSNx7X;45OFSxMI%B zmh-JSvwF9KrX8nYNeaKuRfzp`&%1)4dTHb6+B$H9zwgD7{G0(3a-b2%#*jXVc%Y33 zY^UDE_t}rc{UA80{_u87t7;a;SY6&EU5heLY|Edm*S2XAAXv)@U0+W?C3ivSegeFE zsxX>Byee#m@c;>+wr=Jk$tU%DS^Q-%F;>MO}75V#Ch4u zUhPSJJ+yi!SP6o%ou(IqQ*XyHr(>Q!%sOIf5yUaRoa|mjg%{S-AWS9kFl6RA1G(^> zdb5PcJx45LQr6QT$Sc^>#v#Qyvtum-ZXVtKlO`&|I&~x2X(XdOzvE8IT{E9uGRlkh(V;4v-gB%CI|6!jP(~0 zvJAsd@i60(_RT%8^}n|Vt==M=Vdk?0^LGeH0kv}QUu1l-FpIt=kZbY|Wn=PY5f;M= zS`M9mKI3bk)cC^M53mS~$i482JpR&R%uMfk(=mI_`4|!IolGO5`W81Tsuk zmDXj45#YxTk2YiSeZ&C9uk@uCFl3a%VcZpSEr^z~TIREL!%@D|(Ck;A7nn;4ybzO@ z;)xiZXhlp2M(p{eO4m2GZ`4By*;ezp#lX1m)|&N25>tr4t}AW!RhwFa!!>G#0nOzk z@I+=#!BdbNl$pxHfw28V++D;e{VY(B?b3Yf3vSZBaoxiY%+c>-N5wPym*|3mr|&)i ziO!cMXNJBkO{Dt-B5p~LRUWq|yQd|?ik~qQUmyl|aHf z9L}+p7XfPUCb7TnUvWq#xIqu6pg!dG8Culb)R@nLbo01GQT`#%bQQiyg8$Aoc*ys@@M*jdYffPm1I2UF;Z@Gz(I=`G0?N{}{-FFCvM14{<&N+^a@K ztkpHU$L~5}nvSyBtcbc>azG@?P0Z4w5D(Ki*ZZ%*ATf7KOAtkD?+-tH2e zPGpkH2d(!)Vf={!WKF|B)bwV4<+I=U3`mU09#2na;8Ltt%$rX`?<1@UJ4$HSVY42@ zr!^P+{oSGnU-Xf|c|tj4v)sA`%}n`c3js+M4Oz$Vk7z`PfjC@8b9@38YI0`+GS73j zVrhhQ_%>;R_P0(fL6Gbvi{}!^8xHPDzO5^9le4n2o0c~t0fPo$zsZRIa5e>4%}8!K zL;CvIJD#$^opbCZ29<53-lf{DH6IFC2*zeQb&XadP+s|hSSEu;ZfN0>e6iN(X;>CR zj8y0j{$N>#5CXA<{B{3eN^ol6BnnE#=F6fqCn8!4c{XiQ4cEd9V_n(IHjcFI+pY70 zu}XP)ukXog7&?fK;l7X%Wr!iV9H|kYo{2bUWB|G2Ra2na;DAFvrjj8|8neXg#jZ)5 zewT4sNP6>!6YRR!?=duQ<*8^DvbwL_&xJ2n>L?|}RUza7H{eBU2ePZXphz->T&Z#+ z0(xb|3^PkZSk|QYF9MZtExZy2kJjx}qhWK(k4Th^4Kqn!RqsU0=j$({aTv7fg$+Vp z;(x%;;Bp4~VrZ4Ptg<6`H6v=TE)*FQp`>{e3q3POEh!2*^ks0ydW42d;7YHK)2AWM9z+Pez%+Zdrro03&Q5;aX9UNVgS@1M?4`vNRVY`ILmNs_D1U}$(OkF zz1jA06g5n+g)UcH_GPiK>E9376XhJ;rH5F(CJU(6FiyBSS=GGUMtsDtr;BMPf_B80 z>wpR1Q)^28n;++00+WfBJYw|VScM8G!dDQ{wcnzNHT@YjdsQt5SnSe;m?^VH--Yp#MO;5@HR2{}Xll?v% zi6sv&hhu3L#$*2++{Kc32u^b9m$Tq1;_B;{U5}q*oXkc^GyF6t9appJd4HRnYw~zm zD?aWG41?5NYQR$f3qDy*iRe4UF#%S|n3n|uYVMQrJ=zU z*RDSiDx-w)(H!w7GvSt=x(+|7ta()Q2~&Y2Qbhu=QyQZB5QHr;($WcMd`zLghNQ%4;Wn)Y7Kna-K5&rml(>I%&=5Q+yP(WmGM6F@UAHVX{^d26g?1aY4#q*^_hG_N1K=9g8yj^C~ zrRTAQT{Go&$o?N0*&c>tm@_WL81wn*956LovZNu}DBA0#TgL_jZ@a#1d@uQQfgE@t zeZB8C%m>j<63A0v1{qy<&Z7=WG}K(7bYD%i=B^k)DNE>(_;Hrxl$jYJe&2+N?|@pi zin}dVEB+f6(*UVtjqcPR9jW~zJkG(%oQCF)gj3ld<2X^_nSK7Sqk4`b3WN`+g7h8@-&ssoL!@4lQHOI zWcp~-V`@84T%6xgyIy(5LbJep6Aac$`E4?876N!IK4 zFDkOaxmf(fRp_cEy0`!~{VbxIFtYUbEeeNLHI>B~kE!-O@h+-P66!I_Yjve@hew4T zMU+2AChj&xFixp?xt0m3nh`8(RWYuZ{0reE?(Oo-8|1!BqMeH*rxM+XgMd%YLrXY;gubu0wjG}7fw5Q z>$M+d_u0>?&uT|38-)BKxgIFMh~JL?L@KVW%P$j62~CTrhO+wPhX~<0WxiZjTz{|G zZ-X~L8w2hUD|BmP88B1ST2-;}`6xC^9X?E(S(>eidfNVkPeur@lde)