diff --git a/DESCRIPTION b/DESCRIPTION index 784b8f4..edc18c7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -34,6 +34,8 @@ Imports: data.table, dplyr, duckdb, + furrr, + future, glue, purrr, readr, @@ -46,5 +48,6 @@ VignetteBuilder: knitr Suggests: knitr, rmarkdown, + writexl, testthat (>= 3.0.0) Config/roxygen2/version: 8.0.0 diff --git a/R/data_curation.R b/R/data_curation.R index f056c37..3ce7a8f 100644 --- a/R/data_curation.R +++ b/R/data_curation.R @@ -33,12 +33,12 @@ #' 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, # end if speed_time - speed_limit = 2048L, # B/s - 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) exts <- c(".fna", ".PATRIC.faa", ".PATRIC.gff") @@ -46,7 +46,6 @@ 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 +58,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,32 +68,29 @@ 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) } -#' 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 = 4L, - workers_second = 4L, - 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)) @@ -107,56 +103,60 @@ ) } - # 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)) - res1 <- future.apply::future_lapply( + future::plan(future::multisession, workers = max(1L, workers_first)) + + res1 <- furrr::future_map( genome_ids, function(gid) { - ok <- .ftpes_download_one(gid, out_dir, + ok <- .ftps_download_one( + gid, out_dir, connect_timeout = 10L, max_time = 45L, speed_time = 30L, speed_limit = 2048L ) 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 + 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) } - # 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)) - res2 <- future.apply::future_lapply( + future::plan(future::multisession, workers = max(1L, workers_second)) + + res2 <- furrr::future_map( fail_ids, function(gid) { - ok <- .ftpes_download_one(gid, out_dir, + ok <- .ftps_download_one( + gid, out_dir, connect_timeout = 10L, max_time = 120L, speed_time = 30L, speed_limit = 2048L ) 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) + 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)), @@ -234,36 +234,309 @@ 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 max_checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%). +#' @param min_checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%). +#' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content. +#' @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, + max_checkm_contam = 5, + min_checkm_complete = 95, + 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) -# Make sure the BV-BRC metadata live where they're supposed to + 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 > max_checkm_contam, + flag_checkm_complete = !is.na(.data$checkm_completeness_qc) & .data$checkm_completeness_qc < min_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_rejects <- 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" + ) + } + + reject_list <- list( + make_rejects( + qc$flag_checkm_missing, + "checkm_missing", + ifelse(qc$flag_checkm_missing, "missing", NA_character_), + "present", + "is present" + ), + make_rejects( + qc$flag_checkm_contam, + "checkm_contamination", + qc$checkm_contamination_qc, + max_checkm_contam, + ">" + ), + make_rejects( + qc$flag_checkm_complete, + "checkm_completeness", + qc$checkm_completeness_qc, + min_checkm_complete, + "<" + ) + ) + + if (!is.null(gc_deviations)) { + reject_list[[length(reject_list) + 1L]] <- make_rejects( + qc$flag_gc, + "gc_deviation_sd", + z_gc, + gc_deviations, + ">" + ) + } + + if (!is.null(length_deviations)) { + reject_list[[length(reject_list) + 1L]] <- make_rejects( + qc$flag_length, + "length_deviation_sd", + z_len, + length_deviations, + ">" + ) + } + + if (!is.null(cds_deviations)) { + reject_list[[length(reject_list) + 1L]] <- make_rejects( + qc$flag_cds, + "cds_deviation_sd", + z_cds, + cds_deviations, + ">" + ) + } + + rejections <- dplyr::bind_rows(reject_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 + ) +} + +# 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, and are fresh .ensure_bvbrc_cache <- function(base_dir = ".", verbose = TRUE, + max_age_days = 30L, cache_rel = file.path("data", "bvbrc", "bvbrcData.duckdb"), cache_table = "bvbrc_bac_data") { base_dir <- normalizePath(base_dir, mustWork = FALSE) cache_db <- file.path(base_dir, cache_rel) - need_build <- !file.exists(cache_db) - con_cache <- NULL + # Always delegate to .updateBVBRCdata() so its max_age_days staleness check + # actually runs. Previously this only rebuilt when the cache file/table was + # entirely absent, so an existing-but-stale cache was silently reused + # forever, and results (e.g. genome counts) drifted across machines/sessions + # depending on whenever each cache happened to be built. + .updateBVBRCdata(base_dir = base_dir, max_age_days = max_age_days, verbose = verbose) - if (!need_build) { - con_cache <- DBI::dbConnect(duckdb::duckdb(), dbdir = cache_db) - on.exit(try(DBI::dbDisconnect(con_cache, shutdown = TRUE), silent = TRUE), add = TRUE) - need_build <- !(cache_table %in% DBI::dbListTables(con_cache)) - } - - if (need_build) { - if (isTRUE(verbose)) message("BV-BRC cache missing or incomplete. Building via .updateBVBRCdata(). Please wait.") - .updateBVBRCdata(base_dir = base_dir, verbose = verbose) + if (!file.exists(cache_db)) stop("After .updateBVBRCdata(), cache DB still missing at: ", cache_db) - if (!is.null(con_cache)) try(DBI::dbDisconnect(con_cache, shutdown = TRUE), silent = TRUE) - if (!file.exists(cache_db)) stop("After .updateBVBRCdata(), cache DB still missing at: ", cache_db) - - con_cache <- DBI::dbConnect(duckdb::duckdb(), dbdir = cache_db) - on.exit(try(DBI::dbDisconnect(con_cache, shutdown = TRUE), silent = TRUE), add = TRUE) - if (!(cache_table %in% DBI::dbListTables(con_cache))) { - stop("After .updateBVBRCdata(), table '", cache_table, "' still not found in ", cache_db) - } + con_cache <- DBI::dbConnect(duckdb::duckdb(), dbdir = cache_db, read_only = TRUE) + on.exit(try(DBI::dbDisconnect(con_cache, shutdown = TRUE), silent = TRUE), add = TRUE) + if (!(cache_table %in% DBI::dbListTables(con_cache))) { + stop("After .updateBVBRCdata(), table '", cache_table, "' still not found in ", cache_db) } invisible(cache_db) @@ -910,7 +1183,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 max_checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%). +#' @param min_checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%). +#' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content. +#' @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 +1202,29 @@ retrieveMetadata <- function(user_bacs, abx = "All", overwrite = FALSE, image = "danylmb/bvbrc:5.3", + max_checkm_contam = 5, + min_checkm_complete = 95, + gc_deviations = NULL, + length_deviations = NULL, + cds_deviations = NULL, + debug = FALSE, + export_tables = FALSE, + load_tables = 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 +1234,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 +1249,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 +1272,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,70 +1291,123 @@ 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 + 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) { - .extractAMRtable( - base_dir = base_dir, - batch_genome_IDs = batch, - abx_filter = abx_filter, - drug_fields = drug_fields, - image = image, - verbose = FALSE - ) - }, BPPARAM = param) - combined_drug_data <- unlist(batch_drug_data, use.names = FALSE) - if (length(combined_drug_data) == 0) { + batch_drug_data <- furrr::future_map( + genome_batches, + function(batch) { + raw <- .extractAMRtable( + base_dir = base_dir, + batch_genome_IDs = batch, + abx_filter = abx_filter, + drug_fields = drug_fields, + image = image, + verbose = FALSE + ) + .parse_bvbrc_tsv(raw) + }, + .options = furrr::furrr_options(seed = TRUE) + ) + + combined_drug_data_tbl <- dplyr::bind_rows(batch_drug_data) |> + 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, - microtrait_fields = microtrait_fields, - image = image, - verbose = FALSE - ) - }, BPPARAM = param) - combined_genome_data <- unlist(batch_genome_data, use.names = FALSE) - if (length(combined_genome_data) == 0) { + batch_genome_data <- furrr::future_map( + genome_batches, + function(batch) { + raw <- .extractGenomeData( + base_dir = base_dir, + batch_genome_IDs = batch, + filter_type = filter_type, + amr_fields = amr_fields, + microtrait_fields = microtrait_fields, + image = image, + verbose = FALSE + ) + .parse_bvbrc_tsv(raw) + }, + .options = furrr::furrr_options(seed = TRUE) + ) + + combined_genome_data_tbl <- dplyr::bind_rows(batch_genome_data) |> + 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, + max_checkm_contam = max_checkm_contam, + min_checkm_complete = min_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 ) @@ -1082,48 +1417,88 @@ retrieveMetadata <- function(user_bacs, 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)) { @@ -1138,6 +1513,24 @@ retrieveMetadata <- function(user_bacs, } } + export_res <- NULL + if (isTRUE(export_tables) || isTRUE(load_tables)) { + export_res <- exportTables( + paths$db_path, + export_tables = export_tables, + load_tables = load_tables, + verbose = verbose + ) + } + + if (isTRUE(load_tables)) { + return(list( + duckdbConnection = con, + table_name = "metadata", + data = if (!is.null(export_res)) export_res$data else NULL + )) + } + list(duckdbConnection = con, table_name = "metadata") } @@ -1546,8 +1939,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. @@ -1557,28 +1950,35 @@ 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 + evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"), 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("Found existing 'filtered' table.") con <- con0 tbl <- "filtered" + on.exit(try(DBI::dbDisconnect(con0, shutdown = TRUE), silent = TRUE), add = TRUE) } else { - if (isTRUE(verbose)) message("No 'filtered' table found; filtering now.") + DBI::dbDisconnect(con0, shutdown = TRUE) + + if (isTRUE(verbose)) + message("No 'filtered' table found; filtering now.") f_out <- .filterGenomes( base_dir = base_dir, user_bacs = user_bacs, @@ -1587,8 +1987,10 @@ 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) } + # What genomes need to be downloaded? Build set as `ids` ids <- tibble::as_tibble(DBI::dbReadTable(con, tbl)) |> dplyr::distinct(`genome.genome_id`) |> dplyr::pull(`genome.genome_id`) @@ -1599,74 +2001,109 @@ retrieveGenomes <- function(base_dir = ".", dir.create(genome_path, recursive = TRUE, showWarnings = FALSE) dir.create(logs_dir, recursive = TRUE, showWarnings = FALSE) + # Checks what is already downloaded vs. the full list needed; takes diff 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) } + # Is diff length 0? If so, all genomes ready to go! 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) + 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`) + return(.list_complete(genome_path, all_ids)) } 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), - BPPARAM = ftp_param + if (isTRUE(verbose)) + message("Trying FTPS download. Workers=", ftp_workers) + + ok_ids <- .ftps_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 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( - 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 + + 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.") - # 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( - 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 + + 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 } @@ -1685,6 +2122,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) @@ -1706,6 +2144,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")) @@ -1763,7 +2206,7 @@ genomeList <- function(base_dir = ".", #' /data//.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 +2219,14 @@ 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 max_checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%). +#' @param min_checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%). +#' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content. +#' @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: @@ -1788,7 +2239,16 @@ 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"), + max_checkm_contam = 5, + min_checkm_complete = 95, + gc_deviations = NULL, + length_deviations = NULL, + cds_deviations = NULL, + export_tables = FALSE, + load_tables = FALSE, + debug = FALSE, verbose = TRUE) { method <- match.arg(method) evidence_mode <- match.arg(evidence_mode) @@ -1804,6 +2264,12 @@ prepareGenomes <- function(user_bacs, base_dir = base_dir, abx = "All", overwrite = overwrite, + max_checkm_contam = max_checkm_contam, + min_checkm_complete = min_checkm_complete, + gc_deviations = gc_deviations, + length_deviations = length_deviations, + cds_deviations = cds_deviations, + debug = debug, verbose = verbose )) @@ -1820,13 +2286,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 +2302,15 @@ 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, + ftp_workers = num_workers, + cli_fasta_workers = num_workers, + cli_gff_workers = num_workers, evidence_mode = evidence_mode, - verbose = verbose + verbose = verbose ) if (length(ids) == 0L) { message("No genomes downloaded.") @@ -1850,13 +2319,196 @@ 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 + 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), '")') + } + + export_res <- NULL + if (isTRUE(export_tables) || isTRUE(load_tables)) { + export_res <- exportTables( + paths$db_path, + export_tables = export_tables, + load_tables = load_tables, + verbose = verbose + ) } - out + + if (isTRUE(load_tables)) { + return(list( + duckdb_path = paths$db_path, + table_name = "files", + data = if (!is.null(export_res)) export_res$data else NULL + )) + } + + invisible(out) +} + +#' Export DuckDB tables and optionally load them into R +#' +#' Writes selected DuckDB tables to CSV files and optionally returns them as +#' in-memory R data frames. +#' +#' @param duckdb_path Character. Path to the DuckDB file. +#' @param output_dir Character or NULL. Directory for exports. Defaults to +#' file.path(dirname(duckdb_path), "exports"). +#' @param tables Character vector or NULL. Tables to export. If NULL, exports +#' all tables in the database. +#' @param skip_tables Character vector of table names to exclude. +#' @param include_summary Logical. If TRUE, writes summary.csv and summary.txt. +#' @param load_tables Logical. If TRUE, also return the exported tables as +#' in-memory R data frames. +#' @param verbose Logical. If TRUE, prints progress messages. +#' +#' @return Invisibly returns a list with exported file paths and, if requested, +#' in-memory tables. +#' @keywords internal +exportTables <- function(duckdb_path, + output_dir = NULL, + tables = NULL, + skip_tables = c(NULL), + include_summary = TRUE, + export_tables = TRUE, + load_tables = FALSE, + verbose = TRUE) { + duckdb_path <- normalizePath(duckdb_path, mustWork = TRUE) + + if (is.null(output_dir)) { + output_dir <- file.path(dirname(duckdb_path), "exports") + } + output_dir <- normalizePath(output_dir, mustWork = FALSE) + dir.create(output_dir, recursive = TRUE, showWarnings = FALSE) + + con <- DBI::dbConnect(duckdb::duckdb(), dbdir = duckdb_path) + on.exit(try(DBI::dbDisconnect(con, shutdown = TRUE), silent = TRUE), add = TRUE) + + available_tables <- DBI::dbListTables(con) + if (!length(available_tables)) { + stop("No tables found in DuckDB: ", duckdb_path) + } + + # More tables will exist in a full DuckDB output after data processing, but + # limiting this to the basic genome stats and metadata you'd get from only + # running data_curation.R + basic_tables <- c( + "bac_data", + "filtered", + "metadata", + "genome_data", + "amr_phenotype", + "metadata_qc", + "metadata_qc_rejections" + ) + + if (is.null(tables)) { + tables <- intersect(basic_tables, available_tables) + } else { + tables <- intersect(as.character(tables), available_tables) + if (!length(tables)) { + stop("None of the requested tables were found in the DuckDB.") + } + } + + if (length(skip_tables)) { + tables <- setdiff(tables, skip_tables) + } + if (!length(tables)) { + stop("No tables left to export after applying skip_tables.") + } + + exported_files <- list() + loaded_tables <- list() + + for (tbl in tables) { + out_file <- file.path(output_dir, paste0(tbl, ".csv")) + df <- DBI::dbReadTable(con, tbl) + + if (isTRUE(export_tables)) { + readr::write_csv(df, out_file, na = "") + exported_files[[tbl]] <- out_file + if (isTRUE(verbose)) { + message("Exported table: ", tbl, " -> ", out_file) + } + } + + if (isTRUE(load_tables)) { + loaded_tables[[tbl]] <- df + } + } + + count_if_present <- function(tbl) { + if (tbl %in% available_tables) { + as.character(DBI::dbGetQuery( + con, + paste0("SELECT COUNT(*) AS n FROM ", DBI::dbQuoteIdentifier(con, tbl)) + )$n[[1]]) + } else { + NA_character_ + } + } + + summary_tbl <- tibble::tibble( + metric = c( + "duckdb_path", + "export_dir", + "tables_exported", + "table_names", + "metadata_rows", + "genome_data_rows", + "amr_phenotype_rows", + "metadata_qc_rows", + "metadata_qc_rejections_rows", + "filtered_rows", + "files_rows", + "bac_data_rows" + ), + value = c( + duckdb_path, + output_dir, + as.character(length(tables)), + paste(tables, collapse = ", "), + count_if_present("metadata"), + count_if_present("genome_data"), + count_if_present("amr_phenotype"), + count_if_present("metadata_qc"), + count_if_present("metadata_qc_rejections"), + count_if_present("filtered"), + count_if_present("files"), + count_if_present("bac_data") + ) + ) + + if (isTRUE(export_tables) && isTRUE(include_summary)) { + readr::write_csv(summary_tbl, file.path(output_dir, "summary.csv"), na = "") + writeLines( + c( + paste0("DuckDB: ", duckdb_path), + paste0("Export directory: ", output_dir), + paste0("Tables exported: ", length(tables)), + paste0("Table names: ", paste(tables, collapse = ", ")) + ), + file.path(output_dir, "summary.txt"), + useBytes = TRUE + ) + if (isTRUE(verbose)) { + message("Exported summary files.") + } + } + + invisible(list( + export_dir = output_dir, + tables = tables, + files = exported_files, + data = if (isTRUE(load_tables)) loaded_tables else NULL, + summary = summary_tbl + )) } diff --git a/R/data_processing.R b/R/data_processing.R index 17b7682..c564769 100644 --- a/R/data_processing.R +++ b/R/data_processing.R @@ -203,21 +203,26 @@ 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( + 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, - 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 +1210,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 +1217,35 @@ 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()) + old_plan <- future::plan() + on.exit(future::plan(old_plan), add = TRUE) + future::plan(future::sequential) + + results <- furrr::future_map( + 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 + }, + .options = furrr::furrr_options(seed = TRUE) + ) # Combine results tsvs <- Filter(function(x) !is.null(x) && file.exists(x), results) @@ -1771,7 +1783,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.") @@ -1801,3 +1813,224 @@ runDataProcessing <- function(duckdb_path, parquet_duckdb_path = normalizePath(parquet_duckdb_path) )) } + +#' Export processed tables from DuckDB database +#' +#' Reads tables from the DuckDB database produced by the `runDataProcessing()` workflow +#' and exports them as CSV, TSV, Parquet, and/or XLSX. This is an optional step +#' that allows users to take their processed data outside our amR workflow for +#' use in their own custom analyses. This is not required to run `amRml`! +#' +#' @param duckdb_path Character. Path to the DuckDB database created by the +#' workflow (for example, `Sar.duckdb`). +#' @param output_path Character or NULL. Directory for exports. Defaults to +#' file.path(dirname(duckdb_path), "processed_exports"). +#' @param amr_phenotype_mode Character. One of "separate" or "append". +#' "separate" exports the AMR labels as a separate wide table. +#' "append" joins those labels onto the main feature tables before export. +#' @param export_formats Character vector. Any of "csv", "tsv", "parquet", "xlsx". +#' @param tables Character vector or NULL. Tables to export. If NULL, exports the +#' standard processed tables present in the database. +#' @param verbose Logical. If TRUE, prints progress messages. +#' +#' @return Invisibly returns a list containing the export path, table names, and mode. +#' @export +exportProcessedData <- function(duckdb_path, + output_path = NULL, + amr_phenotype_mode = c("separate", "append"), + export_formats = c("csv"), + export_sequences = FALSE, + tables = NULL, + verbose = TRUE) { + duckdb_path <- normalizePath(duckdb_path, mustWork = TRUE) + amr_phenotype_mode <- match.arg(amr_phenotype_mode) + + export_formats <- unique(tolower(export_formats)) + allowed_formats <- c("csv", "tsv", "parquet", "xlsx") + unknown_formats <- setdiff(export_formats, allowed_formats) + if (length(unknown_formats)) { + stop("Unsupported export format(s): ", paste(unknown_formats, collapse = ", ")) + } + + if ("xlsx" %in% export_formats) { + message( + "Excel spreadsheets of these features can be extremely large and may not open properly even on powerful hardware." + ) + if (!requireNamespace("writexl", quietly = TRUE)) { + stop("Format 'xlsx' was requested but package 'writexl' is not available.") + } + } + + if (is.null(output_path)) { + output_path <- file.path(dirname(duckdb_path), "processed_exports") + } + output_path <- normalizePath(output_path, mustWork = FALSE) + dir.create(output_path, recursive = TRUE, showWarnings = FALSE) + + con <- DBI::dbConnect(duckdb::duckdb(), dbdir = duckdb_path, read_only = TRUE) + DBI::dbExecute( + con, + sprintf( + "SET file_search_path='%s'", + dirname(normalizePath(duckdb_path)) + ) + ) + on.exit(try(DBI::dbDisconnect(con, shutdown = TRUE), silent = TRUE), add = TRUE) + + available_tables <- DBI::dbListTables(con) + if (!length(available_tables)) { + stop("No tables found in DuckDB: ", duckdb_path) + } + + read_tbl <- function(tbl) { + tibble::as_tibble(DBI::dbReadTable(con, tbl)) + } + + write_one <- function(df, stem) { + if ("csv" %in% export_formats) { + readr::write_csv(df, file.path(output_path, paste0(stem, ".csv")), na = "") + } + if ("tsv" %in% export_formats) { + readr::write_tsv(df, file.path(output_path, paste0(stem, ".tsv")), na = "") + } + if ("parquet" %in% export_formats) { + arrow::write_parquet(df, file.path(output_path, paste0(stem, ".parquet"))) + } + if ("xlsx" %in% export_formats) { + writexl::write_xlsx(list(data = df), file.path(output_path, paste0(stem, ".xlsx"))) + } + } + + build_amr_wide <- function() { + source_tbl <- if ("metadata" %in% available_tables) { + "metadata" + } else if ("amr_phenotype" %in% available_tables) { + "amr_phenotype" + } else { + NULL + } + + if (is.null(source_tbl)) { + return(NULL) + } + + md <- read_tbl(source_tbl) + + needed <- c("genome.genome_id", "genome_drug.antibiotic", "genome_drug.resistant_phenotype") + if (!all(needed %in% names(md))) { + return(NULL) + } + + md |> + dplyr::transmute( + genome_id = `genome.genome_id`, + antibiotic = `genome_drug.antibiotic`, + phenotype = `genome_drug.resistant_phenotype` + ) |> + dplyr::filter(!is.na(genome_id), !is.na(antibiotic), !is.na(phenotype)) |> + dplyr::distinct() |> + dplyr::group_by(genome_id, antibiotic) |> + dplyr::summarise( + phenotype = paste(sort(unique(phenotype)), collapse = ";"), + .groups = "drop" + ) |> + tidyr::pivot_wider( + names_from = antibiotic, + values_from = phenotype, + values_fill = NA_character_ + ) |> + dplyr::arrange(genome_id) + } + + # Appendable here means whether we glue AST phenotypes on the end or not + table_specs <- list( + gene_count = list(source = "gene_count", stem = "gene_count", appendable = TRUE), + protein_count = list(source = "protein_count", stem = "protein_count", appendable = TRUE), + domain_count = list(source = "domain_count", stem = "domain_count", appendable = TRUE), + struct = list(source = "struct", stem = "struct", appendable = TRUE), + gene_names = list(source = "gene_names", stem = "gene_names", appendable = FALSE), + protein_names = list(source = "protein_names", stem = "protein_names", appendable = FALSE), + domain_names = list(source = "domain_names", stem = "domain_names", appendable = FALSE), + metadata = list(source = "metadata", stem = "metadata", appendable = FALSE), + genome_data = list(source = "genome_data", stem = "genome_data", appendable = FALSE), + amr_phenotype_wide = list(source = NULL, stem = "amr_phenotype_wide", appendable = FALSE) + ) + + if (isTRUE(export_sequences)) { + table_specs$gene_seqs <- list(source = "gene_ref_seq", stem = "gene_seqs", appendable = FALSE) + table_specs$protein_seqs <- list(source = "protein_cluster_seq", stem = "protein_seqs", appendable = FALSE) + table_specs$genome_gene_protein <- list(source = "genome_gene_protein", stem = "genome_gene_protein", appendable = FALSE) + } + + if (is.null(tables)) { + selected_keys <- c( + "gene_count", + "protein_count", + "domain_count", + "struct", + "gene_names", + "protein_names", + "domain_names", + "metadata", + "genome_data", + "amr_phenotype_wide" + ) + if (isTRUE(export_sequences)) { + selected_keys <- c(selected_keys, "gene_seqs", "protein_seqs", "genome_gene_protein") + } + } else { + selected_keys <- intersect(as.character(tables), names(table_specs)) + } + + if (!length(selected_keys)) { + stop("No requested tables were found in the DuckDB.") + } + + phenotype_wide <- build_amr_wide() + exported <- character(0) + + for (key in selected_keys) { + spec <- table_specs[[key]] + + if (key == "amr_phenotype_wide") { + if (is.null(phenotype_wide)) { + if (isTRUE(verbose)) message("Skipping amr_phenotype_wide: no AMR source table found.") + next + } + write_one(phenotype_wide, spec$stem) + exported <- c(exported, spec$stem) + if (isTRUE(verbose)) message("Exported: ", spec$stem) + next + } + + if (is.null(spec$source) || !(spec$source %in% available_tables)) { + if (isTRUE(verbose)) message("Skipping missing table: ", key) + next + } + + df <- read_tbl(spec$source) + out_stem <- spec$stem + + if (identical(amr_phenotype_mode, "append") && + isTRUE(spec$appendable) && + !is.null(phenotype_wide) && + "genome_id" %in% names(df)) { + df <- dplyr::left_join(df, phenotype_wide, by = "genome_id") + out_stem <- paste0(spec$stem, "_with_phenotypes") + } + + write_one(df, out_stem) + exported <- c(exported, out_stem) + + if (isTRUE(verbose)) message("Exported: ", out_stem) + } + + invisible(list( + duckdb_path = duckdb_path, + output_path = output_path, + tables = exported, + amr_phenotype_mode = amr_phenotype_mode, + export_formats = export_formats, + export_sequences = isTRUE(export_sequences) + )) +} 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: diff --git a/inst/extdata/Campy_testdata.zip b/inst/extdata/Campy_testdata.zip new file mode 100644 index 0000000..32a4e05 Binary files /dev/null and b/inst/extdata/Campy_testdata.zip differ