Metadata-Version: 2.5
Name: fertilizer-genomics
Version: 0.2.0
Summary: Identify fertile ground in genomes for minimal-edit regulatory design.
Project-URL: Homepage, https://github.com/jmschrei/fertilizer
Project-URL: Issues, https://github.com/jmschrei/fertilizer/issues
Author: Jacob Schreiber
License: MIT License
        
        Copyright (c) 2026 Jacob Schreiber
        
        Permission is hereby granted, free of charge, to any person obtaining a copy
        of this software and associated documentation files (the "Software"), to deal
        in the Software without restriction, including without limitation the rights
        to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the Software is
        furnished to do so, subject to the following conditions:
        
        The above copyright notice and this permission notice shall be included in all
        copies or substantial portions of the Software.
        
        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
        LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
        OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
        SOFTWARE.
License-File: LICENSE
Keywords: atac-seq,bam,bigwig,bioinformatics,chip-seq,deseq2,enrichment,fragments,genomics
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Typing :: Typed
Requires-Python: >=3.10
Requires-Dist: joblib>=1.3
Requires-Dist: numpy>=1.24
Requires-Dist: pandas>=2.0
Requires-Dist: pybigwig>=0.3.22
Requires-Dist: pysam>=0.22
Requires-Dist: scipy>=1.10
Provides-Extra: dev
Requires-Dist: pytest>=8.0; extra == 'dev'
Requires-Dist: ruff>=0.5; extra == 'dev'
Description-Content-Type: text/markdown

# fertilizer

**What this does (concretely):** `fertilizer` takes one bigWig, BAM or 10x fragment file per condition (or one fragment file plus a table assigning cell barcodes to conditions) and a set of BED regions, computes a summary statistic or read count per region per condition, and calls regions where one condition has significantly higher signal than the others. Output is a TSV with effect size, p-value, q-value, and the name of the enriched condition. The statistical model is a DESeq2-inspired negative-binomial GLM likelihood-ratio test adapted to the one-replicate-per-condition setting.

**Why "fertilizer".** The **Fertile Ground Hypothesis** is that genomes are full of "almost-regulatory" regions — sequences that do not do anything on their own, but can be minimally edited to achieve subtle and precise activity. Many near-motifs, for example, sit one or two substitutions away from binding a transcription factor and recruiting its downstream regulatory activity. `fertilizer` helps identify the fertile ground in a genome that is most useful for your design task by flagging regions where signal in one condition stands out from the others.

> **Status.** v0.2.0 — API is unstable until 1.0. Please report issues. Most
> of this package (including the test suite) was drafted with Claude
> assistance; the statistical methodology and calibration are exercised by
> the simulation tests in `tests/test_enrichment.py` and we are continuing
> to validate against external benchmarks.

## Installation

