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. Drop an inversion-prone label (e.g. Fibroblast on a mesenchymal tumor) or add one (e.g. Plasma_cell where plasma cells are reactive infiltrate). 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. 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.

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

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_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). Two gates downgrade untrustworthy tumor calls to uncertain rather than asserting them: 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. Both gates are subtractive (they never promote 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 and n_low_complexity are recorded in qc.json.

Caveat on tumor_score: it is a sum across segments, so its absolute magnitude scales with the number of segments (n_segments in qc.json) and the strength of the CN signal. There is no fixed unit or universal cutoff; a value is only meaningful relative to other cells in the same run. Normal cells sit near 0 (or a low baseline) and tumor cells form a separate higher cluster, but the actual numbers (e.g. a tumor cluster at 5-40) are not comparable across samples. For any cross-sample comparison, rely on class / confidence or n_segments_altered instead of the raw score.

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 and the 1.0-centered chr_cnv_matrix.csv output.

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