# RBPNet/eCLIP data workflow

```{warning}
Everything in this guide is experimental: eCLIP preprocessing, window scanning
and selection, bundle construction, RBPNet modeling, training, and evaluation.
The workflow has been minimally tested and has only been confirmed to process
data successfully and train reasonable models on PUM2 eCLIP data. It needs
substantially more validation should
not be treated as production-ready or broadly validated.
```

TranscriptML includes an eCLIP path from ordinary alignments through fixed-shape,
memory-mappable arrays and structured RBPNet training. The scanner is
descriptive and the selectors prepare model examples; none is intended as a
general-purpose peak caller. GC matching and post-selection sampling are not
implemented.

```text
FASTA + one-transcript-per-gene GTF + IP BAM(s) + SMInput BAM
                              |
                              v
                         preprocess
                              |
                              v
       canonical locus-coordinate HDF5 + FASTA + metadata
                              |
                              v
                     descriptive window scan
                              |
                              v
                       region selection
                              |
                              v
                 versioned selection manifest
                              |
                              v
                     RBPNet DatasetBundle
                              |
                              v
              fixed-shape, memory-mappable .npy arrays
                              |
                              v
           structured RBPNet profile/enrichment training
```

Install the optional assay dependencies with:

```bash
python -m pip install -e ".[rbpnet]"
```

## 1. Canonical preprocessing

The GTF must select one transcript per gene. GTF coordinates are converted to
zero-based, half-open intervals. `--coordinate-space mature_transcript` is the
default and preserves the original behavior: exons are spliced in transcript
5′→3′ order. `--coordinate-space gene` instead retains the contiguous genomic
span from the first through last selected exon, including introns. Both spaces
are represented in annotated RNA 5′→3′ orientation. Thus, position zero on a
minus-strand locus is its highest-genomic-coordinate base and gene-space
sequence is reverse-complemented.

Protein-coding exon sequence is partitioned into `5putr`, `cds`, and `3putr`;
exons without a CDS are `noncoding_exon`. Gene space additionally labels the
gaps between exons as `intron`. The labels exhaustively partition each stored
locus, so a boundary-crossing window is reported as `mixed`.

```bash
transcriptml rbpnet preprocess \
  --genome-fasta ../RBPNet2/Data/chr21_test/chr21.fa \
  --gtf ../RBPNet2/Data/chr21_test/gencode_v50_MANE_select_chr21.gtf \
  --ip-bam ip1=../RBPNet2/Data/chr21_test/ip1_chr21.bam \
  --ip-bam ip2=../RBPNet2/Data/chr21_test/ip2_chr21.bam \
  --sminput-bam sminput=../RBPNet2/Data/chr21_test/sminput_chr21.bam \
  --coordinate-space mature_transcript \
  --output-dir processed/chr21
```

Use `--coordinate-space gene --output-dir processed/chr21_gene` to build the
corresponding full-gene experiment. Coordinate space is recorded in
`manifest.json`, `signals.h5`, scan metadata, selection provenance, and bundle
metadata. Downstream stages reject mismatched inputs.

For backward compatibility, table columns such as `transcript_id`, `tx_start`,
`tx_end`, and `transcript_anchor` retain their established names. In a
gene-space experiment, their numeric values are gene-space coordinates in the
selected transcript's 5′→3′ orientation; `coordinate_space` removes the
ambiguity.

Only read1 alignments are considered. The crosslink-position signal is the
aligned 5′ reference base: `reference_start` for forward alignments and
`reference_end - 1` for reverse alignments. In the eCLIP libraries used during
development, read1 aligns opposite the RNA strand, so `--read1-rna-strand
opposite` is the default. Use `same` or `unstranded` for libraries with another
convention. This read1 convention describes the supplied BAMs and should not be
confused with papers that name the crosslink-bearing FASTQ mate R2 before BAM
construction.

Unmapped, secondary, supplementary, QC-failed, duplicate (by default), and
low-MAPQ alignments are filtered. Assignment then depends on coordinate space:

- In mature-transcript space, the 5′ base must map uniquely to a
  strand-compatible selected exon. Aligned/deleted reference segments must
  remain in exons and each CIGAR `N` must exactly match an adjacent selected
  exon junction. Intronic/pre-mRNA alignments are not retained.
