# RiboCode P-site Detection Fix Proposal

## Problem Summary

RiboCode's `metaplots` command applies three stringent filters to determine which read lengths have valid 3-nucleotide periodicity:

| Filter | Default | Description |
|--------|---------|-------------|
| `--frame0_percent` | 0.65 | Requires 65% of reads in frame 0 |
| `--pvalue1_cutoff` | 0.001 | Wilcoxon test: frame0 > frame1 |
| `--pvalue2_cutoff` | 0.001 | Wilcoxon test: frame0 > frame2 |

When no read lengths pass all three tests, the config file is written with only comment lines and no data. RiboCode then exits with:

> Error, can not determine the P-site locations. Please check your input parameters and try again.

This is a known upstream issue ([xryanglab/RiboCode#62](https://github.com/xryanglab/RiboCode/issues/62), [#65](https://github.com/xryanglab/RiboCode/issues/65)) with no developer fix. It affects the nf-core/riboseq pipeline reliably on the full test dataset, and has been reported by multiple users.

## Root Cause (in source)

In `metaplots.py`, the function `_predict_psite()` (line 159) loops over candidate P-site positions for each read length. A read length is rejected if:

1. Fewer than 10 reads in frame 0 (`f0.sum() < 10`)
2. Frame 0 fraction below threshold (`f0_percent < frame0_percent_cutoff`)
3. Either p-value is not significant (`pv1.pvalue >= pv1_cutoff` or `pv2.pvalue >= pv2_cutoff`)

When every read length is rejected, `pre_psite_dict` is empty (line 216). The config file gets no data line, and `loadconfig.py` line 66 triggers the fatal error (`if i == 0: sys.exit()`).

Critically, `metaplots` already computes the peak position for each read length at line 202 (`mask_max_psite = d[:50].argmax()`) but discards it when the statistical tests fail.

## Proposed Fix: Pipeline-level Fallback Config Generation

Since we cannot modify RiboCode's source code, the fix lives entirely in the nf-core/riboseq pipeline. The strategy has three layers.

### Layer 1: Relax metaplots defaults

Change the default `ext.args` for `RIBOCODE_METAPLOTS` in `conf/modules.config`:

```groovy
withName: 'RIBOCODE_METAPLOTS' {
    ext.args = { params.extra_ribocode_metaplots_args ?: '-pv1 0.05 -pv2 0.05 -f0_percent 0.5' }
}
```

This relaxes the thresholds enough to rescue many borderline samples without fully disabling the tests (which [issue #65](https://github.com/xryanglab/RiboCode/issues/65) showed can shift P-sites to non-peak positions). Users can still override via `--extra_ribocode_metaplots_args`.

This alone will not fix all cases - some samples still fail even with `f0_percent` set to 0.

### Layer 2: Use riboWaltz P-site offsets as fallback

The pipeline already runs riboWaltz by default, which computes per-read-length P-site offsets via its own method (`riboWaltz::psite()`). The output `*.psite_offset.tsv` contains:

| Column | Description |
|--------|-------------|
| `length` | Read length |
| `corrected_offset_from_5` | P-site offset from 5' end (after correction) |
| `total_percentage` | Percentage of reads at this length |
| `sample` | Sample name |

This is exactly the information RiboCode's config file needs. When metaplots fails to detect periodicity, we can use riboWaltz's offsets to generate a valid config file instead. This directly addresses the concern raised by Ira - no one needs to guess read lengths a priori, because riboWaltz has already determined them from the data.

#### Workflow restructuring

Currently riboWaltz and RiboCode run independently. The fix adds a config validation/generation process between metaplots and RiboCode that can optionally consume riboWaltz output:

```
riboWaltz ─────────────────────────┐
                                   v
metaplots ──> CONFIG_FALLBACK ──> RiboCode
```

#### New local process: RIBOCODE_CONFIG_FALLBACK

```groovy
process RIBOCODE_CONFIG_FALLBACK {
    tag "$meta.id"
    label 'process_single'

    input:
    tuple val(meta), path(config)
    tuple val(meta2), path(bam)
    tuple val(meta3), path(ribowaltz_offsets)   // optional: from RIBOWALTZ
    val(default_offset)                          // static fallback P-site offset

    output:
    tuple val(meta), path("*_validated_config.txt"), emit: config

    script:
    def prefix = meta.id
    """
    #!/usr/bin/env python3
    import csv
    import os
    import sys

    config_file = "${config}"
    bam_file = "${bam}"
    ribowaltz_file = "${ribowaltz_offsets}"
    default_offset = ${default_offset}
    prefix = "${prefix}"

    # Check if metaplots config has valid data lines
    has_data = False
    with open(config_file) as f:
        for line in f:
            if line.strip() and not line.startswith("#"):
                has_data = True
                break

    if has_data:
        # Metaplots config is valid - pass through
        import shutil
        shutil.copy(config_file, f"{prefix}_validated_config.txt")
        sys.exit(0)

    # Metaplots failed - try riboWaltz offsets
    print("WARNING: metaplots failed to detect periodicity for ${prefix}.", file=sys.stderr)

    read_lengths = []
    psite_offsets = []

    if ribowaltz_file != "NO_RIBOWALTZ" and os.path.exists(ribowaltz_file):
        print("Using riboWaltz P-site offsets as fallback.", file=sys.stderr)
        with open(ribowaltz_file) as f:
            reader = csv.DictReader(f, delimiter="\\t")
            for row in reader:
                length = int(row["length"])
                offset = int(row["corrected_offset_from_5"])
                read_lengths.append(length)
                psite_offsets.append(offset)

    if not read_lengths:
        # No riboWaltz data either - use static defaults
        print(f"Using static fallback: read lengths 28-30 with P-site offset {default_offset}.", file=sys.stderr)
        read_lengths = [28, 29, 30]
        psite_offsets = [default_offset] * 3

    sample_name = os.path.splitext(os.path.basename(bam_file))[0]
    plen_str = ",".join(map(str, read_lengths))
    psite_str = ",".join(map(str, psite_offsets))

    with open(f"{prefix}_validated_config.txt", "w") as f:
        f.write(f"# Fallback config: metaplots periodicity detection failed.\\n")
        if ribowaltz_file != "NO_RIBOWALTZ" and os.path.exists(ribowaltz_file):
            f.write(f"# P-site offsets derived from riboWaltz.\\n")
        else:
            f.write(f"# Using static default offsets.\\n")
        f.write(f"{sample_name}\\t{bam_file}\\tyes\\t{plen_str}\\t{psite_str}\\n")
    """
}
```

#### Workflow changes in `workflows/riboseq/main.nf`

```groovy
if (!params.skip_ribocode){
    // ... (steps 1-3 unchanged) ...

    // Step 3: Generate metaplots and config for each sample
    RIBOCODE_METAPLOTS(
        ch_transcriptome_bams_for_ribocode,
        RIBOCODE_PREPARE.out.annotation
    )

    // Prepare riboWaltz offsets channel (empty if riboWaltz was skipped)
    ch_ribowaltz_offsets = params.skip_ribowaltz
        ? ch_transcriptome_bams_for_ribocode.map { meta, bam -> [ meta, "NO_RIBOWALTZ" ] }
        : RIBOWALTZ.out.offset

    // Step 3b: Validate config, use riboWaltz offsets as fallback if needed
    RIBOCODE_CONFIG_FALLBACK(
        RIBOCODE_METAPLOTS.out.config,
        ch_transcriptome_bams_for_ribocode,
        ch_ribowaltz_offsets,
        params.ribocode_default_psite_offset
    )

    // Step 4: Run RiboCode ORF detection (now uses validated config)
    ch_ribocode_inputs = ch_transcriptome_bams_for_ribocode
        .join(RIBOCODE_CONFIG_FALLBACK.out.config)

    RIBOCODE_RIBOCODE(
        ch_ribocode_inputs.map { meta, bam, config -> [ meta, bam ] },
        RIBOCODE_PREPARE.out.annotation,
        ch_ribocode_inputs.map { meta, bam, config -> [ meta, config ] }
    )
}
```

### Layer 3: Static fallback (safety net)

If riboWaltz was also skipped (`--skip_ribowaltz`), the process falls back to static defaults:

| Parameter | Default | Description |
|-----------|---------|-------------|
| `ribocode_default_psite_offset` | 12 | Default P-site offset (mirrors `plastid_default_psite_offset`) |

The static fallback uses read lengths 28-30 with the default offset. This is a last resort - the riboWaltz-informed fallback is strongly preferred.

### Fallback priority

```
1. metaplots auto-detection (relaxed defaults)   -- best case
2. riboWaltz-derived offsets                      -- data-driven fallback
3. Static defaults (28-30, offset 12)             -- last resort
```

## Why this approach

- **No RiboCode source changes required** - works with the tool as-is
- **Data-driven fallback** - riboWaltz offsets are computed from the actual data, not guessed a priori (addresses Ira's concern)
- **Three layers of defence** - relaxed defaults, riboWaltz fallback, static fallback
- **Consistent with pipeline conventions** - mirrors the existing `plastid_default_psite_offset` pattern
- **User-controllable** - all behaviour is exposed as pipeline parameters
- **Transparent** - fallback emits clear warnings so users know auto-detection failed

## Risks and mitigations

| Risk | Mitigation |
|------|------------|
| Relaxed p-values may shift P-sites off-peak ([#65](https://github.com/xryanglab/RiboCode/issues/65)) | We only relax to 0.05, not 1.0; still provides meaningful filtering |
| riboWaltz and RiboCode use different offset calculation methods | Both compute 5'-end offsets; riboWaltz's method is well-validated ([Lauria et al. 2018](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1006169)) |
| Adding a riboWaltz dependency makes RiboCode wait | Only in the fallback path; metaplots and riboWaltz already run in parallel |
| Static fallback lengths (28-30) may not match actual data | Last resort only; riboWaltz fallback is preferred and data-driven |
| Downstream ORF calls less accurate with fallback | Still better than `--skip_ribocode` which gives no results |

## Alternatives considered

### Only relaxing p-value cutoffs

Setting `-pv1 1 -pv2 1` (suggested in [issue #62](https://github.com/xryanglab/RiboCode/issues/62)) effectively disables the statistical test. This helps some samples but the PR comment from Sebastian Uhrig confirmed that even `f0_percent 0` does not fix all samples. A fallback is needed as a safety net regardless.

### User-supplied config file

Adding a parameter for users to provide a pre-made config file would work but puts the burden on users to understand RiboCode's config format and know their read lengths a priori. The riboWaltz-based fallback automates this.

### Static-only fallback (no riboWaltz)

This was the original proposal but doesn't address Ira's question about how to know which read lengths to trust. The riboWaltz integration solves this by using already-computed, data-driven offsets.