Add mutation spectra analysis#37
Merged
Merged
Conversation
Introduce derip2-spectra, a new entry point that builds SigProfiler SBS-96 (pyrimidine-collapsed) and SBS-192 (strand-resolved) trinucleotide mutation spectra by calling every substitution of each aligned sequence against the deRIP'd ancestral consensus. - src/derip2/spectra/: channel ordering + pyrimidine folding + nearest non-gap context (channels.py), SigProfiler-compliant matrix I/O (matrix_io.py) - src/derip2/stats/mutation_spectra.py: SpectraResult + vectorised compute_spectra, with a multi-hit-column homoplasy proxy - src/derip2/plotting/spectra.py: native SBS-96/192, strand-asymmetry (binomial screen) and homoplasy figures, reusing the strand-bias palette - DeRIP.calculate_spectra/write_spectra_matrix/plot_spectra - app_spectra.py CLI + derip2-spectra script entry point - tests: channels, spectra (synthetic + golden), plots, CLI (32 tests) Pure-Python, no new dependencies. Full suite green (234 passed). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Merging this PR will not alter performance
Performance Changes
Comparing 
Footnotes
|
Layer the rigorous IQ-TREE path on top of the tree-free baseline. Instead of comparing every sequence to one reference, it reconstructs ancestral sequences on a rooted tree and calls each substitution per parent->child branch, reading context from the parent. This counts recurrent (homoplasic) deamination correctly: on mintest, a site the baseline scores as 4 hits is 2 independent branch events. - src/derip2/spectra/tree_asr.py: first subprocess wrapper in the codebase. Runs IQ-TREE --ancestral (iqtree3/2 auto-detected), parses .state, loads the tree with ete4, roots it (outgroup/midpoint/none) by orienting edges from an existing node so every node keeps a sequence. reconstruct(), assign_clades(), orientation_flip_fraction() for root sensitivity. - src/derip2/spectra/call_mutations.py: per-branch event stream -> shared assembler; true per-branch homoplasy; parent/child event provenance; min_prob posterior filtering. - src/derip2/stats/mutation_spectra.py: extract assemble_matrices() + IDX96/IDX192 code tables shared by both paths; generalise SpectraResult with method + parent/child event names. - src/derip2/spectra/qc.py: per-column gap/ambiguity profile + reports. - app_spectra.py: --method phylo, --tree, --iqtree-model, --threads, --rooting, --outgroup, --min-prob, --root-sensitivity, --partition-by clade; run manifest. - environment.yml/pyproject: iqtree + ete4 (optional 'spectra' group). - tests: hand-built trees (polarity, homoplasy, inheritance, min_prob, clades), committed .state/.treefile fixtures, QC, IQ-TREE runs guarded by shutil.which. 253 passed. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…o spectra_profiling
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #37 +/- ##
==========================================
+ Coverage 86.99% 87.50% +0.51%
==========================================
Files 12 22 +10
Lines 1853 3034 +1181
Branches 306 488 +182
==========================================
+ Hits 1612 2655 +1043
- Misses 151 252 +101
- Partials 90 127 +37 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…tone 3) - docs/tutorials/mutation-spectra.md: full tutorial with real example images from the Sahana transposon, added to mkdocs nav. Covers baseline vs phylogenetic, supplying a precalculated phylogeny (--tree with iqtree3), and the recommended RIP-masked-topology workflow (deRIP2 --mask -> iqtree3 -> --tree so the masked topology carries ancestral states recomputed for the UNMASKED sequences). README + CLI docs updated with the same workflows. - run_iqtree now always passes -m model, so --iqtree-model applies to fixed-tree (-te) runs — required for the masked-topology workflow. - tree_asr: map IQ-TREE-sanitised tree leaf names ([^A-Za-z0-9._-] -> _) back to alignment ids, so realistic FASTA headers (e.g. 'scf:1-9(+)') work. New _sanitize_name + tests. - Benchmarks for compute_spectra and compute_spectra_from_tree. - Reconcile pyproject: ete4 is a core dep; drop the redundant spectra extra (IQ-TREE binary still installed separately). Baseline-vs-phylo contrast on 60 Sahana copies: 38,727 vs 7,560 events, 1,163 sites hit >=3x independently. 256 tests pass. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ilestone 4) Feature: per-group spectra - --groups FILE maps sequence names to group labels (e.g. species) and reports one spectrum per group. Works for baseline (per-row samples) and phylo (a branch is attributed to a group only when its whole descendant clade shares it; spanning branches -> 'mixed'). Name matching tolerates IQ-TREE sanitisation, so the same label file works for MSA ids and tree tips. assign_groups() in tree_asr, _load_group_lookup() in app_spectra, samples= override on DeRIP.calculate_spectra. Tests for both methods. Plot fixes - Strand asymmetry: bars now encode strand by colour (blue=coding, orange= template) so the legend matches; added a footnote explaining '*' = binomial p<0.05 vs 50:50. - Homoplasy: translucent stacked markers (overlaps stay visible) + a substitution-class legend; y-axis names the recurrence unit per method. - homoplasy_table now counts recurrence over the actual per-event ref/alt (fixes a phylo KeyError where the root base equalled a derived base). Docs (docs/tutorials/mutation-spectra.md) - Real COSMIC decomposition of the full Sahana spectrum via SigProfilerAssignment (cosine 0.80) with a prominent caveat that human-cancer signatures do not apply directly to fungi and a fungal reference library is the right path. - Per-clade and per-group example figures; masked-topology spectrum figure; baseline/phylo homoplasy figures. - Full-Sahana (n=396) baseline and phylogenetic spectra (n in titles); honest note that ML ASR's majority bias inflates T>C on heavily-RIP'd families, so the consensus baseline is more robust for direction. - CLI docs + README note --groups; masked workflow gains -st DNA (heavily masked alignments defeat IQ-TREE's sequence-type auto-detection). 260 tests pass. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Statistical comparison of spectra (new stats/spectra_compare.py, no SciPy): - cosine_similarity (effect size) and a dependency-free chi-squared test of homogeneity (regularised incomplete gamma; validated to match SciPy exactly), with per-channel standardised residuals and Cramer's V. - compare_spectra (two vectors -> cosine, chi2, top differentiating channels) and pairwise_compare (all group pairs, Bonferroni). Works on group columns or two precalculated matrices. Tutorial adds worked examples for both, with the caveat to read the p-value beside the effect size (huge n over-powers chi2). Strand-asymmetry plot: only star a class when both strands carry >= min_count events (default 10) so a near-empty strand is not flagged; asterisk sits on the bar top (hugs small bars at the axis) not floating; footnote states the rule. Homoplasy plot: legend moved to a horizontal row below the axes so it no longer overlaps the densely-packed data points. Baseline: full worked example in the tutorial explaining the deRIP-consensus ancestor - the default internal reconstruction and how to compute/pass it explicitly via `derip2` + `--ancestor` (CLI and Python). 269 tests pass; docs build clean. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…+tests derip2-spectra baseline can now reuse an ancestor already in the input MSA: - new --reference-tag (default deRIPseq, exact id) detects a consensus row, uses it as the ancestor, and excludes it from the counted sequences - precedence: --ancestor file > in-MSA tag > recomputed deRIP consensus - --ancestor validated for length == alignment width up front; read via new ao.loadFirstSequence so single-sequence FASTAs work Input validation: - ao.validate_no_degenerate gates input to A/C/G/T/- (lowercase accepted); N and IUPAC codes raise a clear, located error - ao.uppercase_alignment normalises soft-masking so lowercase input no longer silently changes the deRIP consensus (RIP detection is upper-case only) strandbias plot: named ARM_LABEL_HEADROOM (1.22) replaces the inline 1.06 so the forward/reverse strand arm labels clear the tallest bars. Tests: new tests/conftest.py (data-path fixtures, make_reconstruction factory, have_iqtree); offline coverage for tree_asr guards, app_spectra new features and edge cases, spectra-plot branches, strandbias headroom. 307 tests pass; plotting/spectra 100%, tree_asr/app_spectra 97%. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…igures
- README + CLI docs + mutation-spectra tutorial: --reference-tag / in-MSA
deRIPseq reuse, --ancestor length validation, and the A/C/G/T/- input
requirement (IUPAC rejected). Corrected the worked example so the reused
row id matches derip2's --prefix (default deRIPseq).
- Regenerated strand_bias_{derip,logo,noise}.png so the forward/reverse arm
labels clear the tallest bars (ARM_LABEL_HEADROOM change). Content identical,
same dimensions, only added vertical headroom.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…o spectra_profiling
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
1 participant
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.
Introduce derip2-spectra, a new entry point that builds SigProfiler SBS-96 (pyrimidine-collapsed) and SBS-192 (strand-resolved) trinucleotide mutation spectra by calling every substitution of each aligned sequence against the deRIP'd ancestral consensus.
Pure-Python, no new dependencies. Full suite green (234 passed).