From 7202f92b07dc6fddefd9caf74dc93375ca880c51 Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Thu, 30 Jul 2026 11:16:32 -0600
Subject: [PATCH 01/17] Update runHMMER.R

---
 R/runHMMER.R | 93 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 93 insertions(+)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index 6d77565..520bcd9 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -1,3 +1,96 @@
+#' Title
+#'
+#' @param db_dir
+#' @param db_name
+#'
+#' @returns
+#'
+#' @export
+#' @examples
+.checkHmmerDatabase <- function(hmmer_db_dir, db_name) {
+
+  dir.create(db_dir, recursive = TRUE, showWarnings = FALSE)
+
+  dbs <- list(
+
+    Pfam = list(
+      hmm = file.path(db_dir, "Pfam-A", "Pfam-A.hmm"),
+      url = "https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz"
+    ),
+
+    COG = list(
+      hmm = file.path(db_dir, "COG", "COG.hmm"),
+      url = "http://boabio.belozersky.msu.ru/media/COG_database2024.zip"
+    ),
+
+    AMRFinder = list(
+      hmm = file.path(db_dir, "AMRFinder", "AMRFinder.hmm"),
+      url = "https://ftp.ncbi.nlm.nih.gov/hmm/NCBIfam-AMRFinder/latest/NCBIfam-AMRFinder.HMM.tar.gz"
+    )
+  )
+
+  db <- dbs[[db_name]]
+
+  if (!file.exists(db$hmm)) {
+
+    if (is.na(db$url)) {
+      stop(
+        db_name,
+        " database not found in ",
+        db_dir,
+        ". Please place ",
+        basename(db$hmm),
+        " in this directory."
+      )
+    }
+
+    message("Downloading ", db_name, " database")
+
+    tmp <- tempfile(fileext = ".gz")
+
+    utils::download.file(
+      db$url,
+      tmp,
+      mode = "wb"
+    )
+
+    R.utils::gunzip(
+      tmp,
+      destname = db$hmm,
+      overwrite = TRUE,
+      remove = FALSE
+    )
+  }
+
+  pressed_files <- paste0(
+    db$hmm,
+    c(".h3m", ".h3i", ".h3f", ".h3p")
+  )
+
+  if (!all(file.exists(pressed_files))) {
+
+    message("Running hmmpress on ", basename(db$hmm))
+
+    output <- system2(
+      "hmmpress",
+      db$hmm,
+      stdout = TRUE,
+      stderr = TRUE
+    )
+
+    if (!all(file.exists(pressed_files))) {
+      stop(
+        "hmmpress failed for ",
+        db$hmm,
+        "\n",
+        paste(output, collapse = "\n")
+      )
+    }
+  }
+
+  normalizePath(db$hmm)
+}
+
 #' Write a data frame to a compressed Parquet file
 #'
 #' @param df A data frame or tibble to write.

From 5b2de7b2413c530f169468f1a087cfa282a34cd9 Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Thu, 30 Jul 2026 15:07:56 -0600
Subject: [PATCH 02/17] Modify HMM download and press

---
 R/runHMMER.R | 379 +++++++++++++++++++++++++++++++++++++++++----------
 1 file changed, 307 insertions(+), 72 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index 520bcd9..41d33f2 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -7,90 +7,232 @@
 #'
 #' @export
 #' @examples
