Bio Applied Ribo Seq

Ribosome Profiling (Ribo-seq) Analysis

When to Use

• Processing raw ribosome-profiling FASTQ files: 3' adapter trimming and rRNA/tRNA contaminant depletion before genome/transcriptome alignment
• Calibrating P-site offsets per read length from a Ribo-seq BAM and checking 3-nt (triplet) reading-frame periodicity as a QC gate
• Detecting actively translated ORFs (annotated CDS, uORFs, novel/non-canonical ORFs) with RiboCode, ribotricer, or plastid
• Computing translation efficiency (TE = Ribo-seq footprint density / matched RNA-seq mRNA abundance) per gene or transcript
• Deciding whether a Ribo-seq library passed QC (periodicity, footprint length distribution, frame distribution) before downstream ORF/TE analysis

Version Compatibility

• cutadapt >= 4.4, bowtie2 >= 2.5.x (rRNA/tRNA depletion index), STAR >= 2.7.11a (splice-aware genome alignment)
• plastid >= 0.6.1 (Python 3.9-3.11; metagene/P-site tooling, psite/metagene CLI and library)
• RiboCode >= 1.2.31 (annotation-based ORF calling with permutation-test periodicity)
• ribotricer >= 1.3.3 (reference-free-ish periodicity-score ORF detection, prepare-orfs + detect-orfs)
• pysam >= 0.22, numpy >= 1.26, pandas >= 2.0, scipy >= 1.11, statsmodels >= 0.14

Prerequisites

• pip install pysam pandas numpy scipy statsmodels plastid
• pip install RiboCode ribotricer; cutadapt, bowtie2/STAR on PATH
• A genome FASTA + GTF, an rRNA/tRNA-only bowtie2 index, and (for TE) a matched RNA-seq count matrix from the same samples
• Familiarity with bio-applied-rna-seq-analysis (general RNA-seq counting/DE) and bio-applied-ngs-fundamentals (FASTQ/BAM basics)

Adapter Trimming and rRNA/tRNA Depletion

