Adding the sylph coverage model to yacht

[rtraborn]

Hi @dkoslicki and team!
I created a sylph coverage model from Shaw and Yu, 2024 and added it to yacht, in a branch I named superyacht just for fun.
This is a draft that I'm still testing, so that and other caveats still apply. A few notes:

• The coverage code is in a function called cov_calc, which calculates lambda and ani according as specified by the sylph paper.
• For engineering reasons, I decided to put cov_calc inside get_exclusive_hashes, given that that function provides us with the signature objects needed to make the calculations.
• Because of this, I am passing the output of cov_calc, a pandas dataframe, along with hypothesis_recovery. There are probably good ways to integrate this, and I'll give this some more thought.
• Also, I decided to not incorporate the output of cov_calc more deeply into hypothesis_recovery for now. I have some ideas on what might be the best approach that we could discuss if you'd like. I thought it would be best to share this new branch while I look into this more deeply.
• The script internal_superyacht_test.py is just a script that I have been using to test the new branch, and this can be ignored; I'll remove it once we move towards publication.
• I plan to update the way I instantiated the AdjustStatusLambda enum in a more idiomatic python way this week. It should be a relatively quick fix.
• I did not incorporate the taxonomic reassignment/winner_map routine from sylph, but it's something I would like to add.

I'm going to do more testing this week on additional datasets. Happy to discuss here or via email/video!

[rtraborn]

After some more testing, I just pushed some additional updates to this branch.

• Fixed a typo: corrected to logger.warning in cov_calc.py
• Added missing scipy.special.gamma import to utils.py
• Fixed a few bugs I discovered in binary_search_lambda()
• Replaced print statements with logger.info() for consistency
• Consolidated duplicate constants into utils.py
• Removed a few unused local variables from hypothesis_recovery_src.py

[rtraborn]

A small update, but with my most recent commit from last night I made the promised change to the AdjustStatusLambda enum in cov_calc to make it more idiomatically python-like. I think it looks cleaner- thanks to @standage for the flagging this!

[rtraborn]

Hi All!

I've made some more changes over the past week. Here's an overview of my most recent updates (Part 1 of 2):

• Formally adds the Winner takes all strategy from sylph, which performs abundance calculation and k-mer reassignment. This helps prevent double-counting of k-mers, assigning shared k-mers to the taxon with the highest calculated ANI. I'll note that in the interest of performance this procedure is only being done once per instance, rather than twice (as in sylph). We can discuss this in the future.
• This procedure tracks k-mers lost to reassignment (kmers_lost column) and filters orgs with final_est_ani < 0.90 (90% ANI threshold)
• Coverage results are now merged into the overall Excel output of yacht run. Previously the cov_calc output was passed along as a separate dataframe. I tried to do this without touching too much of the original yacht code; let me know what you think and if we need to make any tweaks.
• The columns added to this new output are naive_ani, final_est_ani, final_est_cov, mean_cov, median_cov, lambda_status, ani_ci, lambda_ci, rel_abund, kmers_lost
• Fixed a bug created by the refactoring described above 🙃 (MIN_ANI_THRESHOLD is defined).
• I fixed a bug that I encountered after more extensive testing. What happened was that ani_from_lambda function had a ZeroDivisionError when lambda_val was 0 or very close to 0. It took my testing on a pretty diverse metagenomic sample (more on this test later) for this to crop up, but I'm glad it did.

[sonarqubecloud]

Quality Gate passed

Issues
20 New issues
0 Accepted issues

Measures
0 Security Hotspots
0.0% Coverage on New Code
0.0% Duplication on New Code

See analysis details on SonarQube Cloud

[rtraborn]

Hi all! I made more updates to this branch (along with a decent amount of testing) as we discussed and the code is ready for your review.
Updates include the following:

• added new CLI argument (--winner-takes-all) for the winner-map functionality. Added two-pass winner-takes all.
• added ANI capping- prevents lambda-corrected ANI values from being higher than 1.
• added test that skipped correction under high-coverage conditions (coverage >= 3.0) as per sylph
• Most notably, added coverage calculation for yacht using the top-level --calculate-coverage CLI flag that replaces the user-supplied values for c. When run, instead of n worksheets, a single worksheet is produced in the (*.xlsx) output.
• (Minor) this adds a small fix to the gzip functionality when running yacht train- replaced python zgip with the system gzip. This corrects an error with incomplete uncompression that I encountered on Mac OSX.
  This was tested on about 5 standard metagenomic datasets with no errors or taxa misidentifications. I look forward to your comments and thoughts!
  [P.S. As an aside, the reason that the most recent SonarCloud scan failed appears to be because of regular yacht code that preexists this branch, so I refrained from pushing the fixes here. I can make those changes separately if desired.]

