From f708e7a81c403ecc974dbf46e77377bf3b9bfc10 Mon Sep 17 00:00:00 2001 From: Evan Pierce Brenner <108823789+epbrenner@users.noreply.github.com> Date: Wed, 12 Aug 2026 16:05:01 -0600 Subject: [PATCH] Major workflow refactor 1. Moved many common and utility functions to helpers.R 2. Added new manifest feature, logging all steps of workflow - Manifest is also used to identify which HMMER databases have run on a bug. - Adds extensive data provenance logging + moves databases to R user cache. 3. Removes InterProScan entirely 4. Adds Abhirupa's HMMER functionality to data_processing.R, removed runHMMER.R 5. Replaced "domain" tables/features with Pfam/COG/AMRFinder/DefenseCas 6. Added jsonlite to Imports 7. stuff Co-Authored-By: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> --- DESCRIPTION | 1 + R/data_curation.R | 254 +++-- R/data_processing.R | 2409 +++++++++++++++++++++++++++++++------------ R/helpers.R | 939 +++++++++++++++++ R/runHMMER.R | 1352 ------------------------ 5 files changed, 2867 insertions(+), 2088 deletions(-) create mode 100644 R/helpers.R delete mode 100644 R/runHMMER.R diff --git a/DESCRIPTION b/DESCRIPTION index 64266d6..ee2f1c5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -42,6 +42,7 @@ Imports: glue, grid, gridExtra, + jsonlite, knitr, purrr, readr, diff --git a/R/data_curation.R b/R/data_curation.R index 6df70cc..d5a5833 100644 --- a/R/data_curation.R +++ b/R/data_curation.R @@ -1,34 +1,3 @@ -#' Helps ensure trailing 0s are retained in genome IDs for proper downloading -#' @keywords internal -.id_checker <- function(x) { - # Taxon IDs are just numbers, genome IDs have decimals, this tells them apart - grepl("^[0-9]+$", x) -} - -#' A helper used in data_curation.R and data_processing.R to ensure exported tables -#' don't lose trailing zeroes. Should be relocated into a common helpers/utilities -#' script later. -#' @keywords internal -.preserve_export_id_text <- function(df) { - df <- tibble::as_tibble(df) - - id_pattern <- paste0( - "(^|[._])(", - "genome(_drug)?_id|taxon_id|", - "assembly_accession|bioproject_accession|biosample_accession|", - "refseq_accessions?|genbank_accessions?|sra_accession|pmid|", - "gene_id|protein_id|domain_id|cluster_id|AccNum|id", - ")$" - ) - - id_cols <- names(df)[grepl(id_pattern, names(df), ignore.case = TRUE)] - if (length(id_cols)) { - df[id_cols] <- lapply(df[id_cols], as.character) - } - - df -} - #' Helps tag genomes with their AMR evidence for parsing #' @keywords internal .create_amr_tagged_view <- function(con) { @@ -1559,53 +1528,6 @@ retrieveMetadata <- function(user_bacs, list(duckdbConnection = con, table_name = "metadata") } -# FASTA sanitizer to ensure Panaroo compatibility with BV-BRC CLI downloads -.strip_fasta_preamble <- function(fna_path) { - if (!file.exists(fna_path)) { - return(invisible(FALSE)) - } - txt <- readLines(fna_path, warn = FALSE) - first <- which(grepl("^\\s*>", txt))[1] - if (is.na(first)) { - return(invisible(FALSE)) - } - if (first > 1L) { - txt <- txt[first:length(txt)] - txt[1] <- sub("^\\ufeff", "", txt[1]) - writeLines(txt, fna_path, sep = "\n", useBytes = TRUE) - return(invisible(TRUE)) - } - invisible(FALSE) -} - -# GFF sanitizer to ensure Panaroo compatibility with BV-BRC CLI downloads -.sanitize_gff <- function(gff_path) { - if (!file.exists(gff_path)) { - return(invisible(FALSE)) - } - lines <- readLines(gff_path, warn = FALSE) - if (length(lines) == 0L) { - return(invisible(FALSE)) - } - if (!grepl("^##gff-version\\s*3", lines[1])) { - lines <- c("##gff-version 3", lines) - } - out <- purrr::map_chr(lines, function(line) { - if (grepl("^#", line)) { - return(line) - } - parts <- strsplit(line, "[\t ]", perl = TRUE)[[1]] - if (length(parts) >= 9) { - paste(c(parts[1:8], paste(parts[9:length(parts)], collapse = " ")), collapse = "\t") - } else { - line - } - }) - writeLines(out, gff_path, sep = "\n", useBytes = TRUE) - invisible(TRUE) -} - - #' Filter genomes by AMR phenotype and metadata, and store results in DuckDB #' #' Preferred path: use per-selection DB "metadata" table (from retrieveMetadata()) and @@ -1819,24 +1741,6 @@ retrieveMetadata <- function(user_bacs, } ### BV-BRC CLI downloader [slower by comparison, but does not need FTP server] - -#' Helps normalize Docker paths -#' @keywords internal -.docker_path <- function(p) gsub("\\\\", "/", normalizePath(p, mustWork = FALSE)) - -#' Helps run a shell inside a container, and prefers bash (don't we all?) -#' @keywords internal -.pick_shell <- function(image) { - chk <- suppressWarnings(system2("docker", - c( - "run", "--rm", image, "sh", "-lc", - "command -v bash >/dev/null || echo NOBASH" - ), - stdout = TRUE, stderr = TRUE - )) - if (length(chk) && any(grepl("NOBASH", chk))) "sh" else "bash" -} - #' Using p3-dump-genomes in CLI to fetch FASTA and .gto files #' @keywords internal .cli_dump_fastas_gto_chunk <- function(image, out_dir, genome_ids, tag, tries = 3L) { @@ -2279,9 +2183,99 @@ prepareGenomes <- function(user_bacs, evidence_mode <- match.arg(evidence_mode) base_dir <- normalizePath(base_dir, mustWork = FALSE) - .ensure_bvbrc_cache(base_dir = base_dir, verbose = verbose) + paths <- .buildDBpath( + base_dir = base_dir, + user_bacs = user_bacs + ) + + manifest_path <- file.path( + dirname(paths$db_path), + paste0("manifest_", .manifest_run_id(), ".json") + ) + + manifest <- .manifest_start( + manifest_path = manifest_path, + dataset_id = .generateDBname(user_bacs), + duckdb_path = normalizePath(paths$db_path, mustWork = FALSE), + base_dir = base_dir, + selection = list( + user_bacs = as.character(user_bacs), + genome_id_file = if (is.null(genome_id_file)) { + NULL + } else { + normalizePath(genome_id_file, mustWork = FALSE) + } + ), + hash_files = FALSE + ) + + manifest <- .manifest_event( + manifest, + message = "Started genome curation run.", + details = list( + method = method, + evidence_mode = evidence_mode, + overwrite = overwrite + ) + ) + + run_failed <- TRUE + on.exit( + if (run_failed) { + .manifest_finish( + manifest, + status = "failed", + error = "prepareGenomes() exited before successful completion." + ) + }, + add = TRUE + ) + + manifest <- .manifest_stage( + manifest, + name = "prepare_bvbrc_cache", + status = "success", + parameters = list( + max_age_days = 30L + ), + outputs = file.path(base_dir, "data", "bvbrc", "bvbrcData.duckdb"), + tool = list( + name = "BV-BRC", + interface = "p3-all-genomes" + ) + ) + + .ensure_bvbrc_cache( + base_dir = base_dir, + verbose = verbose + ) + + if (isTRUE(verbose)) { + message("Step 0: Building AMR metadata (retrieveMetadata)") + } + + manifest <- .manifest_stage( + manifest, + name = "retrieve_metadata", + status = "running", + parameters = list( + filter_type = "AMR", + abx = "All", + 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, + overwrite = overwrite + ), + inputs = if (!is.null(genome_id_file)) genome_id_file else character(), + tool = list( + name = "BV-BRC", + docker_image = "danylmb/bvbrc:5.3" + ) + ) - if (isTRUE(verbose)) message("Step 0: Building AMR metadata (retrieveMetadata)") invisible(retrieveMetadata( user_bacs = user_bacs, genome_id_file = genome_id_file, @@ -2298,6 +2292,28 @@ prepareGenomes <- function(user_bacs, verbose = verbose )) + manifest <- .manifest_stage( + manifest, + name = "retrieve_metadata", + status = "success", + parameters = list( + filter_type = "AMR", + abx = "All", + 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, + overwrite = overwrite + ), + outputs = normalizePath(paths$db_path, mustWork = FALSE), + tool = list( + name = "BV-BRC", + docker_image = "danylmb/bvbrc:5.3" + ) + ) + if (isTRUE(verbose)) message("Step 1: Filtering genomes for download by evidence: ", evidence_mode) f_out <- .filterGenomes( base_dir = base_dir, @@ -2342,6 +2358,22 @@ prepareGenomes <- function(user_bacs, return(NULL) } + manifest <- .manifest_stage( + manifest, + name = "download_genomes", + status = "success", + parameters = list( + method = method, + workers = num_workers, + evidence_mode = evidence_mode, + overwrite = overwrite + ), + outputs = file.path(paths$db_dir, "genomes"), + metrics = list( + genomes_downloaded = length(ids) + ) + ) + if (isTRUE(verbose)) message("Step 3: Formatting data into a database for further processing") out <- genomeList( base_dir = base_dir, @@ -2355,8 +2387,27 @@ prepareGenomes <- function(user_bacs, message("") message("Continue with downstream processing using:") message('runDataProcessing("', normalizePath(paths$db_path), '")') + message("") + message("Provenance manifest saved to:") + message(" ", normalizePath(manifest$path)) } + manifest <- .manifest_stage( + manifest, + name = "build_genome_file_table", + status = "success", + outputs = c( + paths$db_path, + file.path( + paths$db_dir, + paste0(.generateDBname(user_bacs), ".txt") + ) + ), + metrics = list( + genomes = length(ids) + ) + ) + export_res <- NULL if (isTRUE(export_tables) || isTRUE(load_tables)) { export_res <- exportTables( @@ -2375,6 +2426,13 @@ prepareGenomes <- function(user_bacs, )) } + run_failed <- FALSE + + .manifest_finish( + manifest, + status = "success" + ) + invisible(out) } diff --git a/R/data_processing.R b/R/data_processing.R index 31190f2..ae21216 100644 --- a/R/data_processing.R +++ b/R/data_processing.R @@ -1,32 +1,6 @@ #' @importFrom data.table := NULL -#' Normalize a host filesystem path for use in Docker -#' -#' Converts Windows and mixed-separator paths to forward slashes -#' and applies `normalizePath()` without requiring the path to exist. -#' -#' @param p Character scalar. A filesystem path on the host OS. -#' -#' @return A normalized path string. -#' -#' @keywords internal -.docker_path <- function(p) gsub("\\\\", "/", normalizePath(p, mustWork = FALSE)) - -# Map host paths under mounted root to container path -#' .to_container() -#' -#' Used for OS-agnostic mapping of Docker directories and mount paths -#' -#' @keywords internal -#' @examples NULL -.to_container <- function(x, host_root, container_root = "/work") { - host_root_unix <- .docker_path(host_root) - x_unix <- .docker_path(x) - pattern <- paste0("^", gsub("([\\^$.|?*+(){}\\[\\]\\\\])", "\\\\\\\\\\1", host_root_unix)) - sub(pattern, container_root, x_unix) -} - # Launch Panaroo to build a pangenome (per batch) #' processPanaroo() #' @@ -139,126 +113,7 @@ NULL invisible(res) } -#' Remove pseudogene annotations from Panaroo input GFF files -#' -#' Cleans GFF annotation files of `pseudogene` feature records only. Cleaned -#' GFFs are written to a subdirectory under `output_path` and swapped into the -#' Panaroo input list, leaving the original genome annotations alone. -#' -#' This optional preprocessing step can reduce weird runtime stalls during -#' Panaroo graph construction for some BV-BRC/PATRIC genome annotations that -#' contain troublesome pseudogene features. -#' -#' @param panaroo_input_files Character vector of `"gff fna"` input lines used -#' by Panaroo. -#' @param output_path Character scalar. Base directory for temporary cleaned -#' GFF files and audit outputs. -#' @param clean_dir Character scalar. Name of the subdirectory created beneath -#' `output_path` to store cleaned GFF files. Default `"gff_clean"`. -#' -#' @return A list containing: -#' \itemize{ -#' \item `panaroo_input_files` — rewritten Panaroo input lines pointing to the -#' cleaned GFF files. -#' \item `audit` — a tibble summarizing, for each genome, the total number of -#' annotated features, the number of pseudogenes removed, and the number of -#' remaining features. -#' } -#' -#' @details -#' This performs lightweight preprocessing only, removing feature records whose -#' third GFF column is exactly `"pseudogene"` and does not otherwise modify -#' annotation coordinates, attributes, or sequence files. FASTA paths are unchanged. -#' -#' @keywords internal -.stripPseudogeneGFFs <- function(panaroo_input_files, - output_path, - clean_dir = "gff_clean") { - # Normalize our paths - panaroo_input_files <- as.character(panaroo_input_files) - output_path <- .docker_path(output_path) - - # Set the directory to place cleaned GFFs into - clean_root <- file.path(output_path, clean_dir) - dir.create(clean_root, recursive = TRUE, showWarnings = FALSE) - - # Where the cleaned up Panaroo input and clean audit is stored - out_lines <- character(length(panaroo_input_files)) - audit <- vector("list", length(panaroo_input_files)) - - for (i in seq_along(panaroo_input_files)) { - - # Read the Panaroo gff + fna input lines - line <- panaroo_input_files[[i]] - parts <- strsplit(line, "\\s+")[[1]] - - # If you're missing either a gff or an fna file in there somehow - if (length(parts) < 2L) { - stop("Broken Panaroo input line: ", line) - } - - # Read in the files parsed above - gff_in <- .docker_path(parts[1]) - fna_in <- .docker_path(parts[2]) - - # If it didn't read in - if (!file.exists(gff_in)) { - stop("Missing GFF file: ", gff_in) - } - if (!file.exists(fna_in)) { - stop("Missing FNA file: ", fna_in) - } - - # What we're saving out - gff_out <- file.path(clean_root, basename(gff_in)) - - # Read in the GFF lines and fine the comment lined headers - gff_lines <- readLines(gff_in, warn = FALSE) - is_header <- startsWith(gff_lines, "#") - body <- gff_lines[!is_header] - - # If there's nothing in there to parse - if (length(body) == 0L) { - writeLines(gff_lines, gff_out, useBytes = TRUE) - n_total <- 0L - n_pseudogene <- 0L - n_kept <- 0L - } else { - # Otherwise, find the pseudogene lines and save everything but those - # Now with 100% more purrr - fields <- strsplit(body, "\t", fixed = TRUE) - types <- purrr::map_chr( - fields, - \(x) if (length(x) >= 3L) x[[3]] else NA_character_ - ) - keep <- !is.na(types) & types != "pseudogene" - - cleaned <- c(gff_lines[is_header], body[keep]) - writeLines(cleaned, gff_out, useBytes = TRUE) - - # We love stats - n_total <- length(body) - n_pseudogene <- sum(!keep, na.rm = TRUE) - n_kept <- sum(keep, na.rm = TRUE) - } - - # Record what we did and save it into the audit log - audit[[i]] <- tibble::tibble( - gff_in = gff_in, - gff_out = gff_out, - n_total_features = n_total, - n_pseudogene = n_pseudogene, - n_kept = n_kept - ) - - out_lines[[i]] <- paste(gff_out, fna_in) - } - list( - panaroo_input_files = out_lines, - audit = dplyr::bind_rows(audit) - ) -} #' Run Panaroo for Pangenome Analysis in Parallel Batches @@ -802,7 +657,7 @@ NULL #' * `.runPanaroo()` — core Panaroo execution #' * `.mergePanaroo()` — merge multiple Panaroo batches #' * `.panaroo2duckdb()` — import Panaroo results into DuckDB -#' * [runDataProcessing()] — full pipeline including CD-HIT & InterProScan +#' * [runDataProcessing()] — full pipeline including CD-HIT & HMMER #' #' @examples #' \dontrun{ @@ -1105,436 +960,1293 @@ CDHIT2duckdb <- function(duckdb_path, invisible(TRUE) } - -#' Check or install InterProScan data bundle -#' -#' Ensures that the InterProScan data directory exists locally, downloading -#' and verifying the appropriate tarball when necessary. +#' Download and prepare HMMER databases for generating new file types. #' -#' @param version InterProScan version string. -#' @param dest_dir Directory where data should be installed. -#' @param docker_image Docker image string for InterProScan. -#' @param platform Character indicating Docker platform (e.g. `"linux/amd64"`). -#' @param curl_bin Path to curl executable. -#' @param verbose Logical; print status messages. +#' @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. #' -#' @return A list containing `data_dir` and `ready` status. +#' @returns A list of paths to the database hmm files. #' #' @keywords internal -.checkInterProData <- function( - version = "5.76-107.0", - dest_dir = "inst/extdata/interpro", - docker_image = sprintf("interpro/interproscan:%s", version), - platform = "linux/amd64", - curl_bin = "curl", - verbose = TRUE +#' @examples +.prepareHmmerDatabases <- function( + hmmer_db_dir, + databases = c("Pfam", "COG", "AMRFinder"), + docker_image = "staphb/hmmer", + hmmer_db_url = NULL, + verbose = TRUE ) { - msg <- function(...) if (verbose) message(sprintf(...)) - if (!dir.exists(dest_dir)) dir.create(dest_dir, recursive = TRUE, showWarnings = FALSE) - dest_dir <- normalizePath(dest_dir, mustWork = TRUE) + hmmer_db_dir <- normalizePath( + hmmer_db_dir, + mustWork = FALSE + ) + + dir.create( + hmmer_db_dir, + recursive = TRUE, + showWarnings = FALSE + ) - root_dir <- file.path(dest_dir, sprintf("interproscan-%s", version)) - data_dir <- file.path(root_dir, "data") + options(timeout = max(3600, getOption("timeout"))) - # Simple existence check - if (dir.exists(data_dir) && length(list.files(data_dir, recursive = TRUE)) > 0) { - msg("InterProScan data already present at: %s", data_dir) - return(list(data_dir = normalizePath(data_dir), ready = TRUE)) - } + dbs <- list( + Pfam = list( + dir = file.path(hmmer_db_dir, "Pfam"), + hmm_name = "Pfam-A.hmm", + url = "https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz", + type = "gz" + ), - # Download bundle if needed - tar_url <- sprintf( - "http://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/%s/alt/interproscan-data-%s.tar.gz", - version, version + COG = list( + dir = file.path(hmmer_db_dir, "COG"), + hmm_name = "COG_database2024.hmm", + url = "http://boabio.belozersky.msu.ru/media/COG_database2024.zip", + type = "zip" + ), + + AMRFinder = list( + dir = file.path(hmmer_db_dir, "AMRFinder"), + hmm_name = NULL, + url = "https://ftp.ncbi.nlm.nih.gov/hmm/NCBIfam-AMRFinder/latest/NCBIfam-AMRFinder.HMM.tar.gz", + type = "tar.gz" + ) ) - md5_url <- paste0(tar_url, ".md5") - tar_path <- file.path(dest_dir, basename(tar_url)) - md5_path <- paste0(tar_path, ".md5") - if (!file.exists(tar_path)) { - msg("Downloading InterProScan data bundle.") - status_tar <- system2(curl_bin, c("-L", "-o", tar_path, tar_url)) - status_md5 <- system2(curl_bin, c("-L", "-o", md5_path, md5_url)) - if (status_tar != 0 || status_md5 != 0) { - stop("Failed to download InterProScan data bundle.") + # Add custom database(s) + missing_dbs <- setdiff(databases, names(dbs)) + + if (length(missing_dbs) > 0) { + + if (is.null(hmmer_db_url)) { + stop( + "hmmer_db_url must be supplied when using custom databases" + ) + } + + get_db_type <- function(url) { + + file <- basename(url) + + if (grepl("\\.(tar\\.gz|tgz)$", file, ignore.case = TRUE)) { + return("tar.gz") + } else if (grepl("\\.zip$", file, ignore.case = TRUE)) { + return("zip") + } else if (grepl("\\.gz$", file, ignore.case = TRUE)) { + return("gz") + } else { + stop( + "Unsupported archive type: ", + file + ) + } } - } - msg("Verifying MD5 checksum.") - md5_expected <- sub("\\s+.*$", "", readLines(md5_path)[1]) - md5_actual <- tools::md5sum(tar_path)[[1]] - if (!identical(tolower(md5_expected), tolower(md5_actual))) { - stop("MD5 checksum mismatch for InterProScan data bundle.") + for (db_name in missing_dbs) { + dbs[[db_name]] <- list( + dir = file.path(hmmer_db_dir, db_name), + hmm_name = NULL, + url = hmmer_db_url, + type = get_db_type(hmmer_db_url) + ) + } } - msg("Extracting InterProScan data bundle.") - utils::untar(tar_path, exdir = dest_dir, tar = "internal") + dbs <- dbs[databases] + db_paths <- list() - msg("Data unpacked successfully.") - return(list(data_dir = normalizePath(data_dir), ready = TRUE)) -} + for (db_name in names(dbs)) { + db <- dbs[[db_name]] -#' Internal helpers for reading InterProScan TSV outputs -#' -#' Provide standardized column names, types, and a reader wrapper for the -#' InterProScan tab-delimited output format. -#' -#' @param filepath Path to a `.tsv` or `.tsv.gz` InterProScan result file. -#' -#' @return A tibble of parsed InterProScan output. -#' -#' @keywords internal -.getDfIPRColNames <- function() { - c( - "AccNum", "SeqMD5Digest", "SLength", "Analysis", - "DB.ID", "SignDesc", "StartLoc", "StopLoc", "Score", - "Status", "RunDate", "IPRAcc", "IPRDesc", "placeholder" - ) -} + dir.create( + db$dir, + recursive = TRUE, + showWarnings = FALSE + ) -#' Internal helpers for reading InterProScan TSV outputs -#' -#' Provide standardized column names, types, and a reader wrapper for the -#' InterProScan tab-delimited output format. -#' -#' @param filepath Path to a `.tsv` or `.tsv.gz` InterProScan result file. -#' -#' @return A tibble of parsed InterProScan output. -#' -#' @keywords internal + if (verbose) { + message("Checking ", db_name) + } -.getDfIPRColTypes <- function() { - readr::cols( - "AccNum" = readr::col_character(), - "SeqMD5Digest" = readr::col_character(), - "SLength" = readr::col_integer(), - "Analysis" = readr::col_character(), - "DB.ID" = readr::col_character(), - "SignDesc" = readr::col_character(), - "StartLoc" = readr::col_integer(), - "StopLoc" = readr::col_integer(), - "Score" = readr::col_double(), - "Status" = readr::col_character(), - "RunDate" = readr::col_character(), - "IPRAcc" = readr::col_character(), - "IPRDesc" = readr::col_character(), - "placeholder" = readr::col_character() - ) -} + hmm_files <- list.files( + db$dir, + pattern = "\\.hmm$", + recursive = TRUE, + full.names = TRUE, + ignore.case = TRUE + ) -#' Internal helpers for reading InterProScan TSV outputs -#' -#' Provide standardized column names, types, and a reader wrapper for the -#' InterProScan tab-delimited output format. -#' -#' @param filepath Path to a `.tsv` or `.tsv.gz` InterProScan result file. -#' -#' @return A tibble of parsed InterProScan output. -#' -#' @keywords internal -.readIPRscanTsv <- function(filepath) { - readr::read_tsv(filepath, - col_types = .getDfIPRColTypes(), - col_names = .getDfIPRColNames() - ) -} + if (length(hmm_files) == 0) { + if (verbose) { + message("Downloading ", db_name) + } -#' Run InterProScan on a sequence chunk inside Docker -#' -#' Executes InterProScan on a subset of protein sequences, writing temporary -#' FASTA and reading back `.tsv` or `.tsv.gz` results. -#' -#' @param chunk A tibble with columns `name` and `sequence`. -#' @param path Working directory used for temporary files. -#' @param ipr_data_path Path to InterProScan data directory. -#' @param out_file_base Output prefix for chunk results. -#' @param appl Character vector of InterProScan applications (e.g. `"Pfam"`). -#' @param chunk_id Integer chunk index. -#' @param threads Number of CPUs for InterProScan container. -#' @param file_format Output format (`"TSV"`). -#' @param docker_image InterProScan Docker image. -#' -#' @return Path to a `.tsv` or `.tsv.gz` InterProScan output file. -#' -#' @keywords internal -.process_chunk <- function(chunk, - path, - ipr_data_path = "inst/extdata/interpro/data", - out_file_base, - appl, - chunk_id, - threads, - file_format, - docker_image = sprintf("interpro/interproscan:%s", "5.76-107.0")) { - # Normalize and mount paths - dir.create(file.path(path, "tmp", "iprscan"), recursive = TRUE, showWarnings = FALSE) - path <- .docker_path(path) - bind_data <- .docker_path(ipr_data_path) - - fasta_sequences <- Biostrings::AAStringSet(chunk$sequence) - names(fasta_sequences) <- chunk$name - temp_fasta_file <- tempfile(tmpdir = path, fileext = ".fa") - Biostrings::writeXStringSet(fasta_sequences, temp_fasta_file) - - chunk_out_file_base_host <- file.path(path, sprintf("%s_chunk_%d", out_file_base, chunk_id)) - chunk_out_file_base_cont <- .to_container(chunk_out_file_base_host, path, "/work") - - # Pull image (best-effort) - try(suppressWarnings(system2("docker", args = c("pull", docker_image))), silent = TRUE) - - appl_str <- paste(appl, collapse = ",") + tmp <- tempfile() - cmd_args <- c( - "run", "--rm", - "-v", paste0(path, ":", "/work"), - "-v", paste0(bind_data, ":/opt/interproscan/data"), - "-w", "/work", - docker_image, - "--input", .to_container(temp_fasta_file, path, "/work"), - "--cpu", as.character(threads), - "-f", file_format, - "--appl", appl_str, - "-b", chunk_out_file_base_cont - ) + utils::download.file( + url = db$url, + destfile = tmp, + mode = "wb", + method = "libcurl" + ) + switch( + db$type, + + gz = { + hmm_file <- file.path( + db$dir, + db$hmm_name %||% basename( + sub( + "\\.gz$", + "", + basename(db$url), + ignore.case = TRUE + ) + ) + ) + + R.utils::gunzip( + filename = tmp, + destname = hmm_file, + overwrite = TRUE, + remove = FALSE + ) + }, + + zip = { + utils::unzip( + zipfile = tmp, + exdir = db$dir + ) + }, + + `tar.gz` = { + utils::untar( + tarfile = tmp, + exdir = db$dir + ) + } + ) - status <- tryCatch( - { - system2( + unlink(tmp) + + hmm_files <- list.files( + db$dir, + pattern = "\\.hmm$", + recursive = TRUE, + full.names = TRUE, + ignore.case = TRUE + ) + } + + if (length(hmm_files) == 0) { + stop( + "No .hmm file found for ", + db_name + ) + } + + if (length(hmm_files) == 1) { + + hmm_file <- hmm_files[[1]] + + } else { + + hmm_file <- file.path( + db$dir, + paste0(db_name, ".hmm") + ) + + source_hmms <- setdiff( + normalizePath(hmm_files), + normalizePath(hmm_file, mustWork = FALSE) + ) + + valid_hmms <- purrr::map_lgl( + source_hmms, + .isValidHmmFile + ) + + if (any(!valid_hmms)) { + + bad_files <- basename( + source_hmms[!valid_hmms] + ) + + if (isTRUE(verbose)) { + warning( + "Ignoring ", + length(bad_files), + " invalid HMM file(s):\n", + paste(bad_files, collapse = "\n"), + call. = FALSE + ) + } + + source_hmms <- source_hmms[valid_hmms] + } + + if (length(source_hmms) == 0) { + stop( + "No valid HMM files found for ", + db_name + ) + } + + if (!file.exists(hmm_file)) { + + if (verbose) { + message( + "Combining ", + length(source_hmms), + " HMM files for ", + db_name + ) + } + + file.create(hmm_file) + + for (f in sort(source_hmms)) { + file.append(hmm_file, f) + } + } + } + + pressed_files <- paste0( + hmm_file, + c(".h3m", ".h3i", ".h3f", ".h3p") + ) + + if (!all(file.exists(pressed_files))) { + + if (verbose) { + message( + "Running hmmpress for ", + basename(hmm_file) + ) + } + + output <- system2( "docker", args = c( "run", "--rm", - "--platform", "linux/amd64", # force amd64 for ARM hosts - "-v", paste0(path, ":", "/work"), - "-v", paste0(bind_data, ":/opt/interproscan/data"), - "-w", "/work", + "-v", + paste0(dirname(hmm_file), ":/db"), docker_image, - "--input", .to_container(temp_fasta_file, path, "/work"), - "--cpu", as.character(threads), - "-f", file_format, - "--appl", appl_str, - "-b", chunk_out_file_base_cont + "hmmpress", + file.path("/db", basename(hmm_file)) ), stdout = TRUE, stderr = TRUE ) - }, - error = function(e) { - stop(sprintf("InterProScan execution failed for chunk %d: %s", chunk_id, e$message)) + + if (!all(file.exists(pressed_files))) { + stop( + "hmmpress failed for ", + db_name, + "\n", + paste(output, collapse = "\n") + ) + } } - ) - out_tsv <- paste0(chunk_out_file_base_host, ".tsv") - out_tsvgz <- paste0(chunk_out_file_base_host, ".tsv.gz") + db_paths[[db_name]] <- list( + hmm = hmm_file, + source = db$url, + type = db$type, + pressed = pressed_files + ) - if (file.exists(out_tsv)) { - return(out_tsv) - } else if (file.exists(out_tsvgz)) { - return(out_tsvgz) - } else { - stop(sprintf( - "InterProScan produced no output for chunk %d. Checked: %s and %s.\nLast message:\n%s", - chunk_id, out_tsv, out_tsvgz, paste(status, collapse = "\n") - )) + if (verbose) { + message( + db_name, + " ready: ", + hmm_file + ) + } } -} -#' Derive protein domain presence/absence and counts via InterProScan and write to DuckDB -domainFromIPR <- function(duckdb_path, - path, - out_file_base = "iprscan", - appl = c("Pfam"), - ipr_version = "5.76-107.0", - ipr_dest_dir = "inst/extdata/interpro", - ipr_platform = "linux/amd64", - auto_prepare_data = TRUE, - threads = 8, - file_format = "TSV", - docker_repo = "interpro/interproscan") { - duckdb_path <- normalizePath(duckdb_path) - if (missing(path) || path %in% c(".", "results", "results/")) { - path <- dirname(duckdb_path) - } - dir.create(path, recursive = TRUE, showWarnings = FALSE) - path <- normalizePath(path) + db_paths +} - ipr_image <- sprintf("%s:%s", docker_repo, ipr_version) - # Prepare data if needed - ipr_info <- if (isTRUE(auto_prepare_data)) { - .checkInterProData( - version = ipr_version, - dest_dir = ipr_dest_dir, - docker_image = ipr_image, - platform = ipr_platform, - verbose = TRUE - ) - } else { - list( - data_dir = file.path(ipr_dest_dir, sprintf("interproscan-%s", ipr_version), "data"), - ready = NA - ) - } - ipr_data_path <- ipr_info$data_dir - # Pull image once - try(suppressWarnings(system2("docker", args = c("pull", ipr_image))), silent = TRUE) - con <- DBI::dbConnect(duckdb::duckdb(), duckdb_path) - on.exit(try(DBI::dbDisconnect(con, shutdown = FALSE), silent = TRUE), add = TRUE) - sequences_df <- dplyr::tbl(con, "protein_cluster_seq") |> tibble::as_tibble() - if (nrow(sequences_df) == 0L) { - stop("No sequences found in 'protein_cluster_seq'. Please run CDHIT2duckdb() first.") - } +#' The function to run HMMER with docker +#' +#' @param JOB_NAME +#' @param FASTA +#' @param DB +#' @param Total_proteins +#' @param output_path +#' @param db_paths +#' @param docker_image +#' @param threads +#' @param n_workers +#' +#' @returns +#' +#' @keywords internal +.runHmmerJob <- function(JOB_NAME, FASTA, DB, Total_proteins, + output_path = NULL, db_paths, + docker_image = "staphb/hmmer", threads = 8L, + n_workers = 8L, + verbose = TRUE +) { + hmmer_input <- file.path(output_path, FASTA) + hmmer_output <- file.path(output_path, paste0(JOB_NAME, ".tbl")) - # Chunking for parallel (not currently implemented due to memory limits) - chunks <- list(sequences_df) # Force 1 chunk for RAM limits + # database paths + database_path <- db_paths[[DB]]$hmm + db_host_dir <- dirname(database_path) + db_filename <- basename(database_path) + db_cont_dir <- "/opt/hmmer/data" + db_cont_path <- file.path(db_cont_dir, db_filename) - # Forcing 1 container operation for RAM limits - cpu_per_container <- threads + # mounts + mount_host <- output_path + mount_cont <- "/work" - message(sprintf( - "InterPro: running in single-container mode with %d CPU(s).", - cpu_per_container - )) + threads_per_job <- max( + 1L, + floor(threads / n_workers) + ) - old_plan <- future::plan() - on.exit(future::plan(old_plan), add = TRUE) - future::plan(future::sequential) + cmd_args <- c( + "run", "--rm", + "-v", paste0(mount_host, ":", mount_cont), + "-v", paste0(db_host_dir, ":", db_cont_dir), + docker_image, + "hmmsearch", + "--notextw", + "--cpu", as.character(threads_per_job), + "-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) + ) - 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 + if(verbose) message("Running hmmsearch via Docker...") + output <- tryCatch( + { + system2("docker", args = cmd_args, stdout = TRUE, stderr = TRUE) }, - .options = furrr::furrr_options(seed = TRUE) + error = function(e) { + stop("hmmsearch execution failed: ", e$message) + } ) - # Combine results - tsvs <- Filter(function(x) !is.null(x) && file.exists(x), results) - if (length(tsvs) == 0L) { - stop("InterProScan produced no usable outputs. Check Docker logs above.") + if (!file.exists(hmmer_output)) { + stop("hmmsearch failed: output file not found. Check stderr:\n", paste(output, collapse = "\n")) } - df_iprscan <- purrr::map(tsvs, .readIPRscanTsv) |> purrr::list_rbind() - - # Load processed tables (unchanged) - DBI::dbWriteTable(con, "domain_names", - df_iprscan |> - dplyr::select(AccNum, DB.ID, SignDesc, IPRAcc, IPRDesc, StartLoc, StopLoc), - overwrite = TRUE - ) - - df_protein_domain_pa <- df_iprscan |> - dplyr::select(AccNum, DB.ID, IPRAcc, placeholder) |> - dplyr::mutate(domain_ID = stringr::str_glue("{DB.ID}_{IPRAcc}")) |> - dplyr::distinct() |> - dplyr::mutate(placeholder = stringr::str_replace_all(placeholder, "-", "1")) |> - tidyr::pivot_wider( - id_cols = AccNum, names_from = domain_ID, values_from = placeholder, - values_fill = "0" - ) |> - dplyr::group_by(AccNum) |> - dplyr::summarize(across(everything(), ~ ifelse(any(. == "1"), "1", "0")), .groups = "drop") |> - dplyr::mutate(across(-AccNum, as.numeric)) + if(verbose) message("hmmsearch completed successfully.") - protein_filter <- dplyr::tbl(con, "protein_count") |> tibble::as_tibble() - accs <- unique(df_protein_domain_pa$AccNum) - accs_in_matrix <- intersect(accs, colnames(protein_filter)) - if (length(accs_in_matrix) == 0L) { - stop("No InterPro accessions match protein_count columns.") - } + # Adding an E value cutoff here + hmmer_tbl <- .parseHMMEROutput(hmmer_output) |> + dplyr::filter(i_evalue <= 1e-5) |> + dplyr::select( + protein, + query_name, + query_accession, + target_description, + i_evalue, + domain_score + ) - protein_filter <- protein_filter |> dplyr::select(genome_id, dplyr::all_of(accs_in_matrix)) - df_protein_domain_pa <- df_protein_domain_pa |> - dplyr::filter(AccNum %in% accs_in_matrix) |> - dplyr::arrange(match(AccNum, accs_in_matrix)) + hmmer_tbl_filename <- file.path( + dirname(hmmer_output), + paste0(tools::file_path_sans_ext(basename(hmmer_output)), ".parquet") + ) - domain_count <- as.matrix(protein_filter |> dplyr::select(-genome_id)) %*% - as.matrix(df_protein_domain_pa |> dplyr::select(-AccNum)) |> - tibble::as_tibble() |> - dplyr::mutate(genome_id = protein_filter |> dplyr::pull(genome_id)) |> - dplyr::relocate(genome_id, .before = dplyr::everything()) + .write_compressed_parquet(hmmer_tbl, hmmer_tbl_filename) - DBI::dbWriteTable(conn = con, name = "domain_count", domain_count, overwrite = TRUE) - invisible(TRUE) + hmmer_tbl_filename } -# Clean BV-BRC metadata, then save as Parquet files + +#' Wrapper for preparing HMM databases and running HMMER on protein sequences from duckdb and writing them. #' #' @param duckdb_path -#' @param path -#' @param ref_file_path +#' @param output_path +#' @param threads +#' @param hmmer_db_dir +#' @param databases +#' @param docker_image +#' @param num_of_splits +#' @param n_workers #' #' @returns -#' @export #' +#' @keywords internal #' @examples -cleanMetaData <- function(duckdb_path, path, ref_file_path = "data_raw/") { - duckdb_path <- normalizePath(duckdb_path) - # If no explicit path is provided (or a generic one), choose results// when - # the DuckDB lives under data//, or else fall back to the DuckDB directory. - if (missing(path) || path %in% c(".", "results", "results/")) { - bug_dir <- dirname(duckdb_path) - mapped_results <- sub( - paste0(.Platform$file.sep, "data", .Platform$file.sep), - paste0(.Platform$file.sep, "results", .Platform$file.sep), - bug_dir, - fixed = TRUE +.runHMMER <- function(duckdb_path, + output_path, + threads = 8L, + hmmer_db_dir, + databases = c("Pfam", "COG", "AMRFinder"), + docker_image = "staphb/hmmer", + num_of_splits = 8L, + n_workers = 8L, + verbose = TRUE +) { + # Fail fast if Docker is missing + if (!nzchar(Sys.which("docker"))) { + stop("Docker is not available on your PATH but is required to run HMMER.") + } + + # But also check if Docker is on the PATH but isn't running + docker_ok <- system2( + "docker", + "info", + stdout = FALSE, + stderr = FALSE + ) == 0L + + if (!docker_ok) { + stop( + "Docker is installed but is not running or cannot be reached. ", + "Please (re)start Docker Desktop and try again." ) - path <- if (!identical(mapped_results, bug_dir)) mapped_results else bug_dir } - path <- normalizePath(path, mustWork = FALSE) - if (!dir.exists(path)) dir.create(path, recursive = TRUE) + duckdb_path <- .docker_path(duckdb_path) + if (missing(output_path) || output_path %in% c(".", "results", "results/")) { + output_path <- dirname(duckdb_path) + } + output_path <- .docker_path(output_path) + if (!dir.exists(output_path)) dir.create(output_path, recursive = TRUE) con <- DBI::dbConnect(duckdb::duckdb(), duckdb_path) on.exit(try(DBI::dbDisconnect(con, shutdown = FALSE), silent = TRUE), add = TRUE) - ref_file_path <- normalizePath(ref_file_path) - clean_drug <- readr::read_tsv(file.path(ref_file_path, "clean_drug.tsv")) - drug_class <- readr::read_tsv(file.path(ref_file_path, "drug_class.tsv")) - drug_abbr <- readr::read_tsv(file.path(ref_file_path, "drug_abbr.tsv")) - class_abbr <- readr::read_tsv(file.path(ref_file_path, "class_abbr.tsv")) - clean_countries <- readr::read_tsv(file.path(ref_file_path, "cleaned_bvbrc_countries.tsv")) |> - dplyr::select("raw_entry", "clean_name", "short_name") |> - dplyr::distinct() + prot_seqs <- DBI::dbReadTable(con, "protein_cluster_seq") |> + tibble::as_tibble() - # Define lab methods - lab_methods <- c("Disk diffusion", "MIC", "Broth dilution", "Agar dilution", "Biofosun Gram-positive panels broth dilution", - "Vitek_2-P607_card", "cation-adjusted Mueller-Hinton broth", "gradient_diffusion", "kirby-bauer_disc_diffusion") + # Just in case CD-HIT failed to generate sequences somehow + if (nrow(prot_seqs) == 0L) { + stop("No sequences found in 'protein_cluster_seq'. Please run CDHIT2duckdb() first.") + } - dplyr::tbl(con, "filtered") |> - tibble::as_tibble() |> - dplyr::select("genome.genome_id") |> - dplyr::left_join(dplyr::tbl(con, "metadata") |> - tibble::as_tibble(), by = dplyr::join_by("genome.genome_id" == "genome_drug.genome_id")) |> - dplyr::select( - "genome.genome_id", "genome_drug.antibiotic", - "genome_drug.genome_name", "genome_drug.evidence", "genome_drug.laboratory_typing_method", + # required to define the database size for hmmsearch --Z and --domZ parameters + Total_proteins <- nrow(prot_seqs) + + if (is.null(hmmer_db_dir)) { + hmmer_db_dir <- .defaultHmmerDbDir() + } + + dir.create( + hmmer_db_dir, + recursive = TRUE, + showWarnings = FALSE + ) + + # database paths + if(verbose) message ("Preparing HMM databases") + db_paths <- .prepareHmmerDatabases( + hmmer_db_dir = hmmer_db_dir, + databases = databases, + docker_image = docker_image, + verbose = verbose + ) + + db_paths <- db_paths[databases] + + # clamp splits to the number of sequences available + chunk_count <- min(as.integer(num_of_splits), nrow(prot_seqs)) + + split_fasta <- function(seqs, prefix) { + records <- paste0(">", seqs$name, "\n", seqs$sequence) + chunk_size <- ceiling(length(records) / chunk_count) + chunks <- split(records, ceiling(seq_along(records) / chunk_size)) + + purrr::walk2(chunks, seq_along(chunks), function(chunk, i) { + chunk_path <- file.path(output_path, sprintf("%s_chunk_%02d.fasta", prefix, i)) + readr::write_lines(chunk, chunk_path) + }) + } + + split_fasta(prot_seqs, "protein") + + job_list <- expand.grid( + chunk = sprintf("%02d", seq_len(chunk_count)), + db = databases, + stringsAsFactors = FALSE + ) |> + dplyr::mutate( + JOB_NAME = paste0("protein_chunk_", chunk, "_", db), + FASTA = paste0("protein_chunk_", chunk, ".fasta"), + DB = db + ) |> + dplyr::select(JOB_NAME, FASTA, DB) + + old_plan <- future::plan() + on.exit(future::plan(old_plan), add = TRUE) + + future::plan( + future::multisession, + workers = max(1L, n_workers) + ) + + if (verbose) message("Running HMMER jobs") + parquet_files <- furrr::future_map_chr( + seq_len(nrow(job_list)), + function(i) { + + .runHmmerJob( + JOB_NAME = job_list$JOB_NAME[i], + FASTA = job_list$FASTA[i], + DB = job_list$DB[i], + Total_proteins = Total_proteins, + output_path = output_path, + db_paths = db_paths, + docker_image = docker_image, + threads = threads, + n_workers = n_workers, + verbose = verbose + ) + } + ) + + parquet_tbl <- tibble::tibble( + parquet = parquet_files, + db = job_list$DB + ) + + final_parquets <- list() + + for (database_name in databases) { + + if(verbose) message("Combining ", database_name) + + db_files <- parquet_tbl |> + dplyr::filter( + db == database_name + ) |> + dplyr::pull(parquet) + + combined_tbl <- purrr::map( + db_files, + arrow::read_parquet + ) |> + dplyr::bind_rows() |> + dplyr::left_join(.parse_hmmer_profiles(db_paths[[database_name]]$hmm) |> + dplyr::select(query_name = profile_name, query_accession = profile_accession, description = profile_description), + by = "query_name") + + final_parquet <- file.path( + output_path, + paste0( + "protein_", + database_name, + ".parquet" + ) + ) + + .write_compressed_parquet( + combined_tbl, + final_parquet + ) + + DBI::dbWriteTable( + con, + name = paste0( + "protein_", + database_name + ), + value = combined_tbl, + overwrite = TRUE + ) + + final_parquets[[database_name]] <- final_parquet + + message( + "Created ", + basename(final_parquet) + ) + } + + unlink( + list.files( + output_path, + pattern = "^protein_chunk_.*\\.(fasta|tbl|parquet)$", + full.names = TRUE + ) + ) + + invisible(list( + databases = db_paths, + outputs = final_parquets + )) + + # purrr::map(parquet_files, arrow::read_parquet) |> + # dplyr::bind_rows() |> + # .write_compressed_parquet(final_parquet) + + # message("Combined parquet written.") + + # arrow::read_parquet(final_parquet) |> + # DBI::dbWriteTable(conn = con, name = tools::file_path_sans_ext(basename(final_parquet)), overwrite = TRUE) +} + +#' Map HMMER protein annotations to genome-level count matrix and load into DuckDB +#' +#' Reads a Parquet file of HMMER hits (produced by [.runHMMER()]), joins the +#' annotations to the protein-cluster count matrix already in DuckDB, aggregates +#' 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()]). +#' +#' @return Invisibly returns the path to the written count Parquet file. +#' +#' @seealso [CDHIT2duckdb()], [runDataProcessing()] +#' +#' @keywords internal +.proteinAnnotations2Duckdb <- function( + duckdb_path, + databases, + output_path = dirname(duckdb_path) +) { + + duckdb_path <- .docker_path(duckdb_path) + + con <- DBI::dbConnect( + duckdb::duckdb(), + duckdb_path + ) + + on.exit( + try( + DBI::dbDisconnect( + con, + shutdown = FALSE + ), + silent = TRUE + ), + add = TRUE + ) + + protein_long <- DBI::dbReadTable( + con, + "protein_count" + ) |> + tibble::as_tibble() |> + tidyr::pivot_longer( + cols = -genome_id, + names_to = "protein", + values_to = "count" + ) |> + dplyr::filter(count > 0) |> + dplyr::mutate( + protein = stringr::str_replace( + protein, + "^fig\\.", + "fig|" + ) + ) + + count_paths <- list() + + for (database in databases) { + + annotation_table <- paste0( + "protein_", + database + ) + + if (!DBI::dbExistsTable(con, annotation_table)) { + + warning( + annotation_table, + " not found in DuckDB. Skipping." + ) + + next + } + + message( + "Processing ", + annotation_table + ) + + annotation <- DBI::dbReadTable( + con, + annotation_table + ) |> + tibble::as_tibble() |> + dplyr::distinct( + protein, + query_name + ) + + genome_annot_matrix <- protein_long |> + dplyr::inner_join( + annotation |> + dplyr::select( + protein, + query_name + ), + by = "protein", + relationship = "many-to-many" + ) |> + dplyr::group_by( + genome_id, + query_name + ) |> + dplyr::summarise( + count = sum(count), + .groups = "drop" + ) |> + tidyr::pivot_wider( + names_from = query_name, + values_from = count, + values_fill = 0 + ) + + count_table <- paste0( + annotation_table, + "_count" + ) + + count_path <- file.path( + output_path, + paste0( + count_table, + ".parquet" + ) + ) + + arrow::write_parquet( + genome_annot_matrix, + count_path + ) + + DBI::dbWriteTable( + con, + count_table, + genome_annot_matrix, + overwrite = TRUE + ) + + count_paths[[database]] <- count_path + + message( + "Created ", + count_table + ) + } + + invisible(count_paths) +} + +#' Annotate proteins using DefenseFinder + CasFinder HMMs +#' Will add to the duckdb + create the parquet file. +#' +#' @param defense_db_dir Directory used to store downloaded HMMs +#' @param docker_image Docker image containing HMMER +#' @param duckdb_path DuckDB database path +#' @param output_path Output directory +#' @param threads Number of HMMER threads +#' +#' @returns Path to annotation parquet +#' @keywords internal +.defenseHMMER <- function( + defense_db_dir, + docker_image = "staphb/hmmer", + duckdb_path = "inst/extdata/Sfl.duckdb", + output_path = NULL, + threads = 8L, + verbose = TRUE +) { + + if (!nzchar(Sys.which("docker"))) { + stop("Docker is required.") + } + + defense_db_dir <- normalizePath( + defense_db_dir, + mustWork = FALSE + ) + + if (is.null(output_path)) { + output_path <- dirname( + normalizePath( + duckdb_path, + mustWork = FALSE + ) + ) + } + + dir.create( + defense_db_dir, + recursive = TRUE, + showWarnings = FALSE + ) + + dir.create( + output_path, + recursive = TRUE, + showWarnings = FALSE + ) + + #################################################################### + # download repositories + #################################################################### + + defense_dir <- file.path( + defense_db_dir, + "DefenseFinder" + ) + + cas_dir <- file.path( + defense_db_dir, + "CasFinder" + ) + + if (!dir.exists(defense_dir)) { + + if(verbose) message( + "Downloading DefenseFinder models" + ) + + tmp <- tempfile(fileext = ".zip") + + utils::download.file( + "https://github.com/mdmparis/defense-finder-models/archive/refs/heads/master.zip", + tmp, + mode = "wb", + method = "libcurl" + ) + + utils::unzip( + tmp, + exdir = defense_dir + ) + + unlink(tmp) + } + + if (!dir.exists(cas_dir)) { + + if(verbose) message( + "Downloading CasFinder models" + ) + + tmp <- tempfile(fileext = ".zip") + + utils::download.file( + "https://github.com/macsy-models/CasFinder/archive/refs/heads/main.zip", + tmp, + mode = "wb", + method = "libcurl" + ) + + utils::unzip( + tmp, + exdir = cas_dir + ) + + unlink(tmp) + } + + #################################################################### + # helper + #################################################################### + + build_database <- function( + repo_dir, + db_name + ) { + + profile_dirs <- list.dirs( + repo_dir, + recursive = TRUE, + full.names = TRUE + ) + + profile_dirs <- profile_dirs[ + basename(profile_dirs) == "profiles" + ] + + # moving to purrr implementation + hmm_files <- profile_dirs |> + purrr::map(\(x) list.files(x, + pattern = "\\.hmm$", + recursive = TRUE, + full.names = TRUE, + ignore.case = TRUE)) |> + purrr::flatten_chr() |> + unique() + + if (length(hmm_files) == 0) { + + stop( + "No HMM files found for ", + db_name + ) + } + + valid_hmms <- purrr::map_lgl( + hmm_files, + .isValidHmmFile + ) + + if (any(!valid_hmms)) { + bad_files <- basename(hmm_files[!valid_hmms]) + + if (isTRUE(verbose)) { + warning( + "Ignoring ", + length(bad_files), + " invalid HMM file(s):\n", + paste(bad_files, collapse = "\n"), + call. = FALSE + ) + } + + hmm_files <- hmm_files[valid_hmms] + } + + if (length(hmm_files) == 0) { + stop( + "No valid HMM files found for ", + db_name + ) + } + + combined_hmm <- file.path( + repo_dir, + paste0( + db_name, + ".hmm" + ) + ) + + if (file.exists(combined_hmm)) { + unlink(combined_hmm) + } + + file.create(combined_hmm) + + for (f in sort(hmm_files)) { + + file.append( + combined_hmm, + f + ) + } + + pressed_files <- paste0( + combined_hmm, + c( + ".h3m", + ".h3i", + ".h3f", + ".h3p" + ) + ) + + if (!all(file.exists(pressed_files))) { + + if(verbose) message( + "Running hmmpress for ", + db_name + ) + + output <- system2( + "docker", + args = c( + "run", + "--rm", + "-v", + paste0( + dirname(combined_hmm), + ":/db" + ), + docker_image, + "hmmpress", + file.path( + "/db", + basename(combined_hmm) + ) + ), + stdout = TRUE, + stderr = TRUE + ) + + if (!all(file.exists(pressed_files))) { + + stop( + "hmmpress failed for ", + db_name, + "\n", + paste(output, + collapse = "\n") + ) + } + } + + combined_hmm + } + + #################################################################### + # build separate databases + #################################################################### + + defense_hmm <- build_database( + defense_dir, + "DefenseFinder" + ) + + cas_hmm <- build_database( + cas_dir, + "CasFinder" + ) + + #################################################################### + # load proteins + #################################################################### + + con <- DBI::dbConnect( + duckdb::duckdb(), + duckdb_path + ) + + on.exit( + try( + DBI::dbDisconnect( + con, + shutdown = FALSE + ), + silent = TRUE + ), + add = TRUE + ) + + prot_seqs <- DBI::dbReadTable( + con, + "protein_cluster_seq" + ) |> + tibble::as_tibble() + + fasta_file <- file.path( + output_path, + "protein_DefenseCas.faa" + ) + + # required to define the database size for hmmsearch --Z and --domZ parameters + Total_proteins <- nrow(prot_seqs) + + readr::write_lines( + paste0( + ">", + prot_seqs$name, + "\n", + prot_seqs$sequence + ), + fasta_file + ) + + #################################################################### + # run hmmsearch separately + #################################################################### + + databases <- list( + DefenseFinder = defense_hmm, + CasFinder = cas_hmm + ) + + all_hits <- list() + + for (db_name in names(databases)) { + + if(verbose) message( + "Running ", + db_name + ) + + hmm_file <- databases[[db_name]] + + tbl_file <- file.path( + output_path, + paste0( + "protein_", + db_name, + ".tbl" + ) + ) + + output <- system2( + "docker", + args = c( + "run", + "--rm", + "-v", + paste0(output_path, ":/work"), + "-v", + paste0(dirname(hmm_file), ":/db"), + docker_image, + "hmmsearch", + "--notextw", + "--cpu", + as.character(threads), + "-Z", Total_proteins, + "--domZ", Total_proteins, + "--domtblout", + file.path( + "/work", + basename(tbl_file) + ), + file.path( + "/db", + basename(hmm_file) + ), + "/work/protein_DefenseCas.faa" + ), + stdout = TRUE, + stderr = TRUE + ) + + if (!file.exists(tbl_file)) { + stop( + "hmmsearch failed for ", + db_name, + "\n", + paste(output, collapse = "\n") + ) + } + + hits <- .parseHMMEROutput( + tbl_file + ) |> + dplyr::select( + protein, + query_name + ) |> + dplyr::mutate( + database = db_name + )|> + dplyr::left_join(.parse_hmmer_profiles(hmm_file) |> + dplyr::select(query_name = profile_name, query_accession = profile_accession, description = profile_description), + by = "query_name") + + all_hits[[db_name]] <- hits + } + + #################################################################### + # merge at parquet stage + #################################################################### + + combined_tbl <- dplyr::bind_rows( + all_hits + ) + + parquet_file <- file.path( + output_path, + "protein_DefenseCas.parquet" + ) + + .write_compressed_parquet( + combined_tbl, + parquet_file + ) + + DBI::dbWriteTable( + con, + "protein_DefenseCas", + combined_tbl, + overwrite = TRUE + ) + + message( + "Created protein_DefenseCas" + ) + + unlink( + c( + fasta_file, + file.path( + output_path, + paste0("protein_", names(databases), ".tbl") + ) + ) + ) + + invisible(list( + databases = list( + DefenseFinder = defense_hmm, + CasFinder = cas_hmm + ), + output = parquet_file + )) +} + + +# Clean BV-BRC metadata, then save as Parquet files +#' +#' @param duckdb_path +#' @param path +#' @param ref_file_path +#' +#' @returns +#' @export +#' +#' @examples +cleanMetaData <- function(duckdb_path, path, ref_file_path = "data_raw/") { + duckdb_path <- normalizePath(duckdb_path) + # If no explicit path is provided (or a generic one), choose results// when + # the DuckDB lives under data//, or else fall back to the DuckDB directory. + if (missing(path) || path %in% c(".", "results", "results/")) { + bug_dir <- dirname(duckdb_path) + mapped_results <- sub( + paste0(.Platform$file.sep, "data", .Platform$file.sep), + paste0(.Platform$file.sep, "results", .Platform$file.sep), + bug_dir, + fixed = TRUE + ) + path <- if (!identical(mapped_results, bug_dir)) mapped_results else bug_dir + } + + path <- normalizePath(path, mustWork = FALSE) + if (!dir.exists(path)) dir.create(path, recursive = TRUE) + + con <- DBI::dbConnect(duckdb::duckdb(), duckdb_path) + on.exit(try(DBI::dbDisconnect(con, shutdown = FALSE), silent = TRUE), add = TRUE) + ref_file_path <- normalizePath(ref_file_path) + + clean_drug <- readr::read_tsv(file.path(ref_file_path, "clean_drug.tsv")) + drug_class <- readr::read_tsv(file.path(ref_file_path, "drug_class.tsv")) + drug_abbr <- readr::read_tsv(file.path(ref_file_path, "drug_abbr.tsv")) + class_abbr <- readr::read_tsv(file.path(ref_file_path, "class_abbr.tsv")) + clean_countries <- readr::read_tsv(file.path(ref_file_path, "cleaned_bvbrc_countries.tsv")) |> + dplyr::select("raw_entry", "clean_name", "short_name") |> + dplyr::distinct() + + # Define lab methods + lab_methods <- c("Disk diffusion", "MIC", "Broth dilution", "Agar dilution", "Biofosun Gram-positive panels broth dilution", + "Vitek_2-P607_card", "cation-adjusted Mueller-Hinton broth", "gradient_diffusion", "kirby-bauer_disc_diffusion") + + dplyr::tbl(con, "filtered") |> + tibble::as_tibble() |> + dplyr::select("genome.genome_id") |> + dplyr::left_join(dplyr::tbl(con, "metadata") |> + tibble::as_tibble(), by = dplyr::join_by("genome.genome_id" == "genome_drug.genome_id")) |> + dplyr::select( + "genome.genome_id", "genome_drug.antibiotic", + "genome_drug.genome_name", "genome_drug.evidence", "genome_drug.laboratory_typing_method", "genome_drug.resistant_phenotype", "genome_drug.taxon_id", "genome_drug.pmid", "genome.collection_year", "genome.isolation_country", "genome.host_common_name", @@ -1665,9 +2377,64 @@ cleanData <- function(duckdb_path, path) { path <- normalizePath(path, mustWork = FALSE) if (!dir.exists(path)) dir.create(path, recursive = TRUE) + # Fun new manifest action allows cleanData to find applicable database names + 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 + ) + + hmmer_stage <- NULL + + for (run in rev(manifest$runs)) { + 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 <- unlist( + hmmer_stage$parameters$databases + ) + + if (!length(hmmer_databases)) { + stop( + "HMMER stage in manifest does not contain any databases." + ) + } + con <- DBI::dbConnect(duckdb::duckdb(), duckdb_path) on.exit(try(DBI::dbDisconnect(con, shutdown = FALSE), silent = TRUE), add = TRUE) + .proteinAnnotations2Duckdb( + duckdb_path = duckdb_path, + databases = hmmer_databases, + output_path = path + ) + # Parquet output paths genes_parquet <- file.path(path, "gene_count.parquet") gene_names_parquet <- file.path(path, "gene_names.parquet") @@ -1676,9 +2443,7 @@ cleanData <- function(duckdb_path, path) { struct_parquet <- file.path(path, "struct.parquet") proteins_parquet <- file.path(path, "protein_count.parquet") - domains_parquet <- file.path(path, "domain_count.parquet") - domain_names_parquet <- file.path(path, "domain_names.parquet") protein_names_parquet <- file.path(path, "protein_names.parquet") protein_cluster_seq_parquet <- file.path(path, "protein_seqs.parquet") @@ -1720,13 +2485,39 @@ cleanData <- function(duckdb_path, path) { writeCompressedParquet(proteins_parquet) DBI::dbExecute(con_new, sprintf("CREATE OR REPLACE VIEW protein_count AS SELECT * FROM read_parquet('%s')", basename(proteins_parquet))) - # domain_count -> long parquet + view - DBI::dbReadTable(con, "domain_count") |> - tidyr::pivot_longer(-genome_id, names_to = "domain", values_to = "value") |> - dplyr::filter(!is.na(value) & value != "") |> - dplyr::mutate(value = as.integer(value)) |> - writeCompressedParquet(domains_parquet) - DBI::dbExecute(con_new, sprintf("CREATE OR REPLACE VIEW domain_count AS SELECT * FROM read_parquet('%s')", basename(domains_parquet))) + # HMMER annotation counts -> long Parquet + views per database in manifest + for (database in hmmer_databases) { + + count_table <- paste0( + "protein_", + database, + "_count" + ) + + count_parquet <- file.path( + path, + paste0(count_table, ".parquet") + ) + + DBI::dbReadTable(con, count_table) |> + tidyr::pivot_longer( + -genome_id, + names_to = "annotation", + values_to = "value" + ) |> + dplyr::filter(!is.na(value) & value != "") |> + dplyr::mutate(value = as.integer(value)) |> + writeCompressedParquet(count_parquet) + + DBI::dbExecute( + con_new, + sprintf( + "CREATE OR REPLACE VIEW %s AS SELECT * FROM read_parquet('%s')", + count_table, + basename(count_parquet) + ) + ) + } # gene_struct -> long parquet + view DBI::dbReadTable(con, "gene_struct") |> @@ -1745,10 +2536,31 @@ cleanData <- function(duckdb_path, path) { writeCompressedParquet(protein_names_parquet) DBI::dbExecute(con_new, sprintf("CREATE OR REPLACE VIEW protein_names AS SELECT * FROM read_parquet('%s')", basename(protein_names_parquet))) - DBI::dbReadTable(con, "domain_names") |> - dplyr::select(-c(IPRAcc, IPRDesc)) |> - writeCompressedParquet(domain_names_parquet) - DBI::dbExecute(con_new, sprintf("CREATE OR REPLACE VIEW domain_names AS SELECT * FROM read_parquet('%s')", basename(domain_names_parquet))) + # Parsing through the different HMMER result Parquets + for (database in hmmer_databases) { + + annotation_table <- paste0( + "protein_", + database + ) + + annotation_parquet <- file.path( + path, + paste0(annotation_table, ".parquet") + ) + + DBI::dbReadTable(con, annotation_table) |> + writeCompressedParquet(annotation_parquet) + + DBI::dbExecute( + con_new, + sprintf( + "CREATE OR REPLACE VIEW %s AS SELECT * FROM read_parquet('%s')", + annotation_table, + basename(annotation_parquet) + ) + ) + } DBI::dbReadTable(con, "gene_ref_seq") |> writeCompressedParquet(gene_ref_seq_parquet) DBI::dbExecute(con_new, sprintf("CREATE OR REPLACE VIEW gene_seqs AS SELECT * FROM read_parquet('%s')", basename(gene_ref_seq_parquet))) @@ -1766,137 +2578,143 @@ cleanData <- function(duckdb_path, path) { } -#' Run the full amRdata processing pipeline (Panaroo → CD-HIT → InterProScan → Parquet) +#' Run the full amRdata processing pipeline (Panaroo -> CD-HIT -> HMMER -> Parquet) #' #' @description #' `runDataProcessing()` orchestrates the complete feature-extraction pipeline for a -#' BV-BRC selection, starting from a **per-selection DuckDB** (created by -#' [prepareGenomes()] and populated by downstream steps). It: +#' BV-BRC selection, starting from a **per-selection DuckDB** created by +#' [prepareGenomes()] and populated by downstream genome processing steps. It: #' 1. Runs **Panaroo** to build the pangenome and writes gene/struct outputs into DuckDB. #' 2. Runs **CD-HIT** to cluster proteins and writes protein outputs into DuckDB. -#' 3. Runs **InterProScan** (Pfam) to annotate protein domains and writes domain outputs into DuckDB. +#' 3. Runs **HMMER** against the requested protein databases and writes annotation +#' tables into DuckDB. #' 4. **Cleans BV-BRC metadata** (drug names/classes, countries, years) and -#' exports all feature/metadata tables as **compressed Parquet** files, then creates -#' a **Parquet-backed DuckDB** with read-only views of those Parquets for downstream ML. +#' exports feature and metadata tables as compressed Parquet files, then creates +#' a **Parquet-backed DuckDB** with read-only views for downstream ML. #' #' The function is a thin controller that delegates each stage to the corresponding -#' internal helpers (Dockerized tools where applicable) and ensures consistent -#' output locations and table schemas across stages. +#' internal helpers (Dockerized tools where applicable) and records processing +#' parameters and provenance in the dataset manifest. #' #' @section Pipeline Steps: #' \enumerate{ -#' \item **Panaroo** via runPanaroo2Duckdb() → writes: +#' \item **Panaroo** via [runPanaroo2Duckdb()] -> writes: #' \itemize{ -#' \item `gene_count` (genome × gene counts)\cr +#' \item `gene_count` (genome x gene counts)\cr #' \item `gene_names`\cr #' \item `gene_struct` (structural variants)\cr #' \item `gene_ref_seq`, `genome_gene_protein` #' } -#' \item **CD-HIT** via CDHIT2duckdb() (calls internal `.runCDHIT()`) → writes: +#' \item **CD-HIT** via [CDHIT2duckdb()] -> writes: #' \itemize{ -#' \item `protein_count` (genome × protein-cluster counts)\cr +#' \item `protein_count` (genome x protein-cluster counts)\cr #' \item `protein_names`\cr -#' \item `protein_cluster_seq` (representative sequences) +#' \item `protein_cluster_seq` (representative sequences)\cr +#' \item `protein_members` #' } -#' \item **InterProScan (Pfam)** via domainFromIPR() → writes: +#' \item **HMMER** via the configured HMMER databases -> writes: #' \itemize{ -#' \item `domain_names`\cr -#' \item `domain_count` (genome × domain-family matrix) -#' } -#' \item **Metadata cleaning + Parquet export** via cleanData() → writes Parquet -#' files to `output_path`, and builds a **Parquet-backed DuckDB** -#' (`*_parquet.duckdb`) with views: -#' \itemize{ -#' \item `gene_count`, `protein_count`, `domain_count`, `struct`\cr -#' \item `metadata` (cleaned), plus `amr_phenotype`, `genome_data`, `original_metadata`\cr -#' \item `gene_names`, `protein_names`, `domain_names`\cr -#' \item `gene_seqs`, `protein_seqs`\cr -#' \item `genome_gene_protein` +#' \item `protein_` annotation tables\cr +#' \item `protein__count` genome-by-annotation count tables +#' \item The default databases are `Pfam`, `COG`, `AMRFinder`, and `DefenseCas`. #' } +#' \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. #' } #' #' @param duckdb_path Character. Path to the **per-selection DuckDB** produced by #' [prepareGenomes()] (e.g., `"data//.duckdb"`). This DB must -#' already contain at least the tables written by `prepareGenomes()` and subsequent -#' download steps (e.g., `files`, `filtered`, and metadata tables). -#' @param output_path Character or `NULL`. Base directory for writing Panaroo/CD-HIT/InterProScan -#' outputs and final Parquet files. If `NULL`, defaults to `dirname(duckdb_path)`. -#' -#' @param threads Integer. Shared concurrency budget used across tools (Panaroo, CD-HIT, -#' InterProScan). Passed through to each stage as appropriate. Defaults to `8`. -#' -#' @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`. +#' already contain the tables written by [prepareGenomes()] and the upstream +#' genome-processing steps. +#' @param output_path Character or `NULL`. Base directory for writing Panaroo, +#' CD-HIT, HMMER, and final Parquet outputs. If `NULL`, defaults to +#' `dirname(duckdb_path)`. +#' +#' @param threads Integer. Shared concurrency budget used across Panaroo, CD-HIT, +#' and HMMER. Defaults to `8`. +#' +#' @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`. #' @param panaroo_core_threshold Numeric. Panaroo `--core_threshold`. Default: `0.90`. #' @param panaroo_len_dif_percent Numeric. Panaroo `--len_dif_percent`. Default: `0.95`. #' @param panaroo_cluster_threshold Numeric. Panaroo `--threshold`. Default: `0.95`. -#' @param panaroo_family_seq_identity Numeric. Panaroo `-f` (gene family identity). Default: `0.5`. -#' @param panaroo_refind_mode Character. Panaroo's `--refind-mode` (`"off"`, `"default"`, -#' or `"strict"`). See [.processPanaroo()] for what refinding does and the runtime -#' caveat behind the current default. Default `"off"`. +#' @param panaroo_family_seq_identity Numeric. Panaroo `-f` gene family identity. +#' Default: `0.5`. +#' @param panaroo_refind_mode Character. Panaroo's `--refind-mode` (`"off"`, +#' `"default"`, or `"strict"`). See [.processPanaroo()] for the runtime caveat +#' behind refinding. Default: `"off"`. +#' @param panaroo_strip_pseudogenes Logical. If `TRUE`, remove pseudogene feature +#' records from Panaroo input GFF files before running Panaroo. Default: `FALSE`. +#' @param panaroo_pseudogene_clean_dir Character. Directory name for cleaned GFF +#' files. Default: `"gff_clean"`. +#' @param panaroo_write_pseudogene_audit Logical. If `TRUE`, write a pseudogene +#' cleaning audit file. Default: `TRUE`. #' #' @param cdhit_identity Numeric. CD-HIT `-c` identity threshold. Default: `0.9`. #' @param cdhit_word_length Integer. CD-HIT `-n` word length. Default: `5`. -#' @param cdhit_memory Integer. CD-HIT `-M` memory limit (MB). Use `0` for unlimited. Default: `0`. -#' @param cdhit_extra_args Character vector. Extra arguments forwarded to `cd-hit` -#' (e.g., `c("-g","1")`). Default: `c("-g","1")`. -#' @param cdhit_output_prefix Character. Prefix for CD-HIT output files. Default: `"cdhit_out"`. -#' -#' @param ipr_appl Character vector. InterProScan applications to run; typically `c("Pfam")`. -#' Default: `c("Pfam")`. -#' @param ipr_threads_unused Deprecated/unused. Kept for backward compatibility; ignored. -#' @param ipr_version Character. InterProScan image tag (e.g., `"5.76-107.0"`). Default: `"5.76-107.0"`. -#' @param ipr_dest_dir Character. Local destination for InterProScan data bundle -#' (used by `.checkInterProData()`). Default: `"inst/extdata/interpro"`. -#' @param ipr_platform Character. Docker platform string for InterProScan containers, -#' e.g., `"linux/amd64"`. Default: `"linux/amd64"`. -#' @param auto_prepare_data Logical. If `TRUE`, ensure InterProScan data are present -#' (download/verify if missing). Default: `TRUE`. -#' -#' @param ref_file_path Character. Directory containing reference TSVs used by cleanData() -#' for metadata harmonization (e.g., `"data_raw/"`). **Required**; defaults to `"data_raw/"`. +#' @param cdhit_memory Integer. CD-HIT `-M` memory limit in MB. Use `0` for +#' unlimited. Default: `0`. +#' @param cdhit_extra_args Character vector. Extra arguments forwarded to +#' `cd-hit`. Default: `c("-g", "1")`. +#' @param cdhit_output_prefix Character. Prefix for CD-HIT output files. +#' Default: `"cdhit_out"`. +#' +#' @param hmmer_databases Character vector. HMMER annotation databases to run. +#' Default: `c("Pfam", "COG", "AMRFinder", "DefenseCas")`. +#' @param hmmer_db_dir Character or `NULL`. Directory containing the shared HMMER +#' database cache. If `NULL`, uses the amRdata user cache. +#' @param hmmer_docker_image Character. Docker image containing HMMER. +#' Default: `"staphb/hmmer"`. +#' @param hmmer_num_splits Integer. Number of protein-sequence chunks for HMMER. +#' Default: `8`. +#' @param hmmer_workers Integer. Number of parallel HMMER workers. Default: `8`. #' +#' @param ref_file_path Character. Directory containing reference TSVs used by +#' [cleanMetaData()] and [cleanData()] for metadata harmonization. +#' Default: `"data_raw/"`. #' @param verbose Logical. Print progress messages. Default: `TRUE`. #' #' @return #' Invisibly returns a list with: #' \itemize{ -#' \item `duckdb_path` – input DuckDB path -#' \item `panaroo_output` – path to the selected Panaroo output directory used for import -#' \item `parquet_duckdb_path` – absolute path to the created Parquet-backed DuckDB +#' \item `duckdb_path` - input DuckDB path +#' \item `panaroo_output` - path to the selected Panaroo output directory used for import +#' \item `parquet_duckdb_path` - absolute path to the created Parquet-backed DuckDB #' } #' #' @details #' **Docker & Platform Notes** -#' * All heavy tools (Panaroo, CD-HIT, InterProScan) run inside Docker containers. -#' * On Apple Silicon/ARM hosts, images are forced to `--platform linux/amd64` to ensure compatibility. -#' * Ensure Docker Desktop is running and has sufficient memory/CPUs configured. +#' * Panaroo, CD-HIT, and HMMER run inside Docker containers. +#' * HMMER databases are stored separately from individual bug directories and +#' are reused across datasets unless a custom `hmmer_db_dir` is supplied. +#' * Ensure Docker Desktop is running and has sufficient memory and CPU resources. #' #' **Input Requirements** -#' * The `duckdb_path` must reference a per-selection DuckDB that contains: -#' `files` (paths to `.gff`, `.fna`, `.PATRIC.faa`), -#' `filtered` (genomes selected for download/filtering), and -#' BV-BRC metadata tables written by earlier steps. +#' * `duckdb_path` must reference a per-selection DuckDB containing the genome +#' file table, filtered genome selection, and BV-BRC metadata produced by the +#' upstream curation workflow. #' #' **Outputs & Side Effects** -#' * Writes tool-specific intermediate outputs under `output_path` (e.g., `panaroo_out_*`, CD-HIT files). -#' * Writes Parquet files to `output_path`: -#' `gene_count.parquet`, `protein_count.parquet`, `domain_count.parquet`, `struct.parquet`, -#' `gene_names.parquet`, `protein_names.parquet`, `domain_names.parquet`, -#' `gene_seqs.parquet`, `protein_seqs.parquet`, `genome_gene_protein.parquet`, -#' `metadata.parquet`, `amr_phenotype.parquet`, `genome_data.parquet`, `original_metadata.parquet`. -#' * Creates a new Parquet-backed DuckDB (`*_parquet.duckdb`) with read-only views pointing to those Parquets. +#' * Writes tool-specific intermediate outputs under `output_path`. +#' * Writes feature and metadata Parquet files under `output_path`. +#' * Creates a new Parquet-backed DuckDB (`*_parquet.duckdb`) with read-only views +#' over the generated Parquet files. +#' * Records processing parameters, software versions, database selections, and +#' other provenance information in the dataset manifest. #' #' **Threading** -#' * `threads` is a shared budget; each stage uses a portion or all of it. -#' * InterProScan can be memory-intensive; on laptops, single-container mode is used internally. +#' * `threads` provides the shared CPU budget for the major processing stages. +#' * Panaroo, CD-HIT, and HMMER allocate that budget according to their respective +#' stage parameters. #' #' @seealso -#' prepareGenomes(), runPanaroo2Duckdb(), CDHIT2duckdb(), domainFromIPR(), cleanData() +#' [prepareGenomes()], [runPanaroo2Duckdb()], [CDHIT2duckdb()], [cleanMetaData()], +#' [cleanData()] #' #' @examples #' \dontrun{ -#' # Paths below are illustrative; adapt to your project layout. #' runDataProcessing( #' duckdb_path = "data/Shigella_flexneri/Sfl.duckdb", #' output_path = "data/Shigella_flexneri", @@ -1905,46 +2723,119 @@ cleanData <- function(duckdb_path, path) { #' ) #' #' # After completion: -#' # data/Shigella_flexneri/Sfl_parquet.duckdb +#' # data/Shigella_flexneri/Sfl_parquet.duckdb #' # will contain views over the Parquet files for downstream ML. #' } #' -#' @export -runDataProcessing <- function(duckdb_path, - output_path = NULL, - # unified threads for all tools - threads = 8, - # Panaroo - panaroo_split_jobs = FALSE, - panaroo_core_threshold = 0.90, - panaroo_len_dif_percent = 0.95, - panaroo_cluster_threshold = 0.95, - panaroo_family_seq_identity = 0.5, - panaroo_refind_mode = c("off", "default", "strict"), - panaroo_strip_pseudogenes = FALSE, - panaroo_pseudogene_clean_dir = "gff_clean", - panaroo_write_pseudogene_audit = TRUE, - # CD-HIT - cdhit_identity = 0.9, - cdhit_word_length = 5, - cdhit_memory = 0, - cdhit_extra_args = c("-g", "1"), - cdhit_output_prefix = "cdhit_out", - # InterPro - ipr_appl = c("Pfam"), - ipr_threads_unused = NULL, - ipr_version = "5.76-107.0", - ipr_dest_dir = "inst/extdata/interpro", - ipr_platform = "linux/amd64", - auto_prepare_data = TRUE, - # Metadata cleaning - ref_file_path = "data_raw/", - verbose = TRUE) { +runDataProcessing <- function( + duckdb_path, + output_path = NULL, + threads = 8, + + # Panaroo + panaroo_split_jobs = FALSE, + panaroo_core_threshold = 0.90, + panaroo_len_dif_percent = 0.95, + panaroo_cluster_threshold = 0.95, + panaroo_family_seq_identity = 0.5, + panaroo_refind_mode = c("off", "default", "strict"), + panaroo_strip_pseudogenes = FALSE, + panaroo_pseudogene_clean_dir = "gff_clean", + panaroo_write_pseudogene_audit = TRUE, + + # CD-HIT + cdhit_identity = 0.9, + cdhit_word_length = 5, + cdhit_memory = 0, + cdhit_extra_args = c("-g", "1"), + cdhit_output_prefix = "cdhit_out", + + # HMMER + hmmer_databases = c( + "Pfam", + "COG", + "AMRFinder", + "DefenseCas" + ), + hmmer_db_dir = NULL, + hmmer_docker_image = "staphb/hmmer", + hmmer_num_splits = 8L, + hmmer_workers = 8L, + + # Metadata cleaning + ref_file_path = "data_raw/", + verbose = TRUE +) { panaroo_refind_mode <- match.arg(panaroo_refind_mode) duckdb_path <- normalizePath(duckdb_path) out_dir <- if (is.null(output_path)) dirname(duckdb_path) else normalizePath(output_path) + # Find the latest manifest + manifest_path <- .manifest_find_latest(duckdb_path) + + if (is.null(manifest_path)) { + stop( + "No provenance manifest found for: ", + duckdb_path, + "\nRun prepareGenomes() first or provide a dataset with an existing manifest." + ) + } + + # Append a new processing run to the existing manifest + manifest <- .manifest_resume( + manifest_path = manifest_path, + base_dir = dirname(dirname(dirname(duckdb_path))), + hash_files = FALSE + ) + + run_failed <- TRUE + + on.exit( + if (run_failed) { + .manifest_finish( + manifest, + status = "failed", + error = "runDataProcessing() exited before successful completion." + ) + }, + add = TRUE + ) + + # Record the start of this processing run + manifest <- .manifest_event( + manifest, + message = "Started data-processing run.", + details = list( + duckdb_path = duckdb_path, + output_path = out_dir + ) + ) + # 1) Panaroo (run + optional merge) -> write Panaroo tables + if (isTRUE(verbose)) message("Running Panaroo and writing gene & struct tables to DuckDB.") + + # Log! + manifest <- .manifest_stage( + manifest, + name = "panaroo", + status = "running", + parameters = list( + core_threshold = panaroo_core_threshold, + len_dif_percent = panaroo_len_dif_percent, + cluster_threshold = panaroo_cluster_threshold, + family_seq_identity = panaroo_family_seq_identity, + threads = threads, + split_jobs = panaroo_split_jobs, + refind_mode = panaroo_refind_mode, + strip_pseudogenes = panaroo_strip_pseudogenes + ), + inputs = duckdb_path, + tool = list( + name = "Panaroo", + docker_image = "staphb/panaroo:1.7.0" + ) + ) + pan_dir <- runPanaroo2Duckdb( duckdb_path = duckdb_path, output_path = out_dir, @@ -1961,8 +2852,56 @@ runDataProcessing <- function(duckdb_path, verbose = verbose ) + manifest <- .manifest_stage( + manifest, + name = "panaroo", + status = "success", + parameters = list( + core_threshold = panaroo_core_threshold, + len_dif_percent = panaroo_len_dif_percent, + cluster_threshold = panaroo_cluster_threshold, + family_seq_identity = panaroo_family_seq_identity, + threads = threads, + split_jobs = panaroo_split_jobs, + refind_mode = panaroo_refind_mode, + strip_pseudogenes = panaroo_strip_pseudogenes + ), + inputs = duckdb_path, + outputs = c( + pan_dir, + duckdb_path + ), + tool = list( + name = "Panaroo", + version = "1.7.0", + docker_image = "staphb/panaroo:1.7.0" + ) + ) + # 2) CD-HIT -> write `protein` tables if (isTRUE(verbose)) message("Running CD-HIT and writing protein tables to DuckDB.") + + # Log! + manifest <- .manifest_stage( + manifest, + name = "cdhit", + status = "running", + parameters = list( + identity = cdhit_identity, + word_length = cdhit_word_length, + memory = cdhit_memory, + threads = threads, + extra_args = cdhit_extra_args, + output_prefix = cdhit_output_prefix + ), + inputs = duckdb_path, + tool = list( + name = "CD-HIT", + version = "4.8.1", + docker_image = "weizhongli1987/cdhit:4.8.1" + ) + ) + CDHIT2duckdb( duckdb_path = duckdb_path, output_path = out_dir, @@ -1974,19 +2913,182 @@ runDataProcessing <- function(duckdb_path, extra_args = cdhit_extra_args ) - # 3) InterProScan -> write `domain` tables - if (isTRUE(verbose)) message("Running InterProScan and writing domain tables to DuckDB.") - domainFromIPR( - duckdb_path = duckdb_path, - path = out_dir, - out_file_base = "iprscan", - appl = ipr_appl, - ipr_version = ipr_version, - ipr_dest_dir = ipr_dest_dir, - ipr_platform = ipr_platform, - auto_prepare_data = auto_prepare_data, - threads = threads, - file_format = "TSV" + manifest <- .manifest_stage( + manifest, + name = "cdhit", + status = "success", + parameters = list( + identity = cdhit_identity, + word_length = cdhit_word_length, + memory = cdhit_memory, + threads = threads, + extra_args = cdhit_extra_args, + output_prefix = cdhit_output_prefix + ), + inputs = duckdb_path, + outputs = c( + file.path(out_dir, paste0(cdhit_output_prefix, "_input.fa")), + file.path(out_dir, paste0(cdhit_output_prefix, "_proteins")), + file.path(duckdb_path) + ), + tool = list( + name = "CD-HIT", + version = "4.8.1", + docker_image = "weizhongli1987/cdhit:4.8.1" + ) + ) + + # 3) HMMER -> write HMM-based match tables for desired databases + if (isTRUE(verbose)) { + message( + "Running HMMER with databases: ", + paste(hmmer_databases, collapse = ", ") + ) + } + + hmmer_db_dir <- if (is.null(hmmer_db_dir)) { + .defaultHmmerDbDir() + } else { + normalizePath(hmmer_db_dir, mustWork = FALSE) + } + + dir.create( + hmmer_db_dir, + recursive = TRUE, + showWarnings = FALSE + ) + + manifest <- .manifest_stage( + manifest, + name = "hmmer", + status = "running", + parameters = list( + databases = hmmer_databases, + database_dir = hmmer_db_dir, + docker_image = hmmer_docker_image, + threads = threads, + num_of_splits = hmmer_num_splits, + workers = hmmer_workers + ), + inputs = duckdb_path, + tool = list( + name = "HMMER", + version = .hmmer_version(hmmer_docker_image), + docker_image = hmmer_docker_image + ) + ) + + generic_databases <- intersect( + hmmer_databases, + c("Pfam", "COG", "AMRFinder") + ) + + if (length(generic_databases)) { + hmmer_result <- .runHMMER( + duckdb_path = duckdb_path, + output_path = out_dir, + threads = threads, + hmmer_db_dir = hmmer_db_dir, + databases = generic_databases, + docker_image = hmmer_docker_image, + num_of_splits = hmmer_num_splits, + n_workers = hmmer_workers, + verbose = verbose + ) + } + + + + if ("DefenseCas" %in% hmmer_databases) { + defense_result <- .defenseHMMER( + defense_db_dir = if (is.null(hmmer_db_dir)) { + .defaultHmmerDbDir() + } else { + file.path(hmmer_db_dir, "DefenseCas") + }, + docker_image = hmmer_docker_image, + duckdb_path = duckdb_path, + output_path = out_dir, + threads = threads, + verbose = verbose + ) + } + + expected_outputs <- file.path( + out_dir, + paste0("protein_", hmmer_databases, ".parquet") + ) + + missing_outputs <- expected_outputs[!file.exists(expected_outputs)] + + if (length(missing_outputs)) { + stop( + "HMMER did not produce all expected outputs:\n", + paste(missing_outputs, collapse = "\n") + ) + } + + con <- DBI::dbConnect( + duckdb::duckdb(), + duckdb_path + ) + on.exit( + DBI::dbDisconnect(con, shutdown = FALSE), + add = TRUE + ) + + expected_tables <- paste0("protein_", hmmer_databases) + + missing_tables <- expected_tables[ + !vapply( + expected_tables, + DBI::dbExistsTable, + logical(1), + conn = con + ) + ] + + if (length(missing_tables)) { + stop( + "HMMER did not produce all expected DuckDB tables:\n", + paste(missing_tables, collapse = "\n") + ) + } + + manifest <- .manifest_stage( + manifest, + name = "hmmer", + status = "success", + parameters = list( + databases = hmmer_databases, + database_dir = hmmer_db_dir, + docker_image = hmmer_docker_image, + threads = threads, + num_of_splits = hmmer_num_splits, + workers = hmmer_workers + ), + inputs = duckdb_path, + outputs = c( + purrr::map( + hmmer_databases, + ~ file.path(out_dir, paste0("protein_", .x, ".parquet")) + ), + duckdb_path + ), + metrics = list( + annotation_tables = paste0( + "protein_", + hmmer_databases + ), + database_provenance = list( + generic = if (!is.null(hmmer_result)) hmmer_result$databases else NULL, + DefenseCas = if (!is.null(defense_result)) defense_result$databases else NULL + ) + ), + tool = list( + name = "HMMER", + docker_image = hmmer_docker_image + ) ) # 4) Clean metadata and export Parquet + Parquet-backed DuckDB @@ -2004,7 +3106,7 @@ runDataProcessing <- function(duckdb_path, if (isTRUE(verbose)) { message("\n============================================") - message("Completed data-processing pipeline successfully.") + message("Completed data-processing workflow successfully.") message("Parquet-backed DuckDB created at:") message(" ", normalizePath(parquet_duckdb_path)) message("\nYou can use the amRml package to train machine") @@ -2014,6 +3116,37 @@ runDataProcessing <- function(duckdb_path, message("============================================\n") } + # Log! + manifest <- .manifest_stage( + manifest, + name = "clean_metadata_and_export", + status = "success", + parameters = list( + reference_path = normalizePath( + ref_file_path, + mustWork = FALSE + ) + ), + inputs = c( + duckdb_path, + ref_file_path + ), + outputs = c( + out_dir, + parquet_duckdb_path + ), + metrics = list( + parquet_duckdb = parquet_duckdb_path + ) + ) + + run_failed <- FALSE + + .manifest_finish( + manifest, + status = "success" + ) + invisible(list( duckdb_path = duckdb_path, panaroo_output = pan_dir, diff --git a/R/helpers.R b/R/helpers.R new file mode 100644 index 0000000..6069462 --- /dev/null +++ b/R/helpers.R @@ -0,0 +1,939 @@ +### Helpers for amRdata live in this script +######################### +# Data curation helpers # +######################### +#' Helps ensure trailing 0s are retained in genome IDs for proper downloading +#' @keywords internal +.id_checker <- function(x) { + # Taxon IDs are just numbers, genome IDs have decimals, this tells them apart + grepl("^[0-9]+$", x) +} + +#' A helper used in data_curation.R and data_processing.R to ensure exported tables +#' don't lose trailing zeroes. Should be relocated into a common helpers/utilities +#' script later. +#' @keywords internal +.preserve_export_id_text <- function(df) { + df <- tibble::as_tibble(df) + + id_pattern <- paste0( + "(^|[._])(", + "genome(_drug)?_id|taxon_id|", + "assembly_accession|bioproject_accession|biosample_accession|", + "refseq_accessions?|genbank_accessions?|sra_accession|pmid|", + "gene_id|protein_id|domain_id|cluster_id|AccNum|id", + ")$" + ) + + id_cols <- names(df)[grepl(id_pattern, names(df), ignore.case = TRUE)] + if (length(id_cols)) { + df[id_cols] <- lapply(df[id_cols], as.character) + } + + df +} + +#' Helps normalize Docker paths +#' @keywords internal +.docker_path <- function(p) gsub("\\\\", "/", normalizePath(p, mustWork = FALSE)) + +#' Helps run a shell inside a container, and prefers bash (don't we all?) +#' @keywords internal +.pick_shell <- function(image) { + chk <- suppressWarnings(system2("docker", + c( + "run", "--rm", image, "sh", "-lc", + "command -v bash >/dev/null || echo NOBASH" + ), + stdout = TRUE, stderr = TRUE + )) + if (length(chk) && any(grepl("NOBASH", chk))) "sh" else "bash" +} + +# FASTA sanitizer to ensure Panaroo compatibility with BV-BRC CLI downloads +.strip_fasta_preamble <- function(fna_path) { + if (!file.exists(fna_path)) { + return(invisible(FALSE)) + } + txt <- readLines(fna_path, warn = FALSE) + first <- which(grepl("^\\s*>", txt))[1] + if (is.na(first)) { + return(invisible(FALSE)) + } + if (first > 1L) { + txt <- txt[first:length(txt)] + txt[1] <- sub("^\\ufeff", "", txt[1]) + writeLines(txt, fna_path, sep = "\n", useBytes = TRUE) + return(invisible(TRUE)) + } + invisible(FALSE) +} + +# GFF sanitizer to ensure Panaroo compatibility with BV-BRC CLI downloads +.sanitize_gff <- function(gff_path) { + if (!file.exists(gff_path)) { + return(invisible(FALSE)) + } + lines <- readLines(gff_path, warn = FALSE) + if (length(lines) == 0L) { + return(invisible(FALSE)) + } + if (!grepl("^##gff-version\\s*3", lines[1])) { + lines <- c("##gff-version 3", lines) + } + out <- purrr::map_chr(lines, function(line) { + if (grepl("^#", line)) { + return(line) + } + parts <- strsplit(line, "[\t ]", perl = TRUE)[[1]] + if (length(parts) >= 9) { + paste(c(parts[1:8], paste(parts[9:length(parts)], collapse = " ")), collapse = "\t") + } else { + line + } + }) + writeLines(out, gff_path, sep = "\n", useBytes = TRUE) + invisible(TRUE) +} + +######################### +# Manifest helpers # +######################### + +#' Returns the basics about a file for manifest logging +#' +#' @param path Character vector of file paths. +#' @param hash Logical. If TRUE, calculate SHA-256 checksums. +#' +#' @return A list of file records. +#' @keywords internal +.manifest_file_info <- function(path, hash = FALSE) { + path <- unique(as.character(path)) + path <- path[nzchar(path)] + + if (!length(path)) { + return(list()) + } + + # See what exists + purrr::map(path, function(x) { + exists <- file.exists(x) + + out <- list( + path = x, + exists = exists, + size_bytes = if (exists) file.info(x)$size else NA_real_, + modified_at = if (exists) as.character(file.info(x)$mtime) else NA_character_ + ) + + # Hash what exists, if desired + if (isTRUE(hash) && exists && !dir.exists(x)) { + out$sha256 <- unname(tools::sha256(x)) + } + + out + }) +} + + +#' Capture basic GitHub repo state for manifest provenance +#' +#' @param base_dir Character. Project root. +#' +#' @return A named list. +#' @keywords internal +.manifest_git_info <- function(base_dir = ".") { + base_dir <- normalizePath(base_dir, mustWork = FALSE) + + # Find Git + git <- Sys.which("git") + + if (!nzchar(git)) { + return(list( + available = FALSE + )) + } + + # Run Git through system commands + run_git <- function(args) { + tryCatch( + system2( + git, + args = args, + stdout = TRUE, + stderr = FALSE + ), + error = function(e) character() + ) + } + + inside <- run_git(c("-C", shQuote(base_dir), "rev-parse", "--is-inside-work-tree")) + + if (!length(inside) || !identical(trimws(inside[[1]]), "true")) { + return(list( + available = TRUE, + repository = FALSE + )) + } + + commit <- run_git(c("-C", shQuote(base_dir), "rev-parse", "HEAD")) + branch <- run_git(c("-C", shQuote(base_dir), "rev-parse", "--abbrev-ref", "HEAD")) + dirty <- run_git(c("-C", shQuote(base_dir), "status", "--porcelain")) + + list( + available = TRUE, + repository = TRUE, + commit = if (length(commit)) trimws(commit[[1]]) else NA_character_, + branch = if (length(branch)) trimws(branch[[1]]) else NA_character_, + dirty = length(dirty) > 0L + ) +} + + +#' Capture package versions currently loaded in the R session +#' +#' @return Named character vector of package versions. +#' @keywords internal +.manifest_package_versions <- function() { + pkgs <- sort(loadedNamespaces()) + + stats::setNames( + as.list( + purrr::map_chr( + pkgs, + function(pkg) { + tryCatch( + as.character(utils::packageVersion(pkg)), + error = function(e) NA_character_ + ) + } + ) + ), + pkgs + ) +} + + +#' Generate a unique manifest run identifier +#' +#' @return Character scalar. +#' @keywords internal +.manifest_run_id <- function() { + paste0( + "run_", + format(Sys.time(), "%Y%m%dT%H%M%OS3", tz = "UTC"), + "_pid", + Sys.getpid() + ) |> + gsub("[^A-Za-z0-9_]", "", x = _) +} + + +#' Start or load a dataset provenance manifest +#' +#' @param manifest_path Character. Path to the JSON manifest. +#' @param dataset_id Character scalar. +#' @param duckdb_path Character scalar. +#' @param base_dir Character scalar. +#' @param selection Optional named list describing the dataset selection. +#' @param hash_files Logical. Calculate SHA-256 for manifest-recorded files. +#' +#' @return A manifest object with `path` and `run_index`. +#' @keywords internal +.manifest_start <- function( + manifest_path, + dataset_id, + duckdb_path, + base_dir = ".", + selection = list(), + hash_files = FALSE +) { + if (!requireNamespace("jsonlite", quietly = TRUE)) { + stop("Package 'jsonlite' is required for manifest generation.") + } + + manifest_path <- normalizePath( + manifest_path, + mustWork = FALSE + ) + + dir.create( + dirname(manifest_path), + recursive = TRUE, + showWarnings = FALSE + ) + + manifest <- list( + schema_version = 1L, + manifest_created_at = as.character(Sys.time()), + manifest_updated_at = as.character(Sys.time()), + dataset_id = dataset_id, + dataset = list( + duckdb = duckdb_path, + selection = selection + ), + runs = list() + ) + + run <- list( + run_id = .manifest_run_id(), + status = "running", + started_at = as.character(Sys.time()), + finished_at = NA_character_, + command = commandArgs(trailingOnly = FALSE), + working_directory = getwd(), + host = as.list(Sys.info()), + r = list( + version = R.version.string, + platform = R.version$platform + ), + git = .manifest_git_info(base_dir), + packages = .manifest_package_versions(), + stages = list(), + events = list() + ) + + if (is.null(manifest$runs)) { + manifest$runs <- list() + } + + manifest$runs[[length(manifest$runs) + 1L]] <- run + manifest$manifest_updated_at <- as.character(Sys.time()) + + run_index <- length(manifest$runs) + + jsonlite::write_json( + manifest, + manifest_path, + auto_unbox = TRUE, + pretty = TRUE, + null = "null" + ) + + structure( + list( + manifest = manifest, + path = manifest_path, + run_index = run_index, + hash_files = isTRUE(hash_files) + ), + class = "amr_manifest" + ) +} + + +#' Update a manifest stage +#' +#' @param manifest_state Manifest state returned by [.manifest_start()]. +#' @param name Character stage name. +#' @param status Character stage status. +#' @param parameters Optional named list. +#' @param inputs Optional character vector of input paths. +#' @param outputs Optional character vector of output paths. +#' @param tool Optional named list describing the tool. +#' @param metrics Optional named list of metrics. +#' @param message Optional log message. +#' +#' @return Updated manifest state. +#' @keywords internal +.manifest_stage <- function( + manifest_state, + name, + status = "success", + parameters = list(), + inputs = character(), + outputs = character(), + tool = list(), + metrics = list(), + message = NULL +) { + if (!inherits(manifest_state, "amr_manifest")) { + stop("Invalid manifest state.") + } + + stage_index <- which( + purrr::map_lgl( + manifest_state$manifest$runs[[manifest_state$run_index]]$stages, + ~ identical(.x$name, name) && identical(.x$status, "running") + ) + ) + + stage <- list( + name = name, + status = status, + started_at = as.character(Sys.time()), + parameters = parameters, + inputs = .manifest_file_info(inputs, hash = manifest_state$hash_files), + outputs = .manifest_file_info(outputs, hash = manifest_state$hash_files), + tool = tool, + metrics = metrics + ) + + if (!is.null(message)) { + stage$message <- as.character(message) + } + + if (length(stage_index) == 1L) { + existing <- manifest_state$manifest$runs[[manifest_state$run_index]]$stages[[stage_index]] + + stage$started_at <- existing$started_at + stage$finished_at <- if (status != "running") { + as.character(Sys.time()) + } else { + NULL + } + + manifest_state$manifest$runs[[manifest_state$run_index]]$stages[[stage_index]] <- stage + } else { + if (status != "running") { + stage$finished_at <- as.character(Sys.time()) + } + + manifest_state$manifest$runs[[manifest_state$run_index]]$stages <- + append( + manifest_state$manifest$runs[[manifest_state$run_index]]$stages, + list(stage) + ) + } + + manifest_state$manifest$manifest_updated_at <- as.character(Sys.time()) + + jsonlite::write_json( + manifest_state$manifest, + manifest_state$path, + auto_unbox = TRUE, + pretty = TRUE, + null = "null" + ) + + manifest_state +} + + +#' Append a provenance event to the active manifest run +#' +#' @param manifest_state Manifest state returned by [.manifest_start()]. +#' @param level Character event level. +#' @param message Character message. +#' @param details Optional named list. +#' +#' @return Updated manifest state. +#' @keywords internal +.manifest_event <- function( + manifest_state, + level = "info", + message, + details = list() +) { + manifest_state$manifest$runs[[manifest_state$run_index]]$events <- + append( + manifest_state$manifest$runs[[manifest_state$run_index]]$events, + list( + list( + timestamp = as.character(Sys.time()), + level = level, + message = message, + details = details + ) + ) + ) + + manifest_state$manifest$manifest_updated_at <- as.character(Sys.time()) + + jsonlite::write_json( + manifest_state$manifest, + manifest_state$path, + auto_unbox = TRUE, + pretty = TRUE, + null = "null" + ) + + manifest_state +} + + +#' Finish an active provenance manifest run +#' +#' @param manifest_state Manifest state returned by [.manifest_start()]. +#' @param status Final run status. +#' @param error Optional error message. +#' +#' @return Invisibly returns the final manifest state. +#' @keywords internal +.manifest_finish <- function( + manifest_state, + status = "success", + error = NULL +) { + manifest_state$manifest$runs[[manifest_state$run_index]]$status <- status + manifest_state$manifest$runs[[manifest_state$run_index]]$finished_at <- + as.character(Sys.time()) + + if (!is.null(error)) { + manifest_state$manifest$runs[[manifest_state$run_index]]$error <- as.character(error) + } + + manifest_state$manifest$manifest_updated_at <- as.character(Sys.time()) + + jsonlite::write_json( + manifest_state$manifest, + manifest_state$path, + auto_unbox = TRUE, + pretty = TRUE, + null = "null" + ) + + invisible(manifest_state) +} + +# To distinguish multiple manifests in the same bug directory +.manifest_find_latest <- function(duckdb_path) { + manifest_dir <- dirname(normalizePath( + duckdb_path, + mustWork = FALSE + )) + + manifests <- list.files( + manifest_dir, + pattern = "^manifest_.*\\.json$", + full.names = TRUE + ) + + if (!length(manifests)) { + return(NULL) + } + + manifests[which.max(file.info(manifests)$mtime)] +} + +#' Resume provenance logging in an existing manifest +#' +#' Loads an existing manifest and appends a new run. +#' +#' @param manifest_path Character. Path to an existing JSON manifest. +#' @param base_dir Character. Project root. +#' @param hash_files Logical. Calculate SHA-256 checksums for manifest-recorded files. +#' +#' @return A manifest object with `path` and `run_index`. +#' @keywords internal +.manifest_resume <- function( + manifest_path, + base_dir = ".", + hash_files = FALSE +) { + if (!requireNamespace("jsonlite", quietly = TRUE)) { + stop("Package 'jsonlite' is required for manifest generation.") + } + + manifest_path <- normalizePath( + manifest_path, + mustWork = TRUE + ) + + manifest <- jsonlite::read_json( + manifest_path, + simplifyVector = FALSE + ) + + if (is.null(manifest$runs)) { + manifest$runs <- list() + } + + run <- list( + run_id = .manifest_run_id(), + status = "running", + started_at = as.character(Sys.time()), + finished_at = NA_character_, + command = commandArgs(trailingOnly = FALSE), + working_directory = getwd(), + host = as.list(Sys.info()), + r = list( + version = R.version.string, + platform = R.version$platform + ), + git = .manifest_git_info(base_dir), + packages = .manifest_package_versions(), + stages = list(), + events = list() + ) + + manifest$runs[[length(manifest$runs) + 1L]] <- run + manifest$manifest_updated_at <- as.character(Sys.time()) + + run_index <- length(manifest$runs) + + jsonlite::write_json( + manifest, + manifest_path, + auto_unbox = TRUE, + pretty = TRUE, + null = "null" + ) + + structure( + list( + manifest = manifest, + path = manifest_path, + run_index = run_index, + hash_files = isTRUE(hash_files) + ), + class = "amr_manifest" + ) +} + +########################### +# Data processing helpers # +########################### + +# Map host paths under mounted root to container path +#' .to_container() +#' +#' Used for OS-agnostic mapping of Docker directories and mount paths +#' +#' @keywords internal +#' @examples NULL +.to_container <- function(x, host_root, container_root = "/work") { + host_root_unix <- .docker_path(host_root) + x_unix <- .docker_path(x) + pattern <- paste0("^", gsub("([\\^$.|?*+(){}\\[\\]\\\\])", "\\\\\\\\\\1", host_root_unix)) + sub(pattern, container_root, x_unix) +} + +#' Remove pseudogene annotations from Panaroo input GFF files +#' +#' Cleans GFF annotation files of `pseudogene` feature records only. Cleaned +#' GFFs are written to a subdirectory under `output_path` and swapped into the +#' Panaroo input list, leaving the original genome annotations alone. +#' +#' This optional preprocessing step can reduce weird runtime stalls during +#' Panaroo graph construction for some BV-BRC/PATRIC genome annotations that +#' contain troublesome pseudogene features. +#' +#' @param panaroo_input_files Character vector of `"gff fna"` input lines used +#' by Panaroo. +#' @param output_path Character scalar. Base directory for temporary cleaned +#' GFF files and audit outputs. +#' @param clean_dir Character scalar. Name of the subdirectory created beneath +#' `output_path` to store cleaned GFF files. Default `"gff_clean"`. +#' +#' @return A list containing: +#' \itemize{ +#' \item `panaroo_input_files` — rewritten Panaroo input lines pointing to the +#' cleaned GFF files. +#' \item `audit` — a tibble summarizing, for each genome, the total number of +#' annotated features, the number of pseudogenes removed, and the number of +#' remaining features. +#' } +#' +#' @details +#' This performs lightweight preprocessing only, removing feature records whose +#' third GFF column is exactly `"pseudogene"` and does not otherwise modify +#' annotation coordinates, attributes, or sequence files. FASTA paths are unchanged. +#' +#' @keywords internal +.stripPseudogeneGFFs <- function(panaroo_input_files, + output_path, + clean_dir = "gff_clean") { + # Normalize our paths + panaroo_input_files <- as.character(panaroo_input_files) + output_path <- .docker_path(output_path) + + # Set the directory to place cleaned GFFs into + clean_root <- file.path(output_path, clean_dir) + dir.create(clean_root, recursive = TRUE, showWarnings = FALSE) + + # Where the cleaned up Panaroo input and clean audit is stored + out_lines <- character(length(panaroo_input_files)) + audit <- vector("list", length(panaroo_input_files)) + + for (i in seq_along(panaroo_input_files)) { + + # Read the Panaroo gff + fna input lines + line <- panaroo_input_files[[i]] + parts <- strsplit(line, "\\s+")[[1]] + + # If you're missing either a gff or an fna file in there somehow + if (length(parts) < 2L) { + stop("Broken Panaroo input line: ", line) + } + + # Read in the files parsed above + gff_in <- .docker_path(parts[1]) + fna_in <- .docker_path(parts[2]) + + # If it didn't read in + if (!file.exists(gff_in)) { + stop("Missing GFF file: ", gff_in) + } + if (!file.exists(fna_in)) { + stop("Missing FNA file: ", fna_in) + } + + # What we're saving out + gff_out <- file.path(clean_root, basename(gff_in)) + + # Read in the GFF lines and fine the comment lined headers + gff_lines <- readLines(gff_in, warn = FALSE) + is_header <- startsWith(gff_lines, "#") + body <- gff_lines[!is_header] + + # If there's nothing in there to parse + if (length(body) == 0L) { + writeLines(gff_lines, gff_out, useBytes = TRUE) + n_total <- 0L + n_pseudogene <- 0L + n_kept <- 0L + } else { + # Otherwise, find the pseudogene lines and save everything but those + # Now with 100% more purrr + fields <- strsplit(body, "\t", fixed = TRUE) + types <- purrr::map_chr( + fields, + \(x) if (length(x) >= 3L) x[[3]] else NA_character_ + ) + keep <- !is.na(types) & types != "pseudogene" + + cleaned <- c(gff_lines[is_header], body[keep]) + writeLines(cleaned, gff_out, useBytes = TRUE) + + # We love stats + n_total <- length(body) + n_pseudogene <- sum(!keep, na.rm = TRUE) + n_kept <- sum(keep, na.rm = TRUE) + } + + # Record what we did and save it into the audit log + audit[[i]] <- tibble::tibble( + gff_in = gff_in, + gff_out = gff_out, + n_total_features = n_total, + n_pseudogene = n_pseudogene, + n_kept = n_kept + ) + + out_lines[[i]] <- paste(gff_out, fna_in) + } + + list( + panaroo_input_files = out_lines, + audit = dplyr::bind_rows(audit) + ) +} + +######################### +# HMMER helpers # +######################### + +#' Validate if a HMM file has old HMMER3 format +#' +#' @param hmm_file +#' +#' @returns +#' +#' @keywords internal +#' @examples +.isValidHmmFile <- function(hmm_file) { + + lines <- tryCatch( + readLines(hmm_file, warn = FALSE), + error = function(e) character(0) + ) + + if (length(lines) == 0) { + return(FALSE) + } + + first_line <- trimws(lines[1]) + last_line <- trimws(tail(lines, 1)) + + starts_ok <- grepl("^HMMER3/f", first_line) + ends_ok <- identical(last_line, "//") + + starts_ok && ends_ok +} + +#' Parsing HMM database to extract profile names, accessions and descriptions +#' +#' @param hmm_file path to the HMM database file (`.hmm`) +#' +#' @returns a tibble +#' +#' @keywords internal +.parse_hmmer_profiles <- function(hmm_file) { + + lines <- readLines(hmm_file, warn = FALSE) + + starts <- c( + which(grepl("^NAME\\s+", lines)), + length(lines) + 1L + ) + + blocks <- purrr::map2( + starts[-length(starts)], + starts[-1L] - 1L, + ~ lines[.x:.y] + ) + + extract_field <- function(block, pattern) { + + hit <- stringr::str_subset(block, pattern) + + if (length(hit) == 0) { + return(NA_character_) + } + + stringr::str_remove(hit[[1]], pattern) + } + + purrr::map_dfr( + blocks, + ~ tibble::tibble( + profile_name = extract_field(.x, "^NAME\\s+"), + profile_accession = extract_field(.x, "^ACC\\s+"), + profile_description = extract_field(.x, "^DESC\\s+") + ) + ) +} + +#' Parse HMMER tabular output into a tibble +#' +#' Reads a HMMER `--domtblout` file and returns a tidy tibble with one row per +#' target-query hit. Comment lines are stripped and the free-text description +#' field is reunited from the remaining whitespace-delimited columns. +#' +#' @param file Path to a HMMER `.tbl` output file produced with `--domtblout`. +#' +#' @return A tibble with 19 columns matching the HMMER per-sequence hit table +#' +#' @references Adapted from the rhmmer package +#' (). +#' +#' @examples +#' \dontrun{ +#' hits <- .parseHMMEROutput("results/Ecoli/protein_chunk_01_COG.tbl") +#' hits |> dplyr::filter(sequence_evalue < 1e-5) +#' } +#' +#' @keywords internal +.parseHMMEROutput <- function(file) { + + # target name accession tlen query name accession qlen E-value score bias # of c-Evalue i-Evalue score bias from to from to from to acc description of target + col_types <- readr::cols( + protein = readr::col_character(), # target name + protein_accession = readr::col_character(), + tlen = readr::col_integer(), + + query_name = readr::col_character(), # query name + query_accession = readr::col_character(), + qlen = readr::col_integer(), + + sequence_evalue = readr::col_double(), + sequence_score = readr::col_double(), + sequence_bias = readr::col_double(), + + domain_num = readr::col_integer(), + domain_of = readr::col_integer(), + + c_evalue = readr::col_double(), + i_evalue = readr::col_double(), + + domain_score = readr::col_double(), + domain_bias = readr::col_double(), + + hmm_from = readr::col_integer(), + hmm_to = readr::col_integer(), + + ali_from = readr::col_integer(), + ali_to = readr::col_integer(), + + env_from = readr::col_integer(), + env_to = readr::col_integer(), + + acc = readr::col_double(), + + target_description = readr::col_character() + ) + # the line delimiter should always be just "\n", even on Windows + lines <- readr::read_lines(file, lazy = FALSE, progress = FALSE) + + # drop comment lines + data_lines <- lines[!grepl("^#", lines)] + + # split: whitespace-separated fields + split_fields <- strsplit(data_lines, "\\s+", perl = TRUE) + + # count space separated fields + N <- max(sapply(split_fields, length)) + + # Parsing differently to avoid fussy read_tsv() warnings + txt <- sub( + pattern = sprintf("(%s).*", paste0(rep("\\S+", N), collapse = " +")), + replacement = "\\1", + x = lines, + perl = TRUE + ) |> + gsub(pattern = " *", replacement = "\t") |> + paste0(collapse = "\n") + + table <- readr::read_tsv( + I(txt), + col_names = names(col_types$cols), + comment = "#", + na = "-", + col_types = col_types, + lazy = FALSE, + progress = FALSE + ) + + table +} + +#' Write a data frame to a compressed Parquet file +#' +#' @param df A data frame or tibble to write. +#' @param path Output file path (`.parquet` extension). +#' +#' @keywords internal +.write_compressed_parquet <- function(df, path) { + arrow::write_parquet( + df, + path, + compression = "zstd", + compression_level = 9, + use_dictionary = TRUE + ) +} + +# Default persistent cache for shared HMMER databases +.defaultHmmerDbDir <- function() { + file.path( + tools::R_user_dir("amRdata", "cache"), + "hmmer" + ) +} + +.hmmer_version <- function(docker_image = "staphb/hmmer") { + output <- system2( + "docker", + args = c( + "run", + "--rm", + docker_image, + "hmmsearch", + "-h" + ), + stdout = TRUE, + stderr = TRUE + ) + + version <- stringr::str_match( + paste(output, collapse = "\n"), + "HMMER ([0-9.]+)" + )[, 2] + + if (is.na(version)) { + return(NA_character_) + } + + version +} diff --git a/R/runHMMER.R b/R/runHMMER.R deleted file mode 100644 index a463c8f..0000000 --- a/R/runHMMER.R +++ /dev/null @@ -1,1352 +0,0 @@ -#' Validate if a HMM file has old HMMER3 format and remove them -#' -#' @param hmm_file -#' -#' @returns -#' -#' @keywords internal -#' @examples -.isValidHmmFile <- function(hmm_file) { - - lines <- tryCatch( - readLines(hmm_file, warn = FALSE), - error = function(e) character(0) - ) - - if (length(lines) == 0) { - return(FALSE) - } - - first_line <- trimws(lines[1]) - last_line <- trimws(tail(lines, 1)) - - starts_ok <- grepl("^HMMER3/f", first_line) - ends_ok <- identical(last_line, "//") - - starts_ok && ends_ok -} - -#' 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. -#' -#' @returns A list of paths to the database hmm files. -#' -#' @keywords internal -#' @examples -.prepareHmmerDatabases <- function( - hmmer_db_dir, - databases = c("Pfam", "COG", "AMRFinder"), - docker_image = "staphb/hmmer", - hmmer_db_url = NULL, - verbose = TRUE -) { - - hmmer_db_dir <- normalizePath(hmmer_db_dir) - - options(timeout = max(3600, getOption("timeout"))) - - dbs <- list( - - Pfam = list( - dir = file.path(hmmer_db_dir, "Pfam"), - hmm_name = "Pfam-A.hmm", - url = "https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz", - type = "gz" - ), - - COG = list( - dir = file.path(hmmer_db_dir, "COG"), - hmm_name = "COG_database2024.hmm", - url = "http://boabio.belozersky.msu.ru/media/COG_database2024.zip", - type = "zip" - ), - - AMRFinder = list( - dir = file.path(hmmer_db_dir, "AMRFinder"), - hmm_name = NULL, - url = "https://ftp.ncbi.nlm.nih.gov/hmm/NCBIfam-AMRFinder/latest/NCBIfam-AMRFinder.HMM.tar.gz", - type = "tar.gz" - ) - ) - - ## Add custom database(s) - missing_dbs <- setdiff(databases, names(dbs)) - - if (length(missing_dbs) > 0) { - - if (is.null(hmmer_db_url)) { - stop( - "hmmer_db_url must be supplied when using custom databases" - ) - } - - get_db_type <- function(url) { - - file <- basename(url) - - if (grepl("\\.(tar\\.gz|tgz)$", - file, - ignore.case = TRUE)) { - - return("tar.gz") - - } else if (grepl("\\.zip$", - file, - ignore.case = TRUE)) { - - return("zip") - - } else if (grepl("\\.gz$", - file, - ignore.case = TRUE)) { - - return("gz") - - } else { - - stop( - "Unsupported archive type: ", - file - ) - } - } - - for (db_name in missing_dbs) { - - dbs[[db_name]] <- list( - dir = file.path(hmmer_db_dir, db_name), - hmm_name = NULL, - url = hmmer_db_url, - type = get_db_type(hmmer_db_url) - ) - } - } - - dbs <- dbs[databases] - - db_paths <- list() - - for (db_name in names(dbs)) { - - db <- dbs[[db_name]] - - dir.create( - db$dir, - recursive = TRUE, - showWarnings = FALSE - ) - - if(verbose) message("Checking ", db_name) - - hmm_files <- list.files( - db$dir, - pattern = "\\.hmm$", - recursive = TRUE, - full.names = TRUE, - ignore.case = TRUE - ) - - if (length(hmm_files) == 0) { - - message("Downloading ", db_name) - - tmp <- tempfile() - - utils::download.file( - url = db$url, - destfile = tmp, - mode = "wb", - method = "libcurl" - ) - - switch( - - db$type, - - gz = { - - hmm_file <- file.path( - db$dir, - db$hmm_name %||% basename( - sub("\\.gz$", - "", - basename(db$url), - ignore.case = TRUE) - ) - ) - - R.utils::gunzip( - filename = tmp, - destname = hmm_file, - overwrite = TRUE, - remove = FALSE - ) - }, - - zip = { - - utils::unzip( - zipfile = tmp, - exdir = db$dir - ) - }, - - `tar.gz` = { - - utils::untar( - tarfile = tmp, - exdir = db$dir - ) - } - ) - - hmm_files <- list.files( - db$dir, - pattern = "\\.hmm$", - recursive = TRUE, - full.names = TRUE, - ignore.case = TRUE - ) - } - - hmm_files <- list.files( - db$dir, - pattern = "\\.hmm$", - recursive = TRUE, - full.names = TRUE, - ignore.case = TRUE - ) - - if (length(hmm_files) == 0) { - - stop( - "No .hmm file found for ", - db_name - ) - - } else if (length(hmm_files) == 1) { - - hmm_file <- hmm_files[1] - - } else { - - hmm_file <- file.path( - db$dir, - paste0(db_name, ".hmm") - ) - - source_hmms <- setdiff( - normalizePath(hmm_files), - normalizePath(hmm_file, mustWork = FALSE) -) - - # purrr implementation - valid_hmms <- purrr::map_lgl(source_hmms, .isValidHmmFile) - -if (any(!valid_hmms)) { - bad_files <- basename(hmm_files[!valid_hmms]) - - if (isTRUE(verbose)) { - warning( - "Ignoring ", - length(bad_files), - " invalid HMM file(s):\n", - paste(bad_files, collapse = "\n"), - call. = FALSE - ) - } - - hmm_files <- hmm_files[valid_hmms] -} - -if (length(source_hmms) == 0) { - - stop( - "No valid HMM files found for ", - db_name - ) -} - - if (!file.exists(hmm_file)) { - - if(verbose) message( - "Combining ", - length(source_hmms), - " HMM files for ", - db_name - ) - - file.create(hmm_file) - - for (f in sort(source_hmms)) { - file.append(hmm_file, f) - } - } - } - - pressed_files <- paste0( - hmm_file, - c(".h3m", ".h3i", ".h3f", ".h3p") - ) - - if (!all(file.exists(pressed_files))) { - - if(verbose) message( - "Running hmmpress for ", - basename(hmm_file) - ) - - output <- system2( - "docker", - args = c( - "run", - "--rm", - "-v", - paste0(dirname(hmm_file), ":/db"), - docker_image, - "hmmpress", - file.path("/db", basename(hmm_file)) - ), - stdout = TRUE, - stderr = TRUE - ) - - if (!all(file.exists(pressed_files))) { - - stop( - "hmmpress failed for ", - db_name, - "\n", - paste(output, collapse = "\n") - ) - } - } - - db_paths[[db_name]] <- hmm_file - - message( - db_name, - " ready: ", - hmm_file - ) - } - - db_paths -} - -#' Write a data frame to a compressed Parquet file -#' -#' @param df A data frame or tibble to write. -#' @param path Output file path (`.parquet` extension). -#' -#' @keywords internal -.write_compressed_parquet <- function(df, path) { - arrow::write_parquet( - df, - path, - compression = "zstd", - compression_level = 9, - use_dictionary = TRUE - ) -} - -#' Parsing HMM database to extract profile names, accessions and descriptions -#' -#' @param hmm_file path to the HMM database file (`.hmm`) -#' -#' @returns a tibble -#' -#' @keywords internal -.parse_hmmer_profiles <- function(hmm_file) { - - lines <- readLines(hmm_file, warn = FALSE) - - starts <- c( - which(grepl("^NAME\\s+", lines)), - length(lines) + 1L - ) - - blocks <- purrr::map2( - starts[-length(starts)], - starts[-1L] - 1L, - ~ lines[.x:.y] - ) - - extract_field <- function(block, pattern) { - - hit <- stringr::str_subset(block, pattern) - - if (length(hit) == 0) { - return(NA_character_) - } - - stringr::str_remove(hit[[1]], pattern) - } - - purrr::map_dfr( - blocks, - ~ tibble::tibble( - profile_name = extract_field(.x, "^NAME\\s+"), - profile_accession = extract_field(.x, "^ACC\\s+"), - profile_description = extract_field(.x, "^DESC\\s+") - ) - ) -} - -#' The function to run HMMER with docker -#' -#' @param JOB_NAME -#' @param FASTA -#' @param DB -#' @param Total_proteins -#' @param output_path -#' @param db_paths -#' @param docker_image -#' @param threads -#' @param n_workers -#' -#' @returns -#' -#' @keywords internal -.runHmmerJob <- function(JOB_NAME, FASTA, DB, Total_proteins, - output_path = NULL, db_paths, - docker_image = "staphb/hmmer", threads = 8L, - n_workers = 8L, - verbose = TRUE -) { - hmmer_input <- file.path(output_path, FASTA) - hmmer_output <- file.path(output_path, paste0(JOB_NAME, ".tbl")) - - # database paths - database_path <- db_paths[[DB]] - db_host_dir <- dirname(database_path) - db_filename <- basename(database_path) - db_cont_dir <- "/opt/hmmer/data" - db_cont_path <- file.path(db_cont_dir, db_filename) - - # mounts - mount_host <- output_path - mount_cont <- "/work" - - threads_per_job <- max( - 1L, - floor(threads / n_workers) -) - - cmd_args <- c( - "run", "--rm", - "-v", paste0(mount_host, ":", mount_cont), - "-v", paste0(db_host_dir, ":", db_cont_dir), - docker_image, - "hmmsearch", - "--notextw", - "--cpu", as.character(threads_per_job), - "-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) - ) - - if(verbose) message("Running hmmsearch via Docker...") - output <- tryCatch( - { - system2("docker", args = cmd_args, stdout = TRUE, stderr = TRUE) - }, - error = function(e) { - stop("hmmsearch execution failed: ", e$message) - } - ) - - if (!file.exists(hmmer_output)) { - stop("hmmsearch failed: output file not found. Check stderr:\n", paste(output, collapse = "\n")) - } - - if(verbose) message("hmmsearch completed successfully.") - - hmmer_tbl <- .parseHMMEROutput(hmmer_output) |> - dplyr::select("protein", "query_name") - - hmmer_tbl_filename <- file.path( - dirname(hmmer_output), - paste0(tools::file_path_sans_ext(basename(hmmer_output)), ".parquet") - ) - - .write_compressed_parquet(hmmer_tbl, hmmer_tbl_filename) - - hmmer_tbl_filename - } - - -#' Wrapper for preparing HMM databases and running HMMER on protein sequences from duckdb and writing them. -#' -#' @param duckdb_path -#' @param output_path -#' @param threads -#' @param hmmer_db_dir -#' @param databases -#' @param docker_image -#' @param num_of_splits -#' @param n_workers -#' -#' @returns -#' -#' @keywords internal -#' @examples -.runHMMER <- function(duckdb_path, - output_path, - threads = 8L, - hmmer_db_dir, - databases = c("Pfam", "COG", "AMRFinder"), - docker_image = "staphb/hmmer", - num_of_splits = 8L, - n_workers = 8L, - verbose = TRUE - ) { - # Fail fast if Docker is missing - if (!nzchar(Sys.which("docker"))) { - stop("Docker is not available on your PATH but is required to run HMMER.") - } - - # But also check if Docker is on the PATH but isn't running - docker_ok <- system2( - "docker", - "info", - stdout = FALSE, - stderr = FALSE - ) == 0L - - if (!docker_ok) { - stop( - "Docker is installed but is not running or cannot be reached. ", - "Please (re)start Docker Desktop and try again." - ) - } - - duckdb_path <- .docker_path(duckdb_path) - if (missing(output_path) || output_path %in% c(".", "results", "results/")) { - output_path <- dirname(duckdb_path) - } - output_path <- .docker_path(output_path) - if (!dir.exists(output_path)) dir.create(output_path, recursive = TRUE) - - con <- DBI::dbConnect(duckdb::duckdb(), duckdb_path) - on.exit(try(DBI::dbDisconnect(con, shutdown = FALSE), silent = TRUE), add = TRUE) - - prot_seqs <- DBI::dbReadTable(con, "protein_cluster_seq") |> - tibble::as_tibble() - - # required to define the database size for hmmsearch --Z and --domZ parameters - Total_proteins <- nrow(prot_seqs) - - if(is.null(hmmer_db_dir)) { - hmmer_db_dir <- output_path - } - - # database paths - if(verbose) message ("Preparing HMM databases") - db_paths <- .prepareHmmerDatabases( - hmmer_db_dir = hmmer_db_dir, - databases = databases, - docker_image = docker_image, - verbose = verbose -) - -db_paths <- db_paths[databases] - - # clamp splits to the number of sequences available - chunk_count <- min(as.integer(num_of_splits), nrow(prot_seqs)) - - split_fasta <- function(seqs, prefix) { - records <- paste0(">", seqs$name, "\n", seqs$sequence) - chunk_size <- ceiling(length(records) / chunk_count) - chunks <- split(records, ceiling(seq_along(records) / chunk_size)) - - purrr::walk2(chunks, seq_along(chunks), function(chunk, i) { - chunk_path <- file.path(output_path, sprintf("%s_chunk_%02d.fasta", prefix, i)) - readr::write_lines(chunk, chunk_path) - }) - } - - split_fasta(prot_seqs, "protein") - - job_list <- expand.grid( - chunk = sprintf("%02d", seq_len(chunk_count)), - db = databases, - stringsAsFactors = FALSE - ) |> - dplyr::mutate( - JOB_NAME = paste0("protein_chunk_", chunk, "_", db), - FASTA = paste0("protein_chunk_", chunk, ".fasta"), - DB = db - ) |> - dplyr::select(JOB_NAME, FASTA, DB) - - future::plan( - future::multisession, - workers = max(1L, n_workers) -) - - if (verbose) message("Running HMMER jobs") -parquet_files <- furrr::future_map_chr( - seq_len(nrow(job_list)), - function(i) { - - .runHmmerJob( - JOB_NAME = job_list$JOB_NAME[i], - FASTA = job_list$FASTA[i], - DB = job_list$DB[i], - Total_proteins = Total_proteins, - output_path = output_path, - db_paths = db_paths, - docker_image = docker_image, - threads = threads, - n_workers = n_workers, - verbose = verbose - ) - } -) - -future::plan(future::sequential) - - parquet_tbl <- tibble::tibble( - parquet = parquet_files, - db = job_list$DB -) - - final_parquets <- list() - -for (database_name in databases) { - - if(verbose) message("Combining ", database_name) - - db_files <- parquet_tbl |> - dplyr::filter( - db == database_name - ) |> - dplyr::pull(parquet) - - combined_tbl <- purrr::map( - db_files, - arrow::read_parquet - ) |> - dplyr::bind_rows() |> -dplyr::left_join(.parse_hmmer_profiles(db_paths[[database_name]]) |> - dplyr::select(query_name = profile_name, query_accession = profile_accession, description = profile_description), -by = "query_name") - - final_parquet <- file.path( - output_path, - paste0( - "protein_", - database_name, - ".parquet" - ) - ) - - .write_compressed_parquet( - combined_tbl, - final_parquet - ) - - DBI::dbWriteTable( - con, - name = paste0( - "protein_", - database_name - ), - value = combined_tbl, - overwrite = TRUE - ) - - final_parquets[[database_name]] <- final_parquet - - message( - "Created ", - basename(final_parquet) - ) -} - -invisible(final_parquets) - - # purrr::map(parquet_files, arrow::read_parquet) |> - # dplyr::bind_rows() |> - # .write_compressed_parquet(final_parquet) - - # message("Combined parquet written.") - - # arrow::read_parquet(final_parquet) |> - # DBI::dbWriteTable(conn = con, name = tools::file_path_sans_ext(basename(final_parquet)), overwrite = TRUE) -} - -#' Parse HMMER tabular output into a tibble -#' -#' Reads a HMMER `--domtblout` file and returns a tidy tibble with one row per -#' target-query hit. Comment lines are stripped and the free-text description -#' field is reunited from the remaining whitespace-delimited columns. -#' -#' @param file Path to a HMMER `.tbl` output file produced with `--domtblout`. -#' -#' @return A tibble with 19 columns matching the HMMER per-sequence hit table -#' -#' @references Adapted from the rhmmer package -#' (). -#' -#' @examples -#' \dontrun{ -#' hits <- .parseHMMEROutput("results/Ecoli/protein_chunk_01_COG.tbl") -#' hits |> dplyr::filter(sequence_evalue < 1e-5) -#' } -#' -#' @keywords internal -.parseHMMEROutput <- function(file) { - - # target name accession tlen query name accession qlen E-value score bias # of c-Evalue i-Evalue score bias from to from to from to acc description of target - col_types <- readr::cols( - protein = readr::col_character(), # target name - protein_accession = readr::col_character(), - tlen = readr::col_integer(), - - query_name = readr::col_character(), # query name - query_accession = readr::col_character(), - qlen = readr::col_integer(), - - sequence_evalue = readr::col_double(), - sequence_score = readr::col_double(), - sequence_bias = readr::col_double(), - - domain_num = readr::col_integer(), - domain_of = readr::col_integer(), - - c_evalue = readr::col_double(), - i_evalue = readr::col_double(), - - domain_score = readr::col_double(), - domain_bias = readr::col_double(), - - hmm_from = readr::col_integer(), - hmm_to = readr::col_integer(), - - ali_from = readr::col_integer(), - ali_to = readr::col_integer(), - - env_from = readr::col_integer(), - env_to = readr::col_integer(), - - acc = readr::col_double(), - - target_description = readr::col_character() -) - # the line delimiter should always be just "\n", even on Windows - lines <- readr::read_lines(file, lazy = FALSE, progress = FALSE) - - # drop comment lines - data_lines <- lines[!grepl("^#", lines)] - - # split: whitespace-separated fields - split_fields <- strsplit(data_lines, "\\s+", perl = TRUE) - - # count space separated fields - N <- max(sapply(split_fields, length)) - - # Parsing differently to avoid fussy read_tsv() warnings - txt <- sub( - pattern = sprintf("(%s).*", paste0(rep("\\S+", N), collapse = " +")), - replacement = "\\1", - x = lines, - perl = TRUE - ) |> - gsub(pattern = " *", replacement = "\t") |> - paste0(collapse = "\n") - - table <- readr::read_tsv( - I(txt), - col_names = names(col_types$cols), - comment = "#", - na = "-", - col_types = col_types, - lazy = FALSE, - progress = FALSE - ) - - table -} - -#' Map HMMER protein annotations to genome-level count matrix and load into DuckDB -#' -#' Reads a Parquet file of HMMER hits (produced by [.runHMMER()]), joins the -#' annotations to the protein-cluster count matrix already in DuckDB, aggregates -#' 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()]). -#' -#' @return Invisibly returns the path to the written count Parquet file. -#' -#' @seealso [CDHIT2duckdb()], [runDataProcessing()] -#' -#' @keywords internal -.proteinAnnotations2Duckdb <- function( - duckdb_path, - databases = c("Pfam", "COG", "AMRFinder") -) { - - duckdb_path <- .docker_path(duckdb_path) - - con <- DBI::dbConnect( - duckdb::duckdb(), - duckdb_path - ) - - on.exit( - try( - DBI::dbDisconnect( - con, - shutdown = FALSE - ), - silent = TRUE - ), - add = TRUE - ) - - protein_long <- DBI::dbReadTable( - con, - "protein_count" - ) |> - tibble::as_tibble() |> - tidyr::pivot_longer( - cols = -genome_id, - names_to = "protein", - values_to = "count" - ) |> - dplyr::filter(count > 0) |> - dplyr::mutate( - protein = stringr::str_replace( - protein, - "^fig\\.", - "fig|" - ) - ) - - count_paths <- list() - - for (database in databases) { - - annotation_table <- paste0( - "protein_", - database - ) - - if (!DBI::dbExistsTable(con, annotation_table)) { - - warning( - annotation_table, - " not found in DuckDB. Skipping." - ) - - next - } - - message( - "Processing ", - annotation_table - ) - - annotation <- DBI::dbReadTable( - con, - annotation_table - ) |> - tibble::as_tibble() - - genome_annot_matrix <- protein_long |> - dplyr::inner_join( - annotation |> - dplyr::select( - protein, - query_name - ), - by = "protein" - ) |> - dplyr::group_by( - genome_id, - query_name - ) |> - dplyr::summarise( - count = sum(count), - .groups = "drop" - ) |> - tidyr::pivot_wider( - names_from = query_name, - values_from = count, - values_fill = 0 - ) - - count_table <- paste0( - annotation_table, - "_count" - ) - - count_path <- file.path( - dirname(duckdb_path), - paste0( - count_table, - ".parquet" - ) - ) - - arrow::write_parquet( - genome_annot_matrix, - count_path - ) - - DBI::dbWriteTable( - con, - count_table, - genome_annot_matrix, - overwrite = TRUE - ) - - count_paths[[database]] <- count_path - - message( - "Created ", - count_table - ) - } - - invisible(count_paths) -} - -#' Annotate proteins using DefenseFinder + CasFinder HMMs -#' Will add to the duckdb + create the parquet file. -#' -#' @param defense_db_dir Directory used to store downloaded HMMs -#' @param docker_image Docker image containing HMMER -#' @param duckdb_path DuckDB database path -#' @param output_path Output directory -#' @param threads Number of HMMER threads -#' -#' @returns Path to annotation parquet -#' @keywords internal -.defenseHMMER <- function( - defense_db_dir, - docker_image = "staphb/hmmer", - duckdb_path = "inst/extdata/Sfl.duckdb", - output_path = NULL, - threads = 8L, - verbose = TRUE -) { - - if (!nzchar(Sys.which("docker"))) { - stop("Docker is required.") - } - - defense_db_dir <- normalizePath( - defense_db_dir, - mustWork = FALSE - ) - - if (is.null(output_path)) { - output_path <- dirname( - normalizePath( - duckdb_path, - mustWork = FALSE - ) - ) - } - - dir.create( - defense_db_dir, - recursive = TRUE, - showWarnings = FALSE - ) - - dir.create( - output_path, - recursive = TRUE, - showWarnings = FALSE - ) - - #################################################################### - # download repositories - #################################################################### - - defense_dir <- file.path( - defense_db_dir, - "DefenseFinder" - ) - - cas_dir <- file.path( - defense_db_dir, - "CasFinder" - ) - - if (!dir.exists(defense_dir)) { - - if(verbose) message( - "Downloading DefenseFinder models" - ) - - tmp <- tempfile(fileext = ".zip") - - utils::download.file( - "https://github.com/mdmparis/defense-finder-models/archive/refs/heads/master.zip", - tmp, - mode = "wb", - method = "libcurl" - ) - - utils::unzip( - tmp, - exdir = defense_dir - ) - - unlink(tmp) - } - - if (!dir.exists(cas_dir)) { - - if(verbose) message( - "Downloading CasFinder models" - ) - - tmp <- tempfile(fileext = ".zip") - - utils::download.file( - "https://github.com/macsy-models/CasFinder/archive/refs/heads/main.zip", - tmp, - mode = "wb", - method = "libcurl" - ) - - utils::unzip( - tmp, - exdir = cas_dir - ) - - unlink(tmp) - } - - #################################################################### - # helper - #################################################################### - - build_database <- function( - repo_dir, - db_name - ) { - - profile_dirs <- list.dirs( - repo_dir, - recursive = TRUE, - full.names = TRUE - ) - - profile_dirs <- profile_dirs[ - basename(profile_dirs) == "profiles" - ] - - # moving to purrr implementation - hmm_files <- profile_dirs |> - purrr::map(\(x) list.files(x, - pattern = "\\.hmm$", - recursive = TRUE, - full.names = TRUE, - ignore.case = TRUE)) |> - purrr::flatten_chr() |> - unique() - - if (length(hmm_files) == 0) { - - stop( - "No HMM files found for ", - db_name - ) - } - - valid_hmms <- purrr::map_lgl(hmm_files, .isValidHmmFile) - - if (any(!valid_hmms)) { - bad_files <- basename(hmm_files[!valid_hmms]) - - if (isTRUE(verbose)) { - warning( - "Ignoring ", - length(bad_files), - " invalid HMM file(s):\n", - paste(bad_files, collapse = "\n"), - call. = FALSE - ) - } - - hmm_files <- hmm_files[valid_hmms] - } - -if (length(hmm_files) == 0) { - stop( - "No valid HMM files found for ", - db_name - ) -} - - combined_hmm <- file.path( - repo_dir, - paste0( - db_name, - ".hmm" - ) - ) - - if (file.exists(combined_hmm)) { - unlink(combined_hmm) - } - - file.create(combined_hmm) - - for (f in sort(hmm_files)) { - - file.append( - combined_hmm, - f - ) - } - - pressed_files <- paste0( - combined_hmm, - c( - ".h3m", - ".h3i", - ".h3f", - ".h3p" - ) - ) - - if (!all(file.exists(pressed_files))) { - - if(verbose) message( - "Running hmmpress for ", - db_name - ) - - output <- system2( - "docker", - args = c( - "run", - "--rm", - "-v", - paste0( - dirname(combined_hmm), - ":/db" - ), - docker_image, - "hmmpress", - file.path( - "/db", - basename(combined_hmm) - ) - ), - stdout = TRUE, - stderr = TRUE - ) - - if (!all(file.exists(pressed_files))) { - - stop( - "hmmpress failed for ", - db_name, - "\n", - paste(output, - collapse = "\n") - ) - } - } - - combined_hmm - } - - #################################################################### - # build separate databases - #################################################################### - - defense_hmm <- build_database( - defense_dir, - "DefenseFinder" - ) - - cas_hmm <- build_database( - cas_dir, - "CasFinder" - ) - - #################################################################### - # load proteins - #################################################################### - - con <- DBI::dbConnect( - duckdb::duckdb(), - duckdb_path - ) - - on.exit( - try( - DBI::dbDisconnect( - con, - shutdown = FALSE - ), - silent = TRUE - ), - add = TRUE - ) - - prot_seqs <- DBI::dbReadTable( - con, - "protein_cluster_seq" - ) |> - tibble::as_tibble() - - fasta_file <- file.path( - output_path, - "protein_DefenseCas.faa" - ) - - # required to define the database size for hmmsearch --Z and --domZ parameters - Total_proteins <- nrow(prot_seqs) - - readr::write_lines( - paste0( - ">", - prot_seqs$name, - "\n", - prot_seqs$sequence - ), - fasta_file - ) - - #################################################################### - # run hmmsearch separately - #################################################################### - - databases <- list( - DefenseFinder = defense_hmm, - CasFinder = cas_hmm - ) - - all_hits <- list() - - for (db_name in names(databases)) { - - if(verbose) message( - "Running ", - db_name - ) - - hmm_file <- databases[[db_name]] - - tbl_file <- file.path( - output_path, - paste0( - "protein_", - db_name, - ".tbl" - ) - ) - - output <- system2( - "docker", - args = c( - "run", - "--rm", - "-v", - paste0(output_path, ":/work"), - "-v", - paste0(dirname(hmm_file), ":/db"), - docker_image, - "hmmsearch", - "--notextw", - "--cpu", - as.character(threads), - "-Z", Total_proteins, - "--domZ", Total_proteins, - "--domtblout", - file.path( - "/work", - basename(tbl_file) - ), - file.path( - "/db", - basename(hmm_file) - ), - "/work/protein_DefenseCas.faa" - ), - stdout = TRUE, - stderr = TRUE - ) - - if (!file.exists(tbl_file)) { - - warning( - "hmmsearch failed for ", - db_name, - "\n", - paste(output, collapse = "\n") - ) - - next -} - - hits <- .parseHMMEROutput( - tbl_file - ) |> - dplyr::select( - protein, - query_name - ) |> - dplyr::mutate( - database = db_name - )|> -dplyr::left_join(.parse_hmmer_profiles(hmm_file) |> - dplyr::select(query_name = profile_name, query_accession = profile_accession, description = profile_description), -by = "query_name") - - all_hits[[db_name]] <- hits - } - - #################################################################### - # merge at parquet stage - #################################################################### - - combined_tbl <- dplyr::bind_rows( - all_hits - ) - - parquet_file <- file.path( - output_path, - "protein_DefenseCas.parquet" - ) - - .write_compressed_parquet( - combined_tbl, - parquet_file - ) - - DBI::dbWriteTable( - con, - "protein_DefenseCas", - combined_tbl, - overwrite = TRUE - ) - - message( - "Created protein_DefenseCas" - ) - - invisible(parquet_file) -}