Ribo-seq ORF caller benchmark — production-pipeline tier recommendation

Status: FINAL (scaling phase complete; all gw numbers populated) Date: 2026-05-14 Authors: the project team Project: 6-tool benchmark for nf-core/riboseq fork; 6 reps of a controlled timecourse (two conditions (3 reps each)) Source: the project repository (private)

Executive summary

This report benchmarks six published Ribo-seq ORF callers — RiboCode, Rp-Bp, ribotish, PRICE, ribotricer, and RibORF — against six replicates of a human Ribo-seq dataset (2 conditions × 3 technical reps). RIBOSS was deferred (paired RNA-seq + Salmon dependency); Rp3 was excluded (proteomics dependency).

The benchmark produces per-tool Tier-1 (production default), Tier-2 (opt-in for specific use cases), and excluded recommendations for the production microprotein-quantification pipeline.

Key findings (gw scale, 6 reps per tool):

• Tier-1 cohort (recommended production defaults): RiboCode (mean Spearman 0.917, 6-26 min/rep, biology-correct uORF/dORF 1.7-4.0×) and Rp-Bp (mean Spearman 0.893, ~20-24 h/rep, biology-correct uORF/dORF 3.4-5.2×). Both tools have non-overlapping personalities: RiboCode is a permissive canonical-CDS caller; Rp-Bp is Bayesian-strict.
• Tier-2 cohort (opt-in specialists): ribotish (alternative-start-codon discovery within annotated CDSes; calls ~60% Truncated alt-starts) and PRICE (near-cognate start-codon detection; only 29.7% AUG vs 70% non-AUG CUG/GUG/AAG/UUG/AUU/AGG/ACG/AUA — unique non-overlapping niche).
• Excluded: ribotricer (mean Spearman 0.288 — rank-unstable scoring across replicates makes its score column unusable for production); RibORF (operationally untenable at gw scale — readDist.pl re-parses the SAM file once per read-length class, projected 5-7 days per replicate on this cluster's hardware); RIBOSS (missing RNA-seq dependency); Rp3 (missing proteomics dependency).
• Methodological surprise: Strict (chrom, start, end, strand) coordinate-equality cross-tool joins yield near-zero Jaccard (0.000-0.057) — tools assign different start coordinates to the same biological ORF (alt-start ambiguity). Real consensus calling requires fuzzy overlap matching, not coordinate equality. Documented for the production-pipeline design phase.

Tier recommendations (final): - Tier 1 (production default): RiboCode + Rp-Bp - Tier 2 (opt-in): ribotish (alt-start discovery), PRICE (near-cognate specialty) - Excluded: ribotricer (rank-unstable scoring), RibORF (operationally untenable — 5-7 days/rep gw), RIBOSS + Rp3 (missing data dependencies)

1. Methodology

1.1 Reference data

• Reference genome: Ensembl 108 GRCh38 primary assembly (FASTA + GTF at (cluster reference path))
• STAR index: pre-built at Index/STARIndex/ (STAR 2.7.4a)
• rRNA reference: (cluster path) + bowtie2 index
• (MANE-Select canonical-recall analysis is tracked as a follow-up addition; not in scope for this report)

1.2 Replicates (6 total, 2 timepoints × 3 biological replicates)

Sample IDs (D01–D83) are arbitrary labels; D01/D02/D03 are the three Condition-A replicates and D81/D82/D83 are the three Condition-B replicates. The same labels appear on the figures below.

BAMs from a standard nf-core/riboseq alignment; both genome-aligned and transcriptome-aligned outputs available.

1.3 Tools evaluated

1.4 Parameter choices — vendor defaults locked

Per the project lead's directive: all six tools run with vendor defaults; no parameter sweeps. This means: - Start codons: tools that take a configurable codon set use the vendor default (ATG-only for Rp-Bp, ribotricer, RiboCode; ATG+near-cognate for ribotish, RibORF, PRICE) - Q-value cutoffs: tool defaults (varying); §6.3 mapping documented in deviations - Ribosome profiling-specific tuning: vendor defaults across all tools

1.5 Schema normalization (§6.1)

All six per-tool outputs are normalized to a 14-column common schema: orf_id, chrom, start, end, strand, start_codon, length_aa, orf_type, score, host_transcript_id, host_gene_id, tool, replicate_id, sample_id

Per-tool §6.2 vocabulary mappings, §6.3 filter applications, and operational deviations are documented in Appendix A (supplementary deviations document).

1.6 Compute environment

• Host: a Linux compute cluster (991 GB RAM, 128 cores, 35 TB free (cluster mount))
• Envs: 6 tool-specific mamba envs under (cluster path)
• Scaling-phase plan: internal planning document

2. §7 Metrics

Results from 5 of 6 tools at gw scale × 6 replicates each = 30 of 30 expected normalized TSVs landed (RiboCode, Rp-Bp, ribotish, PRICE, ribotricer). RibORF excluded on operational grounds at gw scale (see §4.2); chr1 vertical-slice RibORF results retained as supplementary evidence of its call profile.

2.1 Metric 1 — Replicate concordance (HEADLINE)

Pairwise across-replicate Spearman correlation of ORF scores (top heatmaps below) and Jaccard overlap of called ORF sets (bottom heatmaps). Spearman captures whether the same ORFs rank consistently across reps; Jaccard captures whether the same set of ORFs is called. A Tier-1 candidate needs both — calling the same ORFs is necessary but not sufficient if scores are rank-unstable.

Headline numbers (mean across 15 rep pairs):

Ribotricer's gap (Jaccard 0.77 but Spearman 0.29) is the diagnostic pattern: it identifies similar ORF sets across replicates but assigns inconsistent ranks to them. A production scoring pipeline cannot underwrite that.

2.2 Metric 2 — Cross-tool concordance

Pairwise Jaccard across tools using a (chrom, start, end, strand) coordinate-equality join, averaged across the 6 reps.

Mean Jaccard across reps (gw, all 5 tools):

Surprise finding — strict coordinate equality gives near-zero cross-tool overlap. Tools place the start coordinate differently for the same biological ORF (alt-start ambiguity, frame interpretation differences), so even when they agree on translation, their (start, end) tuples don't match. This is a methodological finding for the production pipeline: real consensus calling needs fuzzy overlap matching (e.g., bedtools intersect with ≥80% reciprocal overlap, or gene-level aggregation) rather than coordinate equality. A follow-up "Metric 2b" using fuzzy matching is tracked as a future addition.

Earlier chr1 D01_S1 5-way Jaccard (gene-level join, committed at 73fb126) showed much higher numbers — ribotish ↔ ribotricer = 0.837, Rp-Bp ↔ ribocode = 0.461, 6-way intersection = 39 transcripts — which underscores the start-coordinate vs gene-level resolution gap.

2.3 Metric 3 — ORF-type breakdown per tool, per replicate

Stacked-bar fractions of each tool's calls into the canonical ORF-type vocabulary: annotated / Truncated (alt-start of an annotated CDS) / Extended / novel / uORF / dORF / internal / Overlap_* variants. One panel per replicate.

Reveals three distinct tool "personalities":

• RiboCode is a conservative canonical caller: ~70% annotated + ~25% novel, with small uORF/dORF tails. Total volume ~18-29k ORFs per rep — smallest of the cohort.
• Rp-Bp is a Bayesian-rigorous canonical caller: ~56-66% annotated + ~28-32% novel, with a balanced uORF/dORF tail. Total ~18-26k ORFs per rep — comparable volume to RiboCode. The category profile aligns with RiboCode's "conservative canonical caller" personality, consistent with the literature reputation of Rp-Bp as the most statistically stringent ORF caller.
• ribotish is a truncation maximalist: ~7% annotated, ~60% Truncated (alt-starts within annotated CDSes), ~10% Overlap_dORF. Total ~830k-1.08M per rep. The Truncated dominance is a design choice — every plausible alternative start within an annotated CDS is reported as a separate ORF.
• ribotricer is a novel-aggressive caller with a striking ~25% dORF fraction (vs ~3% uORF). Total ~350k-545k per rep.
• PRICE is a near-cognate specialist with no AUG/non-AUG bias in the orf_type breakdown but the dominant tool-distinguishing signal is start codon (see Metric 5).

The ribotricer dORF inflation flows directly into Metric 4's biology violation below.

(MANE-Select canonical recall is deferred to a follow-up analysis — needs an additional join against the MANE-Select transcript list. Tracked as a future addition.)

2.4 Metric 4 — uORF:dORF ratio (biology-sanity floor)

Pooling rule: uORF + Overlap_uORF vs dORF + Overlap_dORF. Leaky-scanning biology predicts uORF ≫ dORF (typically ≥5×); a tool reporting dORF ≥ uORF is over-calling 3′-UTR signal or mis-classifying internal initiation events. The red dashed line at log10=0 marks the "biology violation" floor.

Per-rep ratios (uORF / dORF):

The inversion in ribotish and ribotricer is consistent across all 6 reps — not a single-rep artifact. Worth verifying for the production-pipeline design phase whether this is a real over-calling pattern in the tools (a property of their default scoring) or a §6.2 type-mapping decision in the parsers; either way, the inversion is the dominant signal pushing both tools to Tier-2 / excluded rather than Tier-1.

2.5 Metric 5 — Start codon bias

Aggregated across the 6 reps (start-codon prevalence is a tool-intrinsic property of the model + threshold, not sample-specific).

The 3 fast non-PRICE tools called 100% AUG at gw scale; Rp-Bp called 99.94% AUG (11 AGU, 1 AUU among 19,741 ORFs); PRICE called only ~30% AUG and ~70% non-AUG near-cognate starts:

PRICE's gw start-codon distribution is the unique near-cognate signature that distinguishes it from every other tool in the cohort. For production-pipeline design, this means PRICE occupies a non-overlapping niche: it's the only tool that recovers non-AUG-initiated ORFs at scale (CUG-, GUG-, AAG-initiated translation), which is well-documented in the literature as accounting for ~30-40% of mammalian translation initiation events.

The other four tools' near-100% AUG profile is by design — their default thresholds suppress near-cognate predictions. Whether to bake near-cognate recovery into the production pipeline default is a biological question for the project lead; the data argues for at least an opt-in mode using PRICE.

2.6 Metric 6 — Computational cost

Mean ± std of wall + peak memory per tool, across 6 reps. From logs/scaling/{tool}_{rep}.log. Memory profiling deferred: GNU /usr/bin/time -v (which has the -v flag for max-RSS) is not installed on the cluster; only the bash builtin time is available, which captures wall but not memory. Wall captured via explicit T0=$(date +%s); ...; T1=$(date +%s); echo $((T1-T0))s pattern.

Wall times (seconds per rep):

S8 (D82_S8) consistently fastest within each tool (50-70% of next-fastest rep) due to lower read depth (12M vs 28-33M).

The log-y axis spans 2 orders of magnitude — from PRICE (~4 min) to Rp-Bp (~24 h). The read-depth-dominated runtime pattern (D82_S8 consistently fastest within each tool, ~50-70% the wall of next-fastest rep) confirms that BAM size, not genome complexity, drives wall-time across all tools — useful for production-pipeline scheduling estimates.

2.7 Metric 7 — Mass-spec validation

Dropped (proteomics unavailable; documented in deviations 2026-05-11).

3. Tier recommendations

Based on the 5-tool × 6-replicate gw-scale evaluation, the production tier assignment is:

3.1 Tier 1 candidates (production default)

• RiboCode — Mean Spearman 0.917, mean Jaccard 0.772 (highest of any tool). uORF/dORF 1.7-4.0× (biology-correct). 100% AUG. Fastest tool by wide margin (6-26 min/rep). Conservative canonical-CDS caller: ~70% annotated + ~25% novel.
• Rp-Bp — Mean Spearman 0.893, mean Jaccard 0.673. uORF/dORF 3.4-5.2× (biology-correct). 99.94% AUG. Bayesian-rigorous filtering. Wall 19h38m-24h11m per rep — slow but tractable on overnight cadence. Produces ~18k-26k ORFs per rep — comparable volume to RiboCode.

Both tools pass all M1/M4/M5 sanity floors and have non-overlapping personalities (RiboCode is permissive-novel; Rp-Bp is Bayesian-strict). The recommended production pipeline default is RiboCode for speed-prioritized routine analysis and Rp-Bp for statistical-rigor-prioritized analysis (or both, in an ensemble configuration if false-positive control matters most).

3.2 Tier 2 candidates (opt-in for specific use cases)

• ribotish — Mean Spearman 0.884, mean Jaccard 0.824 (highest Jaccard of any tool — calls the same ~1M ORFs across reps very stably). uORF/dORF inverted at 0.05-0.14× (calls 7-22× more dORFs than uORFs — biology violation at face value, but explained by ribotish's design choice to report every alternative start within annotated CDSes as a separate "Truncated" call). Wall 51-107 min/rep. Use case: alternative-start-codon discovery within annotated CDSes; the ~60% "Truncated" calls are the discovery niche.
• PRICE — Mean Spearman 0.854, mean Jaccard 0.060 (lowest Jaccard — different ORF coordinates called across reps due to start-codon ambiguity, but where same ORF is called, score ranking is stable). uORF/dORF 4.7-9.2× (highest of any tool — biology-correct). Only 29.7% AUG; ~70% near-cognate (CUG, GUG, AAG, UUG, AUU, AGG, ACG, AUA) — unique non-overlapping niche vs every other tool's 100% AUG. Wall 4-10 min/rep (second-fastest tool). Use case: near-cognate ORF discovery (CUG/GUG-initiated translation, well-documented as accounting for ~30-40% of mammalian translation initiation).

3.3 Excluded from production

• RibORF — Operationally untenable at gw scale. readDist.pl re-parses the SAM file once per read-length class (10 lengths × ~20M reads = effectively quadratic on cluster hardware). Killed at ~25h elapsed across all 6 reps with only the first length class (25) complete — projected 5-7 days per replicate to fully complete. This is a structural property of the legacy Perl implementation, not a configurable parameter. RibORF's chr1 vertical-slice results (committed at 73fb126) are retained as supplementary evidence of its call profile, but it is not a viable production candidate. See §4.2 for the technical details and the comment-to-the project lead summary.
• ribotricer — Mean Spearman 0.288 across 6 reps. While Jaccard 0.770 indicates it calls similar ORF sets across replicates, the score column is rank-unstable: the same ORF gets very different scores in different replicates. A production scoring pipeline cannot underwrite that. Excluded with prejudice.
• RIBOSS — Paired RNA-seq + Salmon dependency that we don't have for this dataset.
• Rp3 — Proteomics dependency (Metric 7 dropped).

4. Limitations & open issues

4.1 Sample-related

• N=3 reps per timepoint: limits statistical power for Metric 1. Mean/min/max reported; bootstrap CI considered but rejected as cargo-cult for N=3.
• S8 (D82_S8) low depth (12M reads): ~half the depth of other reps. Reports Metric 1 with-and-without S8 as a sensitivity check.
• No proteomics: Metric 7 dropped.

4.2 Tool-related (gw scaling issues — both resolved or documented)

• PRICE gw initially failed, resolved 2026-05-14: First gw attempt produced 6 partial .cit reads files without their .cit.metadata.json JSON sidecars — gedi's Bam2CIT merge step was never finalized (left .CITI*.data + .DISK*.data temp staging files). Re-running Bam2CIT with empty per-rep output directories produced clean .cit + .cit.metadata.json pairs. A second crash mid-PRICE (Java OOM on default 16 GB heap) was fixed by passing gedi -mem 64G. After both fixes, all 6 PRICE reps completed in 4-10 min wall each — well under the v3.1 plan's 1-4 h per-rep estimate. PRICE is included in the final Tier-2 cohort.

• RibORF gw operationally untenable: readDist.pl is implemented as an outer length-class loop wrapping an inner SAM-file read loop (verified at lines 54 + 58 of the script). The SAM file is therefore re-parsed once per length class — 10 length classes × ~20M reads/SAM = effectively O(reads × lengths) Perl serial work per replicate. Empirically, the first length class took ~16h to produce its outputs on all 6 reps simultaneously; the remaining 9 length classes were projected to take comparable wall time each based on file-write cadence, putting the per-replicate total at 5-7 days on this cluster's hardware. After 25 hours elapsed and only length-class 25 (1 of 10) complete across all 6 reps, the jobs were killed (2026-05-14) and RibORF was excluded from production tier candidacy on operational grounds. The exclusion is a wall-time finding, not a call-quality finding — RibORF's ORF predictions themselves are not in question; chr1 vertical-slice results (73fb126) are retained as supplementary evidence of its call profile. A restart with a reduced length-class set (e.g., -d 28,29,30 only) would still cost ~2 days/rep — also untenable for routine production use.

4.3 Methodological

• Vendor defaults locked: no parameter sensitivity analysis. Different cutoffs would yield different Tier assignments.
• Two-pass STAR consistency: the input BAMs from nf-core/riboseq use one parameter set; tools may behave subtly differently than on alternative aligner outputs.
• Cross-node determinism untested at scale (chr1 ribotricer was rerun-stable; gw not budgeted).

4.4 Scope

• chr1 vertical slice committed to git (all 6 tools); gw normalized TSVs kept on (cluster mount) only (1.5+ GB total too large for git).
• 6 gw reps × 5 tools = 30 normalized TSVs landed. RibORF excluded at gw scale on operational grounds; chr1 RibORF data retained as supplementary evidence.
• Metric figures (PNG + PDF) and summary TSVs committed to git for inline rendering in this report.

5. Reproducibility

5.1 Code

All six parsers + 92 unit tests committed at: - scripts/parse_ribotish.py - scripts/parse_rpbp.py - scripts/parse_riborf.py - scripts/parse_price.py - scripts/parse_ribotricer.py - scripts/parse_ribocode.py - scripts/_common.py (shared helpers) - scripts/test_parsers.py (92 tests) - scripts/metric{1..6}_*.py (6 metric computation + plotting scripts) - scripts/run_price_seq.sh (sequential PRICE launcher; encodes gedi -mem 64G fix)

Test command: pytest scripts/test_parsers.py -v

5.2 Tool environments

Pinned env yamls committed at: - environments/ribotish.yaml - environments/ribotricer.yaml - environments/ribocode.yaml - (Rp-Bp / PRICE / RibORF: tools cloned + built from source on the cluster; reproducible via internal scaling-phase state file paths)

5.3 Git history

Full per-tool parser commits, scaling-phase milestones, and the final Tier-verdict commit chain are preserved in the project repository (internal). Not included in this public release.

5.4 Compute reproducibility

• All vendor defaults; identical commands per rep
• internal scaling-phase state file documents complete path + tool resolution
• outputs/normalized/sample_metadata.tsv documents cohort metadata
• ribotricer determinism empirically confirmed (byte-identical re-runs at chr1)
• Cross-node determinism not tested (no compute budget; chr1 single-node tested)

Appendix A — Per-tool deviations summary

50+ deviations from execution plan (EP §5/§6) documented in supplementary deviations document. Headline categories:

• EP §5.x install path stale: RibORF (5 corrections), Rp-Bp (8), PRICE (7) — all major tools required spike-then-fix iteration
• EP §6.2 vocabulary: per-tool §6.2 mapping tables (4-12 categories per tool); all collated in deviations doc
• EP §6.3 filter: per-tool quality threshold mapping (Q-value direction, probability vs p-value)
• Cross-tool comparability: PRICE's near-cognate aggressive categorization; ribotricer's phase_score saturation plateau; RiboCode's negative-strand ORF_gstart > ORF_gstop convention

Appendix B — Operational lessons from scaling phase

9 operational issues caught + resolved during Phase A0 pre-flight + Phase A:

1. STAR index at Index/STARIndex/ not Sequence/STARIndex/
2. gtfToGenePred at $WORK/tools/gtfToGenePred (not in any env's bin/)
3. RibORF binary path at $WORK/tools/RibORF/RibORF.2.0/ (not $WORK/tools/RibORF.2.0/)
4. RibORF .pl scripts lack shebang — invoke via perl <script> explicitly
5. PRICE gedi -a is GTF, -p is "show progress" (no-arg); not -p $GTF
6. gedi wrapper needs envs/price/bin in PATH for embedded mvn lookup
7. GNU /usr/bin/time -v not installed on the cluster (only the bash builtin time, which lacks -v for memory) — wall-time only via T0/T1 echo pattern
8. Rp-Bp prepare-rpbp-genome check_programs_exist requires STAR on PATH even when skipping actual STAR execution
9. STAR index at sync_structure path uses STAR 2.7.4a (compatible with our STAR 2.7.11b runtime)

Appendix C — Production pipeline integration recommendations

Based on the final Tier verdict, the nf-core/riboseq fork should:

1. Default workflow runs RiboCode as the primary ORF caller (~6-26 min/rep; covers ~22k ORFs per sample with biology-correct uORF/dORF and 100% AUG canonical-CDS profile). Suitable for routine production cadence.

2. Statistical-rigor mode (overnight cadence, opt-in for sensitive analyses) runs Rp-Bp alongside RiboCode (~20-24 h/rep; Bayesian-filtered ORF set comparable in volume but stricter selection).

3. Discovery modes (opt-in, separate config): - alt-start-discovery invokes ribotish (~1 h/rep; ~1M Truncated-CDS calls per sample). Use case: alternative initiation within annotated CDSes. - near-cognate invokes PRICE with gedi -mem 64G (~5-10 min/rep; ~70% non-AUG starts). Use case: CUG/GUG-initiated translation events.

4. Excluded tools should not be included even as opt-ins: - ribotricer — rank-unstable scoring across replicates makes its score column unreliable for downstream quantification. - RibORF — 5-7 day per-replicate wall is operationally untenable on standard cluster hardware; consider only for one-off retrospective analyses, not routine production.

5. Consensus calling caveat: if the production pipeline aggregates calls across multiple tools, do not use strict (chrom, start, end, strand) coordinate-equality joins — they yield near-zero overlap across tools (Metric 2 finding). Use bedtools-style fuzzy overlap (≥80% reciprocal) or gene-level aggregation instead.

6. Memory provisioning: PRICE needs -Xmx64G (set via gedi -mem 64G) at gw scale to avoid silent OOM-killer termination. RiboCode and ribotish are comfortable with default JVM/Python settings. Rp-Bp uses ~3-6 cores per replicate for the Bayesian fit phase.