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..e53ffc5 100644 --- a/R/data_curation.R +++ b/R/data_curation.R @@ -1177,6 +1177,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. @@ -1196,6 +1198,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 +1209,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 +1221,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 +1320,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 +2169,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 +2197,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 +2211,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 +2246,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 +2264,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 +2310,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 +2333,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 +2359,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/dev/bvbrc_api_prototype.R b/dev/bvbrc_api_prototype.R new file mode 100644 index 0000000..181542d --- /dev/null +++ b/dev/bvbrc_api_prototype.R @@ -0,0 +1,159 @@ +#!/usr/bin/env Rscript +# BV-BRC Data API prototype for amRdata +# ------------------------------------- +# Exploratory/benchmark script — NOT wired into the package namespace. +# Demonstrates the three things that fix the current Docker bottleneck: +# 1. bvbrc_count() -- total matches WITHOUT downloading (Content-Range) +# 2. bvbrc_fetch_all() -- keyset pagination past the 25k cap, optional +# id-space partitioning across a parallel worker pool +# 3. bvbrc_rank_genomes_by_drug_coverage() -- faceted subsampling +# +# Design rationale + limitations: docs/bvbrc-api-feasibility.md +# Requires: httr2, jsonlite, data.table (+ future, future.apply for parallel) + +suppressPackageStartupMessages({ + library(httr2) + library(jsonlite) +}) + +BVBRC_BASE <- "https://www.bv-brc.org/api" +PAGE_MAX <- 25000L # hard per-request AND offset ceiling + +# --- helpers ---------------------------------------------------------------- + +# URL-encode ONE RQL value (BV-BRC wants field names/values encoded individually) +.enc <- function(x) utils::URLencode(as.character(x), reserved = TRUE) + +.bvbrc_req <- function(collection, rql, accept = "application/json") { + url <- sprintf("%s/%s/?%s", BVBRC_BASE, collection, rql) + request(url) |> + req_headers(Accept = accept) |> + req_timeout(120) |> + req_retry( + max_tries = 4L, + is_transient = function(resp) resp_status(resp) %in% c(429, 500, 502, 503, 504) + ) +} + +# --- 1. count() ------------------------------------------------------------- +# Total matching rows without downloading them. Reads "items 0-1/". +bvbrc_count <- function(collection, rql_filter) { + resp <- .bvbrc_req(collection, paste0(rql_filter, "&limit(1)")) |> req_perform() + as.integer(sub(".*/", "", resp_header(resp, "Content-Range"))) +} + +# --- one keyset page -------------------------------------------------------- +# `key` is the unique sort/seek field: "id" for genome_amr/feature/sequence, +# "genome_id" for the genome collection (which has no "id" field). +.bvbrc_page <- function(collection, rql_filter, select, last_id, page_size, + key = "id") { + seek <- if (nzchar(last_id)) { + sprintf("and(%s,gt(%s,%s))", rql_filter, key, .enc(last_id)) + } else { + rql_filter + } + rql <- sprintf("%s&select(%s,%s)&sort(+%s)&limit(%d)", + seek, key, select, key, page_size) + resp <- .bvbrc_req(collection, rql) |> req_perform() + fromJSON(resp_body_string(resp), simplifyDataFrame = TRUE) +} + +# Walk one result set to exhaustion via keyset (seek) pagination. +# gt(key, X) is INCLUSIVE on this API, so we dedupe the boundary row and judge +# termination on the RAW fetched size, not the post-dedupe size. +.keyset_walk <- function(collection, rql_filter, select, page_size = PAGE_MAX, + key = "id") { + out <- list(); last_id <- ""; k <- 0L + repeat { + pg <- .bvbrc_page(collection, rql_filter, select, last_id, page_size, key) + raw_n <- if (is.data.frame(pg)) nrow(pg) else 0L + if (raw_n == 0L) break + if (nzchar(last_id)) pg <- pg[pg[[key]] != last_id, , drop = FALSE] + if (nrow(pg) > 0L) { + k <- k + 1L; out[[k]] <- pg + last_id <- pg[[key]][nrow(pg)] + } + if (raw_n < page_size) break + } + if (k == 0L) return(NULL) + data.table::rbindlist(out, fill = TRUE, use.names = TRUE) +} + +# --- 2. fetch_all() --------------------------------------------------------- +# Pull an entire result set of ANY size. +# parallel = FALSE : single keyset walk (simple, good for interactive subsets) +# parallel = TRUE : partition the id space into hex buckets and keyset-walk +# each in parallel. UUID ids => naturally balanced buckets. +bvbrc_fetch_all <- function(collection, rql_filter, select, + page_size = PAGE_MAX, parallel = FALSE, workers = 6L, + key = "id") { + if (!parallel) { + return(.keyset_walk(collection, rql_filter, select, page_size, key)) + } + if (key != "id") { + stop("parallel id-partitioning requires key = 'id' (UUID collections). ", + "Use parallel = FALSE for the genome collection (key = 'genome_id').") + } + # force promises so the VALUES (not unevaluated args) are exported to workers + force(rql_filter); force(select); force(page_size); force(key) + hex <- c(0:9, letters[1:6]) # '0'..'f' (UUID first char) + edges <- c(hex, "g") # 'g' > 'f' as an exclusive upper bound + buckets <- Map(c, edges[-length(edges)], edges[-1]) + + worker <- function(b) { + filt <- sprintf("and(%s,ge(id,%s),lt(id,%s))", rql_filter, b[1], b[2]) + .keyset_walk(collection, filt, select, page_size, key) + } + + future::plan(future::multisession, workers = workers) + on.exit(future::plan(future::sequential), add = TRUE) + parts <- future.apply::future_lapply( + buckets, worker, + future.packages = c("httr2", "jsonlite", "data.table"), + future.globals = c(".keyset_walk", ".bvbrc_page", ".bvbrc_req", ".enc", + "BVBRC_BASE", "PAGE_MAX", "rql_filter", "select", + "page_size", "key") + ) + data.table::rbindlist(Filter(Negate(is.null), parts), fill = TRUE, use.names = TRUE) +} + +# --- 3. faceted subsampling ------------------------------------------------- +# "Give me the N genomes with the most coverage across these drugs." +# Facet fields come back sorted by count desc, so the top-N buckets ARE the +# best-covered genomes. +bvbrc_rank_genomes_by_drug_coverage <- function(species, antibiotics, n = 500L) { + ab <- paste(vapply(antibiotics, .enc, character(1)), collapse = ",") + rql <- sprintf("and(eq(genome_name,%s),in(antibiotic,(%s)))", .enc(species), ab) + q <- sprintf( + "%s&limit(1)&facet((field,genome_id),(mincount,1),(limit,%d))&json(nl,map)", + rql, n + ) + resp <- .bvbrc_req("genome_amr", q, accept = "application/solr+json") |> req_perform() + ff <- fromJSON(resp_body_string(resp))$facet_counts$facet_fields$genome_id + data.frame(genome_id = names(ff), n_drugs = as.integer(unlist(ff)), + row.names = NULL) +} + +# --- example / benchmark ---------------------------------------------------- +# Call bvbrc_demo() explicitly (sourcing this file no longer auto-runs it). +bvbrc_demo <- function() { + cat("count (S. aureus AMR rows):", + bvbrc_count("genome_amr", "eq(genome_name,Staphylococcus%20aureus)"), "\n") + + t <- system.time( + d <- bvbrc_fetch_all("genome_amr", "eq(genome_name,Staphylococcus%20aureus)", + "genome_id,antibiotic,resistant_phenotype") + ) + cat(sprintf("sequential fetch_all: %d rows in %.1fs\n", nrow(d), t["elapsed"])) + + tp <- system.time( + dp <- bvbrc_fetch_all("genome_amr", "eq(genome_name,Staphylococcus%20aureus)", + "genome_id,antibiotic,resistant_phenotype", + parallel = TRUE, workers = 6L) + ) + cat(sprintf("parallel fetch_all: %d rows in %.1fs\n", nrow(dp), tp["elapsed"])) + + top <- bvbrc_rank_genomes_by_drug_coverage( + "Staphylococcus aureus", c("ciprofloxacin", "gentamicin", "oxacillin"), n = 5L) + cat("top genomes by drug coverage:\n"); print(top) +} diff --git a/dev/bvbrc_pipeline_prototype.R b/dev/bvbrc_pipeline_prototype.R new file mode 100644 index 0000000..2603105 --- /dev/null +++ b/dev/bvbrc_pipeline_prototype.R @@ -0,0 +1,81 @@ +#!/usr/bin/env Rscript +# BV-BRC API pipeline prototype (Phase 1) +# --------------------------------------- +# Standalone demo that the direct-API path is a viable, RESILIENT replacement +# for the Docker/p3-* metadata+AMR download in amRdata (issue #30: stochastic +# BV-BRC downloads). NOT wired into the package. +# +# For a species it: +# 1. pulls Good-quality genome metadata (chosen columns) -- key = genome_id +# 2. pulls AMR phenotype rows (paginated past 25k) -- key = id +# 3. joins on genome_id and writes parquet +# All requests inherit retry-with-backoff on 429/5xx (dev/bvbrc_api_prototype.R), +# so a transient 503 is retried, never parsed as data -> no "more columns" crash. +# +# Requires: dev/bvbrc_api_prototype.R + arrow + +source("dev/bvbrc_api_prototype.R") + +# minimal genome column set (see docs §9); drop the rest +GENOME_COLS <- paste( + "genome_id", "assembly_accession", "genbank_accessions", "genome_quality", + "genome_status", "checkm_completeness", "checkm_contamination", "cds", + "genome_length", "gc_content", "host_name", "isolation_country", + "geographic_group", "species", "taxon_id", sep = "," +) +AMR_COLS <- paste( + "genome_id", "genome_name", "antibiotic", "resistant_phenotype", + "measurement", "measurement_unit", "laboratory_typing_method", "evidence", + sep = ",") + +# Pull Good-quality genome metadata for a species (sequential; genome_id key). +pull_genome_metadata <- function(species, cols = GENOME_COLS, verbose = TRUE) { + filt <- sprintf("and(eq(species,%s),eq(genome_quality,Good))", .enc(species)) + if (verbose) message(" genome metadata: ", bvbrc_count("genome", filt), " Good genomes") + bvbrc_fetch_all("genome", filt, cols, key = "genome_id") +} + +# Pull AMR phenotype rows for a species (parallel; UUID id key). +pull_amr <- function(species, cols = AMR_COLS, parallel = TRUE, verbose = TRUE) { + filt <- sprintf("eq(genome_name,%s)", .enc(species)) + if (verbose) message(" AMR rows: ", bvbrc_count("genome_amr", filt)) + bvbrc_fetch_all("genome_amr", filt, cols, parallel = parallel, key = "id") +} + +# End-to-end: pull both, join, write parquet, return a small summary. +run_pipeline <- function(species, out_dir = "dev/out", parallel_amr = TRUE) { + dir.create(out_dir, showWarnings = FALSE, recursive = TRUE) + message("Species: ", species) + + t_meta <- system.time(meta <- pull_genome_metadata(species)) + t_amr <- system.time(amr <- pull_amr(species, parallel = parallel_amr)) + + # keep AMR only for the Good genomes we pulled, then join + amr_good <- amr[amr$genome_id %in% meta$genome_id, ] + joined <- merge(meta, amr_good, by = "genome_id", suffixes = c("", ".amr")) + + slug <- gsub("[^A-Za-z0-9]+", "_", species) + arrow::write_parquet(meta, file.path(out_dir, paste0(slug, "_genome.parquet"))) + arrow::write_parquet(amr, file.path(out_dir, paste0(slug, "_amr.parquet"))) + arrow::write_parquet(joined, file.path(out_dir, paste0(slug, "_joined.parquet"))) + + summary <- data.frame( + species = species, + good_genomes = nrow(meta), + amr_rows = nrow(amr), + good_genomes_w_amr = length(unique(amr_good$genome_id)), + joined_rows = nrow(joined), + meta_secs = round(t_meta["elapsed"], 1), + amr_secs = round(t_amr["elapsed"], 1), + row.names = NULL + ) + message(" wrote parquet to ", out_dir, "/") + summary +} + +if (!interactive()) { + # Enterococcus faecium: ~7.7k Good genomes (1 page) + ~69k AMR rows (keyset, + # multi-page) -> exercises pagination + parallel without a huge runtime. + s <- run_pipeline("Enterococcus faecium") + cat("\n=== pipeline summary ===\n"); print(s, row.names = FALSE) +} diff --git a/dev/bvbrc_species_roster.R b/dev/bvbrc_species_roster.R new file mode 100644 index 0000000..c9f6aad --- /dev/null +++ b/dev/bvbrc_species_roster.R @@ -0,0 +1,66 @@ +#!/usr/bin/env Rscript +# Per-species roster for the "raise the 25k row cap" request to the bvbrc +# Python package. For each species, count (all via header-only requests): +# total_genomes : genome collection, eq(species,X) +# clean_genomes : + genome_quality = Good (BV-BRC QC; uses CheckM) +# clean_with_amr : + antimicrobial_resistance present +# amr_rows : genome_amr rows, eq(genome_name,) [text match; +# aggregates all strains -- genome_amr.taxon_id is unreliable +# (mixes species- and strain-level taxa), so name match is +# the robust per-species key] +# and flag whether any of those exceed the 25,000 per-request ceiling. +# +# See docs/bvbrc-api-feasibility.md. Requires dev/bvbrc_api_prototype.R. + +source("dev/bvbrc_api_prototype.R") + +CAP <- 25000L + +bvbrc_species_row <- function(species) { + s <- .enc(species) + total <- bvbrc_count("genome", sprintf("eq(species,%s)", s)) + clean <- bvbrc_count("genome", + sprintf("and(eq(species,%s),eq(genome_quality,Good))", s)) + clean_amr <- bvbrc_count("genome", + sprintf("and(eq(species,%s),eq(genome_quality,Good),eq(antimicrobial_resistance,*))", s)) + amr_rows <- bvbrc_count("genome_amr", sprintf("eq(genome_name,%s)", s)) + data.frame(species = species, total_genomes = total, + clean_genomes = clean, clean_with_amr = clean_amr, + amr_rows = amr_rows, + genomes_over_cap = total > CAP, amr_rows_over_cap = amr_rows > CAP) +} + +bvbrc_roster <- function(species_vec) { + rows <- lapply(species_vec, function(sp) { + message(" ", sp); tryCatch(bvbrc_species_row(sp), + error = function(e) { message(" ! ", conditionMessage(e)); NULL }) + }) + do.call(rbind, Filter(Negate(is.null), rows)) +} + +# WHO Bacterial Priority Pathogens List (2024) + CDC AR Threats (2019), bacterial +WHO_CDC_PRIORITY <- c( + "Acinetobacter baumannii", "Pseudomonas aeruginosa", "Escherichia coli", + "Klebsiella pneumoniae", "Enterobacter cloacae", "Staphylococcus aureus", + "Enterococcus faecium", "Enterococcus faecalis", "Streptococcus pneumoniae", + "Salmonella enterica", "Shigella flexneri", "Shigella sonnei", + "Neisseria gonorrhoeae", "Neisseria meningitidis", "Campylobacter jejuni", + "Campylobacter coli", "Haemophilus influenzae", "Helicobacter pylori", + "Mycobacterium tuberculosis", "Clostridioides difficile", + "Streptococcus pyogenes", "Streptococcus agalactiae", "Serratia marcescens", + "Proteus mirabilis", "Morganella morganii" +) + +if (!interactive()) { + r <- bvbrc_roster(WHO_CDC_PRIORITY) + r <- r[order(-r$amr_rows), ] + out <- "dev/bvbrc_species_roster.csv" + write.csv(r, out, row.names = FALSE) + cat("\n=== Per-species roster (WHO/CDC priority) ===\n") + print(r, row.names = FALSE) + cat(sprintf("\nSpecies with >25k genomes: %d / %d\n", + sum(r$genomes_over_cap), nrow(r))) + cat(sprintf("Species with >25k genome_amr rows: %d / %d\n", + sum(r$amr_rows_over_cap), nrow(r))) + cat("Saved:", out, "\n") +} diff --git a/dev/bvbrc_species_roster.csv b/dev/bvbrc_species_roster.csv new file mode 100644 index 0000000..9927cd0 --- /dev/null +++ b/dev/bvbrc_species_roster.csv @@ -0,0 +1,26 @@ +"species","total_genomes","clean_genomes","clean_with_amr","amr_rows","genomes_over_cap","amr_rows_over_cap" +"Escherichia coli",118185,97915,8623,7388629,TRUE,TRUE +"Mycobacterium tuberculosis",46979,42536,39253,2257046,TRUE,TRUE +"Klebsiella pneumoniae",42836,33618,32134,2006673,TRUE,TRUE +"Salmonella enterica",43005,36716,7285,1969740,TRUE,TRUE +"Streptococcus pneumoniae",32218,29630,25855,1060510,TRUE,TRUE +"Staphylococcus aureus",29360,25585,19069,590284,TRUE,TRUE +"Acinetobacter baumannii",17684,15112,13777,515917,FALSE,TRUE +"Campylobacter jejuni",39644,28123,96,478916,TRUE,TRUE +"Pseudomonas aeruginosa",16370,13282,11585,328979,FALSE,TRUE +"Neisseria gonorrhoeae",15004,13962,5527,264833,FALSE,TRUE +"Shigella sonnei",6381,6135,1965,128515,FALSE,TRUE +"Enterococcus faecium",9645,7742,2920,69474,FALSE,TRUE +"Campylobacter coli",6457,5142,140,55414,FALSE,TRUE +"Enterobacter cloacae",2822,2474,273,31616,FALSE,TRUE +"Clostridioides difficile",4622,4008,827,16580,FALSE,FALSE +"Shigella flexneri",5988,5487,379,14933,FALSE,FALSE +"Serratia marcescens",2219,1869,25,4205,FALSE,FALSE +"Enterococcus faecalis",6117,5043,137,3499,FALSE,FALSE +"Neisseria meningitidis",2988,2778,337,1867,FALSE,FALSE +"Proteus mirabilis",1782,1532,22,1783,FALSE,FALSE +"Morganella morganii",668,469,4,1432,FALSE,FALSE +"Streptococcus agalactiae",2864,2456,160,1039,FALSE,FALSE +"Haemophilus influenzae",1104,953,87,705,FALSE,FALSE +"Streptococcus pyogenes",2868,2624,1,294,FALSE,FALSE +"Helicobacter pylori",3367,3142,54,271,FALSE,FALSE 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..4eae20a --- /dev/null +++ b/docs/bvbrc-api-feasibility.md @@ -0,0 +1,302 @@ +# 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 + +A runnable prototype lives at [`dev/bvbrc_api_prototype.R`](../dev/bvbrc_api_prototype.R) +(standalone; not wired into the package). It implements `bvbrc_count()`, +`bvbrc_fetch_all()` (keyset + optional parallel id-partitioning), and +`bvbrc_rank_genomes_by_drug_coverage()`. Run `bvbrc_demo()` after sourcing. + +**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. + +Roster produced by [`dev/bvbrc_species_roster.R`](../dev/bvbrc_species_roster.R) +→ [`dev/bvbrc_species_roster.csv`](../dev/bvbrc_species_roster.csv). Counts are +header-only (`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. Full table in the CSV. + +**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..f46bd21 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, @@ -42,6 +43,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 (\%).} 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) +})