Genome metadata QC filtering

[epbrenner]

Adding new functionality.

1. .parse_bvbrc_tsv() is a helper to more safely merge TSV data from parallel workers after I noticed a tab shift in data somewhere yielded strain names in the accession ID column and such after a failed merge.
2. .apply_metadata_qc is the big one, which allows you to screen BV-BRC genome metadata for CheckM completeness and CheckM contamination values. By default, I set these values very loose at 90% complete and 10% contamination being the cutoffs. We can adjust as needed. Additionally, you can enable filters based on CDS counts, genome length, and GC content, which are set based on deviations away from the median value of the distribution of those values for the input taxon. This will not work well for multi-taxon bug sets so this is set to NULL by default!
3. Updated the retrieveMetadata() and prepareGenomes() functions to accept the apply_metadata_qc steps.
4. Updated .ftpes_download_one to incorporate the old .is_complete helper inside the loop. This avoids a BiocParallel issue where the workers did not inherit the helpers.
5. Updated a spot where an old ftp_download_one function was called instead of ftpes_download_one

@amcim Most importantly for you, I know you've been working on the BiocParallel stuff and had since moved on to another implementation (futures?). I just want to confirm that this script is working for me locally after my patches, but this is just based off main. If we're changing implementations, some of my updates aren't necessary.

For testers:
Open GitHub Desktop or your preferred method of fetching a PR, fetch this PR, source data_curation.R into your favorite CLI, and try retrieveMetadata(user_bacs = "Shigella flexneri") for a quick bug, or prepareGenomes(user_bacs = "Shigella flexneri") to also test the genome downloader. I have tried this out locally for Enterobacter spp., K. pnumoniae, and Campylobacter jejuni. Give it a spin!

[epbrenner]

Update: user_bacs = "Staphylococcus argenteus" has a mere 73 genomes post-filtering, so it's good for very fast download tests!

[jananiravi]

Few comments on chunks that didn't need a review!

[eboyer221]

future, future.apply, and furrr are used in data_curation.R and data_processing.R (e.g. future::plan(), future.apply::future_lapply(), furrr::future_map()) but aren't declared in DESCRIPTION. All the calls in this PR are already ::-qualified, so I don't think any NAMESPACE/roxygen changes are needed. The suggestion is just to add these three to Imports:

​ Imports: arrow, BiocParallel, Biostrings, DBI, data.table, dplyr, duckdb, furrr, future, future.apply, glue, purrr, readr, reshape2, stringr, tibble, tidyr ​

While testing, I saw prepareGenomes("Staphylococcus argenteus") resolve to different genome counts across runs (248 genomes / 1 taxon ID vs. 21 genomes / 3 taxon IDs) for the identical input. The difference in numbers only showed up after an earlier run had hit a future::multisession error (from testing via devtools::load_all() before fully installing the branch. Restarting R fully resolved the issue, so it seems like an orphaned worker process left a stale lock on bvbrcData.duckdb that a later run then hit.

Might be worth a line in the docs/README recommending a full R session restart (not just load_all()) after any future::multisession error before re-testing, so this doesn't get mistaken for a genuine QC/parallelism bug by the next person who hits it.

A few smaller things I noticed while reviewing, none blocking:

• No test coverage yet for .apply_metadata_qc() or .parse_bvbrc_tsv() which are the two new core functions this PR adds. So we may want to add a couple of testthat cases for those in a separate PR.
• retrieveMetadata() defaults to checkm_contam = 10, checkm_complete = 90, while prepareGenomes() defaults to checkm_contam = 5, checkm_complete = 95 (stricter). Is that intentional, or should these match?
• The new is_complete() closure inside .ftpes_download_one() duplicates the existing .is_complete_set() helper (still used by .list_complete()/.missing_any()/.audit_gaps()/retrieveGenomes()). Could just call .is_complete_set() directly instead of re-implementing it.

[jananiravi]

good catch, Emily. @epbrenner happy to re-review once you make these changes.

[eboyer221]

good catch, Emily. @epbrenner happy to re-review once you make these changes.

Changes I pushed to this PR:
-Added future and furrr to DESCRIPTION Imports (they were used throughout but never declared so this would've broken on a fresh install)
-Fixed a scoping bug in .with_future_plan() (data_processing.R) where its own on.exit() reverted the plan before the caller's parallel map ever ran, silently killing parallelism in the Panaroo batch step and InterPro chunk processing
-Consolidated all remaining future.apply calls onto furrr::future_map()/future_map2_lgl() for consistency with the rest of the codebase's purrr usage, and to match the existing future.apply → dropped as a dependency
-Added a note in the README about the multisession error I encountered previously.

Verified: package loads clean, all existing tests pass.

[jananiravi]

I didn't notice anything from my quick manual inspection, but here are 3 from C -- flagging for @epbrenner @eboyer221. Please check to make sure these won't cause immediate or downstream issues.

1. prepareGenomes()'s broken return value (wrong key, no live connection, connection leak)
2. .ftps_download_two_pass()'s missing early-return + lost Pass-1 logging when everything succeeds first try
3. retrieveMetadata()'s missing "no drug data returned" guard.

Detailed version:
Three regressions in d658224:

1. prepareGenomes() return value is now broken. It used to return out (the actual result of genomeList(), containing a live duckdbConnection). It now returns a hand-built invisible(list(duckdb_path = paths$db_path, table_name = "files")) instead — wrong key name (duckdb_path vs. the documented duckdbConnection), and a path string instead of an open connection. prepareGenomes()'s own roxygen @return still promises "Active DBI connection to the per-bug DuckDB" — this commit breaks that contract. Any caller doing out$duckdbConnection (which the README's own workflow examples imply) will get NULL. It also silently leaks the connection genomeList() opened, since it's never captured or disconnected now. Looks like an overcorrection for "reduced console clutter for final 'out' print" — the fix should have been invisible(out), not replacing out's contents.
2. .ftps_download_two_pass() lost its "all passed in pass 1" early return. The old code had if (!length(fail_ids)) { ...; return(ok_ids_1) } right after Pass 1. That's gone now — so even when every genome succeeds in pass 1, the function still unconditionally spins up a second future::multisession worker pool and calls furrr::future_map() over an empty fail_ids. Not fatal (mapping over a length-0 vector is a no-op), but it's wasted worker-pool startup/teardown on every single successful run, and — given eboyer's finding that orphaned future::multisession workers are exactly what caused the stale-DuckDB-lock flakiness — adding an unnecessary extra worker-pool spin-up on every call cuts against the fix this whole PR is for. Same edit also dropped the "Pass 1: ok=X fail=Y" console message and its log_file line, so that pass's outcome is no longer logged at all.
3. retrieveMetadata() lost its "no drug data returned" guard. The old code had if (nrow(combined_drug_data_tbl) == 0L) { message("No drug data returned."); return(NULL) } right after building combined_drug_data_tbl. That block is gone (only the equivalent check for combined_genome_data_tbl survived). Now if BV-BRC returns zero AMR rows for a query, the function proceeds anyway, writes an empty amr_phenotype table, and the downstream INNER JOIN into metadata_full silently produces an empty result instead of a clear early message — trading a clear failure for a silent empty one.