diff --git a/NAMESPACE b/NAMESPACE index aedaaa0..e43ed5d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,6 +5,7 @@ export(.updateBVBRCdata) export(CDHIT2duckdb) export(buildClusterFeatureMap) export(buildDyadFeatureMap) +export(checkDataAvailability) export(cleanData) export(cleanMetaData) export(exportProcessedData) diff --git a/R/data_curation.R b/R/data_curation.R index 189555d..0b95aee 100644 --- a/R/data_curation.R +++ b/R/data_curation.R @@ -2645,3 +2645,95 @@ exportTables <- function(duckdb_path, summary = summary_tbl )) } + +#' Check BV-BRC data availability for selected bacterial taxa +#' +#' Performs a lightweight metadata-only query to summarize the data available +#' for one or more bacterial taxa before running the genome-download and +#' feature-processing workflows. +#' +#' By default, genome IDs and metadata are resolved through the BV-BRC Data API. +#' The legacy CLI path can be selected with `metadata_method = "cli"`. +#' +#' @param user_bacs Character vector. Mixed taxon IDs and/or species names. +#' @param base_dir Character. Project root. Used by the legacy CLI path. +#' Default `"."`. +#' @param metadata_method Character. Metadata backend: `"api"` (default) or +#' `"cli"`. +#' @param max_checkm_contam Numeric. Maximum allowed CheckM contamination +#' (%) for the QC summary. Default `5`. +#' @param min_checkm_complete Numeric. Minimum allowed CheckM completeness +#' (%) for the QC summary. Default `95`. +#' @param gc_deviations Numeric. Maximum SDs from the median GC content. +#' Optional. Default `NULL`. +#' @param length_deviations Numeric. Maximum SDs from the median genome +#' length. Optional. Default `NULL`. +#' @param cds_deviations Numeric. Maximum SDs from the median CDS count. +#' Optional. Default `NULL`. +#' @param verbose Logical. If TRUE, print progress messages. Default `TRUE`. +#' +#' @return A tibble with one row per requested taxon containing summary +#' statistics describing genome availability, sequencing status, AMR data +#' availability, and the number of genomes passing metadata QC. +#' +#' @examples +#' \dontrun{ +#' checkDataAvailability( +#' c("Staphylococcus argenteus", "Streptococcus suis") +#' ) +#' +#' checkDataAvailability( +#' c("1280", "28901"), +#' metadata_method = "api" +#' ) +#' } +#' +#' @export +checkDataAvailability <- function( + user_bacs, + base_dir = ".", + metadata_method = c("api", "cli"), + max_checkm_contam = 5, + min_checkm_complete = 95, + gc_deviations = NULL, + length_deviations = NULL, + cds_deviations = NULL, + verbose = TRUE +) { + metadata_method <- match.arg(metadata_method) + base_dir <- normalizePath(base_dir, mustWork = FALSE) + + if (missing(user_bacs) || length(user_bacs) == 0L) { + stop("`user_bacs` must contain at least one taxon ID or species name.") + } + + user_bacs <- unique(trimws(as.character(user_bacs))) + user_bacs <- user_bacs[nzchar(user_bacs)] + + if (!length(user_bacs)) { + stop("`user_bacs` must contain at least one non-empty taxon ID or species name.") + } + + if (isTRUE(verbose)) { + message( + "Checking BV-BRC data availability using metadata_method = \"", + metadata_method, + "\"." + ) + } + + dplyr::bind_rows( + purrr::map( + user_bacs, + .checkDataPerTaxon, + base_dir = base_dir, + metadata_method = metadata_method, + max_checkm_contam = max_checkm_contam, + min_checkm_complete = min_checkm_complete, + gc_deviations = gc_deviations, + length_deviations = length_deviations, + cds_deviations = cds_deviations, + verbose = verbose + ) + ) +} diff --git a/R/data_processing.R b/R/data_processing.R index 45875c2..4825c16 100644 --- a/R/data_processing.R +++ b/R/data_processing.R @@ -962,10 +962,17 @@ CDHIT2duckdb <- function(duckdb_path, #' Download and prepare HMMER databases for generating new file types. #' -#' @param hmmer_db_dir Directory to store HMMER databases -#' @param databases List of databases to prepare (default: c("Pfam", "COG", "AMRFinder")) -#' @param docker_image Docker image containing HMMER (default: "staphb/hmmer") -#' @param hmmer_db_url If the databases contain custom database(s), the url is required to download the database. +#' @param hmmer_db_dir Character. Directory where HMMER databases are cached. +#' @param databases Character vector of database names to prepare. Supported +#' built-in databases are `Pfam`, `COG`, and `AMRFinder`; `DefenseCas` is +#' prepared separately by the DefenseFinder/CasFinder workflow. +#' @param docker_image Character. Docker image containing HMMER tools used to +#' press the prepared databases. Default: `"staphb/hmmer"`. +#' @param hmmer_db_url NON-FUNCTIONAL. Character or `NULL`. URL used to download +#' a custom HMMER database when `databases` contains names not covered by the +#' built-in database definitions. This function is not currently active! +#' @param verbose Logical. Print status messages while checking, downloading, +#' combining, and pressing databases. Default: `TRUE`. #' #' @returns A list of paths to the database hmm files. #' @@ -1280,20 +1287,28 @@ CDHIT2duckdb <- function(duckdb_path, #' The function to run HMMER with docker #' -#' @param JOB_NAME protein_chunk id -#' @param FASTA fasta sequences -#' @param DB HMM database -#' @param Total_proteins protein sequence count -#' @param output_path path for saving hmmer outputs -#' @param db_paths path to HMM database -#' @param docker_image hmmer docker image (ideally from dockerhub) -#' @param threads number of threads -#' @param n_workers number of parallel workers +##' @param JOB_NAME Character. Identifier used for the HMMER job and output +#' filename. +#' @param FASTA Character. File name of the protein FASTA chunk to search. +#' @param DB Character. Name of the HMMER database to search against. +#' @param total_proteins Integer. Total number of proteins in the full input +#' dataset, used to set HMMER's `-Z` and `--domZ` values. +#' @param output_path Character. Directory containing the FASTA input, HMMER +#' output, and final Parquet result. +#' @param db_paths List. Prepared HMMER database metadata indexed by database +#' name. +#' @param docker_image Character. Docker image containing HMMER. Default: +#' `"staphb/hmmer"`. +#' @param threads Integer. Total CPU budget used when calculating the number +#' of threads allocated to this job. Default: `8`. +#' @param n_workers Integer. Number of HMMER jobs being run in parallel. +#' Used to divide the CPU budget among jobs. Default: `8`. +#' @param verbose Logical. Print progress messages. Default: `TRUE`. #' #' @returns the filename of the parquet file with hmmer output post parsing #' #' @keywords internal -.runHmmerJob <- function(JOB_NAME, FASTA, DB, Total_proteins, +.runHmmerJob <- function(JOB_NAME, FASTA, DB, total_proteins, output_path = NULL, db_paths, docker_image = "staphb/hmmer", threads = 8L, n_workers = 8L, @@ -1326,8 +1341,8 @@ CDHIT2duckdb <- function(duckdb_path, "hmmsearch", "--notextw", "--cpu", as.character(threads_per_job), - "-Z", Total_proteins, - "--domZ", Total_proteins, + "-Z", total_proteins, + "--domZ", total_proteins, "--domtblout", .to_container(hmmer_output, mount_host, mount_cont), db_cont_path, .to_container(hmmer_input, mount_host, mount_cont) @@ -1374,14 +1389,23 @@ CDHIT2duckdb <- function(duckdb_path, #' Wrapper for preparing HMM databases and running HMMER on protein sequences from duckdb and writing them. #' -#' @param duckdb_path path to the duckdb with protein sequences and list -#' @param output_path path where HMMER output will be saved -#' @param threads number of threads -#' @param hmmer_db_dir path to the directory where HMM databases are/will be downloaded -#' @param databases list of HMM databases -#' @param docker_image the docker image of HMMER -#' @param num_of_splits The number of splits of the protein sequence file for parallel processing -#' @param n_workers The number of parallel runs +#' @param duckdb_path Character. Path to the DuckDB database containing +#' `protein_cluster_seq`, which provides the protein sequences to analyze. +#' @param output_path Character. Directory for HMMER intermediate and final +#' Parquet outputs. +#' @param threads Integer. Total CPU budget used by HMMER jobs. Default: `8`. +#' @param hmmer_db_dir Character. Directory containing the prepared HMMER +#' databases. If `NULL`, the default `amRdata` HMMER database cache is used. +#' @param databases Character vector of HMMER databases to run. +#' @param docker_image Character. Docker image containing HMMER. Default: +#' `"staphb/hmmer"`. +#' @param num_of_splits Integer. Number of chunks into which the protein +#' sequences should be divided. Must be a positive integer. The requested +#' value is automatically reduced when fewer protein sequences are available. +#' Default: `8`. +#' @param n_workers Integer. Number of parallel HMMER jobs to run. Default: `8`. +#' @param verbose Logical. Print progress messages. Default: `TRUE`.#' +#' @returns #' #' @keywords internal .runHMMER <- function(duckdb_path, @@ -1433,7 +1457,7 @@ CDHIT2duckdb <- function(duckdb_path, } # required to define the database size for hmmsearch --Z and --domZ parameters - Total_proteins <- nrow(prot_seqs) + total_proteins <- nrow(prot_seqs) if (is.null(hmmer_db_dir)) { hmmer_db_dir <- .defaultHmmerDbDir() @@ -1509,7 +1533,7 @@ CDHIT2duckdb <- function(duckdb_path, JOB_NAME = job_list$JOB_NAME[i], FASTA = job_list$FASTA[i], DB = job_list$DB[i], - Total_proteins = Total_proteins, + total_proteins = total_proteins, output_path = output_path, db_paths = db_paths, docker_image = docker_image, @@ -1608,11 +1632,13 @@ CDHIT2duckdb <- function(duckdb_path, #' counts per genome and annotation, and writes the result both as a Parquet file #' and as a new table in the DuckDB database. #' -#' @param annotated_parquet Path to the combined HMMER results Parquet file -#' (e.g. `"results/Ecoli/protein_COG.parquet"`). The filename stem is used as -#' the table name in DuckDB. -#' @param duckdb_path Path to the per-selection DuckDB database containing a -#' `protein_count` table (created by [CDHIT2duckdb()]). +#' @param duckdb_path Character. Path to the per-selection DuckDB database +#' containing the `protein_count` table created by [CDHIT2duckdb()]. +#' @param databases Character vector of HMMER database names to process. +#' Each database must correspond to a `protein_` annotation table +#' already present in the DuckDB. +#' @param output_path Character. Directory where the genome-by-annotation +#' Parquet files will be written. Defaults to `dirname(duckdb_path)`. #' #' @return Invisibly returns the path to the written count Parquet file. #' @@ -2052,7 +2078,7 @@ CDHIT2duckdb <- function(duckdb_path, ) # required to define the database size for hmmsearch --Z and --domZ parameters - Total_proteins <- nrow(prot_seqs) + total_proteins <- nrow(prot_seqs) readr::write_lines( paste0( @@ -2107,8 +2133,8 @@ CDHIT2duckdb <- function(duckdb_path, "--notextw", "--cpu", as.character(threads), - "-Z", Total_proteins, - "--domZ", Total_proteins, + "-Z", total_proteins, + "--domZ", total_proteins, "--domtblout", file.path( "/work", @@ -2628,6 +2654,9 @@ cleanData <- function(duckdb_path, path) { #' \item **Metadata cleaning + Parquet export** via [cleanData()] -> writes #' Parquet files to `output_path`, and builds a **Parquet-backed DuckDB** #' (`*_parquet.duckdb`) with views over those Parquets. +#' +#' \item **Dyad feature mapping** via [buildDyadFeatureMap()] -> writes the +#' protein-gene dyad feature network used for downstream graph analysis. #' } #' #' @param duckdb_path Character. Path to the **per-selection DuckDB** produced by @@ -2641,6 +2670,10 @@ cleanData <- function(duckdb_path, path) { #' @param threads Integer. Shared concurrency budget used across Panaroo, CD-HIT, #' and HMMER. Defaults to `8`. #' +#' @param export_tabular_data Logical. If TRUE, automatically call +#' [exportProcessedData()] after the processing pipeline completes to create +#' the default human-readable exports. Default: `FALSE`. +#' #' @param panaroo_split_jobs Logical. If `TRUE`, Panaroo runs in multiple batches #' that can be merged by [.mergePanaroo()]. If `FALSE`, Panaroo runs once on all #' isolates. Default: `FALSE`. @@ -2710,6 +2743,10 @@ cleanData <- function(duckdb_path, path) { #' over the generated Parquet files. #' * Records processing parameters, software versions, database selections, and #' other provenance information in the dataset manifest. +#' * Creates the protein-gene dyad feature map for downstream graph analysis. +#' * If `export_tabular_data = TRUE`, exports features in human-readable CSVs. +#' Additional export options are available by calling `exportProcessedData()` +#' after `runDataProcessing()` completes. #' #' **Threading** #' * `threads` provides the shared CPU budget for the major processing stages. @@ -2739,6 +2776,7 @@ runDataProcessing <- function( duckdb_path, output_path = NULL, threads = 8, + export_tabular_data = FALSE, # Panaroo panaroo_split_jobs = FALSE, @@ -3115,6 +3153,24 @@ runDataProcessing <- function( "_parquet.duckdb" ) + # With all features completed, create the dyad feature map +if (isTRUE(verbose)) message("Building the mapping of protein|gene dyad to all features and exporting Parquet-backed views.") + buildDyadFeatureMap(duckdb_path = duckdb_path, output_path = out_dir) + + # And if the user wants to export processed data + if (export_tabular_data == TRUE) { + if (isTRUE(verbose)) + message("\n============================================") + message("Exporting human-readable processed data tables.") + message("Additional export options are available through `exportProcessedData()`.") + message("\n============================================") + + # Export with default parameters + exportProcessedData(duckdb_path = duckdb_path, + output_path = out_dir, + verbose = verbose) + } + if (isTRUE(verbose)) { message("\n============================================") message("Completed data-processing workflow successfully.") @@ -3127,6 +3183,45 @@ runDataProcessing <- function( message("============================================\n") } + # Final Parquets generated by this run + parquet_files <- c( + # Metadata outputs + file.path(out_dir, "metadata.parquet"), + file.path(out_dir, "amr_phenotype.parquet"), + file.path(out_dir, "genome_data.parquet"), + file.path(out_dir, "original_metadata.parquet"), + + # Core feature outputs + file.path(out_dir, "gene_count.parquet"), + file.path(out_dir, "gene_names.parquet"), + file.path(out_dir, "gene_seqs.parquet"), + file.path(out_dir, "genome_gene_protein.parquet"), + file.path(out_dir, "struct.parquet"), + file.path(out_dir, "protein_count.parquet"), + file.path(out_dir, "protein_names.parquet"), + file.path(out_dir, "protein_seqs.parquet"), + file.path(out_dir, "protein_members.parquet"), + + # HMMER outputs + file.path( + out_dir, + paste0("protein_", hmmer_databases, ".parquet") + ), + file.path( + out_dir, + paste0("protein_", hmmer_databases, "_count.parquet") + ), + + # Dyad feature map + file.path(out_dir, "dyad_feature.parquet") + ) + + # Record only files that exist, so manifest reflects successful runs + parquet_files <- normalizePath( + parquet_files[file.exists(parquet_files)], + mustWork = TRUE + ) + # Log! manifest <- .manifest_stage( manifest, @@ -3143,11 +3238,12 @@ runDataProcessing <- function( ref_file_path ), outputs = c( - out_dir, + parquet_files, parquet_duckdb_path ), metrics = list( - parquet_duckdb = parquet_duckdb_path + parquet_duckdb = parquet_duckdb_path, + parquet_files = parquet_files ) ) @@ -3180,8 +3276,17 @@ runDataProcessing <- function( #' "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 export_sequences Logical. If TRUE, also exports gene and protein +#' sequence tables and the genome-to-gene-to-protein mapping. Default FALSE. +#' @param export_dyads Logical. If TRUE, exports the optional dyad annotation +#' table. Each row represents a protein-gene dyad with semicolon-separated +#' mapped feature values for structural annotations and the HMMER databases +#' recorded in the dataset manifest. The export uses `export_formats`. +#' Default FALSE. +#' @param tables Character vector or NULL. Tables to export. If NULL, exports +#' the standard processed tables plus HMMER tables recorded in the manifest. +#' @param export_tables Logical. If TRUE, write the selected tables to disk. +#' Default TRUE. #' @param verbose Logical. If TRUE, prints progress messages. #' #' @return Invisibly returns a list containing the export path, table names, and mode. @@ -3191,6 +3296,7 @@ exportProcessedData <- function(duckdb_path, amr_phenotype_mode = c("separate", "append"), export_formats = c("csv"), export_sequences = FALSE, + export_dyads = FALSE, tables = NULL, export_tables = TRUE, verbose = TRUE) { @@ -3264,14 +3370,14 @@ exportProcessedData <- function(duckdb_path, df, file = file.path(output_path, paste0(stem, ".csv")), sep = ",", row.names = FALSE, col.names = TRUE, quote = TRUE, na = "", qmethod = "double", fileEncoding = "UTF-8" - ) + ) } if ("tsv" %in% export_formats) { utils::write.table( df, file = file.path(output_path, paste0(stem, ".tsv")), sep = "\t", row.names = FALSE, col.names = TRUE, quote = TRUE, na = "", qmethod = "double", fileEncoding = "UTF-8" - ) + ) } if ("parquet" %in% export_formats) { arrow::write_parquet(df, file.path(output_path, paste0(stem, ".parquet"))) @@ -3281,6 +3387,56 @@ exportProcessedData <- function(duckdb_path, } } + # Determine which HMMER databases were actually run from the latest + # successful `runDataProcessing()` manifest + manifest_path <- .manifest_find_latest(duckdb_path) + + hmmer_databases <- character() + + if (!is.null(manifest_path)) { + manifest <- jsonlite::read_json( + manifest_path, + simplifyVector = FALSE + ) + + successful_hmmer <- list() + + for (run in rev(manifest$runs %||% list())) { + stages <- run$stages %||% list() + + matches <- purrr::keep( + stages, + ~ identical(.x$name, "hmmer") && + identical(.x$status, "success") + ) + + if (length(matches)) { + successful_hmmer <- matches[[1]] + break + } + } + + if (length(successful_hmmer)) { + hmmer_databases <- unlist( + successful_hmmer$parameters$databases %||% character(), + use.names = FALSE + ) + hmmer_databases <- unique(as.character(hmmer_databases)) + } + } + + if (!length(hmmer_databases) && isTRUE(verbose)) { + message( + "No successful HMMER runs were found in the manifest. ", + "HMMER tables will not be selected automatically." + ) + } else if (isTRUE(verbose)) { + message( + "Successful HMMER runs were identified in the manifest: ", + paste(hmmer_databases, collapse = ", ") + ) + } + build_amr_wide <- function() { source_tbl <- if ("metadata" %in% available_tables) { "metadata" @@ -3325,35 +3481,75 @@ exportProcessedData <- function(duckdb_path, 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 = "gene_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) ) + # Add only the HMMER databases recorded in the manifest. + for (database in hmmer_databases) { + annotation_key <- paste0("protein_", database) + count_key <- paste0(annotation_key, "_count") + + table_specs[[annotation_key]] <- list( + source = annotation_key, + stem = annotation_key, + appendable = FALSE + ) + + table_specs[[count_key]] <- list( + source = count_key, + stem = count_key, + appendable = TRUE + ) + } + 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) + 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" + "amr_phenotype_wide", + paste0( + "protein_", + hmmer_databases, + "_count" + ), + paste0( + "protein_", + hmmer_databases + ) ) + + selected_keys <- selected_keys[selected_keys %in% names(table_specs)] + if (isTRUE(export_sequences)) { selected_keys <- c(selected_keys, "gene_seqs", "protein_seqs", "genome_gene_protein") } @@ -3372,6 +3568,30 @@ exportProcessedData <- function(duckdb_path, } exported <- character(0) + # Optionally export a mapping table rooted on the protein-gene dyads + if (isTRUE(export_dyads)) { + dyad_tbl <- .exportDyadAnnotations( + duckdb_path = duckdb_path, + verbose = verbose + ) + + if (isTRUE(export_tables)) { + write_one( + dyad_tbl, + "dyad_annotations" + ) + + exported <- c( + exported, + "dyad_annotations" + ) + + if (isTRUE(verbose)) { + message("Exported: dyad_annotations") + } + } + } + for (key in selected_keys) { spec <- table_specs[[key]] @@ -3403,10 +3623,14 @@ exportProcessedData <- function(duckdb_path, out_stem <- paste0(spec$stem, "_with_phenotypes") } - write_one(df, out_stem) - exported <- c(exported, out_stem) + if (isTRUE(export_tables)) { + write_one(df, out_stem) + exported <- c(exported, out_stem) - if (isTRUE(verbose)) message("Exported: ", out_stem) + if (isTRUE(verbose)) { + message("Exported: ", out_stem) + } + } } invisible(list( @@ -3415,6 +3639,8 @@ exportProcessedData <- function(duckdb_path, tables = exported, amr_phenotype_mode = amr_phenotype_mode, export_formats = export_formats, - export_sequences = isTRUE(export_sequences) + export_sequences = isTRUE(export_sequences), + hmmer_databases = hmmer_databases, + manifest_path = manifest_path )) } diff --git a/R/helpers.R b/R/helpers.R index 76c3287..7ad8842 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -96,6 +96,547 @@ invisible(TRUE) } +#' Check BV-BRC data availability for a single taxon +#' +#' Internal worker used by `checkDataAvailability()`. Each call resolves and +#' summarizes a taxon independently. +#' +#' @param user_bac Character scalar. Taxon ID or species name. +#' @inheritParams checkDataAvailability +#' +#' @return A one-row tibble containing genome and AMR availability statistics. +#' @keywords internal +.checkDataPerTaxon <- function( + user_bac, + base_dir = ".", + metadata_method = c("api", "cli"), + max_checkm_contam = 5, + min_checkm_complete = 95, + gc_deviations = NULL, + length_deviations = NULL, + cds_deviations = NULL, + verbose = TRUE +) { + metadata_method <- match.arg(metadata_method) + base_dir <- normalizePath(base_dir, mustWork = FALSE) + + # Little cache of internal helper helpers to help the helper + empty_result <- function() { + tibble::tibble( + query = user_bac, + total_genomes = 0L, + wgs_genomes = 0L, + complete_genomes = 0L, + amr_genomes = 0L, + amr_records = 0L, + unique_antibiotics = 0L, + antibiotics = NA_character_, + drug_classes = NA_character_, + checkm_available = 0L, + qc_pass_genomes = 0L, + qc_fail_genomes = 0L, + median_genome_length = NA_real_, + median_gc_content = NA_real_, + median_cds = NA_real_, + median_checkm_completeness = NA_real_, + median_checkm_contamination = NA_real_, + collection_year_min = NA_integer_, + collection_year_max = NA_integer_ + ) + } + + safe_median <- function(x) { + x <- suppressWarnings(as.numeric(x)) + x <- x[is.finite(x)] + + if (!length(x)) { + return(NA_real_) + } + + stats::median(x) + } + + collapse_unique <- function(x) { + x <- trimws(as.character(x)) + x <- sort(unique(x[!is.na(x) & nzchar(x)])) + + if (!length(x)) { + return(NA_character_) + } + + paste(x, collapse = ", ") + } + + # Resolve this taxon to genome IDs + genome_ids <- if (identical(metadata_method, "api")) { + .resolveGenomeIDsApi( + base_dir = base_dir, + user_bacs = user_bac, + verbose = verbose + ) + } else { + # Legacy CLI shenanigans + bac_input_data <- .retrieveCustomQuery( + base_dir = base_dir, + user_bacs = user_bac + ) + + if (is.null(bac_input_data) || nrow(bac_input_data) == 0L) { + character(0) + } else { + cache_db <- file.path( + base_dir, + "data", + "bvbrc", + "bvbrcData.duckdb" + ) + + if (!file.exists(cache_db)) { + stop( + "BV-BRC cache not found at: ", + cache_db, + ". Run .updateBVBRCdata() first." + ) + } + + 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 + ) + + taxon_ids <- unique(bac_input_data$genome.taxon_id) + taxon_sql <- paste( + DBI::dbQuoteString(con_cache, taxon_ids), + collapse = ", " + ) + + query <- sprintf( + paste0( + "SELECT DISTINCT \"genome.genome_id\" AS genome_id ", + "FROM bvbrc_bac_data ", + "WHERE \"genome.taxon_id\" IN (%s)" + ), + taxon_sql + ) + + result <- DBI::dbGetQuery(con_cache, query) + + if (nrow(result)) { + unique(as.character(result$genome_id)) + } else { + character(0) + } + } + } + + genome_ids <- unique(as.character(genome_ids)) + genome_ids <- genome_ids[ + !is.na(genome_ids) & nzchar(genome_ids) + ] + + if (!length(genome_ids)) { + if (isTRUE(verbose)) { + message("No genomes matched '", user_bac, "'.") + } + return(empty_result()) + } + + # Fetch genome metadata + genome_fields <- paste( + c( + "genome_id", + "genome_name", + "species", + "taxon_id", + "genome_quality", + "genome_status", + "collection_year", + "genome_length", + "gc_content", + "cds", + "checkm_completeness", + "checkm_contamination" + ), + collapse = "," + ) + + genome_data <- if (identical(metadata_method, "api")) { + .extractGenomeDataApi( + genome_ids = genome_ids, + fields = genome_fields, + verbose = verbose + ) + } else { + raw <- .extractGenomeData( + base_dir = base_dir, + batch_genome_IDs = genome_ids, + filter_type = "AMR", + amr_fields = genome_fields, + microtrait_fields = genome_fields, + verbose = verbose + ) + + .parse_bvbrc_tsv(raw) + } + + genome_data <- tibble::as_tibble(genome_data) + + if (!nrow(genome_data)) { + if (isTRUE(verbose)) { + message( + "No genome metadata were returned for '", + user_bac, + "'." + ) + } + return(empty_result()) + } + + id_col <- dplyr::case_when( + "genome.genome_id" %in% names(genome_data) ~ + "genome.genome_id", + "genome_id" %in% names(genome_data) ~ + "genome_id", + TRUE ~ NA_character_ + ) + + if (is.na(id_col)) { + stop("Genome metadata did not contain a proper genome ID column.") + } + + genome_data <- genome_data |> + dplyr::mutate( + .genome_id = as.character(.data[[id_col]]) + ) |> + dplyr::filter( + !is.na(.data$.genome_id), + nzchar(.data$.genome_id) + ) |> + dplyr::distinct(.data$.genome_id, .keep_all = TRUE) + + # Retrieve AMR phenotype records + if (isTRUE(verbose)) { + message("Checking AMR phenotype availability for '", user_bac, "'.") + } + + amr_data <- if (identical(metadata_method, "api")) { + .extractAMRtableApi( + genome_ids = genome_ids, + abx = "All", + verbose = verbose + ) + } else { + drug_fields <- paste( + c( + "genome_id", + "antibiotic", + "evidence", + "laboratory_typing_method", + "resistant_phenotype" + ), + collapse = "," + ) + + raw <- .extractAMRtable( + base_dir = base_dir, + batch_genome_IDs = genome_ids, + abx_filter = "--required antibiotic", + drug_fields = drug_fields, + verbose = verbose + ) + + .parse_bvbrc_tsv(raw) + } + + amr_data <- tibble::as_tibble(amr_data) + + amr_id_col <- dplyr::case_when( + "genome_drug.genome_id" %in% names(amr_data) ~ + "genome_drug.genome_id", + "genome_id" %in% names(amr_data) ~ + "genome_id", + TRUE ~ NA_character_ + ) + + amr_ids <- if (!is.na(amr_id_col) && nrow(amr_data)) { + unique(as.character(amr_data[[amr_id_col]])) + } else { + character(0) + } + + amr_ids <- amr_ids[ + !is.na(amr_ids) & nzchar(amr_ids) + ] + + antibiotic_col <- dplyr::case_when( + "genome_drug.antibiotic" %in% names(amr_data) ~ + "genome_drug.antibiotic", + "antibiotic" %in% names(amr_data) ~ + "antibiotic", + TRUE ~ NA_character_ + ) + + genome_length_col <- dplyr::case_when( + "genome.genome_length" %in% names(genome_data) ~ + "genome.genome_length", + "genome_length" %in% names(genome_data) ~ + "genome_length", + TRUE ~ NA_character_ + ) + + gc_content_col <- dplyr::case_when( + "genome.gc_content" %in% names(genome_data) ~ + "genome.gc_content", + "gc_content" %in% names(genome_data) ~ + "gc_content", + TRUE ~ NA_character_ + ) + + cds_col <- dplyr::case_when( + "genome.cds" %in% names(genome_data) ~ + "genome.cds", + "cds" %in% names(genome_data) ~ + "cds", + TRUE ~ NA_character_ + ) + + # Apply the same metadata QC used by retrieveMetadata() + qc_out <- .apply_metadata_qc( + genome_tbl = genome_data, + max_checkm_contam = max_checkm_contam, + min_checkm_complete = min_checkm_complete, + gc_deviations = gc_deviations, + length_deviations = length_deviations, + cds_deviations = cds_deviations + ) + + qc_tbl <- tibble::as_tibble(qc_out$qc_tbl) + + total_genomes <- nrow(genome_data) + + wgs_genomes <- if ( + "genome.genome_status" %in% names(genome_data) + ) { + sum( + genome_data$genome.genome_status == "WGS", + na.rm = TRUE + ) + } else if ("genome_status" %in% names(genome_data)) { + sum(genome_data$genome_status == "WGS", na.rm = TRUE) + } else { + NA_integer_ + } + + complete_genomes <- if ( + "genome.genome_status" %in% names(genome_data) + ) { + sum( + genome_data$genome.genome_status == "Complete", + na.rm = TRUE + ) + } else if ("genome_status" %in% names(genome_data)) { + sum( + genome_data$genome_status == "Complete", + na.rm = TRUE + ) + } else { + NA_integer_ + } + + checkm_complete_col <- if ( + "genome.checkm_completeness" %in% names(genome_data) + ) { + "genome.checkm_completeness" + } else if ("checkm_completeness" %in% names(genome_data)) { + "checkm_completeness" + } else { + NA_character_ + } + + checkm_contam_col <- if ( + "genome.checkm_contamination" %in% names(genome_data) + ) { + "genome.checkm_contamination" + } else if ("checkm_contamination" %in% names(genome_data)) { + "checkm_contamination" + } else { + NA_character_ + } + + checkm_available <- if ( + !is.na(checkm_complete_col) && + !is.na(checkm_contam_col) + ) { + sum( + !is.na(suppressWarnings( + as.numeric(genome_data[[checkm_complete_col]]) + )) & + !is.na(suppressWarnings( + as.numeric(genome_data[[checkm_contam_col]]) + )) + ) + } else { + NA_integer_ + } + + qc_pass <- sum(qc_tbl$qc_keep %in% TRUE, na.rm = TRUE) + qc_fail <- sum(qc_tbl$qc_keep %in% FALSE, na.rm = TRUE) + + collection_year_col <- if ( + "genome.collection_year" %in% names(genome_data) + ) { + "genome.collection_year" + } else if ("collection_year" %in% names(genome_data)) { + "collection_year" + } else { + NA_character_ + } + + collection_year <- if (!is.na(collection_year_col)) { + suppressWarnings( + as.integer(genome_data[[collection_year_col]]) + ) + } else { + integer(0) + } + + antibiotics <- if (!is.na(antibiotic_col)) { + collapse_unique(amr_data[[antibiotic_col]]) + } else { + NA_character_ + } + + unique_antibiotics <- if (!is.na(antibiotic_col)) { + x <- trimws(as.character(amr_data[[antibiotic_col]])) + length(unique(x[!is.na(x) & nzchar(x)])) + } else { + 0L + } + + drug_classes <- NA_character_ + + if (!is.na(antibiotic_col)) { + observed_drugs <- trimws(as.character(amr_data[[antibiotic_col]])) + observed_drugs <- unique( + observed_drugs[!is.na(observed_drugs) & nzchar(observed_drugs)] + ) + + drug_class_file <- file.path( + base_dir, + "data_raw", + "drug_class.tsv" + ) + + if (file.exists(drug_class_file)) { + drug_class_map <- utils::read.delim( + drug_class_file, + stringsAsFactors = FALSE, + check.names = FALSE + ) + + drug_classes <- collapse_unique( + drug_class_map$drug_class[ + drug_class_map$drug %in% observed_drugs + ] + ) + } + } + + median_genome_length <- if (!is.na(genome_length_col)) { + safe_median(genome_data[[genome_length_col]]) + } else { + NA_real_ + } + + median_gc_content <- if (!is.na(gc_content_col)) { + safe_median(genome_data[[gc_content_col]]) + } else { + NA_real_ + } + + median_cds <- if (!is.na(cds_col)) { + safe_median(genome_data[[cds_col]]) + } else { + NA_real_ + } + + median_checkm_completeness <- if (!is.na(checkm_complete_col)) { + safe_median(genome_data[[checkm_complete_col]]) + } else { + NA_real_ + } + + median_checkm_contamination <- if (!is.na(checkm_contam_col)) { + safe_median(genome_data[[checkm_contam_col]]) + } else { + NA_real_ + } + + summary_row <- tibble::tibble( + query = user_bac, + + total_genomes = as.integer(total_genomes), + wgs_genomes = as.integer(wgs_genomes), + complete_genomes = as.integer(complete_genomes), + + amr_genomes = as.integer( + length(intersect(genome_data$.genome_id, amr_ids)) + ), + amr_records = as.integer(nrow(amr_data)), + unique_antibiotics = as.integer(unique_antibiotics), + antibiotics = antibiotics, + drug_classes = drug_classes, + + checkm_available = as.integer(checkm_available), + qc_pass_genomes = as.integer(qc_pass), + qc_fail_genomes = as.integer(qc_fail), + + median_genome_length = median_genome_length, + median_gc_content = median_gc_content, + median_cds = median_cds, + median_checkm_completeness = median_checkm_completeness, + median_checkm_contamination = median_checkm_contamination, + collection_year_min = if ( + length(collection_year) && + any(!is.na(collection_year)) + ) { + min(collection_year, na.rm = TRUE) + } else { + NA_integer_ + }, + collection_year_max = if ( + length(collection_year) && + any(!is.na(collection_year)) + ) { + max(collection_year, na.rm = TRUE) + } else { + NA_integer_ + } + ) + + if (isTRUE(verbose)) { + message( + "Availability summary for '", user_bac, "': ", + total_genomes, " genomes; ", + wgs_genomes, " WGS; ", + complete_genomes, " Complete; ", + length(intersect(genome_data$.genome_id, amr_ids)), + " with AMR records; ", + qc_pass, " pass metadata QC." + ) + } + + summary_row +} + ######################### # Manifest helpers # ######################### @@ -509,7 +1050,10 @@ } # To distinguish multiple manifests in the same bug directory -.manifest_find_latest <- function(duckdb_path) { +.manifest_find_latest <- function( + duckdb_path, + require_success = TRUE +) { manifest_dir <- dirname(normalizePath( duckdb_path, mustWork = FALSE @@ -525,7 +1069,33 @@ return(NULL) } - manifests[which.max(file.info(manifests)$mtime)] + manifests <- manifests[ + order(file.info(manifests)$mtime, decreasing = TRUE) + ] + + if (!isTRUE(require_success)) { + return(manifests[[1]]) + } + + for (path in manifests) { + manifest <- tryCatch( + jsonlite::read_json(path, simplifyVector = FALSE), + error = function(e) NULL + ) + + if (is.null(manifest) || !length(manifest$runs)) { + next + } + + if (any(purrr::map_lgl( + manifest$runs, + ~ identical(.x$status, "success") + ))) { + return(path) + } + } + + NULL } #' Resume provenance logging in an existing manifest @@ -742,6 +1312,395 @@ ) } +#' Export a dyad-centric feature table +#' +#' Builds a one-row-per-dyad table linking protein-gene dyads to structural +#' features and HMMER annotations recorded in the dataset provenance manifest. +#' Feature values are deduplicated and combined into semicolon-separated +#' character fields. +#' +#' @param duckdb_path Character. Path to the source dataset DuckDB. Associated +#' Parquet files and provenance manifest are expected there. +#' @param output_path Character or NULL. Directory where the output Parquet file +#' will be written. Defaults to the DuckDB directory. +#' @param output_stem Character. Output filename stem. Default +#' `"dyad_annotations"`. +#' @param feature_scales Character vector of optional feature types to include. +#' If NULL, includes `struct` plus all HMMER databases recorded in the +#' manifest. +#' @param verbose Logical. Print progress messages. +#' +#' @return Invisibly returns the path to the generated Parquet file. +#' +#' @keywords internal +.exportDyadAnnotations <- function( + duckdb_path, + feature_scales = NULL, + verbose = TRUE +) { + duckdb_path <- normalizePath(duckdb_path, mustWork = TRUE) + parquet_dir <- dirname(duckdb_path) + + manifest_path <- .manifest_find_latest(duckdb_path) + + if (is.null(manifest_path)) { + stop( + "No provenance manifest found for: ", + duckdb_path + ) + } + + manifest <- jsonlite::read_json( + manifest_path, + simplifyVector = FALSE + ) + + # Find latest successful HMMER stage + hmmer_stage <- NULL + + for (run in rev(manifest$runs %||% list())) { + stages <- run$stages %||% list() + + matches <- purrr::keep( + stages, + ~ identical(.x$name, "hmmer") && + identical(.x$status, "success") + ) + + if (length(matches)) { + hmmer_stage <- matches[[1]] + break + } + } + + if (is.null(hmmer_stage)) { + stop( + "No successful HMMER stage found in manifest: ", + manifest_path + ) + } + + hmmer_databases <- unique(as.character( + unlist( + hmmer_stage$parameters$databases %||% character(), + use.names = FALSE + ) + )) + + allowed_features <- c("struct", hmmer_databases) + + if (is.null(feature_scales)) { + feature_scales <- allowed_features + } else { + feature_scales <- unique(as.character(feature_scales)) + + unknown_features <- setdiff( + feature_scales, + allowed_features + ) + + if (length(unknown_features)) { + stop( + "Unsupported feature scale(s): ", + paste(unknown_features, collapse = ", "), + ". Available features: ", + paste(allowed_features, collapse = ", ") + ) + } + } + + con <- DBI::dbConnect( + duckdb::duckdb(), + dbdir = ":memory:" + ) + + duckdb_temp_dir <- file.path( + parquet_dir, + ".duckdb_temp" + ) + + dir.create( + duckdb_temp_dir, + recursive = TRUE, + showWarnings = FALSE + ) + + DBI::dbExecute( + con, + sprintf( + "SET temp_directory=%s", + DBI::dbQuoteString( + con, + normalizePath( + duckdb_temp_dir, + winslash = "/", + mustWork = TRUE + ) + ) + ) + ) + + on.exit( + { + try(DBI::dbDisconnect(con, shutdown = TRUE), silent = TRUE) + unlink( + duckdb_temp_dir, + recursive = TRUE, + force = TRUE + ) + }, + add = TRUE + ) + + parquet_sql <- function(dataset_name) { + path <- file.path( + parquet_dir, + paste0(dataset_name, ".parquet") + ) + + if (!file.exists(path)) { + return(NULL) + } + + normalizePath( + path, + winslash = "/", + mustWork = TRUE + ) + } + + # Initialize using protein-gene dyads + genome_gene_protein_path <- parquet_sql( + "genome_gene_protein" + ) + + if (is.null(genome_gene_protein_path)) { + stop( + "Required Parquet file not found: ", + file.path( + parquet_dir, + "genome_gene_protein.parquet" + ) + ) + } + + DBI::dbExecute( + con, + sprintf( + " + CREATE OR REPLACE VIEW protein_gene AS + SELECT DISTINCT + protein_ids AS protein, + REPLACE(Gene, '~', '.') AS gene + FROM read_parquet('%s') + WHERE protein_ids IS NOT NULL + AND Gene IS NOT NULL + ", + genome_gene_protein_path + ) + ) + + DBI::dbExecute( + con, + " + CREATE OR REPLACE VIEW protein_gene_dyad AS + SELECT DISTINCT + protein, + gene, + CONCAT(protein, '|', gene) AS dyad + FROM protein_gene + " + ) + + # Finish initializing with one row per dyad + feature_select <- character() + feature_joins <- character() + + # Pangenome graph structural variant ('struct') annotations + if ("struct" %in% feature_scales) { + struct_path <- parquet_sql("struct") + + if (is.null(struct_path)) { + if (isTRUE(verbose)) { + message( + "Skipping struct: struct.parquet was not found." + ) + } + } else { + DBI::dbExecute( + con, + sprintf( + " + CREATE OR REPLACE VIEW struct_genes AS + SELECT DISTINCT + s.struct, + t.gene + FROM read_parquet('%s') s + CROSS JOIN UNNEST( + string_split(s.struct, '.') + ) AS t(gene) + WHERE s.value = 1 + ", + struct_path + ) + ) + + DBI::dbExecute( + con, + " + CREATE OR REPLACE VIEW dyad_struct AS + SELECT + pgd.dyad, + string_agg( + DISTINCT sg.struct, + ';' + ORDER BY sg.struct + ) AS struct + FROM protein_gene_dyad pgd + JOIN struct_genes sg + ON pgd.gene = sg.gene + GROUP BY pgd.dyad + " + ) + + feature_select <- c( + feature_select, + "ds.struct" + ) + + feature_joins <- c( + feature_joins, + "LEFT JOIN dyad_struct ds ON b.dyad = ds.dyad" + ) + } + } + + # HMMER feature annotations + for (database in intersect( + hmmer_databases, + feature_scales + )) { + dataset_name <- paste0( + "protein_", + database + ) + + hmmer_path <- parquet_sql(dataset_name) + + if (is.null(hmmer_path)) { + if (isTRUE(verbose)) { + message( + "Skipping ", + database, + ": ", + dataset_name, + ".parquet was not found." + ) + } + next + } + + view_name <- paste0( + "dyad_", + make.names(database) + ) + + DBI::dbExecute( + con, + sprintf( + " + CREATE OR REPLACE VIEW %s AS + SELECT + pgd.dyad, + string_agg( + DISTINCT h.query_name, + ';' + ORDER BY h.query_name + ) AS feature + FROM protein_gene_dyad pgd + JOIN read_parquet('%s') h + ON pgd.protein = h.protein + WHERE h.query_name IS NOT NULL + GROUP BY pgd.dyad + ", + view_name, + hmmer_path + ) + ) + + # Give AMRFinder a better human-readable name (ARG, in this case) + output_column <- if (identical(database, "AMRFinder")) { + "ARG" + } else { + database + } + + alias <- paste0( + "d_", + make.names(database) + ) + + feature_select <- c( + feature_select, + sprintf( + '%s.feature AS "%s"', + alias, + output_column + ) + ) + + feature_joins <- c( + feature_joins, + sprintf( + "LEFT JOIN %s %s ON b.dyad = %s.dyad", + view_name, + alias, + alias + ) + ) + } + + select_features <- if (length(feature_select)) { + paste0( + ",\n ", + paste(feature_select, collapse = ",\n ") + ) + } else { + "" + } + + join_features <- if (length(feature_joins)) { + paste0( + "\n ", + paste(feature_joins, collapse = "\n ") + ) + } else { + "" + } + + result_sql <- paste0( + " + SELECT + b.dyad, + b.protein, + b.gene", + select_features, + " + FROM protein_gene_dyad b", + join_features + ) + + if (isTRUE(verbose)) { + message("Building dyad annotation table.") + } + + DBI::dbGetQuery( + con, + result_sql + ) |> + tibble::as_tibble() +} + ######################### # HMMER helpers # ######################### diff --git a/man/checkDataAvailability.Rd b/man/checkDataAvailability.Rd new file mode 100644 index 0000000..67aa521 --- /dev/null +++ b/man/checkDataAvailability.Rd @@ -0,0 +1,71 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data_curation.R +\name{checkDataAvailability} +\alias{checkDataAvailability} +\title{Check BV-BRC data availability for selected bacterial taxa} +\usage{ +checkDataAvailability( + user_bacs, + base_dir = ".", + metadata_method = c("api", "cli"), + max_checkm_contam = 5, + min_checkm_complete = 95, + gc_deviations = NULL, + length_deviations = NULL, + cds_deviations = NULL, + verbose = TRUE +) +} +\arguments{ +\item{user_bacs}{Character vector. Mixed taxon IDs and/or species names.} + +\item{base_dir}{Character. Project root. Used by the legacy CLI path. +Default \code{"."}.} + +\item{metadata_method}{Character. Metadata backend: \code{"api"} (default) or +\code{"cli"}.} + +\item{max_checkm_contam}{Numeric. Maximum allowed CheckM contamination +(\%) for the QC summary. Default \code{5}.} + +\item{min_checkm_complete}{Numeric. Minimum allowed CheckM completeness +(\%) for the QC summary. Default \code{95}.} + +\item{gc_deviations}{Numeric. Maximum SDs from the median GC content. +Optional. Default \code{NULL}.} + +\item{length_deviations}{Numeric. Maximum SDs from the median genome +length. Optional. Default \code{NULL}.} + +\item{cds_deviations}{Numeric. Maximum SDs from the median CDS count. +Optional. Default \code{NULL}.} + +\item{verbose}{Logical. If TRUE, print progress messages. Default \code{TRUE}.} +} +\value{ +A tibble with one row per requested taxon containing summary +statistics describing genome availability, sequencing status, AMR data +availability, and the number of genomes passing metadata QC. +} +\description{ +Performs a lightweight metadata-only query to summarize the data available +for one or more bacterial taxa before running the genome-download and +feature-processing workflows. +} +\details{ +By default, genome IDs and metadata are resolved through the BV-BRC Data API. +The legacy CLI path can be selected with \code{metadata_method = "cli"}. +} +\examples{ +\dontrun{ +checkDataAvailability( + c("Staphylococcus argenteus", "Streptococcus suis") +) + +checkDataAvailability( + c("1280", "28901"), + metadata_method = "api" +) +} + +} diff --git a/man/dot-checkDataPerTaxon.Rd b/man/dot-checkDataPerTaxon.Rd new file mode 100644 index 0000000..6a6ff3c --- /dev/null +++ b/man/dot-checkDataPerTaxon.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.checkDataPerTaxon} +\alias{.checkDataPerTaxon} +\title{Check BV-BRC data availability for a single taxon} +\usage{ +.checkDataPerTaxon( + user_bac, + base_dir = ".", + metadata_method = c("api", "cli"), + max_checkm_contam = 5, + min_checkm_complete = 95, + gc_deviations = NULL, + length_deviations = NULL, + cds_deviations = NULL, + verbose = TRUE +) +} +\arguments{ +\item{user_bac}{Character scalar. Taxon ID or species name.} + +\item{base_dir}{Character. Project root. Used by the legacy CLI path. +Default \code{"."}.} + +\item{metadata_method}{Character. Metadata backend: \code{"api"} (default) or +\code{"cli"}.} + +\item{max_checkm_contam}{Numeric. Maximum allowed CheckM contamination +(\%) for the QC summary. Default \code{5}.} + +\item{min_checkm_complete}{Numeric. Minimum allowed CheckM completeness +(\%) for the QC summary. Default \code{95}.} + +\item{gc_deviations}{Numeric. Maximum SDs from the median GC content. +Optional. Default \code{NULL}.} + +\item{length_deviations}{Numeric. Maximum SDs from the median genome +length. Optional. Default \code{NULL}.} + +\item{cds_deviations}{Numeric. Maximum SDs from the median CDS count. +Optional. Default \code{NULL}.} + +\item{verbose}{Logical. If TRUE, print progress messages. Default \code{TRUE}.} +} +\value{ +A one-row tibble containing genome and AMR availability statistics. +} +\description{ +Internal worker used by \code{checkDataAvailability()}. Each call resolves and +summarizes a taxon independently. +} +\keyword{internal} diff --git a/man/dot-exportDyadAnnotations.Rd b/man/dot-exportDyadAnnotations.Rd new file mode 100644 index 0000000..61e9597 --- /dev/null +++ b/man/dot-exportDyadAnnotations.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.exportDyadAnnotations} +\alias{.exportDyadAnnotations} +\title{Export a dyad-centric feature table} +\usage{ +.exportDyadAnnotations(duckdb_path, feature_scales = NULL, verbose = TRUE) +} +\arguments{ +\item{duckdb_path}{Character. Path to the source dataset DuckDB. Associated +Parquet files and provenance manifest are expected there.} + +\item{feature_scales}{Character vector of optional feature types to include. +If NULL, includes \code{struct} plus all HMMER databases recorded in the +manifest.} + +\item{verbose}{Logical. Print progress messages.} + +\item{output_path}{Character or NULL. Directory where the output Parquet file +will be written. Defaults to the DuckDB directory.} + +\item{output_stem}{Character. Output filename stem. Default +\code{"dyad_annotations"}.} +} +\value{ +Invisibly returns the path to the generated Parquet file. +} +\description{ +Builds a one-row-per-dyad table linking protein-gene dyads to structural +features and HMMER annotations recorded in the dataset provenance manifest. +Feature values are deduplicated and combined into semicolon-separated +character fields. +} +\keyword{internal} diff --git a/man/dot-filterGenomes.Rd b/man/dot-filterGenomes.Rd index c66b32a..64ccf30 100644 --- a/man/dot-filterGenomes.Rd +++ b/man/dot-filterGenomes.Rd @@ -9,7 +9,7 @@ base_dir = ".", evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"), verbose = TRUE, - fallback_to_bvbrc_cache = TRUE + fallback_to_bvbrc_cache = FALSE ) } \arguments{ diff --git a/man/dot-prepareHmmerDatabases.Rd b/man/dot-prepareHmmerDatabases.Rd index 9f8f0d1..8eed4f3 100644 --- a/man/dot-prepareHmmerDatabases.Rd +++ b/man/dot-prepareHmmerDatabases.Rd @@ -13,13 +13,21 @@ ) } \arguments{ -\item{hmmer_db_dir}{Directory to store HMMER databases} +\item{hmmer_db_dir}{Character. Directory where HMMER databases are cached.} -\item{databases}{List of databases to prepare (default: c("Pfam", "COG", "AMRFinder"))} +\item{databases}{Character vector of database names to prepare. Supported +built-in databases are \code{Pfam}, \code{COG}, and \code{AMRFinder}; \code{DefenseCas} is +prepared separately by the DefenseFinder/CasFinder workflow.} -\item{docker_image}{Docker image containing HMMER (default: "staphb/hmmer")} +\item{docker_image}{Character. Docker image containing HMMER tools used to +press the prepared databases. Default: \code{"staphb/hmmer"}.} -\item{hmmer_db_url}{If the databases contain custom database(s), the url is required to download the database.} +\item{hmmer_db_url}{NON-FUNCTIONAL. Character or \code{NULL}. URL used to download +a custom HMMER database when \code{databases} contains names not covered by the +built-in database definitions. This function is not currently active!} + +\item{verbose}{Logical. Print status messages while checking, downloading, +combining, and pressing databases. Default: \code{TRUE}.} } \value{ A list of paths to the database hmm files. diff --git a/man/dot-proteinAnnotations2Duckdb.Rd b/man/dot-proteinAnnotations2Duckdb.Rd index cfffa09..1d66a0a 100644 --- a/man/dot-proteinAnnotations2Duckdb.Rd +++ b/man/dot-proteinAnnotations2Duckdb.Rd @@ -11,12 +11,15 @@ ) } \arguments{ -\item{duckdb_path}{Path to the per-selection DuckDB database containing a -\code{protein_count} table (created by \code{\link[=CDHIT2duckdb]{CDHIT2duckdb()}}).} +\item{duckdb_path}{Character. Path to the per-selection DuckDB database +containing the \code{protein_count} table created by \code{\link[=CDHIT2duckdb]{CDHIT2duckdb()}}.} -\item{annotated_parquet}{Path to the combined HMMER results Parquet file -(e.g. \code{"results/Ecoli/protein_COG.parquet"}). The filename stem is used as -the table name in DuckDB.} +\item{databases}{Character vector of HMMER database names to process. +Each database must correspond to a \verb{protein_} annotation table +already present in the DuckDB.} + +\item{output_path}{Character. Directory where the genome-by-annotation +Parquet files will be written. Defaults to \code{dirname(duckdb_path)}.} } \value{ Invisibly returns the path to the written count Parquet file. diff --git a/man/dot-runHMMER.Rd b/man/dot-runHMMER.Rd index 72618a7..d693531 100644 --- a/man/dot-runHMMER.Rd +++ b/man/dot-runHMMER.Rd @@ -17,21 +17,30 @@ ) } \arguments{ -\item{duckdb_path}{path to the duckdb with protein sequences and list} +\item{duckdb_path}{Character. Path to the DuckDB database containing +\code{protein_cluster_seq}, which provides the protein sequences to analyze.} -\item{output_path}{path where HMMER output will be saved} +\item{output_path}{Character. Directory for HMMER intermediate and final +Parquet outputs.} -\item{threads}{number of threads} +\item{threads}{Integer. Total CPU budget used by HMMER jobs. Default: \code{8}.} -\item{hmmer_db_dir}{path to the directory where HMM databases are/will be downloaded} +\item{hmmer_db_dir}{Character. Directory containing the prepared HMMER +databases. If \code{NULL}, the default \code{amRdata} HMMER database cache is used.} -\item{databases}{list of HMM databases} +\item{databases}{Character vector of HMMER databases to run.} -\item{docker_image}{the docker image of HMMER} +\item{docker_image}{Character. Docker image containing HMMER. Default: +\code{"staphb/hmmer"}.} -\item{num_of_splits}{The number of splits of the protein sequence file for parallel processing} +\item{num_of_splits}{Integer. Number of chunks into which the protein +sequences should be divided. Must be a positive integer. The requested +value is automatically reduced when fewer protein sequences are available. +Default: \code{8}.} -\item{n_workers}{The number of parallel runs} +\item{n_workers}{Integer. Number of parallel HMMER jobs to run. Default: \code{8}.} + +\item{verbose}{Logical. Print progress messages. Default: \code{TRUE}.#'} } \description{ Wrapper for preparing HMM databases and running HMMER on protein sequences from duckdb and writing them. diff --git a/man/dot-runHmmerJob.Rd b/man/dot-runHmmerJob.Rd index 121d513..13ac6e4 100644 --- a/man/dot-runHmmerJob.Rd +++ b/man/dot-runHmmerJob.Rd @@ -8,7 +8,7 @@ JOB_NAME, FASTA, DB, - Total_proteins, + total_proteins, output_path = NULL, db_paths, docker_image = "staphb/hmmer", @@ -18,23 +18,32 @@ ) } \arguments{ -\item{JOB_NAME}{protein_chunk id} +\item{JOB_NAME}{Character. Identifier used for the HMMER job and output +filename.} -\item{FASTA}{fasta sequences} +\item{FASTA}{Character. File name of the protein FASTA chunk to search.} -\item{DB}{HMM database} +\item{DB}{Character. Name of the HMMER database to search against.} -\item{Total_proteins}{protein sequence count} +\item{total_proteins}{Integer. Total number of proteins in the full input +dataset, used to set HMMER's \code{-Z} and \code{--domZ} values.} -\item{output_path}{path for saving hmmer outputs} +\item{output_path}{Character. Directory containing the FASTA input, HMMER +output, and final Parquet result.} -\item{db_paths}{path to HMM database} +\item{db_paths}{List. Prepared HMMER database metadata indexed by database +name.} -\item{docker_image}{hmmer docker image (ideally from dockerhub)} +\item{docker_image}{Character. Docker image containing HMMER. Default: +\code{"staphb/hmmer"}.} -\item{threads}{number of threads} +\item{threads}{Integer. Total CPU budget used when calculating the number +of threads allocated to this job. Default: \code{8}.} -\item{n_workers}{number of parallel workers} +\item{n_workers}{Integer. Number of HMMER jobs being run in parallel. +Used to divide the CPU budget among jobs. Default: \code{8}.} + +\item{verbose}{Logical. Print progress messages. Default: \code{TRUE}.} } \value{ the filename of the parquet file with hmmer output post parsing diff --git a/man/exportProcessedData.Rd b/man/exportProcessedData.Rd index 07aff2a..f5dec37 100644 --- a/man/exportProcessedData.Rd +++ b/man/exportProcessedData.Rd @@ -10,6 +10,7 @@ exportProcessedData( amr_phenotype_mode = c("separate", "append"), export_formats = c("csv"), export_sequences = FALSE, + export_dyads = FALSE, tables = NULL, export_tables = TRUE, verbose = TRUE @@ -28,8 +29,20 @@ file.path(dirname(duckdb_path), "processed_exports").} \item{export_formats}{Character vector. Any of "csv", "tsv", "parquet", "xlsx".} -\item{tables}{Character vector or NULL. Tables to export. If NULL, exports the -standard processed tables present in the database.} +\item{export_sequences}{Logical. If TRUE, also exports gene and protein +sequence tables and the genome-to-gene-to-protein mapping. Default FALSE.} + +\item{export_dyads}{Logical. If TRUE, exports the optional dyad annotation +table. Each row represents a protein-gene dyad with semicolon-separated +mapped feature values for structural annotations and the HMMER databases +recorded in the dataset manifest. The export uses \code{export_formats}. +Default FALSE.} + +\item{tables}{Character vector or NULL. Tables to export. If NULL, exports +the standard processed tables plus HMMER tables recorded in the manifest.} + +\item{export_tables}{Logical. If TRUE, write the selected tables to disk. +Default TRUE.} \item{verbose}{Logical. If TRUE, prints progress messages.} } diff --git a/man/runDataProcessing.Rd b/man/runDataProcessing.Rd index 4a80dc7..81c327d 100644 --- a/man/runDataProcessing.Rd +++ b/man/runDataProcessing.Rd @@ -8,6 +8,7 @@ runDataProcessing( duckdb_path, output_path = NULL, threads = 8, + export_tabular_data = FALSE, panaroo_split_jobs = FALSE, panaroo_core_threshold = 0.9, panaroo_len_dif_percent = 0.95, @@ -44,6 +45,10 @@ CD-HIT, HMMER, and final Parquet outputs. If \code{NULL}, defaults to \item{threads}{Integer. Shared concurrency budget used across Panaroo, CD-HIT, and HMMER. Defaults to \code{8}.} +\item{export_tabular_data}{Logical. If TRUE, automatically call +\code{\link[=exportProcessedData]{exportProcessedData()}} after the processing pipeline completes to create +the default human-readable exports. Default: \code{FALSE}.} + \item{panaroo_split_jobs}{Logical. If \code{TRUE}, Panaroo runs in multiple batches that can be merged by \code{\link[=.mergePanaroo]{.mergePanaroo()}}. If \code{FALSE}, Panaroo runs once on all isolates. Default: \code{FALSE}.} @@ -153,6 +158,10 @@ upstream curation workflow. over the generated Parquet files. \item Records processing parameters, software versions, database selections, and other provenance information in the dataset manifest. +\item Creates the protein-gene dyad feature map for downstream graph analysis. +\item If \code{export_tabular_data = TRUE}, exports features in human-readable CSVs. +Additional export options are available by calling \code{exportProcessedData()} +after \code{runDataProcessing()} completes. } \strong{Threading} @@ -188,6 +197,9 @@ stage parameters. \item \strong{Metadata cleaning + Parquet export} via \code{\link[=cleanData]{cleanData()}} -> writes Parquet files to \code{output_path}, and builds a \strong{Parquet-backed DuckDB} (\verb{*_parquet.duckdb}) with views over those Parquets. + +\item \strong{Dyad feature mapping} via \code{\link[=buildDyadFeatureMap]{buildDyadFeatureMap()}} -> writes the +protein-gene dyad feature network used for downstream graph analysis. } }