Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 10 additions & 2 deletions inst/rmd/bsdb_signatures.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down Expand Up @@ -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)) {
Expand All @@ -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)
```

Expand Down
108 changes: 85 additions & 23 deletions vignettes/articles/data_integration_and_stats.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down Expand Up @@ -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
}

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Without validation_tsv, the validation table comes from the installed bugphyzz via the internal bugphyzz:::.validationData(), not from the release named by bugphyzz_version. If you pass a DOI and leave validation_tsv empty, Figure 2 and the new Supplementary Table 2 can show MCC/R2 values from a different version than the annotation counts next to them. Nothing warns about this.

Suggested fixes:

  • Make validation_tsv required whenever bugphyzz_version != "devel", or stop() / warning() when the two are inconsistent.
  • Replace ::: with an exported accessor. Bioconductor guidelines discourage :::, and an internal rename would silently break the vignette. If no accessor exists, export one in bugphyzz.

length(bp)
```
```{r}
Expand Down Expand Up @@ -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) |>
Expand All @@ -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
)

Expand Down Expand Up @@ -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)) |>
Expand Down Expand Up @@ -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)
Comment on lines +629 to +637

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The new supplementary table has no test or check. Total is n() while Curated + IBD + ASR counts only specific Evidence codes, so any other or NA Evidence value makes the columns stop summing to Total without notice. The join to validation is by tolower(Attribute), so a mismatch gives silent NA in the validation column.

Suggested fixes:

  • Add stopifnot(with(supp_attr_tbl, all(Curated + IBD + ASR == Total))) after the table is built.
  • Assert that every row has a validation string, or document which attributes are expected to be NA.
  • Note that Sources / Type use unique(...)[1], which hides inconsistent values within an attribute.

)
}) |>
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)
```
Comment on lines +654 to +660

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This chunk writes TSVs to the knit working directory. That is a side effect on the source tree when built with R CMD build/BiocCheck, and the location varies with how the article is rendered (rmarkdown::render vs pkgdown). It also uses an unnamed relative path, and supplementary_table_sources.tsv is not mentioned in the PR's table list.

Fix: take an output directory from a param (e.g. params$out_dir, defaulting to tempdir()), build the path with file.path(), and use readr::write_tsv (readr is already in Imports). Also add *.tsv outputs to .gitignore/.Rbuildignore if they can land in the repo.


# Session information

```{r class.soruce = 'fold-show'}
Expand Down
35 changes: 14 additions & 21 deletions vignettes/articles/enrichment_dbbact.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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)) {
Expand Down
29 changes: 8 additions & 21 deletions vignettes/articles/nychanes.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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")) |>
Expand Down
27 changes: 7 additions & 20 deletions vignettes/articles/subgingival_plaque.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down Expand Up @@ -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'
Expand Down
Loading