-.checkHmmerDatabase <- function(hmmer_db_dir, db_name) {
+.prepareHmmerDatabases <- function(
+    hmmer_db_dir,
+    databases = c("Pfam", "COG", "AMRFinder"),
+    docker_image = "staphb/hmmer"
+  ) {
+  
+  hmmer_db_dir <- path.expand(hmmer_db_dir)
 
-  dir.create(db_dir, recursive = TRUE, showWarnings = FALSE)
+  options(timeout = max(3600, getOption("timeout")))
 
   dbs <- list(
 
     Pfam = list(
-      hmm = file.path(db_dir, "Pfam-A", "Pfam-A.hmm"),
-      url = "https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz"
+      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(
-      hmm = file.path(db_dir, "COG", "COG.hmm"),
-      url = "http://boabio.belozersky.msu.ru/media/COG_database2024.zip"
+      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(
-      hmm = file.path(db_dir, "AMRFinder", "AMRFinder.hmm"),
-      url = "https://ftp.ncbi.nlm.nih.gov/hmm/NCBIfam-AMRFinder/latest/NCBIfam-AMRFinder.HMM.tar.gz"
+      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"
     )
   )
 
-  db <- dbs[[db_name]]
+  dbs <- dbs[databases]
 
-  if (!file.exists(db$hmm)) {
+  db_paths <- list()
 
-    if (is.na(db$url)) {
-      stop(
-        db_name,
-        " database not found in ",
-        db_dir,
-        ". Please place ",
-        basename(db$hmm),
-        " in this directory."
+  for (db_name in names(dbs)) {
+
+    db <- dbs[[db_name]]
+
+    dir.create(
+      db$dir,
+      recursive = TRUE,
+      showWarnings = FALSE
+    )
+
+    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
+          )
+
+          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
       )
     }
 
-    message("Downloading ", db_name, " database")
+   if (db_name == "AMRFinder") {
 
-    tmp <- tempfile(fileext = ".gz")
+  hmm_files <- list.files(
+    db$dir,
+    pattern = "\\.hmm$",
+    recursive = TRUE,
+    full.names = TRUE,
+    ignore.case = TRUE
+  )
 
-    utils::download.file(
-      db$url,
-      tmp,
-      mode = "wb"
-    )
+  hmm_file <- file.path(
+    db$dir,
+    "AMRFinder.hmm"
+  )
+
+  source_hmms <- setdiff(
+    normalizePath(hmm_files),
+    normalizePath(hmm_file, mustWork = FALSE)
+  )
+
+  if (!file.exists(hmm_file)) {
 
-    R.utils::gunzip(
-      tmp,
-      destname = db$hmm,
-      overwrite = TRUE,
-      remove = FALSE
+    message(
+      "Combining ",
+      length(source_hmms),
+      " AMRFinder HMM files"
     )
-  }
 
-  pressed_files <- paste0(
-    db$hmm,
-    c(".h3m", ".h3i", ".h3f", ".h3p")
-  )
+    file.create(hmm_file)
 
-  if (!all(file.exists(pressed_files))) {
+    for (f in sort(source_hmms)) {
+      file.append(hmm_file, f)
+    }
+  }
 
-    message("Running hmmpress on ", basename(db$hmm))
+} else {
+      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
+        )
+      }
+      
+      hmm_file <- hmm_files[1]
+    }
 
-    output <- system2(
-      "hmmpress",
-      db$hmm,
-      stdout = TRUE,
-      stderr = TRUE
+    pressed_files <- paste0(
+      hmm_file,
+      c(
+        ".h3m",
+        ".h3i",
+        ".h3f",
+        ".h3p"
+      )
     )
 
     if (!all(file.exists(pressed_files))) {
-      stop(
-        "hmmpress failed for ",
-        db$hmm,
-        "\n",
-        paste(output, collapse = "\n")
+
+      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
+    )
   }
 
-  normalizePath(db$hmm)
+  db_paths
 }
-
 #' Write a data frame to a compressed Parquet file
 #'
 #' @param df A data frame or tibble to write.
@@ -107,19 +249,37 @@
   )
 }
 
+#' Title
+#'
+#' @param duckdb_path
+#' @param output_path
+#' @param threads
+#' @param hmmer_db_dir
+#' @param databases
+#' @param docker_image
+#' @param split_jobs
+#' @param num_of_splits
+#' @param n_workers
+#'
+#' @returns
+#'
+#' @export
+#' @examples
 .runHMMER <- function(duckdb_path,
                       output_path,
                       threads = 8L,
-                      database_path,
+                      hmmer_db_dir,
+                      databases = c("Pfam", "COG", "AMRFinder"),
                       docker_image = "staphb/hmmer",
                       split_jobs = TRUE,
                       num_of_splits = 20L,
-                      n_workers = 4L) {
+                      n_workers = 4L
+                    ) {
   # 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.")
   }
-
+  
   duckdb_path <- .docker_path(duckdb_path)
   if (missing(output_path) || output_path %in% c(".", "results", "results/")) {
     output_path <- dirname(duckdb_path)
@@ -133,8 +293,14 @@
   prot_seqs <- DBI::dbReadTable(con, "protein_cluster_seq") |>
     tibble::as_tibble()
 
-  # derive a clean label from the database filename
-  database <- tools::file_path_sans_ext(basename(database_path))
+  # database paths 
+ db_paths <- .prepareHmmerDatabases(
+  hmmer_db_dir = hmmer_db_dir,
+  databases = databases,
+  docker_image = docker_image
+)
+
+db_paths <- db_paths[databases]
 
   # clamp splits to the number of sequences available
   chunk_count <- min(as.integer(num_of_splits), nrow(prot_seqs))
@@ -154,7 +320,7 @@
 
   job_list <- expand.grid(
     chunk = sprintf("%02d", seq_len(chunk_count)),
-    db = database,
+    db = databases,
     stringsAsFactors = FALSE
   ) |>
     dplyr::mutate(
@@ -169,6 +335,7 @@
     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"
@@ -178,13 +345,18 @@
     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,
       "hmmscan",
-      "--cpu", as.character(threads),
+      "--cpu", as.character(threads_per_job),
       "--tblout", .to_container(hmmer_output, mount_host, mount_cont),
       db_cont_path,
       .to_container(hmmer_input, mount_host, mount_cont)
@@ -219,27 +391,90 @@
     hmmer_tbl_filename
   }
 
-  hmmer_param <- BiocParallel::SnowParam(workers = max(1L, n_workers))
-  parquet_files <- BiocParallel::bpmapply(
-    FUN = .runHmmerJob,
-    JOB_NAME = job_list$JOB_NAME,
-    FASTA = job_list$FASTA,
-    DB = job_list$DB,
-    SIMPLIFY = TRUE,
-    USE.NAMES = FALSE,
-    BPPARAM = hmmer_param
+  future::plan(
+  future::multisession,
+  workers = max(1L, n_workers)
+)
+
+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]
+    )
+  }
+)
+
+future::plan(future::sequential)
+  
+  parquet_tbl <- tibble::tibble(
+  parquet = parquet_files,
+  db = job_list$DB
+)
+  
+ final_parquets <- list()
+
+for (database_name in databases) {
+
+  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()
+
+  final_parquet <- file.path(
+    output_path,
+    paste0(
+      "protein_",
+      database_name,
+      ".parquet"
+    )
+  )
+
+  .write_compressed_parquet(
+    combined_tbl,
+    final_parquet
   )
 
-  final_parquet <- file.path(output_path, paste0("protein_", database, ".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)
+  # purrr::map(parquet_files, arrow::read_parquet) |>
+  #   dplyr::bind_rows() |>
+  #   .write_compressed_parquet(final_parquet)
 
-  message("Combined parquet written.")
+  # message("Combined parquet written.")
 
-  arrow::read_parquet(final_parquet) |>
-    DBI::dbWriteTable(conn = con, name = tools::file_path_sans_ext(basename(final_parquet)), overwrite = TRUE)
+  # arrow::read_parquet(final_parquet) |>
+  #   DBI::dbWriteTable(conn = con, name = tools::file_path_sans_ext(basename(final_parquet)), overwrite = TRUE)
 }
 
 

From fe9bc573584316b572bcce641e3233762a1e729d Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Thu, 30 Jul 2026 18:24:46 -0600
Subject: [PATCH 03/17] default HMM database, start runHMMER from checking for
 database existence.

---
 R/runHMMER.R | 144 ++++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 113 insertions(+), 31 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index 41d33f2..15f8477 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -583,52 +583,134 @@ invisible(final_parquets)
 #' )
 #' }
 #'
-#' @export
-proteinAnnotations2Duckdb <- function(annotated_parquet, duckdb_path) {
-  annotated_parquet <- .docker_path(annotated_parquet)
+#' @internal
+.proteinAnnotations2Duckdb <- function(
+    duckdb_path,
+    databases = c("Pfam", "COG", "AMRFinder")
+) {
+
   duckdb_path <- .docker_path(duckdb_path)
 
-  # derive table name from the annotation filename stem
-  database <- tools::file_path_sans_ext(basename(annotated_parquet))
+  con <- DBI::dbConnect(
+    duckdb::duckdb(),
+    duckdb_path
+  )
 
-  con <- DBI::dbConnect(duckdb::duckdb(), duckdb_path)
-  on.exit(try(DBI::dbDisconnect(con, shutdown = FALSE), silent = TRUE), add = TRUE)
+  on.exit(
+    try(
+      DBI::dbDisconnect(
+        con,
+        shutdown = FALSE
+      ),
+      silent = TRUE
+    ),
+    add = TRUE
+  )
 
-  protein_long <- DBI::dbReadTable(con, "protein_count") |>
+  protein_long <- DBI::dbReadTable(
+    con,
+    "protein_count"
+  ) |>
     tibble::as_tibble() |>
     tidyr::pivot_longer(
       cols = -genome_id,
       names_to = "query_name",
       values_to = "count"
     ) |>
-    dplyr::filter(count > 0)
+    dplyr::filter(count > 0) |>
+    dplyr::mutate(
+      query_name = stringr::str_replace(
+        query_name,
+        "^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."
+      )
 
-  annotation <- arrow::read_parquet(annotated_parquet)
+      next
+    }
 
-  genome_annot_matrix <- protein_long |>
-    # protein IDs are stored with "." separator in DuckDB but "|" in HMMER output
-    dplyr::mutate(query_name = stringr::str_replace(query_name, "^fig\\.", "fig|")) |>
-    dplyr::inner_join(
-      dplyr::select(annotation, name, query_name),
-      by = "query_name"
+    message(
+      "Processing ",
+      annotation_table
+    )
+
+    annotation <- DBI::dbReadTable(
+      con,
+      annotation_table
     ) |>
-    dplyr::group_by(genome_id, name) |>
-    dplyr::summarise(count = sum(count), .groups = "drop") |>
-    tidyr::pivot_wider(
-      names_from = name,
-      values_from = count,
-      values_fill = 0
+      tibble::as_tibble()
+
+    genome_annot_matrix <- protein_long |>
+      dplyr::inner_join(
+        dplyr::select(
+          annotation,
+          name,
+          query_name
+        ),
+        by = "query_name"
+      ) |>
+      dplyr::group_by(
+        genome_id,
+        name
+      ) |>
+      dplyr::summarise(
+        count = sum(count),
+        .groups = "drop"
+      ) |>
+      tidyr::pivot_wider(
+        names_from = name,
+        values_from = count,
+        values_fill = 0
+      )
+
+    count_table <- paste0(
+      annotation_table,
+      "_count"
     )
 
-  count_path <- file.path(dirname(duckdb_path), paste0(database, "_count.parquet"))
-  arrow::write_parquet(genome_annot_matrix, count_path)
+    count_path <- file.path(
+      dirname(duckdb_path),
+      paste0(
+        count_table,
+        ".parquet"
+      )
+    )
 
-  DBI::dbWriteTable(
-    conn = con,
-    name = tools::file_path_sans_ext(basename(count_path)),
-    value = genome_annot_matrix,
-    overwrite = TRUE
-  )
+    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_path)
+  invisible(count_paths)
 }

From 239d71634501a9fbc1e16dae9f2045d1885d1432 Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Tue, 4 Aug 2026 16:17:08 -0600
Subject: [PATCH 04/17] Hmmering

---
 R/runHMMER.R | 615 ++++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 555 insertions(+), 60 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index 15f8477..81dd06b 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -10,10 +10,11 @@
 .prepareHmmerDatabases <- function(
     hmmer_db_dir,
     databases = c("Pfam", "COG", "AMRFinder"),
-    docker_image = "staphb/hmmer"
-  ) {
-  
-  hmmer_db_dir <- path.expand(hmmer_db_dir)
+    docker_image = "staphb/hmmer",
+    hmmer_db_url = NULL
+) {
+
+  hmmer_db_dir <- normalizePath(hmmer_db_dir)
 
   options(timeout = max(3600, getOption("timeout")))
 
@@ -41,6 +42,59 @@
     )
   )
 
+  ## 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()
@@ -79,12 +133,19 @@
       )
 
       switch(
+
         db$type,
 
         gz = {
+
           hmm_file <- file.path(
             db$dir,
-            db$hmm_name
+            db$hmm_name %||% basename(
+              sub("\\.gz$",
+                  "",
+                  basename(db$url),
+                  ignore.case = TRUE)
+            )
           )
 
           R.utils::gunzip(
@@ -96,6 +157,7 @@
         },
 
         zip = {
+
           utils::unzip(
             zipfile = tmp,
             exdir = db$dir
@@ -103,6 +165,7 @@
         },
 
         `tar.gz` = {
+
           utils::untar(
             tarfile = tmp,
             exdir = db$dir
@@ -119,68 +182,57 @@
       )
     }
 
-   if (db_name == "AMRFinder") {
-
-  hmm_files <- list.files(
-    db$dir,
-    pattern = "\\.hmm$",
-    recursive = TRUE,
-    full.names = TRUE,
-    ignore.case = TRUE
-  )
-
-  hmm_file <- file.path(
-    db$dir,
-    "AMRFinder.hmm"
-  )
+      hmm_files <- list.files(
+      db$dir,
+      pattern = "\\.hmm$",
+      recursive = TRUE,
+      full.names = TRUE,
+      ignore.case = TRUE
+    )
 
-  source_hmms <- setdiff(
-    normalizePath(hmm_files),
-    normalizePath(hmm_file, mustWork = FALSE)
-  )
+    if (length(hmm_files) == 0) {
 
-  if (!file.exists(hmm_file)) {
+      stop(
+        "No .hmm file found for ",
+        db_name
+      )
 
-    message(
-      "Combining ",
-      length(source_hmms),
-      " AMRFinder HMM files"
-    )
+    } else if (length(hmm_files) == 1) {
 
-    file.create(hmm_file)
+      hmm_file <- hmm_files[1]
 
-    for (f in sort(source_hmms)) {
-      file.append(hmm_file, f)
-    }
-  }
+    } else {
 
-} else {
-      hmm_files <- list.files(
+      hmm_file <- file.path(
         db$dir,
-        pattern = "\\.hmm$",
-        recursive = TRUE,
-        full.names = TRUE,
-        ignore.case = TRUE
+        paste0(db_name, ".hmm")
       )
-      
-      if (length(hmm_files) == 0) {
-        stop(
-          "No .hmm file found for ",
+
+      source_hmms <- setdiff(
+        normalizePath(hmm_files),
+        normalizePath(hmm_file, mustWork = FALSE)
+      )
+
+      if (!file.exists(hmm_file)) {
+
+        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)
+        }
       }
-      
-      hmm_file <- hmm_files[1]
     }
 
     pressed_files <- paste0(
       hmm_file,
-      c(
-        ".h3m",
-        ".h3i",
-        ".h3f",
-        ".h3p"
-      )
+      c(".h3m", ".h3i", ".h3f", ".h3p")
     )
 
     if (!all(file.exists(pressed_files))) {
@@ -196,16 +248,10 @@
           "run",
           "--rm",
           "-v",
-          paste0(
-            dirname(hmm_file),
-            ":/db"
-          ),
+          paste0(dirname(hmm_file), ":/db"),
           docker_image,
           "hmmpress",
-          file.path(
-            "/db",
-            basename(hmm_file)
-          )
+          file.path("/db", basename(hmm_file))
         ),
         stdout = TRUE,
         stderr = TRUE
@@ -233,6 +279,7 @@
 
   db_paths
 }
+
 #' Write a data frame to a compressed Parquet file
 #'
 #' @param df A data frame or tibble to write.
@@ -293,6 +340,10 @@
   prot_seqs <- DBI::dbReadTable(con, "protein_cluster_seq") |>
     tibble::as_tibble()
 
+  if(is.null(hmmer_db_dir)) {
+  hmmer_db_dir <- output_path
+  }
+
   # database paths 
  db_paths <- .prepareHmmerDatabases(
   hmmer_db_dir = hmmer_db_dir,
@@ -714,3 +765,447 @@ invisible(final_parquets)
 
   invisible(count_paths)
 }
+
+# Annotate the proteins using defense finder and cas finder
+ #' Annotate proteins using DefenseFinder + CasFinder HMMs
+#'
+#' Downloads DefenseFinder and CasFinder model repositories,
+#' extracts all HMMs located within profile directories,
+#' concatenates them into a single DefenseCas.hmm database,
+#' runs hmmpress, performs HMMER annotation against proteins
+#' stored in DuckDB, and stores the results in DuckDB.
+#'
+#' @param defense_db_dir Directory used to store downloaded models.
+#' @param docker_image Docker image containing HMMER.
+#' @param duckdb_path Path to duckdb database.
+#' @param output_path Output directory.
+#' @param threads Number of HMMER threads.
+#' @param split_jobs Split sequences into chunks.
+#' @param num_of_splits Number of fasta chunks.
+#' @param n_workers Parallel workers.
+#'
+#' @returns Invisibly returns parquet file path.
+#'
+#' @export
+#' Annotate proteins using DefenseFinder + CasFinder models
+#'
+#' @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
+#' @export
+.defenseHMMER <- function(
+    defense_db_dir,
+    docker_image = "staphb/hmmer",
+    duckdb_path = "inst/extdata/Sfl.duckdb",
+    output_path = NULL,
+    threads = 8L
+) {
+
+  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)) {
+
+    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)) {
+
+    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"
+    ]
+
+    hmm_files <- unique(
+      unlist(
+        lapply(
+          profile_dirs,
+          function(x) {
+
+            list.files(
+              x,
+              pattern = "\\.hmm$",
+              recursive = TRUE,
+              full.names = TRUE,
+              ignore.case = TRUE
+            )
+
+          }
+        )
+      )
+    )
+
+    if (length(hmm_files) == 0) {
+
+      stop(
+        "No HMM files found for ",
+        db_name
+      )
+    }
+
+    hmm_files <- hmm_files[
+      vapply(
+        hmm_files,
+        function(f) {
+
+          first_line <- tryCatch(
+            readLines(
+              f,
+              n = 1,
+              warn = FALSE
+            ),
+            error = function(e) ""
+          )
+
+          grepl(
+            "^HMMER",
+            first_line
+          )
+
+        },
+        logical(1)
+      )
+    ]
+
+    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))) {
+
+      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_seq"
+  ) |>
+    tibble::as_tibble()
+
+  fasta_file <- file.path(
+    output_path,
+    "protein_DefenseCas.faa"
+  )
+
+  readr::write_lines(
+    paste0(
+      ">",
+      prot_seqs$name,
+      "\n",
+      prot_seqs$sequence
+    ),
+    fasta_file
+  )
+
+  ####################################################################
+  # run hmmscan separately
+  ####################################################################
+
+  databases <- list(
+    DefenseFinder = defense_hmm,
+    CasFinder = cas_hmm
+  )
+
+  all_hits <- list()
+
+  for (db_name in names(databases)) {
+
+    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,
+        "hmmscan",
+        "--cpu",
+        as.character(threads),
+        "--tblout",
+        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(
+        "hmmscan failed for ",
+        db_name,
+        "\n",
+        paste(output,
+              collapse = "\n")
+      )
+    }
+
+    hits <- .parseHMMEROutput(
+      tbl_file
+    ) |>
+      dplyr::select(
+        name,
+        query_name,
+        description
+      ) |>
+      dplyr::mutate(
+        database = db_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)
+}
\ No newline at end of file

