diff --git a/DESCRIPTION b/DESCRIPTION index ee2f1c5..e7e92a2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -56,4 +56,4 @@ Suggests: rmarkdown, writexl, testthat (>= 3.0.0) -Config/roxygen2/version: 8.0.0 +Config/roxygen2/version: 8.1.0 diff --git a/NAMESPACE b/NAMESPACE index cd61148..aedaaa0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ export(.extractAMRtable) export(.updateBVBRCdata) export(CDHIT2duckdb) export(buildClusterFeatureMap) +export(buildDyadFeatureMap) export(cleanData) export(cleanMetaData) export(exportProcessedData) diff --git a/R/data_processing.R b/R/data_processing.R index 12df346..45875c2 100644 --- a/R/data_processing.R +++ b/R/data_processing.R @@ -970,7 +970,6 @@ CDHIT2duckdb <- function(duckdb_path, #' @returns A list of paths to the database hmm files. #' #' @keywords internal -#' @examples .prepareHmmerDatabases <- function( hmmer_db_dir, databases = c("Pfam", "COG", "AMRFinder"), @@ -1281,17 +1280,17 @@ CDHIT2duckdb <- function(duckdb_path, #' 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 +#' @param JOB_NAME protein_chunk id +#' @param FASTA fasta sequences +#' @param DB HMM database +#' @param Total_proteins protein sequence count +#' @param output_path path for saving hmmer outputs +#' @param db_paths path to HMM database +#' @param docker_image hmmer docker image (ideally from dockerhub) +#' @param threads number of threads +#' @param n_workers number of parallel workers #' -#' @returns +#' @returns the filename of the parquet file with hmmer output post parsing #' #' @keywords internal .runHmmerJob <- function(JOB_NAME, FASTA, DB, Total_proteins, @@ -1375,19 +1374,16 @@ CDHIT2duckdb <- function(duckdb_path, #' 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 +#' @param duckdb_path path to the duckdb with protein sequences and list +#' @param output_path path where HMMER output will be saved +#' @param threads number of threads +#' @param hmmer_db_dir path to the directory where HMM databases are/will be downloaded +#' @param databases list of HMM databases +#' @param docker_image the docker image of HMMER +#' @param num_of_splits The number of splits of the protein sequence file for parallel processing +#' @param n_workers The number of parallel runs #' #' @keywords internal -#' @examples .runHMMER <- function(duckdb_path, output_path, threads = 8L, @@ -2203,16 +2199,18 @@ CDHIT2duckdb <- function(duckdb_path, } -# Clean BV-BRC metadata, then save as Parquet files -#' -#' @param duckdb_path -#' @param path -#' @param ref_file_path +#' Clean BV-BRC metadata, then save as Parquet files #' -#' @returns +#' @param duckdb_path Path to the **per-selection DuckDB** produced by +#' [prepareGenomes()] (e.g., `"data//.duckdb"`). This DB must +#' already contain the tables written by [prepareGenomes()] and the upstream +#' genome-processing steps. +#' @param path the path to working directory +#' @param ref_file_path Directory containing reference TSVs used by +#' [cleanMetaData()] and [cleanData()] for metadata harmonization. +#' Default: `"data_raw/"`. +#' #' @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 @@ -2358,15 +2356,15 @@ cleanMetaData <- function(duckdb_path, path, ref_file_path = "data_raw/") { invisible(TRUE) } -# Clean feature matrices, then save as Parquet files +#' Clean feature matrices, then save as Parquet files #' -#' @param duckdb_path -#' @param path +#' @param duckdb_path Path to the **per-selection DuckDB** produced by +#' [prepareGenomes()] (e.g., `"data//.duckdb"`). This DB must +#' already contain the tables written by [prepareGenomes()] and the upstream +#' genome-processing steps. +#' @param path the path to working directory #' -#' @returns #' @export -#' -#' @examples cleanData <- function(duckdb_path, path) { duckdb_path <- normalizePath(duckdb_path) # If no explicit path is provided (or a generic one), choose results// when @@ -2513,6 +2511,7 @@ cleanData <- function(duckdb_path, path) { names_to = "annotation", values_to = "value" ) |> + dplyr::rename(!!database := annotation) |> dplyr::filter(!is.na(value) & value != "") |> dplyr::mutate(value = as.integer(value)) |> writeCompressedParquet(count_parquet) diff --git a/R/feature_to_head.R b/R/feature_to_head.R new file mode 100644 index 0000000..711d796 --- /dev/null +++ b/R/feature_to_head.R @@ -0,0 +1,612 @@ +#' Build a protein-gene dyad feature network using DuckDB +#' +#' Constructs a two-mode (bipartite) network in which every node is either a +#' protein-gene dyad or a biological feature, and every edge links a dyad to one +#' of its features. Features include the dyad's own protein and gene, structural +#' (pangenome graph) gene arrangements, and HMMER annotations: protein domains +#' (Pfam), COGs, phage-defense systems (including Cas), and antimicrobial +#' resistance genes. +#' +#' @param duckdb_path Character. Path to the source per-selection DuckDB. The +#' associated Parquet files and provenance manifest are expected to live +#' alongside it. This is the same `duckdb_path` produced by +#' \code{\link{runDataProcessing}}. +#' @param additional_feature_scales Character vector of optional feature types to +#' include. If `NULL`, all HMMER databases recorded in the manifest are used, +#' plus `struct` to include pangenome graph triplets. If a character vector is +#' supplied, individual feature scales can be selected. +#' @param output_path Character or NULL. Directory where the output Parquet file +#' will be written. If NULL, the output is written alongside the source DuckDB. +#' +#' @details +#' The function performs the following steps: +#' \enumerate{ +#' \item Creates a protein-gene mapping from +#' \code{genome_gene_protein.parquet}. +#' \item Constructs a unique protein-gene dyad identifier of the form +#' \code{"protein|gene"}. +#' \item Optionally loads structural and HMMER annotations. +#' \item Generates a network edge list linking each dyad to its associated +#' features. +#' \item Writes the resulting edge list to a compressed Parquet file. +#' } +#' +#' Every edge runs from a `protein|gene` dyad to a type-prefixed feature node: +#' +#' \preformatted{ +#' protein|gene --> protein:PROTEIN_ID +#' protein|gene --> gene:GENE_ID +#' protein|gene --> pfam:PFXXXXX +#' protein|gene --> cog:COGXXXX +#' protein|gene --> amr:GENE_NAME +#' protein|gene --> defensecas:DEFENSE_SYSTEM +#' protein|gene --> struct:STRUCTURE +#' } +#' +#' The type prefix on each target keeps feature namespaces separate, so a +#' generic feature name cannot collide across scales downstream. +#' +#' The output edge list contains two columns: +#' \describe{ +#' \item{source}{Protein-gene dyad identifier.} +#' \item{target}{Feature node identifier prefixed by feature type.} +#' } +#' +#' @return Invisibly returns the path to the generated +#' \code{dyad_feature.parquet} file. +#' +#' @examples +#' \dontrun{ +#' buildDyadFeatureMap( +#' duckdb_path = "data/Staphylococcus_argenteus/Sar.duckdb" +#' ) +#' } +#' +#' @import DBI duckdb +#' @export +buildDyadFeatureMap <- function( + duckdb_path, + additional_feature_scales = NULL, + output_path = NULL +) { + + duckdb_path <- normalizePath( + duckdb_path, + mustWork = TRUE + ) + + parquet_dir <- dirname(duckdb_path) + + manifest_path <- .manifest_find_latest( + duckdb_path + ) + + if (is.null(manifest_path)) { + stop( + "No provenance manifest found for: ", + duckdb_path + ) + } + + # Read in the `runDataProcessing()` output manifest + manifest <- jsonlite::read_json( + manifest_path, + simplifyVector = FALSE + ) + + hmmer_stage <- NULL + + # Find the most recent run whose HMMER stage completed successfully; its + # recorded databases determine which feature scales are available. + 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 + } + } + + # Without a successful HMMER stage the annotation Parquets we join against + # are not guaranteed to exist, so there is nothing to build from. + if (is.null(hmmer_stage)) { + stop( + "No successful HMMER stage found in manifest: ", + manifest_path + ) + } + + hmmer_databases <- unlist( + hmmer_stage$parameters$databases + ) + + # A successful stage with no databases recorded should not happen; fail loudly. + if (!length(hmmer_databases)) { + stop( + "HMMER stage in manifest does not contain any databases." + ) + } + + # Feature scales the caller is allowed to request: the structural view plus + # every HMMER database from the manifest (so new databases are picked up + # automatically). + allowed_features <- c( + "struct", + hmmer_databases + ) + + # Default to every allowed scale; otherwise reject anything unrecognised. + if (is.null(additional_feature_scales)) { + additional_feature_scales <- allowed_features + } else { + unknown_features <- setdiff( + additional_feature_scales, + allowed_features + ) + + if (length(unknown_features)) { + stop( + "Unsupported feature scale(s): ", + paste(unknown_features, collapse = ", "), + ". Available features: ", + paste(allowed_features, collapse = ", ") + ) + } + } + + out_dir <- if (is.null(output_path)) { + parquet_dir + } else { + normalizePath( + output_path, + winslash = "/", + mustWork = FALSE + ) + } + + dir.create( + out_dir, + recursive = TRUE, + showWarnings = FALSE + ) + + parquet_path <- file.path( + out_dir, + "dyad_feature.parquet" + ) + + con <- DBI::dbConnect( + duckdb::duckdb(), + dbdir = ":memory:" + ) + + on.exit( + DBI::dbDisconnect( + con, + shutdown = TRUE + ), + add = TRUE + ) + + # Local helpers ------------------------------------------------------------- + + # Escape single quotes so a path can be embedded in a SQL string literal. + .sql_escape <- function(x) { + gsub( + "'", + "''", + x, + fixed = TRUE + ) + } + + # Resolve "/.parquet" to an escaped absolute path, + # erroring if the expected Parquet file is missing. + .parquet_dataset_sql <- function( + parquet_dir, + dataset_name + ) { + + path <- file.path( + parquet_dir, + paste0( + dataset_name, + ".parquet" + ) + ) + + if (!file.exists(path)) { + stop( + "Required Parquet file not found: ", + path + ) + } + + .sql_escape( + normalizePath( + path, + winslash = "/", + mustWork = TRUE + ) + ) + } + + # ========================= + # Gene -> protein mapping + # ========================= + + # Base protein/gene pairs. Empty strings are excluded alongside NULLs to match + # the `value != ""` convention in data_processing.R and to avoid emitting + # bogus "protein|" or "|gene" dyads. + sql_path <- .parquet_dataset_sql( + parquet_dir, + "genome_gene_protein" + ) + + DBI::dbExecute( + con, + sprintf( + " + CREATE OR REPLACE VIEW protein_gene AS + SELECT DISTINCT + protein_ids AS protein, + REPLACE(Gene, '~', '.') AS gene + FROM read_parquet('%s') + WHERE protein_ids IS NOT NULL + AND protein_ids <> '' + AND Gene IS NOT NULL + AND Gene <> '' + ", + sql_path + ) + ) + + # ========================= + # Protein-gene -> dyad + # ========================= + + # Collapse each protein/gene pair into a single "protein|gene" dyad id; this + # is the source node for every edge in the output network. + DBI::dbExecute( + con, + " + CREATE OR REPLACE VIEW protein_gene_dyad AS + SELECT DISTINCT + protein, + gene, + CONCAT(protein, '|', gene) AS dyad + FROM protein_gene + " + ) + + # Track which optional feature views were successfully created, so the edge + # queries below only join against views that exist. + has_struct <- FALSE + has_pfam <- FALSE + has_cog <- FALSE + has_amr <- FALSE + has_defensecas <- FALSE + + # ========================= + # Structural gene features + # ========================= + + if ("struct" %in% additional_feature_scales) { + + struct_path <- file.path( + parquet_dir, + "struct.parquet" + ) + + if (!file.exists(struct_path)) { + + # Same skip-and-continue behaviour as .create_feature_view() uses for the + # HMMER datasets, kept inline here because the struct view is built with a + # different (UNNEST) query. + message( + "Skipping struct: no parquet found. Generate struct parquet first." + ) + + } else { + + sql_path <- .parquet_dataset_sql( + parquet_dir, + "struct" + ) + + DBI::dbExecute( + con, + sprintf( + " + CREATE OR REPLACE VIEW v_struct_genes AS + SELECT DISTINCT + struct, + gene + FROM read_parquet('%s') s + CROSS JOIN UNNEST( + string_split(s.struct, '.') + ) AS t(gene) + WHERE s.value = 1 + ", + sql_path + ) + ) + + has_struct <- TRUE + } + } + + # ========================== + # Generic HMMER feature view + # ========================== + + # Local helper: build a `protein -> feature` view from one HMMER annotation + # Parquet. `feature_expr` is the SQL expression that derives the feature label + # from `query_name` (identity for most databases, a REPLACE() for AMRFinder). + # Returns TRUE if the view was created, FALSE if the Parquet was missing. + .create_feature_view <- function( + con, + view_name, + parquet_dir, + dataset_name, + feature_expr + ) { + + sql_path <- tryCatch( + .parquet_dataset_sql( + parquet_dir, + dataset_name + ), + error = function(e) { + message( + "Skipping ", + dataset_name, + ": ", + e$message + ) + NULL + } + ) + + if (is.null(sql_path)) { + return(FALSE) + } + + DBI::dbExecute( + con, + sprintf( + " + CREATE OR REPLACE VIEW %s AS + SELECT DISTINCT + protein, + %s AS feature + FROM read_parquet('%s') + WHERE protein IS NOT NULL + AND query_name IS NOT NULL + ", + view_name, + feature_expr, + sql_path + ) + ) + + TRUE + } + + # ========================= + # HMMER annotation features + # ========================= + + if ("Pfam" %in% additional_feature_scales) { + + has_pfam <- .create_feature_view( + con = con, + view_name = "v_pfam", + parquet_dir = parquet_dir, + dataset_name = "protein_Pfam", + feature_expr = "query_name" + ) + } + + if ("COG" %in% additional_feature_scales) { + + has_cog <- .create_feature_view( + con = con, + view_name = "v_cog", + parquet_dir = parquet_dir, + dataset_name = "protein_COG", + feature_expr = "query_name" + ) + } + + if ("AMRFinder" %in% additional_feature_scales) { + + has_amr <- .create_feature_view( + con = con, + view_name = "v_amrfinder", + parquet_dir = parquet_dir, + dataset_name = "protein_AMRFinder", + feature_expr = "REPLACE(REPLACE(query_name, '-NCBIFAM', ''), '-', '.')" + ) + } + + if ("DefenseCas" %in% additional_feature_scales) { + + has_defensecas <- .create_feature_view( + con = con, + view_name = "v_defensecas", + parquet_dir = parquet_dir, + dataset_name = "protein_DefenseCas", + feature_expr = "query_name" + ) + } + + # ========================= + # Build network edge queries + # ========================= + + # Each entry is a SELECT returning (source, target) rows that are UNIONed into + # the final edge list. In the joined queries below `pgd` aliases the + # `protein_gene_dyad` view, so `pgd.dyad` is the "protein|gene" dyad id. + edge_queries <- c( + " + SELECT DISTINCT + dyad AS source, + CONCAT('protein:', protein) AS target + FROM protein_gene_dyad + ", + " + SELECT DISTINCT + dyad AS source, + CONCAT('gene:', gene) AS target + FROM protein_gene_dyad + " + ) + + # ========================= + # Structure edges + # ========================= + + if (has_struct) { + + edge_queries <- c( + edge_queries, + " + SELECT DISTINCT + pgd.dyad AS source, + CONCAT('struct:', sg.struct) AS target + FROM protein_gene_dyad pgd + JOIN v_struct_genes sg + ON pgd.gene = sg.gene + " + ) + } + + # ========================= + # Pfam edges + # ========================= + + if (has_pfam) { + + edge_queries <- c( + edge_queries, + " + SELECT DISTINCT + pgd.dyad AS source, + CONCAT('pfam:', pf.feature) AS target + FROM protein_gene_dyad pgd + JOIN v_pfam pf + ON pgd.protein = pf.protein + " + ) + } + + # ========================= + # COG edges + # ========================= + + if (has_cog) { + + edge_queries <- c( + edge_queries, + " + SELECT DISTINCT + pgd.dyad AS source, + CONCAT('cog:', cf.feature) AS target + FROM protein_gene_dyad pgd + JOIN v_cog cf + ON pgd.protein = cf.protein + " + ) + } + + # ========================= + # AMRFinder edges + # ========================= + + if (has_amr) { + + edge_queries <- c( + edge_queries, + " + SELECT DISTINCT + pgd.dyad AS source, + CONCAT('amr:', af.feature) AS target + FROM protein_gene_dyad pgd + JOIN v_amrfinder af + ON pgd.protein = af.protein + " + ) + } + + # ========================= + # DefenseCas edges + # ========================= + + if (has_defensecas) { + + edge_queries <- c( + edge_queries, + " + SELECT DISTINCT + pgd.dyad AS source, + CONCAT('defensecas:', df.feature) AS target + FROM protein_gene_dyad pgd + JOIN v_defensecas df + ON pgd.protein = df.protein + " + ) + } + + # ========================= + # Final edge list + # ========================= + + DBI::dbExecute( + con, + paste0( + " + CREATE OR REPLACE VIEW network_edges AS + ", + paste( + edge_queries, + collapse = "\nUNION\n" + ) + ) + ) + + # ========================= + # Export Parquet + # ========================= + + parquet_sql <- DBI::dbQuoteString( + con, + normalizePath( + parquet_path, + winslash = "/", + mustWork = FALSE + ) + ) + + DBI::dbExecute( + con, + paste0( + " + COPY network_edges + TO ", + parquet_sql, + " + (FORMAT PARQUET, COMPRESSION ZSTD) + " + ) + ) + + invisible(parquet_path) +} diff --git a/data/Staphylococcus_epidermidis/dyad_feature.parquet b/data/Staphylococcus_epidermidis/dyad_feature.parquet new file mode 100644 index 0000000..f7951e1 Binary files /dev/null and b/data/Staphylococcus_epidermidis/dyad_feature.parquet differ diff --git a/data/Staphylococcus_epidermidis/genome_gene_protein.parquet b/data/Staphylococcus_epidermidis/genome_gene_protein.parquet new file mode 100644 index 0000000..0873eed Binary files /dev/null and b/data/Staphylococcus_epidermidis/genome_gene_protein.parquet differ diff --git a/data/Staphylococcus_epidermidis/protein_AMRFinder.parquet b/data/Staphylococcus_epidermidis/protein_AMRFinder.parquet new file mode 100644 index 0000000..d7bcd87 Binary files /dev/null and b/data/Staphylococcus_epidermidis/protein_AMRFinder.parquet differ diff --git a/data/Staphylococcus_epidermidis/protein_COG.parquet b/data/Staphylococcus_epidermidis/protein_COG.parquet new file mode 100644 index 0000000..2350f5d Binary files /dev/null and b/data/Staphylococcus_epidermidis/protein_COG.parquet differ diff --git a/data/Staphylococcus_epidermidis/struct.parquet b/data/Staphylococcus_epidermidis/struct.parquet new file mode 100644 index 0000000..129fa72 Binary files /dev/null and b/data/Staphylococcus_epidermidis/struct.parquet differ diff --git a/man/buildDyadFeatureMap.Rd b/man/buildDyadFeatureMap.Rd new file mode 100644 index 0000000..43fcd11 --- /dev/null +++ b/man/buildDyadFeatureMap.Rd @@ -0,0 +1,80 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/feature_to_head.R +\name{buildDyadFeatureMap} +\alias{buildDyadFeatureMap} +\title{Build a protein-gene dyad feature network using DuckDB} +\usage{ +buildDyadFeatureMap( + duckdb_path, + additional_feature_scales = NULL, + output_path = NULL +) +} +\arguments{ +\item{duckdb_path}{Character. Path to the source per-selection DuckDB. The +associated Parquet files and provenance manifest are expected to live +alongside it. This is the same \code{duckdb_path} produced by +\code{\link{runDataProcessing}}.} + +\item{additional_feature_scales}{Character vector of optional feature types to +include. If \code{NULL}, all HMMER databases recorded in the manifest are used, +plus \code{struct} to include pangenome graph triplets. If a character vector is +supplied, individual feature scales can be selected.} + +\item{output_path}{Character or NULL. Directory where the output Parquet file +will be written. If NULL, the output is written alongside the source DuckDB.} +} +\value{ +Invisibly returns the path to the generated +\code{dyad_feature.parquet} file. +} +\description{ +Constructs a two-mode (bipartite) network in which every node is either a +protein-gene dyad or a biological feature, and every edge links a dyad to one +of its features. Features include the dyad's own protein and gene, structural +(pangenome graph) gene arrangements, and HMMER annotations: protein domains +(Pfam), COGs, phage-defense systems (including Cas), and antimicrobial +resistance genes. +} +\details{ +The function performs the following steps: +\enumerate{ +\item Creates a protein-gene mapping from +\code{genome_gene_protein.parquet}. +\item Constructs a unique protein-gene dyad identifier of the form +\code{"protein|gene"}. +\item Optionally loads structural and HMMER annotations. +\item Generates a network edge list linking each dyad to its associated +features. +\item Writes the resulting edge list to a compressed Parquet file. +} + +Every edge runs from a \code{protein|gene} dyad to a type-prefixed feature node: + +\preformatted{ +protein|gene --> protein:PROTEIN_ID +protein|gene --> gene:GENE_ID +protein|gene --> pfam:PFXXXXX +protein|gene --> cog:COGXXXX +protein|gene --> amr:GENE_NAME +protein|gene --> defensecas:DEFENSE_SYSTEM +protein|gene --> struct:STRUCTURE +} + +The type prefix on each target keeps feature namespaces separate, so a +generic feature name cannot collide across scales downstream. + +The output edge list contains two columns: +\describe{ +\item{source}{Protein-gene dyad identifier.} +\item{target}{Feature node identifier prefixed by feature type.} +} +} +\examples{ +\dontrun{ +buildDyadFeatureMap( + duckdb_path = "data/Staphylococcus_argenteus/Sar.duckdb" +) +} + +} diff --git a/man/cleanData.Rd b/man/cleanData.Rd new file mode 100644 index 0000000..ee2ba2b --- /dev/null +++ b/man/cleanData.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data_processing.R +\name{cleanData} +\alias{cleanData} +\title{Clean feature matrices, then save as Parquet files} +\usage{ +cleanData(duckdb_path, path) +} +\arguments{ +\item{duckdb_path}{Path to the \strong{per-selection DuckDB} produced by +\code{\link[=prepareGenomes]{prepareGenomes()}} (e.g., \code{"data//.duckdb"}). This DB must +already contain the tables written by \code{\link[=prepareGenomes]{prepareGenomes()}} and the upstream +genome-processing steps.} + +\item{path}{the path to working directory} +} +\description{ +Clean feature matrices, then save as Parquet files +} diff --git a/man/cleanMetaData.Rd b/man/cleanMetaData.Rd new file mode 100644 index 0000000..b2a0912 --- /dev/null +++ b/man/cleanMetaData.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data_processing.R +\name{cleanMetaData} +\alias{cleanMetaData} +\title{Clean BV-BRC metadata, then save as Parquet files} +\usage{ +cleanMetaData(duckdb_path, path, ref_file_path = "data_raw/") +} +\arguments{ +\item{duckdb_path}{Path to the \strong{per-selection DuckDB} produced by +\code{\link[=prepareGenomes]{prepareGenomes()}} (e.g., \code{"data//.duckdb"}). This DB must +already contain the tables written by \code{\link[=prepareGenomes]{prepareGenomes()}} and the upstream +genome-processing steps.} + +\item{path}{the path to working directory} + +\item{ref_file_path}{Directory containing reference TSVs used by +\code{\link[=cleanMetaData]{cleanMetaData()}} and \code{\link[=cleanData]{cleanData()}} for metadata harmonization. +Default: \code{"data_raw/"}.} +} +\description{ +Clean BV-BRC metadata, then save as Parquet files +} diff --git a/man/domainFromIPR.Rd b/man/domainFromIPR.Rd deleted file mode 100644 index 6b3229a..0000000 --- a/man/domainFromIPR.Rd +++ /dev/null @@ -1,23 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_processing.R -\name{domainFromIPR} -\alias{domainFromIPR} -\title{Derive protein domain presence/absence and counts via InterProScan and write to DuckDB} -\usage{ -domainFromIPR( - 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" -) -} -\description{ -Derive protein domain presence/absence and counts via InterProScan and write to DuckDB -} diff --git a/man/dot-checkInterProData.Rd b/man/dot-checkInterProData.Rd deleted file mode 100644 index 7f4c794..0000000 --- a/man/dot-checkInterProData.Rd +++ /dev/null @@ -1,36 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_processing.R -\name{.checkInterProData} -\alias{.checkInterProData} -\title{Check or install InterProScan data bundle} -\usage{ -.checkInterProData( - 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 -) -} -\arguments{ -\item{version}{InterProScan version string.} - -\item{dest_dir}{Directory where data should be installed.} - -\item{docker_image}{Docker image string for InterProScan.} - -\item{platform}{Character indicating Docker platform (e.g. \code{"linux/amd64"}).} - -\item{curl_bin}{Path to curl executable.} - -\item{verbose}{Logical; print status messages.} -} -\value{ -A list containing \code{data_dir} and \code{ready} status. -} -\description{ -Ensures that the InterProScan data directory exists locally, downloading -and verifying the appropriate tarball when necessary. -} -\keyword{internal} diff --git a/man/dot-defenseHMMER.Rd b/man/dot-defenseHMMER.Rd index 2a3a772..c2ab6d9 100644 --- a/man/dot-defenseHMMER.Rd +++ b/man/dot-defenseHMMER.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/runHMMER.R +% Please edit documentation in R/data_processing.R \name{.defenseHMMER} \alias{.defenseHMMER} \title{Annotate proteins using DefenseFinder + CasFinder HMMs diff --git a/man/dot-docker_path.Rd b/man/dot-docker_path.Rd index d0b6f70..c3385ca 100644 --- a/man/dot-docker_path.Rd +++ b/man/dot-docker_path.Rd @@ -1,21 +1,12 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_curation.R, R/data_processing.R +% Please edit documentation in R/helpers.R \name{.docker_path} \alias{.docker_path} \title{Helps normalize Docker paths} \usage{ .docker_path(p) - -.docker_path(p) -} -\arguments{ -\item{p}{Character scalar. A filesystem path on the host OS.} -} -\value{ -A normalized path string. } \description{ -Converts Windows and mixed-separator paths to forward slashes -and applies \code{normalizePath()} without requiring the path to exist. +Helps normalize Docker paths } \keyword{internal} diff --git a/man/dot-getDfIPRColNames.Rd b/man/dot-getDfIPRColNames.Rd deleted file mode 100644 index 7151d2d..0000000 --- a/man/dot-getDfIPRColNames.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_processing.R -\name{.getDfIPRColNames} -\alias{.getDfIPRColNames} -\title{Internal helpers for reading InterProScan TSV outputs} -\usage{ -.getDfIPRColNames() -} -\arguments{ -\item{filepath}{Path to a \code{.tsv} or \code{.tsv.gz} InterProScan result file.} -} -\value{ -A tibble of parsed InterProScan output. -} -\description{ -Provide standardized column names, types, and a reader wrapper for the -InterProScan tab-delimited output format. -} -\keyword{internal} diff --git a/man/dot-getDfIPRColTypes.Rd b/man/dot-getDfIPRColTypes.Rd deleted file mode 100644 index c827bc1..0000000 --- a/man/dot-getDfIPRColTypes.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_processing.R -\name{.getDfIPRColTypes} -\alias{.getDfIPRColTypes} -\title{Internal helpers for reading InterProScan TSV outputs} -\usage{ -.getDfIPRColTypes() -} -\arguments{ -\item{filepath}{Path to a \code{.tsv} or \code{.tsv.gz} InterProScan result file.} -} -\value{ -A tibble of parsed InterProScan output. -} -\description{ -Provide standardized column names, types, and a reader wrapper for the -InterProScan tab-delimited output format. -} -\keyword{internal} diff --git a/man/dot-id_checker.Rd b/man/dot-id_checker.Rd index 00fd8eb..aa4a11e 100644 --- a/man/dot-id_checker.Rd +++ b/man/dot-id_checker.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_curation.R +% Please edit documentation in R/helpers.R \name{.id_checker} \alias{.id_checker} \title{Helps ensure trailing 0s are retained in genome IDs for proper downloading} diff --git a/man/dot-isValidHmmFile.Rd b/man/dot-isValidHmmFile.Rd index 1268a9b..4ba8824 100644 --- a/man/dot-isValidHmmFile.Rd +++ b/man/dot-isValidHmmFile.Rd @@ -1,12 +1,19 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/runHMMER.R +% Please edit documentation in R/helpers.R \name{.isValidHmmFile} \alias{.isValidHmmFile} -\title{Validate if a HMM file has old HMMER3 format and remove them} +\title{Validate if a HMM file has old HMMER3 format} \usage{ .isValidHmmFile(hmm_file) } +\arguments{ +\item{hmm_file}{Path to a \code{.hmm} file.} +} +\value{ +\code{TRUE} if the file has valid HMMER3 formatting (starts with +\code{HMMER3/f} and ends with \verb{//}), \code{FALSE} otherwise. +} \description{ -Validate if a HMM file has old HMMER3 format and remove them +Validate if a HMM file has old HMMER3 format } \keyword{internal} diff --git a/man/dot-manifest_event.Rd b/man/dot-manifest_event.Rd new file mode 100644 index 0000000..fc2182d --- /dev/null +++ b/man/dot-manifest_event.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.manifest_event} +\alias{.manifest_event} +\title{Append a provenance event to the active manifest run} +\usage{ +.manifest_event(manifest_state, level = "info", message, details = list()) +} +\arguments{ +\item{manifest_state}{Manifest state returned by \code{\link[=.manifest_start]{.manifest_start()}}.} + +\item{level}{Character event level.} + +\item{message}{Character message.} + +\item{details}{Optional named list.} +} +\value{ +Updated manifest state. +} +\description{ +Append a provenance event to the active manifest run +} +\keyword{internal} diff --git a/man/dot-manifest_file_info.Rd b/man/dot-manifest_file_info.Rd new file mode 100644 index 0000000..69e9f77 --- /dev/null +++ b/man/dot-manifest_file_info.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.manifest_file_info} +\alias{.manifest_file_info} +\title{Returns the basics about a file for manifest logging} +\usage{ +.manifest_file_info(path, hash = FALSE) +} +\arguments{ +\item{path}{Character vector of file paths.} + +\item{hash}{Logical. If TRUE, calculate MD5 checksums.} +} +\value{ +A list of file records. +} +\description{ +Returns the basics about a file for manifest logging +} +\keyword{internal} diff --git a/man/dot-manifest_finish.Rd b/man/dot-manifest_finish.Rd new file mode 100644 index 0000000..25f6495 --- /dev/null +++ b/man/dot-manifest_finish.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.manifest_finish} +\alias{.manifest_finish} +\title{Finish an active provenance manifest run} +\usage{ +.manifest_finish(manifest_state, status = "success", error = NULL) +} +\arguments{ +\item{manifest_state}{Manifest state returned by \code{\link[=.manifest_start]{.manifest_start()}}.} + +\item{status}{Final run status.} + +\item{error}{Optional error message.} +} +\value{ +Invisibly returns the final manifest state. +} +\description{ +Finish an active provenance manifest run +} +\keyword{internal} diff --git a/man/dot-manifest_git_info.Rd b/man/dot-manifest_git_info.Rd new file mode 100644 index 0000000..937c627 --- /dev/null +++ b/man/dot-manifest_git_info.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.manifest_git_info} +\alias{.manifest_git_info} +\title{Capture basic GitHub repo state for manifest provenance} +\usage{ +.manifest_git_info(base_dir = ".") +} +\arguments{ +\item{base_dir}{Character. Project root.} +} +\value{ +A named list. +} +\description{ +Capture basic GitHub repo state for manifest provenance +} +\keyword{internal} diff --git a/man/dot-manifest_package_versions.Rd b/man/dot-manifest_package_versions.Rd new file mode 100644 index 0000000..4df2167 --- /dev/null +++ b/man/dot-manifest_package_versions.Rd @@ -0,0 +1,15 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.manifest_package_versions} +\alias{.manifest_package_versions} +\title{Capture package versions currently loaded in the R session} +\usage{ +.manifest_package_versions() +} +\value{ +Named character vector of package versions. +} +\description{ +Capture package versions currently loaded in the R session +} +\keyword{internal} diff --git a/man/dot-manifest_resume.Rd b/man/dot-manifest_resume.Rd new file mode 100644 index 0000000..46f28b5 --- /dev/null +++ b/man/dot-manifest_resume.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.manifest_resume} +\alias{.manifest_resume} +\title{Resume provenance logging in an existing manifest} +\usage{ +.manifest_resume(manifest_path, base_dir = ".", hash_files = FALSE) +} +\arguments{ +\item{manifest_path}{Character. Path to an existing JSON manifest.} + +\item{base_dir}{Character. Project root.} + +\item{hash_files}{Logical. Calculate SHA-256 checksums for manifest-recorded files.} +} +\value{ +A manifest object with \code{path} and \code{run_index}. +} +\description{ +Loads an existing manifest and appends a new run. +} +\keyword{internal} diff --git a/man/dot-manifest_run_id.Rd b/man/dot-manifest_run_id.Rd new file mode 100644 index 0000000..28d860c --- /dev/null +++ b/man/dot-manifest_run_id.Rd @@ -0,0 +1,15 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.manifest_run_id} +\alias{.manifest_run_id} +\title{Generate a unique manifest run identifier} +\usage{ +.manifest_run_id() +} +\value{ +Character scalar. +} +\description{ +Generate a unique manifest run identifier +} +\keyword{internal} diff --git a/man/dot-manifest_stage.Rd b/man/dot-manifest_stage.Rd new file mode 100644 index 0000000..803a3fd --- /dev/null +++ b/man/dot-manifest_stage.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.manifest_stage} +\alias{.manifest_stage} +\title{Update a manifest stage} +\usage{ +.manifest_stage( + manifest_state, + name, + status = "success", + parameters = list(), + inputs = character(), + outputs = character(), + tool = list(), + metrics = list(), + message = NULL +) +} +\arguments{ +\item{manifest_state}{Manifest state returned by \code{\link[=.manifest_start]{.manifest_start()}}.} + +\item{name}{Character stage name.} + +\item{status}{Character stage status.} + +\item{parameters}{Optional named list.} + +\item{inputs}{Optional character vector of input paths.} + +\item{outputs}{Optional character vector of output paths.} + +\item{tool}{Optional named list describing the tool.} + +\item{metrics}{Optional named list of metrics.} + +\item{message}{Optional log message.} +} +\value{ +Updated manifest state. +} +\description{ +Update a manifest stage +} +\keyword{internal} diff --git a/man/dot-manifest_start.Rd b/man/dot-manifest_start.Rd new file mode 100644 index 0000000..078352e --- /dev/null +++ b/man/dot-manifest_start.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.manifest_start} +\alias{.manifest_start} +\title{Start or load a dataset provenance manifest} +\usage{ +.manifest_start( + manifest_path, + dataset_id, + duckdb_path, + base_dir = ".", + selection = list(), + hash_files = FALSE +) +} +\arguments{ +\item{manifest_path}{Character. Path to the JSON manifest.} + +\item{dataset_id}{Character scalar.} + +\item{duckdb_path}{Character scalar.} + +\item{base_dir}{Character scalar.} + +\item{selection}{Optional named list describing the dataset selection.} + +\item{hash_files}{Logical. Calculate SHA-256 for manifest-recorded files.} +} +\value{ +A manifest object with \code{path} and \code{run_index}. +} +\description{ +Start or load a dataset provenance manifest +} +\keyword{internal} diff --git a/man/dot-parseHMMEROutput.Rd b/man/dot-parseHMMEROutput.Rd index 5ace7e1..4b95caf 100644 --- a/man/dot-parseHMMEROutput.Rd +++ b/man/dot-parseHMMEROutput.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/runHMMER.R +% Please edit documentation in R/helpers.R \name{.parseHMMEROutput} \alias{.parseHMMEROutput} \title{Parse HMMER tabular output into a tibble} diff --git a/man/dot-parse_hmmer_profiles.Rd b/man/dot-parse_hmmer_profiles.Rd index 428796e..12f056e 100644 --- a/man/dot-parse_hmmer_profiles.Rd +++ b/man/dot-parse_hmmer_profiles.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/runHMMER.R +% Please edit documentation in R/helpers.R \name{.parse_hmmer_profiles} \alias{.parse_hmmer_profiles} \title{Parsing HMM database to extract profile names, accessions and descriptions} diff --git a/man/dot-pick_shell.Rd b/man/dot-pick_shell.Rd index fb4e8bf..694f49a 100644 --- a/man/dot-pick_shell.Rd +++ b/man/dot-pick_shell.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_curation.R +% Please edit documentation in R/helpers.R \name{.pick_shell} \alias{.pick_shell} \title{Helps run a shell inside a container, and prefers bash (don't we all?)} diff --git a/man/dot-prepareHmmerDatabases.Rd b/man/dot-prepareHmmerDatabases.Rd index f0879b4..9f8f0d1 100644 --- a/man/dot-prepareHmmerDatabases.Rd +++ b/man/dot-prepareHmmerDatabases.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/runHMMER.R +% Please edit documentation in R/data_processing.R \name{.prepareHmmerDatabases} \alias{.prepareHmmerDatabases} \title{Download and prepare HMMER databases for generating new file types.} diff --git a/man/dot-preserve_export_id_text.Rd b/man/dot-preserve_export_id_text.Rd new file mode 100644 index 0000000..fea92d8 --- /dev/null +++ b/man/dot-preserve_export_id_text.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{.preserve_export_id_text} +\alias{.preserve_export_id_text} +\title{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.} +\usage{ +.preserve_export_id_text(df) +} +\description{ +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. +} +\keyword{internal} diff --git a/man/dot-processPanaroo.Rd b/man/dot-processPanaroo.Rd index aa42f1f..f3aa5b7 100644 --- a/man/dot-processPanaroo.Rd +++ b/man/dot-processPanaroo.Rd @@ -12,7 +12,8 @@ cluster_threshold, family_seq_identity, panaroo_threads_per_job, - refind_mode = c("off", "default", "strict") + refind_mode = c("off", "default", "strict"), + verbose = TRUE ) } \arguments{ diff --git a/man/dot-process_chunk.Rd b/man/dot-process_chunk.Rd deleted file mode 100644 index b6c9c21..0000000 --- a/man/dot-process_chunk.Rd +++ /dev/null @@ -1,45 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_processing.R -\name{.process_chunk} -\alias{.process_chunk} -\title{Run InterProScan on a sequence chunk inside Docker} -\usage{ -.process_chunk( - 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") -) -} -\arguments{ -\item{chunk}{A tibble with columns \code{name} and \code{sequence}.} - -\item{path}{Working directory used for temporary files.} - -\item{ipr_data_path}{Path to InterProScan data directory.} - -\item{out_file_base}{Output prefix for chunk results.} - -\item{appl}{Character vector of InterProScan applications (e.g. \code{"Pfam"}).} - -\item{chunk_id}{Integer chunk index.} - -\item{threads}{Number of CPUs for InterProScan container.} - -\item{file_format}{Output format (\code{"TSV"}).} - -\item{docker_image}{InterProScan Docker image.} -} -\value{ -Path to a \code{.tsv} or \code{.tsv.gz} InterProScan output file. -} -\description{ -Executes InterProScan on a subset of protein sequences, writing temporary -FASTA and reading back \code{.tsv} or \code{.tsv.gz} results. -} -\keyword{internal} diff --git a/man/dot-proteinAnnotations2Duckdb.Rd b/man/dot-proteinAnnotations2Duckdb.Rd index e59930d..cfffa09 100644 --- a/man/dot-proteinAnnotations2Duckdb.Rd +++ b/man/dot-proteinAnnotations2Duckdb.Rd @@ -1,12 +1,13 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/runHMMER.R +% Please edit documentation in R/data_processing.R \name{.proteinAnnotations2Duckdb} \alias{.proteinAnnotations2Duckdb} \title{Map HMMER protein annotations to genome-level count matrix and load into DuckDB} \usage{ .proteinAnnotations2Duckdb( duckdb_path, - databases = c("Pfam", "COG", "AMRFinder") + databases, + output_path = dirname(duckdb_path) ) } \arguments{ diff --git a/man/dot-readIPRscanTsv.Rd b/man/dot-readIPRscanTsv.Rd deleted file mode 100644 index d4503e5..0000000 --- a/man/dot-readIPRscanTsv.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_processing.R -\name{.readIPRscanTsv} -\alias{.readIPRscanTsv} -\title{Internal helpers for reading InterProScan TSV outputs} -\usage{ -.readIPRscanTsv(filepath) -} -\arguments{ -\item{filepath}{Path to a \code{.tsv} or \code{.tsv.gz} InterProScan result file.} -} -\value{ -A tibble of parsed InterProScan output. -} -\description{ -Provide standardized column names, types, and a reader wrapper for the -InterProScan tab-delimited output format. -} -\keyword{internal} diff --git a/man/dot-runHMMER.Rd b/man/dot-runHMMER.Rd index 2ead2ac..72618a7 100644 --- a/man/dot-runHMMER.Rd +++ b/man/dot-runHMMER.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/runHMMER.R +% Please edit documentation in R/data_processing.R \name{.runHMMER} \alias{.runHMMER} \title{Wrapper for preparing HMM databases and running HMMER on protein sequences from duckdb and writing them.} @@ -16,6 +16,23 @@ verbose = TRUE ) } +\arguments{ +\item{duckdb_path}{path to the duckdb with protein sequences and list} + +\item{output_path}{path where HMMER output will be saved} + +\item{threads}{number of threads} + +\item{hmmer_db_dir}{path to the directory where HMM databases are/will be downloaded} + +\item{databases}{list of HMM databases} + +\item{docker_image}{the docker image of HMMER} + +\item{num_of_splits}{The number of splits of the protein sequence file for parallel processing} + +\item{n_workers}{The number of parallel runs} +} \description{ Wrapper for preparing HMM databases and running HMMER on protein sequences from duckdb and writing them. } diff --git a/man/dot-runHmmerJob.Rd b/man/dot-runHmmerJob.Rd index 5cc8340..121d513 100644 --- a/man/dot-runHmmerJob.Rd +++ b/man/dot-runHmmerJob.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/runHMMER.R +% Please edit documentation in R/data_processing.R \name{.runHmmerJob} \alias{.runHmmerJob} \title{The function to run HMMER with docker} @@ -17,6 +17,28 @@ verbose = TRUE ) } +\arguments{ +\item{JOB_NAME}{protein_chunk id} + +\item{FASTA}{fasta sequences} + +\item{DB}{HMM database} + +\item{Total_proteins}{protein sequence count} + +\item{output_path}{path for saving hmmer outputs} + +\item{db_paths}{path to HMM database} + +\item{docker_image}{hmmer docker image (ideally from dockerhub)} + +\item{threads}{number of threads} + +\item{n_workers}{number of parallel workers} +} +\value{ +the filename of the parquet file with hmmer output post parsing +} \description{ The function to run HMMER with docker } diff --git a/man/dot-runPanaroo.Rd b/man/dot-runPanaroo.Rd index ed083fa..c302921 100644 --- a/man/dot-runPanaroo.Rd +++ b/man/dot-runPanaroo.Rd @@ -16,7 +16,8 @@ refind_mode = c("off", "default", "strict"), strip_pseudogenes = FALSE, pseudogene_clean_dir = "gff_clean", - write_pseudogene_audit = TRUE + write_pseudogene_audit = TRUE, + verbose = TRUE ) } \arguments{ diff --git a/man/dot-stripPseudogeneGFFs.Rd b/man/dot-stripPseudogeneGFFs.Rd index f2e4471..3853720 100644 --- a/man/dot-stripPseudogeneGFFs.Rd +++ b/man/dot-stripPseudogeneGFFs.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_processing.R +% Please edit documentation in R/helpers.R \name{.stripPseudogeneGFFs} \alias{.stripPseudogeneGFFs} \title{Remove pseudogene annotations from Panaroo input GFF files} diff --git a/man/dot-to_container.Rd b/man/dot-to_container.Rd index ce8b31c..c12cc21 100644 --- a/man/dot-to_container.Rd +++ b/man/dot-to_container.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_processing.R +% Please edit documentation in R/helpers.R \name{.to_container} \alias{.to_container} \title{.to_container()} diff --git a/man/dot-write_compressed_parquet.Rd b/man/dot-write_compressed_parquet.Rd index b9e45de..586cf50 100644 --- a/man/dot-write_compressed_parquet.Rd +++ b/man/dot-write_compressed_parquet.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/runHMMER.R +% Please edit documentation in R/helpers.R \name{.write_compressed_parquet} \alias{.write_compressed_parquet} \title{Write a data frame to a compressed Parquet file} diff --git a/man/exportProcessedData.Rd b/man/exportProcessedData.Rd index c8f16d4..07aff2a 100644 --- a/man/exportProcessedData.Rd +++ b/man/exportProcessedData.Rd @@ -11,6 +11,7 @@ exportProcessedData( export_formats = c("csv"), export_sequences = FALSE, tables = NULL, + export_tables = TRUE, verbose = TRUE ) } diff --git a/man/exportTables.Rd b/man/exportTables.Rd index 4fffb78..98618a7 100644 --- a/man/exportTables.Rd +++ b/man/exportTables.Rd @@ -8,8 +8,9 @@ exportTables( duckdb_path, output_dir = NULL, tables = NULL, - skip_tables = c(NULL), + skip_tables = NULL, include_summary = TRUE, + export_formats = c("csv"), export_tables = TRUE, load_tables = FALSE, verbose = TRUE diff --git a/man/runDataProcessing.Rd b/man/runDataProcessing.Rd index fa2eda4..4a80dc7 100644 --- a/man/runDataProcessing.Rd +++ b/man/runDataProcessing.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/data_processing.R \name{runDataProcessing} \alias{runDataProcessing} -\title{Run the full amRdata processing pipeline (Panaroo → CD-HIT → InterProScan → Parquet)} +\title{Run the full amRdata processing pipeline (Panaroo -> CD-HIT -> HMMER -> Parquet)} \usage{ runDataProcessing( duckdb_path, @@ -22,12 +22,11 @@ runDataProcessing( cdhit_memory = 0, cdhit_extra_args = c("-g", "1"), cdhit_output_prefix = "cdhit_out", - 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, + hmmer_databases = c("Pfam", "COG", "AMRFinder", "DefenseCas"), + hmmer_db_dir = NULL, + hmmer_docker_image = "staphb/hmmer", + hmmer_num_splits = 8L, + hmmer_workers = 8L, ref_file_path = "data_raw/", verbose = TRUE ) @@ -35,17 +34,19 @@ runDataProcessing( \arguments{ \item{duckdb_path}{Character. Path to the \strong{per-selection DuckDB} produced by \code{\link[=prepareGenomes]{prepareGenomes()}} (e.g., \code{"data//.duckdb"}). This DB must -already contain at least the tables written by \code{prepareGenomes()} and subsequent -download steps (e.g., \code{files}, \code{filtered}, and metadata tables).} +already contain the tables written by \code{\link[=prepareGenomes]{prepareGenomes()}} and the upstream +genome-processing steps.} -\item{output_path}{Character or \code{NULL}. Base directory for writing Panaroo/CD-HIT/InterProScan -outputs and final Parquet files. If \code{NULL}, defaults to \code{dirname(duckdb_path)}.} +\item{output_path}{Character or \code{NULL}. Base directory for writing Panaroo, +CD-HIT, HMMER, and final Parquet outputs. If \code{NULL}, defaults to +\code{dirname(duckdb_path)}.} -\item{threads}{Integer. Shared concurrency budget used across tools (Panaroo, CD-HIT, -InterProScan). Passed through to each stage as appropriate. Defaults to \code{8}.} +\item{threads}{Integer. Shared concurrency budget used across Panaroo, CD-HIT, +and HMMER. Defaults to \code{8}.} -\item{panaroo_split_jobs}{Logical. If \code{TRUE}, Panaroo runs in multiple batches that can be -merged by \code{\link[=.mergePanaroo]{.mergePanaroo()}}. If \code{FALSE}, Panaroo runs once on all isolates. Default: \code{FALSE}.} +\item{panaroo_split_jobs}{Logical. If \code{TRUE}, Panaroo runs in multiple batches +that can be merged by \code{\link[=.mergePanaroo]{.mergePanaroo()}}. If \code{FALSE}, Panaroo runs once on all +isolates. Default: \code{FALSE}.} \item{panaroo_core_threshold}{Numeric. Panaroo \code{--core_threshold}. Default: \code{0.90}.} @@ -53,139 +54,145 @@ merged by \code{\link[=.mergePanaroo]{.mergePanaroo()}}. If \code{FALSE}, Panaro \item{panaroo_cluster_threshold}{Numeric. Panaroo \code{--threshold}. Default: \code{0.95}.} -\item{panaroo_family_seq_identity}{Numeric. Panaroo \code{-f} (gene family identity). Default: \code{0.5}.} +\item{panaroo_family_seq_identity}{Numeric. Panaroo \code{-f} gene family identity. +Default: \code{0.5}.} -\item{panaroo_refind_mode}{Character. Panaroo's \code{--refind-mode} (\code{"off"}, \code{"default"}, -or \code{"strict"}). See \code{\link[=.processPanaroo]{.processPanaroo()}} for what refinding does and the runtime -caveat behind the current default. Default \code{"off"}.} +\item{panaroo_refind_mode}{Character. Panaroo's \code{--refind-mode} (\code{"off"}, +\code{"default"}, or \code{"strict"}). See \code{\link[=.processPanaroo]{.processPanaroo()}} for the runtime caveat +behind refinding. Default: \code{"off"}.} + +\item{panaroo_strip_pseudogenes}{Logical. If \code{TRUE}, remove pseudogene feature +records from Panaroo input GFF files before running Panaroo. Default: \code{FALSE}.} + +\item{panaroo_pseudogene_clean_dir}{Character. Directory name for cleaned GFF +files. Default: \code{"gff_clean"}.} + +\item{panaroo_write_pseudogene_audit}{Logical. If \code{TRUE}, write a pseudogene +cleaning audit file. Default: \code{TRUE}.} \item{cdhit_identity}{Numeric. CD-HIT \code{-c} identity threshold. Default: \code{0.9}.} \item{cdhit_word_length}{Integer. CD-HIT \code{-n} word length. Default: \code{5}.} -\item{cdhit_memory}{Integer. CD-HIT \code{-M} memory limit (MB). Use \code{0} for unlimited. Default: \code{0}.} +\item{cdhit_memory}{Integer. CD-HIT \code{-M} memory limit in MB. Use \code{0} for +unlimited. Default: \code{0}.} -\item{cdhit_extra_args}{Character vector. Extra arguments forwarded to \code{cd-hit} -(e.g., \code{c("-g","1")}). Default: \code{c("-g","1")}.} +\item{cdhit_extra_args}{Character vector. Extra arguments forwarded to +\code{cd-hit}. Default: \code{c("-g", "1")}.} -\item{cdhit_output_prefix}{Character. Prefix for CD-HIT output files. Default: \code{"cdhit_out"}.} +\item{cdhit_output_prefix}{Character. Prefix for CD-HIT output files. +Default: \code{"cdhit_out"}.} -\item{ipr_appl}{Character vector. InterProScan applications to run; typically \code{c("Pfam")}. -Default: \code{c("Pfam")}.} +\item{hmmer_databases}{Character vector. HMMER annotation databases to run. +Default: \code{c("Pfam", "COG", "AMRFinder", "DefenseCas")}.} -\item{ipr_threads_unused}{Deprecated/unused. Kept for backward compatibility; ignored.} +\item{hmmer_db_dir}{Character or \code{NULL}. Directory containing the shared HMMER +database cache. If \code{NULL}, uses the amRdata user cache.} -\item{ipr_version}{Character. InterProScan image tag (e.g., \code{"5.76-107.0"}). Default: \code{"5.76-107.0"}.} +\item{hmmer_docker_image}{Character. Docker image containing HMMER. +Default: \code{"staphb/hmmer"}.} -\item{ipr_dest_dir}{Character. Local destination for InterProScan data bundle -(used by \code{.checkInterProData()}). Default: \code{"inst/extdata/interpro"}.} +\item{hmmer_num_splits}{Integer. Number of protein-sequence chunks for HMMER. +Default: \code{8}.} -\item{ipr_platform}{Character. Docker platform string for InterProScan containers, -e.g., \code{"linux/amd64"}. Default: \code{"linux/amd64"}.} +\item{hmmer_workers}{Integer. Number of parallel HMMER workers. Default: \code{8}.} -\item{auto_prepare_data}{Logical. If \code{TRUE}, ensure InterProScan data are present -(download/verify if missing). Default: \code{TRUE}.} - -\item{ref_file_path}{Character. Directory containing reference TSVs used by cleanData() -for metadata harmonization (e.g., \code{"data_raw/"}). \strong{Required}; defaults to \code{"data_raw/"}.} +\item{ref_file_path}{Character. Directory containing reference TSVs used by +\code{\link[=cleanMetaData]{cleanMetaData()}} and \code{\link[=cleanData]{cleanData()}} for metadata harmonization. +Default: \code{"data_raw/"}.} \item{verbose}{Logical. Print progress messages. Default: \code{TRUE}.} } \value{ Invisibly returns a list with: \itemize{ -\item \code{duckdb_path} – input DuckDB path -\item \code{panaroo_output} – path to the selected Panaroo output directory used for import -\item \code{parquet_duckdb_path} – absolute path to the created Parquet-backed DuckDB +\item \code{duckdb_path} - input DuckDB path +\item \code{panaroo_output} - path to the selected Panaroo output directory used for import +\item \code{parquet_duckdb_path} - absolute path to the created Parquet-backed DuckDB } } \description{ \code{runDataProcessing()} orchestrates the complete feature-extraction pipeline for a -BV-BRC selection, starting from a \strong{per-selection DuckDB} (created by -\code{\link[=prepareGenomes]{prepareGenomes()}} and populated by downstream steps). It: +BV-BRC selection, starting from a \strong{per-selection DuckDB} created by +\code{\link[=prepareGenomes]{prepareGenomes()}} and populated by downstream genome processing steps. It: \enumerate{ \item Runs \strong{Panaroo} to build the pangenome and writes gene/struct outputs into DuckDB. \item Runs \strong{CD-HIT} to cluster proteins and writes protein outputs into DuckDB. -\item Runs \strong{InterProScan} (Pfam) to annotate protein domains and writes domain outputs into DuckDB. +\item Runs \strong{HMMER} against the requested protein databases and writes annotation +tables into DuckDB. \item \strong{Cleans BV-BRC metadata} (drug names/classes, countries, years) and -exports all feature/metadata tables as \strong{compressed Parquet} files, then creates -a \strong{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 \strong{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. } \details{ \strong{Docker & Platform Notes} \itemize{ -\item All heavy tools (Panaroo, CD-HIT, InterProScan) run inside Docker containers. -\item On Apple Silicon/ARM hosts, images are forced to \verb{--platform linux/amd64} to ensure compatibility. -\item Ensure Docker Desktop is running and has sufficient memory/CPUs configured. +\item Panaroo, CD-HIT, and HMMER run inside Docker containers. +\item HMMER databases are stored separately from individual bug directories and +are reused across datasets unless a custom \code{hmmer_db_dir} is supplied. +\item Ensure Docker Desktop is running and has sufficient memory and CPU resources. } \strong{Input Requirements} \itemize{ -\item The \code{duckdb_path} must reference a per-selection DuckDB that contains: -\code{files} (paths to \code{.gff}, \code{.fna}, \code{.PATRIC.faa}), -\code{filtered} (genomes selected for download/filtering), and -BV-BRC metadata tables written by earlier steps. +\item \code{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. } \strong{Outputs & Side Effects} \itemize{ -\item Writes tool-specific intermediate outputs under \code{output_path} (e.g., \verb{panaroo_out_*}, CD-HIT files). -\item Writes Parquet files to \code{output_path}: -\code{gene_count.parquet}, \code{protein_count.parquet}, \code{domain_count.parquet}, \code{struct.parquet}, -\code{gene_names.parquet}, \code{protein_names.parquet}, \code{domain_names.parquet}, -\code{gene_seqs.parquet}, \code{protein_seqs.parquet}, \code{genome_gene_protein.parquet}, -\code{metadata.parquet}, \code{amr_phenotype.parquet}, \code{genome_data.parquet}, \code{original_metadata.parquet}. -\item Creates a new Parquet-backed DuckDB (\verb{*_parquet.duckdb}) with read-only views pointing to those Parquets. +\item Writes tool-specific intermediate outputs under \code{output_path}. +\item Writes feature and metadata Parquet files under \code{output_path}. +\item Creates a new Parquet-backed DuckDB (\verb{*_parquet.duckdb}) with read-only views +over the generated Parquet files. +\item Records processing parameters, software versions, database selections, and +other provenance information in the dataset manifest. } \strong{Threading} \itemize{ -\item \code{threads} is a shared budget; each stage uses a portion or all of it. -\item InterProScan can be memory-intensive; on laptops, single-container mode is used internally. +\item \code{threads} provides the shared CPU budget for the major processing stages. +\item Panaroo, CD-HIT, and HMMER allocate that budget according to their respective +stage parameters. } } \section{Pipeline Steps}{ \enumerate{ -\item \strong{Panaroo} via runPanaroo2Duckdb() → writes: +\item \strong{Panaroo} via \code{\link[=runPanaroo2Duckdb]{runPanaroo2Duckdb()}} -> writes: \itemize{ -\item \code{gene_count} (genome × gene counts)\cr +\item \code{gene_count} (genome x gene counts)\cr \item \code{gene_names}\cr \item \code{gene_struct} (structural variants)\cr \item \code{gene_ref_seq}, \code{genome_gene_protein} } -\item \strong{CD-HIT} via CDHIT2duckdb() (calls internal \code{.runCDHIT()}) → writes: +\item \strong{CD-HIT} via \code{\link[=CDHIT2duckdb]{CDHIT2duckdb()}} -> writes: \itemize{ -\item \code{protein_count} (genome × protein-cluster counts)\cr +\item \code{protein_count} (genome x protein-cluster counts)\cr \item \code{protein_names}\cr -\item \code{protein_cluster_seq} (representative sequences) -} -\item \strong{InterProScan (Pfam)} via domainFromIPR() → writes: -\itemize{ -\item \code{domain_names}\cr -\item \code{domain_count} (genome × domain-family matrix) +\item \code{protein_cluster_seq} (representative sequences)\cr +\item \code{protein_members} } -\item \strong{Metadata cleaning + Parquet export} via cleanData() → writes Parquet -files to \code{output_path}, and builds a \strong{Parquet-backed DuckDB} -(\verb{*_parquet.duckdb}) with views: +\item \strong{HMMER} via the configured HMMER databases -> writes: \itemize{ -\item \code{gene_count}, \code{protein_count}, \code{domain_count}, \code{struct}\cr -\item \code{metadata} (cleaned), plus \code{amr_phenotype}, \code{genome_data}, \code{original_metadata}\cr -\item \code{gene_names}, \code{protein_names}, \code{domain_names}\cr -\item \code{gene_seqs}, \code{protein_seqs}\cr -\item \code{genome_gene_protein} +\item \verb{protein_} annotation tables\cr +\item \verb{protein__count} genome-by-annotation count tables +\item The default databases are \code{Pfam}, \code{COG}, \code{AMRFinder}, and \code{DefenseCas}. } +\item \strong{Metadata cleaning + Parquet export} via \code{\link[=cleanData]{cleanData()}} -> writes +Parquet files to \code{output_path}, and builds a \strong{Parquet-backed DuckDB} +(\verb{*_parquet.duckdb}) with views over those Parquets. } } \examples{ \dontrun{ -# Paths below are illustrative; adapt to your project layout. runDataProcessing( duckdb_path = "data/Shigella_flexneri/Sfl.duckdb", output_path = "data/Shigella_flexneri", @@ -194,11 +201,12 @@ runDataProcessing( ) # After completion: -# data/Shigella_flexneri/Sfl_parquet.duckdb +# data/Shigella_flexneri/Sfl_parquet.duckdb # will contain views over the Parquet files for downstream ML. } } \seealso{ -prepareGenomes(), runPanaroo2Duckdb(), CDHIT2duckdb(), domainFromIPR(), cleanData() +\code{\link[=prepareGenomes]{prepareGenomes()}}, \code{\link[=runPanaroo2Duckdb]{runPanaroo2Duckdb()}}, \code{\link[=CDHIT2duckdb]{CDHIT2duckdb()}}, \code{\link[=cleanMetaData]{cleanMetaData()}}, +\code{\link[=cleanData]{cleanData()}} } diff --git a/man/runPanaroo2Duckdb.Rd b/man/runPanaroo2Duckdb.Rd index 48bec57..38f2b51 100644 --- a/man/runPanaroo2Duckdb.Rd +++ b/man/runPanaroo2Duckdb.Rd @@ -131,6 +131,6 @@ runPanaroo2Duckdb( \item \code{.runPanaroo()} — core Panaroo execution \item \code{.mergePanaroo()} — merge multiple Panaroo batches \item \code{.panaroo2duckdb()} — import Panaroo results into DuckDB -\item \code{\link[=runDataProcessing]{runDataProcessing()}} — full pipeline including CD-HIT & InterProScan +\item \code{\link[=runDataProcessing]{runDataProcessing()}} — full pipeline including CD-HIT & HMMER } } diff --git a/tests/testthat/test-buildDyadFeatureMap.R b/tests/testthat/test-buildDyadFeatureMap.R new file mode 100644 index 0000000..04dedd9 --- /dev/null +++ b/tests/testthat/test-buildDyadFeatureMap.R @@ -0,0 +1,23 @@ +test_that("buildDyadFeatureMap() errors when the DuckDB path does not exist", { + expect_error( + buildDyadFeatureMap(duckdb_path = tempfile(fileext = ".duckdb")) + ) +}) + +test_that("buildDyadFeatureMap() errors when no provenance manifest is found", { + tmp_dir <- file.path(tempdir(), paste0("dyad-test-", as.integer(runif(1, 1, 1e6)))) + dir.create(tmp_dir) + on.exit(unlink(tmp_dir, recursive = TRUE), add = TRUE) + + fake_db <- file.path(tmp_dir, "fake.duckdb") + file.create(fake_db) + + expect_error( + buildDyadFeatureMap(duckdb_path = fake_db), + regexp = "provenance manifest" + ) +}) + +# TODO (Bioconductor): add an end-to-end test that runs buildDyadFeatureMap() +# against a small fixture DuckDB + manifest + annotation Parquets and checks the +# resulting edge list (source/target columns, dyad format, type prefixes).