- In gene space, the 5′ base must map uniquely to a strand-compatible selected
  gene span. Every reference-consuming CIGAR operation must remain inside that
  span. Intronic alignments are retained, and `N` operations may represent any
  splice within the selected gene span; they need not match the one selected
  mature isoform. Reads overlapping multiple same-strand gene spans remain
  ambiguous rather than being assigned arbitrarily.

The same read filtering and unique-gene rule apply in both modes. Rejections
are itemized in `qc.json`.

Missing `.fai`/`.bai` indexes are created when the source files and their
directories are writable. GTF transcripts on contigs absent from the analysis
FASTA are reported and skipped; retained FASTA contigs must be represented in
every BAM.

The canonical directory contains:

| File | Contract |
| --- | --- |
| `manifest.json` | Format version, exact sample order/roles, input file stat records, configuration, and effective library sizes |
| `qc.json` | Annotation, sequence, read-filter, assignment, and per-sample totals |
| `transcripts.tsv` | Gene/transcript metadata, coordinate space, genomic span, length, strand, raw counts, SMInput TPM, and compact region annotations |
| `exons.tsv.gz` | Transcript interval ↔ genomic exon block mappings |
| `regions.tsv.gz` | Long-form transcript region intervals |
| `transcripts.fa` + `.fai` | Indexed 5′→3′ selected-locus sequences |
| `signals.h5` | Canonical base-resolution retained read1 5′ counts |

SMInput TPM is calculated from retained locus counts divided by stored locus
length, followed by normalization of those rates to one million. Therefore,
gene-space TPM uses gene-space retained events and full gene-span length; it is
not numerically interchangeable with mature-transcript TPM.
Each sample's `effective_library_size` is exactly the number of retained read1
5′ events used to construct its HDF5 track. That field is the CPM denominator.
Pooled-IP CPM uses the sum of IP counts divided by the sum of IP effective
library sizes.

Signal rows are sparse-aggregated and then written one complete touched HDF5
chunk at a time. Zero-only chunks retain the HDF5 fill value and are not
allocated. The default is shuffled gzip level 1, chosen to reduce write and
slice-decompression time while retaining compact sparse storage. Use
`--signal-compression gzip --signal-compression-level 4` for smaller but slower
files, `--signal-compression lzf` for faster/larger files, or
`--signal-compression none` when storage is unimportant. Compression and chunk
length are recorded as `signals.h5` attributes and in preprocessing provenance.

### HDF5 signal layout

HDF5 is a hierarchical binary container: datasets behave like typed,
multidimensional arrays stored inside a file, can be compressed and chunked,
and can be sliced without loading the entire array. `signals.h5` uses a compact
concatenated-locus layout:

| Dataset | Shape/dtype | Meaning |
| --- | --- | --- |
| `counts` | `(S, total_locus_bases)`, `uint32` | One row per manifest sample |
| `ip_pooled` | `(total_locus_bases,)`, `uint32` | Chunkwise sum of all IP rows |
| `sample_names` / `sample_roles` | `(S,)`, UTF-8 | HDF5 row identity |
| `transcript_ids` | `(T,)`, UTF-8 | Transcript order |
| `transcript_offsets` / `transcript_lengths` | `(T,)`, `int64` | Slice boundaries in concatenated space |

Users normally do not need to calculate flat offsets. The lazy reader handles
FASTA, HDF5, sample order, and exon blocks:

```python
from transcriptml.rbpnet import ProcessedECLIPDataset

with ProcessedECLIPDataset("processed/chr21") as ds:
    print(ds.transcripts)
    print(ds.samples)
    print(ds.coordinate_space)
    seq = ds.get_sequence("ENST...", 100, 400)
    input_profile = ds.get_profile("ENST...", 100, 400, sample="sminput")
    ip1_profile = ds.get_profile("ENST...", 100, 400, sample="ip1")
    pooled = ds.get_pooled_ip_profile("ENST...", 100, 400)
    blocks = ds.get_genomic_blocks("ENST...", 100, 400)
    genomic = ds.coordinate_to_genome("ENST...", 100)
    coordinate = ds.genome_to_coordinate("ENST...", "chr21", genomic[1])
```

## 2. Descriptive window scanning

The scanner reads the canonical experiment, not BAMs. A new window size or
stride therefore requires only a new scan:

