Metadata-Version: 2.4
Name: ebin
Version: 1.2.0
Summary: Posterior bin-expectation normalization for sorting-based MPRA
Requires-Python: >=3.9
Description-Content-Type: text/markdown
Requires-Dist: jax>=0.4
Requires-Dist: numpy
Requires-Dist: scipy
Requires-Dist: pandas
Provides-Extra: netcdf
Requires-Dist: xarray; extra == "netcdf"
Requires-Dist: h5netcdf; extra == "netcdf"

# EBin: Normalization of MPRA data (FACS bins × replicates × sequences)

Raw per-bin read counts are turned into one number per sequence per cell line —
the **posterior-mean expected bin**, `E[bin]` — together with the posterior SD
that states how far it can be trusted. The counts enter a compound NB–Poisson
likelihood; what is estimated is the sequence's latent expression distribution.



<img src="http://data.georgy.top/media/ebin_model_layers.png" alt="EBin model"/>


```MPRA 
X[n,s,b] | T[n,s,b] ~ NB(R[s,b] · T[n,s,b], P[s,b])       reads from cells
T[n,s,b]            ~ Poisson(a_n · Pi[n,b] · lambda_s)   cells in bin b
Pi[n,b]             = Phi((q_b − mu_n)/sigma_n) − Phi((q_{b−1} − mu_n)/sigma_n)
```
Each sequence `n` carries a latent Gaussian effect law `N(mu_n, sigma_n²)` — its
expression distribution over cells — and an abundance `a_n`. The library's
marginal effect law is the **`N`-component mixture** of those per-sequence
normals, `G(x) = Σ_n w_n Phi((x − mu_n)/sigma_n)` with `w_n = 1/N`, and the `B`
sorting bins are its quantile cells: the cut points `q_1 < … < q_{B−1}` solve
`G(q_j) = j/B`. They are shared, and re-solved at every step, so every sequence
is tied to the marginal through them. This approach has 2 benefits: first one is that
it allows us to cut the number of params down to 2 instead of `B-1`, and the second one is that those params are further regularized, both by the normality assumption and the fact that they are all should form the global marginal distribution. 

The activity that ships is the posterior **mean** of `E[bin] = Σ_b b·Pi[n,b]`,
integrated over the posterior of `(mu_n, sigma_n)` with the abundance profiled
per grid point. The plug-in point estimate is a posterior *mode*, and for a
low-count sequence the mode sits at a simplex vertex — activity exactly 1 or
`B`, regardless of depth.

## Install
```
pip install ebin
```
or

```bash
pip install -e /path/to/package        # or place this directory on PYTHONPATH
```

`jax` (GPU-version strongly recommended: `pip install jax[cuda]`), `numpy`, `scipy` and `pandas` are required.
`xarray` + `h5netcdf` are needed only for the netCDF export and `matplotlib`
only for the diagnostics. Everything runs in float64 as some parts of the algorithm require high precision, namely the tridiagonal eigenvalue solver for orthogonal polynomials.

## Input

One row per sequence, a 3-level column header `(cell line, replicate, bin)`:

```
seq,       A549,A549,A549,A549,A549,...
           rep1,rep1,rep1,rep1,rep2,...
           bin1,bin2,bin3,bin4,bin1,...
ACGT...,   112, 340, 88,  15,  97,  ...
```
(in other words, the tabular should be readable with `pd.read_csv(..., header=[0,1,2])`)


Bins are ordered by the number in their name. Cell lines are fitted
independently. **Empty fields are read as zero counts, not as missing data** —
some exports write every zero as a blank, and treating those as missin

## Use

```bash
python -m ebin counts.csv -o activity.csv
```

```python
from ebin import activity_table

table, results = activity_table("counts.csv", out="activity.csv")
table[("A549", "activity")]        # posterior-mean E[bin] per sequence
```

The CSV has two header rows, `(cell line, field)`, one row per sequence, and
carries everything the model estimates per sequence:

| field | meaning |
| --- | --- |
| `activity` | posterior-mean `E[bin]` — the normalized value intended for downstream use |
| `activity_sd` | its posterior SD — per-sequence confidence |
| `activity_map` | plug-in mode, kept for reference only |
| `mu`, `mu_sd` | posterior-mean effect location and its standard error |
| `sigma`, `sigma_sd` | posterior-mean effect scale and its standard error |
| `abundance` | fitted size factor `a_n` |
| `n_eff` | posterior grid support — large means the estimate is diffuse |
| `tot` | total observed reads |

All of it comes from the same posterior weights, so the columns are mutually
consistent: `activity` is `Σ_b b·bin_prob`. `mu` and `sigma` are in gauge units
(median cut 0, first cut −1), so they are comparable across sequences within a
cell line. A subset is selected with `--fields activity,activity_sd,tot` (or
`fields=[...]` in the API).

The `B−1` cut points are a property of the cell line, not of a sequence, so
they are not columns: `cuts_table(results)` returns them as one row per cell
line, and the CLI writes them next to the activity table as
`<out>_cuts.csv`. The gauge fixes the first two at −1 and 0, so what varies
between lines is where the upper cuts land (A549 0.990, MDA-MB-231 0.964).