Goal: turn a raw Ribo-seq FASTQ (17-34 nt footprints, usually 3'-adapter-ligated) into a clean, rRNA-free FASTQ ready for splice-aware alignment. Approach: trim the 3' adapter with size selection matching monosome footprint length, then align to a combined rRNA+tRNA bowtie2 index and keep only the unaligned reads (--un-gz) — rRNA/tRNA make up 50-90% of raw Ribo-seq reads and must be removed before genome alignment or they dominate coverage and periodicity signal.

import subprocess
import re


def remove_adapters_and_rrna(fastq_in: str, rrna_index: str, adapter: str = "AGATCGGAAGAGCACACGTCT",
                              min_len: int = 20, max_len: int = 38, threads: int = 8) -> dict:
    """Trim 3' adapter and deplete rRNA/tRNA reads from a raw Ribo-seq FASTQ.

    fastq_in: raw single-end Ribo-seq FASTQ(.gz).
    rrna_index: bowtie2 index built from species rRNA+tRNA sequences (e.g. from
        UCSC/RefSeq rRNA repeat annotations plus a tRNA database).
    adapter: 3' adapter (default is the standard Illumina/TruSeq small-RNA adapter
        used by most Ribo-seq kits).
    min_len/max_len: post-trim size window; 20-38 nt keeps monosome footprints
        (~28-32 nt) while allowing disome/short-footprint variation.
    Returns a dict of read counts through each step (n_input, n_trimmed, n_rrna,
    n_retained) parsed from cutadapt/bowtie2 stderr.
    """
    trimmed_fq = fastq_in.replace(".fastq", ".trimmed.fastq")
    cutadapt_cmd = ["cutadapt", "-a", adapter, "-m", str(min_len), "-M", str(max_len),
                     "--discard-untrimmed", "-j", str(threads), "-o", trimmed_fq, fastq_in]
    cutadapt_res = subprocess.run(cutadapt_cmd, check=True, capture_output=True, text=True)
    n_input = int(re.search(r"Total reads processed:\s*([\d,]+)", cutadapt_res.stdout).group(1).replace(",", ""))
    n_trimmed = int(re.search(r"Reads written \(passing filters\):\s*([\d,]+)", cutadapt_res.stdout).group(1).replace(",", ""))

    clean_fq = fastq_in.replace(".fastq", ".norrna.fastq.gz")
    bowtie2_cmd = ["bowtie2", "-x", rrna_index, "-U", trimmed_fq, "-p", str(threads),
                   "--un-gz", clean_fq, "-S", "/dev/null"]
    bowtie2_res = subprocess.run(bowtie2_cmd, check=True, capture_output=True, text=True)
    aligned_pct = float(re.search(r"([\d.]+)% overall alignment rate", bowtie2_res.stderr).group(1))
    n_rrna = round(n_trimmed * aligned_pct / 100)

    return {"n_input": n_input, "n_trimmed": n_trimmed, "n_rrna": n_rrna,
            "n_retained": n_trimmed - n_rrna, "clean_fastq": clean_fq}

P-site Offset Calibration and 3-nt Periodicity

Goal: convert raw footprint 5' alignment positions into P-site (decoding-center) positions, then confirm the library shows the expected triplet periodicity around start codons. Approach: for reads mapped to a transcriptome/genome BAM near annotated start codons, find the modal distance from each read's 5' end to the start codon, stratified by read length (footprint offset is length-dependent, typically 12-15 nt in eukaryotes); apply that offset to every read of the same length to get P-site positions, then check what fraction fall in frame 0 vs 1/2 — a healthy library shows >60-70% of P-sites in frame 0.

import pysam
from collections import defaultdict, Counter


def calibrate_psite_offsets(bam_path: str, start_codon_pos: dict, read_length_range: tuple = (25, 35)) -> dict:
    """Calibrate per-read-length P-site offsets from 5' footprint ends near start codons.

    bam_path: transcriptome-coordinate BAM (splice-aware genome BAMs need CDS
        coordinates converted to transcript space first, e.g. with plastid).
    start_codon_pos: dict transcript_id -> 0-based transcript-coordinate position
        of the start codon's first base (from a GTF/plastid Transcript object).
    read_length_range: footprint lengths to calibrate over (monosome range).
    Returns {read_length: offset_nt}, the modal 5'-end-to-start-codon distance
    per length — the offset to add to a read's 5' position to get its P-site.
    """
    bam = pysam.AlignmentFile(bam_path, "rb")
    offset_counts = defaultdict(Counter)
    for tx_id, start_pos in start_codon_pos.items():
        if tx_id not in bam.references:
            continue
        for read in bam.fetch(tx_id):
            if read.is_unmapped or read.is_reverse:
                continue
            rl = read.query_length
            if not (read_length_range[0] <= rl <= read_length_range[1]):
                continue
            offset = start_pos - read.reference_start
            if 0 <= offset < 20:
                offset_counts[rl][offset] += 1
    return {rl: counts.most_common(1)[0][0] for rl, counts in offset_counts.items() if counts}


def frame_periodicity(bam_path: str, start_codon_pos: dict, psite_offsets: dict) -> dict:
    """Compute the fraction of calibrated P-sites falling in frame 0/1/2 across all CDS.

    Uses the offsets from calibrate_psite_offsets to shift each read's 5' end to
    its P-site, then bins (P-site - start_codon_pos) % 3. Frame 0 should dominate
    (>60-70%) for a usable library; a flat ~33/33/33 split means periodicity was
    lost (over-digestion, poor size selection, or wrong offsets).
    """
    bam = pysam.AlignmentFile(bam_path, "rb")
    frame_counts = Counter()
    for tx_id, start_pos in start_codon_pos.items():
        if tx_id not in bam.references:
            continue
        for read in bam.fetch(tx_id):
            if read.is_unmapped or read.is_reverse:
                continue
            offset = psite_offsets.get(read.query_length)
            if offset is None:
                continue
            psite = read.reference_start + offset
            frame_counts[(psite - start_pos) % 3] += 1
    total = sum(frame_counts.values()) or 1
    return {f"frame_{f}_pct": round(100 * frame_counts.get(f, 0) / total, 1) for f in (0, 1, 2)}

ORF Detection and Translation Efficiency

Goal: identify actively translated ORFs and quantify translation efficiency (TE) relative to mRNA abundance. Approach: for ORF calling, use a dedicated periodicity-aware tool rather than hand-rolled logic — RiboCode tests annotated-plus-novel ORFs (uORFs, dORFs, ncRNA ORFs) with a permutation test on triplet periodicity (RiboCode_prepare_transcripts then RiboCode -c config.txt -l no -g); ribotricer scores arbitrary ORF candidates (ribotricer prepare-orfs + ribotricer detect-orfs) using a multitaper/periodicity metric, tool-agnostic to species annotation quality; plastid provides the underlying P-site/metagene machinery (plastid.plotting, psite, metagene) both tools build on, and is the right choice when you need custom metagene profiles rather than a full ORF caller. Once footprint counts per gene are in hand, TE is simply the library-size-normalized Ribo-seq/RNA-seq ratio:

import numpy as np
import pandas as pd
from scipy import stats


def translation_efficiency(ribo_counts: pd.DataFrame, rna_counts: pd.DataFrame, pseudocount: float = 1.0) -> pd.DataFrame:
    """Compute per-gene log2 translation efficiency (Ribo-seq / RNA-seq) across matched samples.

    ribo_counts, rna_counts: genes x samples raw count DataFrames with identical
        column order (same samples, footprint counts vs. mRNA read counts).
    pseudocount: added before logging to avoid log(0) for low-count genes.
    Returns a DataFrame with per-sample log2(TE), the across-sample mean log2(TE),
    and a one-sample t-test p-value against log2(TE)=0 (tests whether a gene is
    translationally up/down-regulated beyond what mRNA level predicts).
    """
    ribo_cpm = ribo_counts / ribo_counts.sum(axis=0) * 1e6
    rna_cpm = rna_counts / rna_counts.sum(axis=0) * 1e6
    log2_te = np.log2(ribo_cpm + pseudocount) - np.log2(rna_cpm + pseudocount)

    result = log2_te.copy()
    result["mean_log2_TE"] = log2_te.mean(axis=1)
    result["pvalue"] = [stats.ttest_1samp(row.values, popmean=0.0).pvalue for _, row in log2_te.iterrows()]
    return result.sort_values("pvalue")

Pitfalls

• Genomic vs. transcript coordinates: P-site offset calibration and frame periodicity must be computed in transcript (spliced) coordinates — a genomic BAM with introns will corrupt the modulo-3 frame calculation across splice junctions
• Skipping rRNA/tRNA depletion: raw Ribo-seq libraries are often 50-90% rRNA/tRNA; aligning without depletion wastes compute and can leave contaminant reads polluting coverage-based ORF calls
• One offset for all read lengths: the 5'-end-to-P-site distance depends on footprint length (nuclease digestion varies); applying a single fixed offset instead of a per-length table destroys periodicity
• TE from unmatched libraries: translation efficiency requires Ribo-seq and RNA-seq from the same samples/conditions with independent library-size normalization (CPM/TMM each); comparing to a public RNA-seq reference confounds TE with batch/condition effects
• Trusting low-count genes: log2(TE) is very noisy below ~1-2 reads/kb in either assay; filter genes on a minimum RNA-seq and Ribo-seq count/CPM threshold before ranking or testing TE changes

See Also

• bio-applied-rna-seq-analysis — general RNA-seq quantification/DE, needed for the matched mRNA side of TE
• bio-applied-ngs-fundamentals — FASTQ/BAM/GTF basics underlying alignment and coordinate handling here
• bio-applied-advanced-ngs — splice-aware alignment (STAR) and coverage-track workflows used upstream of P-site calibration
• bio-applied-mirna-seq-pipeline — parallel short-read adapter-trimming/size-selection workflow for another small-RNA-length library type