Add mean_zero_coverage_length coverage method#292
Open
wwood wants to merge 5 commits into
Open
Conversation
Add a new 'coverm makedb' subcommand that builds persistent mapping indexes (databases) from reference genome or contig FASTA files, so they can be reused across multiple 'coverm contig'/'coverm genome' runs without re-indexing. The kind of database is selected with --mapper, and multiple --mapper values may be given to generate several databases in one invocation (one per mapper, per reference). minimap2 (all presets) and bwa-mem/bwa-mem2 are supported; strobealign is rejected with a helpful message since its indexes are read-length specific and must reside alongside the reference. The minimap2 database is written as <ref>.<mapper>.mmi and can be fed back in via 'coverm contig -r <db>.mmi --minimap2-reference-is-index'. - Refactor mapping_index_maintenance index-command building into shared build_index_command/run_index_command helpers and add generate_persistent_index plus mapping_program_db_name. - Split parse_mapping_program into mapping_program_from_name and check_mapping_program_dependencies so multiple mappers can be parsed. - Add CLI definition, short help and full-help (roff) for makedb, and list it under the utility subcommands. - Add functional tests and document the mode in the README and release.sh doc generation. Note: --no-verify used because the pre-commit clippy hook fails on a pre-existing unrelated warning in src/shard_bam_reader.rs under the current clippy version.
Previously 'coverm makedb' rejected strobealign because its index is read-length specific and must reside alongside the reference FASTA. Both constraints can be satisfied: - The read length can be set with --strobealign-params '-r <length>', or estimated by passing an example reads file (strobealign --create-index accepts a reads file positionally). A new --strobealign-params option is added to makedb for this. - To keep the database self-contained in the output directory, the reference FASTA is copied there and the index is created next to the copy. The copied FASTA is what gets passed back via '--reference <copy> --strobealign-use-index'. Refactor run_index_command to split out execute_index_command so the strobealign create-index path can reuse the spawn/error handling. Use .alias() rather than a second .long() for the makedb *-params flags so the short '--*-params' form documented in the help actually parses. Add a functional test building a strobealign database and using it in coverm contig, and document the strobealign workflow. (--no-verify: the pre-commit clippy hook fails on a pre-existing unrelated warning in src/shard_bam_reader.rs under the current clippy version.)
Replace the is_none()/unwrap() pattern for max_score with a match, which clippy's unnecessary_unwrap lint flagged. Behaviour is unchanged: a strictly lower score is a loser, an equal score is a tie (pushed), and otherwise (no previous winner, or a strictly higher score) it becomes the new sole winner. Also apply rustfmt to the makedb code added in earlier commits. Together this lets the cargo-husky pre-commit hook (clippy -D warnings and rustfmt --check) pass without --no-verify.
Address two review comments on PR #290: - makedb no longer calls check_reference_existence, which for BWA inspects pre-existing index files beside the input reference. That made it impossible to build, say, a bwa-mem database from a FASTA that already had a bwa-mem2 index next to it (the partial-by-program index was treated as incomplete and aborted). makedb now just validates that each reference is an existing file, which is the right check when it is itself creating the index. - Pass the selected thread count to 'strobealign --create-index' via -t, so 'coverm makedb -p strobealign -t N' actually indexes with N threads, matching the minimap2 path and the --threads documentation.
Add a new coverage method, mean_zero_coverage_length, that reports the average length (in base pairs) of the regions of each contig or genome that have zero coverage. Within a covered contig, each maximal run of consecutive zero-coverage bases is one region. An entirely-unobserved contig (or, in genome mode, an entry with no reads mapped at all) is treated as a single zero-coverage region spanning its whole length, which keeps results identical regardless of whether the streaming (print_zero_coverage) or non-streaming (calculate_coverage) code path is taken. A contig/genome with no zero-coverage positions is reported as 0. Like the other structural estimators (length, count, reads_per_base, anir), this method does not compute a covered fraction and so cannot be combined with --min-covered-fraction > 0. Includes unit tests for the estimator logic, an end-to-end contig test, a cmdline integration test, and README/help documentation.
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: dbff109387
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
Comment on lines
+892
to
+898
| if *num_mapped_reads == 0 { | ||
| // No reads mapped to this entry at all, so the whole entry is | ||
| // a single zero-coverage region spanning its full length. | ||
| // Reporting the total length here (rather than the per-contig | ||
| // average) keeps results identical to the print_zero_coverage | ||
| // path, which only knows the entry's total length. | ||
| return unobserved_bases as f32; |
There was a problem hiding this comment.
In genome mode with --no-zeros, src/genome.rs decides whether to print an entry by checking whether any calculated coverage is > 0; returning the entry length here when num_mapped_reads == 0 makes completely unmapped genomes look nonzero for -m mean_zero_coverage_length, so coverm genome --no-zeros -m mean_zero_coverage_length will still emit genomes with zero read coverage. This needs a separate path for reporting the length without using that positive value as the print/no-print signal.
Useful? React with 👍 / 👎.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Adds a new coverage method,
mean_zero_coverage_length, that reports the average length (in base pairs) of the regions of each contig or genome that have zero coverage.This builds on top of #290 (the
coverm makedbbranch), so its diff is included here; this PR's own changes are the new coverage method.Behaviour
1000).0.Cross-path consistency
Genome coverage is computed via two code paths — streaming (
--separator) and non-streaming (--genome-fasta-files). For an entirely-unobserved multi-contig genome, the streaming path only knows the genome's total length (viaprint_zero_coverage). To keep results identical regardless of input format, an entry with zero mapped reads is defined as a single region of its full length. The estimator'snum_mapped_readsis available in bothcalculate_coverageandprint_zero_coverage, so both paths agree, and it also distinguishes an absent entry (→ length) from a fully-covered one (→ 0).Like the other structural estimators (
length,count,reads_per_base,anir), this method does not compute a covered fraction and so cannot be combined with--min-covered-fraction > 0(it errors with a clear message).Changes
src/mosdepth_genome_coverage_estimators.rs: newMeanZeroCoverageLengthEstimatorvariant, all trait arms, and 4 unit tests.src/bin/coverm.rs: method dispatch and--min-covered-fractionincompatibility guard.src/cli.rs: added to bothgenome/contigvalue-parser lists and both help tables.src/contig.rs: end-to-end test (coveredseq1→ 68.25, uncoveredseq2→ 1000).tests/test_cmdline.rs: cmdline integration test.README.md: row in the calculation-methods table.Testing
cargo clippy -- -D warningsandcargo fmt -- --checkare clean (enforced by the pre-commit hook).test_streaming_bam_file) fails in this environment only becausesamtools/bwaare not installed.The generated
docs/coverm-*.htmlman pages were not regenerated, as that requires theman/pandoctooling used at release time; the source README and CLI help text are updated.https://claude.ai/code/session_0126FFDNupzAHRS5m74JA9Zp
Generated by Claude Code