From cdc9fa435ad5d9380a78555e2a39d7498042ab27 Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Tue, 4 Aug 2026 16:54:34 -0600
Subject: [PATCH 05/17] change the stop with warning

---
 R/runHMMER.R | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index 81dd06b..b9e8efe 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -1152,16 +1152,17 @@ invisible(final_parquets)
       stderr = TRUE
     )
 
-    if (!file.exists(tbl_file)) {
+   if (!file.exists(tbl_file)) {
 
-      stop(
-        "hmmscan failed for ",
-        db_name,
-        "\n",
-        paste(output,
-              collapse = "\n")
-      )
-    }
+  warning(
+    "hmmscan failed for ",
+    db_name,
+    "\n",
+    paste(output, collapse = "\n")
+  )
+
+  next
+}
 
     hits <- .parseHMMEROutput(
       tbl_file

From 9fe0af081cbf5d3df9dc5346f0857e6e39b2d1e7 Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Tue, 4 Aug 2026 17:16:33 -0600
Subject: [PATCH 06/17] protein_cluster_seq

---
 R/runHMMER.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index b9e8efe..9f6511e 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -1076,7 +1076,7 @@ invisible(final_parquets)
 
   prot_seqs <- DBI::dbReadTable(
     con,
-    "protein_seq"
+    "protein_cluster_seq"
   ) |>
     tibble::as_tibble()
 

From a8cc93ffbbc12153ab100988e06b4ce8843866cc Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Tue, 4 Aug 2026 17:22:54 -0600
Subject: [PATCH 07/17] find the HMMs with old version and remove

---
 R/runHMMER.R | 56 ++++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 52 insertions(+), 4 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index 9f6511e..ffddcfa 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -1,3 +1,23 @@
+.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
+}
+
 #' Title
 #'
 #' @param db_dir
