Conversation
Add a --gff option to `coverm contig` and `coverm genome` that accepts one or more GFF3 files. When provided, CoverM reports coverage statistics for each feature (gene) defined in the GFF files rather than per-contig or per-genome. Coverage for each gene is computed over its region on the parent contig (GFF seqid), reusing the existing coverage estimators, so all coverage methods (mean, trimmed_mean, covered_fraction, variance, length, count, reads_per_base, rpkm, tpm, anir, ...) are supported. Reads are attributed to genes by reference overlap, giving per-gene read counts and identities. - New src/gene.rs: GFF3 parser and gene_coverage() function. - Wire --gff into the contig and genome subcommands and full help. - Output uses a "Gene" header and one row per feature. - Add test GFF files (tests/data/genes/) plus unit and assert_cli tests. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013JBSD2WGRUBCErFRkjjeCB
![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: 8f99a1f578
ℹ️ 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".
| for estimator in coverage_estimators.iter_mut() { | ||
| estimator.setup(); | ||
| estimator.add_contig( | ||
| &gene_ups_and_downs, |
There was a problem hiding this comment.
When --gff is used with the default --contig-end-exclusion 75, each gene is converted into its own gene_ups_and_downs slice and passed to CoverageEstimator::add_contig, whose mean/trimmed/variance estimators drop the first/last contig_end_exclusion bases of the slice. That makes any covered gene ≤150 bp report 0 coverage and trims 150 bp from every longer gene in the README's default --gff example; the exclusion should stay tied to contig ends or be forced to 0 in gene mode.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Fixed in e65462b. Contig-end exclusion is now disabled for the coverage estimators when --gff is used (new CoverageEstimator::disable_contig_end_exclusion()), since gene boundaries are not contig ends. Added a unit test asserting a whole-contig gene reports the same value with --contig-end-exclusion 75 as with 0.
| if is_primary { | ||
| gene.num_mapped_reads += 1; | ||
| if aligned_len > 0 { | ||
| gene.sum_identity += identity; |
There was a problem hiding this comment.
For reads that span multiple genes, the code computes identity/NM over the entire alignment and then adds that same value to every overlapped gene. With --gff -m anir, a mismatch or indel in one gene will lower the ANIr reported for an adjacent gene that the read also overlaps, even when the bases inside that gene match perfectly; compute identity/edit distance over the read's overlap with each gene instead.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Considered, but keeping whole-read attribution deliberately. The per-read edit distance comes from the NM tag, which is a whole-read statistic and cannot be split positionally between genes without parsing MD tags / the reference. Prorating NM by overlap length reduces mathematically back to the whole-read identity, so it isn't genuinely more accurate. The alternative (only counting reads fully contained in a gene) avoids cross-gene contamination but silently drops boundary-spanning reads from ANIr. Since genes are typically much longer than reads, we've opted to keep the simple whole-read attribution and documented the limitation in a code comment.
Addresses PR review: in gene mode each gene is a sub-region of a contig, so the default --contig-end-exclusion (75) was trimming the ends of every gene and zeroing out any covered gene shorter than 150 bp. Contig-end exclusion is now disabled for the coverage estimators when --gff is used, since gene boundaries are not contig ends. Add CoverageEstimator::disable_contig_end_exclusion() and a unit test that a whole-contig gene reports the same coverage with contig-end-exclusion 75 as with 0. Per-gene edit distance / ANIr continue to use whole-read NM statistics (NM cannot be split positionally between genes), documented in a comment. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013JBSD2WGRUBCErFRkjjeCB
1 participant
Summary
Adds a
--gffoption tocoverm contigandcoverm genome. When one or more GFF3 files are supplied, CoverM reports coverage statistics per feature (gene) rather than per-contig (contig mode) or per-genome (genome mode).How it works
src/gene.rscontains a small GFF3 parser (read_gff_files) andgene_coverage(), which mirrorscontig_coveragebut computes coverage over each gene's region on its parent contig.CoverageEstimators, so all coverage methods work (mean, trimmed_mean, covered_fraction, covered_bases, variance, length, count, reads_per_base, rpkm, tpm, anir, ...).IDattribute (falling back toName, thenseqid_start_end). The output uses aGeneheader with one row per feature.##FASTAsection of a GFF3 file is ignored.Tests
src/gene.rsfor the GFF parser and gene coverage (a gene spanning a whole contig reproduces the whole-contig coverage; per-gene read counts).assert_cliintegration tests intests/test_cmdline.rsfor bothcontig --gffandgenome --gff(mean, length and count methods).tests/data/genes/.All new tests pass;
cargo clippy --fixandcargo fmtapplied.🤖 Generated with Claude Code