```bash
transcriptml rbpnet scan-windows \
  --processed-dir processed/chr21 \
  --window-size 100 \
  --stride 50 \
  --min-sminput-tpm 0 \
  --pseudocount 1 \
  --output-prefix processed/chr21_windows_100nt
```

It writes equivalent `*.tsv.gz` and `*.parquet` tables plus `*.scan.json`.
Incomplete terminal windows are omitted unless
`--include-incomplete-terminal-windows` is given. Per-transcript cumulative
sums make count aggregation efficient even for stride-1 scans.

Rows include transcript/genomic coordinates and exon blocks; exact overlap
counts/fractions for every region class; GC; SMInput TPM; dynamic per-sample
counts, CPM, and maximum positional counts; pooled-IP values; combined
coverage; and a CPM-scale pooled-IP/SMInput log ratio. Boundary-crossing
windows are `mixed`, never silently assigned to one region.

The scanner is deliberately descriptive. It does not label peaks, negatives,
or training examples.

## 3. Region selection

Selection asks which experimental loci are eligible and why. It writes a
versioned `*.parquet` manifest, equivalent `*.tsv.gz`, and a
`*.selection.json` provenance sidecar.

Every selector can restrict its candidate universe by scanner annotation
without regenerating the descriptive window table:

```bash
transcriptml rbpnet select-regions \
  --processed-dir processed/chr21 \
  --windows processed/chr21_windows_100nt.parquet \
  --strategy peak_gray_negative \
  --region-types 5putr,cds,3putr \
  --output-prefix processed/chr21_exonic_selection
```

Valid values are the biological classes `5putr`, `cds`, `3putr`,
`noncoding_exon`, and `intron`. `mixed` is an annotation status, not a region
class, and therefore is not a valid value for `--region-types`.

By default, a requested class includes both pure windows and mixed windows
with a positive overlap with that class. Thus `--region-types 3putr` includes
windows wholly contained in 3' UTR and CDS/3' UTR boundary windows. Two flags
control boundary-crossing windows:

- `--discard-mixed` retains only pure windows of the requested classes.
- `--only-mixed` retains only boundary-crossing windows overlapping the
  requested classes. It requires `--region-types`.

The flags are mutually exclusive. Omitting `--region-types` preserves the
complete window universe; `--discard-mixed` may still be used to remove all
boundary-crossing windows.

The overlap restriction and mixed-window policy are applied before each
strategy's signal/statistical rules. Consequently, excluded windows do not affect coverage eligibility, the
original selector's testing/50-nt advance, or the `peak_gray_negative` BH
correction universe. The published IP locus-density Poisson null remains based
on the complete locus; this option restricts which windows are tested, not how
that published null is defined. The selection provenance records the requested
types and source/eligible window counts, while the sidecar also reports
selected-example counts by region type.

### Published RBPNet v1

First make the published 100-nt, stride-1 descriptive scan, then select:

```bash
transcriptml rbpnet scan-windows \
  --processed-dir processed/chr21 --window-size 100 --stride 1 \
  --output-prefix processed/chr21_windows_v1

transcriptml rbpnet select-regions \
  --processed-dir processed/chr21 \
  --windows processed/chr21_windows_v1.parquet \
  --strategy original_rbpnet \
  --poisson-null ip_locus_density \
  --output-prefix processed/chr21_original
```