@@ -208,10 +228,38 @@
         paste0(db_name, ".hmm")
       )
 
-      source_hmms <- setdiff(
-        normalizePath(hmm_files),
-        normalizePath(hmm_file, mustWork = FALSE)
-      )
+     source_hmms <- setdiff(
+  normalizePath(hmm_files),
+  normalizePath(hmm_file, mustWork = FALSE)
+)
+
+valid_hmms <- vapply(
+  source_hmms,
+  .isValidHmmFile,
+  logical(1)
+)
+
+if (any(!valid_hmms)) {
+
+  bad_files <- source_hmms[!valid_hmms]
+
+  warning(
+    "Ignoring ",
+    length(bad_files),
+    " invalid HMM file(s):\n",
+    paste(basename(bad_files), collapse = "\n")
+  )
+
+  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)) {
 

From a9974dadbc9f87b181aff8dbd0ea758c7fdfab0f Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Tue, 4 Aug 2026 17:25:02 -0600
Subject: [PATCH 08/17] Roxygen skeleton for .defenseHMMER

---
 R/runHMMER.R | 23 +----------------------
 1 file changed, 1 insertion(+), 22 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index ffddcfa..e322188 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -814,27 +814,6 @@ invisible(final_parquets)
   invisible(count_paths)
 }
 