Sequences with no reads at all in a cell line keep an activity (they are fitted
as structural zeros and land near the neutral prior mean, ~2.5 for `B = 4`);
sequences with no observed cell at all are `NaN`.

Useful CLI options:

```
--groups A549,HepG2        cell lines to run, in output order
--abundance gamma          integrate a_n out instead of fitting it
--zero-truncation          drop all-zero rows instead of fitting them
--grid 61x25               readout grid (n_mu x n_log_sigma)
--a-max 20                 pin the abundance ceiling (default: auto, see below)
--a-max-tol 0.001          sequences allowed to sit at the ceiling before it doubles
--tau-buckets 6            latent-count grid lengths (1 = one grid for everybody)
--save-fits fits.pkl       keep the fitted parameters
--warm-from fits.pkl       start from an earlier fit of the same data
--netcdf activity.nc       also write the full posterior (see below)
--plots DIR                also write the diagnostic plots
```

### Per-bin distribution

The one field too wide for a flat table is `bin_prob`, the posterior-mean
probability of each bin — an `(N, B)` array per cell line. It is returned in
`results`, and `--netcdf` writes it alongside every scalar field:

```python
from ebin import activity_table, write_netcdf

table, results = activity_table("counts.csv")
results["A549"]["bin_prob"]                    # (N, B), rows sum to 1
write_netcdf("activity.nc", results, table.index)
```

Readouts can be redone from saved fits without refitting, which is the cheap
way to try a different grid or prior:

```python
from ebin import readout_table
table, results = readout_table("counts.csv", "fits.pkl")
```

## Diagnostics

Three plots, one panel per cell line, drawn from the table alone (`--plots DIR`
on the CLI writes all three):

```python
from ebin import (plot_marginal_distribution, plot_activity_vs_raw,
                  plot_gc_vs_abundance)

plot_marginal_distribution(table, out="marginal.png")
plot_activity_vs_raw(table, "counts.csv", out="vs_raw.png")
plot_gc_vs_abundance(table, out="gc.png")
```

**`plot_marginal_distribution`** — the marginal effect law: the `N`-component
mixture `G(x) = Σ_n w_n Phi((x − mu_n)/sigma_n)`, one normal per sequence,
drawn as the density it defines and sliced into the sorting bins, each in its
own colour. The x ticks are the cut points with their values. Weights are
uniform by default, as in the fit; `weights=` takes a per-sequence weight (the
fitted abundances, say) to see the marginal over cells rather than over
sequences. Every slice holds exactly `1/B` of the mass by construction, so the
gauge (cuts at −1 and 0) and any unusual tightness or skew in a cell line's
effect law are visible here.

**`plot_activity_vs_raw`** — fitted activity against the depth-normalized raw
mass centre (each replicate×bin column divided by its own library size, then
`Σ_b b·p_b`), coloured by fitted abundance. Well-measured sequences sit on the
diagonal; the horizontal band at the centre is the low-abundance tail, where the
raw estimate is noise and the model shrinks toward neutral. A cell line whose
*high*-abundance sequences leave the diagonal is a real warning sign.

**`plot_gc_vs_abundance`** — GC content against the fitted size factor, with the
binned median trend. A strong trend means `a_n` is absorbing a
cloning/sequencing bias rather than biological abundance alone. GC is taken from
the table index when the index is the sequence; otherwise it is supplied through
`sequences=`.

`raw_mass_center(counts)` and `gc_content(seqs)` are exported separately for
cases where the raw numbers are more appropriate than a picture.

## Configuration

Defaults are the shipped configuration; they are not arbitrary and are worth
understanding before they are changed.

