Integrate clustering track and parser loader reliability updates - #121
Open
ngamarra wants to merge 322 commits into
Open
Integrate clustering track and parser loader reliability updates#121ngamarra wants to merge 322 commits into
ngamarra wants to merge 322 commits into
Conversation
…nsemble control-calibration
Phase 0
- plotting: prepare_dmr_site_data / prepare_dmr_segment_data /
prepare_dmr_multi_summary_data (backend-agnostic, flexible modkit-column detection)
+ plot_dmr_site/segment/multi_summary_matplotlib renderers.
- dmr_multi_sample.ipynb: demonstrate the new DMR helpers on the demo BAMs.
- region_contrasts: remove dead _estimate_beta_binomial_prior /
_beta_binomial_two_sided_p_value (+ their tests); import numpy for the LRT
annotations; split semicolon statements.
- __init__: surface ContrastSpec/SampleSpec/CohortSpec/RegionContrastResult/
RegionDiscoveryResult (additive, back-compat).
- tests: new test_dmr_plotting, test_dmr_integration (parse->pileup->dmr, skips
without modkit/reference), test_plotting_matplotlib, test_utils.
- lint: workflows.py logger below imports; cluster.py zip strict= (repo-wide ruff clean).
S1 — Q2 ensemble control-calibration
- score_regions: normalization_mode ('none'|'global') + global_normalization_factors
kwarg; _apply_global_normalization subtracts the per-(sample,motif) global_offset,
clamps, and re-derives modified_count (valid_count unchanged; NaN/absent = identity);
applied on the ensemble path before pooling and the LRT. Guard: normalization is
ensemble-only (fail loud for single_read/cluster_occupancy).
- background_adjusted mode: _pool_background_rate (coverage-weighted b) +
_apply_background_correction (clamp((p-b)/(1-b)); b>=1 -> 0); reuses the existing
replicate-level beta-binomial LRT + BH unchanged; adds background_fraction/
background_present; global normalization runs first.
- plotting: prepare_region_contrast_correction_overlay_data +
plot_region_contrast_correction_overlay_matplotlib (region-level raw-vs-corrected).
- tests: test_s1_control_calibration (engine math, ordering, precedence, guard, overlay).
Reviewed via two adversarial workflows; all confirmed findings fixed.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…background subtraction) New region_contrast_control_calibration.ipynb demonstrating the S1 API on the CTCF demo BAMs: build pileups + SampleSpecs, compute global normalization factors and run score_regions(normalization_mode="global") (runnable), plus a background_adjusted + raw-vs-corrected overlay template (point background=[...] at a negative control). The global-normalization chain is validated end-to-end on the demo data. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ignal calling (Q3/G1) New dimelo/background.py: fit a beta-binomial background null from negative-control per-read (k, n) counts, then call each target read as true signal at a controlled FDR. - fit_beta_binomial: MLE (alpha, beta) via (mu, s) with multi-start on concentration (generalizes the polii read-level prototype). - build_control_null: per-site null where control coverage >= min_control_reads, else a pooled fallback null. (Context-conditioned pooled fallback is a planned refinement.) - call_true_signal_reads: one-sided upper-tail P(K>=k | n, null) per target read -> BH q_value -> is_true_signal at FDR alpha; valid_count==0 reads are untestable. Consumes the region_contrasts.build_single_read_mod_fraction_evidence_table schema. - summarize_site_occupancy: per-site control-calibrated occupancy rate (feeds Q6). - Exported dimelo.background; tests in test_background.py. First vertical slice of S2; occupancy_posterior (soft mixture), background-removed pileup, and plotting/tutorial to follow. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
… pileup, and plots - background.call_true_signal_reads now also returns occupancy_posterior: a soft per-read P(signal | k, n) from a 2-component beta-binomial mixture (background fixed to each read's assigned control null; signal component + mixing weight fit by EM, signal mean constrained above background). valid_count==0 reads get NaN. - background.background_removed_pileup: per-site pileup from true-signal reads (weighting="hard") or posterior-weighted (weighting="posterior") -> the background-removed signal track (links to the Q2 ensemble correction). - plotting: prepare_true_signal_read_data + plot_true_signal_read_raster_matplotlib (per-read raster colored by is_true_signal or posterior, background-rate line); prepare_occupancy_rate_track_data + plot_occupancy_rate_track_matplotlib; prepare_background_removed_pileup_overlay_data + plot_background_removed_pileup_overlay_matplotlib (raw vs removed). - tests extended in test_background.py. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ol-calibration notebook New section demonstrating dimelo.background: call_true_signal_reads (per-read FDR + occupancy_posterior), summarize_site_occupancy, background_removed_pileup, and the three S2 plots. Runs on a small inline evidence table (validated end-to-end); notes how to build real evidence via build_single_read_mod_fraction_evidence_table. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…test-hardening Bayesian methods (user request): summarize_site_occupancy now reports a Beta-Binomial conjugate posterior for occupancy theta (default Jeffreys Beta(0.5, 0.5) prior) alongside the raw point estimate: - count_mode="posterior" (default when occupancy_posterior is present) uses the soft per-read posteriors as expected bound counts (propagates read-level uncertainty); count_mode="hard" uses the binary is_true_signal calls. - returns occupancy_posterior_mean + occupancy_ci_low/high (credible_mass, default 0.95). - plot_occupancy_rate_track_matplotlib draws the posterior mean with credible-interval error bars when those columns are present. S2 adversarial review (0 code bugs) test-hardening: - background_removed_pileup posterior branch now unit-tested with hand-built reads where is_true_signal and occupancy_posterior differ, pinning exact hard vs posterior values. - added a one-sidedness unit test for _upper_tail_pvalue (k below the null -> p~1). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
binding_strength: per-site binding strength = single-molecule occupancy (bound/total from Q3 is_true_signal), reported with the Bayesian Beta-posterior mean + credible interval, and control-calibrated by comparing target vs control occupancy (control reads scored against their own null = false-positive baseline): occupancy_excess + a two-proportion test of bound counts (BH across sites) + significance + rank. Deliberately excludes SM methylation magnitude / footprint from the strength metric (those are Q5). binding_strength_consistency_qc: separate validation view — per-site mean read_mod_fraction of true-signal vs background reads (signal_gap), confirming bound reads carry more signal without feeding it into the strength metric. plotting: prepare_binding_strength_data + plot_binding_strength_matplotlib (ranked occupancy with credible-interval error bars, control occupancy overlaid, significant sites highlighted). Tests in test_binding_strength.py. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…criminating tests - summarize_site_occupancy count_mode="auto" now resolves on occupancy_posterior column PRESENCE, not on whether any value is finite, so a site's Beta-posterior occupancy no longer depends on the other sites in the batch (a zero-coverage-only site was flipping between posterior and hard modes by batch composition). - test_binding_strength: pin occupancy to the hard n_bound/n_reads (distinct from the soft posterior mean), replace the tautological CI-containment check with a signal/noise separation (ci_low > 0.5), and add a control-calibration test with a "leaky" site whose control carries spurious high signal — high raw occupancy but ~0 excess and NOT significant — which fails if the control-subtraction arm is dropped. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ion notebook Demonstrates background.binding_strength (control-calibrated occupancy + Beta-posterior CI + excess + significance), binding_strength_consistency_qc, and plot_binding_strength on the inline single-molecule evidence from section 3. Validated end-to-end. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
New dimelo/footprint.py: a lightweight, hand-rolled 2-state Bernoulli HMM (no new dependency; hmmlearn absent on the cluster) for single-molecule footprint detection. - BernoulliHMM: log-space forward-backward, Viterbi, and Baum-Welch fit over pooled read sequences (nan = missing/uninformative); states re-ordered so state 0 is protected. - call_footprints: per read, Viterbi-decode the accessible/protected path and call a footprint = a protected run near the motif anchor within [min_width, max_width]; reports footprint_present, protected coords/width, and a posterior footprint_score. - aggregate_footprint_profile: per-position methylation density + mean protected posterior (the footprint shows as a methylation dip / posterior rise at the anchor). - plotting: prepare_footprint_profile_data + plot_footprint_profile_matplotlib. Footprint is a Q5 chromatin-state feature and a confidence input to the Q3 bound/unbound call -- deliberately NOT part of the Q6 binding-strength metric. Tests in test_footprint.py. Follow-ups: wire footprint into cluster state annotation; per-read state raster overlay; tutorial. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Replace the hand-rolled Bernoulli HMM with hmmlearn's CategoricalHMM (two symbols: 0 = unmethylated, 1 = methylated), per request. BernoulliHMM keeps the same public API (fit / viterbi / posterior / p_emit) so call_footprints, aggregate_footprint_profile, the plotting, and the tests are unchanged. nan observations are dropped for fit/decode and the decoded path is carried forward to stay aligned to the input columns. fit uses an informed initialization (init_params="", emissions spread low->high) so Baum-Welch separates a protected from an accessible state rather than collapsing from a random init. hmmlearn added to install_requires (footprint is imported by dimelo/__init__). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…anchors + harden tests - call_footprints now requires observed methylation data near the anchor and that a candidate protected run contains >=1 observed position; otherwise Viterbi's prior fallback over all-nan / anchor-missing reads produced spurious full-width footprints (S4 review finding). Adds an n_observed column per read. - tests: add a smoothing/Viterbi test on ambiguous data (sticky HMM keeps an isolated methylation flip inside one protected run), an all-nan / anchor-unobserved guard test, and max_width, anchor_tolerance-boundary, and equal-width tie-break tests. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
New dimelo/joint_occupancy.py: for user-supplied anchor pairs, quantify whether two anchors are co-occupied on the same single molecule more than expected under independence. - joint_occupancy: for reads spanning BOTH anchors, per-read joint state (both/a_only/b_only/neither) from the Q3 is_true_signal call at each anchor; per pair a 2x2 contingency -> Fisher exact + odds ratio + log2 observed/expected co-occupancy, BH across pairs, plus a Spearman correlation of the soft occupancy_posterior between anchors. Reports pair_distance_bp for the Q4 trans-contact confound. - summarize_joint_states (stacked-bar counts/fractions) and anchor_pairs_from_sequence. - plotting: prepare_joint_state_distribution_data / prepare_joint_occupancy_distance_data + plot_joint_state_distribution_matplotlib / plot_joint_occupancy_distance_matplotlib. - tests in test_joint_occupancy.py (coupled vs independent pairs, spanning-read filter, min-spanning-reads gating, plotting). Generalizes the polii promoter-enhancer joint-occupancy prototype; reuses Q3 per-read calls. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…-read re-lookup) Iterate the merged spanning frame once instead of re-looking-up occupancy_posterior by read_id per row (O(n^2) and ambiguous under duplicate read_ids). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…th DiMeLo signal New dimelo/tracks.py (coverage-first modality integration): - bin_regions: tile regions into equal-width bins. - import_bigwig_signal: per-bin bigWig coverage aligned to regions (via pyBigWig .stats). - import_bedgraph_signal: per-bin bedGraph coverage (overlap-length-weighted mean). - correlate_binned_signals: Pearson + Spearman of DiMeLo vs external binned signal on matched (region_id, bin), NaN-safe with constant/insufficient-data guards. - plotting: prepare_track_correlation_data + plot_track_correlation_matplotlib (scatter annotated with the correlation). tracks exported from dimelo. Covers MNase / ATAC / CUT&Tag / ChIP coverage comparison. Hi-C (cooler) import and the Hi-C-vs-joint-occupancy view are the next slice (heavier dependency). Tests in test_tracks.py (bigWig/bedGraph fixtures, binning, overlap weighting, correlation). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Extend the single-molecule notebook with: - Q5: footprint.call_footprints + aggregate_footprint_profile on a synthetic footprinted window, with the aggregate footprint plot. - Q7: joint_occupancy on per-read Q3 calls at two anchors (2x2 Fisher + log2 O/E + posterior correlation) with the joint-state distribution plot. Both validated end-to-end. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…coverage tests - import_bigwig_signal now returns NaN for zero-width bins (bins >= region length collide integer edges) instead of crashing pyBigWig.stats with 'Invalid interval bounds!' (matches the bedGraph path). - tests: bedGraph test now uses an ASYMMETRIC straddling bin (20bp value-1 + 10bp value-5 -> 70/30) so overlap-length weighting is distinguishable from a plain mean; add a non-divisible bin_regions remainder test; add a zero-width-bin regression test; strengthen the correlation-plot test (labels, scatter count, annotation text) and cover the empty-paired and NaN-annotation branches. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…finite reads Spearman posterior correlation now masks to reads with a finite occupancy_posterior at BOTH anchors before correlating (default nan_policy='propagate' otherwise let a single low-coverage/NaN-posterior read null the coupling for a well-sampled pair); the n>=2 and std>0 guards now apply to the filtered reads. Tests: add an asymmetric coupled pair (n_a_only != n_b_only) pinning all four contingency cells + per-read a_only/b_only counts (catches an a/b transpose), and a NaN-posterior test. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
… (Q8/Q7 bridge) - import_hic_contacts: per anchor-pair Hi-C contact frequency from a .cool file (sum of the inter-anchor sub-matrix via cooler; NaN on empty/absent), pairs given as region strings. - correlate_hic_vs_joint_occupancy: joins Hi-C contact with a joint_occupancy pair-summary on pair_id and correlates contact vs log2 observed/expected co-occupancy (Pearson + Spearman) -- the architecture bridge and a Q4 trans-contact flag (high contact + high co-occupancy at large 1D distance). - cooler added as an optional 'hic' extra (imported lazily). Tests build a small .cool fixture (skipped if cooler is absent). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…earch spike) New dimelo/trans_contact.py: separate in-site (cis) from contact-driven (trans) methylation using two independent pieces of evidence, since the negative control (Q3) does not catch target-driven trans methylation. - estimate_cis_reach: empirical enzyme 1D reach from a methylation-density decay profile around isolated bound sites (half-decay offset; user-overridable). - classify_cis_trans: per position, cis if within cis_reach_bp of a Q3-bound site; else trans if its Hi-C contact to a bound site clears a threshold (3D-proximal, 1D-distal); else distal_background. - apply_trans_correction / summarize_trans_correction: in-site-only corrected signal + a separate trans track (feeds Q10 architecture; down-weights Hi-C-explained Q7 pairs). Exploratory (spike): validate that trans-flagging reduces spurious Q7 joint occupancy at high-Hi-C distal pairs before relying on it. Depends on S6b Hi-C + Q3 calls + Q5 footprint. Tests in test_trans_contact.py. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ach branch coverage - assert nearest_bound_distance numerically and add a position exactly at cis_reach_bp (inclusive <=) so a reach off-by-one is caught. - cover estimate_cis_reach's no-crossing (len(below)==0) branch with an explicit background_level, and a non-default decay_fraction; rename the flat-profile test to reflect the peak==background branch it exercises. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ule state composition New dimelo/integration.py (the capstone synthesis of Q2-Q8): - build_integrated_site_table: outer-join the per-site metrics from occupancy (S2), binding strength (S3), footprint summary (S4), joint-occupancy summary (S5), external-track/Hi-C signal (S6/S6b) and trans-correction (S7) into one denormalized per-site table; columns colliding across sources are prefixed with the source name, unique columns kept as-is. - single_molecule_state_composition: cross-tabulate per-read annotations (bound/unbound, footprint present, co-occupancy state) into the fraction of single molecules in each combinatorial chromatin state per site -- the cell-population heterogeneity view. - plotting: prepare_state_composition_data + plot_state_composition_matplotlib (stacked bar). Completes the S1-S8 roadmap (Q2-Q8 + Q4 trans + Q10 synthesis). Tests in test_integration.py. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ollision guard - single_molecule_state_composition: build combined_state via elementwise str() (map) instead of astype(str), so a NaN/None in a state column (common when joining per-read outputs) becomes a literal 'nan'/'None' label rather than crashing '|'.join under pandas 3.x (which preserves NaN through astype(str)). Env is pandas 3.0.3. - build_integrated_site_table: after on-collision prefixing, fail loud if any output column name still collides across sources (e.g. a source column literally named like another source's prefixed output), instead of letting the outer merge emit _x/_y. - tests: NaN-state composition + residual-collision cases. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…iter `modkit extract full` emits a row for every read position, including soft-clipped / inserted bases that are not aligned to the reference (ref_position == -1). read_by_base_txt_to_hdf5 updated read_start/read_end from every row, so a -1 row drove start_candidate negative: read_start collapsed far below the true leftmost aligned position and each read's stored per-base vector (sized read_end - read_start) ballooned to the genomic-coordinate scale (observed: median read span ~49 Mb, dense vectors up to 180 MB/read). This made downstream per-read analyses intractably slow and their fractions meaningless. Skip rows with pos_in_genome < 0 at the top of the parse loop. Read spans return to sane values (median ~1.6 kb, none >1 Mb) and h5s shrink from GB-scale to MB-scale. Fixes the E-P joint-occupancy extraction and the shared S7 clustering extraction path. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…tion (Q8 connectors) Public-reference-data connector slice (fetch/cache separated from offline compare): - jaspar.py: fetch_jaspar_matrix (cached JASPAR REST client; requests lazily imported) + jaspar_matrix_from_json (offline PFM parser). - motifs.py: counts_to_probability, information_content (logo bits), scan_sequences (log-odds PSSM best hit), observed_pfm, compare_motifs (observed vs annotated), and extract_site_sequences (pysam) -- "do DiMeLo sites carry the annotated motif / match the logo?". Generalizes the notebook CTCF PWM prototype. - plotting: prepare_sequence_logo_data + plot_sequence_logo_matplotlib (information-content logo). - cluster_reference.py: characterize what a cluster's sites are associated with via reference site-groups -- site_overlaps_group + cluster_feature_enrichment (per-cluster Fisher vs the rest + BH), annotate_clusters_by_site_groups (peak overlap per group), annotate_clusters_by_motif, and summarize_cluster_associations (per-cluster associated groups). Answers "what is this cluster's sites associated with using these site groups?". Modules exported; tests in test_motifs.py / test_cluster_reference.py (offline; JASPAR fetch uses a cache fixture, no live network). See the reference-connectors vision for the roadmap (ENCODE, gget, scRNA next). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…n overlap tests - jaspar.fetch_jaspar_matrix: validate the payload (parse) BEFORE writing the cache so a malformed HTTP-200 response never poisons it, and treat a corrupt/unparseable cache file as a miss (refetch) instead of hard-failing. - cluster_reference.cluster_feature_enrichment: drop sites with a NaN cluster label before computing totals, so unassigned/noise sites no longer leak into every cluster's 2x2 background row (sum of per-cluster n_sites now equals the tested total). - tests: half-open overlap boundary cases (touching = no overlap, 1 bp = overlap) and an unlabeled-site exclusion test. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
New dimelo/encode.py: search the ENCODE portal and fetch files, routing into the existing comparison machinery. Network fetch is separated from offline parsing/loading (testable without network); downloads are cached. - parse_encode_search / parse_encode_files: offline parsers (experiments table; downloadable files with absolute URLs + format/output_type/assembly filters). - search_encode_experiments / fetch_encode_experiment_files: network queries (requests lazy). - download_encode_file: streamed, cached download (temp-file then rename; cache hit = no network). - load_encode_peaks: BED/narrowPeak (.gz ok) -> chromosome/start/end for overlap. Feeds ChIP/TF peaks into cluster_reference/regulatory_enrichment (validate binding sites, cluster association), histone-mark bigWigs into tracks correlation, and Hi-C .cool into tracks.import_hic_contacts. Tests offline (fixture JSON + a cached-download path) incl. an end-to-end ENCODE-peaks -> cluster association wiring test. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Coerce start/end to numeric and drop non-data rows so a BED/narrowPeak with a UCSC 'track'/'browser' header (not '#'-commented) is parsed cleanly instead of crashing astype(int) with IntCastingNaNError. + regression test. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
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.
Summary:
Validation:
Notes: