Metadata-Version: 2.5
Name: valius-kopya
Version: 2.1.0
Summary: Expression-only single-cell copy-number-variation caller.
Project-URL: Homepage, https://github.com/ValiusSciences/kopya
Project-URL: Repository, https://github.com/ValiusSciences/kopya
Project-URL: Issues, https://github.com/ValiusSciences/kopya/issues
Author-email: Bulent Arman Aksoy <arman@valiussciences.com>
License-Expression: MIT
License-File: LICENSE
Keywords: CNV,bioinformatics,cancer,copy-number,genomics,scRNA-seq,single-cell
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.11
Requires-Dist: anndata>=0.10
Requires-Dist: click>=8.1
Requires-Dist: igraph>=0.11
Requires-Dist: leidenalg>=0.10
Requires-Dist: numpy<2,>=1.26
Requires-Dist: pandas>=2.2
Requires-Dist: pyarrow>=15
Requires-Dist: pyucell>=0.6
Requires-Dist: ruptures>=1.1.9
Requires-Dist: scanpy>=1.10
Requires-Dist: scikit-learn>=1.4
Requires-Dist: scipy>=1.11
Requires-Dist: tqdm>=4.66
Provides-Extra: dev
Requires-Dist: matplotlib>=3.7; extra == 'dev'
Requires-Dist: pytest-cov>=5; extra == 'dev'
Requires-Dist: pytest>=8; extra == 'dev'
Provides-Extra: plot
Requires-Dist: matplotlib>=3.7; extra == 'plot'
Description-Content-Type: text/markdown

# kopya

[![PyPI](https://img.shields.io/pypi/v/valius-kopya)](https://pypi.org/project/valius-kopya/)
[![Python](https://img.shields.io/pypi/pyversions/valius-kopya)](https://pypi.org/project/valius-kopya/)
[![License](https://img.shields.io/pypi/l/valius-kopya)](https://github.com/ValiusSciences/kopya/blob/main/LICENSE)
[![tests](https://github.com/ValiusSciences/kopya/actions/workflows/test.yml/badge.svg)](https://github.com/ValiusSciences/kopya/actions/workflows/test.yml)

Expression-only single-cell CNV caller: a pure-Python, pip-installable tool that finds the diploid baseline, segments the genome, and calls tumor vs normal automatically. An end-to-end, R-free alternative to CopyKAT, SCEVAN, and inferCNV.

## Motivation

Copy-number calling from scRNA-seq is already well served by several mature tools: CopyKAT, SCEVAN, and inferCNV in R, and infercnvpy in Python. kopya fills a specific gap. It is a single pip-installable Python package and an **end-to-end automated caller**, not just a CNV-signal generator: it finds the diploid baseline, segments the genome into discrete CNV regions, and emits per-cell tumor / normal / uncertain calls with QC gates and CopyKAT-compatible outputs, sparse-first and fast enough for large cohorts.

For the detailed, tool-by-tool comparison, algorithm, scope, runtime, and exactly where kopya differs from each (including the closest peer, infercnvpy), see [`docs/comparison.md`](https://github.com/ValiusSciences/kopya/blob/main/docs/comparison.md). Allele-aware calling (LOH, copy-neutral) stays the domain of Numbat.

## What we borrowed from the field

The algorithm is not novel; it is a careful synthesis of the best ideas from the published tools:

| Idea | Taken from |
|------|-----------|
| CP10k + log1p normalization, then project onto GENCODE coordinates | Standard scRNA-seq practice; CopyKAT's preprocessing |
| UCell scoring against non-malignant signatures to find the diploid baseline | SCEVAN |
| Lowest-variance Leiden cluster as a diploid fallback | CopyKAT |
| Pool tumor signal across cells before segmentation, then fill per-cell CN inside the shared segment table | SCEVAN's multichannel trick |
| PELT changepoint detection (L2 cost) with a BIC-style penalty | ruptures library; used by SCEVAN |
| 2-component GMM on L1 deviation from baseline, giving tumor / normal | CopyKAT |
| Per-segment CN matrix, Leiden-clustered into subclones | SCEVAN |

Our contribution is twofold. First, **an end-to-end automated caller**: kopya finds the diploid baseline, segments the genome into discrete CNV regions, and calls tumor / normal / uncertain with QC gates, rather than handing you a signal matrix to interpret yourself. Second, **the engineering**: a single sparse-first Python package that runs the same ideas roughly ~100× faster than CopyKAT's reported 1-2 h runtime and in a fraction of the memory (kopya does 10k cells in ~12 s at ~1.8 GB; see the scaling table below), installs in one step, and writes CopyKAT-style output files (`prediction.csv`, `chr_cnv_matrix.csv`) that drop into a CopyKAT-based workflow with minimal changes.

## How it works

Five-step pipeline:

| Step | Module | What it does |
|------|--------|-------------|
| 1. Normalize | `normalize.py` | CP10k + log1p; project genes onto GENCODE v49 hg38 coordinates; drop chrX / chrY / chrM, `MT-*`, `HLA-*`, immunoglobulin (IGH/IGK/IGL) and cell-cycle genes. Per-chromosome outputs are autosome-only (chr1-22) |
| 2. Pick baseline | `baseline.py` | 4-mode cascade: **supervised** (barcodes you supply), then **signature** (UCell on bundled immune/stromal/endothelial gene sets, à la SCEVAN), then **variance cluster** (Leiden + lowest-variance cluster, à la CopyKAT), then **GMM fallback** |
| 3. Smooth + segment | `smooth.py`, `segment.py` | Center each cell against the per-gene normal-pool median; per-chromosome moving average; PELT changepoint detection on the **pooled tumor signal**; per-cell CN computed inside the shared segment table |
| 4. Classify + subclones | `classify.py` | 2-component GMM on per-cell L1 distance from baseline, giving tumor / normal / uncertain; Leiden on tumor cells' per-segment CN matrix, giving subclones (max 5) |
| 5. Write outputs | `outputs.py` | `prediction.csv`, `chr_cnv_matrix.csv` (CopyKAT-style columns), IGV `.seg`, `qc.json` |

Performance scales near-linearly with cell count. The table below is 10 replicate runs per size of the full `kopya run` pipeline end to end (load through outputs) on an Apple M5 Pro (15 cores, 24 GB RAM, macOS), peak resident memory via `/usr/bin/time`; values are mean ± sd. Absolute times are hardware-dependent; the near-linear scaling is the point. Reproduce with [`scripts/benchmark_scaling.py`](https://github.com/ValiusSciences/kopya/blob/main/scripts/benchmark_scaling.py):

| Cells | Wall time (s) | Peak RAM (GB) |
|------:|--------------:|--------------:|
| 5,000 | 10.9 ± 0.3 | 1.17 ± 0.01 |
| 10,000 | 12.2 ± 0.3 | 1.82 ± 0.01 |
| 25,000 | 15.8 ± 0.4 | 3.54 ± 0.11 |
| 50,000 | 22.0 ± 0.7 | 4.22 ± 0.11 |
| 100,000 | 39.7 ± 0.6 | 4.55 ± 0.07 |
| 200,000 | 70.1 ± 10.2 | 7.28 ± 0.43 |
| 400,000 | 119.3 ± 3.8 | 11.52 ± 0.53 |
| 800,000 | 242.1 ± 4.4 | 14.91 ± 0.57 |

At 10k cells that is ~12 s and ~1.8 GB peak RSS. CopyKAT's own reports put it at 1-2 h and 30-60 GB on comparable data, so kopya is roughly two orders of magnitude faster in a much smaller footprint. Those R-tool figures are the tools' published/reported numbers, not measured here; the only tool benchmarked head-to-head is infercnvpy (see [`docs/comparison.md`](https://github.com/ValiusSciences/kopya/blob/main/docs/comparison.md)). kopya also scales to ~800k cells in about 4 minutes at ~15 GB, beyond where the R tools stay practical on commodity hardware.

## Install

```bash
# From PyPI. The distribution is `valius-kopya`; the import name and CLI are `kopya`.
pip install valius-kopya

# With the optional plotting stack (matplotlib, for `plot-heatmap` and the Python-API heatmap):
pip install "valius-kopya[plot]"
```

Development install from a clone:

```bash
# Create the conda env, run from the repo root (Python + numpy/scipy/scanpy/sklearn/ruptures/leidenalg/pyucell).
# The yaml's `-e .` editable-installs this package, and it resolves against the
# directory you invoke mamba from, so the working directory must be the repo root.
mamba env create -f envs/kopya.yaml
mamba activate kopya

# The env already editable-installs this package. Re-run to add the optional
# plotting stack (matplotlib, for `plot-heatmap` and the Python-API heatmap)
# and the test tools (pytest):
pip install -e ".[plot,dev]"

# Confirm the install (full unit-test suite, ~2 min).
pytest
```

The package ships a bundled GENCODE v49 hg38 gene-order table (~77k genes, 2.6 MB) and a Tirosh 2016 cell-cycle gene list; no external resources required to run.

## Quickstart

A runnable first pass on the committed 500-cell demo fixture, no data to download:

```bash
# Run the pipeline, then render the inferCNV-style heatmap.
kopya run --anndata tests/fixtures/tiny_simulated.h5ad --out-dir /tmp/kopya-demo --sample demo
kopya plot-heatmap --run-dir /tmp/kopya-demo --out /tmp/kopya-demo/heatmap.png --scale linear
```

This writes `prediction.csv`, `chr_cnv_matrix.csv`, `segments.parquet`, `qc.json`, an IGV `.seg`, and the heatmap into `/tmp/kopya-demo/`. The fixture has three planted events (chr7 gain, chr10 loss, focal MYC amp), so it's a good place to build intuition. See [`examples/tiny_simulated/`](https://github.com/ValiusSciences/kopya/tree/main/examples/tiny_simulated) for a pre-computed copy of these outputs and [`notebooks/working_with_results.ipynb`](https://github.com/ValiusSciences/kopya/blob/main/notebooks/working_with_results.ipynb) for a walkthrough of every file.

## Inputs

Four accepted entry forms, all producing identical downstream results, listed most-common first.

### CellRanger directory (Form C)

The standard external format: the `filtered_feature_bc_matrix/` directory that CellRanger, STARsolo, alevin-fry, and most droplet-sequencing pipelines write:

```
outs/
└── filtered_feature_bc_matrix/
    ├── barcodes.tsv.gz
    ├── features.tsv.gz      # genes.tsv.gz for CellRanger v2
    └── matrix.mtx.gz
```

```bash
kopya run --cellranger-dir outs/filtered_feature_bc_matrix/ --out-dir results/ --sample my-sample
```

### CellRanger H5 (Form D)

The single-file HDF5 alternative from the same CellRanger run:

```bash
kopya run --cellranger-h5 outs/filtered_feature_bc_matrix.h5 --out-dir results/ --sample my-sample
```

> **Loupe files (`.cloupe`) are not supported.** They are a proprietary format with no public parser. Use the `.h5` or `filtered_feature_bc_matrix/` output from the same CellRanger run.

### AnnData h5ad (Form A)

`adata.X` (or `layers["counts"]`) holds raw integer UMI counts. `var_names` are HGNC gene symbols. `obs_names` are cell barcodes.

```bash
kopya run --anndata processed_ad.h5ad --out-dir results/ --sample my-sample
```

### Matrix trio (Form B)

A Matrix Market counts file alongside cell- and gene-metadata CSVs, in this layout:

```
results/
├── normalization/
│   └── counts.mtx
└── metadata/
    ├── obs.csv
    └── var.csv
```

### Finding inputs automatically

`find-inputs` searches a directory tree for any of the four forms and prints the ready-to-run command. Detects the CellRanger `outs/` layout, the matrix-trio layout, and any `.h5ad`. Warns clearly if a `.cloupe` file is found.

```bash
kopya find-inputs --patient-dir /path/to/patient --sample sample-01
```

### Optional inputs

- `--norm-cell-names path.txt`: one barcode per line; bypasses the baseline cascade and uses these cells as the diploid pool. Useful when you have manual immune-cell annotations. **For some tumors this is not optional; see [When to use supervised mode](#when-to-use-supervised-mode).**
- `--signatures path.json`: override the bundled non-malignant signature library (same schema as the built-in `normal_signatures.json`: an optional `_meta.non_malignant_labels` allow-list plus one `"<label>": [genes...]` entry per signature). Retunes the **signature** baseline mode for tissues the defaults don't fit.
- `--non-malignant-labels A,B,C`: override just the allow-list of cell types treated as "normal", without authoring a whole JSON. Add back a label the defaults exclude (e.g. `Fibroblast` on a sample whose normal stroma is genuinely fibroblast-rich) or drop one your tumor happens to score on. Validated against the active signature set, so a typo fails fast.
- `--gene-order path.tsv`: overrides the bundled GENCODE table (4 columns: `gene_symbol chr start end`).

> **The bundled signatures and allow-list are defaults, not universal truth.** They're curated for immune-infiltrated solid tumors. On mesenchymal (sarcoma/osteosarcoma) or hematologic samples the defaults can mislabel the malignant population as "normal": `Plasma_cell` is intentionally **excluded** from the allow-list, for instance, because plasma cells are the malignancy in myeloma, and `Fibroblast` is excluded because mesenchymal / ECM-like tumor cells score on the fibroblast program and would otherwise seed their own baseline. Both remain defined as signatures, so such cells still resolve to them and are kept out of the seed rather than falling through to a label that is allow-listed. If a run looks off (see [When to use supervised mode](#when-to-use-supervised-mode)), retune with `--non-malignant-labels` / `--signatures`, or pin the baseline outright with `--norm-cell-names`. The effective allow-list is recorded in `qc.json` (`baseline_diagnostics.non_malignant_labels`).

## Running

### Quick run, all defaults

```bash
kopya run \
    --counts  /path/to/counts.mtx \
    --obs     /path/to/obs.csv \
    --var     /path/to/var.csv \
    --out-dir results/ \
    --sample  sample-01
```

### All the knobs

```bash
kopya run \
    --anndata /path/to/processed_ad.h5ad \
    --out-dir results/ \
    --sample  sample-01 \
    --threads 8 \
    \
    # cell + gene filtering
    --min-genes 200       \
    --low-dr    0.05      \
    --up-dr     1.0       \
    \
    # baseline picking
    --normal-score-threshold 0.3   \
    --min-normal-cells       50    \
    --norm-cell-names        immune_barcodes.txt   `# optional supervised override` \
    --non-malignant-labels   T_cell,B_cell,Endothelial   `# optional allow-list override` \
    --signatures             custom_signatures.json       `# optional signature library` \
    \
    # smoothing + segmentation
    --smooth-window  100   \
    --penalty-coef   1.0   \
    --min-seg-genes  25    \
    \
    # classification
    --call-confidence      0.5  \
    --max-subclones        5    \
    --coherence-gate       0.45 `# downgrade scattered-signal tumor calls to uncertain; 0 disables, lower for focal-amp tumors` \
    --low-complexity-frac  0.5  `# flag ambient cells below this fraction of the median detected-gene count; 0 disables` \
    --subclones            `# or --no-subclones`
```

Every knob has a sensible default; the bare-minimum invocation above is all most runs need. See `kopya run --help` for the full surface.

## Python API (scanpy-style)

Prefer to work in memory alongside scanpy / [infercnvpy](https://infercnvpy.readthedocs.io/)? Run on an AnnData and the results are written back onto it, no files, no CLI:

```python
import scanpy as sc
import kopya as kp

adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")   # raw counts in .X
# ... your usual cell-type annotation pass populates adata.obs["cell_type"] ...

kp.tl.cnv(adata, reference_key="cell_type",
          reference_cat=["T cell", "B cell", "Macrophage", "Endothelial"])

adata.obs["cnv_class"]          # 'tumor' / 'normal' / 'uncertain'  (NaN = QC-filtered)
adata.obs["cnv_score"]          # per-cell aneuploidy score
adata.obs["cnv_subclone"]       # Leiden subclone among tumor cells ('' for non-tumor)
adata.obs["cnv_low_complexity"] # bool: ambient / low-gene-count cell flag
adata.obsm["cnv_chr"]           # cells × chromosomes, 1.0-centered (>1 gain, <1 loss)
adata.uns["cnv"]                # chroms, baseline method, params, segment table

kp.pl.chromosome_heatmap(adata, groupby="cell_type")   # needs the [plot] extra (matplotlib)
```

Pass `reference_key`/`reference_cat` (your normal cell types) to run supervised; recommended, and required for mesenchymal tumors (see below). Omit them for the unsupervised cascade. Raw counts are required (the tool does its own CP10k+log1p); if they live in a layer, pass `layer="counts"`. The heatmap needs the `[plot]` extra (`pip install -e ".[plot]"`).

**Coming from Seurat?** See the round-trip walkthrough (annotate once in Seurat, run kopya, bring the calls back) in the executed notebook [`docs/tutorials/seurat.ipynb`](https://github.com/ValiusSciences/kopya/blob/main/docs/tutorials/seurat.ipynb) (previews inline on GitHub, with a rendered CNV heatmap).

For ~1M-cell cohorts, the CLI's `run` command has a memory-lean streaming path; `tl.cnv` is the in-memory convenience path for typical objects.

## When to use supervised mode

The unsupervised baseline cascade (signature, then variance, then GMM) works out of the box for most samples, but it has a failure mode worth knowing about: **on some tumors it inverts.** It picks the *tumor* itself as the diploid baseline, so every call comes out backwards: tumor cells look normal, normal cells look aneuploid, and the per-chromosome profile anti-correlates with matched bulk.

This happens when the tumor's expression profile matches one of the bundled "normal" signatures:

- **Mesenchymal / non-immune tumors** (osteosarcoma, other sarcomas). The malignant cells express fibroblast / osteoblast markers, so the signature step labels the tumor population "normal."
- **Full-length protocols (SMART-seq2).** UCell rank-scoring misfires on read-length counts; see the Patel 2014 XFAIL under [Validation](#validation).
- More generally, any sample where the genuinely-diploid cells are too small a minority for the cascade to lock onto.

**How to spot an inverted baseline:**

- `qc.json` reports a `baseline_method` of `signature`/`variance` with an implausibly large `n_normal_seed` (most of the sample called "normal") and a compressed `tumor_signal_max_abs`.
- Known-normal/immune cells (T cells, plasma cells) score *highest*, not lowest: the tell-tale sign the reference frame is flipped.
- If you have matched bulk, the single-cell profile anti-correlates with it.

**The fix: seed the baseline yourself.**

```bash
kopya run --anndata processed_ad.h5ad --out-dir results/ --sample my-sample \
    --norm-cell-names normal_barcodes.txt
```

`--norm-cell-names` (one barcode per line) bypasses the cascade and uses exactly those cells as the diploid pool (`baseline_method: supervised`, the top tier of the cascade). Source the barcodes from any upstream cell-typing; a T-cell / immune cluster from your Leiden annotation is the usual choice, since lymphocytes are reliably diploid.

> **Worked example: a mesenchymal tumor (e.g. osteosarcoma).** Unsupervised, the osteoblast/fibroblast signatures capture the mesenchymal tumor itself as the baseline, so the single-cell profile anti-correlates with matched bulk (negative concordance) and plasma cells come out as the highest-scoring "tumor" population. Re-seeding the baseline with an annotated T-cell barcode list flips the sign back to positive concordance (strong-event recall goes from 0 to complete) and clears the plasma-cell false positives. The lesson: for mesenchymal and other inversion-prone tumors, **supervised mode is not a fallback, it's the right call.**

**Lighter-weight alternative.** If you don't have curated barcodes but you *do* know the cascade is picking the wrong cell type as normal, you can keep the unsupervised cascade and just retune which signatures count as "normal", e.g. `--non-malignant-labels T_cell,B_cell,Endothelial,NK_cell` to exclude the mesenchymal labels, or supply a whole custom `--signatures` library. See [Optional inputs](#optional-inputs).

## High ambient RNA

Ambient RNA — free-floating transcripts counted into every droplet — is the same background added to every cell, and kopya reads copy number as a *ratio* against the diploid pool. Adding it to both sides pushes that ratio toward 1: a real deviation is attenuated by roughly `(1 - rho)`, where `rho` is the ambient fraction, estimated here with [SoupX](https://github.com/constantAmateur/SoupX)'s `autoEstCont`. **The CNV is still there, just fainter**, and normalization will not recover it. **In practice this is only worth acting on above `rho` ≈ 0.5** — routine contamination attenuates the signal too, but correcting it changed nothing measurable on our cohort, so treat the rest of this section as a checklist for a suspected-dirty sample rather than a routine step.

**How to spot it** from `qc.json` alone, with no ambient estimate:

- A compressed `tumor_signal_max_abs` alongside an implausibly large `n_normal_seed`.
- `baseline_diagnostics.per_label_counts` dominated by a single label.
- A high `baseline_diagnostics.median_top_score` — cells match one signature strongly because they all look alike.

These overlap with the inverted-baseline tells in [When to use supervised mode](#when-to-use-supervised-mode); heavy ambient makes an inversion more likely, so check both. (Both `baseline_diagnostics` fields need `baseline_method: signature`.)

**What to do:**

- **`rho` up to ~0.2: don't pre-correct.** Check instead that no mesenchymal label is seeding the baseline; `Fibroblast` is no longer allow-listed by default (see [Optional inputs](#optional-inputs)).
- **`rho` above ~0.5: correct**, and don't expect supervised mode alone to rescue it — the ambient sits in the reference cells and the tumor alike, so subtracting even a perfect reference subtracts mostly ambient.
- **If you correct, check what was delivered rather than the reported `rho`.** Rounding corrected counts back to integers erases sub-0.5 per-entry subtractions; across our sub-0.2 samples the delivered fraction was a median ~60% of the estimate, and as little as a third.

**Not the same as the `low_complexity` flag.** Both get called "ambient" here. `low_complexity` asks whether *this barcode* caught a real cell: it fires below `--low-complexity-frac` × the sample's median detected-gene count (0.5 by default) and bars that cell from seeding the baseline (see [Interpreting `prediction.csv`](#interpreting-predictioncsv)). `rho` asks how much background is in the cells you *keep* — a normal-looking cell sitting in `rho` 0.4 of soup is never flagged. And because the threshold tracks the sample's own median, heavy ambient drags it down: **a low `n_low_complexity` is not evidence the sample is clean.**

**What we measured.** On a 12-sample internal cohort (SoupX `rho` 0.033 to 0.737), ambient-corrected counts changed median AUC by +0.003 and improved matched-bulk concordance in only 3 of 11 scorable samples (median -0.007). Only the one sample above `rho` 0.5 gained materially (+0.10 AUC), and it carries the cohort's weakest ground truth — treat ~0.5 as indicative, not calibrated. This cohort is internal and is not reproducible from this repository.

## Outputs

All artifacts land under `--out-dir`:

| File | Shape / Format | Purpose |
|------|----------------|---------|
| `prediction.csv` | cells × `{class, confidence, tumor_score, subclone, n_segments_altered, low_complexity}` | Per-cell call, in CopyKAT-style column semantics (a `prediction`/`class` column plus per-cell scores); `low_complexity` is appended last. Aligns with CopyKAT's per-cell prediction table so downstream consumers need little adaptation |
| `chr_cnv_matrix.csv` | cells × chromosomes | Per-chromosome CN summary, centered on 1.0 (`>1` gain, `<1` loss), in the same cells × chromosomes layout CopyKAT emits. Calibrated **up to a per-cell scale factor**: each cell's whole row carries an offset that cancels within a cell, leaves a pseudobulk factor-weighted rather than equally weighted, and does not cancel at all when cells are compared to each other at one chromosome — centre first for per-cell work (see [Per-cell offset](#the-per-cell-offset-in-chr_cnv_matrixcsv)) |
| `chr_cnv_matrix_centered.csv` | cells × chromosomes | **Optional**, written only with `--centered-chr-matrix`. The same matrix with each cell's row divided by its own median chromosome, so cells are comparable at a fixed chromosome. `1.0` = that cell's median chromosome, **not** diploid — for visualization and relative per-cell reads, not absolute ploidy. Additive: the raw `chr_cnv_matrix.csv` is always written unchanged |
| `{sample}_clones.seg` | IGV `.seg` | Per-clone consensus segments, loadable directly into IGV. Clone IDs prefixed with the sample name so multi-sample sessions stay disambiguated |
| `segments.parquet` | per-segment summary (one row per segment) | The **segment-to-genome map**. Original columns `chr`, `start_idx`/`end_idx` (gene-axis indices, start inclusive/end exclusive), `n_genes`, `tumor_mean` (pooled-tumor log-deviation: `>0` gain, `<0` loss), plus three appended: `segment_id` (0-based key), `start_bp`/`end_bp` (genomic span of the segment, the min gene start and max gene end across its genes). `segment_id`/`start_bp`/`end_bp` make this self-contained: you can locate any segment without the source `adata.var`. Segment row `i` equals column `i` of `cn_per_segment.npz` |
| `cn_per_segment.npz` | dense (cells × segments) | Per-cell × per-segment CN matrix (`cn`) plus `cell_barcodes`; column `i` corresponds to `segment_id == i` in `segments.parquet`. Stored **raw** (log-space, carries a per-gene pedestal + per-cell offset); re-center before plotting (see [How the signal is recovered](#visualize-the-calls-infercnv-style-heatmap)). Optimized for numpy consumers |
| `cn_per_segment_denoised.npz` | dense (cells × segments) | **Optional**, written only with `--denoise-outputs`. The **recentered + inferCNV-denoised** signal (the exact transform `plot-heatmap` draws), for a clean, ready-to-read gain/loss matrix: `0` ≈ diploid, `>0` gain, `<0` loss. Strong hallmark events survive; the normal noise floor collapses to ~0. Additive: the raw `cn_per_segment.npz` is always written unchanged. Also stores `sd_amplifier` |
| `qc.json` | dict | Run-level diagnostics: cell/gene counts, baseline method, subclone counts, per-step timings |

> **New to the output files?** Start with the worked example and notebook:
> [`examples/tiny_simulated/`](https://github.com/ValiusSciences/kopya/tree/main/examples/tiny_simulated) holds a full set of
> results from a tiny 500-cell run, and
> [`notebooks/working_with_results.ipynb`](https://github.com/ValiusSciences/kopya/blob/main/notebooks/working_with_results.ipynb)
> walks through every file: what each column means, the numeric scale of each
> value (log-space vs 1.0-centered, when negatives are expected), and how to
> re-center `cn_per_segment.npz` into gain/loss.

#### Value scales at a glance

Different files center their values differently; this trips people up, so keep it straight:

| File | Scale | Diploid value | Negatives? |
|------|-------|---------------|-----------|
| `chr_cnv_matrix.csv` | 1.0-centered multiplicative ratio | `1.0`, up to a per-cell factor | no (always > 0) |
| `chr_cnv_matrix_centered.csv` (opt-in) | 1.0-centered multiplicative ratio, per-cell centered | `1.0` = the cell's median chromosome | no (always > 0) |
| `segments.parquet` `tumor_mean` | natural-log deviation | `0` | **yes** |
| `cn_per_segment.npz` | natural-log signal (raw; re-center first) | `~0` after re-centering | **yes** |
| `{sample}_clones.seg` `seg.mean` | natural-log ratio vs normal reference | `0` | **yes** |
| `cn_per_segment_denoised.npz` (opt-in) | natural-log deviation, denoised | `0` | **yes** |
| `prediction.csv` `tumor_score` | signed consensus-template projection, in log-deviation amplitude units (within-run rank, not portable) | `~0` | **yes** |
| `prediction.csv` `cn_burden` | gene-weighted mean \|deviation\| across the genome (within-run rank, not portable) | `~0` | no |

#### Denoised outputs: scope and rationale

`--denoise-outputs` writes **one** extra file, `cn_per_segment_denoised.npz`, and nothing else. Denoise is a *per-cell, per-segment* soft-threshold against the reference noise band, so it only helps at the full resolution where that noise actually lives. It is applied at plot/output time only; the raw `cn_per_segment.npz` (and every other file) is always written unchanged. Use the denoised matrix for a clean "is the hallmark event there?" read (strong events survive, the normal floor collapses to ~0); keep the raw matrix as the source of truth for magnitude, since denoise is lossy and can attenuate weak/focal events.

The other outputs deliberately have **no** `_denoised` twin; they don't need one:

- **`chr_cnv_matrix.csv`** is already aggregated to chromosome level, which averages per-segment noise toward 0 (it's why raw chr7 already reads ~1.8 in tumor vs ~1.0 in normal). Chromosome aggregation is itself a form of denoising, and this file uses a scalar pedestal to preserve bulk concordance. Separately from denoise, it does keep the **per-cell offset** (a scalar pedestal is pool-wide, not per-cell); `--centered-chr-matrix` removes that, see [Per-cell offset](#the-per-cell-offset-in-chr_cnv_matrixcsv).
- **`{sample}_clones.seg`** is a per-subclone *median* across cells (with the per-segment baseline already subtracted); a cross-cell median is a strong denoiser, so an explicit denoise would be roughly a no-op.
- **`segments.parquet` `tumor_mean`** is a pooled summary vector, not a per-cell matrix; there is nothing per-cell to threshold.
- **`prediction.csv`** holds calls/scores, not a signal matrix. (`tumor_score` is computed on the *raw* matrix, by design.)

Want a denoised **chromosome-level** view? Aggregate the denoised per-segment matrix: denoise first, then length-weight-average per chromosome. Aggregating the raw matrix and denoising afterward is not equivalent; the order matters.

#### The per-cell offset in `chr_cnv_matrix.csv`

`chr_cnv_matrix.csv` is calibrated only **up to a per-cell scale factor**. The per-gene centering in step 3 leaves a positive pedestal on every cell — genes detected in ≤50% of the normal pool have a normal-median of 0 and pass through uncentered, so a cell's pedestal grows with how many genes it detected — and the scalar baseline subtracted at write time removes only the pool-wide part of it.

What is left is a factor on a cell's **whole row**. It cancels exactly in one place: **within one cell** — a ratio between two chromosomes, or a ranking of a cell's own chromosomes.

It does **not** cancel when cells are compared to each other at a fixed chromosome. On real data that offset correlates strongly with sequencing depth (Pearson `r ≈ 0.91` against `log10(n_counts)` on a 191k-cell glioblastoma sample) and can be larger than the per-chromosome biology, so a heatmap, ranking or clustering built on the raw values may be showing depth rather than copy number. Depth is a strong correlate and the detected-gene mechanism above explains it, but this is not a controlled demonstration that depth is the sole cause.

**An average over cells does not cancel it either.** Write a raw value as `s_c × t_ck` (the cell's factor × the true ratio); a pseudobulk is then

```
mean_c(s_c · t_ck) = mean(s) · mean(t_k) + cov_c(s, t_k)
```

— a factor-weighted average, not the plain mean, and the two terms fail differently. `mean(s)` is one constant shared by every chromosome, so it rescales the profile without changing its shape: the per-chromosome Pearson this file is validated on is scale-invariant and therefore blind to it, which is why the matched-bulk concordance never flagged the offset — though it does mean a pseudobulk's absolute level is not diploid-calibrated. The `cov` term is the one that can bite: it bends the profile chromosome by chromosome whenever the per-cell factor correlates with a chromosome's value across the pooled cells — high-depth cells over-represented in the clone that carries an event, say. That has not shown up against matched bulk, but it is an empirical observation about these samples, not a mathematical guarantee.

So: the raw file is right for **within-cell** work, and it remains the file validated against matched bulk — just read those pseudobulks as factor-weighted, not equally weighted. Centre first for a pseudobulk that should weight every cell equally, and for anything that ranks, sorts, colours or clusters **individual cells**. Either pass `--centered-chr-matrix` to get `chr_cnv_matrix_centered.csv`, or do it in one line:

```python
m = pd.read_csv("chr_cnv_matrix.csv", index_col="barcode")
centered = m.div(m.median(axis=1), axis=0)   # what --centered-chr-matrix writes
lg = np.log2(centered)                       # optional: symmetric scale for plotting
```

The trade-off is why this is opt-in rather than the default: centering makes `1.0` mean *this cell's median chromosome* instead of *diploid*. For a cell whose median chromosome is genuinely altered — a whole-genome doubling, or a tumor cell with most chromosomes gained — that misstates ploidy. The raw file keeps the absolute frame and stays the source of truth for copy-number level; the centered file is for visualization and relative per-cell interpretation. Note the classifier, `plot-heatmap` and `cn_per_segment_denoised.npz` all remove this offset internally, so calls in `prediction.csv` are unaffected by it.

### Example `qc.json`

```json
{
  "sample": "sample-01",
  "n_cells_loaded": 77418,
  "n_genes_loaded": 18487,
  "n_cells_post_m1": 77418,
  "n_genes_post_m1": 11379,
  "baseline_method": "signature",
  "n_normal_seed": 40063,
  "n_segments": 182,
  "n_tumor": 18096,
  "n_normal_called": 56981,
  "n_uncertain": 2441,
  "n_low_complexity": 9663,
  "n_outlier_fenced": 0,
  "n_anti_aligned": 0,
  "anti_alignment_tumor_n": 18096,
  "n_subclones_observed": 5,
  "subclone_counts": {"subclone_1": 10097, "subclone_2": 3045, "subclone_3": 2606, "subclone_4": 2049, "subclone_5": 299},
  "timings_secs": {
    "load": 0.6,
    "m1_filter_normalize_project": 5.0,
    "m2_pick_baseline": 6.2,
    "m3_smooth_segment_cn": 19.1,
    "m4_classify_write_outputs": 10.4
  }
}
```

`baseline_method` is the most important diagnostic; it tells downstream how confident the baseline call is (`supervised` > `signature` > `variance` > `gmm_fallback`).

### Interpreting `prediction.csv`

The intended per-cell call is the **`class`** column (`tumor` / `normal` / `uncertain`); it is the GMM-derived label, not a threshold you apply yourself to a raw score. The other columns are supporting signals for triage:

| Column | Meaning | How to read it |
|--------|---------|----------------|
| `class` | Tumor / normal / uncertain call from the 2-component GMM on `tumor_score` | The call. `uncertain` = GMM posterior below `--call-confidence` (default 0.5); filter these out for clean downstream analysis |
| `confidence` | Max GMM posterior in `[0, 1]` | How sure the call is. Cells near 0.5 are genuinely ambiguous; sort tumor calls by this to find borderline ones |
| `tumor_score` | **Signed** alignment between the cell's CN deviation and this sample's own consensus CN profile, gene-count weighted (log-space) | Higher = deviating *with* the clone. `~0` = no CN. **Negative = deviating against it**, which is a transcriptional outlier, not a "very normal" cell. See caveat below |
| `cn_burden` | Non-negative magnitude companion: gene-count-weighted mean \|per-segment deviation\| across the genome | "How much CN does this cell carry", direction ignored. Rank on this — not on `tumor_score` — when you want the cells *closest to diploid* |
| `n_segments_altered` | Count of segments deviating >0.2 log-units from baseline | Often more interpretable than `tumor_score`: "how much of the genome is rearranged" rather than a magnitude. A real tumor cell lights up many segments; a high score from one or two noisy focal events will have a low count here |
| `subclone` | Leiden cluster among tumor cells (`subclone_1…`, `""` for non-tumor) | Largest clone is `subclone_1`; capped at `--max-subclones` (default 5) |
| `low_complexity` | bool: ambient / empty-droplet-like cell (few detected genes) | Flagged, never silently dropped. A tumor call on such a cell is untrustworthy (noise-dominated CN signal), so the classifier downgrades it to `uncertain`. Filter these out for a clean tumor set |

**On `uncertain` (flag, don't force).** Three gates downgrade untrustworthy calls to `uncertain` rather than asserting them. Two act on `tumor` calls: the **low-complexity** gate (above) and a **coherence** gate, where a high score built from scattered, single-segment spikes rather than contiguous chromosome-arm-scale runs is more likely a transcriptionally-"loud" normal cell than copy number. The third acts on `normal` calls: the **anti-alignment** gate, where a cell carrying a large, coherent deviation pointing *against* the consensus profile is not a diploid cell even though the signed score puts it next to one — most often a transcriptional outlier, occasionally a subclone whose profile opposes the one the template locked onto (see *Known limitations* in `CHANGELOG.md`). All three gates are *subtractive* — none ever promotes a cell to tumor — so they cannot inflate the tumor set. Genuinely ambiguous cells, including same-lineage normals whose expression mimics CNV (e.g. normal astrocytes vs AC-like glioma, which expression-only CNV cannot separate), land in `uncertain` by design. `n_uncertain`, `n_low_complexity` and `n_anti_aligned` are recorded in `qc.json`. The anti-alignment gate's thresholds are set relative to the malignant population's own score and burden, so it **declines to run at all when fewer than 30 cells are called tumor** — below that the median it needs is one or two cells and the gate would move with their noise. `anti_alignment_tumor_n` in `qc.json` reports how many tumor calls the estimate rested on, and is `0` when the gate was skipped.

**Caveat on `tumor_score`:** it is *genome-fraction weighted and normalized*, so unlike a raw sum it does not scale with the number of segments (`n_segments` in `qc.json`). It does scale with the amplitude of the sample's own consensus profile: a cell carrying the full consensus scores roughly the consensus's own mean |deviation|, so a sample with ±0.5 events puts its tumor mode near 0.5 and a sample with ±0.2 events near 0.2. There is still **no fixed unit or universal cutoff**; a value is only meaningful *relative to other cells in the same run*. Normal cells sit near 0 and tumor cells form a separate higher cluster, but the numbers are **not comparable across samples**. For any cross-sample comparison, rely on `class` / `confidence` or `n_segments_altered` instead of the raw score.

**The score is signed, and the sign matters.** The two ends of the range are not "most tumor" and "most normal" — they are "aligned with the clone" and "aligned against it". A strongly negative cell has a *large* real deviation pointing the opposite way, which usually means a transcriptionally extreme cell type (erythrocytes, platelets), not a confidently diploid one. Cells closest to diploid sit near **zero**, not at the minimum. If you are selecting a diploid reference set, rank on `cn_burden` (or `|tumor_score|`), never on `tumor_score` itself. The anti-alignment gate catches the sharpest version of this — a strongly negative cell whose deviation is also large and chromosome-arm-scale is labelled `uncertain` rather than `normal` — but the gate is deliberately conservative, so ranking on the signed score is still the wrong move for anything milder.

A good triage filter is to cross-tab `tumor_score` against `n_segments_altered`: confident tumor cells score high on *both*, while a high score with few altered segments usually signals a noisy or transcriptionally extreme cell rather than a true malignant one.

## Visualize the calls (inferCNV-style heatmap)

`plot-heatmap` renders a genome-ordered CNV heatmap from a completed run's
out-dir, in the familiar inferCNV layout: up to three stacked panels.
**References** (confident-diploid cells) on top, **Observations** (tumor cells,
grouped by subclone) below, and an **Uncertain** panel of flagged cells when any
exist. Genome runs left-to-right with chromosome widths proportional to gene
content; a blue-to-white-to-red diverging colour maps loss to diploid to gain.

The `normal` class on a real run is a mixed bag (genuinely-diploid cells, ambient
junk, and tumor cells under-called as normal), so plotting all of it as
"references" reads as noisy structure. Two curation steps fix that, matching how
inferCNV uses a hand-picked reference:

- **Low-complexity (ambient) cells are dropped** from every panel; they are not
  trustworthy CNV signal. Pass `--show-low-complexity` to keep them.
- **The reference panel (and the baseline the whole figure is recentered
  against) is the *confident-diploid* subset**: non-junk normals whose
  `tumor_score` is in the lowest `--ref-score-quantile` (default 0.5) of that
  pool. This keeps under-called tumor cells out of the reference frame, which
  also sharpens the observation panel. `--ref-score-quantile 1.0` keeps every
  non-junk normal.

```bash
# needs the plotting extra: pip install -e ".[plot]"
kopya plot-heatmap --run-dir results/ --out results/sample-01_cnv.png
```

Two colour scales:

- `--scale log` (default): per-segment deviation in log-space, centered on 0.
- `--scale linear`: copy-number **ratio** centered on 1.0 (diploid), matching
  inferCNV's "Modified Expression" 0.8-1.2 convention. Note this is *not* the
  same frame as `chr_cnv_matrix.csv`, despite both reading "centered on 1.0":
  the heatmap subtracts each cell's own median and the CSV does not.
  `chr_cnv_matrix_centered.csv` (`--centered-chr-matrix`) is the chromosome-level
  equivalent.

On a glioblastoma sample the tumor panel shows the textbook GBM signature:
solid **chr7 gain** and **chr10 loss**, plus chr12/18/19/20 gains, over a clean
confident-diploid reference panel.
(A residual note: any structure left in the reference panel is real biology/QC,
under-called tumor or same-lineage normals that mimic CNV, not a plotting
artifact. Curating the reference with `--norm-cell-names` removes it entirely.)

> **How the signal is recovered.** The stored `cn_per_segment.npz` is the
> per-segment mean of the median-centered smoothed signal; it carries a per-gene
> positive pedestal and a per-cell magnitude offset, so plotted raw it reads as an
> all-gain wash. `plot-heatmap` recovers the signed gain/loss signal the way
> inferCNV displays it: subtract each segment's median over the confident-diploid
> reference (references land near 0), then subtract each cell's own genome-wide median
> (removes the per-cell offset), then **denoise** by soft-thresholding each cell's
> per-segment deviation toward 0 by `--sd-amplifier` reference-cell spreads
> (default 1.0), collapsing the normal noise floor to white while leaving real CNV
> events intact. The band is estimated from the *cleanest* reference cells so
> residual noisy cells don't widen it and erase weak tumor signal. Set
> `--sd-amplifier 0` to disable denoise and see the raw reference-relative signal.

Common knobs: `--ref-score-quantile` (reference strictness, default 0.5),
`--show-low-complexity` (keep ambient cells, dropped by default),
`--sd-amplifier` (denoise strength, default 1.0; 0 disables),
`--vlim` (colour clip, default 0.15 log-units), `--n-ref` /
`--max-obs` (display subsampling caps for the reference / observation panels),
`--sample` / `--subtitle` (header text), `--dpi`. Cells are ordered within each
panel by hierarchical clustering. See
`kopya plot-heatmap --help` for the full surface.

## Validation

We test at two levels. Each level is self-contained and can be run independently.

### Level 1: Synthetic unit tests (always run in CI)

Every algorithm module has a dedicated test file with synthetic fixtures that have known answers. The end-to-end test uses a 500-cell AnnData with three planted CNV events:

- **chr7 broad gain** (200 genes × 2.5×): EGFR amplification, the textbook GBM hallmark
- **chr10 broad loss** (200 genes × 0.4×): PTEN deletion, the other GBM hallmark
- **focal MYC amp** (~30 genes on chr8, × 5.0): a focal event to probe sub-chromosomal resolution

The test asserts ground-truth recovery: tumor/normal recall ≥ 80% both directions, chr7/chr10 signals clearly separated, focal MYC visible.

```bash
pytest                          # full unit suite, ~2 min
pytest tests/test_segment.py -v # just segmentation
pytest tests/test_cli.py -v     # includes CellRanger format tests
```

These run on every push and pull request via CI. **If this passes, the pipeline recovers the planted events on the synthetic fixture** (a regression gate, not a proof of correctness).

The CLI tests cover all four input forms end-to-end on the same 500-cell fixture: `--anndata`, `--cellranger-dir`, `--cellranger-h5`, and `--counts/--obs/--var`. The CellRanger fixtures are built programmatically; no 10x Genomics tools required.

Regenerate the fixture only when intentionally changing the planted events:

```bash
python -m tests.fixtures.build_tiny_simulated
```

### Level 2: External gold-standard datasets

We test against publicly available datasets whose CNV ground truth is established independently, not derived from this tool. The tests skip automatically when the data files are absent.

**Current results:**

| Dataset | Source | Ground truth | Result |
|---|---|---|---|
| Patel 2014 GBM | GSE57872 (SMART-seq2) | Published chr7 gain / chr10 loss | 3 PASS, 4 XFAIL¹ |
| DCIS1 breast cancer | GSE148673 (CopyKAT paper) | Bulk WGS, digitized from Gao et al. 2021 (arm-level) | 5 PASS |
| SCEVAN synthetic matrices | Zenodo 6628423 | Planted tumor/normal labels | 3 PASS |
| 10X ovarian scFFPE (HGSOC) | 10x Genomics 17k dataset | FLEX cell-type annotations (expression-derived) | 4 PASS |
| Maynard 2020 lung | `maynard2020_3k` (infercnvpy) | Cell-type labels (expression-derived) + kopya-vs-inferCNV | 6 PASS |
| UCSF osteosarcoma T1 | osteosarc.com (IPISRC044_T1) | Matched bulk WES (CNVkit) | 5 PASS |

¹ The Patel SMART-seq2 dataset triggers a known baseline inversion: the unsupervised UCell cascade misfires on full-length read counts, so the four hallmark-direction tests are `xfail(strict)`. Fix path: supervised mode with cell-type annotations. Full diagnosis in [`tests/external/GOLD_STANDARD_TESTING.md §6.1`](https://github.com/ValiusSciences/kopya/blob/main/tests/external/GOLD_STANDARD_TESTING.md).

**What "ground truth" means here.** Matched DNA-level CNV truth exists only for DCIS1 and the osteosarcoma sample (each a single patient; the DCIS1 profile is digitized at chromosome-arm resolution from a published figure). The ovarian and Maynard checks compare kopya's calls against expression-derived cell-type annotations, so they measure agreement between two expression-based classifications, not DNA-level CNV accuracy. No head-to-head *accuracy* comparison against CopyKAT or SCEVAN was run; the only cross-tool check is kopya-vs-inferCNV per-chromosome concordance on Maynard.

**To reproduce:**

```bash
# One-time: acquire and convert each dataset (exact commands in the doc below)
# See tests/external/GOLD_STANDARD_TESTING.md §2 for download URLs and conversion scripts.

# Then run:
pytest tests/external/ -v
```

Full documentation (acquisition commands, acceptance criteria, per-metric results, known limitations, improvement history):
[`tests/external/GOLD_STANDARD_TESTING.md`](https://github.com/ValiusSciences/kopya/blob/main/tests/external/GOLD_STANDARD_TESTING.md)

## Documentation

- [`docs/tutorials/seurat.ipynb`](https://github.com/ValiusSciences/kopya/blob/main/docs/tutorials/seurat.ipynb): executed end-to-end tutorial: annotate in Seurat, run kopya, bring the CNV calls back (with a rendered heatmap; previews on GitHub)
- [`notebooks/working_with_results.ipynb`](https://github.com/ValiusSciences/kopya/blob/main/notebooks/working_with_results.ipynb): a hands-on walkthrough of every output file: what each column means, the numeric scale of each value, and how to re-center the CN matrix into gain/loss
- [`examples/tiny_simulated/`](https://github.com/ValiusSciences/kopya/tree/main/examples/tiny_simulated): a complete, committed result set (500-cell demo) to explore or diff against
- [`docs/comparison.md`](https://github.com/ValiusSciences/kopya/blob/main/docs/comparison.md): detailed positioning vs CopyKAT, SCEVAN, inferCNV, CONICSmat, Numbat, CaSpER
- [`tests/external/GOLD_STANDARD_TESTING.md`](https://github.com/ValiusSciences/kopya/blob/main/tests/external/GOLD_STANDARD_TESTING.md): external dataset acquisition, acceptance criteria, results, known limitations
- `kopya run --help`: full CLI surface with every knob