| option | default | why |
| --- | --- | --- |
| `conditional` | `False` | Zero **inflation**: all-zero rows are fitted with a structural-zero probability `phi`, so they enter the mixture that sets the cuts. Zero **truncation** (`True`) drops them, which sets the cuts on a zero-excluded distribution — measurably worse, both in-library and on transfer to an independent library. |
| `abundance` | `True` | A free size factor `a_n` per sequence. Without one, every sequence has the same expected latent total, so real depth variation has nowhere to go but `Pi` — which tilts toward the bins with the largest channel scale, and `Pi` stops being a profile. |
| `abundance_prior` | `None` | `'gamma'` integrates `a_n` out under `Gamma(kappa, mean=1)` instead: 2 parameters per sequence instead of 3, low-count sequences regularized, and the abundance uncertainty propagated into the effect. The activity profile is essentially unchanged (ρ ≈ 0.9996 against free `a_n`), so this is a parameter-economy. Costs ~2× runtime. |
| `mu_prior` | `1.3` | SD of the prior on the effect location. Bounding the effect *spread* alone still lets a single-bin low-count sequence run `mu_n` large and walk `Pi` to a vertex. The gauge pins Q50 = 0 and Q25 = −1, so the population SD is ≈ 1.5 and the between-sequence SD ≈ 1.1 — the gauge-consistent `tau` is ~1.3. A larger value (2–3) is inconsistent with the gauge and under-regularizes. |
| `sigma_prior` | `(0.0, 0.5)` | Prior on `log sigma_n`. The per-sequence MLE is unbounded: a sequence whose reads land only in the outer bins wants a two-point mass, reachable only as `sigma → ∞`. Without this, `sigma` hits its bound in several cell lines; those are bad optima, so the prior often *improves* the data likelihood. `(0, 1)` fits marginally better on well-behaved lines but leaves stuck ones stuck. |
| `lambda_fix` | `50.0` | The latent Poisson level is not identified upward — it rides the `R·lambda` ridge, and only products like `R(1−P)/P · Pi · lambda` are identified. It is held fixed and `R` absorbs the scale. |
| `a_max` | `"auto"` | Ceiling on `a_n`. A fixed `15` was too low: it pinned 0.9–1.4 % of every lib2 cell line at the bound, and moving those sequences off it shifts the whole line's activity (median 0.03, up to 0.5), because the cuts are shared. `"auto"` starts at 20 and doubles the bound of whatever sits on it until under `bound_frac` do — the bound is per sequence, so only the deep ones ever carry a large one. A number pins it and reproduces the old behaviour. |
| `bound_frac` | `0.001` | When `a_max="auto"`: the share of sequences allowed to remain at the bound. An inactive constraint does not move an optimum, so this is how close to "no ceiling at all" the fit is required to get. 0.5 % is measurably too loose — see below. |
| `cap_headroom` | `3.0` | Where each sequence's bound starts, as a multiple of its own depth-derived abundance. Smaller is faster (it sizes that sequence's latent-count grid) and anything too small is caught and doubled, so this only trades one extra warm refit against grid length. |
| `tau_buckets` | `6` | Distinct latent-count grid lengths the sequences are bucketed onto. `1` puts every sequence on a grid sized by the deepest one, which is what the model used to do: on lib1 that is 3× the arithmetic and a 2.1× slower gradient at `a_max = 20`, and 14× the arithmetic at `a_max = 100`, where it no longer fits in 16 GB. Each bucket is one more XLA kernel to compile, so this is a padding-vs-compile knob and 4 is nearly as good as 8. |
| readout grid | `121 × 33` | `(mu, log sigma)` grid for the posterior mean. Converged in the bulk: doubling to `181 × 49` moves the median activity by 2e-6 (a handful of extreme low-count sequences move by up to 0.04), while halving to `61 × 25` moves it by 5e-3. |

The readout re-uses the fit's own priors by default, so the posterior mean is
coherent with the fitted model. A flat effect location is obtained with
`readout=dict(tau_within=None)`, and `mu_center=<array>` shrinks toward a
per-sequence target (for instance the same sequence's consensus across cell
lines).

## Runtime

One cell line of 30k sequences × 4 bins takes ~2–4 min on a GPU (fit + readout)
and roughly 20× longer on CPU; the Gamma-abundance variant costs ~2× the
free-`a_n` fit and does not bucket (its abundance nodes are shared across
sequences, so there is nothing per-sequence to size a grid from). Memory is
dominated by the latent-count grid; the likelihood scans over abundance/rate
components so only one is live at a time. `XLA_PYTHON_CLIENT_PREALLOCATE=false`
is set on import so JAX does not claim 75% of VRAM.


## Parameters


| symbol | shape | role | what it is | example (A549) |
| --- | --- | --- | --- | --- |
| `X[n,s,b]` | N×S×B | observed | reads counted for sequence n in channel (s,b) | 1686 median total |
| `T[n,s,b]` | N×S×B | latent | cells of n that landed in bin b — never observed, summed out | ~63 cells |
| `mu_n`, `sigma_n` | N each | per sequence | the sequence's expression law over cells, in gauge units | σ median 1.46 |
| `a_n` | N | per sequence | size factor: how many cells n contributes at all. Geometric mean pinned to 1, capped at `a_max` | 1.42 median |
| `Pi[n,b]` | N×B | derived | bin profile, rows sum to 1. Not free — read off `(mu_n, sigma_n)` and the cuts | — |
| `q_j` | B−1 | shared | cut points = quantiles of the mixture `G` at `j/B`. Gauge pins `q_1 = −1`, `q_2 = 0` | −1, 0, 0.99 |
| `lambda_s` | S | fixed | cell level of replicate s. Not identified — rides the `R·lambda` ridge, so held fixed while `R` absorbs the scale | 50 |
| `R[s,b]` | S×B | per channel | NB size per latent cell: `R·T` is the shape of the read distribution | 12.7 11.4 11.6 20.1 |
| `P[s,b]` | S×B | per channel | NB probability. Small `P` = a noisier channel at the same depth | .33 .32 .31 .52 |
| `phi` | 1 | shared | probability a sequence is a structural zero (dropout, never a cell) | 8e−7 |
| `kappa` | 1 | optional | Gamma shape when `a_n` is integrated out instead of fitted; CV = 1/√κ | 0.9 (variant) |

so `X, T, R, P` are 3D arrays:
<center>
<img src="http://data.georgy.top/media/ebin_index_cube.png" alt="Index cube: n sequences, s replicates, b bins" width="196">
</center>