[rtraborn]

Hi team!
I just updated the code so that rel_abun is renormalized so that the column sums to 1. This required small changes to the appropriate locations of both hypothesis_recovery_src.py and run_YACHT.py.
While I was there I also corrected a variable name (it's now temp_manifest), which looked to be a typo.
This is now working on my tests. I'll run a few more today just to be sure.
[Update at 3:05pm: I just made another commit that renormalizes every worksheet when --winner_takes_all is on but --calculate_coverage is not. I encountered this during my testing.]

[rtraborn]

Good morning all!

Over the past few days I have made updates to standardize_yacht_output.py to make yacht convert compatible with the new "superyacht" --calculate_coverage and --winner_takes_call output.
Briefly, these changes include:

1. Small tweak to how organism_id_list is parsed to only grab the genome accession from gtdb genomes (where the accession and metadata are outputted into the same column).
2. Changing percentages to weights in yacht convert output (if the rel_abund column is present).
3. Rounding the contents of the PERCENTAGE column in the .cami output to 6 digits.
4. Adding a --single_sheet argument to yacht convert for handling superyacht output, as well as help statements. --sheet_name was made optional if --single_sheet is present.
5. Improvements to error handling if rel_abund is absent.
6. Adding appropriate warnings and info statements.

Note that I didn't change anything about the format of the .cami output file, including the PERCENTAGE column, as this appears to be required.

These updates to yacht convert work in my hands on different outputs.

Please let me know what you think, and if you have any questions or concerns.

[rtraborn]

Good afternoon all! I made a few more updates to a new version of superyacht, 2.2.0.
They include:

• adding a --min_ani flag to set MIN_ANI_THRESHOLD
• adding a Newton-Raphson convergence argument, --convergence_nr (the program still does not converge by default).
• Updated fallback coverage to include sample-wide median lambda to avoid being overly stringent (formerly used a coverage of 1)
• Fixed cov_max handling in one the branches (rarely encountered in my hands) in cov_calc.py

These updates give us more flexibility for usage and hopefully some more sensitivity.

[sonarqubecloud]

Quality Gate failed

Failed conditions
1 Security Hotspot

See analysis details on SonarQube Cloud

[rtraborn]

Good afternoon!

I completed a thorough optimization of the "Superyacht" codebase to improve performance, plus a few refinements.

The top-line summary is that this update (v 2.2.1) provides a significant speedup and requires less memory compared to its predecessors.

This optimization took the form of: a performance overhaul of winner-takes-all (WTA), which included:

1. parallelization of recalculate_ani
2. vectorization of build_winner_map, which replaced a per-k-mer Python loop that had slowed things down with a grouped array.
3. an array-based winner assignment, which replaced the massive (!) kmer: (ani, org) dictionary. This uses per-organism-won k-mer arrays with binary search. The result is a speedup and also reduced memory demands.
4. [Minor but meaningful]: an update in how kmers_lost is calculated by only considering organisms which participated in the winner map by having a valid (ie non-NaN) ANI. This results in a small but significant speedup.

Worth noting: instead of each genome's k-mers being extracted from the sketch four times, this extraction is only done once.

Overall, these changes resulted in a ~3.4X speedup in the back-half (that is, the Superyacht-derived part) of yacht run during my tests. This version also required less memory during my tests: from 55GB to 38.6GB running on a diverse CamiII metagenome.

One caveat to mention: to break ties, when two genomes have an exactly equal ANI for a given shared k-mer, the winner is selected arbitrarily via array order. This hasn't changed in this code, but it isn't deterministic because the results table is generated from separate workers. I can update it if we need to by making tie-breaking deterministic.

Finally, I left the timings I used to measure performance in the codebase but commented them out. This is only 20 lines of code overall and I can (and will) remove them once we're satisfied with the update.