`fertilizer` is installable with [uv](https://docs.astral.sh/uv/). The PyPI
distribution is named `fertilizer-genomics` (the `fertilizer` name was taken),
but the importable module is still `fertilizer`.

```bash
# install from PyPI
uv pip install fertilizer-genomics

# or install directly from GitHub
uv pip install git+https://github.com/jmschrei/fertilizer.git

# or install from a local clone
git clone https://github.com/jmschrei/fertilizer.git
cd fertilizer
uv pip install -e .
```

> `pyBigWig` needs `libcurl` and `libssl` headers at install time on Linux.
> If pip fails to build it, install them first:
> `sudo apt-get install libcurl4-openssl-dev libssl-dev zlib1g-dev` (Debian/Ubuntu)
> or `brew install curl openssl` (macOS). `pysam`, used to read BAM files,
> ships prebuilt wheels for Linux and macOS.

For a development environment with test and lint tooling:

```bash
uv venv
uv pip install -e ".[dev]"
pytest tests/
```

## Quickstart — try it on a tiny demo

```bash
python examples/make_demo_data.py
fertilizer extract -w examples/A.bw examples/B.bw examples/C.bw \
    -b examples/regions.bed -o examples/signals.tsv -s sum
fertilizer enrich -i examples/signals.tsv -c A B C \
    -o examples/enrichment.tsv --q-threshold 0.05
```

`examples/C.bw` has signal enriched at 10 randomly chosen regions; the
resulting `enrichment.tsv` should contain ~10 rows, each with
`enriched_condition == C`. See `examples/README.md` for a walkthrough.

## When NOT to use this

`fertilizer` is designed for the **one-sample-per-condition** setting with
many loci, most of which are not differentially enriched. It is the wrong
tool when:

- **You have replicates per condition.** Prefer DESeq2 / edgeR / csaw —
  they estimate per-locus dispersion from within-condition variance and
  give better power. `fertilizer` collapses biological + technical
  variability into a single across-condition α.
- **Most loci genuinely do change between conditions.** Median-of-ratios
  size factors and the across-loci dispersion estimator both assume that
  the null *majority* anchors normalization and variance. Pass a
  background-matched / genome-wide region set, not a candidate-only set.
- **Your experiment causes a global shift in the mark** (e.g. an EZH2
  knockdown that collapses H3K27me3 genome-wide). The null-majority
  assumption is violated by construction; use spike-in normalization
  and pass it via `--size-factors`.
- **Your bigWig aggregates are not count-like.** `enrich` refuses input
  produced by `extract --stat mean / max / min / std / coverage`
  (override at your own risk with `--allow-non-sum`). The NB-GLM
  assumes the variance-mean relationship of count data. If you have the
  BAM or fragment files, count those instead (see
  [BAM and fragment input](#bam-and-fragment-input)).
- **Your regions overlap densely** (sliding/tiling windows). BH controls
  FDR under independence or PRDS; overlapping windows violate this and
  q-values will be optimistic. Thin to non-overlapping regions, or use a
  method that models autocorrelation explicitly.

## Usage

`fertilizer` exposes two analysis subcommands (plus `install-skill`, below). Typical workflow:

```bash
fertilizer extract -w A.bw B.bw C.bw -b regions.bed -o signals.tsv -s sum
fertilizer enrich  -i signals.tsv -c A B C -o enrichment.tsv
```

`extract` also counts reads from BAM files (`-a`) or fragment ends from 10x
fragment files (`-f`), optionally splitting one fragment file into
per-group pseudobulks with a barcode table (`-g`):

```bash
fertilizer extract -a A.bam B.bam C.bam -b regions.bed -o counts.tsv -ps 4 -ns -5
fertilizer extract -f fragments.tsv.gz -g cells.tsv --group-column cluster -b regions.bed -o counts.tsv
```

> The NB-GLM in `enrich` assumes **count-like** input. Use `extract --stat sum`
> (the total signal over each region) for bigWigs, or count BAM/fragment
> input, which is always counted — bigWig `mean`/`max`/`min`/`std`/`coverage`
> are not counts and `enrich` will refuse them unless `--allow-non-sum` is
> passed. `extract` writes a metadata header (`# fertilizer-extract stat=...`)
> that `enrich` reads to enforce this.

### `fertilizer extract` — signal aggregation

For bigWigs, compute a per-region summary statistic (mean by default; `-s` chooses among `mean`/`max`/`min`/`sum`/`std`/`coverage`) for each file over each region in the concatenated BED input. For BAM/SAM/CRAM and fragment files, count reads or fragment ends per region instead ([BAM and fragment input](#bam-and-fragment-input)). One row per region, one column per input file (or per barcode group). Input row order is preserved.

| flag | description |
| --- | --- |
| `-w`, `--bigwigs` | one or more bigWig signal tracks. Exactly one of `-w`, `-a` or `-f` is given per run |
| `-a`, `--bams` | one or more BAM/SAM/CRAM files, counted as described in [BAM and fragment input](#bam-and-fragment-input) |
| `-f`, `--fragments` | one or more 10x fragment files (plain or gzipped), counted as described in [BAM and fragment input](#bam-and-fragment-input) |
| `-b`, `--beds` | one or more BED region files. Columns 1-3 are required (`chrom`/`start`/`end`); columns 4-6 are passed through as `name`/`score`/`strand`; any further columns are passed through as `bed_col_<i>` (BED12 and narrowPeak disagree on the meaning of columns 7+, so generic names are used to avoid mislabeling). `#` comment lines are skipped, as are UCSC `track` and `browser` lines at the top of a file. |
| `-o`, `--output` | path to the output TSV; gzip-compressed when the name ends in `.gz` (the metadata header is kept, and `enrich` reads it from the compressed file) |
| `-s`, `--stat` | bigWig input only. Per-region summary statistic: `mean` (default), `max`, `min`, `sum`, `std`, `coverage`. Maps to pyBigWig's `stats(type=..., exact=True)`, so values come from the full-resolution data rather than the bigWig's zoom levels. **Use `sum` if the output will be passed to `fertilizer enrich`** — the NB-GLM assumes count-like input. `extract` writes a `# fertilizer-extract stat=...` header line so `enrich` can verify this. |
| `-n`, `--names` | optional explicit column names, one per input file. Overrides the default of using each file's name without its extension (`A.bw` → `A`, `C1.fragments.tsv.gz` → `C1.fragments`). Useful when two paths share a basename (e.g. `RNAseq/A.bw` and `ATACseq/A.bw`). A name that matches a BED column present in the input (`chrom`, `start`, `end`, `name`, `score`, `strand`, `bed_col_<i>`) is rejected. |
| `-j`, `--n-jobs` | parallel workers (default `-1`, one per core). On a shared machine set this explicitly. bigWigs use threads, splitting the regions into one chunk per worker per file; BAM/CRAM and fragment input use processes and split within a file ([BAM and fragment input](#bam-and-fragment-input)) |

**Coordinates are 0-based half-open**, matching the standard BED/UCSC bigWig convention. A region `chr1 100 200` covers bases 100..199 inclusive (length 100). If your input is a 1-based file (UCSC table dumps, some BED-like exports), subtract 1 from `start` before running `extract`.

**Partially covered regions.** bigWigs do not store uncovered bases, and every statistic except `coverage` is computed over the covered bases only. For a 200 bp region with 50 bp at 5.0, 50 bp at 1.0 and 100 bp uncovered, `mean` is 3.0 (not 1.5), `min` is 1.0 (not 0.0), `max` is 5.0, `std` is the sample standard deviation of the 100 covered values, `sum` is 300 (uncovered bases add nothing), and `coverage` is 0.5. If your bigWig writes zeros explicitly, those bases are covered and do count. Use `sum` or `coverage` when the uncovered fraction matters.

Zeros never mean "missing" — the output is always numeric, never `NaN`. A region whose summary statistic genuinely resolves to zero (empty bigWig, uncovered span) is reported as `0.0` silently. A region with a locus-level problem (unknown chromosome, coordinates past the end of the chromosome, zero-length interval, negative start) is also reported as `0.0` but triggers a single `FertilizerWarning` — one warning per distinct issue type per run, regardless of how many rows or bigWigs were affected.

Example `signals.tsv`:

```
chrom   start   end     A       B       C
chr1    0       500     2.0     7.0     3.1
chr1    500     1000    6.0     7.0     5.8
chr1    1500    2000    4.0     7.0     4.4
chr2    0       100     0.0     0.0     0.0
```

Column names come from each input file's name without its extension (override with `-n/--names`), so passing two files with the same basename (even from different directories) without `--names` is rejected up front. `extract` also emits a `FertilizerWarning` if more than 95% of the output cells are exactly zero — almost always a wrong path or a chromosome-naming mismatch (e.g. `chr1` in the BED but `1` in the bigWig).

#### BAM and fragment input

With `-a` or `-f`, each region's value is a **count** of positions falling in `[start, end)`:

- **BAM/SAM/CRAM (`-a`):** the 5′ end of each read — the leftmost aligned base of a forward read, the rightmost of a reverse read. CRAM is decoded without read sequences, which the counts never use, so no reference FASTA is needed. Each mate of a pair is counted separately, so for paired-end ATAC-seq this counts both Tn5 insertions of every fragment.
- **Fragment files (`-f`):** both ends of every fragment, i.e. its two Tn5 insertions, at `start` and `end − 1`. The expected layout is 10x's `chrom, start, end, barcode, count`, with optional `#` header lines. Each line counts once; the duplicate count in column 5 is ignored.

| flag | description |
| --- | --- |
| `-ps`, `--pos-shift` | added to each read's or fragment's start coordinate before counting (default 0) |
| `-ns`, `--neg-shift` | added to each read's or fragment's end coordinate before counting (default 0). The two shifts behave exactly as bam2bw's `-ps`/`-ns`; `-ps 4 -ns -5` applies the standard Tn5 offset. 10x fragment files are already shifted, so leave both at 0 for them |
| `--min-mapq` | BAM/CRAM only: skip reads with mapping quality below this (default 30) |
| `--include-flagged` | BAM/CRAM only: count reads carrying these flags, which are skipped by default: `duplicate`, `secondary`, `supplementary`, `qcfail`. Unmapped reads are always skipped |
| `-g`, `--groups` | fragment files only: a tab-separated table with a header row (plain or gzipped) assigning cell barcodes to groups. Writes one column per group, in order of first appearance, summed over every `-f` file; barcodes not in the table are ignored and a barcode listed under two groups is an error. Cannot be combined with `-n` |
| `--barcode-column`, `--group-column` | the `--groups` columns holding barcodes and group labels (defaults `barcode` and `group`) |

Counting reads directly gives `enrich` true counts. A bigWig `sum` over a coverage track is roughly reads × fragment length, so the Poisson part of the NB variance understates sampling noise on that scale, most at low counts.

**Parallelism.** An indexed BAM or CRAM (`.bai`, `.csi` or `.crai`) is split into chromosome pieces, about four per `-j` worker, run in separate processes. Pieces are sized by how many reads the index records on each chromosome (mapped reads for a BAM index, compressed data for a `.crai`), and chromosomes with no reads are skipped, so a single-chromosome file also uses every worker. Each read is counted by the piece where its leftmost aligned base lies. An unindexed BAM or CRAM, or a SAM, is streamed as one task.

Fragment files are streamed in chunks, so memory does not grow with file size and no index is needed. A BGZF-compressed file (what 10x and `bgzip` write) or an uncompressed one is also split into byte ranges so that a single large file uses all `-j` workers: roughly `-j` divided by the number of files ranges per file, each at least 16 MB. Ranges are cut at BGZF block starts and every line is counted by exactly one range. Plain (non-BGZF) gzip cannot be split and is read by one worker. Each worker parses 16 MB of text at a time; on a 3.8 GB fragment file counted into 70 groups, peak memory was about 0.25 GB per worker (8 GB at `-j 32`).

The output header records the input and shifts (`# fertilizer-extract stat=count source=fragments pos_shift=0 neg_shift=0`), and `enrich` accepts it without `--allow-non-sum`. Locus-level problems are reported as for bigWigs, with two differences for fragment files, which carry no chromosome lengths: a chromosome counts as missing when it never appears in the file, and regions past a chromosome's end are not detected.

### `fertilizer enrich` — enrichment analysis

Identify loci where one condition is **enriched** relative to the others. Takes a TSV whose columns include (at minimum) the non-negative, finite numeric columns named via `-c` — the output of `fertilizer extract`, which also carries `chrom`/`start`/`end`, is the canonical input and is passed through verbatim — names the columns to compare, and writes a **filtered** TSV containing only loci that pass the significance threshold, with effect size, p-value, q-value, and the name of the enriched condition. Loci where one condition is *depleted* relative to the others are **not** called. Any columns present in the input beyond the tested conditions are preserved in the output.

The test is a one-sided negative-binomial GLM likelihood-ratio test per locus on `k*` and a **background condition at a user-chosen rank**, inspired by [DESeq2](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8) (Love, Huber, Anders, 2014) and adapted to the **one-replicate-per-condition** setting this package targets. Steps:

1. **Size factors.** DESeq2's median-of-ratios, computed jointly across all conditions from loci with positive signal in every condition.
2. **Dispersion.** With no within-condition replicates, per-locus variance cannot be estimated from repeated observations. A dispersion parameter `α` is estimated *across* loci under the null-majority assumption. Default (`--fit-type common`) uses the median of per-locus method-of-moments estimates, scaled by a closed-form correction for the small-df median-of-χ² bias; an alternative parametric trend `α(μ) = a/μ + b` (`--fit-type parametric`) is available, fit by robust (MAD-trimmed) weighted least squares and subject to the same scale correction.
3. **Pairwise NB-GLM LRT + Bonferroni + BH.** Per locus, identify `k*` = condition with the highest size-factor-normalized signal and `k_bg` = condition at rank `--background-rank` (default 3, so the third-highest condition by normalized signal; 1 = `k*`, 2 = runner-up). The null (`μ_{k*} = μ_{k_bg}`) is fit by intercept-only NB MLE on this pair; the alternative is saturated on this pair with the constraint `μ_{k*} > μ_{k_bg}` (always satisfied by construction). The remaining K−2 conditions enter both models as saturated nuisance and cancel from the likelihood ratio. Under the null with `k*` fixed, the statistic `T = 2·(ℓ_alt − ℓ_null)` follows the chi-bar-squared `½·χ²(0) + ½·χ²(1)` distribution — a standard result for LRTs on a boundary constraint (Self & Liang, *JASA* 1987; Silvapulle & Sen, *Constrained Statistical Inference*, 2005). The per-locus one-sided p-value is multiplied by K (Bonferroni) to correct for selecting `k*` by argmax; the choice of `k_bg` is deterministic given the ordering and adds no extra Bonferroni cost. Benjamini–Hochberg q-values then control the FDR across loci.

The rank knob controls robustness to "competing peaks" — loci where more than one condition is active. At `--background-rank 2`, the test compares `k*` to the runner-up: any second condition that is also elevated shrinks the gap and the LRT collapses to ~0. The default of 3 compares against the third-ranked condition, so a single competing peak does not depress the test; values >3 tolerate more competing peaks at the cost of comparing `k*` against an increasingly-low background. `--background-rank` is silently capped to K when larger, so the default just works at K = 2. Restricting the test to a single pair of conditions is what keeps it *enrichment-only*: a depletion pattern (one condition far below an otherwise-uniform set) has its top conditions all at the high level, so the test does not reject — exactly the behavior we want when looking for "fertile ground" loci.

> **FDR assumes independent (or positively dependent) loci.** BH controls FDR under independence or positive regression dependency (PRDS). BigWig signal on adjacent windows is correlated — overlapping tiling windows, shared peaks, or bin sizes smaller than the underlying signal's autocorrelation length all introduce dependence. For typical "one row per peak / per gene" inputs this is fine. For dense sliding-window inputs (e.g. 100 bp windows stepped every 50 bp), q-values will be optimistic; prefer non-overlapping windows or thin by autocorrelation length before trusting the FDR.

Columns appended to the output:

| column | meaning |
| --- | --- |
| `effect_size` | log2((X_{k\*}/s_{k\*}) + pc) − log2(mean_{j≠k\*}(X_j/s_j) + pc), i.e. the log2 fold change of the enriched condition vs the mean of the other K−1 conditions on the size-factor-normalized scale (computed as a log-difference; equivalent to a log2 ratio when pc is small). Always non-negative by construction. This is a user-facing summary; the test statistic itself is computed only from the `(k*, k_bg)` pair. |
| `p_value` | Bonferroni-corrected one-sided LRT p-value: `min(K · ½ · χ²(1).sf(T), 1)` |
| `q_value` | Benjamini–Hochberg q-value |
| `enriched_condition` | name of the condition column with the highest normalized signal (k\*). Always populated (it is just `argmax(X/s)`), so it is only meaningful for rows that pass a significance threshold — on a row with `q_value ≈ 1` it is whichever column happened to be highest under noise, not a call. |
| `effect_size_pc_dominated` | `True` when some `X_j/s_j < pseudocount` for that locus, meaning the log2 effect size is dominated by `pc` rather than data. Treat the effect size as a lower bound. The LRT itself is unaffected. |
| `lrt_zero_dominated` | `True` when `X_{k*} == 0` or `X_{k_bg} == 0`, i.e. the LRT pair contains a zero. The reported p-value for these loci is driven by the internal `mu_alt` floor (1e-20) rather than by data, and they tend to dominate the top of sparse-data output as spurious hits. Treat with skepticism — common causes are regions of poor mappability or chromosome-naming mismatches in some tracks. |
| `lrt_convergence_failed` | `True` when the intercept-only NB MLE for the null fit did not converge for that locus. Its `p_value` has been set to 1.0; the column is present so users can audit how many loci hit this case. |

The LRT statistic `T` and the per-locus dispersion `α` are not written to the TSV. They are available as `lrt_stat` and `per_locus_dispersion` on the `EnrichmentResult` returned by `enrichment_analysis` (see [Python API](#python-api)).

Size factors, the number of loci that contributed to the size-factor estimate, the estimated dispersion fit, its trend coefficients, the test's effective conservativeness at the current K, and the number of kept/total loci are printed to stderr. The dispersion-fit label is one of:

| label | meaning |
| --- | --- |
| `common` | `--fit-type common` succeeded (single α = median of per-locus MoM estimates) |
| `parametric` | `--fit-type parametric` succeeded (fitted `α(μ) = a/μ + b`) |
| `common-fallback` | `--fit-type parametric` was requested but failed; fell back to the common fit |
| `override` | `--dispersion` was supplied; the fixed α was used at every locus |
| `zero` | `--fit-type zero` was supplied; Poisson was forced at every locus |

Size-factor spread, Poisson fallbacks, `--fit-type zero`, and `common-fallback` all emit a `FertilizerEnrichmentWarning` on stderr.

**Most users should keep the defaults.** `--fit-type common` with the default `--min-signal 5.0` and `--pseudocount 0.5` works well across most datasets. Touch the dispersion knobs only if (a) you have a strong prior that dispersion trends with mean signal (try `--fit-type parametric`), (b) you have an external estimate of α (pass `--dispersion`), or (c) you want a sensitivity analysis (`--dispersion 0.05` and `--dispersion 0.10`).

| CLI flag | effect |
| --- | --- |
| `-i`, `--input` | input TSV; read as gzip when the name ends in `.gz` |
| `-c`, `--conditions` | two or more column names to compare. Names that match an output column (`effect_size`, `p_value`, `q_value`, `enriched_condition`, or one of the three flag columns) are rejected |
| `-o`, `--output` | output TSV, filtered to loci passing the threshold, with extra columns appended; gzip-compressed when the name ends in `.gz` |
| `--q-threshold` | keep loci with q ≤ this (default `0.05`; set to `1.0` to keep all rows) |
| `--p-threshold` | additionally keep only loci with raw p ≤ this (default: off) |
| `--fit-type` | dispersion model: `common` (default, median of MoM estimates, bias-corrected), `parametric` (fits `α(μ) = a/μ + b`), or `zero` (forces Poisson — diagnostic only, strictly anti-conservative if real overdispersion exists) |
| `--min-signal` | minimum mean normalized signal for loci included in dispersion estimation (default `5.0`) |
| `--dispersion` | override the fitted α with a fixed value (finite, ≥ 0; 0 is Poisson) applied to every locus; bypasses `--fit-type` and `--min-signal` entirely and sets `dispersion_fit = override`. Useful for sensitivity analyses (e.g. re-run at 0.05 and 0.10 to see how much calls depend on α) |
| `--size-factors` | externally-supplied size factors, one positive value per `-c` entry in the same order. Bypasses median-of-ratios. Use when you have an external normalization you trust more (RPM/RPKM, spike-in). Pass `1 1 1 ...` to disable normalization entirely. |
| `--background-rank` | rank of the condition compared against `k*` in the LRT pair (default `3` — tolerates one competing peak; `2` compares against the runner-up; larger values tolerate more competing peaks). Capped to K when larger; the default therefore works at K = 2 without special-casing. |
| `--pseudocount` | pseudocount for the effect-size log2 transform only; does not affect the LRT. Must be > 0 (default `0.5`) |
| `--allow-non-sum` | bypass the check that the input was produced by `fertilizer extract --stat sum` or by counting BAM/CRAM/fragment input (`stat=count`). The NB-GLM assumes count-like input; `mean`/`max`/`min`/`std`/`coverage` are not counts, so p-values may be miscalibrated. Use only after empirically verifying calibration on your data. |

**Differences from DESeq2** (non-exhaustive):

- **Enrichment-only, rank-pair LRT.** DESeq2's LRT is two-sided across an arbitrary `full` vs `reduced` design and fires on both enrichment and depletion. We hard-code a 1-df one-sided LRT comparing `k*` (the argmax) against the condition at a user-chosen rank (`--background-rank`, default 3, which tolerates one competing peak), with Bonferroni × K for picking `k*` as the empirical argmax. Loci where one condition is *depleted* relative to the rest are not called by construction.
- **Per-locus dispersion MLE.** DESeq2 uses Cox-Reid adjusted profile likelihood; we use method-of-moments. MoM is less efficient per locus but is consistent and does not fail to converge.
- **Shrinkage.** DESeq2 shrinks per-locus dispersion toward the trend via a log-normal empirical-Bayes prior and retains dispersion outliers. We use the trend value directly for every locus (equivalent to infinite shrinkage, no outlier retention) — robust for the small-K setting, but cannot capture genuinely heterogeneous per-locus dispersion.
- **Bias correction.** Because median(χ²(K−1))/(K−1) < 1 for small K (0.69 at K=3, 0.84 at K=5, 0.91 at K=8), median-of-MoM underestimates α. We apply the closed-form correction.
- **Log2 fold change shrinkage.** DESeq2 optionally shrinks LFC estimates (apeglm / ashr); we report a raw log2 fold change of the enriched condition vs the mean of the others as the effect size.
- **Observation-level outliers.** DESeq2 uses Cook's distance to flag and optionally refit without outliers. We don't.
- **Independent filtering.** DESeq2 filters low-signal loci out of multiple-testing correction to maximize power at a given FDR. We don't — use `--min-signal` (dispersion-only) or pre-filter the input TSV if you want this.
- **Integer counts.** DESeq2 is designed for integer RNA-seq counts; the NB likelihood here is evaluated with `scipy.special.gammaln` and is numerically correct for any non-negative float input. (This matches the common practice of passing fractional RSEM/salmon expected counts to DESeq2 via `tximport`, and is required here because bigWig sums are real-valued.)

**Known calibration behavior.** The Bonferroni × K correction for the data-driven argmax makes this test **conservative under the null** at K ≥ 4 with the default rank, increasingly so as K grows. Empirical Type-I rates at nominal α = 0.05 on Poisson nulls under the default `--background-rank 3`: K = 2 → ≈0.05 (rank capped to 2; `k*` vs runner-up); K = 3 → ≈0.07–0.09, above nominal (rank 3 is the lowest of three conditions, where the order-statistic gap is widest; `enrich` warns about this on every K = 3 run at the default rank); K = 4 → ≈0.015; K ≥ 5 → well under 0.005. With `--background-rank 2` the test is uniformly conservative across all K: K = 3 → ≈0.008, K ≥ 5 → well under 0.001. The test suite verifies Type-I error at α = 0.05 stays under 0.10 at K = 3 and under 0.05 for K ≥ 4 under the default, and under 0.08 across the same grid at `--background-rank 2`. Power is preserved against strong effects (≥80% at 2× fold change, ≥99% at 3× fold change across K ∈ {2…8} in simulations). At low μ (< ~5) the delta-method and the median-bias correction both degrade, and low-μ loci are excluded from dispersion estimation via `--min-signal`. Depletion-only patterns (one condition low, the rest uniform) are simulated in the test suite and confirmed *not* to be called.

> **Supply a mix of positive and negative loci.** The size-factor and dispersion steps both assume that *most* loci are **not** enriched — the "null majority" is what anchors the normalization and the variance estimate. If you run `enrich` on a set of regions pre-filtered to be those you expect to change, you will get sub-optimal results: the estimated library-size differences will absorb real biological differences, and the dispersion estimate will be inflated by the true positives. For best results pass a genome-wide or background-matched set of loci containing both putative-enriched and expected-stable regions.

Example (three conditions):

```bash
fertilizer enrich -i signals.tsv -c A B C -o enrichment.tsv --q-threshold 0.05
```

Output `enrichment.tsv` (filtered, numbers illustrative; `effect_size_pc_dominated`, `lrt_zero_dominated`, `lrt_convergence_failed` columns omitted from the example for brevity):

```
chrom   start   end     A     B     C       effect_size     p_value     q_value     enriched_condition
chr1    0       1000    2.0   2.1   8.4     2.07            0.0004      0.0032      C
chr3    1200    2200    15.2  3.6   4.1     2.08            0.0011      0.0060      A
chr7    900     1900    8.8   8.9   30.1    1.78            0.0018      0.0090      C
```

## Python API

Both subcommands are thin wrappers around library functions. The same work can be done from Python:

```python
import numpy as np
import pandas as pd

from fertilizer.extract import bigwig_region_means, load_regions, FertilizerWarning
from fertilizer.enrichment import (
    enrichment_analysis,
    size_factors,
    bh_qvalues,
    EnrichmentResult,
    FertilizerEnrichmentWarning,
)

# --- extract: region means for one bigWig -------------------------
regions = load_regions(["regions.bed"])                  # chrom/start/end
values, issues = bigwig_region_means(regions, "A.bw")    # np.ndarray, set[str]
# `issues` is a subset of {"missing_chrom", "out_of_bounds", "invalid_region"}

# --- extract: counts from a BAM/CRAM or fragment file (serial) -----
from fertilizer.counting import BarcodeGroups, count_bam, count_fragments
arrays = (regions["chrom"].to_numpy(), regions["start"].to_numpy(), regions["end"].to_numpy())
bam_counts = count_bam("A.bam", *arrays, pos_shift=4, neg_shift=-5)   # (n_regions, 1)
groups = BarcodeGroups.from_table("cells.tsv", group_column="cluster")
frag_counts, seen = count_fragments("fragments.tsv.gz", *arrays, groups=groups)
# frag_counts: (n_regions, n_groups), columns in groups.names order;
# seen: chromosomes that appeared in the file

# --- enrich: enrichment NB-GLM LRT on a (n_loci, n_conditions) array
counts = pd.read_csv("signals.tsv", sep="\t", comment="#")[["A", "B", "C"]].to_numpy(float)
res: EnrichmentResult = enrichment_analysis(counts, fit_type="common")
# Per-locus arrays (length n_loci):
#   res.p_value, res.q_value, res.effect_size, res.lrt_stat,
#   res.per_locus_dispersion, res.enriched_condition_idx,
#   res.effect_size_pc_dominated, res.lrt_zero_dominated,
#   res.lrt_convergence_failed
# Per-condition array (length n_conditions):
#   res.size_factors
# Scalars:
#   res.n_loci_for_size_factors (int), res.dispersion_fit (str),
#   res.dispersion_trend (tuple[float, float]), res.background_rank (int)
```

`FertilizerWarning` (locus-level issues from `extract`) and `FertilizerEnrichmentWarning` (size-factor spread, Poisson and `common-fallback` fallbacks, `--fit-type zero`, the K = 3 calibration note, region overlap and MLE non-convergence from `enrich`) are both `UserWarning` subclasses and are unrelated to each other, so filtering one does not filter the other. From the CLI they are printed to stderr. In Python:

```python
import warnings
from fertilizer.extract import FertilizerWarning
from fertilizer.enrichment import FertilizerEnrichmentWarning, enrichment_analysis

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always", FertilizerWarning)
    warnings.simplefilter("always", FertilizerEnrichmentWarning)
    res = enrichment_analysis(counts)          # `counts` as in the example above
for w in caught:
    print(w.category.__name__, w.message)

# or silence both (simplefilter takes one category per call)
warnings.simplefilter("ignore", FertilizerWarning)
warnings.simplefilter("ignore", FertilizerEnrichmentWarning)
```

## Claude Code skill

`fertilizer` ships an agent skill for [Claude Code](https://claude.com/claude-code)
that teaches the assistant to run `extract` and `enrich` correctly: choosing
between bigWig, BAM/CRAM and fragment input (including per-cluster pseudobulks
of one fragment file), choosing a background region set, picking `--background-rank` for the
question being asked, reading the stderr diagnostics and output flags, mapping
common questions ("regions specific to X", "higher in A than B", "starting
regions for design") to commands, and diagnosing the usual failures. It is a
short router plus topic files that load only when needed.

```bash
fertilizer install-skill
```

This copies the skill into `~/.claude/skills/fertilizer`. Options:

- `-d, --directory DIR` — install into a different skills directory (default `~/.claude/skills`).
- `--symlink` — symlink the packaged skill instead of copying it, so an editable install is reflected without reinstalling.
- `-f, --force` — overwrite an existing installation. **Re-run with `--force` after upgrading `fertilizer`**; without it the command refuses to replace the old copy.

Restart Claude Code to pick it up.

## Project layout

```
fertilizer/
├── pyproject.toml             # package metadata + uv/hatchling build config
├── README.md
├── LICENSE
├── src/
│   └── fertilizer/
│       ├── __init__.py
│       ├── cli.py             # top-level argparse dispatcher
│       ├── extract.py         # signal aggregation + `extract` subcommand
│       ├── counting.py        # BAM/SAM/CRAM and fragment-file counting for `extract`
│       ├── enrichment.py      # enrichment analysis + `enrich` subcommand
│       ├── install_skill.py   # `install-skill` subcommand
│       └── skills/fertilizer/ # bundled Claude Code skill (SKILL.md + references/)
├── examples/
│   ├── make_demo_data.py     # generates a synthetic end-to-end demo dataset
│   └── README.md             # walkthrough of `extract` + `enrich` on demo data
└── tests/
    ├── test_extract.py
    ├── test_counting.py
    ├── test_enrichment.py
    └── test_install_skill.py
```

## Troubleshooting / FAQ

**`pyBigWig` won't install.** It needs system `libcurl` and `libssl` headers.
On Debian/Ubuntu: `sudo apt-get install libcurl4-openssl-dev libssl-dev zlib1g-dev`.
On macOS: `brew install curl openssl`.

**`extract` output is all zeros.** Most often a chromosome-naming mismatch
between the BED and the bigWig (`chr1` vs `1`). `extract` warns when more
than 95% of cells are exactly zero — re-check the inputs. To inspect a
bigWig's chromosome names: `python -c "import pyBigWig; print(pyBigWig.open('A.bw').chroms())"`;
a BAM's or CRAM's: `samtools view -H A.bam | grep '^@SQ'`; a fragment file's:
`zcat fragments.tsv.gz | grep -v '^#' | cut -f1 | uniq`.

**`enrich` gives `q_value` near 1 for everything.** Four common causes,
in order of likelihood:
1. K is large (≥ 5) and the Bonferroni × K argmax correction makes the test
   very conservative. The stderr output prints the effective Type-I rate
   at α=0.05 for your K and rank; if it is far below nominal, low power
   is by construction. Lower K (combine biologically equivalent conditions)
   or accept that strong effects only will be called.
2. Your enriched loci have a **competing peak** (a second condition also
   elevated). The default `--background-rank 3` tolerates one such peak,
   but if you have two or more, raise it (e.g. `--background-rank 4`).
   `--background-rank 2` is the most fragile choice in this regard.
3. Pre-filtering. If you passed a region set already enriched for the
   conditions you care about, the null-majority assumption is violated;
   the size factors absorb real differences and dispersion is inflated.
   Pass a genome-wide / background-matched region set instead.
4. The dispersion estimate fell back to Poisson (`alpha=0`) because too
   few loci passed `--min-signal`. Stderr will say so. Lower `--min-signal`
   or supply more loci.

**I have replicates per condition.** This package targets the
one-replicate-per-condition setting. With replicates, prefer DESeq2 or
edgeR — they estimate per-locus dispersion from within-condition variance
and give better power. `fertilizer` is the tool to reach for when you
have one bigWig per condition (a common setup for ChIP-seq, ATAC-seq, and
many predictive models) and want enrichment calls without making up
fake replicates.

**Which conditions should I use as "background"?** None — the test is
symmetric and compares each condition against the others. Just supply
all of your conditions as `-c`.

**Can I pre-normalize my data?** Yes. Either (a) pass `--size-factors V1 V2 ...`
matched to your `-c` entries to inject an external normalization (RPM/RPKM,
spike-in, etc.), or (b) pass `--size-factors 1 1 1 ...` to disable
normalization entirely if your inputs are already on the same scale.

**`effect_size_pc_dominated == True` for some loci.** Those loci have at
least one condition with normalized signal below `--pseudocount`, so the
log2 ratio is dominated by the pseudocount rather than data — treat the
effect size as a lower bound only. The p-value/q-value are unaffected
(they don't use the pseudocount).

## Citation

If you use `fertilizer` in published work, please cite:

> Schreiber, J. *fertilizer: per-region enrichment from single-replicate
> bigWig signal across conditions.* https://github.com/jmschrei/fertilizer,
> v0.2.0 (2026).

## License

MIT — see [LICENSE](LICENSE).

## Changelog

See [CHANGELOG.md](CHANGELOG.md).
