diff --git a/inst/rmd/bsdb_signatures.Rmd b/inst/rmd/bsdb_signatures.Rmd index 53cdbf4..5d9b809 100644 --- a/inst/rmd/bsdb_signatures.Rmd +++ b/inst/rmd/bsdb_signatures.Rmd @@ -3,6 +3,9 @@ title: "High-throughput ORA with BSDB signatures" output: html_document: toc: true +params: + bugphyzz_version: "devel" + uberon_obo: "https://github.com/obophenotype/uberon/releases/download/v2024-05-13/uberon-basic.obo" --- ```{r, include = FALSE} @@ -84,7 +87,12 @@ dim(bsdb) BugSigDB subsets by body site: ```{r, message=FALSE} -uberon <- getOntology(onto = "uberon") +## Uberon, pinned to the release current when these analyses were first run +## (May 2024). bugsigdbr::getOntology() fetches the latest release and its +## download step fails on some networks, so the OBO is read directly. +uberon <- ontologyIndex::get_ontology( + params$uberon_obo, propagate_relationships = c("is_a", "part_of") +) bsdb_subsets_by_bodysite <- vector("list", length(body_sites)) names(bsdb_subsets_by_bodysite) <- body_sites for (i in seq_along(bsdb_subsets_by_bodysite)) { @@ -102,7 +110,7 @@ for (i in seq_along(bsdb_subsets_by_bodysite)) { Bugphyzz: ```{r, message=FALSE} -bp <- importBugphyzz(v = 0.8) +bp <- importBugphyzz(version = params$bugphyzz_version) # manuscript thresholds (MCC 0.7, R2 0.5) are the defaults names(bp) ``` diff --git a/vignettes/articles/data_integration_and_stats.Rmd b/vignettes/articles/data_integration_and_stats.Rmd index e87673e..f508711 100644 --- a/vignettes/articles/data_integration_and_stats.Rmd +++ b/vignettes/articles/data_integration_and_stats.Rmd @@ -4,6 +4,10 @@ output: html_document: code_folding: hide toc: true +params: + bugphyzz_version: "devel" # Zenodo DOI of the release the paper cites, or "devel" + validation_tsv: "" # path or URL of a taxPProValidation validation_summary.tsv; "" uses the copy in the installed bugphyzz + out_dir: "" # where the supplementary TSVs are written; "" means tempdir() --- ```{r script options, include = FALSE} @@ -38,26 +42,26 @@ library(cowplot) ## Bugphyzz data ```{r import bp, message=FALSE} -discrete_types <- c( - "binary", "multistate-intersection", "multistate-union" -) -bp <- importBugphyzz(v = 0.5) |> - map(~ mutate(.x, NCBI_ID = as.character(NCBI_ID))) |> - map(~ { - attr_type <- unique(.x$Attribute_type) - if (attr_type %in% discrete_types) { - df <- .x |> - filter( - !(Validation < 0.7 & Evidence == "asr") - ) - } else if (attr_type == "numeric") { - df <- .x |> - filter( - !(Validation < 0.5 & Evidence == "asr") - ) - } - df - }) +## importBugphyzz() keeps imputed (ASR) annotations with MCC >= 0.7 +## (discrete) or R2 >= 0.5 (numeric) by default, the manuscript thresholds; +## since bugphyzzExports 2026-11 the exported files are already filtered. +bp <- importBugphyzz(version = params$bugphyzz_version) |> + map(~ mutate(.x, NCBI_ID = as.character(NCBI_ID))) +## Cross-validation table (taxPProValidation): the installed package's copy +## unless params$validation_tsv points to another version +validation <- if (nzchar(params$validation_tsv)) { + readr::read_tsv(params$validation_tsv, show_col_types = FALSE) +} else { + if (params$bugphyzz_version != "devel") { + warning( + "validation_tsv is empty, so Figure 2 and the attributes table use ", + "the validation_summary.tsv shipped with bugphyzz ", + packageVersion("bugphyzz"), ", which may not be the version the ", + params$bugphyzz_version, " release was filtered with." + ) + } + bugphyzz:::.validationData() # bugphyzz exports no accessor for this table +} length(bp) ``` ```{r} @@ -131,7 +135,7 @@ attr_asr_default_names <- map(bp, ~ { attr_asr_default <- length(attr_asr_default_names) -attr_asr_all <- bugphyzz:::.validationData() |> +attr_asr_all <- validation |> filter(rank == "all") |> select(physiology, mcc_mean, r2_mean) |> pull(physiology) |> @@ -148,7 +152,7 @@ count_summary_tbl <- tibble::tribble( "Number of annoations from sources", n_annotations_source, "Number of IBD annotations", n_annotations_ibd, "Number of annotations imputed through ASR", n_annotations_asr, - "Number of attributes with ASR data (default > 0.5 val)", attr_asr_default, + "Number of attributes with ASR data passing validation", attr_asr_default, "Number of attributes with ASR data (all)", attr_asr_all ) @@ -235,7 +239,7 @@ Create plot of MMC values for discrete attributes: ``` ```{r validation discrete} -val <- bugphyzz:::.validationData() |> +val <- validation |> filter(rank == "all") # Validation including taxa of mixed ranks as input (keyword is "all") dis_phys_names <- val |> filter(!is.na(mcc_mean)) |> @@ -597,6 +601,64 @@ sources <- read.table( myDataTable(sources, page_len = nrow(sources)) ``` +# Supplementary Table: attributes + +One row per attribute: type, unit, number of sources, annotations by evidence +(curated from a source, IBD = inherited from strains by taxonomy, ASR), taxa, +and the cross-validation score of each attribute value at rank "all" (MCC for +discrete, R2 for numeric; those at or above the export threshold, 0.7 and +0.5, have imputed annotations in the release). + +```{r supp table attributes} +val_str <- validation |> + filter(rank == "all") |> + mutate( + score = ifelse(is.na(mcc_mean), r2_mean, mcc_mean), + label = ifelse( + physiology == attribute, sprintf("%.2f", score), + sprintf("%s %.2f", attribute, score) + ) + ) |> + group_by(physiology) |> + summarise(`Validation (rank all)` = paste(label, collapse = "; "), .groups = "drop") |> + mutate(Attribute = tolower(physiology)) |> # the export lowercases attribute names + select(-physiology) +supp_attr_tbl <- bp |> + map(~ { + .x |> + summarise( + Type = unique(Attribute_type)[1], + Unit = if ("Unit" %in% names(.x)) paste(unique(na.omit(Unit)), collapse = ";") else "", + Sources = n_distinct(Attribute_source[!is.na(Attribute_source)]), + Curated = sum(Evidence %in% c("exp", "igc", "tas", "nas")), + IBD = sum(Evidence == "tax"), + ASR = sum(Evidence == "asr"), + Total = n(), + Taxa = n_distinct(NCBI_ID) + ) + }) |> + bind_rows(.id = "Attribute") |> + mutate(.key = tolower(Attribute)) |> + left_join(val_str, by = c(".key" = "Attribute")) |> + select(-.key) |> + arrange(desc(Total)) +## every annotation is curated, IBD or ASR, and every attribute with ASR +## annotations has a validation score to show for them +stopifnot( + all(supp_attr_tbl$Curated + supp_attr_tbl$IBD + supp_attr_tbl$ASR == supp_attr_tbl$Total), + !any(supp_attr_tbl$ASR > 0 & is.na(supp_attr_tbl$`Validation (rank all)`)) +) +myDataTable(supp_attr_tbl, page_len = nrow(supp_attr_tbl)) +``` + +```{r write supp tables, include = FALSE} +## TSV copies for the manuscript supplement +out_dir <- if (nzchar(params$out_dir)) params$out_dir else tempdir() +readr::write_tsv(supp_attr_tbl, file.path(out_dir, "supplementary_table_attributes.tsv")) +readr::write_tsv(sources, file.path(out_dir, "supplementary_table_sources.tsv")) +message("Supplementary tables written to ", out_dir) +``` + # Session information ```{r class.soruce = 'fold-show'} diff --git a/vignettes/articles/enrichment_dbbact.Rmd b/vignettes/articles/enrichment_dbbact.Rmd index 4f1ba72..77f7b26 100644 --- a/vignettes/articles/enrichment_dbbact.Rmd +++ b/vignettes/articles/enrichment_dbbact.Rmd @@ -3,6 +3,9 @@ title: "Enrichment with dbBact approach" output: html_document: toc: yes +params: + bugphyzz_version: "devel" # Zenodo DOI of the release the paper cites, or "devel" + uberon_obo: "https://github.com/obophenotype/uberon/releases/download/v2024-05-13/uberon-basic.obo" --- ```{r, include = FALSE} @@ -41,26 +44,11 @@ freq_var <- 1 Import bugphyzz and create signatures at the genus and species levels: ```{r, warning=FALSE, message=FALSE} -discrete_types <- c( - "binary", "multistate-intersection", "multistate-union" -) -bp <- importBugphyzz(v = 0.5) |> - map(~ mutate(.x, NCBI_ID = as.character(NCBI_ID))) |> - map(~ { - attr_type <- unique(.x$Attribute_type) - if (attr_type %in% discrete_types) { - df <- .x |> - filter( - !(Validation < 0.7 & Evidence == "asr") - ) - } else if (attr_type == "numeric") { - df <- .x |> - filter( - !(Validation < 0.5 & Evidence == "asr") - ) - } - df - }) +## importBugphyzz() keeps imputed (ASR) annotations with MCC >= 0.7 +## (discrete) or R2 >= 0.5 (numeric) by default, the manuscript thresholds; +## since bugphyzzExports 2026-11 the exported files are already filtered. +bp <- importBugphyzz(version = params$bugphyzz_version) |> + map(~ mutate(.x, NCBI_ID = as.character(NCBI_ID))) bpSigs_g <- map(bp, ~ { makeSignatures( dat = .x, taxIdType = "NCBI_ID", taxLevel = "genus", @@ -96,7 +84,12 @@ dim(bsdb) Subset by body site: ```{r} -uberon <- getOntology(onto = "uberon") +## Uberon, pinned to the release current when these analyses were first run +## (May 2024). bugsigdbr::getOntology() fetches the latest release and its +## download step fails on some networks, so the OBO is read directly. +uberon <- ontologyIndex::get_ontology( + params$uberon_obo, propagate_relationships = c("is_a", "part_of") +) bsdb_subsets_by_bodysite <- vector("list", length(body_sites)) names(bsdb_subsets_by_bodysite) <- body_sites for (i in seq_along(bsdb_subsets_by_bodysite)) { diff --git a/vignettes/articles/nychanes.Rmd b/vignettes/articles/nychanes.Rmd index 0ee9009..f1cf54a 100644 --- a/vignettes/articles/nychanes.Rmd +++ b/vignettes/articles/nychanes.Rmd @@ -4,6 +4,8 @@ output: html_document: toc: true +params: + bugphyzz_version: "devel" # Zenodo DOI of the release the paper cites, or "devel" --- ```{r, include = FALSE} @@ -89,26 +91,11 @@ tse_gn bugphyzz data: ```{r, message=FALSE} -discrete_types <- c( - "binary", "multistate-intersection", "multistate-union" -) -bp <- importBugphyzz(v = 0.5) |> - map(~ mutate(.x, NCBI_ID = as.character(NCBI_ID))) |> - map(~ { - attr_type <- unique(.x$Attribute_type) - if (attr_type %in% discrete_types) { - df <- .x |> - filter( - !(Validation < 0.7 & Evidence == "asr") - ) - } else if (attr_type == "numeric") { - df <- .x |> - filter( - !(Validation < 0.5 & Evidence == "asr") - ) - } - df - }) +## importBugphyzz() keeps imputed (ASR) annotations with MCC >= 0.7 +## (discrete) or R2 >= 0.5 (numeric) by default, the manuscript thresholds; +## since bugphyzzExports 2026-11 the exported files are already filtered. +bp <- importBugphyzz(version = params$bugphyzz_version) |> + map(~ mutate(.x, NCBI_ID = as.character(NCBI_ID))) aer <- bp$aerophilicity dim(aer) @@ -366,7 +353,7 @@ tibble( ### dbBact - genus - BugSigDB 727 ```{r} -bsdb <- importBugSigDB() |> +bsdb <- importBugSigDB(version = "10.5281/zenodo.10627578") |> # v1.2.1, as in the other articles as_tibble() tob_sigs <- bsdb |> filter(`BSDB ID` %in% c("bsdb:727/1/1", "bsdb:727/1/2")) |> diff --git a/vignettes/articles/subgingival_plaque.Rmd b/vignettes/articles/subgingival_plaque.Rmd index 967b29a..e60b6cb 100644 --- a/vignettes/articles/subgingival_plaque.Rmd +++ b/vignettes/articles/subgingival_plaque.Rmd @@ -3,6 +3,8 @@ title: "BSEA gingival plaque - MicrobiomeBenchmarkData - 35S" output: html_document: toc: true +params: + bugphyzz_version: "devel" # Zenodo DOI of the release the paper cites, or "devel" --- ```{r include = FALSE} @@ -51,26 +53,11 @@ table(tse_subset$body_subsite) Signatures at the genus level: ```{r bp signatures genus, message = FALSE} -discrete_types <- c( - "binary", "multistate-intersection", "multistate-union" -) -bp <- importBugphyzz(v = 0.5) |> - map(~ mutate(.x, NCBI_ID = as.character(NCBI_ID))) |> - map(~ { - attr_type <- unique(.x$Attribute_type) - if (attr_type %in% discrete_types) { - df <- .x |> - filter( - !(Validation < 0.7 & Evidence == "asr") - ) - } else if (attr_type == "numeric") { - df <- .x |> - filter( - !(Validation < 0.5 & Evidence == "asr") - ) - } - df - }) +## importBugphyzz() keeps imputed (ASR) annotations with MCC >= 0.7 +## (discrete) or R2 >= 0.5 (numeric) by default, the manuscript thresholds; +## since bugphyzzExports 2026-11 the exported files are already filtered. +bp <- importBugphyzz(version = params$bugphyzz_version) |> + map(~ mutate(.x, NCBI_ID = as.character(NCBI_ID))) aer <- bp$aerophilicity gn_sigs <- makeSignatures( dat = aer, taxIdType = 'Taxon_name', taxLevel = 'genus'