From ee9bfcb3619e60d31f55f8a4236fe7a03e87ff90 Mon Sep 17 00:00:00 2001
From: Evan Pierce Brenner <108823789+epbrenner@users.noreply.github.com>
Date: Wed, 22 Jul 2026 11:17:14 -0600
Subject: [PATCH 01/14] Add more genome quality filters to data_curation.R

Added ability to filter BV-BRC genomes based on genome stat metadata, including CheckM completeness and CheckM contamination as defaults. Also enabled optional checking of genome length, CDS count, and GC content based on deviations from median distribution per bug. Updated genome downloader functions to use ftpes_download_one instead of old ftp_download_one, and moved a dot function helper to inside the ftpes_download_one loop to avoid issues with BiocParallel workers not understanding what that helper is.
---
 R/data_curation.R | 518 ++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 428 insertions(+), 90 deletions(-)

diff --git a/R/data_curation.R b/R/data_curation.R
index f056c37..acd6bd1 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -36,17 +36,26 @@
 .ftpes_download_one <- function(genomeID, out_dir,
                                 connect_timeout = 10L,
                                 max_time = 30L,
-                                speed_time = 30L, # end if <speed_limit for >speed_time
-                                speed_limit = 2048L, # B/s
+                                speed_time = 30L,
+                                speed_limit = 2048L,
                                 min_bytes = 100L) {
   dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
 
+  # Adding this prior dot function helper into this loop for parallel workers
+  is_complete <- function(dir, genomeID, min_bytes = 100L) {
+    fna <- file.path(dir, paste0(genomeID, ".fna"))
+    faa <- file.path(dir, paste0(genomeID, ".PATRIC.faa"))
+    gff <- file.path(dir, paste0(genomeID, ".PATRIC.gff"))
+    paths <- c(fna, faa, gff)
+    all(file.exists(paths)) &&
+      all(vapply(paths, function(x) file.info(x)$size, numeric(1)) > min_bytes)
+  }
+
   exts <- c(".fna", ".PATRIC.faa", ".PATRIC.gff")
   for (ext in exts) {
     dest <- file.path(out_dir, paste0(genomeID, ext))
     dest_tmp <- paste0(dest, ".tmp")
 
-    # Skip if we already completed these
     if (file.exists(dest) && file.info(dest)$size > min_bytes) next
     if (file.exists(dest) && file.info(dest)$size == 0) try(unlink(dest), silent = TRUE)
     if (file.exists(dest_tmp)) try(unlink(dest_tmp), silent = TRUE)
@@ -59,7 +68,7 @@
       "--max-time", as.character(max_time),
       "--speed-time", as.character(speed_time),
       "--speed-limit", as.character(speed_limit),
-      "--ftp-ssl", # AUTH TLS on port 21 works in testing
+      "--ftp-ssl",
       "--ftp-pasv", "--disable-epsv", "--ipv4",
       "--user", "anonymous:",
       "-o", shQuote(dest_tmp), shQuote(url)
@@ -69,22 +78,19 @@
     status <- attr(res, "status")
     if (is.null(status)) status <- 0L
 
-    # Avoiding 0B files cluttering up results -- atomic rename on successful tmp DL
     ok <- (status == 0L && file.exists(dest_tmp) && file.info(dest_tmp)$size > min_bytes)
     if (ok) {
       if (!file.rename(dest_tmp, dest)) {
-        # If rename fails for some reason, copy + unlink
         file.copy(dest_tmp, dest, overwrite = TRUE)
         unlink(dest_tmp)
       }
     } else {
       try(unlink(dest_tmp), silent = TRUE)
-      return(FALSE) # Abort early, don't accept partial set
+      return(FALSE)
     }
   }
 
-  # Do we have all 3 present?
-  .is_complete_set(out_dir, genomeID, min_bytes = min_bytes)
+  is_complete(out_dir, genomeID, min_bytes = min_bytes)
 }
 
 #' Helps manage FTPS downloading from BV-BRC, tryng a quick download first, and
@@ -234,6 +240,230 @@
   return(df)
 }
 
