Bio Ribo Seq Ribosome Stalling

Version Compatibility

Reference examples tested with: plastid 0.6+, numpy 1.26+, scipy 1.12+, biopython 1.83+, twobitreader 3.1+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Ribosome Stalling Detection

"Find ribosome pause sites in my data" -> Detect codon positions where ribosomes dwell longer than the local average, attribute them to A-site decoding or nascent-chain effects, and judge whether the signal is real or a drug artifact.

• Python: plastid for A-site codon density, local-relative pause scoring, and motif context

Read this first: cycloheximide destroys pause signal

A pause is only meaningful when footprint positions reflect in-vivo dwell times. Cycloheximide (CHX) pre-treatment of live cells violates this: arrest is not instantaneous, ribosomes run on after the drug, density redistributes downstream, codon-specific pausing is attenuated, and an artifactual start-codon peak appears. Hussmann 2015 showed CHX data report a WEAK NEGATIVE correlation between codon rate and tRNA abundance while flash-frozen data report a STRONG POSITIVE one -- the drug flips the conclusion. A pause analysis on CHX data largely measures the drug.

Decision rule before any dwell/pause analysis:

• Flash-frozen, no drug (or CHX only in lysis at high concentration) -> codon-resolution dwell-time analysis is valid.
• CHX pre-incubation of live cells -> restrict to qualitative/gene-level statements; do not report codon dwell times or tRNA correlations as biology.
• Harringtonine/lactimidomycin data are for initiation-site mapping, not elongation pausing (see initiation-site-mapping).

A-site vs P-site: the offset choice changes the biology

The P-site offset (~12 nt from the 5' end for canonical 28-30 nt footprints) must be calibrated per read length, not hardcoded (see ribosome-periodicity). The relevant site depends on the mechanism: tRNA-availability/decoding pauses register at the A-SITE (A-site = P-site + 3), so codon-occupancy and tRNA work assign to the A-site. Nascent-chain effects (polyproline, charge) act at the P-site/exit tunnel and upstream. State which site is used; the peak position relative to A/P/E is itself diagnostic.

Pause-metric selection

The two normalization rules the naive z-score violates: never z-score across positions of differently-expressed genes (high-expression genes dominate) -- normalize each gene to its own mean first; and require a real per-gene coverage floor (a few hundred in-frame footprints), far above a sum > 100 cutoff, or per-position metrics are noise.

Calculate A-site codon density (plastid)

Goal: Get a per-codon occupancy vector for each CDS at the A-site.

Approach: Map footprints to the A-site offset, fetch the CDS count vector, and reduce each codon to its summed in-frame count.

from plastid import BAMGenomeArray, GTF2_TranscriptAssembler, FivePrimeMapFactory
import numpy as np

def asite_codon_occupancy(bam_path, gtf_path, asite_offset=15):
    '''Per-codon A-site occupancy per CDS. A-site offset = P-site (~12) + 3.

    A single fixed offset is a simplification valid only when one read length
    dominates. For production, calibrate per length (ribosome-periodicity) and
    map with VariableFivePrimeMapFactory.from_file using A-site = P-site + 3.
    '''
    alignments = BAMGenomeArray(bam_path, mapping=FivePrimeMapFactory(offset=asite_offset))
    out = {}
    for tx in GTF2_TranscriptAssembler(gtf_path):
        if tx.cds_start is None:
            continue
        cds = tx.get_cds()
        counts = cds.get_counts(alignments)          # numpy vector over the CDS
        n_codons = len(counts) // 3
        # Sum the 3 positions of each codon into a SCALAR (one value per codon)
        per_codon = np.array([counts[i*3:i*3+3].sum() for i in range(n_codons)])
        out[tx.get_name()] = per_codon
    return out

cds.get_counts(alignments) is the count method on the SegmentChain; BAMGenomeArray has no count_in_region/get_density. Reducing each codon to a scalar (sum of its three positions) is essential -- storing the whole vector at each codon makes every downstream metric garbage.

Score pauses with a local-relative metric

Goal: Flag codons where occupancy exceeds the gene's own average.

Approach: Divide each position by the gene mean (a pause score), require adequate coverage, and threshold.

def pause_scores(per_codon_occupancy, min_total=500, score_threshold=5.0):
    '''Pause score = codon occupancy / gene-mean occupancy (local-relative).

    min_total: per-gene footprint floor; below this, scores are noise.
    score_threshold: fold-over-gene-mean to call a pause (tune per dataset).
    '''
    pauses = []
    for tx, occ in per_codon_occupancy.items():
        if occ.sum() < min_total:
            continue
        mean = occ.mean()
        if mean == 0:
            continue
        scores = occ / mean
        for pos in np.where(scores > score_threshold)[0]:
            pauses.append({'transcript': tx, 'codon': int(pos),
                           'pause_score': float(scores[pos])})
    return pauses

Codon occupancy across genes

Goal: Estimate per-codon-type dwell, averaged across the transcriptome.

Approach: Normalize each gene to its own mean BEFORE pooling, then average per codon identity (the A-site codon).

Pool mean-of-ratios, not ratio-of-means: a raw average across genes is dominated by highly expressed genes. The per-codon occupancy is then a relative dwell estimate -- and only on no-drug data. The tRNA-availability correlation (codon occupancy vs tRNA adaptation index) is modest, sign- and protocol-dependent, and reflects charged-tRNA levels rather than gene copy number; report the effect size, not a presumed strong negative correlation.