The default preset follows [Horlacher et al. 2023](https://doi.org/10.1186/s13059-023-03015-7)
and its [reference implementation](https://github.com/mhorlacher/rbpnet): a
one-sided Poisson test against the whole-locus pooled-IP rate. For a window of
length `W` in a locus of length `T`, the published/default null is:

```text
mu = pooled_IP_locus_total / T * W
```

A candidate requires uncorrected `p < 0.01`, at least 8 pooled counts, and a
maximum positional count of at least 2. After accepting one, scanning advances
50 nt. No multiple-testing correction is applied: this is the intentionally
lenient published candidate generator, not a calibrated peak caller. The
selected interval remains 100 nt; later model context is independent.

An explicitly experimental alternative uses the matched input:

```bash
transcriptml rbpnet select-regions \
  --processed-dir processed/chr21 \
  --windows processed/chr21_windows_v1.parquet \
  --strategy original_rbpnet \
  --poisson-null sminput \
  --sminput-poisson-pseudocount 1 \
  --output-prefix processed/chr21_original_sminput_null
```

Its exact expectation is:

```text
mu = (SMInput_window_count + 1)
     * (pooled_IP_effective_library_size / SMInput_effective_library_size)
```

The add-one term is applied in SMInput count space and then exposure-scaled to
the pooled-IP library. This asks whether the observed pooled-IP count exceeds
the background expected at the IP sequencing depth. It retains the same
uncorrected `p < 0.01`, count/height requirements, and 50-nt advance. The
manifest records `poisson_null`, the pseudocount, formula, and each selected
window's `selection_null_mean`. This mode is not attributed to Horlacher et al.

### Broad measured windows (`broad_coverage`)

This strategy applies coverage thresholds without a peak test:

```bash
transcriptml rbpnet select-regions \
  --processed-dir processed/chr21 \
  --windows processed/chr21_windows_100nt.parquet \
  --strategy broad_coverage \
  --min-total-count 6 --min-sminput-count 0 --min-ip-count 0 \
  --replicate-mode per_ip \
  --output-prefix processed/chr21_measured
```

The defaults implement `IP + SMInput > 5`, including cases where either track
is zero. The default `--replicate-mode per_ip` applies that criterion independently as
`IP_replicate_j + SMInput >= 6` and emits a replicate-identified row.
`combined` instead applies it to pooled IP and emits one row while preserving
every replicate column. This selector is inspired by the broad-coverage
training philosophy of Yeo et al.; it works on arbitrary TranscriptML scan
tables and is **not** an exact reproduction of their non-overlapping,
annotation-aware Skipper window generation. A locked exact-window preset can
be added separately later.

### Peak / gray / confident negative

```bash
transcriptml rbpnet select-regions \
  --processed-dir processed/chr21 \
  --windows processed/chr21_windows_100nt.parquet \
  --strategy peak_gray_negative \
  --min-total-count 8 --min-sminput-count 0 --min-ip-count 0 \
  --peak-fdr 0.05 --peak-min-log2-ratio 1 \
  --negative-fdr 0.05 --negative-max-log2-ratio -0.5 \
  --stitch-gap 0 \
  --output-prefix processed/chr21_peak_gray_negative
```

By default, combined pooled-IP + SMInput count of at least 8 is sufficient to
enter classification; neither individual track must be nonzero. Optional
`--min-ip-count` and `--min-sminput-count` knobs remain available. Thus both
`IP=0, SMInput=N` and `IP=N, SMInput=0` are valid, informative cases.
Adequately measured windows are tested by conditioning on pooled-IP + SMInput
counts. The null IP probability is determined by effective library sizes.
One-sided exact binomial enrichment/depletion p-values naturally handle both
extremes and are Benjamini–Hochberg corrected. The log2 effect uses the
scanner's explicit CPM-scale pseudocount, so it remains finite at zero. Peaks
require enrichment plus a minimum log2 effect; confident negatives require
depletion plus a maximum log2 effect; the remaining adequate windows are gray.
Low-total-information windows are omitted.
Overlapping/nearby windows are stitched only when transcript, state, and region
type agree. Peak anchors are pooled-signal maxima; negative and gray anchors
are interval midpoints. These defaults are transparent starting choices, not a
definitive CLIP peak caller.

### Selection manifest v1

Every row has a deterministic content-derived `example_id`; gene, transcript,
chromosome, strand, coordinate space, anchor, and half-open selection interval; region overlap;
strategy/state and optional replicate identity; sample/pooled signal summaries;
statistical fields; and explicit gene/transcript/chromosome grouping columns.
Parquet metadata and the sidecar preserve the complete selection and scan
configuration. Consumers must key by `example_id`, not row order.

## 4. Materialized RBPNet bundle

Bundle construction separates the selected biological interval from future
model context:

```bash
transcriptml rbpnet make-bundle \
  --processed-dir processed/chr21 \
  --selection-manifest processed/chr21_measured.parquet \
  --output-dir data/rbpnet_chr21 \
  --input-length 300 \
  --profile-length 300 \
  --max-jitter 32 \
  --transcript-end-policy shift_to_fit
```

For requested input length `L`, profile length `P`, and maximum future jitter
`J`, the bundle writes:

| File | Shape/dtype |
| --- | --- |
| `X.npy` | `(N, 4, L + 2J)`, `uint8`, TranscriptML RNA4 A/C/G/U encoding |
| `sminput_profiles.npy` | `(N, P + 2J)`, `uint32` |
| `ip_profiles.npy` | `(N, R, P + 2J)`, `uint32`; `R` follows `ip_axis_order` |
| `sequence_valid_mask.npy` | `(N, L + 2J)`, `uint8` |
| `profile_valid_mask.npy` | `(N, P + 2J)`, `uint8` |
| `profile_sminput_totals.npy` | `(N,)`, `uint64` |
| `profile_ip_totals.npy` | `(N, R)`, `uint64` |
| `selection_sminput_counts.npy` | `(N,)`, `uint64` |
| `selection_ip_counts.npy` | `(N, R)`, `uint64` |

Pooled IP is `ip_profiles.sum(axis=1)` and is not the only stored
representation. The ordinary TranscriptML sidecars (`ids.txt`,
`metadata.json`, `schema.json`, `config.json`) accompany the arrays, and the
full selected-example metadata is copied to `examples.parquet`. All arrays can
be loaded with NumPy `mmap_mode="r"`.

Boundary handling is configurable:

- `shift_to_fit` is the default. It first centers the full stored width around
  the anchor, then shifts the interval right or left until it lies entirely in
  the locus. The anchor need not remain centered, but every stored position is
  real sequence/signal. A locus shorter than either requested stored width is
  genuinely insufficient; its examples are dropped and counted as
  `n_dropped_short_loci`, never silently padded.
- `drop` preserves the old behavior of discarding any example whose centered
  stored sequence or profile crosses a boundary.
- `pad` preserves fixed widths with all-zero sequence/profile padding. Validity
  masks distinguish padding from ambiguous sequence or true zero signal.

Each example's metadata records `locus_length`, biological anchor,
`sequence_materialized_start/end`, `profile_materialized_start/end`, anchor
offsets, padding amounts, and the crop offsets at both jitter extremes.

For `shift_to_fit`, a future jitter shift `s` in `[-J,+J]` must be resolved in
biological coordinates—not universally as `J+s`:

```text
desired_crop_start = anchor - crop_length//2 + s
actual_crop_start  = clip(desired_crop_start, 0, locus_length-crop_length)
crop_offset        = actual_crop_start - materialized_start
```

Near a boundary, multiple requested shifts can collapse to the same closest
legal crop. The shifted stored interval of width `L+2J` contains every such
legal `L`-nt crop when the locus is long enough. Sequence and profile starts
are recorded separately because their requested lengths may differ. The helper
`transcriptml.rbpnet.bundle.jitter_crop_offset` implements this formula.
Legacy `pad` bundles retain their centered `J+s` crop convention. Bundle
construction does not itself perform random augmentation. Candidate-scan
advance and training-time jitter remain unrelated.

```python
from transcriptml.rbpnet.bundle import load_rbpnet_bundle

bundle = load_rbpnet_bundle("data/rbpnet_chr21", mmap_mode="r")
X = bundle.X
ip = bundle.arrays["ip_profiles"]
input_profile = bundle.arrays["sminput_profiles"]
print(bundle.config["sample_metadata"]["ip_axis_order"])
```

The canonical experiment uses HDF5 because it provides compressed lazy slicing
over an entire transcriptome. The model bundle uses separate `.npy` files
because its selected fixed-shape arrays are simple to inspect and memory-map.
Changing selection, context, or jitter does not require reprocessing BAMs.

(rbpnet-model-and-training)=
## 5. RBPNet model and training

Create a native starter config and train it with the same TranscriptML command
used by other registered models:

```bash
transcriptml init-run --workflow rbpnet --out-dir configs/rbpnet
# Edit dataset/output paths and profile_length, then:
transcriptml train configs/rbpnet/train_config.json

transcriptml evaluate \
  --checkpoint runs/rbpnet/model/best.pt \
  --dataset data/rbpnet_chr21 \
  --out-dir runs/rbpnet/evaluation \
  --split test \
  --save-profiles
```

`transcriptml models show rbpnet --json` prints every architectural default.
The first model family intentionally requires equal sequence and profile crop
lengths. It uses no observed-control input, pooling, reverse-complement
augmentation, valid convolution, absolute-count head, or larger sequence
context than prediction context.

### Default architecture

The RNA4 sequence alone enters a same-padded 1D convolution with 128 filters
and kernel 12 followed by ReLU, then five residual blocks. Each block is a dilated
kernel-6 convolution, BatchNorm, ReLU, dropout 0.25, and residual addition. The
dilations are `[2, 4, 8, 16, 32]`; no positional pooling occurs. Two independent
kernel-25, stride-one transposed-convolution heads produce target and control
positional logits. A global-average-pooled linear head produces the scalar
mixture logit. Important dimensions, normalization, biases, kernels, dilation
schedule, dropout, head kinds, and profile length are configurable.

The trunk receptive field is reported in checkpoints and `summary.json` and is

```text
RF = 1 + (initial_kernel - 1)
       + sum((residual_kernel - 1) * dilation)
```

which is 322 bases for the defaults (160 indexed positions to the left and 161
to the right under the documented asymmetric even-kernel padding). Explicit
left/right padding preserves output index `i` as index `i`; the extra base of
an even effective kernel is placed on the right.

| Model parameter | Default | Meaning |
| --- | --- | --- |
| `in_ch` | `4` | RNA4 input channels. |
| `n_filters` | `128` | Shared positional hidden width. |
| `initial_kernel_size` | `12` | Initial same-padded convolution kernel. |
| `n_residual_blocks` | `5` | Number of residual convolutions. |
| `residual_kernel_size` | `6` | Residual convolution kernel. |
| `dilations` | `null` | Explicit schedule; `null` resolves powers of two starting at 2. |
| `normalization` | `batch` | `batch`, `layer`, or `none`. |
| `dropout` | `0.25` | Dropout in each residual branch. |
| `profile_head_type` | `transpose_conv` | `transpose_conv` or ordinary same-padded `conv`. |
| `profile_head_kernel_size` | `25` | Kernel shared by the two separately parameterized profile heads. |
| `profile_head_bias` | `true` | Whether profile heads include a bias. |
| `enrichment_head_type` | `none` | `none`, `linear`, or `mlp`. |
| `enrichment_hidden` | `64` | Hidden width for the optional MLP only. |
| `enrichment_dropout` | `0` | Optional MLP dropout. |
| `profile_length` | `300` | Required input/output crop length, or `null` to accept any length. |

### Profile model: target, control, and pi

The two heads define independently normalized distributions
`p_target=softmax(target_logits)` and
`p_control=softmax(control_logits)`. With global mixing logit `a`,
`pi=sigmoid(a)` and the predicted IP distribution is

```text
p_IP = pi * p_target + (1 - pi) * p_control
```

The mixture is evaluated with `logsigmoid` and `logaddexp` for stability. `pi`
is the latent fraction of the **positional IP profile** assigned to the target
component. It is not IP/SMInput enrichment and is never given an IP-vs-SMInput
binomial loss. SMInput is a profile training target, not a neural-network
input.

By default, individual IP profiles are summed once across their replicate axis
and the complete multinomial NLL is calculated for the pooled IP counts under
`p_IP`. A second complete multinomial NLL compares SMInput counts with
`p_control`. The `lgamma` combinatorial constant is included by default and can
be disabled. A zero-total profile has no positional information, so that locus
is excluded from that profile component's mean instead of producing a NaN.
Components are reduced over informative loci and weighted by
`lambda_ip_profile` and `lambda_sm_profile` (both 1 by default).

### Optional enrichment model

Set `"enrichment_head_type": "linear"` to enable the default enrichment head,
or `"mlp"` for a configurable two-layer head. The head average-pools the shared
hidden representation only over the biological selection interval, using the
coordinate-derived mask for the current jittered crop. It supports variable
measurement widths. The result `eta_i` is a sequence-predicted log enrichment,
independent of `pi`.

For replicate `j`, effective retained-event library sizes supply the known
offset and the observed selection-interval counts supply the binomial data:

```text
depth_offset_j = log(L_IP_j / L_SM)
logit(p_ij)    = eta_i + depth_offset_j
N_ij           = IP_ij + SM_i
IP_ij          ~ Binomial(N_ij, p_ij)
```

The logits-based complete binomial NLL is evaluated independently for each
valid locus-replicate pair and averaged over those pairs. One sequence row and
one `eta_i` therefore use every IP replicate without duplicating the locus.
`IP=0` and `SMInput=0` edge cases are exact and require no pseudocount; only a
pair with both counts zero is excluded because it has no information. Enabling
the head adds `lambda_enrichment * L_enrichment`, with weight 1 by default.

### Jitter-ready structured batches

`RBPNetDataset` memory-maps the bundle arrays and returns sequence, pooled and
individual IP profiles, SMInput profile, exact selection counts, effective
library sizes/depth offsets, selection mask, valid-position masks, coordinates,
and identifiers. With `max_train_jitter=J`, a deterministic RNG keyed by
seed/epoch/example samples a shift from `[-J,+J]` for training. Sequence and all
profiles use the same biological crop. The crop start is derived from anchor,
actual materialized start, and locus bounds, so boundary-shifted contexts do
not incorrectly assume offset `J+s`. Evaluation always uses shift zero.

When enrichment is enabled, every allowed jittered crop must fully contain its
selection interval; invalid bundle/context combinations fail before training.
`max_train_jitter` cannot exceed the materialized bundle margin.

### Splits, optimization, and outputs

RBPNet starter configs use a `group` split on `group_gene_id`, keeping all
overlapping loci from one gene together. Transcript, chromosome, metadata, and
explicit predefined groups are also usable through their metadata columns.
Every non-random split is checked for group overlap. Row-random splitting is
rejected unless `allow_random_window_split=true` explicitly acknowledges the
leakage risk. Replicate-specific selection rows describing an identical locus
are deduplicated by default while retaining the complete replicate axis.

For chromosome cross-validation, create one plan from the final bundle and
reuse it for every training job:

```bash
transcriptml cv create-chromosome-plan \
  --dataset data/rbpnet \
  --group-col group_chromosome \
  --n-folds 5 \
  --output runs/rbpnet/cv5.json

transcriptml train configs/rbpnet/train_config.json \
  --dataset data/rbpnet \
  --cv-plan runs/rbpnet/cv5.json \
  --fold 0 \
  --output-dir runs/rbpnet/fold0/model
```

Chromosomes are sorted by decreasing example count and greedily assigned to
the currently smallest fold group, with deterministic ties. Run `k` uses group
`k` for test, `(k+1) mod N` for validation, and every remaining group for
training. Resolution verifies that the dataset's chromosome membership and
counts still exactly match the immutable, content-hashed plan. Training records
both the plan path and its validated `plan_id` in summaries and checkpoints.

AdamW, Adam, and SGD; plateau, cosine, and step schedulers; clipping; early
stopping; device selection; DataLoader workers; seeds; and mixed precision are
configurable. `history.json` logs total, pooled-IP profile, SMInput profile, and
enrichment losses independently. `best.pt` and `last.pt` retain model, loss,
optimizer, samples, coordinate space, split, receptive-field, and training
provenance. Evaluation CSVs contain `pi`, optional `eta`, and each replicate's
depth-adjusted predicted IP fraction. The raw structured tensors remain
available through the Python model output for future attribution work.

### Scientific evaluation reports

For an RBPNet checkpoint, `transcriptml evaluate` writes a structured report
when given `--out-dir`. The default split is `test`; `--split` accepts `train`,
`val`, `test`, or `all`. These indices are resolved **only** from `splits` stored
in the checkpoint used for training. RBPNet evaluation deliberately does not
fall back to `bundle.splits`, because a reused or edited bundle must not change
which observations are considered held out. Evaluation always uses jitter
shift zero and is deterministic apart from platform-level floating-point
details.

```text
evaluation/
  summary.json
  examples.parquet
  stratified_metrics.parquet
  calibration.parquet
  plots/
  predicted_target_profiles.npy   # only with --save-profiles
  predicted_control_profiles.npy
  predicted_ip_profiles.npy
```

The optional profile arrays are normalized positional probabilities (`float32`)
with shape `(N_evaluated, profile_length)`; each row sums to one over valid
positions. Their first axis is exactly `evaluation_row` in
`examples.parquet`; they can be opened without loading them into memory using
`np.load(path, mmap_mode="r")`. `predicted_ip_profiles.npy` stores the final
target/control mixture, not just the latent target component.

For observed counts `y`, total `N`, empirical distribution `q=y/N`, predicted
distribution `p`, uniform distribution `u` over valid positions, and predicted
control distribution `p_control`, profile metrics use natural logarithms:

| Metric | Definition |
| --- | --- |
| complete multinomial NLL | `-log Multinomial(y | N, p)`, including the count combinatorial constant |
| KL/read (saturated gap) | `(NLL_model - NLL_saturated) / N = KL(q || p)` |
| JSD | `0.5 KL(q || m) + 0.5 KL(p || m)`, where `m=(q+p)/2` |
| information gain over uniform/read | `(LL_model - LL_uniform) / N` |
| pooled-IP information gain over control/read | `(LL_predicted_IP - LL_predicted_control) / N` |
| Wasserstein | one-dimensional earth-mover distance between `q` and `p`, in nucleotides |

Empirical-profile metrics are undefined and recorded as null/NaN when the
observed profile total is zero. Counts outside a validity mask are rejected;
all probability distributions are restricted and normalized over valid
positions. The report always exposes the scientifically comparable complete
NLL, while its checkpoint-objective reconstruction honors the checkpoint's
`include_multinomial_constant` and `include_binomial_constant` settings.

With the enrichment head enabled, the report retains the complete
replicate-aware binomial NLL and compares it with the depth-only null
`eta=0`. Its information gain is
`(LL_model - LL_eta=0) / (IP+SMInput)`. The descriptive empirical enrichment is

```text
eta_hat = log((IP + c) / (SMInput + c)) - log(L_IP / L_SM)
```

where `c=0.5` by default and is configurable with
`--enrichment-pseudocount`. This pseudocount is used only for Pearson/Spearman
diagnostics, never for either likelihood. `calibration.parquet` contains both
locus rows and fixed-width predicted-probability bins. Each bin uses
read-weighted summaries:

```text
predicted_bin = sum(N * p) / sum(N)
observed_bin  = sum(IP) / sum(N)
```

With multiple IP replicates, each observed replicate is also compared with the
pooled profile of all other replicates. The resulting leave-one-replicate-out
JSD and Wasserstein values are an experimental reproducibility reference, not
a guaranteed upper bound on every model metric.

`stratified_metrics.parquet` is long-form. It reports locus macro (replicates
within a locus averaged first), gene macro (loci within a gene averaged first),
and read micro summaries. Read micro likelihood and information values sum the
appropriate likelihood/information numerators and divide by contributing
reads; JSD and Wasserstein use read-count-weighted means. Every row records its
informative observation, locus, gene, and read counts. Summaries include
pooled-IP and SMInput read-depth bins, chromosome, and optional
`selection_state` and `region_type` strata.

`examples.parquet` is the inspectable per-locus table. It includes bundle and
evaluation indices, stable example and biological identifiers, selection
metadata where available, profile and selection-interval counts, effective
library sizes/depth offsets, `pi`, optional `eta` and replicate predicted IP
fractions, every per-example profile metric, and replicate-ceiling metrics.
`summary.json` records metric definitions, sample depths, aggregate metrics,
correlations, representative-example sampling, plot provenance, and all output
paths.

The plot collection includes performance versus depth, metric distributions,
eta agreement, read-weighted calibration, fixed-seed representative profiles
(default: at least 10 pooled-IP and 10 SMInput profile reads),
selection/region stratification, and `pi` diagnostics. Missing optional
metadata, a disabled enrichment head, or a single IP replicate causes only the
inapplicable plot/metric to be skipped; `summary.json` records why. The legacy
`--out-csv` route remains available for compact `pi`/`eta` predictions, but the
structured report is preferred for scientific evaluation.

Profile-only model block:

```json
{
  "model": {
    "name": "rbpnet",
    "params": {"profile_length": 300, "enrichment_head_type": "none"}
  },
  "loss": {"name": "rbpnet"},
  "max_train_jitter": 0
}
```

To train profiles plus enrichment, change only the head and, if desired, the
independent component weights:

```json
{
  "model": {
    "name": "rbpnet",
    "params": {"profile_length": 300, "enrichment_head_type": "linear"}
  },
  "loss": {
    "name": "rbpnet",
    "lambda_ip_profile": 1.0,
    "lambda_sm_profile": 1.0,
    "lambda_enrichment": 1.0
  },
  "max_train_jitter": 32
}
```