+#' Parse BV-BRC TSV output
+#'
+#' Safely cleans blank links and weird formatting for a safe merge of the BV-BRC
+#' TSV data into a tibble.
+#'
+#' @param x Character vector of output lines from a BV-BRC CLI command.
+#'
+#' @return A tibble with minimal name repair.
+#' @keywords internal
+.parse_bvbrc_tsv <- function(x) {
+  if (!length(x)) return(tibble::tibble())
+
+  x <- as.character(x)
+  x <- x[nzchar(trimws(x))]
+
+  # Keep only lines that look like TSV records/headers
+  x <- x[grepl("\t", x, fixed = TRUE)]
+
+  if (!length(x)) return(tibble::tibble())
+
+  txt <- paste(x, collapse = "\n")
+  df <- utils::read.table(
+    text = txt,
+    sep = "\t",
+    header = TRUE,
+    fill = TRUE,
+    quote = "",
+    check.names = FALSE,
+    comment.char = "",
+    colClasses = "character"
+  )
+
+  tibble::as_tibble(df, .name_repair = "minimal")
+}
+
+
+
+#' Apply QC filters to BV-BRC genome metadata
+#'
+#' Quality control step with default values for contamination and completeness.
+#' Optional filters for other genome stats to find genomes that are X deviations
+#' away from the median values for those stats. This shouldn't be applied to
+#' jobs that use more than one diverse taxon at a time though! Requires a fairly
+#' normal distribution like you'd see for (most) single species
+#' Returns the cleaned table, retained genome IDs, and a rejection log.
+#'
+#' CheckM contamination and completeness are treated as the primary QC gate.
+#' Deviation-based filters are optional and are best suited to single-taxon runs.
+#'
+#' @param genome_tbl A tibble of genome metadata containing BV-BRC genome columns.
+#' @param checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%).
+#' @param checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%).
+#' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content.
+#' @param length_deviations Optional numeric scalar. Maximum SDs from the median genome length.
+#' @param cds_deviations Optional numeric scalar. Maximum SDs from the median CDS count.
+#'
+#' @return A list with components:
+#' \describe{
+#'   \item{qc_tbl}{QC-annotated genome table}
+#'   \item{keep_ids}{Character vector of retained genome IDs}
+#'   \item{rejections}{Tibble describing dropped genomes and reasons}
+#' }
+#' @keywords internal
+.apply_metadata_qc <- function(genome_tbl,
+                               checkm_contam = 10,
+                               checkm_complete = 90,
+                               gc_deviations = NULL,
+                               length_deviations = NULL,
+                               cds_deviations = NULL) {
+  numify <- function(x) suppressWarnings(as.numeric(x))
+
+  safe_median <- function(x) {
+    x <- x[is.finite(x)]
+    if (!length(x)) return(NA_real_)
+    stats::median(x)
+  }
+
+  safe_sd <- function(x) {
+    x <- x[is.finite(x)]
+    if (length(x) < 2L) return(NA_real_)
+    stats::sd(x)
+  }
+
+  qc <- tibble::as_tibble(genome_tbl) |>
+    dplyr::mutate(
+      genome_length_qc = numify(.data$`genome.genome_length`),
+      gc_content_qc = numify(.data$`genome.gc_content`),
+      cds_qc = numify(.data$`genome.cds`),
+      checkm_completeness_qc = numify(.data$`genome.checkm_completeness`),
+      checkm_contamination_qc = numify(.data$`genome.checkm_contamination`)
+    )
+
+  med_len <- safe_median(qc$genome_length_qc)
+  sd_len <- safe_sd(qc$genome_length_qc)
+  med_gc <- safe_median(qc$gc_content_qc)
+  sd_gc <- safe_sd(qc$gc_content_qc)
+  med_cds <- safe_median(qc$cds_qc)
+  sd_cds <- safe_sd(qc$cds_qc)
+
+  z_len <- if (!is.null(length_deviations) && is.finite(sd_len) && sd_len > 0) {
+    abs(qc$genome_length_qc - med_len) / sd_len
+  } else {
+    rep(NA_real_, nrow(qc))
+  }
+
+  z_gc <- if (!is.null(gc_deviations) && is.finite(sd_gc) && sd_gc > 0) {
+    abs(qc$gc_content_qc - med_gc) / sd_gc
+  } else {
+    rep(NA_real_, nrow(qc))
+  }
+
+  z_cds <- if (!is.null(cds_deviations) && is.finite(sd_cds) && sd_cds > 0) {
+    abs(qc$cds_qc - med_cds) / sd_cds
+  } else {
+    rep(NA_real_, nrow(qc))
+  }
+
+  qc <- qc |>
+    dplyr::mutate(
+      flag_checkm_missing = is.na(.data$checkm_completeness_qc) | is.na(.data$checkm_contamination_qc),
+      flag_checkm_contam = !is.na(.data$checkm_contamination_qc) & .data$checkm_contamination_qc > checkm_contam,
+      flag_checkm_complete = !is.na(.data$checkm_completeness_qc) & .data$checkm_completeness_qc < checkm_complete,
+      flag_checkm = .data$flag_checkm_missing | .data$flag_checkm_contam | .data$flag_checkm_complete,
+      flag_length = if (!is.null(length_deviations) && is.finite(sd_len) && sd_len > 0) {
+        !is.na(z_len) & z_len > length_deviations
+      } else FALSE,
+      flag_gc = if (!is.null(gc_deviations) && is.finite(sd_gc) && sd_gc > 0) {
+        !is.na(z_gc) & z_gc > gc_deviations
+      } else FALSE,
+      flag_cds = if (!is.null(cds_deviations) && is.finite(sd_cds) && sd_cds > 0) {
+        !is.na(z_cds) & z_cds > cds_deviations
+      } else FALSE,
+      qc_keep = !(.data$flag_checkm | .data$flag_length | .data$flag_gc | .data$flag_cds)
+    )
+
+  make_rej <- function(flag, rule, observed, threshold, comparator) {
+    idx <- which(flag %in% TRUE)
+    if (!length(idx)) return(NULL)
+    tibble::tibble(
+      genome_id = qc$`genome.genome_id`[idx],
+      genome_name = qc$`genome.genome_name`[idx],
+      species = qc$`genome.species`[idx],
+      failed_rule = rule,
+      observed = as.character(observed[idx]),
+      threshold = as.character(threshold),
+      comparator = comparator,
+      action_taken = "dropped"
+    )
+  }
+
+  rej_list <- list(
+    make_rej(
+      qc$flag_checkm_missing,
+      "checkm_missing",
+      ifelse(qc$flag_checkm_missing, "missing", NA_character_),
+      "present",
+      "is present"
+    ),
+    make_rej(
+      qc$flag_checkm_contam,
+      "checkm_contamination",
+      qc$checkm_contamination_qc,
+      checkm_contam,
+      ">"
+    ),
+    make_rej(
+      qc$flag_checkm_complete,
+      "checkm_completeness",
+      qc$checkm_completeness_qc,
+      checkm_complete,
+      "<"
+    )
+  )
+
+  if (!is.null(gc_deviations)) {
+    rej_list[[length(rej_list) + 1L]] <- make_rej(
+      qc$flag_gc,
+      "gc_deviation_sd",
+      z_gc,
+      gc_deviations,
+      ">"
+    )
+  }
+
+  if (!is.null(length_deviations)) {
+    rej_list[[length(rej_list) + 1L]] <- make_rej(
+      qc$flag_length,
+      "length_deviation_sd",
+      z_len,
+      length_deviations,
+      ">"
+    )
+  }
+
+  if (!is.null(cds_deviations)) {
+    rej_list[[length(rej_list) + 1L]] <- make_rej(
+      qc$flag_cds,
+      "cds_deviation_sd",
+      z_cds,
+      cds_deviations,
+      ">"
+    )
+  }
+
+  rejections <- dplyr::bind_rows(rej_list)
+  if (nrow(rejections) == 0L) {
+    rejections <- tibble::tibble(
+      genome_id = character(),
+      genome_name = character(),
+      species = character(),
+      failed_rule = character(),
+      observed = character(),
+      threshold = character(),
+      comparator = character(),
+      action_taken = character()
+    )
+  }
+
+  list(
+    qc_tbl = qc,
+    keep_ids = unique(qc$`genome.genome_id`[qc$qc_keep]),
+    rejections = rejections
+  )
+}
 
 # Make sure the BV-BRC metadata live where they're supposed to
 .ensure_bvbrc_cache <- function(base_dir = ".",
@@ -910,7 +1140,13 @@
 #' @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 verbose Logical. If TRUE, prints concise messages.
+#' @param checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%).
+#' @param checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%).
+#' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content.
+#' @param length_deviations Optional numeric scalar. Maximum SDs from the median genome length.
+#' @param cds_deviations Optional numeric scalar. Maximum SDs from the median CDS count.
+#' @param debug Logical. If TRUE, retain `metadata_full` and QC columns for inspection.
+#' @param verbose Logical. If TRUE, print progress messages.
 #'
 #' @return A list with:
 #'   - duckdbConnection: live DBI connection to the created DuckDB
@@ -923,27 +1159,27 @@ retrieveMetadata <- function(user_bacs,
                              abx = "All",
                              overwrite = FALSE,
                              image = "danylmb/bvbrc:5.3",
+                             checkm_contam = 10,
+                             checkm_complete = 90,
+                             gc_deviations = NULL,
+                             length_deviations = NULL,
+                             cds_deviations = NULL,
+                             debug = FALSE,
                              verbose = TRUE) {
   base_dir <- normalizePath(base_dir, mustWork = FALSE)
 
-  # -------------------------------
-  # GENOME ID RESOLUTION (UPDATED)
-  # -------------------------------
   if (!is.null(genome_id_file)) {
     if (!file.exists(genome_id_file)) {
       stop("Provided genome_id_file does not exist.")
     }
-
     if (isTRUE(verbose)) {
       message("Using genome IDs from file: ", genome_id_file)
     }
-
     genome_ids <- readLines(genome_id_file, warn = FALSE)
     genome_ids <- trimws(genome_ids)
     genome_ids <- genome_ids[genome_ids != ""]
   } else {
     if (isTRUE(verbose)) message("Resolving genome IDs for user inputs.")
-
     genome_ids <- .retrieveQueryIDs(
       base_dir = base_dir,
       user_bacs = user_bacs,
@@ -953,13 +1189,11 @@ retrieveMetadata <- function(user_bacs,
   }
 
   genome_ids <- unique(as.character(genome_ids))
-
-  if (length(genome_ids) == 0) {
+  if (length(genome_ids) == 0L) {
     message("No genome IDs available for the specified inputs.")
     return(NULL)
   }
 
-  # Desired fields from BV-BRC
   drug_fields <- paste0(
     "antibiotic,computational_method,",
     "evidence,genome_name,id,",
@@ -970,11 +1204,13 @@ retrieveMetadata <- function(user_bacs,
     "pmid,resistant_phenotype,",
     "source,taxon_id,testing_standard"
   )
+
   abx_filter <- if (identical(abx, "All")) {
     "--required antibiotic"
   } else {
     paste0("--in antibiotic,", paste(abx, collapse = ","))
   }
+
   amr_fields <- paste0(
     "assembly_accession,assembly_method,",
     "bioproject_accession,biosample_accession,",
@@ -991,6 +1227,7 @@ retrieveMetadata <- function(user_bacs,
     "refseq_accessions,",
     "refseq_project_id,sra_accession,species,taxon_id"
   )
+
   microtrait_fields <- paste0(
     amr_fields, ",genome_length,gram_stain,",
     "habitat,",
@@ -1009,121 +1246,202 @@ retrieveMetadata <- function(user_bacs,
     "trna"
   )
 
-  # Batching downloads and parallel implementation
+  qc_fields <- c(
+    "genome_length",
+    "gc_content",
+    "contigs",
+    "cds",
+    "rrna",
+    "trna",
+    "checkm_completeness",
+    "checkm_contamination"
+  )
+
+  amr_fields <- paste(
+    unique(c(strsplit(amr_fields, ",", fixed = TRUE)[[1]], qc_fields)),
+    collapse = ","
+  )
+  microtrait_fields <- paste(
+    unique(c(strsplit(microtrait_fields, ",", fixed = TRUE)[[1]], qc_fields)),
+    collapse = ","
+  )
+
   batch_size <- 500L
-  genome_ids <- as.character(genome_ids)
   genome_batches <- split(genome_ids, ceiling(seq_along(genome_ids) / batch_size))
 
   n_cores <- max(1L, parallel::detectCores(logical = TRUE) - 1L)
   param <- BiocParallel::SnowParam(workers = n_cores)
 
-  # Pull AMR metadata
   if (isTRUE(verbose)) message("Retrieving AMR phenotype data in batches.")
   batch_drug_data <- BiocParallel::bplapply(genome_batches, function(batch) {
-    .extractAMRtable(
-      base_dir          = base_dir,
-      batch_genome_IDs  = batch,
-      abx_filter        = abx_filter,
-      drug_fields       = drug_fields,
-      image             = image,
-      verbose           = FALSE
+    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)
   }, BPPARAM = param)
-  combined_drug_data <- unlist(batch_drug_data, use.names = FALSE)
-  if (length(combined_drug_data) == 0) {
+
+  combined_drug_data_tbl <- dplyr::bind_rows(batch_drug_data) |>
+    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)
   }
-  combined_drug_data_tbl <- tibble::as_tibble(utils::read.table(
-    text = combined_drug_data,
-    sep = "\t", header = TRUE, fill = TRUE,
-    quote = "", check.names = FALSE, comment.char = "", colClasses = "character"
-  )) |>
-    dplyr::mutate(dplyr::across(dplyr::everything(), ~ iconv(.x, from = "", to = "UTF-8", sub = ""))) |>
-    dplyr::mutate(`genome_drug.genome_id` = as.character(`genome_drug.genome_id`))
 
-  # Pull genome metadata
   if (isTRUE(verbose)) message("Retrieving genome metadata in batches.")
   batch_genome_data <- BiocParallel::bplapply(genome_batches, function(batch) {
-    .extractGenomeData(
-      base_dir          = base_dir,
-      batch_genome_IDs  = batch,
-      filter_type       = filter_type,
-      amr_fields        = amr_fields,
+    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
+      image = image,
+      verbose = FALSE
     )
+    .parse_bvbrc_tsv(raw)
   }, BPPARAM = param)
-  combined_genome_data <- unlist(batch_genome_data, use.names = FALSE)
-  if (length(combined_genome_data) == 0) {
+
+  combined_genome_data_tbl <- dplyr::bind_rows(batch_genome_data) |>
+    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)
   }
-  combined_genome_data_tbl <- tibble::as_tibble(utils::read.table(
-    text = combined_genome_data,
-    sep = "\t", header = TRUE, fill = TRUE,
-    quote = "", check.names = FALSE, comment.char = "", colClasses = "character"
-  )) |>
-    dplyr::mutate(dplyr::across(dplyr::everything(), ~ iconv(.x, from = "", to = "UTF-8", sub = ""))) |>
-    dplyr::mutate(`genome.genome_id` = as.character(`genome.genome_id`))
 
-  # Write & join
+  if ("genome.genome_id" %in% names(combined_genome_data_tbl)) {
+    combined_genome_data_tbl <- combined_genome_data_tbl |>
+      dplyr::filter(grepl("^[0-9]+\\.[0-9]+$", .data$`genome.genome_id`)) |>
+      dplyr::distinct(`genome.genome_id`, .keep_all = TRUE)
+  }
+
+  if ("genome_drug.genome_id" %in% names(combined_drug_data_tbl)) {
+    combined_drug_data_tbl <- combined_drug_data_tbl |>
+      dplyr::filter(grepl("^[0-9]+\\.[0-9]+$", .data$`genome_drug.genome_id`))
+  }
+
+  qc_out <- .apply_metadata_qc(
+    genome_tbl = combined_genome_data_tbl,
+    checkm_contam = checkm_contam,
+    checkm_complete = checkm_complete,
+    gc_deviations = gc_deviations,
+    length_deviations = length_deviations,
+    cds_deviations = cds_deviations
+  )
+
+  keep_ids <- qc_out$keep_ids
+  if (!length(keep_ids)) {
+    message("No genomes passed QC.")
+    return(NULL)
+  }
+
+  combined_genome_data_tbl <- combined_genome_data_tbl |>
+    dplyr::filter(`genome.genome_id` %in% keep_ids)
+
+  combined_drug_data_tbl <- combined_drug_data_tbl |>
+    dplyr::filter(`genome_drug.genome_id` %in% keep_ids)
+
   paths <- .buildDBpath(base_dir = base_dir, user_bacs = user_bacs, overwrite = overwrite)
   db_path <- paths$db_path
   logs_dir <- file.path(base_dir, "data", "logs")
   dir.create(logs_dir, recursive = TRUE, showWarnings = FALSE)
+
   cat(sprintf("[%s] Writing metadata DuckDB: %s\n", Sys.time(), db_path),
-    file = file.path(logs_dir, "bvbrc.log"), append = TRUE
-  )
+      file = file.path(logs_dir, "bvbrc.log"), append = TRUE)
 
   con <- DBI::dbConnect(duckdb::duckdb(), dbdir = db_path)
   on.exit(try(DBI::dbDisconnect(con, shutdown = TRUE), silent = TRUE), add = TRUE)
 
   DBI::dbWriteTable(con, "amr_phenotype", combined_drug_data_tbl, overwrite = TRUE)
   DBI::dbWriteTable(con, "genome_data", combined_genome_data_tbl, overwrite = TRUE)
+  DBI::dbWriteTable(con, "metadata_qc", qc_out$qc_tbl, overwrite = TRUE)
+  DBI::dbWriteTable(con, "metadata_qc_rejections", qc_out$rejections, overwrite = TRUE)
 
   if (isTRUE(verbose)) message("Joining AMR phenotype and genome metadata.")
   DBI::dbExecute(con, '
-    CREATE OR REPLACE TABLE metadata AS
+    CREATE OR REPLACE TABLE metadata_full AS
     SELECT *
     FROM amr_phenotype
     INNER JOIN genome_data
       ON amr_phenotype."genome_drug.genome_id" = genome_data."genome.genome_id"
   ')
 
-  # Debug summary after writes
+  qc_drop_cols <- c(
+    "genome.genome_length",
+    "genome.gc_content",
+    "genome.contigs",
+    "genome.cds",
+    "genome.rrna",
+    "genome.trna",
+    "genome.checkm_completeness",
+    "genome.checkm_contamination"
+  )
+
+  if (isTRUE(debug)) {
+    DBI::dbExecute(con, 'CREATE OR REPLACE TABLE metadata AS SELECT * FROM metadata_full')
+  } else {
+    keep_cols <- setdiff(DBI::dbListFields(con, "metadata_full"), qc_drop_cols)
+    select_sql <- paste(DBI::dbQuoteIdentifier(con, keep_cols), collapse = ", ")
+    DBI::dbExecute(
+      con,
+      paste0("CREATE OR REPLACE TABLE metadata AS SELECT ", select_sql, " FROM metadata_full")
+    )
+    DBI::dbExecute(con, "DROP TABLE metadata_full")
+  }
+
   n_targets <- length(genome_ids)
-  n_amr_ids <- DBI::dbGetQuery(con, 'SELECT COUNT(DISTINCT "genome_drug.genome_id") AS n FROM amr_phenotype')$n
-  n_gmeta_ids <- DBI::dbGetQuery(con, 'SELECT COUNT(DISTINCT "genome.genome_id")      AS n FROM genome_data')$n
+  n_amr_ids <- DBI::dbGetQuery(con, '
+    SELECT COUNT(DISTINCT "genome_drug.genome_id") AS n FROM amr_phenotype
+  ')$n
+  n_gmeta_ids <- DBI::dbGetQuery(con, '
+    SELECT COUNT(DISTINCT "genome.genome_id") AS n FROM genome_data
+  ')$n
+
   if (isTRUE(verbose)) {
     message(
       "Initial summary: targets=", n_targets,
       " | AMR genomes=", n_amr_ids,
       " | genome_data genomes=", n_gmeta_ids
     )
+    message(
+      "QC summary: kept=", length(keep_ids),
+      " | dropped=", nrow(qc_out$rejections)
+    )
   }
 
-  # Tagged view
   .create_amr_tagged_view(con)
 
-  # Final debug summary
   n_bac <- if ("bac_data" %in% DBI::dbListTables(con)) {
-    DBI::dbGetQuery(con, 'SELECT COUNT(DISTINCT "genome.genome_id") AS n FROM bac_data')$n
+    DBI::dbGetQuery(con, '
+      SELECT COUNT(DISTINCT "genome.genome_id") AS n FROM bac_data
+    ')$n
   } else {
     NA_integer_
   }
-  n_amr_ids <- DBI::dbGetQuery(con, 'SELECT COUNT(DISTINCT "genome_drug.genome_id") AS n FROM amr_phenotype')$n
-  n_gmeta_ids <- DBI::dbGetQuery(con, 'SELECT COUNT(DISTINCT "genome.genome_id")      AS n FROM genome_data')$n
+
+  n_amr_ids <- DBI::dbGetQuery(con, '
+    SELECT COUNT(DISTINCT "genome_drug.genome_id") AS n FROM amr_phenotype
+  ')$n
+  n_gmeta_ids <- DBI::dbGetQuery(con, '
+    SELECT COUNT(DISTINCT "genome.genome_id") AS n FROM genome_data
+  ')$n
 
   ids_zero_amr <- DBI::dbGetQuery(
     con,
-    'WITH sel AS (SELECT DISTINCT "genome.genome_id" AS gid FROM genome_data),
-          got AS (SELECT DISTINCT "genome_drug.genome_id" AS gid FROM amr_phenotype)
-     SELECT sel.gid
-     FROM sel LEFT JOIN got USING (gid)
-     WHERE got.gid IS NULL
-     ORDER BY sel.gid'
+    '
+    WITH sel AS (SELECT DISTINCT "genome.genome_id" AS gid FROM genome_data),
+         got AS (SELECT DISTINCT "genome_drug.genome_id" AS gid FROM amr_phenotype)
+    SELECT sel.gid
+    FROM sel LEFT JOIN got USING (gid)
+    WHERE got.gid IS NULL
+    ORDER BY sel.gid
+    '
   )$gid
 
   if (isTRUE(verbose)) {
@@ -1357,16 +1675,18 @@ retrieveMetadata <- function(user_bacs,
   list(duckdbConnection = con, table_name = "filtered")
 }
 
+### This has been deprecated after being moved into the loop for .ftpes_download_one()
+### Can remove once testing is completed
 #' Helps check if a complete set exists after DL (.fna + .PATRIC.faa + .PATRIC.gff)
 #' @keywords internal
-.is_complete_set <- function(dir, genomeID, min_bytes = 100) {
-  fna <- file.path(dir, paste0(genomeID, ".fna"))
-  faa <- file.path(dir, paste0(genomeID, ".PATRIC.faa"))
-  gff <- file.path(dir, paste0(genomeID, ".PATRIC.gff"))
-  paths <- c(fna, faa, gff)
-  all(file.exists(paths)) &&
-    all(purrr::map_dbl(paths, function(x) file.info(x)$size) > min_bytes)
-}
+#.is_complete_set <- function(dir, genomeID, min_bytes = 100) {
+#  fna <- file.path(dir, paste0(genomeID, ".fna"))
+#  faa <- file.path(dir, paste0(genomeID, ".PATRIC.faa"))
+#  gff <- file.path(dir, paste0(genomeID, ".PATRIC.gff"))
+#  paths <- c(fna, faa, gff)
+#  all(file.exists(paths)) &&
+#    all(purrr::map_dbl(paths, function(x) file.info(x)$size) > min_bytes)
+#}
 
 #' Helps collate completed genomes into a set
 #' @keywords internal
@@ -1619,7 +1939,7 @@ retrieveGenomes <- function(base_dir = ".",
   if (identical(method, "ftp")) {
     if (isTRUE(verbose)) message("Trying FTPS download. Workers=", ftp_workers)
     ftp_param <- BiocParallel::SnowParam(workers = max(1L, ftp_workers))
-    ft_ok <- BiocParallel::bplapply(ids, function(gid) .ftp_download_one(gid, genome_path),
+    ft_ok <- BiocParallel::bplapply(ids, function(gid) .ftpes_download_one(gid, genome_path),
       BPPARAM = ftp_param
     )
     ok_ids <- ids[unlist(ft_ok)]
@@ -1763,7 +2083,7 @@ genomeList <- function(base_dir = ".",
 #'   <base_dir>/data/<bug_dir>/<abbrev>.duckdb
 #'
 #' @param user_bacs Character vector. Species and/or taxon IDs (e.g.
-#'   `c("Shigella flexneri", "623")`).
+#'   `c("Shigella sonnei", "624")`).
 #' @param genome_id_file Character or NULL. Optional path to a file listing genome
 #'   IDs (one per line), passed through to `retrieveMetadata()`. If provided, the
 #'   metadata step is restricted to these genome IDs instead of resolving them from
@@ -1776,6 +2096,12 @@ genomeList <- function(base_dir = ".",
 #' @param evidence_mode Character. Sets what types of AMR evidence is acceptable.
 #'    Default `lab_only`. `any` will not require AMR data for downloads. This will
 #'    return very large download lists for many species!
+#' @param checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%).
+#' @param checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%).
+#' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content.
+#' @param length_deviations Optional numeric scalar. Maximum SDs from the median genome length.
+#' @param cds_deviations Optional numeric scalar. Maximum SDs from the median CDS count.
+#' @param debug Logical. If TRUE, keep QC columns in metadata tables.
 #' @param verbose Logical. Print progress messages. Default TRUE.
 #'
 #' @return A list (the output of `genomeList()`), containing:
@@ -1789,6 +2115,12 @@ prepareGenomes <- function(user_bacs,
                            method = c("ftp", "cli"),
                            overwrite = FALSE,
                            evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"),
+                           checkm_contam = 10,
+                           checkm_complete = 90,
+                           gc_deviations = NULL,
+                           length_deviations = NULL,
+                           cds_deviations = NULL,
+                           debug = FALSE,
                            verbose = TRUE) {
   method <- match.arg(method)
   evidence_mode <- match.arg(evidence_mode)
@@ -1804,6 +2136,12 @@ prepareGenomes <- function(user_bacs,
     base_dir = base_dir,
     abx = "All",
     overwrite = overwrite,
+    checkm_contam = checkm_contam,
+    checkm_complete = checkm_complete,
+    gc_deviations = gc_deviations,
+    length_deviations = length_deviations,
+    cds_deviations = cds_deviations,
+    debug = debug,
     verbose = verbose
   ))
 
@@ -1820,13 +2158,13 @@ prepareGenomes <- function(user_bacs,
     return(NULL)
   }
 
-  # A little summary of what's left after filtering (or not)
   paths <- .buildDBpath(base_dir = base_dir, user_bacs = user_bacs)
   con <- DBI::dbConnect(duckdb::duckdb(), dbdir = paths$db_path, read_only = TRUE)
   n_filtered <- DBI::dbGetQuery(con, 'SELECT COUNT(DISTINCT "genome.genome_id") AS n FROM filtered')$n
   n_meta <- if ("genome_data" %in% DBI::dbListTables(con)) DBI::dbGetQuery(con, 'SELECT COUNT(DISTINCT "genome.genome_id") AS n FROM genome_data')$n else NA
   n_amr <- if ("amr_phenotype" %in% DBI::dbListTables(con)) DBI::dbGetQuery(con, 'SELECT COUNT(DISTINCT "genome_drug.genome_id") AS n FROM amr_phenotype')$n else NA
   DBI::dbDisconnect(con, shutdown = TRUE)
+
   if (isTRUE(verbose)) {
     message(sprintf(
       "Evidence filter summary: filtered=%d | genomes with AMR=%s | genomes with genome_data=%s",
@@ -1836,12 +2174,12 @@ prepareGenomes <- function(user_bacs,
 
   if (isTRUE(verbose)) message("Step 2: Downloading genomes from BV-BRC (", method, ")")
   ids <- retrieveGenomes(
-    base_dir      = base_dir,
-    user_bacs     = user_bacs,
-    method        = method,
+    base_dir = base_dir,
+    user_bacs = user_bacs,
+    method = method,
     skip_existing = !overwrite,
     evidence_mode = evidence_mode,
-    verbose       = verbose
+    verbose = verbose
   )
   if (length(ids) == 0L) {
     message("No genomes downloaded.")
@@ -1850,9 +2188,9 @@ prepareGenomes <- function(user_bacs,
 
   if (isTRUE(verbose)) message("Step 3: Formatting data into a database for further processing")
   out <- genomeList(
-    base_dir  = base_dir,
+    base_dir = base_dir,
     user_bacs = user_bacs,
-    verbose   = verbose
+    verbose = verbose
   )
 
   if (isTRUE(verbose)) {

From 1c4e21618bdb00e194814cfbc30a1488abb4c944 Mon Sep 17 00:00:00 2001
From: epbrenner <epbrenner@users.noreply.github.com>
Date: Wed, 22 Jul 2026 17:21:08 +0000
Subject: [PATCH 02/14] Style code (GHA)

---
 R/data_curation.R | 42 +++++++++++++++++++++++++++++-------------
 1 file changed, 29 insertions(+), 13 deletions(-)

diff --git a/R/data_curation.R b/R/data_curation.R
index acd6bd1..039d2ca 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -250,7 +250,9 @@
 #' @return A tibble with minimal name repair.
 #' @keywords internal
 .parse_bvbrc_tsv <- function(x) {
-  if (!length(x)) return(tibble::tibble())
+  if (!length(x)) {
+    return(tibble::tibble())
+  }
 
   x <- as.character(x)
   x <- x[nzchar(trimws(x))]
@@ -258,7 +260,9 @@
   # Keep only lines that look like TSV records/headers
   x <- x[grepl("\t", x, fixed = TRUE)]
 
-  if (!length(x)) return(tibble::tibble())
+  if (!length(x)) {
+    return(tibble::tibble())
+  }
 
   txt <- paste(x, collapse = "\n")
   df <- utils::read.table(
@@ -276,7 +280,6 @@
 }
 
 
-
 #' Apply QC filters to BV-BRC genome metadata
 #'
 #' Quality control step with default values for contamination and completeness.
@@ -313,13 +316,17 @@
 
   safe_median <- function(x) {
     x <- x[is.finite(x)]
-    if (!length(x)) return(NA_real_)
+    if (!length(x)) {
+      return(NA_real_)
+    }
     stats::median(x)
   }
 
   safe_sd <- function(x) {
     x <- x[is.finite(x)]
-    if (length(x) < 2L) return(NA_real_)
+    if (length(x) < 2L) {
+      return(NA_real_)
+    }
     stats::sd(x)
   }
 
@@ -365,19 +372,27 @@
       flag_checkm = .data$flag_checkm_missing | .data$flag_checkm_contam | .data$flag_checkm_complete,
       flag_length = if (!is.null(length_deviations) && is.finite(sd_len) && sd_len > 0) {
         !is.na(z_len) & z_len > length_deviations
-      } else FALSE,
+      } else {
+        FALSE
+      },
       flag_gc = if (!is.null(gc_deviations) && is.finite(sd_gc) && sd_gc > 0) {
         !is.na(z_gc) & z_gc > gc_deviations
-      } else FALSE,
+      } else {
+        FALSE
+      },
       flag_cds = if (!is.null(cds_deviations) && is.finite(sd_cds) && sd_cds > 0) {
         !is.na(z_cds) & z_cds > cds_deviations
-      } else FALSE,
+      } else {
+        FALSE
+      },
       qc_keep = !(.data$flag_checkm | .data$flag_length | .data$flag_gc | .data$flag_cds)
     )
 
   make_rej <- function(flag, rule, observed, threshold, comparator) {
     idx <- which(flag %in% TRUE)
-    if (!length(idx)) return(NULL)
+    if (!length(idx)) {
+      return(NULL)
+    }
     tibble::tibble(
       genome_id = qc$`genome.genome_id`[idx],
       genome_name = qc$`genome.genome_name`[idx],
@@ -1353,7 +1368,8 @@ retrieveMetadata <- function(user_bacs,
   dir.create(logs_dir, recursive = TRUE, showWarnings = FALSE)
 
   cat(sprintf("[%s] Writing metadata DuckDB: %s\n", Sys.time(), db_path),
-      file = file.path(logs_dir, "bvbrc.log"), append = TRUE)
+    file = file.path(logs_dir, "bvbrc.log"), append = TRUE
+  )
 
   con <- DBI::dbConnect(duckdb::duckdb(), dbdir = db_path)
   on.exit(try(DBI::dbDisconnect(con, shutdown = TRUE), silent = TRUE), add = TRUE)
@@ -1384,7 +1400,7 @@ retrieveMetadata <- function(user_bacs,
   )
 
   if (isTRUE(debug)) {
-    DBI::dbExecute(con, 'CREATE OR REPLACE TABLE metadata AS SELECT * FROM metadata_full')
+    DBI::dbExecute(con, "CREATE OR REPLACE TABLE metadata AS SELECT * FROM metadata_full")
   } else {
     keep_cols <- setdiff(DBI::dbListFields(con, "metadata_full"), qc_drop_cols)
     select_sql <- paste(DBI::dbQuoteIdentifier(con, keep_cols), collapse = ", ")
@@ -1679,14 +1695,14 @@ retrieveMetadata <- function(user_bacs,
 ### Can remove once testing is completed
 #' Helps check if a complete set exists after DL (.fna + .PATRIC.faa + .PATRIC.gff)
 #' @keywords internal
-#.is_complete_set <- function(dir, genomeID, min_bytes = 100) {
+# .is_complete_set <- function(dir, genomeID, min_bytes = 100) {
 #  fna <- file.path(dir, paste0(genomeID, ".fna"))
 #  faa <- file.path(dir, paste0(genomeID, ".PATRIC.faa"))
 #  gff <- file.path(dir, paste0(genomeID, ".PATRIC.gff"))
 #  paths <- c(fna, faa, gff)
 #  all(file.exists(paths)) &&
 #    all(purrr::map_dbl(paths, function(x) file.info(x)$size) > min_bytes)
-#}
+# }
 
 #' Helps collate completed genomes into a set
 #' @keywords internal

From 5761630df0b13cf480cd5c5545627e225d9e19db Mon Sep 17 00:00:00 2001
From: Evan Pierce Brenner <108823789+epbrenner@users.noreply.github.com>
Date: Wed, 22 Jul 2026 13:31:32 -0600
Subject: [PATCH 03/14] Fixing helper function bug

Uncommenting a helper that was still needed for several calls. Fixing stuff I broke myself!
---
 R/data_curation.R | 19 +++++++++----------
 1 file changed, 9 insertions(+), 10 deletions(-)

diff --git a/R/data_curation.R b/R/data_curation.R
index 039d2ca..c270ab5 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -1691,18 +1691,17 @@ retrieveMetadata <- function(user_bacs,
   list(duckdbConnection = con, table_name = "filtered")
 }
 
-### This has been deprecated after being moved into the loop for .ftpes_download_one()
-### Can remove once testing is completed
+### This has been duplicated after being moved into the loop for .ftpes_download_one()
 #' Helps check if a complete set exists after DL (.fna + .PATRIC.faa + .PATRIC.gff)
 #' @keywords internal
-# .is_complete_set <- function(dir, genomeID, min_bytes = 100) {
-#  fna <- file.path(dir, paste0(genomeID, ".fna"))
-#  faa <- file.path(dir, paste0(genomeID, ".PATRIC.faa"))
-#  gff <- file.path(dir, paste0(genomeID, ".PATRIC.gff"))
-#  paths <- c(fna, faa, gff)
-#  all(file.exists(paths)) &&
-#    all(purrr::map_dbl(paths, function(x) file.info(x)$size) > min_bytes)
-# }
+ .is_complete_set <- function(dir, genomeID, min_bytes = 100) {
+  fna <- file.path(dir, paste0(genomeID, ".fna"))
+  faa <- file.path(dir, paste0(genomeID, ".PATRIC.faa"))
+  gff <- file.path(dir, paste0(genomeID, ".PATRIC.gff"))
+  paths <- c(fna, faa, gff)
+  all(file.exists(paths)) &&
+    all(purrr::map_dbl(paths, function(x) file.info(x)$size) > min_bytes)
+ }
 
 #' Helps collate completed genomes into a set
 #' @keywords internal

From f00c12aec7a9cb7900dbee38c28220f3bfe3feb9 Mon Sep 17 00:00:00 2001
From: epbrenner <epbrenner@users.noreply.github.com>
Date: Wed, 22 Jul 2026 19:33:00 +0000
Subject: [PATCH 04/14] Style code (GHA)

---
 R/data_curation.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/data_curation.R b/R/data_curation.R
index c270ab5..3bb1aec 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -1694,14 +1694,14 @@ retrieveMetadata <- function(user_bacs,
 ### This has been duplicated after being moved into the loop for .ftpes_download_one()
 #' Helps check if a complete set exists after DL (.fna + .PATRIC.faa + .PATRIC.gff)
 #' @keywords internal
- .is_complete_set <- function(dir, genomeID, min_bytes = 100) {
+.is_complete_set <- function(dir, genomeID, min_bytes = 100) {
   fna <- file.path(dir, paste0(genomeID, ".fna"))
   faa <- file.path(dir, paste0(genomeID, ".PATRIC.faa"))
   gff <- file.path(dir, paste0(genomeID, ".PATRIC.gff"))
   paths <- c(fna, faa, gff)
   all(file.exists(paths)) &&
     all(purrr::map_dbl(paths, function(x) file.info(x)$size) > min_bytes)
- }
+}
 
 #' Helps collate completed genomes into a set
 #' @keywords internal

From b0de152b2721195debf8238e7fb9282204812ef5 Mon Sep 17 00:00:00 2001
From: Evan Pierce Brenner <108823789+epbrenner@users.noreply.github.com>
Date: Wed, 22 Jul 2026 17:46:24 -0600
Subject: [PATCH 05/14] Swapping BiocParallel back to Future

Updating to remove BiocParallel implementation that caused bugs. Patching back to prior Future implementation instead. CLI version of genome downloading appears to have issues that may be separate and requires further troubleshooting.
---
 R/data_curation.R   | 114 +++++++++++++++++++++++++++-----------------
 R/data_processing.R | 100 ++++++++++++++++++++++++--------------
 2 files changed, 133 insertions(+), 81 deletions(-)

diff --git a/R/data_curation.R b/R/data_curation.R
index 3bb1aec..f37537c 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -1281,24 +1281,33 @@ retrieveMetadata <- function(user_bacs,
     collapse = ","
   )
 
+  # In retrieveMetadata() replace the BiocParallel block with this:
+
   batch_size <- 500L
   genome_batches <- split(genome_ids, ceiling(seq_along(genome_ids) / batch_size))
 
   n_cores <- max(1L, parallel::detectCores(logical = TRUE) - 1L)
-  param <- BiocParallel::SnowParam(workers = n_cores)
+
+  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 <- BiocParallel::bplapply(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)
-  }, BPPARAM = param)
+  batch_drug_data <- future.apply::future_lapply(
+    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)
+    },
+    future.seed = TRUE
+  )
 
   combined_drug_data_tbl <- dplyr::bind_rows(batch_drug_data) |>
     dplyr::mutate(dplyr::across(dplyr::everything(), ~ iconv(.x, from = "", to = "UTF-8", sub = "")))
@@ -1309,18 +1318,22 @@ retrieveMetadata <- function(user_bacs,
   }
 
   if (isTRUE(verbose)) message("Retrieving genome metadata in batches.")
-  batch_genome_data <- BiocParallel::bplapply(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)
-  }, BPPARAM = param)
+  batch_genome_data <- future.apply::future_lapply(
+    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)
+    },
+    future.seed = TRUE
+  )
 
   combined_genome_data_tbl <- dplyr::bind_rows(batch_genome_data) |>
     dplyr::mutate(dplyr::across(dplyr::everything(), ~ iconv(.x, from = "", to = "UTF-8", sub = "")))
@@ -1951,12 +1964,21 @@ retrieveGenomes <- function(base_dir = ".",
     return(all_complete)
   }
 
+  # In retrieveGenomes() replace the BiocParallel block with this:
+
   if (identical(method, "ftp")) {
     if (isTRUE(verbose)) message("Trying FTPS download. Workers=", ftp_workers)
-    ftp_param <- BiocParallel::SnowParam(workers = max(1L, ftp_workers))
-    ft_ok <- BiocParallel::bplapply(ids, function(gid) .ftpes_download_one(gid, genome_path),
-      BPPARAM = ftp_param
+
+    old_plan <- future::plan()
+    on.exit(future::plan(old_plan), add = TRUE)
+    future::plan(future::multisession, workers = max(1L, ftp_workers))
+
+    ft_ok <- future.apply::future_lapply(
+      ids,
+      function(gid) .ftpes_download_one(gid, genome_path),
+      future.seed = TRUE
     )
+
     ok_ids <- ids[unlist(ft_ok)]
     if (isTRUE(verbose)) message("Complete file sets for ", length(ok_ids), " genomes (FTP).")
     return(c(ok_ids, .list_complete(genome_path, setdiff(ids, ok_ids))))
@@ -1965,33 +1987,35 @@ retrieveGenomes <- function(base_dir = ".",
   # CLI for FASTA, FAA, and GTO, then GFF from GTO
   chunks <- split(ids, ceiling(seq_along(ids) / chunk_size))
 
-  # Parallel chunk containers
   if (isTRUE(verbose)) {
-    message(
-      "CLI being run in parallel for ", length(chunks),
-      " data chunks."
-    )
+    message("CLI being run in parallel for ", length(chunks), " data chunks.")
   }
-  fasta_param <- BiocParallel::SnowParam(workers = max(1L, cli_fasta_workers))
-  fa_res <- BiocParallel::bpmapply(
+
+  old_plan <- future::plan()
+  on.exit(future::plan(old_plan), add = TRUE)
+  future::plan(future::multisession, workers = max(1L, cli_fasta_workers))
+
+  fa_res <- future.apply::future_mapply(
     FUN = function(vec, tag) .cli_dump_fastas_gto_chunk(image, genome_path, vec, tag),
-    vec = chunks, tag = paste0("fa", seq_along(chunks)),
-    SIMPLIFY = TRUE, BPPARAM = fasta_param
+    vec = chunks,
+    tag = paste0("fa", seq_along(chunks)),
+    SIMPLIFY = TRUE,
+    future.seed = TRUE
   )
   if (!all(fa_res) && isTRUE(verbose)) warning(sum(!fa_res), " data chunks failed.")
 
-  # GFF extraction in parallel containers
   if (isTRUE(verbose)) {
-    message(
-      "GFF extraction being run in parallel for ",
-      length(chunks), " data chunks."
-    )
+    message("GFF extraction being run in parallel for ", length(chunks), " data chunks.")
   }
-  gff_param <- BiocParallel::SnowParam(workers = max(1L, cli_gff_workers))
-  g_res <- BiocParallel::bpmapply(
+
+  future::plan(future::multisession, workers = max(1L, cli_gff_workers))
+
+  g_res <- future.apply::future_mapply(
     FUN = function(vec, tag) .cli_export_gff_chunk(image, genome_path, vec, tag),
-    vec = chunks, tag = paste0("gff", seq_along(chunks)),
-    SIMPLIFY = TRUE, BPPARAM = gff_param
+    vec = chunks,
+    tag = paste0("gff", seq_along(chunks)),
+    SIMPLIFY = TRUE,
+    future.seed = TRUE
   )
   if (!all(g_res) && isTRUE(verbose)) warning(sum(!g_res), " GFF chunks had failures.")
 
diff --git a/R/data_processing.R b/R/data_processing.R
index 17b7682..4ff8a61 100644
--- a/R/data_processing.R
+++ b/R/data_processing.R
@@ -27,6 +27,30 @@ NULL
   sub(pattern, container_root, x_unix)
 }
 
+#' Set a future plan for parallel processing
+#'
+#' Sets a `future` plan for the duration of a block and restores the
+#' previous plan on exit.
+#'
+#' @param workers Integer. Number of workers to use; if <= 1, uses sequential mode.
+#' @param plan Character. Either "multisession" or "sequential".
+#'
+#' @return Invisibly returns TRUE.
+#' @keywords internal
+.with_future_plan <- function(workers, plan = c("multisession", "sequential")) {
+  plan <- match.arg(plan)
+  old_plan <- future::plan()
+  on.exit(future::plan(old_plan), add = TRUE)
+
+  if (is.null(workers) || workers <= 1L || identical(plan, "sequential")) {
+    future::plan(future::sequential)
+  } else {
+    future::plan(future::multisession, workers = workers)
+  }
+
+  invisible(TRUE)
+}
+
 # Launch Panaroo to build a pangenome (per batch)
 #' processPanaroo()
 #'
@@ -203,21 +227,20 @@ NULL
   # Ensure sum of per-job CPUs does not exceed `threads`
   panaroo_threads_per_job <- max(1L, floor(threads / n_jobs))
 
-  param <- BiocParallel::SnowParam(workers = max(1L, n_jobs))
-  batch_panaroo_run <- BiocParallel::bplapply(
+  .with_future_plan(workers = n_jobs)
+
+  batch_panaroo_run <- furrr::future_map(
     panaroo_batches,
-    function(batch) {
-      .processPanaroo(
-        batch_input             = batch,
-        output_path             = output_path,
-        core_threshold          = core_threshold,
-        len_dif_percent         = len_dif_percent,
-        cluster_threshold       = cluster_threshold,
-        family_seq_identity     = family_seq_identity,
-        panaroo_threads_per_job = panaroo_threads_per_job
-      )
-    },
-    BPPARAM = param
+    ~ .processPanaroo(
+      batch_input             = .x,
+      output_path             = output_path,
+      core_threshold          = core_threshold,
+      len_dif_percent         = len_dif_percent,
+      cluster_threshold       = cluster_threshold,
+      family_seq_identity     = family_seq_identity,
+      panaroo_threads_per_job = panaroo_threads_per_job
+    ),
+    .options = furrr::furrr_options(seed = TRUE)
   )
 
   invisible(batch_panaroo_run)
@@ -1205,7 +1228,6 @@ domainFromIPR <- function(duckdb_path,
   chunks <- list(sequences_df) # Force 1 chunk for RAM limits
 
   # Forcing 1 container operation for RAM limits
-  workers <- 1
   cpu_per_container <- threads
 
   message(sprintf(
@@ -1213,27 +1235,33 @@ domainFromIPR <- function(duckdb_path,
     cpu_per_container
   ))
 
-  results <- BiocParallel::bplapply(seq_along(chunks), function(i) {
-    res <- try(
-      .process_chunk(
-        chunk         = chunks[[i]],
-        path          = path,
-        ipr_data_path = ipr_data_path,
-        out_file_base = out_file_base,
-        appl          = appl,
-        chunk_id      = i,
-        threads       = cpu_per_container,
-        file_format   = file_format,
-        docker_image  = ipr_image
-      ),
-      silent = TRUE
-    )
-    if (inherits(res, "try-error")) {
-      message(sprintf("Chunk %d failed: %s", i, as.character(res)))
-      return(NULL)
-    }
-    res
-  }, BPPARAM = BiocParallel::SerialParam())
+  .with_future_plan(workers = 1)
+
+  results <- future.apply::future_lapply(
+    seq_along(chunks),
+    function(i) {
+      res <- try(
+        .process_chunk(
+          chunk         = chunks[[i]],
+          path          = path,
+          ipr_data_path = ipr_data_path,
+          out_file_base = out_file_base,
+          appl          = appl,
+          chunk_id      = i,
+          threads       = cpu_per_container,
+          file_format   = file_format,
+          docker_image  = ipr_image
+        ),
+        silent = TRUE
+      )
+      if (inherits(res, "try-error")) {
+        message(sprintf("Chunk %d failed: %s", i, as.character(res)))
+        return(NULL)
+      }
+      res
+    },
+    future.seed = TRUE
+  )
 
   # Combine results
   tsvs <- Filter(function(x) !is.null(x) && file.exists(x), results)

From eec360f7d4ca5377845550b678a5ae02d3a33a6b Mon Sep 17 00:00:00 2001
From: Evan Pierce Brenner <108823789+epbrenner@users.noreply.github.com>
Date: Thu, 23 Jul 2026 16:42:04 -0600
Subject: [PATCH 06/14] Patching ref_file_path parameter check

data_processing.R runDataProcessing() missing(ref_file_path) check ignores if the default path is used. Updated to accept the default path if a user doesn't specify rather than erroring at the end of the processing steps.
---
 R/data_processing.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/data_processing.R b/R/data_processing.R
index 4ff8a61..6e1fc3e 100644
--- a/R/data_processing.R
+++ b/R/data_processing.R
@@ -1799,7 +1799,7 @@ runDataProcessing <- function(duckdb_path,
   )
 
   # 4) Clean metadata and export Parquet + Parquet-backed DuckDB
-  if (missing(ref_file_path) || is.null(ref_file_path)) {
+  if (is.null(ref_file_path) || !nzchar(ref_file_path)) {
     stop("`ref_file_path` (directory with reference TSVs) must be provided to cleanData().")
   }
   if (isTRUE(verbose)) message("Cleaning metadata and exporting Parquet-backed views.")

From ee3e0e430c17f96973bae9e3bd5d288f00a8b6d4 Mon Sep 17 00:00:00 2001
From: Evan Pierce Brenner <108823789+epbrenner@users.noreply.github.com>
Date: Thu, 23 Jul 2026 18:17:03 -0600
Subject: [PATCH 07/14] Genome download updates in data_curation.R

Updated default workers to 8 for CLI and FTP downloads, made two-try FTP download option the default, added a helper to purge broken file sets (i.e., if .faa is missing for an AccID but .fna and .gff download successfully).
---
 R/data_curation.R | 151 ++++++++++++++++++++++++++++++++--------------
 1 file changed, 107 insertions(+), 44 deletions(-)

diff --git a/R/data_curation.R b/R/data_curation.R
index f37537c..b8cf906 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -98,8 +98,8 @@
 #' hiccup. If 2nd pass fails, log and give up on that file.
 #' @keywords internal
 .ftpes_download_two_pass <- function(genome_ids, out_dir,
-                                     workers_first = 4L,
-                                     workers_second = 4L,
+                                     workers_first = 8L,
+                                     workers_second = 8L,
                                      log_file = NULL) {
   genome_ids <- unique(as.character(genome_ids))
   if (!length(genome_ids)) {
@@ -109,17 +109,20 @@
   if (!is.null(log_file)) {
     dir.create(dirname(log_file), recursive = TRUE, showWarnings = FALSE)
     cat(sprintf("[%s] FTPS run start: %d genomes\n", Sys.time(), length(genome_ids)),
-      file = log_file, append = TRUE
+        file = log_file, append = TRUE
     )
   }
 
-  # Pass 1: 30s per-file cap
+  old_plan <- future::plan()
+  on.exit(future::plan(old_plan), add = TRUE)
+
   message("FTPS pass 1 (45s timeout)")
-  future::plan(future::multisession, workers = max(1, workers_first))
+  future::plan(future::multisession, workers = max(1L, workers_first))
   res1 <- future.apply::future_lapply(
     genome_ids,
     function(gid) {
-      ok <- .ftpes_download_one(gid, out_dir,
+      ok <- .ftpes_download_one(
+        gid, out_dir,
         connect_timeout = 10L, max_time = 45L,
         speed_time = 30L, speed_limit = 2048L
       )
@@ -127,32 +130,34 @@
     },
     future.seed = TRUE
   )
+
   ok1 <- vapply(res1, `[[`, logical(1), "ok")
   ok_ids_1 <- genome_ids[ok1]
   fail_ids <- genome_ids[!ok1]
+
   message(sprintf("Pass 1: ok=%d, fail=%d", length(ok_ids_1), length(fail_ids)))
   if (!is.null(log_file)) {
     cat(sprintf("[%s] Pass1 ok=%d fail=%d\n", Sys.time(), length(ok_ids_1), length(fail_ids)),
-      file = log_file, append = TRUE
+        file = log_file, append = TRUE
     )
   }
 
   if (!length(fail_ids)) {
     if (!is.null(log_file)) {
       cat(sprintf("[%s] FTPS run end: all OK\n", Sys.time()),
-        file = log_file, append = TRUE
+          file = log_file, append = TRUE
       )
     }
     return(ok_ids_1)
   }
 
-  # Pass 2: 60s per-file cap where we retry any failures
   message("FTPS pass 2 (120s timeout) for failed genomes")
-  future::plan(future::multisession, workers = max(1, workers_second))
+  future::plan(future::multisession, workers = max(1L, workers_second))
   res2 <- future.apply::future_lapply(
     fail_ids,
     function(gid) {
-      ok <- .ftpes_download_one(gid, out_dir,
+      ok <- .ftpes_download_one(
+        gid, out_dir,
         connect_timeout = 10L, max_time = 120L,
         speed_time = 30L, speed_limit = 2048L
       )
@@ -160,17 +165,19 @@
     },
     future.seed = TRUE
   )
+
   ok2 <- vapply(res2, `[[`, logical(1), "ok")
   ok_ids_2 <- fail_ids[ok2]
   still_fail <- setdiff(fail_ids, ok_ids_2)
+
   message(sprintf("Pass 2: ok=%d, still_fail=%d", length(ok_ids_2), length(still_fail)))
   if (!is.null(log_file)) {
     cat(sprintf("[%s] Pass2 ok=%d still_fail=%d\n", Sys.time(), length(ok_ids_2), length(still_fail)),
-      file = log_file, append = TRUE
+        file = log_file, append = TRUE
     )
     if (length(still_fail)) {
       cat("Fail IDs (excluded): ", paste(head(still_fail, 50), collapse = ", "), "\n",
-        file = log_file, append = TRUE
+          file = log_file, append = TRUE
       )
     }
   }
@@ -480,6 +487,44 @@
   )
 }
 
+# Get rid of genomes that are missing files
+.purge_genome_files <- function(genome_path, genome_ids, log_file = NULL) {
+  genome_ids <- unique(as.character(genome_ids))
+  if (!length(genome_ids)) return(invisible(0L))
+
+  exts <- c(
+    ".fna",
+    ".PATRIC.faa",
+    ".PATRIC.gff",
+    ".gto",
+    ".orig2id.tsv",
+    ".fna.tmp",
+    ".PATRIC.faa.tmp",
+    ".PATRIC.gff.tmp"
+  )
+
+  n_removed <- 0L
+  for (gid in genome_ids) {
+    paths <- file.path(genome_path, paste0(gid, exts))
+    existing <- paths[file.exists(paths)]
+    if (length(existing)) {
+      unlink(existing, force = TRUE)
+      n_removed <- n_removed + length(existing)
+    }
+  }
+
+  if (!is.null(log_file) && length(genome_ids)) {
+    cat(
+      sprintf("[%s] Purged incomplete genomes: %s\n",
+              Sys.time(), paste(head(genome_ids, 50), collapse = ", ")
+      ),
+      file = log_file, append = TRUE
+    )
+  }
+
+  invisible(n_removed)
+}
+
 # Make sure the BV-BRC metadata live where they're supposed to
 .ensure_bvbrc_cache <- function(base_dir = ".",
                                 verbose = TRUE,
@@ -1281,8 +1326,6 @@ retrieveMetadata <- function(user_bacs,
     collapse = ","
   )
 
-  # In retrieveMetadata() replace the BiocParallel block with this:
-
   batch_size <- 500L
   genome_batches <- split(genome_ids, ceiling(seq_along(genome_ids) / batch_size))
 
@@ -1704,7 +1747,6 @@ retrieveMetadata <- function(user_bacs,
   list(duckdbConnection = con, table_name = "filtered")
 }
 
-### This has been duplicated after being moved into the loop for .ftpes_download_one()
 #' Helps check if a complete set exists after DL (.fna + .PATRIC.faa + .PATRIC.gff)
 #' @keywords internal
 .is_complete_set <- function(dir, genomeID, min_bytes = 100) {
@@ -1894,8 +1936,8 @@ retrieveMetadata <- function(user_bacs,
 #' @param image Docker image for CLI path (default "danylmb/bvbrc:5.3").
 #' @param skip_existing Logical; if TRUE, do not re-download genomes already complete. Default TRUE.
 #' @param ftp_workers Parallel workers for FTP path (default 8).
-#' @param cli_fasta_workers Parallel chunk containers for FASTA+GTO (default 4).
-#' @param cli_gff_workers Parallel chunk containers for GFF export (default 4).
+#' @param cli_fasta_workers Parallel chunk containers for FASTA+GTO (default 8).
+#' @param cli_gff_workers Parallel chunk containers for GFF export (default 8).
 #' @param chunk_size Genomes per chunk container (default 50).
 #' @param verbose Verbose messages.
 #' @return Character vector of genome IDs with complete file sets on disk.
@@ -1905,9 +1947,9 @@ retrieveGenomes <- function(base_dir = ".",
                             method = c("ftp", "cli"),
                             image = "danylmb/bvbrc:5.3",
                             skip_existing = TRUE,
-                            ftp_workers = 4L,
-                            cli_fasta_workers = 4L,
-                            cli_gff_workers = 4L,
+                            ftp_workers = 8L,
+                            cli_fasta_workers = 8L,
+                            cli_gff_workers = 8L,
                             chunk_size = 50L,
                             evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"), # NEW
                             verbose = TRUE) {
@@ -1921,11 +1963,15 @@ retrieveGenomes <- function(base_dir = ".",
   db_path <- paths$db_path
   con0 <- DBI::dbConnect(duckdb::duckdb(), dbdir = db_path)
   has_filtered <- "filtered" %in% DBI::dbListTables(con0)
+
   if (has_filtered) {
     if (isTRUE(verbose)) message("Using existing 'filtered' table (skipping re-filter).")
     con <- con0
     tbl <- "filtered"
+    on.exit(try(DBI::dbDisconnect(con0, shutdown = TRUE), silent = TRUE), add = TRUE)
   } else {
+    DBI::dbDisconnect(con0, shutdown = TRUE)
+
     if (isTRUE(verbose)) message("No 'filtered' table found; filtering now.")
     f_out <- .filterGenomes(
       base_dir = base_dir,
@@ -1935,6 +1981,7 @@ retrieveGenomes <- function(base_dir = ".",
     )
     con <- f_out$duckdbConnection
     tbl <- f_out$table_name
+    on.exit(try(DBI::dbDisconnect(con, shutdown = TRUE), silent = TRUE), add = TRUE)
   }
 
   ids <- tibble::as_tibble(DBI::dbReadTable(con, tbl)) |>
@@ -1955,33 +2002,39 @@ retrieveGenomes <- function(base_dir = ".",
 
   if (length(ids) == 0L) {
     if (isTRUE(verbose)) message("All genomes already complete.")
-    all_complete <- .list_complete(
-      genome_path,
-      tibble::as_tibble(DBI::dbReadTable(con, tbl)) |>
-        dplyr::distinct(`genome.genome_id`) |>
-        dplyr::pull(`genome.genome_id`)
-    )
-    return(all_complete)
+    all_ids <- tibble::as_tibble(DBI::dbReadTable(con, tbl)) |>
+      dplyr::distinct(`genome.genome_id`) |>
+      dplyr::pull(`genome.genome_id`)
+    return(.list_complete(genome_path, all_ids))
   }
 
-  # In retrieveGenomes() replace the BiocParallel block with this:
-
   if (identical(method, "ftp")) {
     if (isTRUE(verbose)) message("Trying FTPS download. Workers=", ftp_workers)
 
-    old_plan <- future::plan()
-    on.exit(future::plan(old_plan), add = TRUE)
-    future::plan(future::multisession, workers = max(1L, ftp_workers))
-
-    ft_ok <- future.apply::future_lapply(
-      ids,
-      function(gid) .ftpes_download_one(gid, genome_path),
-      future.seed = TRUE
+    ok_ids <- .ftpes_download_two_pass(
+      genome_ids = ids,
+      out_dir = genome_path,
+      workers_first = ftp_workers,
+      workers_second = ftp_workers,
+      log_file = file.path(logs_dir, "ftp_download.log")
     )
 
-    ok_ids <- ids[unlist(ft_ok)]
-    if (isTRUE(verbose)) message("Complete file sets for ", length(ok_ids), " genomes (FTP).")
-    return(c(ok_ids, .list_complete(genome_path, setdiff(ids, ok_ids))))
+    dropped_ids <- setdiff(ids, ok_ids)
+    if (length(dropped_ids)) {
+      .purge_genome_files(
+        genome_path = genome_path,
+        genome_ids = dropped_ids,
+        log_file = file.path(logs_dir, "ftp_download.log")
+      )
+    }
+
+    if (isTRUE(verbose)) {
+      message("Complete file sets for ", length(ok_ids), " genomes (FTP).")
+      if (length(dropped_ids)) {
+        message(length(dropped_ids), " genomes were excluded because BV-BRC did not provide a complete file set.")
+      }
+    }
+    return(ok_ids)
   }
 
   # CLI for FASTA, FAA, and GTO, then GFF from GTO
@@ -2044,6 +2097,7 @@ retrieveGenomes <- function(base_dir = ".",
 #' @return A list with duckdbConnection and table_name = "files".
 genomeList <- function(base_dir = ".",
                        user_bacs,
+                       expected_ids = NULL,
                        verbose = TRUE) {
   base_dir <- normalizePath(base_dir, mustWork = FALSE)
   paths <- .buildDBpath(base_dir = base_dir, user_bacs = user_bacs)
@@ -2065,6 +2119,11 @@ genomeList <- function(base_dir = ".",
 
   genome_ids <- unique(c(gff_ids, fna_ids, faa_ids))
 
+  # Make sure all the files we want are present
+  if (!is.null(expected_ids)) {
+    genome_ids <- intersect(genome_ids, unique(as.character(expected_ids)))
+  }
+
   list_of_files <- purrr::map(genome_ids, function(genomeID) {
     gff_path <- file.path(genome_path, paste0(genomeID, ".PATRIC.gff"))
     fna_path <- file.path(genome_path, paste0(genomeID, ".fna"))
@@ -2154,8 +2213,8 @@ prepareGenomes <- function(user_bacs,
                            method = c("ftp", "cli"),
                            overwrite = FALSE,
                            evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"),
-                           checkm_contam = 10,
-                           checkm_complete = 90,
+                           checkm_contam = 5,
+                           checkm_complete = 95,
                            gc_deviations = NULL,
                            length_deviations = NULL,
                            cds_deviations = NULL,
@@ -2229,11 +2288,15 @@ prepareGenomes <- function(user_bacs,
   out <- genomeList(
     base_dir = base_dir,
     user_bacs = user_bacs,
+    expected_ids = ids,
     verbose = verbose
   )
 
   if (isTRUE(verbose)) {
-    message("Done. Files are ready! Continue with downstream processing with runDataProcessing().")
+    message("Done. Files are ready!")
+    message("")
+    message("Continue with downstream processing using:")
+    message('runDataProcessing("', normalizePath(paths$db_path), '")')
   }
   out
 }

From 02d0c9c1c060e7b53585049cce9bf37c89c25c30 Mon Sep 17 00:00:00 2001
From: epbrenner <epbrenner@users.noreply.github.com>
Date: Fri, 24 Jul 2026 00:18:28 +0000
Subject: [PATCH 08/14] Style code (GHA)

---
 R/data_curation.R | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/R/data_curation.R b/R/data_curation.R
index b8cf906..5a664c0 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -109,7 +109,7 @@
   if (!is.null(log_file)) {
     dir.create(dirname(log_file), recursive = TRUE, showWarnings = FALSE)
     cat(sprintf("[%s] FTPS run start: %d genomes\n", Sys.time(), length(genome_ids)),
-        file = log_file, append = TRUE
+      file = log_file, append = TRUE
     )
   }
 
@@ -138,14 +138,14 @@
   message(sprintf("Pass 1: ok=%d, fail=%d", length(ok_ids_1), length(fail_ids)))
   if (!is.null(log_file)) {
     cat(sprintf("[%s] Pass1 ok=%d fail=%d\n", Sys.time(), length(ok_ids_1), length(fail_ids)),
-        file = log_file, append = TRUE
+      file = log_file, append = TRUE
     )
   }
 
   if (!length(fail_ids)) {
     if (!is.null(log_file)) {
       cat(sprintf("[%s] FTPS run end: all OK\n", Sys.time()),
-          file = log_file, append = TRUE
+        file = log_file, append = TRUE
       )
     }
     return(ok_ids_1)
@@ -173,11 +173,11 @@
   message(sprintf("Pass 2: ok=%d, still_fail=%d", length(ok_ids_2), length(still_fail)))
   if (!is.null(log_file)) {
     cat(sprintf("[%s] Pass2 ok=%d still_fail=%d\n", Sys.time(), length(ok_ids_2), length(still_fail)),
-        file = log_file, append = TRUE
+      file = log_file, append = TRUE
     )
     if (length(still_fail)) {
       cat("Fail IDs (excluded): ", paste(head(still_fail, 50), collapse = ", "), "\n",
-          file = log_file, append = TRUE
+        file = log_file, append = TRUE
       )
     }
   }
@@ -490,7 +490,9 @@
 # Get rid of genomes that are missing files
 .purge_genome_files <- function(genome_path, genome_ids, log_file = NULL) {
   genome_ids <- unique(as.character(genome_ids))
-  if (!length(genome_ids)) return(invisible(0L))
+  if (!length(genome_ids)) {
+    return(invisible(0L))
+  }
 
   exts <- c(
     ".fna",
@@ -515,8 +517,9 @@
 
   if (!is.null(log_file) && length(genome_ids)) {
     cat(
-      sprintf("[%s] Purged incomplete genomes: %s\n",
-              Sys.time(), paste(head(genome_ids, 50), collapse = ", ")
+      sprintf(
+        "[%s] Purged incomplete genomes: %s\n",
+        Sys.time(), paste(head(genome_ids, 50), collapse = ", ")
       ),
       file = log_file, append = TRUE
     )

From ec6209f8fe77d7ed643a1f00bf0bd6ac4ba2b2a9 Mon Sep 17 00:00:00 2001
From: Janani Ravi <janani.ravi@cuanschutz.edu>
Date: Mon, 27 Jul 2026 14:03:02 -0600
Subject: [PATCH 09/14] rename functions

ftpes --> ftps
is_complete --> .is_complete
---
 R/data_curation.R | 34 +++++++++++++++++-----------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/R/data_curation.R b/R/data_curation.R
index 5a664c0..dedfdfa 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -33,16 +33,16 @@
 #' Helps appropriately interface with BV-BRC FTPS server, and avoids getting stuck
 #' when malformed files can hang an FTPS connection by introducing safeguards
 #' @keywords internal
-.ftpes_download_one <- function(genomeID, out_dir,
-                                connect_timeout = 10L,
-                                max_time = 30L,
-                                speed_time = 30L,
-                                speed_limit = 2048L,
-                                min_bytes = 100L) {
+.ftps_download_one <- function(genomeID, out_dir,
+                               connect_timeout = 10L,
+                               max_time = 30L,
+                               speed_time = 30L,
+                               speed_limit = 2048L,
+                               min_bytes = 100L) {
   dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
 
   # Adding this prior dot function helper into this loop for parallel workers
-  is_complete <- function(dir, genomeID, min_bytes = 100L) {
+  .is_complete <- function(dir, genomeID, min_bytes = 100L) {
     fna <- file.path(dir, paste0(genomeID, ".fna"))
     faa <- file.path(dir, paste0(genomeID, ".PATRIC.faa"))
     gff <- file.path(dir, paste0(genomeID, ".PATRIC.gff"))
@@ -90,17 +90,17 @@
     }
   }
 
-  is_complete(out_dir, genomeID, min_bytes = min_bytes)
+  .is_complete(out_dir, genomeID, min_bytes = min_bytes)
 }
 
-#' Helps manage FTPS downloading from BV-BRC, tryng a quick download first, and
+#' Helps manage FTPS downloading from BV-BRC, trying a quick download first, and
 #' if that fails, trying a longer timeout 2nd pass at the end in case it was a
-#' hiccup. If 2nd pass fails, log and give up on that file.
+#' hiccup. If the 2nd pass fails, log and give up on that file.
 #' @keywords internal
-.ftpes_download_two_pass <- function(genome_ids, out_dir,
-                                     workers_first = 8L,
-                                     workers_second = 8L,
-                                     log_file = NULL) {
+.ftps_download_two_pass <- function(genome_ids, out_dir,
+                                    workers_first = 8L,
+                                    workers_second = 8L,
+                                    log_file = NULL) {
   genome_ids <- unique(as.character(genome_ids))
   if (!length(genome_ids)) {
     return(character(0))
@@ -121,7 +121,7 @@
   res1 <- future.apply::future_lapply(
     genome_ids,
     function(gid) {
-      ok <- .ftpes_download_one(
+      ok <- .ftps_download_one(
         gid, out_dir,
         connect_timeout = 10L, max_time = 45L,
         speed_time = 30L, speed_limit = 2048L
@@ -156,7 +156,7 @@
   res2 <- future.apply::future_lapply(
     fail_ids,
     function(gid) {
-      ok <- .ftpes_download_one(
+      ok <- .ftps_download_one(
         gid, out_dir,
         connect_timeout = 10L, max_time = 120L,
         speed_time = 30L, speed_limit = 2048L
@@ -2014,7 +2014,7 @@ retrieveGenomes <- function(base_dir = ".",
   if (identical(method, "ftp")) {
     if (isTRUE(verbose)) message("Trying FTPS download. Workers=", ftp_workers)
 
-    ok_ids <- .ftpes_download_two_pass(
+    ok_ids <- .ftps_download_two_pass(
       genome_ids = ids,
       out_dir = genome_path,
       workers_first = ftp_workers,

From e2b88bdab52e67fc9410acdd13ca61429f1e4363 Mon Sep 17 00:00:00 2001
From: Emily Boyer <eboyer8@msudenver.edu>
Date: Tue, 28 Jul 2026 09:38:07 -0600
Subject: [PATCH 10/14] Add future/furrr to DESCRIPTION; fix
 .with_future_plan() scoping bug; consolidate on furrr

---
 DESCRIPTION         |  2 ++
 R/data_curation.R   | 38 ++++++++++++++++++--------------------
 R/data_processing.R | 40 ++++++++++++----------------------------
 3 files changed, 32 insertions(+), 48 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 784b8f4..b5fc9a6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -34,6 +34,8 @@ Imports:
     data.table,
     dplyr,
     duckdb,
+    furrr,
+    future,
     glue,
     purrr,
     readr,
diff --git a/R/data_curation.R b/R/data_curation.R
index dedfdfa..afd5df1 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -118,7 +118,7 @@
 
   message("FTPS pass 1 (45s timeout)")
   future::plan(future::multisession, workers = max(1L, workers_first))
-  res1 <- future.apply::future_lapply(
+  res1 <- furrr::future_map(
     genome_ids,
     function(gid) {
       ok <- .ftps_download_one(
@@ -128,7 +128,7 @@
       )
       list(gid = gid, ok = ok)
     },
-    future.seed = TRUE
+    .options = furrr::furrr_options(seed = TRUE)
   )
 
   ok1 <- vapply(res1, `[[`, logical(1), "ok")
@@ -153,7 +153,7 @@
 
   message("FTPS pass 2 (120s timeout) for failed genomes")
   future::plan(future::multisession, workers = max(1L, workers_second))
-  res2 <- future.apply::future_lapply(
+  res2 <- furrr::future_map(
     fail_ids,
     function(gid) {
       ok <- .ftps_download_one(
@@ -163,7 +163,7 @@
       )
       list(gid = gid, ok = ok)
     },
-    future.seed = TRUE
+    .options = furrr::furrr_options(seed = TRUE)
   )
 
   ok2 <- vapply(res2, `[[`, logical(1), "ok")
@@ -1339,7 +1339,7 @@ retrieveMetadata <- function(user_bacs,
   future::plan(future::multisession, workers = n_cores)
 
   if (isTRUE(verbose)) message("Retrieving AMR phenotype data in batches.")
-  batch_drug_data <- future.apply::future_lapply(
+  batch_drug_data <- furrr::future_map(
     genome_batches,
     function(batch) {
       raw <- .extractAMRtable(
@@ -1352,7 +1352,7 @@ retrieveMetadata <- function(user_bacs,
       )
       .parse_bvbrc_tsv(raw)
     },
-    future.seed = TRUE
+    .options = furrr::furrr_options(seed = TRUE)
   )
 
   combined_drug_data_tbl <- dplyr::bind_rows(batch_drug_data) |>
@@ -1364,7 +1364,7 @@ retrieveMetadata <- function(user_bacs,
   }
 
   if (isTRUE(verbose)) message("Retrieving genome metadata in batches.")
-  batch_genome_data <- future.apply::future_lapply(
+  batch_genome_data <- furrr::future_map(
     genome_batches,
     function(batch) {
       raw <- .extractGenomeData(
@@ -1378,7 +1378,7 @@ retrieveMetadata <- function(user_bacs,
       )
       .parse_bvbrc_tsv(raw)
     },
-    future.seed = TRUE
+    .options = furrr::furrr_options(seed = TRUE)
   )
 
   combined_genome_data_tbl <- dplyr::bind_rows(batch_genome_data) |>
@@ -2051,12 +2051,11 @@ retrieveGenomes <- function(base_dir = ".",
   on.exit(future::plan(old_plan), add = TRUE)
   future::plan(future::multisession, workers = max(1L, cli_fasta_workers))
 
-  fa_res <- future.apply::future_mapply(
-    FUN = function(vec, tag) .cli_dump_fastas_gto_chunk(image, genome_path, vec, tag),
-    vec = chunks,
-    tag = paste0("fa", seq_along(chunks)),
-    SIMPLIFY = TRUE,
-    future.seed = TRUE
+  fa_res <- furrr::future_map2_lgl(
+    chunks,
+    paste0("fa", seq_along(chunks)),
+    function(vec, tag) .cli_dump_fastas_gto_chunk(image, genome_path, vec, tag),
+    .options = furrr::furrr_options(seed = TRUE)
   )
   if (!all(fa_res) && isTRUE(verbose)) warning(sum(!fa_res), " data chunks failed.")
 
@@ -2066,12 +2065,11 @@ retrieveGenomes <- function(base_dir = ".",
 
   future::plan(future::multisession, workers = max(1L, cli_gff_workers))
 
-  g_res <- future.apply::future_mapply(
-    FUN = function(vec, tag) .cli_export_gff_chunk(image, genome_path, vec, tag),
-    vec = chunks,
-    tag = paste0("gff", seq_along(chunks)),
-    SIMPLIFY = TRUE,
-    future.seed = TRUE
+  g_res <- furrr::future_map2_lgl(
+    chunks,
+    paste0("gff", seq_along(chunks)),
+    function(vec, tag) .cli_export_gff_chunk(image, genome_path, vec, tag),
+    .options = furrr::furrr_options(seed = TRUE)
   )
   if (!all(g_res) && isTRUE(verbose)) warning(sum(!g_res), " GFF chunks had failures.")
 
diff --git a/R/data_processing.R b/R/data_processing.R
index 6e1fc3e..8b8d518 100644
--- a/R/data_processing.R
+++ b/R/data_processing.R
@@ -27,30 +27,6 @@ NULL
   sub(pattern, container_root, x_unix)
 }
 
-#' Set a future plan for parallel processing
-#'
-#' Sets a `future` plan for the duration of a block and restores the
-#' previous plan on exit.
-#'
-#' @param workers Integer. Number of workers to use; if <= 1, uses sequential mode.
-#' @param plan Character. Either "multisession" or "sequential".
-#'
-#' @return Invisibly returns TRUE.
-#' @keywords internal
-.with_future_plan <- function(workers, plan = c("multisession", "sequential")) {
-  plan <- match.arg(plan)
-  old_plan <- future::plan()
-  on.exit(future::plan(old_plan), add = TRUE)
-
-  if (is.null(workers) || workers <= 1L || identical(plan, "sequential")) {
-    future::plan(future::sequential)
-  } else {
-    future::plan(future::multisession, workers = workers)
-  }
-
-  invisible(TRUE)
-}
-
 # Launch Panaroo to build a pangenome (per batch)
 #' processPanaroo()
 #'
@@ -227,7 +203,13 @@ NULL
   # Ensure sum of per-job CPUs does not exceed `threads`
   panaroo_threads_per_job <- max(1L, floor(threads / n_jobs))
 
-  .with_future_plan(workers = n_jobs)
+  old_plan <- future::plan()
+  on.exit(future::plan(old_plan), add = TRUE)
+  if (n_jobs <= 1L) {
+    future::plan(future::sequential)
+  } else {
+    future::plan(future::multisession, workers = n_jobs)
+  }
 
   batch_panaroo_run <- furrr::future_map(
     panaroo_batches,
@@ -1235,9 +1217,11 @@ domainFromIPR <- function(duckdb_path,
     cpu_per_container
   ))
 
-  .with_future_plan(workers = 1)
+  old_plan <- future::plan()
+  on.exit(future::plan(old_plan), add = TRUE)
+  future::plan(future::sequential)
 
-  results <- future.apply::future_lapply(
+  results <- furrr::future_map(
     seq_along(chunks),
     function(i) {
       res <- try(
@@ -1260,7 +1244,7 @@ domainFromIPR <- function(duckdb_path,
       }
       res
     },
-    future.seed = TRUE
+    .options = furrr::furrr_options(seed = TRUE)
   )
 
   # Combine results

From d6582240df64f0d313a64674602fc33a37aa739a Mon Sep 17 00:00:00 2001
From: Evan Pierce Brenner <108823789+epbrenner@users.noreply.github.com>
Date: Tue, 28 Jul 2026 09:45:07 -0600
Subject: [PATCH 11/14] Applying fixes from code review

Cleaning up parallel implementation in data_curation.R, standardizing with furrr and dropping old duplicated within-loop helper function that was needed for earlier BiocParallel implementation. Standardized values for contamination and completeness across functions. Added a parameter in prepareGenomes to allow users to set number of workers. Reduced console clutter for final "out" print.
---
 R/data_curation.R | 162 +++++++++++++++++++++++-----------------------
 1 file changed, 81 insertions(+), 81 deletions(-)

diff --git a/R/data_curation.R b/R/data_curation.R
index dedfdfa..33a3943 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -41,16 +41,6 @@
                                min_bytes = 100L) {
   dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
 
-  # Adding this prior dot function helper into this loop for parallel workers
-  .is_complete <- function(dir, genomeID, min_bytes = 100L) {
-    fna <- file.path(dir, paste0(genomeID, ".fna"))
-    faa <- file.path(dir, paste0(genomeID, ".PATRIC.faa"))
-    gff <- file.path(dir, paste0(genomeID, ".PATRIC.gff"))
-    paths <- c(fna, faa, gff)
-    all(file.exists(paths)) &&
-      all(vapply(paths, function(x) file.info(x)$size, numeric(1)) > min_bytes)
-  }
-
   exts <- c(".fna", ".PATRIC.faa", ".PATRIC.gff")
   for (ext in exts) {
     dest <- file.path(out_dir, paste0(genomeID, ext))
@@ -90,7 +80,7 @@
     }
   }
 
-  .is_complete(out_dir, genomeID, min_bytes = min_bytes)
+  .is_complete_set(out_dir, genomeID, min_bytes = min_bytes)
 }
 
 #' Helps manage FTPS downloading from BV-BRC, trying a quick download first, and
@@ -118,7 +108,8 @@
 
   message("FTPS pass 1 (45s timeout)")
   future::plan(future::multisession, workers = max(1L, workers_first))
-  res1 <- future.apply::future_lapply(
+
+  res1 <- furrr::future_map(
     genome_ids,
     function(gid) {
       ok <- .ftps_download_one(
@@ -128,32 +119,17 @@
       )
       list(gid = gid, ok = ok)
     },
-    future.seed = TRUE
+    .options = furrr::furrr_options(seed = TRUE)
   )
 
-  ok1 <- vapply(res1, `[[`, logical(1), "ok")
+  ok1 <- purrr::map_lgl(res1, "ok")
   ok_ids_1 <- genome_ids[ok1]
   fail_ids <- genome_ids[!ok1]
 
-  message(sprintf("Pass 1: ok=%d, fail=%d", length(ok_ids_1), length(fail_ids)))
-  if (!is.null(log_file)) {
-    cat(sprintf("[%s] Pass1 ok=%d fail=%d\n", Sys.time(), length(ok_ids_1), length(fail_ids)),
-      file = log_file, append = TRUE
-    )
-  }
-
-  if (!length(fail_ids)) {
-    if (!is.null(log_file)) {
-      cat(sprintf("[%s] FTPS run end: all OK\n", Sys.time()),
-        file = log_file, append = TRUE
-      )
-    }
-    return(ok_ids_1)
-  }
-
   message("FTPS pass 2 (120s timeout) for failed genomes")
   future::plan(future::multisession, workers = max(1L, workers_second))
-  res2 <- future.apply::future_lapply(
+
+  res2 <- furrr::future_map(
     fail_ids,
     function(gid) {
       ok <- .ftps_download_one(
@@ -163,10 +139,10 @@
       )
       list(gid = gid, ok = ok)
     },
-    future.seed = TRUE
+    .options = furrr::furrr_options(seed = TRUE)
   )
 
-  ok2 <- vapply(res2, `[[`, logical(1), "ok")
+  ok2 <- purrr::map_lgl(res2, "ok")
   ok_ids_2 <- fail_ids[ok2]
   still_fail <- setdiff(fail_ids, ok_ids_2)
 
@@ -314,8 +290,8 @@
 #' }
 #' @keywords internal
 .apply_metadata_qc <- function(genome_tbl,
-                               checkm_contam = 10,
-                               checkm_complete = 90,
+                               checkm_contam = 5,
+                               checkm_complete = 95,
                                gc_deviations = NULL,
                                length_deviations = NULL,
                                cds_deviations = NULL) {
@@ -1222,8 +1198,8 @@ retrieveMetadata <- function(user_bacs,
                              abx = "All",
                              overwrite = FALSE,
                              image = "danylmb/bvbrc:5.3",
-                             checkm_contam = 10,
-                             checkm_complete = 90,
+                             checkm_contam = 5,
+                             checkm_complete = 95,
                              gc_deviations = NULL,
                              length_deviations = NULL,
                              cds_deviations = NULL,
@@ -1339,7 +1315,7 @@ retrieveMetadata <- function(user_bacs,
   future::plan(future::multisession, workers = n_cores)
 
   if (isTRUE(verbose)) message("Retrieving AMR phenotype data in batches.")
-  batch_drug_data <- future.apply::future_lapply(
+  batch_drug_data <- furrr::future_map(
     genome_batches,
     function(batch) {
       raw <- .extractAMRtable(
@@ -1352,19 +1328,14 @@ retrieveMetadata <- function(user_bacs,
       )
       .parse_bvbrc_tsv(raw)
     },
-    future.seed = TRUE
+    .options = furrr::furrr_options(seed = TRUE)
   )
 
   combined_drug_data_tbl <- dplyr::bind_rows(batch_drug_data) |>
     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 <- future.apply::future_lapply(
+  batch_genome_data <- furrr::future_map(
     genome_batches,
     function(batch) {
       raw <- .extractGenomeData(
@@ -1378,7 +1349,7 @@ retrieveMetadata <- function(user_bacs,
       )
       .parse_bvbrc_tsv(raw)
     },
-    future.seed = TRUE
+    .options = furrr::furrr_options(seed = TRUE)
   )
 
   combined_genome_data_tbl <- dplyr::bind_rows(batch_genome_data) |>
@@ -1954,28 +1925,32 @@ retrieveGenomes <- function(base_dir = ".",
                             cli_fasta_workers = 8L,
                             cli_gff_workers = 8L,
                             chunk_size = 50L,
-                            evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"), # NEW
+                            evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"),
+                            # NEW
                             verbose = TRUE) {
   method <- match.arg(method)
   evidence_mode <- match.arg(evidence_mode)
   base_dir <- normalizePath(base_dir, mustWork = FALSE)
 
   # Use 'filtered' if already prepared, or start filtering
-  if (isTRUE(verbose)) message("Preparing download set (checking for existing 'filtered').")
+  if (isTRUE(verbose))
+    message("Preparing download set (checking for existing 'filtered').")
   paths <- .buildDBpath(base_dir = base_dir, user_bacs = user_bacs)
   db_path <- paths$db_path
   con0 <- DBI::dbConnect(duckdb::duckdb(), dbdir = db_path)
   has_filtered <- "filtered" %in% DBI::dbListTables(con0)
 
   if (has_filtered) {
-    if (isTRUE(verbose)) message("Using existing 'filtered' table (skipping re-filter).")
+    if (isTRUE(verbose))
+      message("Using existing 'filtered' table (skipping re-filter).")
     con <- con0
     tbl <- "filtered"
     on.exit(try(DBI::dbDisconnect(con0, shutdown = TRUE), silent = TRUE), add = TRUE)
   } else {
     DBI::dbDisconnect(con0, shutdown = TRUE)
 
-    if (isTRUE(verbose)) message("No 'filtered' table found; filtering now.")
+    if (isTRUE(verbose))
+      message("No 'filtered' table found; filtering now.")
     f_out <- .filterGenomes(
       base_dir = base_dir,
       user_bacs = user_bacs,
@@ -1999,12 +1974,14 @@ retrieveGenomes <- function(base_dir = ".",
 
   if (isTRUE(skip_existing)) {
     already <- .list_complete(genome_path, ids)
-    if (isTRUE(verbose)) message(length(already), " genomes already completed; skipping.")
+    if (isTRUE(verbose))
+      message(length(already), " genomes already completed; skipping.")
     ids <- setdiff(ids, already)
   }
 
   if (length(ids) == 0L) {
-    if (isTRUE(verbose)) message("All genomes already complete.")
+    if (isTRUE(verbose))
+      message("All genomes already complete.")
     all_ids <- tibble::as_tibble(DBI::dbReadTable(con, tbl)) |>
       dplyr::distinct(`genome.genome_id`) |>
       dplyr::pull(`genome.genome_id`)
@@ -2012,7 +1989,8 @@ retrieveGenomes <- function(base_dir = ".",
   }
 
   if (identical(method, "ftp")) {
-    if (isTRUE(verbose)) message("Trying FTPS download. Workers=", ftp_workers)
+    if (isTRUE(verbose))
+      message("Trying FTPS download. Workers=", ftp_workers)
 
     ok_ids <- .ftps_download_two_pass(
       genome_ids = ids,
@@ -2034,7 +2012,10 @@ retrieveGenomes <- function(base_dir = ".",
     if (isTRUE(verbose)) {
       message("Complete file sets for ", length(ok_ids), " genomes (FTP).")
       if (length(dropped_ids)) {
-        message(length(dropped_ids), " genomes were excluded because BV-BRC did not provide a complete file set.")
+        message(
+          length(dropped_ids),
+          " genomes were excluded because BV-BRC did not provide a complete file set."
+        )
       }
     }
     return(ok_ids)
@@ -2044,44 +2025,54 @@ retrieveGenomes <- function(base_dir = ".",
   chunks <- split(ids, ceiling(seq_along(ids) / chunk_size))
 
   if (isTRUE(verbose)) {
-    message("CLI being run in parallel for ", length(chunks), " data chunks.")
+    message("CLI being run in parallel for ",
+            length(chunks),
+            " data chunks.")
   }
 
-  old_plan <- future::plan()
-  on.exit(future::plan(old_plan), add = TRUE)
-  future::plan(future::multisession, workers = max(1L, cli_fasta_workers))
-
-  fa_res <- future.apply::future_mapply(
-    FUN = function(vec, tag) .cli_dump_fastas_gto_chunk(image, genome_path, vec, tag),
-    vec = chunks,
-    tag = paste0("fa", seq_along(chunks)),
-    SIMPLIFY = TRUE,
-    future.seed = TRUE
+  run_chunk_phase <- function(vecs, tags, workers, fun) {
+    old_plan <- future::plan()
+    on.exit(future::plan(old_plan), add = TRUE)
+    future::plan(future::multisession, workers = max(1L, workers))
+    furrr::future_map2(vecs, tags, fun, .options = furrr::furrr_options(seed = TRUE))
+  }
+
+  fa_res <- run_chunk_phase(
+    vecs = chunks,
+    tags = paste0("fa", seq_along(chunks)),
+    workers = cli_fasta_workers,
+    fun = function(vec, tag)
+      .cli_dump_fastas_gto_chunk(image, genome_path, vec, tag)
   )
-  if (!all(fa_res) && isTRUE(verbose)) warning(sum(!fa_res), " data chunks failed.")
+  fa_ok <- purrr::map_lgl(fa_res, identity)
+  if (!all(fa_ok) &&
+      isTRUE(verbose))
+    warning(sum(!fa_ok), " data chunks failed.")
 
   if (isTRUE(verbose)) {
-    message("GFF extraction being run in parallel for ", length(chunks), " data chunks.")
+    message("GFF extraction being run in parallel for ",
+            length(chunks),
+            " data chunks.")
   }
 
-  future::plan(future::multisession, workers = max(1L, cli_gff_workers))
-
-  g_res <- future.apply::future_mapply(
-    FUN = function(vec, tag) .cli_export_gff_chunk(image, genome_path, vec, tag),
-    vec = chunks,
-    tag = paste0("gff", seq_along(chunks)),
-    SIMPLIFY = TRUE,
-    future.seed = TRUE
+  g_res <- run_chunk_phase(
+    vecs = chunks,
+    tags = paste0("gff", seq_along(chunks)),
+    workers = cli_gff_workers,
+    fun = function(vec, tag)
+      .cli_export_gff_chunk(image, genome_path, vec, tag)
   )
-  if (!all(g_res) && isTRUE(verbose)) warning(sum(!g_res), " GFF chunks had failures.")
+  g_ok <- purrr::map_lgl(g_res, identity)
+  if (!all(g_ok) &&
+      isTRUE(verbose))
+    warning(sum(!g_ok), " GFF chunks had failures.")
 
   # Success set: .fna + .PATRIC.faa + .PATRIC.gff all present per isolate
   ok_ids <- ids[purrr::map_lgl(ids, .is_complete_set, dir = genome_path)]
   if (isTRUE(verbose)) {
-    message(
-      "Complete file sets downloaded for ",
-      length(ok_ids), " genomes."
-    )
+    message("Complete file sets downloaded for ",
+            length(ok_ids),
+            " genomes.")
   }
   ok_ids
 }
@@ -2197,6 +2188,8 @@ genomeList <- function(base_dir = ".",
 #' @param evidence_mode Character. Sets what types of AMR evidence is acceptable.
 #'    Default `lab_only`. `any` will not require AMR data for downloads. This will
 #'    return very large download lists for many species!
+#' @param num_workers Integer. Parallel workers used for genome download.
+#'    Applied to both FTP and CLI download branches. Default: 8.
 #' @param checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%).
 #' @param checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%).
 #' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content.
@@ -2215,6 +2208,7 @@ prepareGenomes <- function(user_bacs,
                            base_dir = ".",
                            method = c("ftp", "cli"),
                            overwrite = FALSE,
+                           num_workers = 8L,
                            evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"),
                            checkm_contam = 5,
                            checkm_complete = 95,
@@ -2279,6 +2273,9 @@ prepareGenomes <- function(user_bacs,
     user_bacs = user_bacs,
     method = method,
     skip_existing = !overwrite,
+    ftp_workers = num_workers,
+    cli_fasta_workers = num_workers,
+    cli_gff_workers = num_workers,
     evidence_mode = evidence_mode,
     verbose = verbose
   )
@@ -2301,5 +2298,8 @@ prepareGenomes <- function(user_bacs,
     message("Continue with downstream processing using:")
     message('runDataProcessing("', normalizePath(paths$db_path), '")')
   }
-  out
+  invisible(list(
+    duckdb_path = paths$db_path,
+    table_name = "files"
+  ))
 }

From e9b1d39ad810104f562030d73182f7ef28eb59c1 Mon Sep 17 00:00:00 2001
From: Emily Boyer <eboyer8@msudenver.edu>
Date: Tue, 28 Jul 2026 09:51:59 -0600
Subject: [PATCH 12/14] Document future::multisession testing gotcha in README

---
 README.Rmd | 1 +
 README.md  | 8 ++++++++
 2 files changed, 9 insertions(+)

diff --git a/README.Rmd b/README.Rmd
index daa9d62..48c6a58 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -269,6 +269,7 @@ Processing times vary by species and isolate count:
 
 - These numbers will all vary greatly based on isolate number, genome complexity, and available hardware.
 - Parallelization significantly reduces processing time when multiple cores are available.
+- If a `future::multisession` error occurs mid-run (e.g. while testing via `devtools::load_all()` before installing the package), restart your R session fully before retrying. An orphaned background worker process can leave a stale lock on the local DuckDB caches (e.g. `data/bvbrc/bvbrcData.duckdb`), which can produce inconsistent results on the next run that look like a data or QC bug but are actually just leftover session state.
 
 ### Integration with amR suite
 
diff --git a/README.md b/README.md
index 268d529..a2215e3 100644
--- a/README.md
+++ b/README.md
@@ -283,6 +283,14 @@ Processing times vary by species and isolate count:
 - Parallelization significantly reduces processing time when multiple
   cores are available.
 
+- If a `future::multisession` error occurs mid-run (e.g. while testing
+  via `devtools::load_all()` before installing the package), restart
+  your R session fully before retrying. An orphaned background worker
+  process can leave a stale lock on the local DuckDB caches (e.g.
+  `data/bvbrc/bvbrcData.duckdb`), which can produce inconsistent
+  results on the next run that look like a data or QC bug but are
+  actually just leftover session state.
+
 ### Integration with amR suite
 
 amRdata is designed to work seamlessly with other amR packages:

From 2e6af65cf371a5cff05edc8e44a24d68de65f2c7 Mon Sep 17 00:00:00 2001
From: Evan Pierce Brenner <108823789+epbrenner@users.noreply.github.com>
Date: Tue, 28 Jul 2026 13:36:43 -0600
Subject: [PATCH 13/14] Fixed reversions in data_curation.R

Fixed two code reversions that 1) didn't stop pass 2 from being attempted for the FTPS downloader even if everything worked in pass 1, and 2) caught if no drug data was returned in a query and halted the script. Also updated an output continue providing `out`, just invisibly, rather than custom outputs that were by default hidden. These updates were also applied to a separate branch for table exports.

(Thanks Janani!)

Co-Authored-By: Janani Ravi <8397074+jananiravi@users.noreply.github.com>
---
 R/data_curation.R | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/R/data_curation.R b/R/data_curation.R
index 33a3943..20b76cc 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -126,6 +126,17 @@
   ok_ids_1 <- genome_ids[ok1]
   fail_ids <- genome_ids[!ok1]
 
+  message(sprintf("Pass 1: ok=%d, fail=%d", length(ok_ids_1), length(fail_ids)))
+  if (!is.null(log_file)) {
+    cat(sprintf("[%s] Pass1 ok=%d fail=%d\n", Sys.time(), length(ok_ids_1), length(fail_ids)),
+        file = log_file, append = TRUE
+    )
+  }
+
+  if (!length(fail_ids)) {
+    return(ok_ids_1)
+  }
+
   message("FTPS pass 2 (120s timeout) for failed genomes")
   future::plan(future::multisession, workers = max(1L, workers_second))
 
@@ -1334,6 +1345,11 @@ retrieveMetadata <- function(user_bacs,
   combined_drug_data_tbl <- dplyr::bind_rows(batch_drug_data) |>
     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,
@@ -2298,8 +2314,5 @@ prepareGenomes <- function(user_bacs,
     message("Continue with downstream processing using:")
     message('runDataProcessing("', normalizePath(paths$db_path), '")')
   }
-  invisible(list(
-    duckdb_path = paths$db_path,
-    table_name = "files"
-  ))
+  invisible(out)
 }

From fe37c6867cbf1591122c2e0a9de137e8b52b6b23 Mon Sep 17 00:00:00 2001
From: Emily Boyer <eboyer8@msudenver.edu>
Date: Tue, 28 Jul 2026 15:18:29 -0600
Subject: [PATCH 14/14] Add test coverage for .apply_metadata_qc() and
 .parse_bvbrc_tsv()

---
 tests/testthat/test-apply-metadata-qc.R | 64 +++++++++++++++++++++++++
 tests/testthat/test-parse-bvbrc-tsv.R   | 32 +++++++++++++
 2 files changed, 96 insertions(+)
 create mode 100644 tests/testthat/test-apply-metadata-qc.R
 create mode 100644 tests/testthat/test-parse-bvbrc-tsv.R

diff --git a/tests/testthat/test-apply-metadata-qc.R b/tests/testthat/test-apply-metadata-qc.R
new file mode 100644
index 0000000..620a922
--- /dev/null
+++ b/tests/testthat/test-apply-metadata-qc.R
@@ -0,0 +1,64 @@
+.mock_genome_tbl <- function(genome_id, checkm_completeness, checkm_contamination,
+                              genome_name = "x", species = "x",
+                              genome_length = "4600000", gc_content = "50.5", cds = "4300") {
+  tibble::tibble(
+    `genome.genome_id` = genome_id,
+    `genome.genome_name` = genome_name,
+    `genome.species` = species,
+    `genome.genome_length` = genome_length,
+    `genome.gc_content` = gc_content,
+    `genome.cds` = cds,
+    `genome.checkm_completeness` = checkm_completeness,
+    `genome.checkm_contamination` = checkm_contamination
+  )
+}
+
+test_that(".apply_metadata_qc() keeps a genome that passes the default CheckM gate", {
+  tbl <- .mock_genome_tbl("511145.12", checkm_completeness = "99", checkm_contamination = "1")
+
+  out <- .apply_metadata_qc(tbl)
+
+  expect_identical(out$keep_ids, "511145.12")
+  expect_equal(nrow(out$rejections), 0L)
+})
+
+test_that(".apply_metadata_qc() drops a genome with contamination above the threshold", {
+  tbl <- .mock_genome_tbl("511145.12", checkm_completeness = "99", checkm_contamination = "12")
+
+  out <- .apply_metadata_qc(tbl, checkm_contam = 5, checkm_complete = 95)
+
+  expect_length(out$keep_ids, 0L)
+  expect_identical(out$rejections$failed_rule, "checkm_contamination")
+  expect_identical(out$rejections$observed, "12")
+})
+
+test_that(".apply_metadata_qc() drops a genome with completeness below the threshold", {
+  tbl <- .mock_genome_tbl("511145.12", checkm_completeness = "80", checkm_contamination = "1")
+
+  out <- .apply_metadata_qc(tbl, checkm_contam = 5, checkm_complete = 95)
+
+  expect_length(out$keep_ids, 0L)
+  expect_identical(out$rejections$failed_rule, "checkm_completeness")
+})
+
+test_that(".apply_metadata_qc() drops a genome with missing CheckM fields", {
+  tbl <- .mock_genome_tbl("511145.12", checkm_completeness = NA_character_, checkm_contamination = "1")
+
+  out <- .apply_metadata_qc(tbl)
+
+  expect_length(out$keep_ids, 0L)
+  expect_identical(out$rejections$failed_rule, "checkm_missing")
+})
+
+test_that(".apply_metadata_qc() keeps only the passing genomes in a mixed batch", {
+  tbl <- .mock_genome_tbl(
+    genome_id = c("good.1", "bad_contam.1", "bad_complete.1"),
+    checkm_completeness = c("99", "99", "50"),
+    checkm_contamination = c("1", "20", "1")
+  )
+
+  out <- .apply_metadata_qc(tbl)
+
+  expect_identical(out$keep_ids, "good.1")
+  expect_setequal(out$rejections$genome_id, c("bad_contam.1", "bad_complete.1"))
+})
diff --git a/tests/testthat/test-parse-bvbrc-tsv.R b/tests/testthat/test-parse-bvbrc-tsv.R
new file mode 100644
index 0000000..84cb777
--- /dev/null
+++ b/tests/testthat/test-parse-bvbrc-tsv.R
@@ -0,0 +1,32 @@
+test_that(".parse_bvbrc_tsv() returns an empty tibble for empty input", {
+  expect_identical(.parse_bvbrc_tsv(character(0)), tibble::tibble())
+  expect_identical(.parse_bvbrc_tsv(c("", "   ")), tibble::tibble())
+})
+
+test_that(".parse_bvbrc_tsv() parses well-formed TSV lines into a tibble", {
+  lines <- c(
+    "genome_id\tgenome_name",
+    "511145.12\tEscherichia coli",
+    "83333.111\tEscherichia coli K-12"
+  )
+
+  out <- .parse_bvbrc_tsv(lines)
+
+  expect_s3_class(out, "tbl_df")
+  expect_identical(names(out), c("genome_id", "genome_name"))
+  expect_identical(out$genome_id, c("511145.12", "83333.111"))
+})
+
+test_that(".parse_bvbrc_tsv() drops blank lines and non-TSV noise", {
+  lines <- c(
+    "genome_id\tgenome_name",
+    "",
+    "some CLI log line with no tabs",
+    "511145.12\tEscherichia coli"
+  )
+
+  out <- .parse_bvbrc_tsv(lines)
+
+  expect_equal(nrow(out), 1L)
+  expect_identical(out$genome_id, "511145.12")
+})