kopya

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. 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:

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:

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:

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). kopya also scales to ~800k cells in about 4 minutes at ~15 GB, beyond where the R tools stay practical on commodity hardware.

Install

# 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:

# 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:

# 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/ for a pre-computed copy of these outputs and 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

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:

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.

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.

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.
• --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), 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

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

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? Run on an AnnData and the results are written back onto it, no files, no CLI:

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 (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.
• 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.

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.

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'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; 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).
• 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). 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:

New to the output files? Start with the worked example and notebook: examples/tiny_simulated/ holds a full set of results from a tiny 500-cell run, and 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:

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.
• {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:

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

{
  "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:

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.

# 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.

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:

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:

¹ 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.

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:

# 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

Documentation

• 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: 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/: a complete, committed result set (500-cell demo) to explore or diff against
• docs/comparison.md: detailed positioning vs CopyKAT, SCEVAN, inferCNV, CONICSmat, Numbat, CaSpER
• tests/external/GOLD_STANDARD_TESTING.md: external dataset acquisition, acceptance criteria, results, known limitations
• kopya run --help: full CLI surface with every knob