CNery

breseq copy-number-variation extension. CNery reads per-reference coverage tables produced by breseq bam2cov and predicts copy-number variation (CNV) across the genome. Predictions are corrected for coverage biases introduced by sequencing chemistry (GC-content bias) and prokaryotic replication state during DNA isolation (origin-to-terminus / OTR bias).

Recent updates (latest commits):

• Coverage tables are the only input — CNery reads breseq bam2cov coverage tables and nothing else. It no longer needs a BAM or a reference FASTA, and never runs breseq itself. The reference sequence it needs for GC content is already in the table's own ref_base column.
• CSV or TSV, full or --total-only — all four shapes of table are read without being declared. The delimiter is detected from the file's own header row, and the column schema from its column names. --total-only tables are the recommended input: they carry three coverage columns instead of eight, so they are markedly smaller, and CNery loses nothing by using them.
• Files and folders on the command line — name coverage tables directly, or name folders and let CNery find them by file ending (--file-ending, defaults coverage.csv, coverage.tsv and coverage.tab). Files and folders can be mixed in one command.
• Multi-genome CNV analysis — CNery processes all the coverage tables given in one pass. Each reference (chromosome, plasmid, contig, etc.) is preprocessed separately, pooled for a shared LOWESS GC-bias fit, and then bias-corrected and CN-called independently.
• Cumulative GC skew origin/terminus prediction — CNery now also locates the replication origin and terminus from the reference's own cumulative GC skew (Grigoriev 1998), independently of read depth. Reported in GC_skew/ as a marked-up plot and a JSON summary, and used as the fallback origin/terminus when the coverage profile carries no usable gradient of its own.
• Output flexibility — output prefix defaults to CNV_out/ in the current folder. Output subfolders (CNV_plt/, CNV_csv/, GC_bias/, OTR_corr/, GC_skew/) are created automatically.
• Modular bias correction — the --bias flag lets you choose all (GC + OTR), gc, otr, or none.
• A per-base segment-length prior — --change-rate is the probability per base that copy number changes, so re-tiling a genome with different -w/-s does not restate the biology. Read 1/rate as the expected segment length.
• Boundaries priced by fractional change — a boundary inside an already-amplified region is both rarer (--interior-change-rate) and, when the two states are close in ratio, a claim the coverage cannot support (--fold-change-penalty, charged per doubling as BETA / |log2(l/k)|). Without this, deep coverage resolves 126 copies from 139 on counting noise alone and reports one amplification as a segment plus a shoulder.
• Pip-installable package — requirements.txt and a fixed pyproject.toml allow install directly from GitHub via pip install git+....

Recommended: create a conda/mamba environment from the provided spec.

mamba env create -f environment.yml
mamba activate CNery

Install CNery (a.k.a. breseq-ext-cnv) from GitHub:

pip install git+https://github.com/barricklab/breseq-ext-cnv.git

CNery takes coverage tables — see Generating a coverage table if you do not have them yet. Point it at the folder holding them:

CNery <folder> [-o <output folder>] [-w <window>] [-s <step size>] [-f <fragment length>]

Run with no arguments at all and it reads the current folder (CNery --version prints the installed version). You can also name tables directly, or mix files and folders in one command:

CNery REL606.coverage.csv pPlasmid.coverage.csv -o CNV_out
CNery coverage/ extra/pContig.coverage.csv -o CNV_out

