Compute odds ratios between modification and charging.

For each site, use individual reads as the unit of observation to ask whether a read being modified at that site is associated with that same read being charged. Each site yields a 2x2 table of modified/unmodified against charged/uncharged, an odds ratio, and Fisher's exact test.

Usage

compute_charging_odds_ratios(
  calls,
  charging,
  sites = NULL,
  refs = NULL,
  min_reads = 10,
  min_margin = 1,
  max_p = 1,
  ml_threshold = 200,
  dedupe = TRUE,
  p_method = "BH"
)

Arguments

calls: Per-read calls, either a path to a mod_calls.tsv.gz / mismatch_calls.tsv.gz file or a tibble. Requires columns read_id, chrom (or ref), ref_position (or pos), and either call_code (where "-" means no call) or an integer modified column. A logical within_alignment column, if present, is used to drop calls outside the alignment.
charging: Per-read charging status, either a path to a charging_prob.tsv.gz file or a tibble from read_charging_calls().
sites: Optional tibble of sites to test, with columns ref and pos, typically from call_bcerror_sites(). If NULL (default), all positions in calls are tested.
refs: Optional character vector of reference names to include. If NULL (default), all references are processed.
min_reads: Minimum number of reads with both a call and a charging status required for a site to be tested. Default 10.
min_margin: Minimum count in each margin of the 2x2 table: modified, unmodified, charged and uncharged reads. Default 1, which only requires the odds ratio to be defined. Raising it drops sites too thin to carry power, such as a position where three of ten thousand reads are modified.
max_p: Skip sites where no arrangement of the reads could reach this p-value, given the margins. Default 1, which tests everything. Set it to the significance threshold you intend to use.
ml_threshold: Charging likelihood threshold, used only when charging is a path. Default 200.
dedupe: Whether to collapse repeated calls for the same read and position before counting. Needed for modkit output, where a read with several modification channels at one position contributes a row per channel. Default TRUE. Set FALSE for per-read mismatch calls, which carry one row per read and position, to skip a grouping pass that is costly on large tables.
p_method: Method for p-value adjustment, passed to stats::p.adjust(). Default "BH".

Value

A tibble with one row per tested site and columns: ref, pos, n11 (modified and charged), n10 (modified and uncharged), n01 (unmodified and charged), n00 (unmodified and uncharged), total_obs, mod_freq, charged_freq, odds_ratio, log_odds_ratio, p_value, and p_adjusted.

Details

An odds ratio above 1 means modified reads are more likely to be charged than unmodified reads of the same tRNA; below 1 means the opposite. Because both variables are measured on the same read, this is a within-tRNA comparison and is not confounded by differences in abundance or charging between tRNAs or between samples.

Sites can come from either of two sources, set by sites:

Direct modification calls (sites = NULL, the default): every position present in calls is tested.
Base-calling error: pass a site list from call_bcerror_sites() to restrict testing to positions called from error rates. In this case calls should be a per-read mismatch table rather than modification calls, so that the per-read status and the 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.

Pruning sites before testing

A genome-wide site list is mostly sites that were never testable, and each one still consumes FDR budget. Two arguments drop them up front. min_margin requires a minimum count of modified, unmodified, charged and uncharged reads; max_p drops sites whose margins make the intended significance threshold unreachable no matter how the reads fall.

Both filter on the margins of the 2x2 table, never on the observed odds ratio, and that distinction matters. The margins are ancillary, so filtering on them leaves the null distribution of the surviving p-values intact. Filtering on effect size would enrich for small p-values and invalidate the correction applied afterwards, so do not pre-filter sites on an effect measured from the same reads.

Sites near the 3' end are not interpretable

The amino acid is esterified to the terminal A, and charging is called from the nanopore signal over that CCA end. A site there is therefore not independent of the charging call it is being tested against: the base caller is reading out the same thing the charging model is, and the odds ratio between them is close to tautological.

This dominates the result. Across two zebrafish samples, 82% of terminal-A sites came back significant with a median log odds ratio of 1.9, and the effect decayed monotonically with distance from that residue – roughly 50% of sites one or two nucleotides in, 16-33% from three to ten, and 4.6% beyond fifty, where it flattens out. By region, cca and discriminator ran at 63% and 36% against 3-5% in the D-loop and anticodon arms.

Note that the position number of the terminal A varies between references, since tRNAs differ in length, so filter on distance from the 3' end or on the region column of a Sprinzl coordinate table rather than on an absolute position. Dropping cca, discriminator and the 3' acceptor stem removes the worst of it; the gradient suggests discarding everything within about twenty nucleotides of the end. This caveat applies whichever site source is used, since it concerns the charging call rather than the modification call.

Examples

calls <- clover_example("ecoli/mod_calls.tsv.gz")
charging <- clover_example("ecoli/charging_calls.tsv.gz")
compute_charging_odds_ratios(calls, charging, min_reads = 5)
#> # A tibble: 2 × 13
#>   ref     pos   n11   n10   n01   n00 total_obs mod_freq charged_freq odds_ratio
#>   <chr> <int> <dbl> <dbl> <dbl> <dbl>     <dbl>    <dbl>        <dbl>      <dbl>
#> 1 host…    34     3     0     3     9        15      0.2          0.4       19  
#> 2 host…    46     1     2     5     7        15      0.2          0.4        0.7
#> # ℹ 3 more variables: log_odds_ratio <dbl>, p_value <dbl>, p_adjusted <dbl>

# Restrict to sites called from base-calling error instead.
bcerr <- read_bcerror(clover_example(
  "ecoli/summary/tables/wt-15-ctl-01/wt-15-ctl-01.bcerror.tsv.gz"
))
sites <- call_bcerror_sites(bcerr, min_error = 0.05, min_cov = 10)
mismatches <- clover_example("ecoli/mismatch_calls.tsv.gz")
compute_charging_odds_ratios(
  mismatches,
  charging,
  sites = sites,
  min_reads = 5
)
#> # A tibble: 3 × 13
#>   ref     pos   n11   n10   n01   n00 total_obs mod_freq charged_freq odds_ratio
#>   <chr> <int> <dbl> <dbl> <dbl> <dbl>     <dbl>    <dbl>        <dbl>      <dbl>
#> 1 host…    58     6     1     0     8        15    0.467          0.4     73.7  
#> 2 host…    34     3     0     3     9        15    0.2            0.4     19    
#> 3 host…    55     0     3     6     6        15    0.2            0.4      0.143
#> # ℹ 3 more variables: log_odds_ratio <dbl>, p_value <dbl>, p_adjusted <dbl>