diff --git a/DESCRIPTION b/DESCRIPTION
index b3ff708..64266d6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -53,5 +53,6 @@ biocViews: GenomeAssembly, Annotation, Sequencing
 VignetteBuilder: knitr
 Suggests:
     rmarkdown,
+    writexl,
     testthat (>= 3.0.0)
 Config/roxygen2/version: 8.0.0
diff --git a/NAMESPACE b/NAMESPACE
index 69af5ae..70f9411 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -6,6 +6,8 @@ export(CDHIT2duckdb)
 export(buildClusterFeatureMap)
 export(cleanData)
 export(cleanMetaData)
+export(exportProcessedData)
+export(exportTables)
 export(generatePlots)
 export(generateSummary)
 export(prepareGenomes)
diff --git a/R/data_curation.R b/R/data_curation.R
index 7edc518..6c64a05 100644
--- a/R/data_curation.R
+++ b/R/data_curation.R
@@ -288,8 +288,8 @@
 #' Deviation-based filters are optional and are best suited to single-taxon runs.
 #'
 #' @param genome_tbl A tibble of genome metadata containing BV-BRC genome columns.
-#' @param checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%).
-#' @param checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%).
+#' @param max_checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%).
+#' @param min_checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%).
 #' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content.
 #' @param length_deviations Optional numeric scalar. Maximum SDs from the median genome length.
 #' @param cds_deviations Optional numeric scalar. Maximum SDs from the median CDS count.
