Bio Applied Ribo Seq

---
name: bio-applied-ribo-seq
description: "Ribo-seq: cutadapt/bowtie2 adapter+rRNA removal, plastid P-site calibration, 3-nt periodicity QC, RiboCode/ribotricer ORF calling, translation efficiency. Use when user has ribosome profiling or footprint data."
tool_type: python
primary_tool: Python
---

# 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.

```python
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.

```python
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:

```python
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