diff --git a/.Rbuildignore b/.Rbuildignore index 8c24f36..9b11c66 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -9,3 +9,5 @@ ^\.DS_Store$ ^data_raw$ ^doc$ +^dev$ +^docs$ diff --git a/.gitignore b/.gitignore index b5d06b4..72b4d41 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,5 @@ inst/doc /inst/extdata/interpro /data/tmp +dev/out/ +docs/*.html diff --git a/DESCRIPTION b/DESCRIPTION index ee2f1c5..c392741 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -42,6 +42,7 @@ Imports: glue, grid, gridExtra, + httr2, jsonlite, knitr, purrr, @@ -56,4 +57,4 @@ Suggests: rmarkdown, writexl, testthat (>= 3.0.0) -Config/roxygen2/version: 8.0.0 +Config/roxygen2/version: 8.1.0 diff --git a/R/bvbrc_api.R b/R/bvbrc_api.R new file mode 100644 index 0000000..2ca4675 --- /dev/null +++ b/R/bvbrc_api.R @@ -0,0 +1,223 @@ +# BV-BRC Data API path (additive; opt-in via retrieveMetadata(metadata_method = "api")) +# ----------------------------------------------------------------------------- +# Replaces ONLY the Docker/p3-* download of AMR + genome metadata. It returns the +# same two tibbles the Docker path produces -- columns prefixed `genome_drug.*` +# and `genome.*` -- so all downstream QC/join/DuckDB logic is reused unchanged. +# +# Why: the p3-* CLI path is stochastic (issue #30) -- it merges stderr into the +# data stream and never retries, so a transient BV-BRC 503 corrupts a batch. The +# API path uses typed JSON + retry-with-backoff, so a 503 is retried, never +# parsed as data. +# +# @keywords internal + +.BVBRC_API_BASE <- "https://www.bv-brc.org/api" +.BVBRC_PAGE_MAX <- 25000L + +# One HTTP request with retry on transient failures. +.bvbrcApiReq <- function(collection, rql, accept = "application/json") { + url <- sprintf("%s/%s/?%s", .BVBRC_API_BASE, collection, rql) + httr2::request(url) |> + httr2::req_headers(Accept = accept) |> + httr2::req_timeout(120) |> + httr2::req_retry( + max_tries = 5L, + is_transient = function(resp) httr2::resp_status(resp) %in% + c(429, 500, 502, 503, 504) + ) +} + +.bvbrcEnc <- function(x) utils::URLencode(as.character(x), reserved = TRUE) + +# Keyset-paginate a filtered query past the 25k cap. `key` is the unique +# sort/seek field ("id" for genome_amr, "genome_id" for genome). +.bvbrcApiFetch <- function(collection, rql_filter, select, key = "id") { + out <- list() + last <- "" + i <- 0L + repeat { + seek <- if (nzchar(last)) { + sprintf("and(%s,gt(%s,%s))", rql_filter, key, .bvbrcEnc(last)) + } else { + rql_filter + } + rql <- sprintf( + "%s&select(%s,%s)&sort(+%s)&limit(%d)", + seek, key, select, key, .BVBRC_PAGE_MAX + ) + resp <- httr2::req_perform(.bvbrcApiReq(collection, rql)) + pg <- jsonlite::fromJSON(httr2::resp_body_string(resp), + simplifyDataFrame = TRUE + ) + raw_n <- if (is.data.frame(pg)) nrow(pg) else 0L + if (raw_n == 0L) break + if (nzchar(last)) pg <- pg[pg[[key]] != last, , drop = FALSE] + if (nrow(pg) > 0L) { + i <- i + 1L + out[[i]] <- pg + last <- pg[[key]][nrow(pg)] + } + if (raw_n < .BVBRC_PAGE_MAX) break + } + if (i == 0L) { + return(tibble::tibble()) + } + tibble::as_tibble(data.table::rbindlist(out, fill = TRUE, use.names = TRUE)) +} + +# Split a vector into chunks of size n (keeps in(...) URLs within length limits). +.bvbrcChunk <- function(x, n) { + if (length(x) == 0L) { + return(list()) + } + split(x, ceiling(seq_along(x) / n)) +} + +# Ensure every expected field is present (fill missing with NA), coerce to +# character, order columns, and apply the `prefix.` naming convention that the +# Docker/p3 parser produces. +.bvbrcPrefixFill <- function(df, expected, prefix) { + df <- as.data.frame(df, stringsAsFactors = FALSE) + n <- nrow(df) # rep() keeps length right when the query returned 0 rows + for (f in expected) { + if (!f %in% names(df)) df[[f]] <- rep(NA_character_, n) + } + df <- df[, expected, drop = FALSE] + # coerce to character and use "" for missing, matching the Docker/TSV parser + # (.parse_bvbrc_tsv yields "" for blank fields, not NA) so the two paths agree. + df[] <- lapply(df, function(x) { + x <- as.character(x) + x[is.na(x)] <- "" + x + }) + names(df) <- paste0(prefix, ".", expected) + tibble::as_tibble(df) +} + +# --- AMR phenotype (genome_amr) -> genome_drug.* ------------------------------ +.extractAMRtableApi <- function(genome_ids, abx = "All", + chunk_size = 500L, verbose = TRUE) { + # Full genome_amr field set; computational_method and measurement_unit are + # populated on live BV-BRC (confirmed against the real API), so both are + # fetched -- only `source` is consistently absent and gets filled "". + expected <- c( + "genome_id", "antibiotic", "computational_method", "evidence", + "genome_name", "id", "laboratory_typing_method", + "laboratory_typing_platform", "measurement", "measurement_sign", + "measurement_unit", "measurement_value", "pmid", "resistant_phenotype", + "source", "taxon_id", "testing_standard" + ) + have <- setdiff(expected, "source") + # keyset key for genome_amr is "id"; keep genome_id as a data column. + sel <- paste(setdiff(have, "id"), collapse = ",") + + chunks <- .bvbrcChunk(genome_ids, chunk_size) + if (isTRUE(verbose)) { + message(" [api] AMR: ", length(genome_ids), " genomes in ", + length(chunks), " chunk(s)") + } + parts <- furrr::future_map( + chunks, + function(ids) { + ab <- if (identical(abx, "All")) { + "" + } else { + sprintf(",in(antibiotic,(%s))", + paste(vapply(abx, .bvbrcEnc, ""), collapse = ",")) + } + filt <- sprintf("and(in(genome_id,(%s))%s)", + paste(vapply(ids, .bvbrcEnc, ""), collapse = ","), ab) + .bvbrcApiFetch("genome_amr", filt, sel, key = "id") + }, + .options = furrr::furrr_options(seed = TRUE) + ) + df <- data.table::rbindlist(parts, fill = TRUE, use.names = TRUE) + .bvbrcPrefixFill(df, expected, "genome_drug") +} + +# --- genome-ID resolution (genome) -------------------------------------------- +# API replacement for .retrieveQueryIDs(): resolve species names and/or taxon IDs +# to Good-quality WGS/Complete genome IDs, and write the `bac_data` table (used by +# retrieveMetadata()'s summary). Uses the Data API instead of the Docker-built +# cache, so retrieveMetadata(metadata_method = "api") needs no Docker. +.resolveGenomeIDsApi <- function(base_dir = ".", user_bacs, + overwrite = FALSE, verbose = TRUE) { + sel <- "genome_name,taxon_id,species,strain" + parts <- furrr::future_map( + user_bacs, + function(ub) { + ub <- trimws(as.character(ub)) + key_filter <- if (grepl("^[0-9]+$", ub)) { + sprintf("eq(taxon_lineage_ids,%s)", ub) # taxon ID (any rank) + } else { + sprintf("eq(species,%s)", .bvbrcEnc(ub)) # species name + } + filt <- sprintf( + "and(%s,eq(genome_quality,Good),in(genome_status,(WGS,Complete)))", + key_filter + ) + if (isTRUE(verbose)) message(" [api] resolving genome IDs for '", ub, "'") + res <- .bvbrcApiFetch("genome", filt, sel, key = "genome_id") + if (nrow(res) == 0L) { + warning( + "BV-BRC API resolved 0 genomes for user_bacs entry '", ub, "'. ", + "The API path matches species names exactly and taxon IDs against ", + "the full lineage -- if this input relied on substring matching or ", + "exact-rank taxon matching under the CLI path, results will differ.", + call. = FALSE + ) + } + res + }, + .options = furrr::furrr_options(seed = TRUE) + ) + df <- as.data.frame( + data.table::rbindlist(parts, fill = TRUE, use.names = TRUE), + stringsAsFactors = FALSE + ) + if (nrow(df) == 0L) { + return(character(0)) + } + df <- df[grepl("^[0-9]+[.][0-9]+$", df$genome_id), , drop = FALSE] + df <- df[!duplicated(df$genome_id), , drop = FALSE] + + # write bac_data (genome.* columns), mirroring .retrieveQueryIDs() + paths <- .buildDBpath(base_dir = base_dir, user_bacs = user_bacs, overwrite = overwrite) + con <- DBI::dbConnect(duckdb::duckdb(), dbdir = paths$db_path) + on.exit(try(DBI::dbDisconnect(con, shutdown = TRUE), silent = TRUE), add = TRUE) + bac <- .bvbrcPrefixFill( + df, + c("genome_id", "genome_name", "taxon_id", "species", "strain"), + "genome" + ) + DBI::dbWriteTable(con, "bac_data", as.data.frame(bac), overwrite = TRUE) + if (isTRUE(verbose)) { + message(" [api] resolved ", nrow(df), " genome IDs; wrote bac_data") + } + unique(df$genome_id) +} + +# --- genome metadata (genome) -> genome.* ------------------------------------- +.extractGenomeDataApi <- function(genome_ids, fields, + chunk_size = 500L, verbose = TRUE) { + expected <- strsplit(fields, ",", fixed = TRUE)[[1]] + expected <- unique(c("genome_id", expected)) + sel <- paste(setdiff(expected, "genome_id"), collapse = ",") + + chunks <- .bvbrcChunk(genome_ids, chunk_size) + if (isTRUE(verbose)) { + message(" [api] genome metadata: ", length(genome_ids), " genomes in ", + length(chunks), " chunk(s)") + } + parts <- furrr::future_map( + chunks, + function(ids) { + filt <- sprintf("in(genome_id,(%s))", + paste(vapply(ids, .bvbrcEnc, ""), collapse = ",")) + .bvbrcApiFetch("genome", filt, sel, key = "genome_id") + }, + .options = furrr::furrr_options(seed = TRUE) + ) + df <- data.table::rbindlist(parts, fill = TRUE, use.names = TRUE) + .bvbrcPrefixFill(df, expected, "genome") +} diff --git a/R/data_curation.R b/R/data_curation.R index d5a5833..2a7442e 100644 --- a/R/data_curation.R +++ b/R/data_curation.R @@ -1166,7 +1166,9 @@ #' - genome_data #' - metadata (join on genome IDs returned by BV-BRC) #' -#' @param user_bacs Character vector. Mixed taxon IDs and/or species strings (used for naming). +#' @param user_bacs Character vector. Mixed taxon IDs and/or species strings. +#' Also used for naming the per-selection DuckDB. How entries are resolved to +#' genome IDs depends on `metadata_method` -- see Details. #' @param genome_id_file Character or NULL. Optional path to a file listing genome #' IDs (one per line). If provided, genome IDs are read from this file instead of #' being resolved from `user_bacs`; blank lines and surrounding whitespace are @@ -1177,6 +1179,8 @@ #' @param abx Character or vector. Antibiotic filter. "All" for all antibiotics, else names. #' @param overwrite Logical. If FALSE and DuckDB exists already, abort. Default FALSE. #' @param image Character. Docker image. Default "danylmb/bvbrc:5.3". +#' @param metadata_method Character. Download backend: `"api"` (default) or +#' `"cli"` (Dockerized `BV-BRC p3-* CLI`). #' @param max_checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%). #' @param min_checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%). #' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content. @@ -1185,6 +1189,28 @@ #' @param debug Logical. If TRUE, retain `metadata_full` and QC columns for inspection. #' @param verbose Logical. If TRUE, print progress messages. #' +#' @details +#' `metadata_method` selects how `user_bacs` entries are resolved to genome IDs, +#' and the two backends do not use identical matching rules. For the same input +#' they can return different genome sets: +#' +#' \itemize{ +#' \item \strong{Species strings.} `"cli"` does a case-insensitive substring +#' match against `genome.species` (so `"Escherichia"` matches +#' `"Escherichia coli"`). `"api"` does an exact match on the `species` field; +#' a string that does not match any species exactly resolves to zero genomes +#' and emits a warning. +#' \item \strong{Numeric taxon IDs.} `"cli"` matches only genomes whose own +#' `taxon_id` equals the input exactly. `"api"` matches the ID anywhere in +#' the genome's taxonomic lineage (`taxon_lineage_ids`), so a genus- or +#' family-rank ID pulls every genome beneath it, and a species-rank ID also +#' catches strain-level genomes that `"cli"` would miss. +#' } +#' +#' The `"api"` rules are generally the more complete of the two. If you need the +#' two backends to agree, pass an exact species name and a species- or +#' strain-rank taxon ID, or supply `genome_id_file` directly. +#' #' @return A list with: #' - duckdbConnection: live DBI connection to the created DuckDB #' - table_name: "metadata" @@ -1196,6 +1222,7 @@ retrieveMetadata <- function(user_bacs, abx = "All", overwrite = FALSE, image = "danylmb/bvbrc:5.3", + metadata_method = c("api", "cli"), max_checkm_contam = 5, min_checkm_complete = 95, gc_deviations = NULL, @@ -1206,6 +1233,7 @@ retrieveMetadata <- function(user_bacs, load_tables = FALSE, verbose = TRUE) { base_dir <- normalizePath(base_dir, mustWork = FALSE) + metadata_method <- match.arg(metadata_method) if (!is.null(genome_id_file)) { if (!file.exists(genome_id_file)) { @@ -1217,6 +1245,14 @@ retrieveMetadata <- function(user_bacs, genome_ids <- readLines(genome_id_file, warn = FALSE) genome_ids <- trimws(genome_ids) genome_ids <- genome_ids[genome_ids != ""] + } else if (identical(metadata_method, "api")) { + if (isTRUE(verbose)) message("Resolving genome IDs via BV-BRC API.") + genome_ids <- .resolveGenomeIDsApi( + base_dir = base_dir, + user_bacs = user_bacs, + overwrite = overwrite, + verbose = verbose + ) } else { if (isTRUE(verbose)) message("Resolving genome IDs for user inputs.") genome_ids <- .retrieveQueryIDs( @@ -1308,58 +1344,74 @@ retrieveMetadata <- function(user_bacs, batch_size <- 500L genome_batches <- split(genome_ids, ceiling(seq_along(genome_ids) / batch_size)) + # Set the future plan once, up front, so both the API path (chunk-level + # furrr::future_map in R/bvbrc_api.R) and the CLI path below run in parallel. n_cores <- max(1L, parallel::detectCores(logical = TRUE) - 1L) - old_plan <- future::plan() on.exit(future::plan(old_plan), add = TRUE) future::plan(future::multisession, workers = n_cores) - if (isTRUE(verbose)) message("Retrieving AMR phenotype data in batches.") - batch_drug_data <- furrr::future_map( - genome_batches, - function(batch) { - raw <- .extractAMRtable( - base_dir = base_dir, - batch_genome_IDs = batch, - abx_filter = abx_filter, - drug_fields = drug_fields, - image = image, - verbose = FALSE - ) - .parse_bvbrc_tsv(raw) - }, - .options = furrr::furrr_options(seed = TRUE) - ) + if (identical(metadata_method, "api")) { + # BV-BRC Data API path (Docker-free, resilient; see R/bvbrc_api.R, issue #30) + if (isTRUE(verbose)) message("Retrieving AMR phenotype data via BV-BRC API.") + combined_drug_data_tbl <- .extractAMRtableApi( + genome_ids = genome_ids, abx = abx, verbose = verbose + ) - combined_drug_data_tbl <- dplyr::bind_rows(batch_drug_data) |> - dplyr::mutate(dplyr::across(dplyr::everything(), ~ iconv(.x, from = "", to = "UTF-8", sub = ""))) + if (isTRUE(verbose)) message("Retrieving genome metadata via BV-BRC API.") + gfields <- if (identical(filter_type, "AMR")) amr_fields else microtrait_fields + combined_genome_data_tbl <- .extractGenomeDataApi( + genome_ids = genome_ids, fields = gfields, verbose = verbose + ) + } else { + if (isTRUE(verbose)) message("Retrieving AMR phenotype data in batches.") + batch_drug_data <- furrr::future_map( + genome_batches, + function(batch) { + raw <- .extractAMRtable( + base_dir = base_dir, + batch_genome_IDs = batch, + abx_filter = abx_filter, + drug_fields = drug_fields, + image = image, + verbose = FALSE + ) + .parse_bvbrc_tsv(raw) + }, + .options = furrr::furrr_options(seed = TRUE) + ) + combined_drug_data_tbl <- dplyr::bind_rows(batch_drug_data) + + if (isTRUE(verbose)) message("Retrieving genome metadata in batches.") + batch_genome_data <- furrr::future_map( + genome_batches, + function(batch) { + raw <- .extractGenomeData( + base_dir = base_dir, + batch_genome_IDs = batch, + filter_type = filter_type, + amr_fields = amr_fields, + microtrait_fields = microtrait_fields, + image = image, + verbose = FALSE + ) + .parse_bvbrc_tsv(raw) + }, + .options = furrr::furrr_options(seed = TRUE) + ) + combined_genome_data_tbl <- dplyr::bind_rows(batch_genome_data) + } + # Normalize to UTF-8 for both methods (parity with the Docker parser). + combined_drug_data_tbl <- combined_drug_data_tbl |> + dplyr::mutate(dplyr::across(dplyr::everything(), ~ iconv(.x, from = "", to = "UTF-8", sub = ""))) if (nrow(combined_drug_data_tbl) == 0L) { message("No drug data returned.") return(NULL) } - if (isTRUE(verbose)) message("Retrieving genome metadata in batches.") - batch_genome_data <- furrr::future_map( - genome_batches, - function(batch) { - raw <- .extractGenomeData( - base_dir = base_dir, - batch_genome_IDs = batch, - filter_type = filter_type, - amr_fields = amr_fields, - microtrait_fields = microtrait_fields, - image = image, - verbose = FALSE - ) - .parse_bvbrc_tsv(raw) - }, - .options = furrr::furrr_options(seed = TRUE) - ) - - combined_genome_data_tbl <- dplyr::bind_rows(batch_genome_data) |> + combined_genome_data_tbl <- combined_genome_data_tbl |> dplyr::mutate(dplyr::across(dplyr::everything(), ~ iconv(.x, from = "", to = "UTF-8", sub = ""))) - if (nrow(combined_genome_data_tbl) == 0L) { message("No genome data returned.") return(NULL) @@ -2141,8 +2193,10 @@ genomeList <- function(base_dir = ".", #' metadata step is restricted to these genome IDs instead of resolving them from #' `user_bacs`. Default NULL. #' @param base_dir Character. Project root directory. Default `"."`. -#' @param method Character. Download method passed to `retrieveGenomes()`. +#' @param method Character. Genome download method passed to `retrieveGenomes()`. #' `"ftp"` (default) or `"cli"`. +#' @param metadata_method Character. Metadata download method passed to `retrieveMetadata()`. +#' `"api"` (default) or `"cli"`. #' @param overwrite Logical. Passed to metadata filtering and DuckDB creation. #' Default FALSE. #' @param evidence_mode Character. Sets what types of AMR evidence is acceptable. @@ -2167,6 +2221,7 @@ prepareGenomes <- function(user_bacs, genome_id_file = NULL, base_dir = ".", method = c("ftp", "cli"), + metadata_method = c("api", "cli"), overwrite = FALSE, num_workers = 8L, evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"), @@ -2180,6 +2235,7 @@ prepareGenomes <- function(user_bacs, debug = FALSE, verbose = TRUE) { method <- match.arg(method) + metadata_method <- match.arg(metadata_method) evidence_mode <- match.arg(evidence_mode) base_dir <- normalizePath(base_dir, mustWork = FALSE) @@ -2214,6 +2270,7 @@ prepareGenomes <- function(user_bacs, message = "Started genome curation run.", details = list( method = method, + metadata_method = metadata_method, evidence_mode = evidence_mode, overwrite = overwrite ) @@ -2231,24 +2288,40 @@ prepareGenomes <- function(user_bacs, add = TRUE ) - manifest <- .manifest_stage( - manifest, - name = "prepare_bvbrc_cache", - status = "success", - parameters = list( - max_age_days = 30L - ), - outputs = file.path(base_dir, "data", "bvbrc", "bvbrcData.duckdb"), - tool = list( - name = "BV-BRC", - interface = "p3-all-genomes" + # The Docker/p3-all-genomes cache only backs the "cli" method; the "api" + # method queries BV-BRC directly and has no cache-age concept. Record the + # stage either way so the manifest never has a silent gap here. + if (identical(metadata_method, "cli")) { + manifest <- .manifest_stage( + manifest, + name = "prepare_bvbrc_cache", + status = "success", + parameters = list( + max_age_days = 30L + ), + outputs = file.path(base_dir, "data", "bvbrc", "bvbrcData.duckdb"), + tool = list( + name = "BV-BRC", + interface = "p3-all-genomes" + ) ) - ) - .ensure_bvbrc_cache( - base_dir = base_dir, - verbose = verbose - ) + .ensure_bvbrc_cache( + base_dir = base_dir, + verbose = verbose + ) + } else { + manifest <- .manifest_stage( + manifest, + name = "prepare_bvbrc_cache", + status = "skipped", + message = "metadata_method = \"api\" queries BV-BRC directly; no Docker cache to prepare.", + tool = list( + name = "BV-BRC", + interface = "Data API" + ) + ) + } if (isTRUE(verbose)) { message("Step 0: Building AMR metadata (retrieveMetadata)") @@ -2261,6 +2334,7 @@ prepareGenomes <- function(user_bacs, parameters = list( filter_type = "AMR", abx = "All", + metadata_method = metadata_method, max_checkm_contam = max_checkm_contam, min_checkm_complete = min_checkm_complete, gc_deviations = gc_deviations, @@ -2283,6 +2357,7 @@ prepareGenomes <- function(user_bacs, base_dir = base_dir, abx = "All", overwrite = overwrite, + metadata_method = metadata_method, max_checkm_contam = max_checkm_contam, min_checkm_complete = min_checkm_complete, gc_deviations = gc_deviations, @@ -2308,10 +2383,11 @@ prepareGenomes <- function(user_bacs, overwrite = overwrite ), outputs = normalizePath(paths$db_path, mustWork = FALSE), - tool = list( - name = "BV-BRC", - docker_image = "danylmb/bvbrc:5.3" - ) + tool = if (identical(metadata_method, "api")) { + list(name = "BV-BRC", interface = "Data API") + } else { + list(name = "BV-BRC", interface = "BV-BRC CLI", docker_image = "danylmb/bvbrc:5.3") + } ) if (isTRUE(verbose)) message("Step 1: Filtering genomes for download by evidence: ", evidence_mode) diff --git a/data/bvbrc/bvbrcData.duckdb b/data/bvbrc/bvbrcData.duckdb deleted file mode 100644 index 3f96a3a..0000000 Binary files a/data/bvbrc/bvbrcData.duckdb and /dev/null differ diff --git a/dev/parity_check.R b/dev/parity_check.R new file mode 100644 index 0000000..c9ed419 --- /dev/null +++ b/dev/parity_check.R @@ -0,0 +1,74 @@ +# Parity check: cli vs api on the SAME genomes. +# Prereqs: +# 1. Docker running + `docker pull danylmb/bvbrc:5.3` +# 2. Install the package first (the Docker path uses furrr workers that need +# it INSTALLED, not just load_all'd): R CMD INSTALL . +# Then run from the repo root: Rscript dev/parity_check.R +suppressPackageStartupMessages(library(amRdata)) + +species <- "Morganella morganii" + +# Fix the genome set so BOTH methods pull the same genomes (isolates the +# download comparison from any ID-resolution differences). +ids <- amRdata:::.resolveGenomeIDsApi( + base_dir = tempfile(), user_bacs = species, overwrite = TRUE, verbose = FALSE +) +ids <- utils::head(ids, 30) +gf <- tempfile(fileext = ".txt") +writeLines(ids, gf) + +run <- function(metadata_method) { + td <- file.path(tempdir(), paste0("parity_", metadata_method)) + unlink(td, recursive = TRUE) + dir.create(td, recursive = TRUE) + invisible(retrieveMetadata( + user_bacs = species, genome_id_file = gf, metadata_method = metadata_method, + base_dir = td, overwrite = TRUE, verbose = FALSE + )) + db <- list.files(td, pattern = "[.]duckdb$", recursive = TRUE, full.names = TRUE)[1] + con <- DBI::dbConnect(duckdb::duckdb(), dbdir = db, read_only = TRUE) + on.exit(DBI::dbDisconnect(con, shutdown = TRUE)) + list( + genome = DBI::dbReadTable(con, "genome_data"), + amr = DBI::dbReadTable(con, "amr_phenotype") + ) +} + +d <- run("cli") +a <- run("api") + +cat("\n================ PARITY: cli vs api ================\n") +cat(sprintf("genome_data rows cli=%d api=%d\n", nrow(d$genome), nrow(a$genome))) +cat(sprintf("amr_phenotype rows cli=%d api=%d\n", nrow(d$amr), nrow(a$amr))) +cat("same genome set:", + setequal(d$genome[["genome.genome_id"]], a$genome[["genome.genome_id"]]), "\n") + +# AMR: compare the (genome, antibiotic, phenotype) tuples as sets +tup <- function(x) { + sort(paste( + x[["genome_drug.genome_id"]], + x[["genome_drug.antibiotic"]], + x[["genome_drug.resistant_phenotype"]] + )) +} +cat("AMR (genome,antibiotic,phenotype) identical:", identical(tup(d$amr), tup(a$amr)), "\n") + +# Genome metadata: compare core fields for shared IDs +core <- c("genome.genome_id", "genome.species", "genome.genome_quality", + "genome.genome_status", "genome.checkm_completeness", + "genome.checkm_contamination", "genome.cds", "genome.genome_length", + "genome.gc_content") +ord <- function(x) { + cols <- intersect(core, names(x)) + x <- x[order(x[["genome.genome_id"]]), cols, drop = FALSE] + x[] <- lapply(x, as.character) + as.data.frame(x, stringsAsFactors = FALSE) +} +cat("core genome fields identical:", identical(ord(d$genome), ord(a$genome)), "\n") + +cat("\ncolumns only in DOCKER genome_data:", + paste(setdiff(names(d$genome), names(a$genome)), collapse = ", "), "\n") +cat("columns only in API genome_data:", + paste(setdiff(names(a$genome), names(d$genome)), collapse = ", "), "\n") +cat("\nExpected known diffs: api fills source as \"\" in AMR (consistently absent", + "from genome_amr). Focus on the 'identical' lines.\n") diff --git a/dev/parity_diff.R b/dev/parity_diff.R new file mode 100644 index 0000000..7468c70 --- /dev/null +++ b/dev/parity_diff.R @@ -0,0 +1,50 @@ +# Diagnose the AMR tuple mismatch between cli and api. +# Prereqs: same as dev/parity_check.R (Docker running + package installed). +# Run: Rscript dev/parity_diff.R +suppressPackageStartupMessages(library(amRdata)) + +species <- "Morganella morganii" +ids <- amRdata:::.resolveGenomeIDsApi( + base_dir = tempfile(), user_bacs = species, overwrite = TRUE, verbose = FALSE +) +ids <- utils::head(ids, 30) +gf <- tempfile(fileext = ".txt") +writeLines(ids, gf) + +amr <- function(metadata_method) { + td <- file.path(tempdir(), paste0("pd_", metadata_method)) + unlink(td, recursive = TRUE) + dir.create(td, recursive = TRUE) + invisible(retrieveMetadata( + user_bacs = species, genome_id_file = gf, metadata_method = metadata_method, + base_dir = td, overwrite = TRUE, verbose = FALSE + )) + db <- list.files(td, pattern = "[.]duckdb$", recursive = TRUE, full.names = TRUE)[1] + con <- DBI::dbConnect(duckdb::duckdb(), dbdir = db, read_only = TRUE) + on.exit(DBI::dbDisconnect(con, shutdown = TRUE)) + DBI::dbReadTable(con, "amr_phenotype") +} + +d <- amr("cli") +a <- amr("api") + +tup <- function(x) { + paste(x[["genome_drug.genome_id"]], x[["genome_drug.antibiotic"]], + x[["genome_drug.resistant_phenotype"]], sep = " | ") +} +dt <- tup(d) +at <- tup(a) + +cat("\n--- tuples in DOCKER but not API (", length(setdiff(dt, at)), ") ---\n", sep = "") +print(utils::head(sort(setdiff(dt, at)), 12)) +cat("\n--- tuples in API but not DOCKER (", length(setdiff(at, dt)), ") ---\n", sep = "") +print(utils::head(sort(setdiff(at, dt)), 12)) + +show <- function(x, col) paste(sort(unique(x[[col]])), collapse = " | ") +cat("\ndistinct ANTIBIOTIC — docker:\n ", show(d, "genome_drug.antibiotic"), "\n") +cat("distinct ANTIBIOTIC — api:\n ", show(a, "genome_drug.antibiotic"), "\n") +cat("\ndistinct PHENOTYPE — docker:", show(d, "genome_drug.resistant_phenotype"), "\n") +cat("distinct PHENOTYPE — api: ", show(a, "genome_drug.resistant_phenotype"), "\n") + +# are duplicates the cause? (same count but different multiplicity) +cat("\ndup tuples docker:", sum(duplicated(dt)), " | dup tuples api:", sum(duplicated(at)), "\n") diff --git a/docs/bvbrc-api-feasibility.md b/docs/bvbrc-api-feasibility.md new file mode 100644 index 0000000..fa994a0 --- /dev/null +++ b/docs/bvbrc-api-feasibility.md @@ -0,0 +1,301 @@ +# BV-BRC Data API — Feasibility & Limitations for amRdata + +**Status:** evaluation notes (not yet implemented) +**Date:** 2026-07-24 +**Author:** Emily Boyer +**Scope:** Can we replace the current Docker/`p3-*` CLI download path with direct +calls to the BV-BRC Data API, to fix the download bottleneck and let users pull +genome metadata + AMR phenotype data together, choosing their own columns? + +All numbers below were measured with live `curl` calls against +`https://www.bv-brc.org/api/` on the date above. No authentication is required +for public data. + +--- + +## 1. Bottom line + +**The API is viable and removes the bottleneck. We should call it directly from R +(`httr2`), and NOT adopt the `bvbrc` Python package.** + +- The API is a plain REST/Solr service. Everything we need — column selection, + counting without downloading, faceted subsampling — is a URL query. +- The Python package (`bvbrc` v0.2.1, Beta, 1 maintainer, 9 commits) is a thin + wrapper over this same API. It offers no capability R can't reach with `httr2`, + and it does **not** solve the one hard part (the 25,000-row cap — see §4). It + would add a `reticulate` + Python 3.9–3.11 dependency for no net gain. Use its + client design as a reference, not as a dependency. + +--- + +## 2. What our use case needs → feasibility + +Derived from the AMR team meeting notes. + +| Requirement | Feasible | Mechanism | +|---|---|---| +| User specifies **which columns** to return | ✅ | `select(field1,field2,...)` | +| Count matches **without downloading** ("-n" option) | ✅ | `Content-Range` response header — instant | +| Pick **species/pathogen from a list** (not typed) | ✅ | Solr faceting on `species` (genome collection) | +| Subsample "**500 genomes with most coverage** for drugs of interest" | ✅ | Facet `genome_id` filtered by `in(antibiotic,(...))`, ranked by count | +| Pull **QC-passing snapshot** (Good genomes), relevant columns → Zenodo | ✅ | `eq(genome_quality,Good)` + `select(...)` + keyset paging (§4) | +| Get metadata + AMR phenotype **together** | ⚠️ partial | Two collections (`genome`, `genome_amr`); no server-side join. But each call is sub-second, so it's 2 fast pulls + a local join instead of 2 slow Docker pipelines | + +The "together" goal isn't literally a single call, but it stops being a +bottleneck: today's pain is Docker container startup + serial TSV parsing, not +the join. + +--- + +## 3. Verified capabilities & scale + +| Metric | Value | +|---|---| +| Total `genome_amr` rows (whole DB) | 17,266,649 | +| Total `genome` records | 16,885,561 | +| "Good" quality genomes | 8,117,447 | +| Throughput | 25,000 rows / ~1.9 s / ~2.2 MB (4 columns) | +| Column selection (`select`) | works | +| Count-only (`Content-Range`) | works, instant | +| Faceting (species list; genome ranking by drug coverage) | works | +| Auth for public data | none required | + +**Query interface:** RQL over HTTP GET/POST, e.g. +`.../genome_amr/?and(eq(genome_name,Staphylococcus aureus),in(antibiotic,(ciprofloxacin,gentamicin)))&select(genome_id,antibiotic,resistant_phenotype)&sort(+id)&limit(25000)` +(URL-encode field names and values individually). POST the query body for large +`in(...)` lists to avoid URL-length limits. + +--- + +## 4. The key limitation: the 25,000-row cap — and how to get around it + +**Two hard ceilings, both = 25,000:** + +1. **Per-request cap.** `limit(30000)` returns only 25,000 rows. +2. **Offset ceiling.** `limit(count, start)` with `start >= 25000` returns + **HTTP 400**. So you cannot walk a large result set by increasing the offset. + +`cursorMark` (Solr's deep-scroll) is **not reachable through the RQL layer** — +passing `&cursorMark=*` errors with `undefined field: "cursorMark"`. + +### Answer to "we don't know how much we'll need to download at a given time" + +Use a **two-step pattern** that scales to any size without knowing it in advance: + +**Step 1 — Count first (free, instant).** Read the `Content-Range` header on a +`limit(1)` request. `items 0-1/590284` means 590,284 matching rows. Now you know +the exact size before downloading anything. This directly resolves the "we don't +know how much" problem — you always find out up front, for zero cost. + +**Step 2 — If the count > 25,000, use keyset (seek) pagination.** Instead of an +offset, sort by the unique `id` and, each page, ask for rows *after* the last id +you saw. `start` is always 0, so the offset ceiling never applies. Loop until a +page returns fewer than 25,000 rows. This handles arbitrary/unknown totals. + +``` +# Pseudocode — retrieve ALL rows for a query, any size +total <- GET .../genome_amr/?&limit(1) # read Content-Range -> N +rows <- [] +last_id <- "" # empty sorts before all ids +repeat: + page <- GET .../genome_amr/?and(, gt(id, last_id)) + &select()&sort(+id)&limit(25000) + if last_id != "": # gt(id,...) is INCLUSIVE on this API + drop page[0] if page[0].id == last_id # dedupe boundary row + rows <- rows + page + last_id <- page[last].id + until length(page) < 25000 +``` + +**Gotcha (verified):** `gt(id, X)` behaves as `>=` here — the boundary row +reappears as the first row of the next page. Dedupe on `id` when stitching pages, +or drop the first row of every page after the first. + +**Cost example:** the full 17.3M-row `genome_amr` collection is ~692 pages of +25k. At ~1.9 s/page that's ~22 min sequential for the entire AMR table — or a few +minutes with parallel id-partitioning (see below). And we never need the whole +thing; a QC-filtered, column-selected subset is far smaller. Contrast with the +current per-500-genome Docker-container approach. + +### Alternative: partition the query +If a single query is huge, split it by a natural key (species, taxon, or year) +so each sub-query is < 25k, and page each partition. Keyset paging is simpler and +preferred; partitioning is a fallback when one key range dominates. + +### Looping vs. parallelization (both measured) + +**Looping is the keyset walk itself, and it is inherently sequential *within* a +result set:** page N needs the last `id` from page N-1, so a single keyset walk +cannot be parallelized. + +**Parallelization works by partitioning the `id` space** into disjoint ranges and +giving each worker its own range. It does **not** work by offset — the 25k offset +ceiling (§4) rules out `limit(25000, 25000)`-style parallel offsets entirely. + +Because `id` is a UUID (uniformly distributed hex), splitting on the first hex +character gives 16 disjoint, naturally **balanced** buckets — no knowledge of the +data distribution required. Verified bucket sizes: + +- rows with `id` in `[0,1)`: 1,078,515 +- rows with `id` in `[1,2)`: 1,078,740 +- × 16 buckets ≈ 17.3M total ✓ + +Each bucket is bounded with `and(ge(id,),lt(id,))` and is still >25k, so a +worker **keyset-paginates within its bucket**. Two-level structure: + +``` +partition id space -> [0,1), [1,2), ... [f,g) # 16 disjoint buckets + |-- worker pool (BiocParallel / future) --| + each worker: keyset-paginate its bucket(s) with gt(id,last) + sort(+id) +``` + +**Measured speedup:** 4 workers pulling 25k-row pages from 4 buckets took **5.1 s +vs 18.5 s sequential (~3.6×), with no throttling observed.** This maps directly +onto the parallel backend the package already uses (`BiocParallel` / `future`). + +**Be a good citizen.** Rate limits are undocumented, so keep concurrency modest +(≈4–8 workers, not 16+), reuse connections, and add retry-with-backoff on +transient failures (HTTP 5xx) — the same robustness gap that bites the current +Docker path. Finer partitioning (2 hex chars = 256 buckets) is available if you +want a work queue that load-balances across a fixed worker pool. + +--- + +## 5. Other limitations / open items + +- **`genome` and `genome_amr` are separate collections.** No server-side join; + join locally on `genome_id`. Fine, since each pull is fast. +- **Distinct-genome counts** (e.g. "how many genomes actually have AMR data" — the + true Zenodo snapshot size) need Solr's native `json.facet`/`unique()`, which the + RQL facet layer did not return cleanly. Obtainable via the native Solr interface + (same host, `Accept: application/solr+json`); falls out of the same work as deep + scrolling. Not yet measured. `unique()` is an HLL estimate, not exact. +- **Genome sequence files** (`.fna`/`.faa`/`.gff`) are a separate concern — the + Data API returns metadata/phenotype records, not assembly files. Those still + come from FTP/CLI. This evaluation covers metadata + AMR phenotype only. +- **Rate limits** are not documented; be a good citizen (reuse connections, avoid + hammering, prefer count-first over speculative pulls). + +--- + +## 6. Recommendation + +1. Build a small R module (`httr2`) with: + - a query builder (RQL: `eq`/`in`/`gt`/`select`/`sort`/`limit`), + - a `count()` helper reading `Content-Range`, + - a `fetch_all()` that keyset-paginates past 25k with boundary dedupe, with + optional id-space partitioning across a parallel worker pool for bulk pulls, + - faceting helpers for the "pick a species" and "rank genomes by drug + coverage" subsampling steps. +2. Benchmark it against the current Docker `.extractAMRtable()` / + `.extractGenomeData()` path. +3. Do **not** take on the Python package as a dependency. + +--- + +## 7. Prototype & benchmark + +The findings below came from a standalone prototype (`bvbrc_count()`, +`bvbrc_fetch_all()` with keyset + optional parallel id-partitioning, and +`bvbrc_rank_genomes_by_drug_coverage()`). That prototype has since been +superseded by the package implementation in [`R/bvbrc_api.R`](../R/bvbrc_api.R). + +**Measured (all S. aureus AMR rows):** + +| Operation | Result | +|---|---| +| `bvbrc_count()` | 590,284 rows — instant (header only) | +| `bvbrc_fetch_all()` sequential | 590,284 rows in **72 s** | +| `bvbrc_fetch_all(parallel=TRUE, workers=6)` | 590,284 rows in **27 s (~2.6×)**, identical row count | +| keyset correctness (small set, tiny page size) | rows == count, all ids unique, no boundary duplication | + +Parallel returns the exact same row count as sequential — the id-space +partitioning is complete and non-overlapping. + +--- + +## 8. Other amRdata use cases the API covers + +The API is not just for AMR phenotype — it reaches most of what the package +currently shells out to Docker (`p3-*`) for. All verified live: + +| Need (current source) | API collection | Verified | +|---|---|---| +| AMR phenotype (`genome_amr` via `p3-get-genome-drugs`) | `genome_amr` | ✅ | +| Genome metadata (via `p3-get-genome-data`) | `genome` | ✅ | +| Gene/annotation content — `.PATRIC.gff` (via `p3-dump-genomes`) | `genome_feature` | ✅ 2,583 CDS features w/ `patric_id`, `product`, coords, strand, `aa_sequence_md5` | +| **Protein FASTA — `.PATRIC.faa`** (via `p3-dump-genomes`) | `genome_feature` → `feature_sequence` | ✅ md5 → actual AA sequence returned | +| Contig DNA — `.fna` | `genome_sequence` | ✅ contigs w/ length/GC (+ `sequence` field available) | +| **Genotypic AMR / specialty genes** (new) | `sp_gene` | ✅ 574 hits for one genome; `property=Antibiotic Resistance` filterable | + +**Implications:** + +- **Protein FASTA for CD-HIT clustering can come from the API** — pull + `genome_feature` (coords + `aa_sequence_md5`), then batch-fetch sequences from + `feature_sequence` by md5. Removes Docker from the protein path. For thousands + of genomes, dedupe md5s before fetching (identical proteins share an md5). +- **`sp_gene` is a new, relevant data source** — genotypic AMR determinants + (resistance genes) to complement the phenotypic `genome_amr`. Not currently + used by amRdata; worth considering for AMR modeling features. +- **Whole-assembly bulk files** (`.fna`) are the one case where **FTP is still + the better route** — reconstructing multi-MB assemblies from per-contig API + sequence fields is far heavier than downloading the flat file. FTP + (`ftp.bvbrc.org/genomes//`) was **unreachable from the eval sandbox + (network-blocked, HTTP 000), so verify from a real machine** — but it remains + the documented, efficient path for sequence files. Keep FTP for `.fna`; use the + API for metadata, features, protein sequences, and genotypic AMR. + +**Net:** the API can replace the Docker/`p3-*` path for everything except bulk +assembly (`.fna`) downloads, while adding faceted subsampling and a genotypic-AMR +source the package doesn't currently have. + +--- + +## 9. Per-species roster & the "adjustable row limit" request + +Motivation: the `bvbrc` Python package caps `limit="max"` at 25,000, which is the +**server's** ceiling (see §4 — `limit(30000)` returns 25,000, and `start >= 25000` +is rejected). Goal: gather evidence to ask the package maintainer to stop +silently truncating at 25k. + +**Critical framing (get this right in the request):** you cannot "just raise the +number." 25,000 is enforced by BV-BRC's API, not the package. The correct ask is +**automatic pagination** (keyset — §4) so a query that matches N > 25,000 rows +returns all N instead of a silent first 25k. Expose it as e.g. `limit="all"` or a +`max_rows` parameter that paginates under the hood. A bigger single `limit` value +will be clamped to 25k server-side and change nothing. + +A one-off roster script produced the counts below (header-only, via +`Content-Range`); definitions: +`clean` = `genome_quality = Good`; `amr_rows` = `genome_amr` rows via +`eq(genome_name, )`. + +**Headline:** of 25 WHO(2024)/CDC(2019) bacterial priority species, +**14 exceed the 25k AMR-row cap** and **7 exceed 25k in genome count** (so even +the genome-metadata pull truncates). Worst case **E. coli: 7,388,629 AMR rows — +a capped single pull returns 0.3% of the data.** Other 7-figure species: +M. tuberculosis 2.26M, K. pneumoniae 2.01M, S. enterica 1.97M, S. pneumoniae +1.06M. + +**Caveat on `clean_with_amr`:** it counts genomes whose `genome` record has the +`antimicrobial_resistance` summary field populated — which is **sparsely filled** +and undercounts genomes that actually have `genome_amr` phenotype rows (e.g. +*C. jejuni*: 478,916 AMR rows but only 96 flagged genomes). For an accurate +"clean genomes with phenotype data," intersect Good `genome_id`s with the distinct +`genome_id`s present in `genome_amr` (heavier; via native Solr faceting). Treat the +`clean_with_amr` column as a floor, not a true count. + +### Recommended minimal `genome` column set ("optimize / ditch the rest") + +All verified present on the `genome` collection; use in `select(...)`: + +`genome_id`, `assembly_accession` (NCBI/GCA), `genbank_accessions`, +`genome_quality`, `genome_status`, `checkm_completeness`, `checkm_contamination`, +`cds`, `genome_length`, `gc_content`, `host_name`, `isolation_country`, +`geographic_group`, `species`, `taxon_id`. + +## References +- BV-BRC Data API: https://www.bv-brc.org/api/doc/ +- `bvbrc` Python package: https://pypi.org/project/bvbrc/ · + https://github.com/abates20/bvbrc · + https://bvbrc.readthedocs.io/en/latest/ diff --git a/man/prepareGenomes.Rd b/man/prepareGenomes.Rd index 0f7814a..5701feb 100644 --- a/man/prepareGenomes.Rd +++ b/man/prepareGenomes.Rd @@ -9,6 +9,7 @@ prepareGenomes( genome_id_file = NULL, base_dir = ".", method = c("ftp", "cli"), + metadata_method = c("api", "cli"), overwrite = FALSE, num_workers = 8L, evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"), @@ -34,9 +35,12 @@ metadata step is restricted to these genome IDs instead of resolving them from \item{base_dir}{Character. Project root directory. Default \code{"."}.} -\item{method}{Character. Download method passed to \code{retrieveGenomes()}. +\item{method}{Character. Genome download method passed to \code{retrieveGenomes()}. \code{"ftp"} (default) or \code{"cli"}.} +\item{metadata_method}{Character. Metadata download method passed to \code{retrieveMetadata()}. +\code{"api"} (default) or \code{"cli"}.} + \item{overwrite}{Logical. Passed to metadata filtering and DuckDB creation. Default FALSE.} diff --git a/man/retrieveMetadata.Rd b/man/retrieveMetadata.Rd index 90d8cdb..fb1a1ee 100644 --- a/man/retrieveMetadata.Rd +++ b/man/retrieveMetadata.Rd @@ -12,6 +12,7 @@ retrieveMetadata( abx = "All", overwrite = FALSE, image = "danylmb/bvbrc:5.3", + metadata_method = c("api", "cli"), max_checkm_contam = 5, min_checkm_complete = 95, gc_deviations = NULL, @@ -24,7 +25,9 @@ retrieveMetadata( ) } \arguments{ -\item{user_bacs}{Character vector. Mixed taxon IDs and/or species strings (used for naming).} +\item{user_bacs}{Character vector. Mixed taxon IDs and/or species strings. +Also used for naming the per-selection DuckDB. How entries are resolved to +genome IDs depends on \code{metadata_method} -- see Details.} \item{genome_id_file}{Character or NULL. Optional path to a file listing genome IDs (one per line). If provided, genome IDs are read from this file instead of @@ -42,6 +45,9 @@ Default NULL.} \item{image}{Character. Docker image. Default "danylmb/bvbrc:5.3".} +\item{metadata_method}{Character. Download backend: \code{"api"} (default) or +\code{"cli"} (Dockerized \verb{BV-BRC p3-* CLI}).} + \item{max_checkm_contam}{Numeric scalar. Maximum allowed CheckM contamination (\%).} \item{min_checkm_complete}{Numeric scalar. Minimum allowed CheckM completeness (\%).} @@ -74,3 +80,25 @@ Tables written: \item metadata (join on genome IDs returned by BV-BRC) } } +\details{ +\code{metadata_method} selects how \code{user_bacs} entries are resolved to genome IDs, +and the two backends do not use identical matching rules. For the same input +they can return different genome sets: + +\itemize{ +\item \strong{Species strings.} \code{"cli"} does a case-insensitive substring +match against \code{genome.species} (so \code{"Escherichia"} matches +\code{"Escherichia coli"}). \code{"api"} does an exact match on the \code{species} field; +a string that does not match any species exactly resolves to zero genomes +and emits a warning. +\item \strong{Numeric taxon IDs.} \code{"cli"} matches only genomes whose own +\code{taxon_id} equals the input exactly. \code{"api"} matches the ID anywhere in +the genome's taxonomic lineage (\code{taxon_lineage_ids}), so a genus- or +family-rank ID pulls every genome beneath it, and a species-rank ID also +catches strain-level genomes that \code{"cli"} would miss. +} + +The \code{"api"} rules are generally the more complete of the two. If you need the +two backends to agree, pass an exact species name and a species- or +strain-rank taxon ID, or supply \code{genome_id_file} directly. +} diff --git a/tests/testthat/test-bvbrc_api.R b/tests/testthat/test-bvbrc_api.R new file mode 100644 index 0000000..7a23817 --- /dev/null +++ b/tests/testthat/test-bvbrc_api.R @@ -0,0 +1,93 @@ +# Tests for the BV-BRC Data API download path (R/bvbrc_api.R). +# Pure helpers run offline; the live extractor tests skip when offline/CRAN. + +test_that(".bvbrcChunk splits into size-n groups", { + expect_length(.bvbrcChunk(1:10, 3), 4) + expect_length(.bvbrcChunk(character(0), 5), 0) + expect_identical(unname(.bvbrcChunk(1:3, 10)[[1]]), 1:3) +}) + +test_that(".bvbrcPrefixFill fills missing fields, prefixes, coerces to character", { + df <- data.frame( + genome_id = "1280.1", antibiotic = "ciprofloxacin", + stringsAsFactors = FALSE + ) + out <- .bvbrcPrefixFill(df, c("genome_id", "antibiotic", "source"), "genome_drug") + + expect_identical( + names(out), + c("genome_drug.genome_id", "genome_drug.antibiotic", "genome_drug.source") + ) + expect_identical(out[["genome_drug.source"]], "") # missing field -> "" (Docker convention) + expect_type(out[["genome_drug.genome_id"]], "character") +}) + +test_that(".extractAMRtableApi returns genome_drug.* columns keyed by genome_id (live)", { + skip_on_cran() + skip_if_offline("www.bv-brc.org") + + amr <- .extractAMRtableApi("1280.15865", verbose = FALSE) + expect_s3_class(amr, "data.frame") + expect_true("genome_drug.genome_id" %in% names(amr)) + expect_true(all(grepl("^genome_drug[.]", names(amr)))) + expect_gt(nrow(amr), 0) + expect_true(all(amr[["genome_drug.genome_id"]] == "1280.15865")) +}) + +test_that(".bvbrcPrefixFill handles a zero-row frame (empty query result)", { + out <- .bvbrcPrefixFill(data.frame(), c("genome_id", "antibiotic"), "genome_drug") + expect_identical(names(out), c("genome_drug.genome_id", "genome_drug.antibiotic")) + expect_equal(nrow(out), 0L) +}) + +test_that(".resolveGenomeIDsApi resolves a species to valid genome IDs (live)", { + skip_on_cran() + skip_if_offline("www.bv-brc.org") + + td <- file.path(tempdir(), paste0("res_", as.integer(runif(1, 1, 1e6)))) + dir.create(td, showWarnings = FALSE, recursive = TRUE) + on.exit(unlink(td, recursive = TRUE), add = TRUE) + + ids <- .resolveGenomeIDsApi( + base_dir = td, user_bacs = "Morganella morganii", + overwrite = TRUE, verbose = FALSE + ) + expect_type(ids, "character") + expect_gt(length(ids), 0L) + expect_true(all(grepl("^[0-9]+[.][0-9]+$", ids))) + expect_false(any(duplicated(ids))) +}) + +test_that(".resolveGenomeIDsApi warns (not errors) on a zero-match input", { + skip_on_cran() + skip_if_offline("www.bv-brc.org") + + td <- file.path(tempdir(), paste0("res0_", as.integer(runif(1, 1, 1e6)))) + dir.create(td, showWarnings = FALSE, recursive = TRUE) + on.exit(unlink(td, recursive = TRUE), add = TRUE) + + expect_warning( + ids <- .resolveGenomeIDsApi( + base_dir = td, user_bacs = "Nosuchgenusxyzabc", + overwrite = TRUE, verbose = FALSE + ), + "resolved 0 genomes" + ) + expect_identical(ids, character(0)) +}) + +test_that(".extractGenomeDataApi returns genome.* columns incl. QC fields (live)", { + skip_on_cran() + skip_if_offline("www.bv-brc.org") + + gen <- .extractGenomeDataApi( + "1280.15865", + fields = "species,genome_quality,checkm_completeness,cds", + verbose = FALSE + ) + expect_true(all( + c("genome.genome_id", "genome.species", "genome.checkm_completeness") %in% + names(gen) + )) + expect_equal(nrow(gen), 1L) +})