Folders are searched, top level only, for files ending in coverage.csv, coverage.tsv or coverage.tab — all three are found by default, and they may sit side by side. (.tab is the legacy extension breseq's deprecated --table flag writes; its contents are ordinary TSV.) Use --file-ending if your tables are named otherwise; repeat the flag to accept several. Note that any --file-ending replaces the defaults rather than adding to them:

CNery coverage/ --file-ending cov.txt
CNery coverage/ --file-ending cov.txt --file-ending coverage.csv

A draft genome arrives as hundreds of contigs, and CNery calls copy number per reference: each contig refits its own baseline, so every one of them comes out at copy number 1 however amplified it really is. --group-table says which coverage tables are contigs of one molecule (breseq's -c) so that they share a single background coverage distribution.

CNery --group-table 08_mutation_identification/reference_groups.tsv -o CNV_out

The table has two columns, file and group; any other column is ignored, so it can carry provenance. file names a coverage table exactly — no stem matching — resolved relative to the table's own folder. A blank group means the sequence stands alone.

file	group	length	source
REL606.coverage.tsv		4629812	REL606.gbk
pKAN.coverage.tsv		5988	pKAN.gbk
NZ_MJGT01000001.coverage.tsv	pine3-8-t1	412887	pine3-8-t1.fasta
NZ_MJGT01000002.coverage.tsv	pine3-8-t1	98014	pine3-8-t1.fasta

For a group of two or more sequences:

• they share one baseline depth, one dispersion and one copy-number grid, so a contig at three times the assembly's level is called copy number 3;
• origin-to-terminus correction is declined and GC skew is not measured — contig order and orientation are unknown, so there is no coordinate for a replication ramp or an origin to sit along. Both still write their JSON, saying so;
• "Relative copy number" is measured across the whole group, and the longest group anchors the run at 1.0. A finished chromosome beside a longer draft assembly will therefore no longer read exactly 1.0 — it reads its level relative to the assembly.

A sequence in no group, and every sequence when --group-table is not given, behaves exactly as before. A table in which nothing is grouped is a no-op, which is what lets a caller pass the flag unconditionally.

The table can also supply the inputs. With --group-table and no INPUT argument, the coverage tables read are the ones the table names — which is how breseq keeps its command line to one line on a draft assembly rather than repeating hundreds of paths. Given both, they must agree exactly: every input named by a row, every row naming an input.

Naming a folder is another way to keep the command line short. breseq's 08_mutation_identification/ holds one *.coverage.tsv per contig — hundreds of them on a draft genome — beside its own *.coverage.tab files, which share a default ending but not the schema and would be rejected by name. Restricting the ending picks out exactly the tables CNery wants, in one argument instead of hundreds:

CNery 08_mutation_identification/ --file-ending coverage.tsv -o 09_copy_number_variation/cnery_out

A table's sequence ID comes from its file name, with the matched ending and the . in front of it removed — REL606.coverage.csv becomes REL606, and NC_012967.1.coverage.tsv becomes NC_012967.1. That ID names every output file, and no two inputs may share one.

Everything passed in one command is analyzed together, sharing a single GC-bias fit and one global coverage median. That is what you want for the references of one sample — chromosome, plasmids and contigs. Analyze separate samples with separate commands.

Calculate coverage with a 500 bp window sliding in 250 bp steps; sequencing fragment length is 300 bp:

CNery <inputs> -o <output folder> -w 500 -s 250 -f 300

Analyze coverage across the whole genome, but restrict the CNV plot to a specific genomic segment:

CNery <inputs> -o <output folder> --region REL606:3497890-3955678 -w 1000 -s 500

The sequence ID is the one derived from the table's file name. SEQ_ID: may be omitted when the run has only one input sequence:

CNery REL606.coverage.csv -o CNV_out --region 3497890-3955678

Repeat the flag to plot several sequences, at most once each:

CNery coverage/ -o CNV_out --region REL606:3497890-3955678 --region pPlasmid:1-40000

Open intervals work too — REL606:3497890- runs to the end of that sequence, REL606:-3955678 from its start.

Two things to know. Giving any --region also selects which sequences are plotted: a sequence not named gets no CNV plot. And --region affects plotting only — coverage, bias fitting and copy-number calling always cover every sequence, and the output CSVs always contain every window for every sequence, plotted or not.

Control which bias correction is applied before CN prediction:

# Both GC + OTR corrections (default)
CNery <inputs> -o <output folder> -w 500 -s 250 --bias all

# Only correct OTR (replication) bias
CNery <inputs> -o <output folder> -w 500 -s 250 --bias otr

# Only correct GC-content bias
CNery <inputs> -o <output folder> -w 500 -s 250 --bias gc

# No bias correction before CN prediction
CNery <inputs> -o <output folder> -w 500 -s 250 --bias none

# Polymorphism mode: call a continuous copy number instead of an integer one
CNery <inputs> -o <output folder> -p

# ...at a finer grid (levels 1.02, 1.04, ... instead of 1.05, 1.10, ...)
CNery <inputs> -o <output folder> -p --copy-number-resolution 0.02

By default CNery calls an integer copy number: a region is at one copy, or two, or thirty-four. That is the right answer for a clonal isolate, and it is the mode breseq consumes.

-p/--polymorphism-mode instead calls a continuous level, so a region can come back at 1.30. Read that as a measurement of relative depth and nothing more — it does not say why. A level of 1.30 is consistent with 30% of the population carrying a duplication, with a mixed sample, or with aneuploidy, and CNery does not claim to distinguish them; there is no cell-fraction parameter, deliberately. This is the same stance as the "Relative copy number" already reported per sequence in OTR_corr/*_otr_results.json, which is likewise non-integral on purpose.

--copy-number-resolution sets how fine the grid is (default 0.05, so the smallest change above single copy that can be called is 1.05). It is rounded so that single copy always lands exactly on the grid. Finer is not automatically better: the default is roughly the smallest difference a 10 kb event can support at the default windowing, and a grid finer than the coverage's own scatter will simply cost sensitivity.

Two practical notes:

• Consensus-mode output is completely unchanged by this feature, bit for bit.
• -p output is not readable by breseq. The State column of CNV_csv/*_break_pts.csv becomes fractional, and breseq expects an integer copy number there. Use the default mode when the results are destined for breseq.
• A level one step off single copy is usually residual bias, not biology. On real data the GC and replication corrections leave a slow wander of a few percent, and a fine grid resolves it into long stretches at 0.95 or 1.05. Measured across the test corpus, rounding the continuous calls reproduces the ordinary integer call on 99.6–100% of windows — the mode refines calls rather than changing them — so treat a one-step departure as the resolution limit and look for events that stand clear of it.

When OTR correction is applied, the origin and terminus of replication are automatically inferred — no manual coordinates are required. CNery fits them from the coverage profile, falling back to the reference's own cumulative GC skew when the coverage carries no usable gradient. Either way the correction is applied only if it beats a bootstrap null, so a flat genome is left alone.

Coverage tables are CNery's only input. Generate them once with breseq bam2cov, then run CNery against them as often as you like — no BAM, no reference FASTA, and breseq need not be installed on the machine that runs CNery. The reference sequence CNery needs for GC content is already in each table's ref_base column.

The conventional layout is one table per reference sequence, named <seq_id>.coverage.csv (or .coverage.tsv — both are found by default):

coverage/
├── REL606.coverage.csv
└── pPlasmid.coverage.csv

CNery coverage/ -o CNV_out

Requires pre-release breseq, the Barrick lab channel of development builds auto-built from barricklab/breseq master. Released versions on bioconda do not include the current fixes to how coverage tables are written.

conda install -c https://barricklab.github.io/conda/ -c conda-forge -c bioconda breseq-prerelease

Omit --region and breseq writes one table per reference sequence automatically — no need to look up sequence IDs or lengths. --output is then a directory:

mkdir -p coverage
breseq bam2cov \
  --format CSV \
  --total-only \
  --resolution 0 \
  --output coverage \
  -b data/reference.bam \
  -f data/reference.fasta

To generate a single table instead, name the region using the sequence ID exactly as it appears in the FASTA, and give --output a name ending in .coverage.csv:

breseq bam2cov \
  --format CSV \
  --total-only \
  --region REL606:1-4629812 \
  --resolution 0 \
  --output coverage/REL606.coverage.csv \
  -b data/reference.bam \
  -f data/reference.fasta

Four options matter:

• --resolution 0 outputs every position. The default is 600, which samples only 600 points across the whole region — far too sparse for windowed coverage, and it fails silently by producing a well-formed but nearly empty table.
• --format CSV produces a comma-separated table instead of a plot; TSV gives the same columns tab-separated. CNery reads either and works out which from the file itself, so the choice is yours. (The old -t / --table flag is deprecated: it still selects TSV, but writes the legacy .tab extension. CNery matches that too, but prefer --format TSV.)
• --total-only (short flag -1) writes unique_cov, redundant_cov and total_cov in place of the eight strand-split coverage columns — roughly 2.5× smaller files. breseq sums the strands itself, and those sums are exactly what CNery computes from the wider table, so nothing CNery uses is lost: repeat detection via redundant_cov still works, and copy-number calls are unchanged. Recommended.
• --per-read-group (optional) repeats every coverage column once per read group (@RG) in the BAM, prefixed RG-<n>_ where <n> is the read group's index in the BAM header. Requires a table format; it is rejected with --format PNG. A BAM with no read groups yields a single RG-0 set.

breseq bam2cov \
  --format CSV \
  --total-only \
  --per-read-group \
  --resolution 0 \
  --output coverage \
  -b data/reference.bam \
  -f data/reference.fasta

CNery reads such a table without any change. The aggregate columns keep their names and positions and the per-read-group repeats are appended, so the columns CNery needs are still where it expects them; it selects by name and ignores the rest. The per-read-group summary lines added to the footer are stripped along with the others by their # prefix. Nothing in CNery consumes the per-read-group columns today — they pass through — so the option is safe to enable now if you want per-library coverage available in the same file for other tools.

The output carries a header row, one row per reference position, and a trailing #-commented summary block. With --total-only --format CSV:

position,ref_base,unique_cov,redundant_cov,total_cov
1,G,32,0,32
2,G,32,0,32
...
#,region_unique_average_cov,56.7267
#,region_repeat_average_cov,0
#,region_average_cov,56.7267
#,number_of_positions,4629812

Without --total-only, each count is split by strand and six further columns follow:

position	ref_base	unique_top_cov	unique_bot_cov	redundant_top_cov	redundant_bot_cov	...
1	G	14	18	0	0	...

The delimiter is the only difference between --format CSV and --format TSV; it is used for the header, the data and the footer alike.

That summary block is variable in length — --show-average adds a line, and --per-read-group adds three per group (# RG-0_region_unique_average_cov …) — so anything parsing these tables should skip lines by their # prefix rather than dropping a fixed number from the end.

For the same reason, do not assume a fixed column count: position must be the first column and the named coverage columns must be present, but extra columns to the right are expected and should be ignored rather than treated as an error.

CNery needs position, ref_base, and unique-versus-redundant coverage in one of the two shapes above — unique_cov + redundant_cov, or the four strand-split unique_*/redundant_* columns. It checks for them as soon as a table is opened, so a file with the wrong schema is rejected by name instead of failing later. total_cov is ignored: it is unique_cov + redundant_cov by construction. Note that breseq's own 08_mutation_identification/*.coverage.tab files use a different schema (position last, no ref_base) and are not usable as CNery input.

If you do need the sequence IDs and lengths for a --region argument, read them from the FASTA headers or the .fai index:

grep '^>' data/reference.fasta
cut -f1,2 data/reference.fasta.fai

Given an output folder CNV_out/, CNery writes:

The -f fragment size is optional. GC bias acts at the scale of the sequenced fragment, so that size belongs to the library rather than to the analysis — and it is not something a coverage table shows you. Left unset, CNery scores candidate sizes by how well the GC each implies predicts held-out coverage (with the replication ramp divided out and copy-number variants excluded, so neither can be mistaken for a GC effect) and reports what it chose. The 400 bp default is kept unless a candidate beats it by more than the measurement's own error. Pass -f to pin it.

CNery corrects coverage and calls copy number in two passes. The first runs GC correction, origin-to-terminus correction and the HMM as usual. The second repeats both fits with every window the first pass did not call single-copy excluded from them, then calls copy number again — only the second pass's results are written.

An amplification is invisible to the crude censoring the first pass has available (near-zero depth, repeat overlap), so it otherwise sits in the GC and ramp fits at full weight and distorts them. Measured on synthetic coverage where the truth is known, a real 1.5x replication ramp with an amplification on top is detected 0 times in 60 without this and 60 times in 60 with it, at no cost in false positives. If excluding the non-single-copy windows would leave under half a sequence, it is skipped for that sequence and reported as such.

• CNV_out/CNV_plt/ — per-reference CNV prediction plots.

• CNV_out/CNV_csv/ — per-window coverage + CN calls as CSV.

• The GC correction is a fitted curve, not an exact quantity, so how well it is determined at each window's GC is measured (by resampling the fit) and carried into the copy-number model as extra variance. The effect grows with copy number, because a correction factor's error is multiplied by the number of copies — which is why a window at an extreme GC inside an amplification is no longer able to earn its own copy-number segment on the strength of the correction alone.

• CNV_out/GC_bias/ also holds GC_passes.pdf — the GC correction is fitted in two pooled passes, once on raw coverage and again after OTR correction (which reintroduces a GC trend, because the replication ramp varies with position and position correlates with GC). The second pass additionally excludes every window the first pass's copy-number calls did not put at CN=1. The plot shows both curves and their product, which is the total correction actually applied.

• CNV_out/GC_bias/ — pooled LOWESS GC-bias diagnostic plot, GC_vs_NormRds.pdf. Both files here cover the whole run rather than one reference, which is why neither is named for a sequence.

• CNV_out/corr_plots/ — per-reference before/after diagnostic (*_correction_stages.pdf): one row for each correction step — GC and OTR, in each of the two passes — showing coverage before and after it, with the fitted curve overlaid, and directly beneath each row a strip of which windows that particular fit was allowed to see. Deletions are drawn as spans, repeats as a density track, and the second pass's strips additionally mark everything called CN≠1. One row does not continue from the one above it, and says so: the second OTR fit divides the first pass's ramp back out, because a ramp has to be fitted to coverage that still contains it. Each row is labelled with the fraction of windows within 20% of single copy before and after, reported both over all windows and over uncensored windows only — the two can differ a lot on a repeat-heavy replicon, and the strip below shows why. Produced in every --bias mode.

• CNV_out/OTR_corr/ — per-reference OTR bias plots and a JSON summary (*_otr_results.json) containing the inferred origin window, terminus window, normalized coverage at each, the origin-to-terminus ratio, and the sequence's relative copy number.

  This file is written for every reference in every --bias mode, including ones where no correction was attempted and ones with no usable coverage at all. A reference whose table has no position rows, or whose every window reads zero, is reported rather than treated as an error: it gets the ratio "Not detected", a "No usable coverage reason" saying which of the two it was, and its own (empty) CSVs and plots, and the run continues with the other references and exits 0. A plasmid that got no reads does not cost you the chromosome sequenced alongside it.

  It also reports how well the applied ramp actually fits, beyond r²: "Residual structure score" asks whether what the tent failed to explain is systematically structured rather than just noisy — a fit can explain a lot of variance and still be the wrong shape over a long stretch. It is a z-score against a bootstrap null, so read it as roughly: below 1 unstructured, 1–2 mild, above 2 structured. "Residual decorrelation length (bp)" is the scale of correlation that would be present anyway, published alongside so the score can be judged in context. Both are diagnostic only — nothing in CNery acts on them.

  It also records the evidence behind the decision, whether or not a correction was applied: "Coverage fit r-squared" / "Coverage fit p-value", the same pair for the GC-skew-anchored fit, the "Coverage vs skew likelihood ratio" and its p-value when both candidates were live, "Bootstrap surrogates", and "Breakpoint source" (coverage fit, GC skew, or not corrected). A rejected fit is therefore diagnosable from the file alone. Note the p-values are floored at 1/(surrogates+1), so 0.001 is an upper bound rather than a measurement.

  "Relative copy number" is that sequence's coverage relative to the longest sequence in the run, which reads exactly 1.0. It is not rounded to an integer: a plasmid at 2.95 is a measurement, and it is the only place plasmid copy number is reported — prob_copy_number in the CSVs is called per reference, so a uniformly multi-copy plasmid comes out as 1 there.

• CNV_out/GC_skew/ — per-reference cumulative GC-skew plots with the predicted origin and terminus marked, and a JSON summary (*_gc_skew_results.json).

Each coverage table produces its own set of outputs, named with the sequence ID derived from its file name. The GC-bias plot is the exception: one pooled fit covers every table in the run.

Bacterial genomes are G-rich on the leading strand and C-rich on the lagging strand, so the sign of the GC skew (G−C)/(G+C) flips at the two points where the replication strands switch. Summing the skew along the genome turns those sign changes into extrema: following Grigoriev 1998, the cumulative curve reaches its minimum over the replication origin and its maximum at the terminus.

CNery computes this from the coverage table's ref_base column — no FASTA and no read depth involved — and writes the result for every reference, in all four --bias modes:

{
    "Origin (bp)": 3885501,
    "Terminus (bp)": 1526001,
    "Origin window index": 7771,
    "Terminus window index": 3052,
    "Windows": 9258,
    "Separation (fraction of genome)": 0.4903,
    "Cumulative skew amplitude": 151.2949,
    "Replichore skew t-statistic": 34.77,
    "Replichore skew p-value": 0.001,
    "Bootstrap surrogates": 1000,
    "Prediction confident": true,
    "Prediction method": "Ori-ter coordinates from cumulative GC skew (Grigoriev 1998)"
}

Prediction confident requires two things: the two extrema roughly antipodal (35–65% of the sequence apart, as bidirectional replication implies), and a p-value of 0.01 or better. The coordinates are reported either way — a low-confidence call stays diagnosable from the JSON and the plot rather than being reduced to a flag.

Adjacent windows are not independent — genome composition varies on scales far longer than one window — so an ordinary t-test would badly overstate significance. Replichore skew t-statistic is therefore reported as an effect size only; its magnitude is inflated by an unknown factor and should not be converted to a p-value.

The p-value instead comes from a circular block bootstrap. Contiguous blocks of windows are resampled with replacement around the circle, which preserves local autocorrelation while destroying the long-range two-arm pattern being tested for; the null is "a sequence that wobbles like this one but has no single origin". The full procedure — locating the extrema, then scoring the two arms — is re-run on every surrogate, so choosing the breakpoints by looking at the data is paid for rather than ignored.

Two things to know when reading it:

• It is floored at 1/(B+1). A real chromosome beats all 1,000 surrogates and reads back exactly 0.001. That is an upper bound, not a measurement — read it as "p < 0.001". Bootstrap surrogates is reported so the floor is visible.
• Block length adapts to sequence length. What governs power is the number of blocks rather than their size, so CNery targets ~20 blocks, bounded to 10–200 windows each. Sequences too short to give both long-enough and numerous-enough blocks simply do not reach significance, which is an honest reflection of how little evidence they carry.

The bootstrap is seeded, so repeated runs on the same input give the same p. It adds roughly 3% to the per-sequence cost (~0.15 s on a 4.6 Mb genome).

Expect false on plasmids. They have no bidirectional replication origin, so there is no sign change for the cumulative curve to turn on, and the two extrema you get back are noise — the two in the test data land at p = 0.36 and 0.44. That is the intended answer, not a failure, and chromosomes in the same run are unaffected since the prediction is made per reference.

Two properties worth knowing. The prediction depends only on the reference, so different samples aligned to the same reference give identical answers, whatever their depth or growth phase. And it is invariant to circular permutation: a reference whose coordinates start elsewhere predicts the same locus.

These values feed the OTR correction as its second candidate. CNery prefers the coverage-derived origin and terminus when they clear their own significance test and a likelihood-ratio test says they fit better than a ramp hinged at the GC-skew coordinates; otherwise it uses the skew's, provided the skew prediction is confident and the coverage does not contradict which end is the origin. "Correction type" and "Breakpoint source" in *_otr_results.json say which was used.

The practical consequence: a sequence with no replication gradient in its coverage but a confident skew prediction now receives a small correction where it previously received none. The magnitude is still fitted from the coverage, not imported — on a genuinely flat sequence the fitted ramp is close to 1.0 — but the evidence for such a correction is the reference sequence, not the reads. "GC skew fit p-value" reports what the coverage itself had to say about it.

$ CNery -h

usage: CNery [-h] [--file-ending ENDING] [--region SEQ_ID:START-END] [-o O]
             [-w W] [-s S] [-f F]
             [-z DELETION_COVERAGE_FRACTION] [--change-rate CHANGE_RATE]
             [--interior-change-rate INTERIOR_CHANGE_RATE]
             [--fold-change-penalty FOLD_CHANGE_PENALTY]
             [--max-copy-number MAX_COPY_NUMBER] [-p]
             [--copy-number-resolution CN_RESOLUTION]
             [--bias {all,none,gc,otr}]
             [INPUT ...]

CNery is a Python package extension to breseq that analyzes the sequencing
coverage across the genome to predict copy number variation (CNV).

positional arguments:
  INPUT                 Coverage table files, and/or folders containing them.
                        Folders are searched (top level only) for files
                        ending in --file-ending. Every table given is
                        analyzed together, sharing one GC-bias fit, so these
                        should be the reference sequences of a single sample.
                        Defaults to the current folder.

options:
  -h, --help            show this help message and exit
  --file-ending ENDING  File ending that identifies a coverage table inside
                        an input folder. Repeat the flag to accept more than
                        one. Any --file-ending REPLACES the defaults
                        ('coverage.csv', 'coverage.tsv', 'coverage.tab')
                        rather than adding to them. A file named directly on
                        the command line is always used, whatever it is
                        called.
  --region SEQ_ID:START-END
                        Plot the CNV calls for one sequence over a genomic
                        segment, e.g. 'REL606:3497890-3955678'. The sequence
                        ID is the one derived from the table's file name.
                        Repeat the flag to plot several sequences, at most
                        once each. 'SEQ_ID:' may be omitted when the run has
                        only one input sequence. Open intervals are accepted:
                        'REL606:3497890-' runs to the end of the sequence,
                        'REL606:-3955678' from its start. Giving any --region
                        also selects WHICH sequences are plotted: those not
                        named get no CNV plot. This affects plotting only --
                        coverage, bias fitting and copy-number calling always
                        cover every sequence, and the output CSVs always
                        contain every window.
  -o, --output O        output file prefix / storage location. Defaults to
                        the 'CNV_out' folder in the current dir.
  -w, --window W        Window length used to parse the genome and compute
                        coverage and GC statistics. Default: 100. Wider
                        windows smooth the coverage but lose short events: the
                        window statistic is a per-base median, whose precision
                        grows sublinearly with width, so -w is a resolution
                        knob.
  -s, --step-size S     Step size (<= window size) for each progression of
                        the window across the genome. Set step-size = window
                        size for non-overlapping windows. Default: 100, i.e.
                        non-overlapping. Copy-number calls are near-invariant
                        to this: the state-change prior is per base (see
                        --change-rate) and overlapping windows are
                        down-weighted so they do not count the same bases
                        twice.
  -f, --frag-size F     Average fragment size of the sequencing library. GC%
                        is measured over this many bases centred on each
                        window. Ignored when smaller than -w. Default: 400.
  -z, --deletion-coverage-fraction DELETION_COVERAGE_FRACTION
                        Coverage a deleted region still shows, as a fraction of
                        the single-copy level. Sets the mean of the
                        copy-number-0 emission. Real deletions are not empty --
                        mismapping and repeat spill leave a couple of percent
                        behind. A fraction rather than an absolute depth so
                        that what counts as a deletion does not change with how
                        deeply the sample was sequenced. Default: 0.02.
  --change-rate CHANGE_RATE
                        Prior probability PER BASE that copy number changes.
                        The per-window probability is 1 - exp(-rate * step-
                        size), so changing -w/-s no longer changes the
                        implied biology. Read 1/rate as the expected segment
                        length: the default 1e-06 is one copy-number boundary
                        per megabase. Larger values give more, shorter
                        segments.
  --interior-change-rate INTERIOR_CHANGE_RATE
                        Prior probability PER BASE of a boundary INSIDE an
                        already-altered region -- one amplified state abutting
                        another, with no return to single copy between them.
                        Same units as --change-rate and converted the same way,
                        so -w/-s does not restate it. Rarer than an ordinary
                        boundary because it takes two events rather than one.
                        Default: 1e-08, one per 100 Mb.
  --fold-change-penalty FOLD_CHANGE_PENALTY
                        Extra cost in NATS for an interior boundary, per
                        doubling: moving from copy number k to l costs this
                        divided by |log2(l/k)|. So it is 2 nats at a 2-fold
                        change, 7x that at 10%, 70x at 1% -- scale-free, so
                        1 -> 2 and 100 -> 200 cost the same. It is what stops
                        deep coverage resolving 126 from 139 and splitting one
                        amplification into a segment and a shoulder. Applies
                        only between two amplified states; transitions to or
                        from single copy and the deletion state are untouched.
                        0 disables it. Default: 2.
  --max-copy-number MAX_COPY_NUMBER
                        Highest copy number the HMM can call. This is a
                        CEILING, not the grid: the number of states is sized
                        from the data, so raising this changes nothing unless
                        the sample really does carry a segment above the old
                        value. Memory is linear in the states actually used. A
                        call that lands exactly on the ceiling is reported as
                        such, because it is a clipped value rather than a
                        measurement. Default: 500.
  -p, --polymorphism-mode
                        Call a CONTINUOUS copy number instead of an integer
                        one. The HMM decodes over a grid of coverage levels
                        spaced by --copy-number-resolution, finely between 0
                        and 2 and on the integers above, refined where a
                        segment's own data asks for it. A level is a measured
                        relative depth and asserts no mechanism: 1.3 may be a
                        subpopulation carrying a duplication, a mixed sample,
                        or aneuploidy, and this mode does not claim to tell
                        them apart. NOTE the copy number written to
                        break_pts.csv is fractional in this mode, which breseq
                        cannot read; consensus mode is unchanged.
  --copy-number-resolution CN_RESOLUTION
                        Spacing of the copy-number grid under -p, in copies:
                        0.05 means the finest level above single copy is 1.05.
                        Rounded to one that divides 1.0 exactly, so single copy
                        is always on the grid. The default is about the finest
                        step a 10 kb event supports at the default windowing.
                        Default: 0.05.
  --bias {all,none,gc,otr}
                        Select which bias correction to apply before CN
                        prediction. 'all' applies GC + OTR, 'gc' or 'otr'
                        applies only that one, 'none' skips bias correction.
                        Default: all.

Inputs are breseq 'bam2cov' coverage tables (CSV or TSV). Run with no
arguments in a folder that holds them, or name files and/or folders directly.

The test suite has two tiers, and a bare pytest runs both. Set up the development environment first:

conda env create -f dev-environment.yml --prefix=$PWD/env
conda run -p $PWD/env pytest                 # everything
conda run -p $PWD/env pytest -m synthetic    # fast, offline
conda run -p $PWD/env pytest -m authentic    # real data only

Synthetic tier — DataFrames constructed in tests/conftest.py, staged to match the pipeline's column contract at each stage. Fast, self-contained, no network. Use -m synthetic as the inner loop while editing.

Authentic tier — real breseq coverage tables published as GitHub Release assets: genuine coordinate gaps, repeat regions, per-read-group columns, and copy-number variation that synthetic frames cannot reproduce. On a cold cache the first run downloads ~105 MB, then caches it.

Running both by default is deliberate: an opt-in tier is one people forget, and real-data coverage then lapses without anyone noticing. If you want the fast path, ask for it explicitly.

Each dataset is pinned by sha256 in tests/data/registry.json, so an asset replaced in place fails the hash check rather than silently changing what the tests measure. Downloads are cached by pooch; set CNERY_TESTDATA_DIR to relocate that cache.

Datasets ship coverage tables but no BAM — coverage tables are all CNery reads, which keeps them small.

An unavailable dataset causes those tests to skip rather than fail, so offline work stays possible. That means a default run without network access reports passes and skips together — check for skips rather than assuming green means everything ran.

Adding tests, publishing datasets, and updating golden files are covered in DEVELOPER.