-# Annotate the proteins using defense finder and cas finder
- #' Annotate proteins using DefenseFinder + CasFinder HMMs
-#'
-#' Downloads DefenseFinder and CasFinder model repositories,
-#' extracts all HMMs located within profile directories,
-#' concatenates them into a single DefenseCas.hmm database,
-#' runs hmmpress, performs HMMER annotation against proteins
-#' stored in DuckDB, and stores the results in DuckDB.
-#'
-#' @param defense_db_dir Directory used to store downloaded models.
-#' @param docker_image Docker image containing HMMER.
-#' @param duckdb_path Path to duckdb database.
-#' @param output_path Output directory.
-#' @param threads Number of HMMER threads.
-#' @param split_jobs Split sequences into chunks.
-#' @param num_of_splits Number of fasta chunks.
-#' @param n_workers Parallel workers.
-#'
-#' @returns Invisibly returns parquet file path.
-#'
-#' @export
 #' Annotate proteins using DefenseFinder + CasFinder models
 #'
 #' @param defense_db_dir Directory used to store downloaded HMMs
@@ -844,7 +823,7 @@ invisible(final_parquets)
 #' @param threads Number of HMMER threads
 #'
 #' @returns Path to annotation parquet
-#' @export
+#' @keywords internal
 .defenseHMMER <- function(
     defense_db_dir,
     docker_image = "staphb/hmmer",

From 9ff53b6740561725a07346cba539c038ac01adfe Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Tue, 4 Aug 2026 17:38:45 -0600
Subject: [PATCH 09/17] add the file validation check to .defenseHMMER

---
 R/runHMMER.R | 46 +++++++++++++++++++++++++---------------------
 1 file changed, 25 insertions(+), 21 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index e322188..6a2107c 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -968,29 +968,33 @@ invisible(final_parquets)
       )
     }
 
-    hmm_files <- hmm_files[
-      vapply(
-        hmm_files,
-        function(f) {
-
-          first_line <- tryCatch(
-            readLines(
-              f,
-              n = 1,
-              warn = FALSE
-            ),
-            error = function(e) ""
-          )
+    valid_hmms <- vapply(
+  hmm_files,
+  .isValidHmmFile,
+  logical(1)
+)
 
-          grepl(
-            "^HMMER",
-            first_line
-          )
+if (any(!valid_hmms)) {
 
-        },
-        logical(1)
-      )
-    ]
+  warning(
+    "Ignoring ",
+    sum(!valid_hmms),
+    " invalid HMM file(s):\n",
+    paste(
+      basename(hmm_files[!valid_hmms]),
+      collapse = "\n"
+    )
+  )
+}
+
+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,

