Replacing InterProScan with HMMER + .json manifest

[epbrenner]

Lots of changes.

1. Moved many common and utility functions to helpers.R
2. Added new manifest feature, logging all steps of workflow

• Manifest is also used to identify which HMMER databases have run on a bug.
• Adds extensive data provenance logging + moves databases to R user cache.

3. Removes InterProScan entirely
4. Adds Abhirupa's HMMER functionality to data_processing.R, removed runHMMER.R
5. Replaced "domain" tables/features with Pfam/COG/AMRFinder/DefenseCas
6. Added jsonlite to Imports
7. stuff

To test:
source("R/helpers.R")
source("R/data_curation.R")
source("R/data_processing.R")

prepareGenomes("Staphylococcus epidermidis")
runDataProcessing("data/Staphylococcus_epidermidis/Sep.duckdb")

In the data/Staphylococcus_epidermidis directory, you will also find a .json manifest generated that logs a bunch of stuff. Check it out.

[epbrenner]

To-do:

• Need to update and standardize naming of output HMMER hit tables for downstream use
• Need to implement thorough statistical cutoffs for HMMER hits
• Test alternatives to the default workflow and see what breaks (use CLI, pick HMMER database subsets, etc.)
• Documentation!

[AbhirupaGhosh]

I ran

source("R/data_curation.R")
source("R/data_processing.R")
source("R/data_helpers.R")
prepareGenomes("Staphylococcus epidermidis")
runDataProcessing(duckdb_path = "data/Staphylococcus_epidermidis/Sep.duckdb")

The manifest JSON is generated with all the local and global parameters.

[jananiravi]

@epbrenner, convert this to a 'real' PR? @AbhirupaGhosh has already approved this. @eboyer221 and I can take a quick look, too. @AbhirupaGhosh would like to merge this with main and then move forward with #45.

[epbrenner]

Okay, will merge after a second review. I do think the items on that to-do list above are important, but perhaps not blockers for merge.

[AbhirupaGhosh]

Okay, will merge after a second review. I do think the items on that to-do list above are important, but perhaps not blockers for merge.

NEW PR 🚀

[eboyer221]

Tested this locally against commit f708e7a (PR HEAD) and installed the package fresh, started Docker, ran the test suite, and did a live test pulling real genome data from BV-BRC. The core BV-BRC/manifest work seemed perfect in my test, genome metadata retrieval and the DuckDB writing all worked correctly.

I ran into a few problems in the new HMMER code below.

A few fixes that I believe are high priority before merging:

1. Running with anything other than all four databases at once will crash
R/data_processing.R, around lines 2986–3015 and 3084–3085.

If you call runDataProcessing() with, say, just "DefenseCas" (skipping Pfam/COG/AMRFinder), it crashes with object 'hmmer_result' not found. The same thing happens in reverse if DefenseCas is left out. This happens after the Panaroo/CD-HIT/HMMER steps have already finished.

I think the two variables holding the results only get created inside the if block for their specific database group, but they're both checked later regardless of which group actually ran.

Fix: set both to NULL before those if blocks. I tested this change and confirmed it resolves the crash.

2. Every Pfam/COG/AMRFinder run loses its accession column
R/data_processing.R, around lines 1541–1543.

After a normal run, the Pfam/COG/AMRFinder result tables end up with two separate accession columns instead of one, because two different data sources both have a column with the same name, and they get combined without saying which one should win. This isn't an edge case, it happens on every run.
I recreated the same setup on a small example and confirmed that dropping the duplicate before combining fixes it and gives one column back with the right values.

3. A search with zero hits crashes the whole run
R/helpers.R, line 866.

If a batch of proteins or a database (DefenseCas especially - a lot of genomes just won't have matches for it) comes back with zero hits, the code that reads the results throws an error instead of just returning nothing.
I confirmed the crash by feeding it a results file with no hits, then tested a fix (return an empty result instead of trying to process nothing) and confirmed that works both for the empty case and for a normal file with real hits in it.

4. New functions have no documentation, and a function people already call directly could stop working
man/ folder, and R/data_processing.R around line 2730.

None of the new manifest-related functions have documentation yet. Separately, runDataProcessing() (which the README and the in-app messages tell people to call directly) has lost the tag that makes it available to users. The docs file just hasn't been regenerated yet, so it still works for now, but the moment someone runs the doc-generation step (which needs to happen anyway to add docs for the new functions), runDataProcessing() would stop being callable.

Fix: add the tag back before regenerating docs. I tested this and confirmed it stays available afterward.

Smaller things, not urgent

• On a small dataset, a step could look for a file that was never created. The code sometimes splits data into fewer pieces than it plans for later, so with a small enough dataset, one step can end up looking for a piece that doesn't exist. I confirmed this with a small test case (21 items split into what should've been 8 pieces - only 7 got made). Probably won't come up on a typical full-size genome, but it's an easy fix.
• A checksum feature points at something that doesn't exist yet. There's a call to a function for generating file checksums that isn't actually available in base R. It's not being used anywhere right now, so it's not causing problems today, but it'll break as soon as anyone turns that feature on.
• One function's documentation is blank..isValidHmmFile() has documentation placeholders with nothing filled in.

[eboyer221]

Approving. I went back through this after the recent fixes, and everything I'd flagged earlier is taken care of: the crash that happened when you ran with anything other than all four databases at once, the bug where the accession column got split into two, the crash on results with zero hits, and the mismatch that could point at a file that was never created. I re-tested each one against the current code to make sure, and also re-ran the full test suite - still all passing, and the package installs without issues.

There's still documentation missing for a few of the newer internal functions, but that's already on the to-do list above and doesn't affect whether things work.

One heads-up for later: this now requires a small tracking file to exist next to the database before you can process it, which gets created automatically when you run prepareGenomes(). That means any data processed before this PR won't work with the new code as-is - not a problem for merging, just something to know if anyone tries to reuse older data.

[epbrenner]

Thanks for all your time reviewing this, Emily! Merging now and moving the to-do items to a follow-up PR to avoid holding up others.