Measure modification against charging with odds ratios#36
Open
jayhesselberth wants to merge 5 commits into
Open
Measure modification against charging with odds ratios#36jayhesselberth wants to merge 5 commits into
jayhesselberth wants to merge 5 commits into
Conversation
Add compute_charging_odds_ratios(), which uses individual reads as the unit of observation to test whether a read being modified at a site is associated with that same read being charged. Because both variables are measured on the same read, this is a within-tRNA comparison that differences in abundance or charging between tRNAs or samples cannot confound. Sites can come from either of two signals. By default every position in the calls is tested. Passing a site list from the new call_bcerror_sites() instead restricts testing to positions called from base-calling error, which is useful when the modification-caller channels are not trusted for a dataset: error is model-free, so it gives an independent route to the same question. In that mode the calls should be a per-read mismatch table, so that per-read status and site selection come from the same signal. Reads with no call at a site are excluded from that site's table rather than counted as unmodified, so total_obs varies between sites. read_charging_calls() supplies the per-read charging status these need; the existing read_charging() returns per-tRNA aggregates and cannot. compute_odds_ratios() also gains the sites argument, so the same site list can restrict the pairwise analysis. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The R implementation pivoted each reference's calls into a dense read by position matrix and took a crossprod over it. That is the wrong shape for the problem: the only state a 2x2 table needs is four counters per (reference, position), so the matrix was pure overhead and grew with read depth. On a 100M-row call table it needed tens of gigabytes. Accumulate the counts in a single pass in C++ instead. Memory now scales with the number of sites rather than the number of reads: the same table across 46,000 sites takes about 100 seconds and 10 GB, nearly all of it holding the parsed input rather than the computation. Fisher's exact test matches pairwise_fisher_exact, so the two odds ratio functions report the same statistic. Add min_margin and max_p to drop sites before testing rather than after. A genome-wide site list is mostly sites that were never testable, and each one still consumes FDR budget; max_p drops those whose margins put the intended threshold out of reach whatever the reads did. Both filter on the margins, which are ancillary, so the null distribution of the surviving p-values is unchanged. Filtering on the observed odds ratio would enrich for small p-values and break the correction applied afterwards, so that is deliberately not offered. dedupe = FALSE skips the per-read collapse for call tables that already carry one row per read and position. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Measuring against the Watt zebrafish samples put numbers on the caveat. 82% of terminal-A sites come back significant with a median log odds ratio of 1.9, decaying monotonically with distance from that residue to 4.6% beyond fifty nucleotides. The amino acid sits on the terminal A and the charging model reads the signal there, so the base caller is reporting the same thing the charging call is. Also note that the terminal position differs between references, since tRNAs vary in length, so the filter has to be on distance from the 3' end or on the Sprinzl region rather than on an absolute position number. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The 2x2 table a charging odds ratio needs is a summary of what the pipeline's BAM walk already computes: for each site, how many reads fall in each cell of error by charging. Emitting that summary instead of the per-read calls turns tens of millions of rows into tens of thousands -- 488 KB against 52 MB on a zebrafish sample -- and removes the multi- gigabyte parse that dominated the read-level path. charging_odds_ratios_from_counts() consumes it. Both entry points share one cell-testing routine, so a table assembled from reads and one summarized upstream are treated identically; verified on real data, where the two agree exactly across 23,182 sites in odds ratio, p-value and depth. The per-read path stays for work that needs it, such as re-thresholding the charging call or joining other per-read features, since the summary bakes the threshold in. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
pkgdown fails the build when an exported object is absent from the reference index, which the three new entry points were. Co-Authored-By: Claude Opus 5 <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.
Adds the ability to test, for each site, whether a read being modified is associated with that same read being charged — using individual reads as the unit of observation. Because both variables are measured on the same read, this is a within-tRNA comparison that differences in abundance or charging between tRNAs or samples cannot confound.
What's new
compute_charging_odds_ratios()— builds a 2x2 table of modified/unmodified against charged/uncharged per site, with a Haldane-corrected odds ratio and Fisher's exact test matching the definitioncompute_odds_ratios()already uses.call_bcerror_sites()— selects candidate sites from a bcerror tibble by thresholding error rate and coverage. This is the point of the PR: base-calling error is model-free, so it gives a route to this question that does not depend on the modification caller being trustworthy for a given dataset.read_charging_calls()— per-read charging status.read_charging()returns per-tRNA aggregates and cannot support a read-level test.compute_odds_ratios()gainssites, so the same site list can restrict the pairwise analysis.Reads with no call at a site are excluded from that site's table rather than counted as unmodified, so
total_obsvaries between sites.Implementation
Counting is done in cpp11 in a single pass over the long-form calls. A 2x2 table needs only four counters per (reference, position), so nothing materializes a read-by-position matrix and memory scales with the number of sites rather than read depth. A 100M-row call table across 46,000 sites takes about 100 seconds and 10 GB, nearly all of it holding the parsed input.
min_marginandmax_pprune sites before testing rather than after — a genome-wide site list is mostly sites that were never testable, and each one still consumes FDR budget.max_pdrops those whose margins put the intended threshold out of reach whatever the reads did. Both filter on the margins, which are ancillary, so the null distribution of the surviving p-values is unchanged; filtering on the observed odds ratio would enrich for small p-values and break the correction applied afterwards, so that is deliberately not offered. A test asserts pruning never alters the statistics of surviving sites.The 3' end confound
Worth reading the
Sites near the 3' end are not interpretablesection before using this. The amino acid is esterified to the terminal A and the charging model reads the signal there, so a site at that end is not independent of the charging call it is being tested against.Measured across two zebrafish samples: 82% of terminal-A sites come back significant with a median log odds ratio of 1.9, decaying monotonically to 4.6% beyond fifty nucleotides. By region,
ccaanddiscriminatorrun at 63% and 36% against 3-5% in the D-loop and anticodon arms. The terminal position differs between references since tRNAs vary in length, so the filter has to be on distance from the 3' end or on the Sprinzl region, not an absolute position.Testing
581 tests pass, 8 pre-existing skips for uninstalled optional packages. New coverage includes the 2x2 counts, agreement with
stats::fisher.test, exclusion of reads lacking a call or a charging status, both site sources, the dedupe paths, and the pruning filters.The per-read mismatch tables this consumes are produced by a companion change in
aa-tRNA-seq-pipeline(get_mismatch_calls.py,select_bcerror_sites.py), which emits the same schema asmodkit extract callsso the two are interchangeable here.🤖 Generated with Claude Code