From afd1132fa305df7cedc39315e937b41617b9a6db Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Tue, 4 Aug 2026 17:59:55 -0600
Subject: [PATCH 10/17] Enhance documentation for HMMER functions with detailed
 parameter descriptions and examples

---
 R/runHMMER.R | 29 ++++++++++++++++++++---------
 1 file changed, 20 insertions(+), 9 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index 6a2107c..cecf251 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -1,3 +1,11 @@
+#' Validate a HMM file has old HMMER3 format and remove them
+#'
+#' @param hmm_file
+#'
+#' @returns
+#'
+#' @keywords internal
+#' @examples
 .isValidHmmFile <- function(hmm_file) {
 
   lines <- tryCatch(
@@ -18,14 +26,16 @@
   starts_ok && ends_ok
 }
 
-#' Title
+#' Download and prepare HMMER databases for generating new file types. 
 #'
-#' @param db_dir
-#' @param db_name
+#' @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
+#' @returns A list of paths to the database hmm files. 
 #'
-#' @export
+#' @keywords internal
 #' @examples
 .prepareHmmerDatabases <- function(
     hmmer_db_dir,
@@ -344,7 +354,7 @@ if (length(source_hmms) == 0) {
   )
 }
 
-#' Title
+#' Wrapper for preparing HMM databases and running HMMER on protein sequences from duckdb and writing them.
 #'
 #' @param duckdb_path
 #' @param output_path
@@ -358,7 +368,7 @@ if (length(source_hmms) == 0) {
 #'
 #' @returns
 #'
-#' @export
+#' @keywords internal
 #' @examples
 .runHMMER <- function(duckdb_path,
                       output_path,
@@ -682,7 +692,7 @@ invisible(final_parquets)
 #' )
 #' }
 #'
-#' @internal
+#' @keywords internal
 .proteinAnnotations2Duckdb <- function(
     duckdb_path,
     databases = c("Pfam", "COG", "AMRFinder")
@@ -814,7 +824,8 @@ invisible(final_parquets)
   invisible(count_paths)
 }
 
-#' Annotate proteins using DefenseFinder + CasFinder models
+#' 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

From 8d4cd669c1133675cbed8776846fc06667d54177 Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Wed, 5 Aug 2026 14:13:28 -0600
Subject: [PATCH 11/17] hmmscan to hmmsearch

---
 R/runHMMER.R | 154 ++++++++++++++++++++++++++++++++++++---------------
 1 file changed, 109 insertions(+), 45 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index cecf251..73750ed 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -354,6 +354,51 @@ if (length(source_hmms) == 0) {
   )
 }
 
+#' Title
+#'
+#' @param hmm_file
+#'
+#' @returns
+#'
+#' @export
+#' @examples
+.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+")
+    )
+  )
+}
+
+
 #' Wrapper for preparing HMM databases and running HMMER on protein sequences from duckdb and writing them.
 #'
 #' @param duckdb_path
@@ -398,6 +443,9 @@ if (length(source_hmms) == 0) {
   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
   }
@@ -439,7 +487,7 @@ db_paths <- db_paths[databases]
     ) |>
     dplyr::select(JOB_NAME, FASTA, DB)
 
-  .runHmmerJob <- function(JOB_NAME, FASTA, DB) {
+  .runHmmerJob <- function(JOB_NAME, FASTA, DB, Total_proteins) {
     hmmer_input <- file.path(output_path, FASTA)
     hmmer_output <- file.path(output_path, paste0(JOB_NAME, ".tbl"))
 
@@ -464,31 +512,34 @@ db_paths <- db_paths[databases]
       "-v", paste0(mount_host, ":", mount_cont),
       "-v", paste0(db_host_dir, ":", db_cont_dir),
       docker_image,
-      "hmmscan",
+      "hmmsearch",
+      "--notextw",
       "--cpu", as.character(threads_per_job),
-      "--tblout", .to_container(hmmer_output, mount_host, mount_cont),
+      "-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)
     )
 
-    message("Running hmmscan via Docker...")
+    message("Running hmmsearch via Docker...")
     output <- tryCatch(
       {
         system2("docker", args = cmd_args, stdout = TRUE, stderr = TRUE)
       },
       error = function(e) {
-        stop("hmmscan execution failed: ", e$message)
+        stop("hmmsearch execution failed: ", e$message)
       }
     )
 
     if (!file.exists(hmmer_output)) {
-      stop("hmmscan failed: output file not found. Check stderr:\n", paste(output, collapse = "\n"))
+      stop("hmmsearch failed: output file not found. Check stderr:\n", paste(output, collapse = "\n"))
     }
 
-    message("hmmscan completed successfully.")
+    message("hmmsearch completed successfully.")
 
     hmmer_tbl <- .parseHMMEROutput(hmmer_output) |>
-      dplyr::select("name", "query_name", "description")
+      dplyr::select("protein", "query_name")
 
     hmmer_tbl_filename <- file.path(
       dirname(hmmer_output),
@@ -512,7 +563,8 @@ parquet_files <- furrr::future_map_chr(
     .runHmmerJob(
       JOB_NAME = job_list$JOB_NAME[i],
       FASTA = job_list$FASTA[i],
-      DB = job_list$DB[i]
+      DB = job_list$DB[i],
+      Total_proteins = Total_proteins
     )
   }
 )
@@ -540,7 +592,7 @@ for (database_name in databases) {
     db_files,
     arrow::read_parquet
   ) |>
-    dplyr::bind_rows()
+    dplyr::bind_rows() 
 
   final_parquet <- file.path(
     output_path,
@@ -593,13 +645,9 @@ invisible(final_parquets)
 #' 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 `--tblout`.
+#' @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:
-#'   `name`, `accession`, `query_name`, `query_accession`, `sequence_evalue`,
-#'   `sequence_score`, `sequence_bias`, `best_evalue`, `best_score`,
-#'   `best_bias`, `number_exp`, `number_reg`, `number_clu`, `number_ov`,
-#'   `number_env`, `number_dom`, `number_rep`, `number_inc`, `description`.
+#' @return A tibble with 19 columns matching the HMMER per-sequence hit table
 #'
 #' @references Adapted from the rhmmer package
 #'   (<https://github.com/arendsee/rhmmer>).
@@ -612,27 +660,43 @@ invisible(final_parquets)
 #'
 #' @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(
-    name = readr::col_character(),
-    accession = readr::col_character(),
-    query_name = readr::col_character(),
-    query_accession = readr::col_character(),
-    sequence_evalue = readr::col_double(),
-    sequence_score = readr::col_double(),
-    sequence_bias = readr::col_double(),
-    best_evalue = readr::col_double(),
-    best_score = readr::col_double(),
-    best_bias = readr::col_double(),
-    number_exp = readr::col_double(),
-    number_reg = readr::col_integer(),
-    number_clu = readr::col_integer(),
-    number_ov = readr::col_integer(),
-    number_env = readr::col_integer(),
-    number_dom = readr::col_integer(),
-    number_rep = readr::col_integer(),
-    number_inc = readr::col_character(),
-    description = readr::col_character()
-  )
+  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)
 
@@ -660,9 +724,10 @@ invisible(final_parquets)
       col_types = col_types,
       lazy = FALSE,
       progress = FALSE
-    ) |>
-    tidyr::unite(description, description:last_col(), sep = " ")
-  table$description <- gsub("\t", " ", table$description)
+    ) 
+  # |>
+  #   tidyr::unite(description, description:last_col(), sep = " ")
+  # table$description <- gsub("\t", " ", table$description)
 
   table
 }
@@ -1138,7 +1203,7 @@ if (length(hmm_files) == 0) {
   )
 
   ####################################################################
-  # run hmmscan separately
+  # run hmmsearch separately
   ####################################################################
 
   databases <- list(
@@ -1176,7 +1241,7 @@ if (length(hmm_files) == 0) {
         "-v",
         paste0(dirname(hmm_file), ":/db"),
         docker_image,
-        "hmmscan",
+        "hmmsearch",
         "--cpu",
         as.character(threads),
         "--tblout",
@@ -1197,7 +1262,7 @@ if (length(hmm_files) == 0) {
    if (!file.exists(tbl_file)) {
 
   warning(
-    "hmmscan failed for ",
+    "hmmsearch failed for ",
     db_name,
     "\n",
     paste(output, collapse = "\n")
@@ -1210,9 +1275,8 @@ if (length(hmm_files) == 0) {
       tbl_file
     ) |>
       dplyr::select(
-        name,
-        query_name,
-        description
+        protein,
+        query_name
       ) |>
       dplyr::mutate(
         database = db_name

From 895f53a1349dc35fe71c26c435c189c2a9313751 Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Wed, 5 Aug 2026 15:49:03 -0600
Subject: [PATCH 12/17] incorporate description from HMM files in to the final
 HMMER output. co-authored by @epbrenner

---
 R/runHMMER.R | 146 ++++++++++++++++++++++++++++-----------------------
 1 file changed, 80 insertions(+), 66 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index 73750ed..b40942e 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -398,6 +398,81 @@ if (length(source_hmms) == 0) {
   )
 }
 
+#' Title
+#'
+#' @param JOB_NAME
+#' @param FASTA
+#' @param DB
+#' @param Total_proteins
+#'
+#' @returns
+#'
+#' @export
+#' @examples
+.runHmmerJob <- function(JOB_NAME, FASTA, DB, Total_proteins) {
+    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)
+    )
+
+    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"))
+    }
+
+    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.
 #'
@@ -487,70 +562,6 @@ db_paths <- db_paths[databases]
     ) |>
     dplyr::select(JOB_NAME, FASTA, DB)
 
-  .runHmmerJob <- function(JOB_NAME, FASTA, DB, Total_proteins) {
-    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)
-    )
-
-    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"))
-    }
-
-    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
-  }
-
   future::plan(
   future::multisession,
   workers = max(1L, n_workers)
@@ -592,8 +603,11 @@ for (database_name in databases) {
     db_files,
     arrow::read_parquet
   ) |>
-    dplyr::bind_rows() 
-
+    dplyr::bind_rows() |>
+dplyr::left_join(.parse_hmmer_profiles(db_paths[[database_name]]) |>
+  dplyr::select(query_name = profile_name, profile_accession, profile_description),
+by = "query_name")
+  
   final_parquet <- file.path(
     output_path,
     paste0(

From f054debbce216b81e124d7b0893d5efedf0a6c1a Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Wed, 5 Aug 2026 16:18:25 -0600
Subject: [PATCH 13/17] add description

---
 R/runHMMER.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index b40942e..b984540 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -605,7 +605,7 @@ for (database_name in databases) {
   ) |>
     dplyr::bind_rows() |>
 dplyr::left_join(.parse_hmmer_profiles(db_paths[[database_name]]) |>
-  dplyr::select(query_name = profile_name, profile_accession, profile_description),
+  dplyr::select(query_name = profile_name, profile_accession, description = profile_description),
 by = "query_name")
   
   final_parquet <- file.path(

From ee55b367f80b7f22a801e20480bdf9811c0ac58a Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Wed, 5 Aug 2026 16:23:49 -0600
Subject: [PATCH 14/17] update .runHmmerJob()

---
 R/runHMMER.R | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index b984540..ae343ad 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -409,7 +409,11 @@ if (length(source_hmms) == 0) {
 #'
 #' @export
 #' @examples
-.runHmmerJob <- function(JOB_NAME, FASTA, DB, Total_proteins) {
+.runHmmerJob <- function(JOB_NAME, FASTA, DB, Total_proteins, 
+  output_path = NULL, db_paths, 
+  docker_image = "staphb/hmmer", threads = 8L, 
+  n_workers = 4L
+) {
     hmmer_input <- file.path(output_path, FASTA)
     hmmer_output <- file.path(output_path, paste0(JOB_NAME, ".tbl"))
 
@@ -575,7 +579,12 @@ parquet_files <- furrr::future_map_chr(
       JOB_NAME = job_list$JOB_NAME[i],
       FASTA = job_list$FASTA[i],
       DB = job_list$DB[i],
-      Total_proteins = Total_proteins
+      Total_proteins = Total_proteins,
+      output_path = output_path,
+      db_paths = db_paths,
+      docker_image = docker_image,
+      threads = threads,
+      n_workers = n_workers
     )
   }
 )

From 190f84dcd4a4c723b9a6da73cd19d611f36c549b Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Wed, 5 Aug 2026 16:38:31 -0600
Subject: [PATCH 15/17] update .proteinAnnotations2Duckdb() based on the new
 hmmer output parsing

---
 R/runHMMER.R | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index ae343ad..e50cbcc 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -614,7 +614,7 @@ for (database_name in databases) {
   ) |>
     dplyr::bind_rows() |>
 dplyr::left_join(.parse_hmmer_profiles(db_paths[[database_name]]) |>
-  dplyr::select(query_name = profile_name, profile_accession, description = profile_description),
+  dplyr::select(query_name = profile_name, query_accession = profile_accession, description = profile_description),
 by = "query_name")
   
   final_parquet <- file.path(
@@ -811,13 +811,13 @@ invisible(final_parquets)
     tibble::as_tibble() |>
     tidyr::pivot_longer(
       cols = -genome_id,
-      names_to = "query_name",
+      names_to = "protein",
       values_to = "count"
     ) |>
     dplyr::filter(count > 0) |>
     dplyr::mutate(
-      query_name = stringr::str_replace(
-        query_name,
+      protein = stringr::str_replace(
+        protein,
         "^fig\\.",
         "fig|"
       )
@@ -854,24 +854,24 @@ invisible(final_parquets)
       tibble::as_tibble()
 
     genome_annot_matrix <- protein_long |>
-      dplyr::inner_join(
+      dplyr::inner_join( 
+        annotation |>
         dplyr::select(
-          annotation,
-          name,
+          protein,
           query_name
         ),
-        by = "query_name"
+        by = "protein"
       ) |>
       dplyr::group_by(
         genome_id,
-        name
+        query_name
       ) |>
       dplyr::summarise(
         count = sum(count),
         .groups = "drop"
       ) |>
       tidyr::pivot_wider(
-        names_from = name,
+        names_from = query_name,
         values_from = count,
         values_fill = 0
       )

From d316f2cf81de2e2404c7e22981372f476980dd64 Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Wed, 5 Aug 2026 17:42:35 -0600
Subject: [PATCH 16/17] Update .defenseHMMER()

---
 R/runHMMER.R | 17 +++++++++++++----
 1 file changed, 13 insertions(+), 4 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index e50cbcc..fb0b0a6 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -1215,6 +1215,9 @@ if (length(hmm_files) == 0) {
     "protein_DefenseCas.faa"
   )
 
+   # required to define the database size for hmmsearch --Z and --domZ parameters
+  Total_proteins <- nrow(prot_seqs)
+
   readr::write_lines(
     paste0(
       ">",
@@ -1253,7 +1256,7 @@ if (length(hmm_files) == 0) {
         ".tbl"
       )
     )
-
+    
     output <- system2(
       "docker",
       args = c(
@@ -1265,9 +1268,12 @@ if (length(hmm_files) == 0) {
         paste0(dirname(hmm_file), ":/db"),
         docker_image,
         "hmmsearch",
+         "--notextw",
         "--cpu",
         as.character(threads),
-        "--tblout",
+         "-Z", Total_proteins,
+        "--domZ", Total_proteins,
+        "--domtblout",
         file.path(
           "/work",
           basename(tbl_file)
@@ -1303,7 +1309,10 @@ if (length(hmm_files) == 0) {
       ) |>
       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
   }
@@ -1314,7 +1323,7 @@ if (length(hmm_files) == 0) {
 
   combined_tbl <- dplyr::bind_rows(
     all_hits
-  )
+  ) 
 
   parquet_file <- file.path(
     output_path,

From ce95c97d0a9a7e827657e6b8980ac82937649402 Mon Sep 17 00:00:00 2001
From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com>
Date: Wed, 5 Aug 2026 17:55:40 -0600
Subject: [PATCH 17/17] Update documentations

---
 R/runHMMER.R | 26 ++++++++++++++------------
 1 file changed, 14 insertions(+), 12 deletions(-)

diff --git a/R/runHMMER.R b/R/runHMMER.R
index fb0b0a6..81d5f8c 100644
--- a/R/runHMMER.R
+++ b/R/runHMMER.R
@@ -1,4 +1,4 @@
-#' Validate a HMM file has old HMMER3 format and remove them
+#' Validate if a HMM file has old HMMER3 format and remove them
 #'
 #' @param hmm_file
 #'
@@ -354,14 +354,13 @@ if (length(source_hmms) == 0) {
   )
 }
 
-#' Title
+#' Parsing HMM database to extract profile names, accessions and descriptions
 #'
-#' @param hmm_file
+#' @param hmm_file path to the HMM database file (`.hmm`)
 #'
-#' @returns
+#' @returns a tibble 
 #'
-#' @export
-#' @examples
+#' @keywords internal
 .parse_hmmer_profiles <- function(hmm_file) {
 
   lines <- readLines(hmm_file, warn = FALSE)
@@ -398,17 +397,21 @@ if (length(source_hmms) == 0) {
   )
 }
 
-#' Title
+#' The function to run HMMER with docker
 #'
-#' @param JOB_NAME
+#' @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
 #'
-#' @export
-#' @examples
+#' @keywords internal
 .runHmmerJob <- function(JOB_NAME, FASTA, DB, Total_proteins, 
   output_path = NULL, db_paths, 
   docker_image = "staphb/hmmer", threads = 8L, 
@@ -661,10 +664,9 @@ invisible(final_parquets)
   #   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 `--tblout` file and returns a tidy tibble with one row per
+#' 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.
 #'