@@ -302,8 +302,8 @@
 #' }
 #' @keywords internal
 .apply_metadata_qc <- function(genome_tbl,
-                               checkm_contam = 5,
-                               checkm_complete = 95,
+                               max_checkm_contam = 5,
+                               min_checkm_complete = 95,
                                gc_deviations = NULL,
                                length_deviations = NULL,
                                cds_deviations = NULL) {
@@ -362,8 +362,8 @@
   qc <- qc |>
     dplyr::mutate(
       flag_checkm_missing = is.na(.data$checkm_completeness_qc) | is.na(.data$checkm_contamination_qc),
-      flag_checkm_contam = !is.na(.data$checkm_contamination_qc) & .data$checkm_contamination_qc > checkm_contam,
-      flag_checkm_complete = !is.na(.data$checkm_completeness_qc) & .data$checkm_completeness_qc < checkm_complete,
+      flag_checkm_contam = !is.na(.data$checkm_contamination_qc) & .data$checkm_contamination_qc > max_checkm_contam,
+      flag_checkm_complete = !is.na(.data$checkm_completeness_qc) & .data$checkm_completeness_qc < min_checkm_complete,
       flag_checkm = .data$flag_checkm_missing | .data$flag_checkm_contam | .data$flag_checkm_complete,
       flag_length = if (!is.null(length_deviations) && is.finite(sd_len) && sd_len > 0) {
         !is.na(z_len) & z_len > length_deviations
@@ -383,7 +383,7 @@
       qc_keep = !(.data$flag_checkm | .data$flag_length | .data$flag_gc | .data$flag_cds)
     )
 
-  make_rej <- function(flag, rule, observed, threshold, comparator) {
+  make_rejects <- function(flag, rule, observed, threshold, comparator) {
     idx <- which(flag %in% TRUE)
     if (!length(idx)) {
       return(NULL)
@@ -400,32 +400,32 @@
     )
   }
 
-  rej_list <- list(
-    make_rej(
+  reject_list <- list(
+    make_rejects(
       qc$flag_checkm_missing,
       "checkm_missing",
       ifelse(qc$flag_checkm_missing, "missing", NA_character_),
       "present",
       "is present"
     ),
-    make_rej(
+    make_rejects(
       qc$flag_checkm_contam,
       "checkm_contamination",
       qc$checkm_contamination_qc,
-      checkm_contam,
+      max_checkm_contam,
       ">"
     ),
-    make_rej(
+    make_rejects(
       qc$flag_checkm_complete,
       "checkm_completeness",
       qc$checkm_completeness_qc,
-      checkm_complete,
+      min_checkm_complete,
       "<"
     )
   )
 
   if (!is.null(gc_deviations)) {
-    rej_list[[length(rej_list) + 1L]] <- make_rej(
+    reject_list[[length(reject_list) + 1L]] <- make_rejects(
       qc$flag_gc,
       "gc_deviation_sd",
       z_gc,
@@ -435,7 +435,7 @@
   }
 
   if (!is.null(length_deviations)) {
-    rej_list[[length(rej_list) + 1L]] <- make_rej(
+    reject_list[[length(reject_list) + 1L]] <- make_rejects(
       qc$flag_length,
       "length_deviation_sd",
       z_len,
@@ -445,7 +445,7 @@
   }
 
   if (!is.null(cds_deviations)) {
-    rej_list[[length(rej_list) + 1L]] <- make_rej(
+    reject_list[[length(reject_list) + 1L]] <- make_rejects(
       qc$flag_cds,
       "cds_deviation_sd",
       z_cds,
@@ -454,7 +454,7 @@
     )
   }
 
-  rejections <- dplyr::bind_rows(rej_list)
+  rejections <- dplyr::bind_rows(reject_list)
   if (nrow(rejections) == 0L) {
     rejections <- tibble::tibble(
       genome_id = character(),
@@ -516,35 +516,28 @@
   invisible(n_removed)
 }
 
-# Make sure the BV-BRC metadata live where they're supposed to
+# Make sure the BV-BRC metadata live where they're supposed to, and are fresh
 .ensure_bvbrc_cache <- function(base_dir = ".",
                                 verbose = TRUE,
+                                max_age_days = 30L,
                                 cache_rel = file.path("data", "bvbrc", "bvbrcData.duckdb"),
                                 cache_table = "bvbrc_bac_data") {
   base_dir <- normalizePath(base_dir, mustWork = FALSE)
   cache_db <- file.path(base_dir, cache_rel)
 
-  need_build <- !file.exists(cache_db)
-  con_cache <- NULL
+  # Always delegate to .updateBVBRCdata() so its max_age_days staleness check
+  # actually runs. Previously this only rebuilt when the cache file/table was
+  # entirely absent, so an existing-but-stale cache was silently reused
+  # forever, and results (e.g. genome counts) drifted across machines/sessions
+  # depending on whenever each cache happened to be built.
+  .updateBVBRCdata(base_dir = base_dir, max_age_days = max_age_days, verbose = verbose)
 
-  if (!need_build) {
-    con_cache <- DBI::dbConnect(duckdb::duckdb(), dbdir = cache_db)
-    on.exit(try(DBI::dbDisconnect(con_cache, shutdown = TRUE), silent = TRUE), add = TRUE)
-    need_build <- !(cache_table %in% DBI::dbListTables(con_cache))
-  }
-
-  if (need_build) {
-    if (isTRUE(verbose)) message("BV-BRC cache missing or incomplete. Building via .updateBVBRCdata(). Please wait.")
-    .updateBVBRCdata(base_dir = base_dir, verbose = verbose)
+  if (!file.exists(cache_db)) stop("After .updateBVBRCdata(), cache DB still missing at: ", cache_db)
 
-    if (!is.null(con_cache)) try(DBI::dbDisconnect(con_cache, shutdown = TRUE), silent = TRUE)
-    if (!file.exists(cache_db)) stop("After .updateBVBRCdata(), cache DB still missing at: ", cache_db)
-
-    con_cache <- DBI::dbConnect(duckdb::duckdb(), dbdir = cache_db)
-    on.exit(try(DBI::dbDisconnect(con_cache, shutdown = TRUE), silent = TRUE), add = TRUE)
-    if (!(cache_table %in% DBI::dbListTables(con_cache))) {
-      stop("After .updateBVBRCdata(), table '", cache_table, "' still not found in ", cache_db)
-    }
+  con_cache <- DBI::dbConnect(duckdb::duckdb(), dbdir = cache_db, read_only = TRUE)
+  on.exit(try(DBI::dbDisconnect(con_cache, shutdown = TRUE), silent = TRUE), add = TRUE)
+  if (!(cache_table %in% DBI::dbListTables(con_cache))) {
+    stop("After .updateBVBRCdata(), table '", cache_table, "' still not found in ", cache_db)
   }
 
   invisible(cache_db)
@@ -1191,8 +1184,8 @@
 #' @param abx Character or vector. Antibiotic filter. "All" for all antibiotics, else names.
 #' @param overwrite Logical. If FALSE and DuckDB exists already, abort. Default FALSE.
 #' @param image Character. Docker image. Default "danylmb/bvbrc:5.3".
-#' @param checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%).
-#' @param checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%).
+#' @param max_checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%).
+#' @param min_checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%).
 #' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content.
 #' @param length_deviations Optional numeric scalar. Maximum SDs from the median genome length.
 #' @param cds_deviations Optional numeric scalar. Maximum SDs from the median CDS count.
@@ -1210,12 +1203,14 @@ retrieveMetadata <- function(user_bacs,
                              abx = "All",
                              overwrite = FALSE,
                              image = "danylmb/bvbrc:5.3",
-                             checkm_contam = 5,
-                             checkm_complete = 95,
+                             max_checkm_contam = 5,
+                             min_checkm_complete = 95,
                              gc_deviations = NULL,
                              length_deviations = NULL,
                              cds_deviations = NULL,
                              debug = FALSE,
+                             export_tables = FALSE,
+                             load_tables = FALSE,
                              verbose = TRUE) {
   base_dir <- normalizePath(base_dir, mustWork = FALSE)
 
@@ -1390,8 +1385,8 @@ retrieveMetadata <- function(user_bacs,
 
   qc_out <- .apply_metadata_qc(
     genome_tbl = combined_genome_data_tbl,
-    checkm_contam = checkm_contam,
-    checkm_complete = checkm_complete,
+    max_checkm_contam = max_checkm_contam,
+    min_checkm_complete = min_checkm_complete,
     gc_deviations = gc_deviations,
     length_deviations = length_deviations,
     cds_deviations = cds_deviations
@@ -1519,6 +1514,24 @@ retrieveMetadata <- function(user_bacs,
     }
   }
 
+  export_res <- NULL
+  if (isTRUE(export_tables) || isTRUE(load_tables)) {
+    export_res <- exportTables(
+      paths$db_path,
+      export_tables = export_tables,
+      load_tables = load_tables,
+      verbose = verbose
+    )
+  }
+
+  if (isTRUE(load_tables)) {
+    return(list(
+      duckdbConnection = con,
+      table_name = "metadata",
+      data = if (!is.null(export_res)) export_res$data else NULL
+    ))
+  }
+
   list(duckdbConnection = con, table_name = "metadata")
 }
 
@@ -1943,7 +1956,6 @@ retrieveGenomes <- function(base_dir = ".",
                             cli_gff_workers = 8L,
                             chunk_size = 50L,
                             evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"),
-                            # NEW
                             verbose = TRUE) {
   method <- match.arg(method)
   evidence_mode <- match.arg(evidence_mode)
@@ -1959,7 +1971,7 @@ retrieveGenomes <- function(base_dir = ".",
 
   if (has_filtered) {
     if (isTRUE(verbose))
-      message("Using existing 'filtered' table (skipping re-filter).")
+      message("Found existing 'filtered' table.")
     con <- con0
     tbl <- "filtered"
     on.exit(try(DBI::dbDisconnect(con0, shutdown = TRUE), silent = TRUE), add = TRUE)
@@ -1979,6 +1991,7 @@ retrieveGenomes <- function(base_dir = ".",
     on.exit(try(DBI::dbDisconnect(con, shutdown = TRUE), silent = TRUE), add = TRUE)
   }
 
+  # What genomes need to be downloaded? Build set as `ids`
   ids <- tibble::as_tibble(DBI::dbReadTable(con, tbl)) |>
     dplyr::distinct(`genome.genome_id`) |>
     dplyr::pull(`genome.genome_id`)
@@ -1989,6 +2002,7 @@ retrieveGenomes <- function(base_dir = ".",
   dir.create(genome_path, recursive = TRUE, showWarnings = FALSE)
   dir.create(logs_dir, recursive = TRUE, showWarnings = FALSE)
 
+  # Checks what is already  downloaded vs. the full list needed; takes diff
   if (isTRUE(skip_existing)) {
     already <- .list_complete(genome_path, ids)
     if (isTRUE(verbose))
@@ -1996,6 +2010,7 @@ retrieveGenomes <- function(base_dir = ".",
     ids <- setdiff(ids, already)
   }
 
+  # Is diff length 0? If so, all genomes ready to go!
   if (length(ids) == 0L) {
     if (isTRUE(verbose))
       message("All genomes already complete.")
@@ -2207,8 +2222,8 @@ genomeList <- function(base_dir = ".",
 #'    return very large download lists for many species!
 #' @param num_workers Integer. Parallel workers used for genome download.
 #'    Applied to both FTP and CLI download branches. Default: 8.
-#' @param checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%).
-#' @param checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%).
+#' @param max_checkm_contam Numeric scalar. Maximum allowed CheckM contamination (%).
+#' @param min_checkm_complete Numeric scalar. Minimum allowed CheckM completeness (%).
 #' @param gc_deviations Optional numeric scalar. Maximum SDs from the median GC content.
 #' @param length_deviations Optional numeric scalar. Maximum SDs from the median genome length.
 #' @param cds_deviations Optional numeric scalar. Maximum SDs from the median CDS count.
@@ -2227,11 +2242,13 @@ prepareGenomes <- function(user_bacs,
                            overwrite = FALSE,
                            num_workers = 8L,
                            evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"),
-                           checkm_contam = 5,
-                           checkm_complete = 95,
+                           max_checkm_contam = 5,
+                           min_checkm_complete = 95,
                            gc_deviations = NULL,
                            length_deviations = NULL,
                            cds_deviations = NULL,
+                           export_tables = FALSE,
+                           load_tables = FALSE,
                            debug = FALSE,
                            verbose = TRUE) {
   method <- match.arg(method)
@@ -2248,8 +2265,8 @@ prepareGenomes <- function(user_bacs,
     base_dir = base_dir,
     abx = "All",
     overwrite = overwrite,
-    checkm_contam = checkm_contam,
-    checkm_complete = checkm_complete,
+    max_checkm_contam = max_checkm_contam,
+    min_checkm_complete = min_checkm_complete,
     gc_deviations = gc_deviations,
     length_deviations = length_deviations,
     cds_deviations = cds_deviations,
@@ -2315,5 +2332,185 @@ prepareGenomes <- function(user_bacs,
     message("Continue with downstream processing using:")
     message('runDataProcessing("', normalizePath(paths$db_path), '")')
   }
+
+  export_res <- NULL
+  if (isTRUE(export_tables) || isTRUE(load_tables)) {
+    export_res <- exportTables(
+      paths$db_path,
+      export_tables = export_tables,
+      load_tables = load_tables,
+      verbose = verbose
+    )
+  }
+
+  if (isTRUE(load_tables)) {
+    return(list(
+      duckdb_path = paths$db_path,
+      table_name = "files",
+      data = if (!is.null(export_res)) export_res$data else NULL
+    ))
+  }
+
   invisible(out)
 }
+
+#' Export DuckDB tables and optionally load them into R
+#'
+#' Writes selected DuckDB tables to CSV files and optionally returns them as
+#' in-memory R data frames.
+#'
+#' @param duckdb_path Character. Path to the DuckDB file.
+#' @param output_dir Character or NULL. Directory for exports. Defaults to
+#'   file.path(dirname(duckdb_path), "exports").
+#' @param tables Character vector or NULL. Tables to export. If NULL, exports
+#'   all tables in the database.
+#' @param skip_tables Character vector of table names to exclude.
+#' @param include_summary Logical. If TRUE, writes summary.csv and summary.txt.
+#' @param export_tables Logical. If TRUE, writes tables to CSV files in `output_dir`.
+#' @param load_tables Logical. If TRUE, also return the exported tables as
+#'   in-memory R data frames.
+#' @param verbose Logical. If TRUE, prints progress messages.
+#'
+#' @return Invisibly returns a list with exported file paths and, if requested,
+#'   in-memory tables.
+#' @export
+exportTables <- function(duckdb_path,
+                         output_dir = NULL,
+                         tables = NULL,
+                         skip_tables = NULL,
+                         include_summary = TRUE,
+                         export_tables = TRUE,
+                         load_tables = FALSE,
+                         verbose = TRUE) {
+  duckdb_path <- normalizePath(duckdb_path, mustWork = TRUE)
+
+  if (is.null(output_dir)) {
+    output_dir <- file.path(dirname(duckdb_path), "exports")
+  }
+  output_dir <- normalizePath(output_dir, mustWork = FALSE)
+  dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
+
+  con <- DBI::dbConnect(duckdb::duckdb(), dbdir = duckdb_path)
+  on.exit(try(DBI::dbDisconnect(con, shutdown = TRUE), silent = TRUE), add = TRUE)
+
+  available_tables <- DBI::dbListTables(con)
+  if (!length(available_tables)) {
+    stop("No tables found in DuckDB: ", duckdb_path)
+  }
+
+  # More tables will exist in a full DuckDB output after data processing, but
+  # limiting this to the basic genome stats and metadata you'd get from only
+  # running data_curation.R
+  basic_tables <- c(
+    "bac_data",
+    "filtered",
+    "metadata",
+    "genome_data",
+    "amr_phenotype",
+    "metadata_qc",
+    "metadata_qc_rejections"
+  )
+
+  if (is.null(tables)) {
+    tables <- intersect(basic_tables, available_tables)
+  } else {
+    tables <- intersect(as.character(tables), available_tables)
+    if (!length(tables)) {
+      stop("None of the requested tables were found in the DuckDB.")
+    }
+  }
+
+  if (length(skip_tables)) {
+    tables <- setdiff(tables, skip_tables)
+  }
+  if (!length(tables)) {
+    stop("No tables left to export after applying skip_tables.")
+  }
+
+  exported_files <- list()
+  loaded_tables <- list()
+
+  for (tbl in tables) {
+    out_file <- file.path(output_dir, paste0(tbl, ".csv"))
+    df <- DBI::dbReadTable(con, tbl)
+
+    if (isTRUE(export_tables)) {
+      readr::write_csv(df, out_file, na = "")
+      exported_files[[tbl]] <- out_file
+      if (isTRUE(verbose)) {
+        message("Exported table: ", tbl, " -> ", out_file)
+      }
+    }
+
+    if (isTRUE(load_tables)) {
+      loaded_tables[[tbl]] <- df
+    }
+  }
+
+  count_if_present <- function(tbl) {
+    if (tbl %in% available_tables) {
+      as.character(DBI::dbGetQuery(
+        con,
+        paste0("SELECT COUNT(*) AS n FROM ", DBI::dbQuoteIdentifier(con, tbl))
+      )$n[[1]])
+    } else {
+      NA_character_
+    }
+  }
+
+  summary_tbl <- tibble::tibble(
+    metric = c(
+      "duckdb_path",
+      "export_dir",
+      "tables_exported",
+      "table_names",
+      "metadata_rows",
+      "genome_data_rows",
+      "amr_phenotype_rows",
+      "metadata_qc_rows",
+      "metadata_qc_rejections_rows",
+      "filtered_rows",
+      "files_rows",
+      "bac_data_rows"
+    ),
+    value = c(
+      duckdb_path,
+      output_dir,
+      as.character(length(tables)),
+      paste(tables, collapse = ", "),
+      count_if_present("metadata"),
+      count_if_present("genome_data"),
+      count_if_present("amr_phenotype"),
+      count_if_present("metadata_qc"),
+      count_if_present("metadata_qc_rejections"),
+      count_if_present("filtered"),
+      count_if_present("files"),
+      count_if_present("bac_data")
+    )
+  )
+
+  if (isTRUE(export_tables) && isTRUE(include_summary)) {
+    readr::write_csv(summary_tbl, file.path(output_dir, "summary.csv"), na = "")
+    writeLines(
+      c(
+        paste0("DuckDB: ", duckdb_path),
+        paste0("Export directory: ", output_dir),
+        paste0("Tables exported: ", length(tables)),
+        paste0("Table names: ", paste(tables, collapse = ", "))
+      ),
+      file.path(output_dir, "summary.txt"),
+      useBytes = TRUE
+    )
+    if (isTRUE(verbose)) {
+      message("Exported summary files.")
+    }
+  }
+
+  invisible(list(
+    export_dir = output_dir,
+    tables = tables,
+    files = exported_files,
+    data = if (isTRUE(load_tables)) loaded_tables else NULL,
+    summary = summary_tbl
+  ))
+}
diff --git a/R/data_processing.R b/R/data_processing.R
index 945f59d..cbb43a0 100644
--- a/R/data_processing.R
+++ b/R/data_processing.R
@@ -40,6 +40,15 @@ NULL
 #' @param cluster_threshold Numeric. Sequence identity threshold (`--threshold`). Default `0.95`.
 #' @param family_seq_identity Numeric. Gene family clustering identity (`-f`). Default `0.5`.
 #' @param panaroo_threads_per_job Integer. Number of threads for Panaroo and parallel execution.
+#' @param refind_mode Character. Panaroo's `--refind-mode` (`"off"`, `"default"`, or
+#'   `"strict"`). Refinding searches for and recovers gene calls that annotation
+#'   tools missed, comparing each candidate against the rest of the pangenome.
+#'   Caveat: this search can take substantially longer (and in rare cases fail to
+#'   complete within hours) when a genome carries a cluster of CDS with internal
+#'   stop codons, which existing upstream genome-quality fields do not flag.
+#'   Default `"off"` for now, to avoid that runtime risk; plan to move this back to
+#'   `"default"` once a QC step upstream (e.g. in `.apply_metadata_qc()`) can screen
+#'   out affected genomes before they reach Panaroo.
 #'
 #' @returns A list of results for each Panaroo batch in its output directory.
 #'
@@ -51,9 +60,11 @@ NULL
                             len_dif_percent,
                             cluster_threshold,
                             family_seq_identity,
-                            panaroo_threads_per_job) {
-  output_path <- .docker_path(output_path)
+                            panaroo_threads_per_job,
+                            refind_mode = c("off", "default", "strict")) {
+  refind_mode <- match.arg(refind_mode)
   dir.create(output_path, recursive = TRUE, showWarnings = FALSE)
+  output_path <- .docker_path(output_path)
 
   # Fail fast if Docker is missing
   if (!nzchar(Sys.which("docker"))) {
@@ -93,13 +104,14 @@ NULL
     "--rm",
     "-v", paste0(mount_host, ":", mount_cont),
     "-w", mount_cont,
-    "staphb/panaroo:1.5.1",
+    "staphb/panaroo:1.7.0",
     "panaroo",
     "-i", genome_filepath_cont,
     "-o", output_dir_cont,
     "--clean-mode", "strict",
     "--merge_paralogs",
     "--remove-invalid-genes",
+    "--refind-mode", refind_mode,
     "--core_threshold", as.character(core_threshold),
     "--len_dif_percent", as.character(len_dif_percent),
     "--threshold", as.character(cluster_threshold),
@@ -126,6 +138,123 @@ NULL
   invisible(res)
 }
 
+#' Remove pseudogene annotations from Panaroo input GFF files
+#'
+#' Cleans GFF annotation files of `pseudogene` feature records only. Cleaned
+#' GFFs are written to a subdirectory under `output_path` and swapped into the
+#' Panaroo input list, leaving the original genome annotations alone.
+#'
+#' This optional preprocessing step can reduce weird runtime stalls during
+#' Panaroo graph construction for some BV-BRC/PATRIC genome annotations that
+#' contain troublesome pseudogene features.
+#'
+#' @param panaroo_input_files Character vector of `"gff fna"` input lines used
+#'   by Panaroo.
+#' @param output_path Character scalar. Base directory for temporary cleaned
+#'   GFF files and audit outputs.
+#' @param clean_dir Character scalar. Name of the subdirectory created beneath
+#'   `output_path` to store cleaned GFF files. Default `"gff_clean"`.
+#'
+#' @return A list containing:
+#' \itemize{
+#'   \item `panaroo_input_files` — rewritten Panaroo input lines pointing to the
+#'   cleaned GFF files.
+#'   \item `audit` — a tibble summarizing, for each genome, the total number of
+#'   annotated features, the number of pseudogenes removed, and the number of
+#'   remaining features.
+#' }
+#'
+#' @details
+#' This performs lightweight preprocessing only, removing feature records whose
+#' third GFF column is exactly `"pseudogene"` and does not otherwise modify
+#' annotation coordinates, attributes, or sequence files. FASTA paths are unchanged.
+#'
+#' @keywords internal
+.stripPseudogeneGFFs <- function(panaroo_input_files,
+                                 output_path,
+                                 clean_dir = "gff_clean") {
+  # Normalize our paths
+  panaroo_input_files <- as.character(panaroo_input_files)
+  output_path <- .docker_path(output_path)
+
+  # Set the directory to place cleaned GFFs into
+  clean_root <- file.path(output_path, clean_dir)
+  dir.create(clean_root, recursive = TRUE, showWarnings = FALSE)
+
+  # Where the cleaned up Panaroo input and clean audit is stored
+  out_lines <- character(length(panaroo_input_files))
+  audit <- vector("list", length(panaroo_input_files))
+
+  for (i in seq_along(panaroo_input_files)) {
+
+    # Read the Panaroo gff + fna input lines
+    line <- panaroo_input_files[[i]]
+    parts <- strsplit(line, "\\s+")[[1]]
+
+    # If you're missing either a gff or an fna file in there somehow
+    if (length(parts) < 2L) {
+      stop("Broken Panaroo input line: ", line)
+    }
+
+    # Read in the files parsed above
+    gff_in <- .docker_path(parts[1])
+    fna_in <- .docker_path(parts[2])
+
+    # If it didn't read in
+    if (!file.exists(gff_in)) {
+      stop("Missing GFF file: ", gff_in)
+    }
+    if (!file.exists(fna_in)) {
+      stop("Missing FNA file: ", fna_in)
+    }
+
+    # What we're saving out
+    gff_out <- file.path(clean_root, basename(gff_in))
+
+    # Read in the GFF lines and fine the comment lined headers
+    gff_lines <- readLines(gff_in, warn = FALSE)
+    is_header <- startsWith(gff_lines, "#")
+    body <- gff_lines[!is_header]
+
+    # If there's nothing in there to parse
+    if (length(body) == 0L) {
+      writeLines(gff_lines, gff_out, useBytes = TRUE)
+      n_total <- 0L
+      n_pseudogene <- 0L
+      n_kept <- 0L
+    } else {
+      # Otherwise, find the pseudogene lines and save everything but those
+      fields <- strsplit(body, "\t", fixed = TRUE)
+      types <- vapply(fields, function(x) if (length(x) >= 3L) x[[3]] else NA_character_, character(1))
+      keep <- !is.na(types) & types != "pseudogene"
+
+      cleaned <- c(gff_lines[is_header], body[keep])
+      writeLines(cleaned, gff_out, useBytes = TRUE)
+
+      # We love stats
+      n_total <- length(body)
+      n_pseudogene <- sum(!keep, na.rm = TRUE)
+      n_kept <- sum(keep, na.rm = TRUE)
+    }
+
+    # Record what we did and save it into the audit log
+    audit[[i]] <- tibble::tibble(
+      gff_in = gff_in,
+      gff_out = gff_out,
+      n_total_features = n_total,
+      n_pseudogene = n_pseudogene,
+      n_kept = n_kept
+    )
+
+    out_lines[[i]] <- paste(gff_out, fna_in)
+  }
+
+  list(
+    panaroo_input_files = out_lines,
+    audit = dplyr::bind_rows(audit)
+  )
+}
+
 
 #' Run Panaroo for Pangenome Analysis in Parallel Batches
 #'
@@ -143,6 +272,9 @@ NULL
 #' @param threads Integer. Number of threads for Panaroo and parallel execution. Default `8`.
 #' @param split_jobs Logical. If TRUE, split into multiple smaller pangenome
 #'   generation jobs that can be merged by [.mergePanaroo()]. If FALSE, all isolates in one run.
+#' @param refind_mode Character. Panaroo's `--refind-mode` (`"off"`, `"default"`, or
+#'   `"strict"`). See [.processPanaroo()] for what refinding does and the runtime
+#'   caveat behind the current default. Default `"off"`.
 #'
 #' @return A list of results for each Panaroo batch in its output directory.
 #'
@@ -159,7 +291,12 @@ NULL
                         cluster_threshold = 0.95,
                         family_seq_identity = 0.5,
                         threads = 8,
-                        split_jobs = FALSE) {
+                        split_jobs = FALSE,
+                        refind_mode = c("off", "default", "strict"),
+                        strip_pseudogenes = FALSE,
+                        pseudogene_clean_dir = "gff_clean",
+                        write_pseudogene_audit = TRUE) {
+  refind_mode <- match.arg(refind_mode)
   duckdb_path <- normalizePath(duckdb_path)
   con <- DBI::dbConnect(duckdb::duckdb(), duckdb_path)
   on.exit(try(DBI::dbDisconnect(con, shutdown = FALSE), silent = TRUE), add = TRUE)
@@ -167,6 +304,7 @@ NULL
   if (missing(output_path) || output_path %in% c(".", "results", "results/")) {
     output_path <- dirname(duckdb_path)
   }
+  dir.create(output_path, recursive = TRUE, showWarnings = FALSE)
   output_path <- normalizePath(output_path)
 
   genome_query_output <- DBI::dbGetQuery(con, "SELECT * FROM files ORDER BY genome_id")
@@ -177,6 +315,19 @@ NULL
   # Drop true NAs
   panaroo_input_files <- panaroo_input_files[!is.na(panaroo_input_files)]
 
+  # Cleaning pseudogene lines out of GFF files
+  if (isTRUE(strip_pseudogenes)) {
+    cleaned <- .stripPseudogeneGFFs(
+      panaroo_input_files = panaroo_input_files,
+      output_path = output_path,
+      clean_dir = pseudogene_clean_dir
+    )
+    panaroo_input_files <- cleaned$panaroo_input_files
+    if (isTRUE(write_pseudogene_audit)) {
+      readr::write_csv(cleaned$audit, file.path(output_path, "panaroo_pseudogene_audit.csv"))
+    }
+  }
+
   split_files <- strsplit(panaroo_input_files, " ")
 
   valid_entries <- purrr::map_lgl(split_files, function(paths) {
@@ -220,7 +371,8 @@ NULL
       len_dif_percent         = len_dif_percent,
       cluster_threshold       = cluster_threshold,
       family_seq_identity     = family_seq_identity,
-      panaroo_threads_per_job = panaroo_threads_per_job
+      panaroo_threads_per_job = panaroo_threads_per_job,
+      refind_mode             = refind_mode
     ),
     .options = furrr::furrr_options(seed = TRUE)
   )
@@ -278,7 +430,7 @@ NULL
       "--rm",
       "-v", paste0(mount_host, ":", mount_cont),
       "-w", mount_cont,
-      "staphb/panaroo:1.5.1",
+      "staphb/panaroo:1.7.0",
       "panaroo-merge",
       "-d", dir_args,
       "-o", file.path(mount_cont, "merge_output"),
@@ -480,8 +632,8 @@ NULL
   if (missing(output_path) || output_path %in% c(".", "results", "results/")) {
     output_path <- dirname(duckdb_path)
   }
+  dir.create(output_path, recursive = TRUE, showWarnings = FALSE)
   output_path <- .docker_path(output_path)
-  if (!dir.exists(output_path)) dir.create(output_path, recursive = TRUE)
 
   con <- DBI::dbConnect(duckdb::duckdb(), duckdb_path)
   on.exit(try(DBI::dbDisconnect(con, shutdown = FALSE), silent = TRUE), add = TRUE)
@@ -588,13 +740,17 @@ NULL
 #'
 #' @param threads Integer. Total CPU budget to allocate for Panaroo.
 #'   If `split_jobs = TRUE`, threads are divided across batches.
-#'   Default: `16`.
+#'   Default: `8`.
 #'
 #' @param split_jobs Logical. If `TRUE`, Panaroo is run in multiple parallel
 #'   batches (up to 5, depending on dataset size), and batch outputs are merged
 #'   using `.mergePanaroo()`. If `FALSE`, only one Panaroo invocation is run.
 #'   Default: `FALSE`.
 #'
+#' @param refind_mode Character. Panaroo's `--refind-mode` (`"off"`, `"default"`, or
+#'   `"strict"`). See [.processPanaroo()] for what refinding does and the runtime
+#'   caveat behind the current default. Default `"off"`.
+#'
 #' @param verbose Logical. Print status messages during Panaroo execution,
 #'   merging, and DuckDB import. Default: `TRUE`.
 #'
@@ -635,7 +791,7 @@ NULL
 #' runPanaroo2Duckdb(
 #'   duckdb_path = "data/Shigella_flexneri/Sfl.duckdb",
 #'   output_path = "data/Shigella_flexneri",
-#'   threads     = 12,
+#'   threads     = 8,
 #'   split_jobs  = FALSE
 #' )
 #'
@@ -655,9 +811,14 @@ runPanaroo2Duckdb <- function(duckdb_path,
                               len_dif_percent = 0.95,
                               cluster_threshold = 0.95,
                               family_seq_identity = 0.5,
-                              threads = 16,
+                              threads = 8,
                               split_jobs = FALSE,
+                              refind_mode = c("off", "default", "strict"),
+                              strip_pseudogenes = FALSE,
+                              pseudogene_clean_dir = "gff_clean",
+                              write_pseudogene_audit = TRUE,
                               verbose = TRUE) {
+  refind_mode <- match.arg(refind_mode)
   duckdb_path <- normalizePath(duckdb_path)
   out_dir <- if (is.null(output_path)) dirname(duckdb_path) else normalizePath(output_path)
 
@@ -670,7 +831,11 @@ runPanaroo2Duckdb <- function(duckdb_path,
     cluster_threshold = cluster_threshold,
     family_seq_identity = family_seq_identity,
     threads = threads,
-    split_jobs = split_jobs
+    split_jobs = split_jobs,
+    refind_mode = refind_mode,
+    strip_pseudogenes = strip_pseudogenes,
+    pseudogene_clean_dir = pseudogene_clean_dir,
+    write_pseudogene_audit = write_pseudogene_audit
   )
 
   # Identify Panaroo outputs that contain a final_graph.gml file
@@ -883,6 +1048,7 @@ CDHIT2duckdb <- function(duckdb_path,
   if (missing(output_path) || output_path %in% c(".", "results", "results/")) {
     output_path <- dirname(duckdb_path) # e.g., ./results/<bug>
   }
+  dir.create(output_path, recursive = TRUE, showWarnings = FALSE)
   output_path <- normalizePath(output_path)
 
   cdhit_outputs <- .runCDHIT(duckdb_path,
@@ -1084,11 +1250,10 @@ CDHIT2duckdb <- function(duckdb_path,
                            file_format,
                            docker_image = sprintf("interpro/interproscan:%s", "5.76-107.0")) {
   # Normalize and mount paths
+  dir.create(file.path(path, "tmp", "iprscan"), recursive = TRUE, showWarnings = FALSE)
   path <- .docker_path(path)
   bind_data <- .docker_path(ipr_data_path)
 
-  dir.create(file.path(path, "tmp", "iprscan"), recursive = TRUE, showWarnings = FALSE)
-
   fasta_sequences <- Biostrings::AAStringSet(chunk$sequence)
   names(fasta_sequences) <- chunk$name
   temp_fasta_file <- tempfile(tmpdir = path, fileext = ".fa")
@@ -1174,6 +1339,7 @@ domainFromIPR <- function(duckdb_path,
   if (missing(path) || path %in% c(".", "results", "results/")) {
     path <- dirname(duckdb_path)
   }
+  dir.create(path, recursive = TRUE, showWarnings = FALSE)
   path <- normalizePath(path)
 
   ipr_image <- sprintf("%s:%s", docker_repo, ipr_version)
@@ -1637,7 +1803,7 @@ cleanData <- function(duckdb_path, path) {
 #'   outputs and final Parquet files. If `NULL`, defaults to `dirname(duckdb_path)`.
 #'
 #' @param threads Integer. Shared concurrency budget used across tools (Panaroo, CD-HIT,
-#'   InterProScan). Passed through to each stage as appropriate. Defaults to `16`.
+#'   InterProScan). Passed through to each stage as appropriate. Defaults to `8`.
 #'
 #' @param panaroo_split_jobs Logical. If `TRUE`, Panaroo runs in multiple batches that can be
 #'   merged by [.mergePanaroo()]. If `FALSE`, Panaroo runs once on all isolates. Default: `FALSE`.
@@ -1645,6 +1811,9 @@ cleanData <- function(duckdb_path, path) {
 #' @param panaroo_len_dif_percent Numeric. Panaroo `--len_dif_percent`. Default: `0.95`.
 #' @param panaroo_cluster_threshold Numeric. Panaroo `--threshold`. Default: `0.95`.
 #' @param panaroo_family_seq_identity Numeric. Panaroo `-f` (gene family identity). Default: `0.5`.
+#' @param panaroo_refind_mode Character. Panaroo's `--refind-mode` (`"off"`, `"default"`,
+#'   or `"strict"`). See [.processPanaroo()] for what refinding does and the runtime
+#'   caveat behind the current default. Default `"off"`.
 #'
 #' @param cdhit_identity Numeric. CD-HIT `-c` identity threshold. Default: `0.9`.
 #' @param cdhit_word_length Integer. CD-HIT `-n` word length. Default: `5`.
@@ -1711,7 +1880,7 @@ cleanData <- function(duckdb_path, path) {
 #' runDataProcessing(
 #'   duckdb_path   = "data/Shigella_flexneri/Sfl.duckdb",
 #'   output_path   = "data/Shigella_flexneri",
-#'   threads       = 16,
+#'   threads       = 8,
 #'   ref_file_path = "data_raw/"
 #' )
 #'
@@ -1724,13 +1893,17 @@ cleanData <- function(duckdb_path, path) {
 runDataProcessing <- function(duckdb_path,
                               output_path = NULL,
                               # unified threads for all tools
-                              threads = 16,
+                              threads = 8,
                               # Panaroo
                               panaroo_split_jobs = FALSE,
                               panaroo_core_threshold = 0.90,
                               panaroo_len_dif_percent = 0.95,
                               panaroo_cluster_threshold = 0.95,
                               panaroo_family_seq_identity = 0.5,
+                              panaroo_refind_mode = c("off", "default", "strict"),
+                              panaroo_strip_pseudogenes = FALSE,
+                              panaroo_pseudogene_clean_dir = "gff_clean",
+                              panaroo_write_pseudogene_audit = TRUE,
                               # CD-HIT
                               cdhit_identity = 0.9,
                               cdhit_word_length = 5,
@@ -1747,6 +1920,7 @@ runDataProcessing <- function(duckdb_path,
                               # Metadata cleaning
                               ref_file_path = "data_raw/",
                               verbose = TRUE) {
+  panaroo_refind_mode <- match.arg(panaroo_refind_mode)
   duckdb_path <- normalizePath(duckdb_path)
   out_dir <- if (is.null(output_path)) dirname(duckdb_path) else normalizePath(output_path)
 
@@ -1760,6 +1934,10 @@ runDataProcessing <- function(duckdb_path,
     family_seq_identity    = panaroo_family_seq_identity,
     threads                = threads,
     split_jobs             = panaroo_split_jobs,
+    refind_mode            = panaroo_refind_mode,
+    strip_pseudogenes      = panaroo_strip_pseudogenes,
+    pseudogene_clean_dir   = panaroo_pseudogene_clean_dir,
+    write_pseudogene_audit = panaroo_write_pseudogene_audit,
     verbose                = verbose
   )
 
@@ -1822,3 +2000,224 @@ runDataProcessing <- function(duckdb_path,
     parquet_duckdb_path = normalizePath(parquet_duckdb_path)
   ))
 }
+
+#' Export processed tables from DuckDB database
+#'
+#' Reads tables from the DuckDB database produced by the `runDataProcessing()` workflow
+#' and exports them as CSV, TSV, Parquet, and/or XLSX. This is an optional step
+#' that allows users to take their processed data outside our amR workflow for
+#' use in their own custom analyses. This is not required to run `amRml`!
+#'
+#' @param duckdb_path Character. Path to the DuckDB database created by the
+#'   workflow (for example, `Sar.duckdb`).
+#' @param output_path Character or NULL. Directory for exports. Defaults to
+#'   file.path(dirname(duckdb_path), "processed_exports").
+#' @param amr_phenotype_mode Character. One of "separate" or "append".
+#'   "separate" exports the AMR labels as a separate wide table.
+#'   "append" joins those labels onto the main feature tables before export.
+#' @param export_formats Character vector. Any of "csv", "tsv", "parquet", "xlsx".
+#' @param tables Character vector or NULL. Tables to export. If NULL, exports the
+#'   standard processed tables present in the database.
+#' @param verbose Logical. If TRUE, prints progress messages.
+#'
+#' @return Invisibly returns a list containing the export path, table names, and mode.
+#' @export
+exportProcessedData <- function(duckdb_path,
+                                output_path = NULL,
+                                amr_phenotype_mode = c("separate", "append"),
+                                export_formats = c("csv"),
+                                export_sequences = FALSE,
+                                tables = NULL,
+                                verbose = TRUE) {
+  duckdb_path <- normalizePath(duckdb_path, mustWork = TRUE)
+  amr_phenotype_mode <- match.arg(amr_phenotype_mode)
+
+  export_formats <- unique(tolower(export_formats))
+  allowed_formats <- c("csv", "tsv", "parquet", "xlsx")
+  unknown_formats <- setdiff(export_formats, allowed_formats)
+  if (length(unknown_formats)) {
+    stop("Unsupported export format(s): ", paste(unknown_formats, collapse = ", "))
+  }
+
+  if ("xlsx" %in% export_formats) {
+    message(
+      "Excel spreadsheets of these features can be extremely large and may not open properly even on powerful hardware."
+    )
+    if (!requireNamespace("writexl", quietly = TRUE)) {
+      stop("Format 'xlsx' was requested but package 'writexl' is not available.")
+    }
+  }
+
+  if (is.null(output_path)) {
+    output_path <- file.path(dirname(duckdb_path), "processed_exports")
+  }
+  output_path <- normalizePath(output_path, mustWork = FALSE)
+  dir.create(output_path, recursive = TRUE, showWarnings = FALSE)
+
+  con <- DBI::dbConnect(duckdb::duckdb(), dbdir = duckdb_path, read_only = TRUE)
+  DBI::dbExecute(
+    con,
+    sprintf(
+      "SET file_search_path='%s'",
+      dirname(normalizePath(duckdb_path))
+    )
+  )
+  on.exit(try(DBI::dbDisconnect(con, shutdown = TRUE), silent = TRUE), add = TRUE)
+
+  available_tables <- DBI::dbListTables(con)
+  if (!length(available_tables)) {
+    stop("No tables found in DuckDB: ", duckdb_path)
+  }
+
+  read_tbl <- function(tbl) {
+    tibble::as_tibble(DBI::dbReadTable(con, tbl))
+  }
+
+  write_one <- function(df, stem) {
+    if ("csv" %in% export_formats) {
+      readr::write_csv(df, file.path(output_path, paste0(stem, ".csv")), na = "")
+    }
+    if ("tsv" %in% export_formats) {
+      readr::write_tsv(df, file.path(output_path, paste0(stem, ".tsv")), na = "")
+    }
+    if ("parquet" %in% export_formats) {
+      arrow::write_parquet(df, file.path(output_path, paste0(stem, ".parquet")))
+    }
+    if ("xlsx" %in% export_formats) {
+      writexl::write_xlsx(list(data = df), file.path(output_path, paste0(stem, ".xlsx")))
+    }
+  }
+
+  build_amr_wide <- function() {
+    source_tbl <- if ("metadata" %in% available_tables) {
+      "metadata"
+    } else if ("amr_phenotype" %in% available_tables) {
+      "amr_phenotype"
+    } else {
+      NULL
+    }
+
+    if (is.null(source_tbl)) {
+      return(NULL)
+    }
+
+    md <- read_tbl(source_tbl)
+
+    needed <- c("genome.genome_id", "genome_drug.antibiotic", "genome_drug.resistant_phenotype")
+    if (!all(needed %in% names(md))) {
+      return(NULL)
+    }
+
+    md |>
+      dplyr::transmute(
+        genome_id = `genome.genome_id`,
+        antibiotic = `genome_drug.antibiotic`,
+        phenotype = `genome_drug.resistant_phenotype`
+      ) |>
+      dplyr::filter(!is.na(genome_id), !is.na(antibiotic), !is.na(phenotype)) |>
+      dplyr::distinct() |>
+      dplyr::group_by(genome_id, antibiotic) |>
+      dplyr::summarise(
+        phenotype = paste(sort(unique(phenotype)), collapse = ";"),
+        .groups = "drop"
+      ) |>
+      tidyr::pivot_wider(
+        names_from = antibiotic,
+        values_from = phenotype,
+        values_fill = NA_character_
+      ) |>
+      dplyr::arrange(genome_id)
+  }
+
+  # Appendable here means whether we glue AST phenotypes on the end or not
+  table_specs <- list(
+    gene_count = list(source = "gene_count", stem = "gene_count", appendable = TRUE),
+    protein_count = list(source = "protein_count", stem = "protein_count", appendable = TRUE),
+    domain_count = list(source = "domain_count", stem = "domain_count", appendable = TRUE),
+    struct = list(source = "struct", stem = "struct", appendable = TRUE),
+    gene_names = list(source = "gene_names", stem = "gene_names", appendable = FALSE),
+    protein_names = list(source = "protein_names", stem = "protein_names", appendable = FALSE),
+    domain_names = list(source = "domain_names", stem = "domain_names", appendable = FALSE),
+    metadata = list(source = "metadata", stem = "metadata", appendable = FALSE),
+    genome_data = list(source = "genome_data", stem = "genome_data", appendable = FALSE),
+    amr_phenotype_wide = list(source = NULL, stem = "amr_phenotype_wide", appendable = FALSE)
+  )
+
+  if (isTRUE(export_sequences)) {
+    table_specs$gene_seqs <- list(source = "gene_ref_seq", stem = "gene_seqs", appendable = FALSE)
+    table_specs$protein_seqs <- list(source = "protein_cluster_seq", stem = "protein_seqs", appendable = FALSE)
+    table_specs$genome_gene_protein <- list(source = "genome_gene_protein", stem = "genome_gene_protein", appendable = FALSE)
+  }
+
+  if (is.null(tables)) {
+    selected_keys <- c(
+      "gene_count",
+      "protein_count",
+      "domain_count",
+      "struct",
+      "gene_names",
+      "protein_names",
+      "domain_names",
+      "metadata",
+      "genome_data",
+      "amr_phenotype_wide"
+    )
+    if (isTRUE(export_sequences)) {
+      selected_keys <- c(selected_keys, "gene_seqs", "protein_seqs", "genome_gene_protein")
+    }
+  } else {
+    selected_keys <- intersect(as.character(tables), names(table_specs))
+  }
+
+  if (!length(selected_keys)) {
+    stop("No requested tables were found in the DuckDB.")
+  }
+
+  phenotype_wide <- build_amr_wide()
+  exported <- character(0)
+
+  for (key in selected_keys) {
+    spec <- table_specs[[key]]
+
+    if (key == "amr_phenotype_wide") {
+      if (is.null(phenotype_wide)) {
+        if (isTRUE(verbose)) message("Skipping amr_phenotype_wide: no AMR source table found.")
+        next
+      }
+      write_one(phenotype_wide, spec$stem)
+      exported <- c(exported, spec$stem)
+      if (isTRUE(verbose)) message("Exported: ", spec$stem)
+      next
+    }
+
+    if (is.null(spec$source) || !(spec$source %in% available_tables)) {
+      if (isTRUE(verbose)) message("Skipping missing table: ", key)
+      next
+    }
+
+    df <- read_tbl(spec$source)
+    out_stem <- spec$stem
+
+    if (identical(amr_phenotype_mode, "append") &&
+        isTRUE(spec$appendable) &&
+        !is.null(phenotype_wide) &&
+        "genome_id" %in% names(df)) {
+      df <- dplyr::left_join(df, phenotype_wide, by = "genome_id")
+      out_stem <- paste0(spec$stem, "_with_phenotypes")
+    }
+
+    write_one(df, out_stem)
+    exported <- c(exported, out_stem)
+
+    if (isTRUE(verbose)) message("Exported: ", out_stem)
+  }
+
+  invisible(list(
+    duckdb_path = duckdb_path,
+    output_path = output_path,
+    tables = exported,
+    amr_phenotype_mode = amr_phenotype_mode,
+    export_formats = export_formats,
+    export_sequences = isTRUE(export_sequences)
+  ))
+}
diff --git a/README.Rmd b/README.Rmd
index 413bc1d..ec8e4e9 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -270,6 +270,8 @@ Processing times vary by species and isolate count:
 - These numbers will all vary greatly based on isolate number, genome complexity, and available hardware.
 - Parallelization significantly reduces processing time when multiple cores are available.
 - If a `future::multisession` error occurs mid-run (e.g. while testing via `devtools::load_all()` before installing the package), restart your R session fully before retrying. An orphaned background worker process can leave a stale lock on the local DuckDB caches (e.g. `data/bvbrc/bvbrcData.duckdb`), which can produce inconsistent results on the next run that look like a data or QC bug but are actually just leftover session state.
+- If a `furrr`/`future::multisession` worker fails with `could not find function ".xxx"` for an internal amRdata helper, your installed copy of amRdata is stale relative to the source you're editing. `future::multisession` workers are fresh R processes that resolve `amRdata` by loading the *installed* package from `.libPaths()` — they do not see changes made only via `devtools::load_all()` in your interactive session. Run `devtools::install()` (or `pkgbuild::compile_dll(); devtools::document(); devtools::install()`) before exercising any function that runs work via `future`/`furrr`, or temporarily set `future::plan(future::sequential)` while iterating with `load_all()` alone.
+- `runDataProcessing()`/`runPanaroo2Duckdb()` default `panaroo_refind_mode` to `"off"` rather than Panaroo's own default. Refinding recovers gene calls that annotation tools missed, but its search can take substantially longer (or in rare cases fail to complete within hours) when a genome carries a cluster of CDS with internal stop codons — a condition existing genome-quality metadata (CheckM completeness/contamination, consistency scores, quality flags) does not flag. This default trades some gene-recovery accuracy for predictable runtime until a QC step upstream can screen out affected genomes; set `panaroo_refind_mode = "default"` to restore Panaroo's normal behavior.
 
 ### Integration with amR suite
 
diff --git a/README.md b/README.md
index e28a4fe..369719a 100644
--- a/README.md
+++ b/README.md
@@ -291,6 +291,30 @@ Processing times vary by species and isolate count:
   results on the next run that look like a data or QC bug but are
   actually just leftover session state.
 
+- If a `furrr`/`future::multisession` worker fails with
+  `could not find function ".xxx"` for an internal amRdata helper,
+  your installed copy of amRdata is stale relative to the source
+  you're editing. `future::multisession` workers are fresh R
+  processes that resolve `amRdata` by loading the *installed* package
+  from `.libPaths()` — they do not see changes made only via
+  `devtools::load_all()` in your interactive session. Run
+  `devtools::install()` before exercising any function that runs work
+  via `future`/`furrr`, or temporarily set
+  `future::plan(future::sequential)` while iterating with `load_all()`
+  alone.
+
+- `runDataProcessing()`/`runPanaroo2Duckdb()` default
+  `panaroo_refind_mode` to `"off"` rather than Panaroo's own default.
+  Refinding recovers gene calls that annotation tools missed, but its
+  search can take substantially longer (or in rare cases fail to
+  complete within hours) when a genome carries a cluster of CDS with
+  internal stop codons — a condition existing genome-quality metadata
+  (CheckM completeness/contamination, consistency scores, quality
+  flags) does not flag. This default trades some gene-recovery
+  accuracy for predictable runtime until a QC step upstream can screen
+  out affected genomes; set `panaroo_refind_mode = "default"` to
+  restore Panaroo's normal behavior.
+
 ### Integration with amR suite
 
 amRdata is designed to work seamlessly with other amR packages:
diff --git a/inst/extdata/Campy_testdata.zip b/inst/extdata/Campy_testdata.zip
new file mode 100644
index 0000000..32a4e05
Binary files /dev/null and b/inst/extdata/Campy_testdata.zip differ
diff --git a/man/dot-apply_metadata_qc.Rd b/man/dot-apply_metadata_qc.Rd
new file mode 100644
index 0000000..d7f7622
--- /dev/null
+++ b/man/dot-apply_metadata_qc.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data_curation.R
+\name{.apply_metadata_qc}
+\alias{.apply_metadata_qc}
+\title{Apply QC filters to BV-BRC genome metadata}
+\usage{
+.apply_metadata_qc(
+  genome_tbl,
+  max_checkm_contam = 5,
+  min_checkm_complete = 95,
+  gc_deviations = NULL,
+  length_deviations = NULL,
+  cds_deviations = NULL
+)
+}
+\arguments{
+\item{genome_tbl}{A tibble of genome metadata containing BV-BRC genome columns.}
+
+\item{max_checkm_contam}{Numeric scalar. Maximum allowed CheckM contamination (\%).}
+
+\item{min_checkm_complete}{Numeric scalar. Minimum allowed CheckM completeness (\%).}
+
+\item{gc_deviations}{Optional numeric scalar. Maximum SDs from the median GC content.}
+
+\item{length_deviations}{Optional numeric scalar. Maximum SDs from the median genome length.}
+
+\item{cds_deviations}{Optional numeric scalar. Maximum SDs from the median CDS count.}
+}
+\value{
+A list with components:
+\describe{
+\item{qc_tbl}{QC-annotated genome table}
+\item{keep_ids}{Character vector of retained genome IDs}
+\item{rejections}{Tibble describing dropped genomes and reasons}
+}
+}
+\description{
+Quality control step with default values for contamination and completeness.
+Optional filters for other genome stats to find genomes that are X deviations
+away from the median values for those stats. This shouldn't be applied to
+jobs that use more than one diverse taxon at a time though! Requires a fairly
+normal distribution like you'd see for (most) single species
+Returns the cleaned table, retained genome IDs, and a rejection log.
+}
+\details{
+CheckM contamination and completeness are treated as the primary QC gate.
+Deviation-based filters are optional and are best suited to single-taxon runs.
+}
+\keyword{internal}
diff --git a/man/dot-ftpes_download_two_pass.Rd b/man/dot-ftpes_download_two_pass.Rd
deleted file mode 100644
index e15d299..0000000
--- a/man/dot-ftpes_download_two_pass.Rd
+++ /dev/null
@@ -1,22 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/data_curation.R
-\name{.ftpes_download_two_pass}
-\alias{.ftpes_download_two_pass}
-\title{Helps manage FTPS downloading from BV-BRC, tryng a quick download first, and
-if that fails, trying a longer timeout 2nd pass at the end in case it was a
-hiccup. If 2nd pass fails, log and give up on that file.}
-\usage{
-.ftpes_download_two_pass(
-  genome_ids,
-  out_dir,
-  workers_first = 4L,
-  workers_second = 4L,
-  log_file = NULL
-)
-}
-\description{
-Helps manage FTPS downloading from BV-BRC, tryng a quick download first, and
-if that fails, trying a longer timeout 2nd pass at the end in case it was a
-hiccup. If 2nd pass fails, log and give up on that file.
-}
-\keyword{internal}
diff --git a/man/dot-ftpes_download_one.Rd b/man/dot-ftps_download_one.Rd
similarity index 88%
rename from man/dot-ftpes_download_one.Rd
rename to man/dot-ftps_download_one.Rd
index b89f5d4..af5861d 100644
--- a/man/dot-ftpes_download_one.Rd
+++ b/man/dot-ftps_download_one.Rd
@@ -1,11 +1,11 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/data_curation.R
-\name{.ftpes_download_one}
-\alias{.ftpes_download_one}
+\name{.ftps_download_one}
+\alias{.ftps_download_one}
 \title{Helps appropriately interface with BV-BRC FTPS server, and avoids getting stuck
 when malformed files can hang an FTPS connection by introducing safeguards}
 \usage{
-.ftpes_download_one(
+.ftps_download_one(
   genomeID,
   out_dir,
   connect_timeout = 10L,
diff --git a/man/dot-ftps_download_two_pass.Rd b/man/dot-ftps_download_two_pass.Rd
new file mode 100644
index 0000000..0acd93f
--- /dev/null
+++ b/man/dot-ftps_download_two_pass.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data_curation.R
+\name{.ftps_download_two_pass}
+\alias{.ftps_download_two_pass}
+\title{Helps manage FTPS downloading from BV-BRC, trying a quick download first, and
+if that fails, trying a longer timeout 2nd pass at the end in case it was a
+hiccup. If the 2nd pass fails, log and give up on that file.}
+\usage{
+.ftps_download_two_pass(
+  genome_ids,
+  out_dir,
+  workers_first = 8L,
+  workers_second = 8L,
+  log_file = NULL
+)
+}
+\description{
+Helps manage FTPS downloading from BV-BRC, trying a quick download first, and
+if that fails, trying a longer timeout 2nd pass at the end in case it was a
+hiccup. If the 2nd pass fails, log and give up on that file.
+}
+\keyword{internal}
diff --git a/man/dot-parse_bvbrc_tsv.Rd b/man/dot-parse_bvbrc_tsv.Rd
new file mode 100644
index 0000000..049d0d8
--- /dev/null
+++ b/man/dot-parse_bvbrc_tsv.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data_curation.R
+\name{.parse_bvbrc_tsv}
+\alias{.parse_bvbrc_tsv}
+\title{Parse BV-BRC TSV output}
+\usage{
+.parse_bvbrc_tsv(x)
+}
+\arguments{
+\item{x}{Character vector of output lines from a BV-BRC CLI command.}
+}
+\value{
+A tibble with minimal name repair.
+}
+\description{
+Safely cleans blank links and weird formatting for a safe merge of the BV-BRC
+TSV data into a tibble.
+}
+\keyword{internal}
diff --git a/man/dot-processPanaroo.Rd b/man/dot-processPanaroo.Rd
index ecb74b9..aa42f1f 100644
--- a/man/dot-processPanaroo.Rd
+++ b/man/dot-processPanaroo.Rd
@@ -11,7 +11,8 @@
   len_dif_percent,
   cluster_threshold,
   family_seq_identity,
-  panaroo_threads_per_job
+  panaroo_threads_per_job,
+  refind_mode = c("off", "default", "strict")
 )
 }
 \arguments{
@@ -28,6 +29,16 @@
 \item{family_seq_identity}{Numeric. Gene family clustering identity (\code{-f}). Default \code{0.5}.}
 
 \item{panaroo_threads_per_job}{Integer. Number of threads for Panaroo and parallel execution.}
+
+\item{refind_mode}{Character. Panaroo's \code{--refind-mode} (\code{"off"}, \code{"default"}, or
+\code{"strict"}). Refinding searches for and recovers gene calls that annotation
+tools missed, comparing each candidate against the rest of the pangenome.
+Caveat: this search can take substantially longer (and in rare cases fail to
+complete within hours) when a genome carries a cluster of CDS with internal
+stop codons, which existing upstream genome-quality fields do not flag.
+Default \code{"off"} for now, to avoid that runtime risk; plan to move this back to
+\code{"default"} once a QC step upstream (e.g. in \code{.apply_metadata_qc()}) can screen
+out affected genomes before they reach Panaroo.}
 }
 \value{
 A list of results for each Panaroo batch in its output directory.
diff --git a/man/dot-runPanaroo.Rd b/man/dot-runPanaroo.Rd
index df64874..ed083fa 100644
--- a/man/dot-runPanaroo.Rd
+++ b/man/dot-runPanaroo.Rd
@@ -12,7 +12,11 @@
   cluster_threshold = 0.95,
   family_seq_identity = 0.5,
   threads = 8,
-  split_jobs = FALSE
+  split_jobs = FALSE,
+  refind_mode = c("off", "default", "strict"),
+  strip_pseudogenes = FALSE,
+  pseudogene_clean_dir = "gff_clean",
+  write_pseudogene_audit = TRUE
 )
 }
 \arguments{
@@ -32,6 +36,10 @@
 
 \item{split_jobs}{Logical. If TRUE, split into multiple smaller pangenome
 generation jobs that can be merged by \code{\link[=.mergePanaroo]{.mergePanaroo()}}. If FALSE, all isolates in one run.}
+
+\item{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"}.}
 }
 \value{
 A list of results for each Panaroo batch in its output directory.
diff --git a/man/dot-stripPseudogeneGFFs.Rd b/man/dot-stripPseudogeneGFFs.Rd
new file mode 100644
index 0000000..f2e4471
--- /dev/null
+++ b/man/dot-stripPseudogeneGFFs.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data_processing.R
+\name{.stripPseudogeneGFFs}
+\alias{.stripPseudogeneGFFs}
+\title{Remove pseudogene annotations from Panaroo input GFF files}
+\usage{
+.stripPseudogeneGFFs(panaroo_input_files, output_path, clean_dir = "gff_clean")
+}
+\arguments{
+\item{panaroo_input_files}{Character vector of \code{"gff fna"} input lines used
+by Panaroo.}
+
+\item{output_path}{Character scalar. Base directory for temporary cleaned
+GFF files and audit outputs.}
+
+\item{clean_dir}{Character scalar. Name of the subdirectory created beneath
+\code{output_path} to store cleaned GFF files. Default \code{"gff_clean"}.}
+}
+\value{
+A list containing:
+\itemize{
+\item \code{panaroo_input_files} — rewritten Panaroo input lines pointing to the
+cleaned GFF files.
+\item \code{audit} — a tibble summarizing, for each genome, the total number of
+annotated features, the number of pseudogenes removed, and the number of
+remaining features.
+}
+}
+\description{
+Cleans GFF annotation files of \code{pseudogene} feature records only. Cleaned
+GFFs are written to a subdirectory under \code{output_path} and swapped into the
+Panaroo input list, leaving the original genome annotations alone.
+}
+\details{
+This optional preprocessing step can reduce weird runtime stalls during
+Panaroo graph construction for some BV-BRC/PATRIC genome annotations that
+contain troublesome pseudogene features.
+
+This performs lightweight preprocessing only, removing feature records whose
+third GFF column is exactly \code{"pseudogene"} and does not otherwise modify
+annotation coordinates, attributes, or sequence files. FASTA paths are unchanged.
+}
+\keyword{internal}
diff --git a/man/exportProcessedData.Rd b/man/exportProcessedData.Rd
new file mode 100644
index 0000000..c8f16d4
--- /dev/null
+++ b/man/exportProcessedData.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data_processing.R
+\name{exportProcessedData}
+\alias{exportProcessedData}
+\title{Export processed tables from DuckDB database}
+\usage{
+exportProcessedData(
+  duckdb_path,
+  output_path = NULL,
+  amr_phenotype_mode = c("separate", "append"),
+  export_formats = c("csv"),
+  export_sequences = FALSE,
+  tables = NULL,
+  verbose = TRUE
+)
+}
+\arguments{
+\item{duckdb_path}{Character. Path to the DuckDB database created by the
+workflow (for example, \code{Sar.duckdb}).}
+
+\item{output_path}{Character or NULL. Directory for exports. Defaults to
+file.path(dirname(duckdb_path), "processed_exports").}
+
+\item{amr_phenotype_mode}{Character. One of "separate" or "append".
+"separate" exports the AMR labels as a separate wide table.
+"append" joins those labels onto the main feature tables before export.}
+
+\item{export_formats}{Character vector. Any of "csv", "tsv", "parquet", "xlsx".}
+
+\item{tables}{Character vector or NULL. Tables to export. If NULL, exports the
+standard processed tables present in the database.}
+
+\item{verbose}{Logical. If TRUE, prints progress messages.}
+}
+\value{
+Invisibly returns a list containing the export path, table names, and mode.
+}
+\description{
+Reads tables from the DuckDB database produced by the \code{runDataProcessing()} workflow
+and exports them as CSV, TSV, Parquet, and/or XLSX. This is an optional step
+that allows users to take their processed data outside our amR workflow for
+use in their own custom analyses. This is not required to run \code{amRml}!
+}
diff --git a/man/exportTables.Rd b/man/exportTables.Rd
new file mode 100644
index 0000000..0622fe7
--- /dev/null
+++ b/man/exportTables.Rd
@@ -0,0 +1,45 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data_curation.R
+\name{exportTables}
+\alias{exportTables}
+\title{Export DuckDB tables and optionally load them into R}
+\usage{
+exportTables(
+  duckdb_path,
+  output_dir = NULL,
+  tables = NULL,
+  skip_tables = NULL,
+  include_summary = TRUE,
+  export_tables = TRUE,
+  load_tables = FALSE,
+  verbose = TRUE
+)
+}
+\arguments{
+\item{duckdb_path}{Character. Path to the DuckDB file.}
+
+\item{output_dir}{Character or NULL. Directory for exports. Defaults to
+file.path(dirname(duckdb_path), "exports").}
+
+\item{tables}{Character vector or NULL. Tables to export. If NULL, exports
+all tables in the database.}
+
+\item{skip_tables}{Character vector of table names to exclude.}
+
+\item{include_summary}{Logical. If TRUE, writes summary.csv and summary.txt.}
+
+\item{export_tables}{Logical. If TRUE, writes tables to CSV files in \code{output_dir}.}
+
+\item{load_tables}{Logical. If TRUE, also return the exported tables as
+in-memory R data frames.}
+
+\item{verbose}{Logical. If TRUE, prints progress messages.}
+}
+\value{
+Invisibly returns a list with exported file paths and, if requested,
+in-memory tables.
+}
+\description{
+Writes selected DuckDB tables to CSV files and optionally returns them as
+in-memory R data frames.
+}
diff --git a/man/genomeList.Rd b/man/genomeList.Rd
index dfa6f39..ea36ef2 100644
--- a/man/genomeList.Rd
+++ b/man/genomeList.Rd
@@ -4,7 +4,7 @@
 \alias{genomeList}
 \title{Build a table of local genome file paths and write to DuckDB}
 \usage{
-genomeList(base_dir = ".", user_bacs, verbose = TRUE)
+genomeList(base_dir = ".", user_bacs, expected_ids = NULL, verbose = TRUE)
 }
 \arguments{
 \item{base_dir}{Character. Project root.}
diff --git a/man/prepareGenomes.Rd b/man/prepareGenomes.Rd
index 29edd07..0f7814a 100644
--- a/man/prepareGenomes.Rd
+++ b/man/prepareGenomes.Rd
@@ -10,13 +10,22 @@ prepareGenomes(
   base_dir = ".",
   method = c("ftp", "cli"),
   overwrite = FALSE,
+  num_workers = 8L,
   evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"),
+  max_checkm_contam = 5,
+  min_checkm_complete = 95,
+  gc_deviations = NULL,
+  length_deviations = NULL,
+  cds_deviations = NULL,
+  export_tables = FALSE,
+  load_tables = FALSE,
+  debug = FALSE,
   verbose = TRUE
 )
 }
 \arguments{
 \item{user_bacs}{Character vector. Species and/or taxon IDs (e.g.
-\code{c("Shigella flexneri", "623")}).}
+\code{c("Shigella sonnei", "624")}).}
 
 \item{genome_id_file}{Character or NULL. Optional path to a file listing genome
 IDs (one per line), passed through to \code{retrieveMetadata()}. If provided, the
@@ -31,10 +40,25 @@ metadata step is restricted to these genome IDs instead of resolving them from
 \item{overwrite}{Logical. Passed to metadata filtering and DuckDB creation.
 Default FALSE.}
 
+\item{num_workers}{Integer. Parallel workers used for genome download.
+Applied to both FTP and CLI download branches. Default: 8.}
+
 \item{evidence_mode}{Character. Sets what types of AMR evidence is acceptable.
 Default \code{lab_only}. \code{any} will not require AMR data for downloads. This will
 return very large download lists for many species!}
 
+\item{max_checkm_contam}{Numeric scalar. Maximum allowed CheckM contamination (\%).}
+
+\item{min_checkm_complete}{Numeric scalar. Minimum allowed CheckM completeness (\%).}
+
+\item{gc_deviations}{Optional numeric scalar. Maximum SDs from the median GC content.}
+
+\item{length_deviations}{Optional numeric scalar. Maximum SDs from the median genome length.}
+
+\item{cds_deviations}{Optional numeric scalar. Maximum SDs from the median CDS count.}
+
+\item{debug}{Logical. If TRUE, keep QC columns in metadata tables.}
+
 \item{verbose}{Logical. Print progress messages. Default TRUE.}
 }
 \value{
diff --git a/man/retrieveGenomes.Rd b/man/retrieveGenomes.Rd
index 1800b5d..f9d202b 100644
--- a/man/retrieveGenomes.Rd
+++ b/man/retrieveGenomes.Rd
@@ -10,9 +10,9 @@ retrieveGenomes(
   method = c("ftp", "cli"),
   image = "danylmb/bvbrc:5.3",
   skip_existing = TRUE,
-  ftp_workers = 4L,
-  cli_fasta_workers = 4L,
-  cli_gff_workers = 4L,
+  ftp_workers = 8L,
+  cli_fasta_workers = 8L,
+  cli_gff_workers = 8L,
   chunk_size = 50L,
   evidence_mode = c("lab_only", "lab_or_comp", "comp_only", "any"),
   verbose = TRUE
@@ -31,9 +31,9 @@ retrieveGenomes(
 
 \item{ftp_workers}{Parallel workers for FTP path (default 8).}
 
-\item{cli_fasta_workers}{Parallel chunk containers for FASTA+GTO (default 4).}
+\item{cli_fasta_workers}{Parallel chunk containers for FASTA+GTO (default 8).}
 
-\item{cli_gff_workers}{Parallel chunk containers for GFF export (default 4).}
+\item{cli_gff_workers}{Parallel chunk containers for GFF export (default 8).}
 
 \item{chunk_size}{Genomes per chunk container (default 50).}
 
diff --git a/man/retrieveMetadata.Rd b/man/retrieveMetadata.Rd
index bed8547..90d8cdb 100644
--- a/man/retrieveMetadata.Rd
+++ b/man/retrieveMetadata.Rd
@@ -12,6 +12,14 @@ retrieveMetadata(
   abx = "All",
   overwrite = FALSE,
   image = "danylmb/bvbrc:5.3",
+  max_checkm_contam = 5,
+  min_checkm_complete = 95,
+  gc_deviations = NULL,
+  length_deviations = NULL,
+  cds_deviations = NULL,
+  debug = FALSE,
+  export_tables = FALSE,
+  load_tables = FALSE,
   verbose = TRUE
 )
 }
@@ -34,7 +42,19 @@ Default NULL.}
 
 \item{image}{Character. Docker image. Default "danylmb/bvbrc:5.3".}
 
-\item{verbose}{Logical. If TRUE, prints concise messages.}
+\item{max_checkm_contam}{Numeric scalar. Maximum allowed CheckM contamination (\%).}
+
+\item{min_checkm_complete}{Numeric scalar. Minimum allowed CheckM completeness (\%).}
+
+\item{gc_deviations}{Optional numeric scalar. Maximum SDs from the median GC content.}
+
+\item{length_deviations}{Optional numeric scalar. Maximum SDs from the median genome length.}
+
+\item{cds_deviations}{Optional numeric scalar. Maximum SDs from the median CDS count.}
+
+\item{debug}{Logical. If TRUE, retain \code{metadata_full} and QC columns for inspection.}
+
+\item{verbose}{Logical. If TRUE, print progress messages.}
 }
 \value{
 A list with:
diff --git a/man/runDataProcessing.Rd b/man/runDataProcessing.Rd
index 4c1f444..fa2eda4 100644
--- a/man/runDataProcessing.Rd
+++ b/man/runDataProcessing.Rd
@@ -7,12 +7,16 @@
 runDataProcessing(
   duckdb_path,
   output_path = NULL,
-  threads = 16,
+  threads = 8,
   panaroo_split_jobs = FALSE,
   panaroo_core_threshold = 0.9,
   panaroo_len_dif_percent = 0.95,
   panaroo_cluster_threshold = 0.95,
   panaroo_family_seq_identity = 0.5,
+  panaroo_refind_mode = c("off", "default", "strict"),
+  panaroo_strip_pseudogenes = FALSE,
+  panaroo_pseudogene_clean_dir = "gff_clean",
+  panaroo_write_pseudogene_audit = TRUE,
   cdhit_identity = 0.9,
   cdhit_word_length = 5,
   cdhit_memory = 0,
@@ -38,7 +42,7 @@ download steps (e.g., \code{files}, \code{filtered}, and metadata tables).}
 outputs and final Parquet files. 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{16}.}
+InterProScan). Passed through to each stage as appropriate. 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}.}
@@ -51,6 +55,10 @@ merged by \code{\link[=.mergePanaroo]{.mergePanaroo()}}. If \code{FALSE}, Panaro
 
 \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{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}.}
@@ -181,7 +189,7 @@ files to \code{output_path}, and builds a \strong{Parquet-backed DuckDB}
 runDataProcessing(
   duckdb_path   = "data/Shigella_flexneri/Sfl.duckdb",
   output_path   = "data/Shigella_flexneri",
-  threads       = 16,
+  threads       = 8,
   ref_file_path = "data_raw/"
 )
 
diff --git a/man/runPanaroo2Duckdb.Rd b/man/runPanaroo2Duckdb.Rd
index 94b0467..48bec57 100644
--- a/man/runPanaroo2Duckdb.Rd
+++ b/man/runPanaroo2Duckdb.Rd
@@ -11,8 +11,12 @@ runPanaroo2Duckdb(
   len_dif_percent = 0.95,
   cluster_threshold = 0.95,
   family_seq_identity = 0.5,
-  threads = 16,
+  threads = 8,
   split_jobs = FALSE,
+  refind_mode = c("off", "default", "strict"),
+  strip_pseudogenes = FALSE,
+  pseudogene_clean_dir = "gff_clean",
+  write_pseudogene_audit = TRUE,
   verbose = TRUE
 )
 }
@@ -38,13 +42,17 @@ Default: \code{0.5}.}
 
 \item{threads}{Integer. Total CPU budget to allocate for Panaroo.
 If \code{split_jobs = TRUE}, threads are divided across batches.
-Default: \code{16}.}
+Default: \code{8}.}
 
 \item{split_jobs}{Logical. If \code{TRUE}, Panaroo is run in multiple parallel
 batches (up to 5, depending on dataset size), and batch outputs are merged
 using \code{.mergePanaroo()}. If \code{FALSE}, only one Panaroo invocation is run.
 Default: \code{FALSE}.}
 
+\item{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{verbose}{Logical. Print status messages during Panaroo execution,
 merging, and DuckDB import. Default: \code{TRUE}.}
 }
@@ -104,7 +112,7 @@ and modeling steps in \code{amRdata} and \code{amRml}.
 runPanaroo2Duckdb(
   duckdb_path = "data/Shigella_flexneri/Sfl.duckdb",
   output_path = "data/Shigella_flexneri",
-  threads     = 12,
+  threads     = 8,
   split_jobs  = FALSE
 )
 
diff --git a/tests/testthat/test-apply-metadata-qc.R b/tests/testthat/test-apply-metadata-qc.R
index 620a922..24dee5e 100644
--- a/tests/testthat/test-apply-metadata-qc.R
+++ b/tests/testthat/test-apply-metadata-qc.R
@@ -25,7 +25,7 @@ test_that(".apply_metadata_qc() keeps a genome that passes the default CheckM ga
 test_that(".apply_metadata_qc() drops a genome with contamination above the threshold", {
   tbl <- .mock_genome_tbl("511145.12", checkm_completeness = "99", checkm_contamination = "12")
 
-  out <- .apply_metadata_qc(tbl, checkm_contam = 5, checkm_complete = 95)
+  out <- .apply_metadata_qc(tbl, max_checkm_contam = 5, min_checkm_complete = 95)
 
   expect_length(out$keep_ids, 0L)
   expect_identical(out$rejections$failed_rule, "checkm_contamination")
@@ -35,7 +35,7 @@ test_that(".apply_metadata_qc() drops a genome with contamination above the thre
 test_that(".apply_metadata_qc() drops a genome with completeness below the threshold", {
   tbl <- .mock_genome_tbl("511145.12", checkm_completeness = "80", checkm_contamination = "1")
 
-  out <- .apply_metadata_qc(tbl, checkm_contam = 5, checkm_complete = 95)
+  out <- .apply_metadata_qc(tbl, max_checkm_contam = 5, min_checkm_complete = 95)
 
   expect_length(out$keep_ids, 0L)
   expect_identical(out$rejections$failed_rule, "checkm_completeness")
diff --git a/vignettes/intro.Rmd b/vignettes/intro.Rmd
index 9d9a36f..59a3827 100644
--- a/vignettes/intro.Rmd
+++ b/vignettes/intro.Rmd
@@ -113,6 +113,7 @@ runDataProcessing(
   panaroo_len_dif_percent = 0.95,
   panaroo_cluster_threshold = 0.95,
   panaroo_family_seq_identity = 0.5,
+  panaroo_refind_mode = "off",
   # CD-HIT
   cdhit_identity = 0.9,
   cdhit_word_length = 5,