bio-copy-number-copy-ratio-segmentation

What this skill does

## Version Compatibility

Reference examples tested with: R 4.3+ with DNAcopy 1.76+, Python 3.10+ with numpy 1.26+, pandas 2.2+; QDNAseq 1.38+ (optional, GC/mappability normalization).

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('DNAcopy')` then `?segment` to confirm arguments
- Python: `pip show numpy pandas`

If code throws an error, introspect the installed package and adapt the example. CBS lives in Bioconductor `DNAcopy`; HMM segmentation is provided by caller-specific backends (CNVkit uses `pomegranate`; HaarSeg has its own R/Python packages).

# Copy-Ratio Segmentation

**"Turn noisy per-bin depth into clean copy-number segments"** -> Two stages, both error-prone. First, normalize the depth profile so the only remaining variation is copy number (not GC, mappability, or replication timing). Second, partition the normalized profile into segments of constant copy number. The segmentation algorithm choice has a *predictable* bias signature, and the diploid-baseline choice can invert every call.

- R: `DNAcopy::segment` (CBS, the reference implementation)
- Python: HMM via `pomegranate`; HaarSeg via `haarseg`
- The output feeds every CNV caller (cnvkit-analysis, gatk-cnv, allele-specific-copy-number)

## Stage 1: Why Depth Is Biased Before It Is Copy Number

Raw read depth confounds copy number with three systematic biases:

| Bias | Cause | Correction |
|------|-------|------------|
| GC content | PCR efficiency and probe hybridization vary with GC | Loess fit of depth vs GC (QDNAseq), or matched normal |
| Mappability | Multi-mapping reads under-counted in repetitive regions | Mappability track filter/weight; exclude low-mappability bins |
| Replication timing | Late-replicating DNA is under-represented — the "wave artifact" | Matched normal or PoN; GC correction alone does NOT remove it |
| Capture efficiency | Per-probe hybridization varies 10-100x (hybrid capture) | Panel of normals — the dominant bias for exomes/panels |

The key postdoc-level point: **GC correction alone is insufficient.** The wave artifact in cancer WGS is driven by replication timing, a biological signal GC normalization cannot flatten. Only a matched normal or a panel of normals removes it. This is why a CNVkit flat reference (GC-only) produces systematic false focal calls and why GATK tangent normalization exists.

## Stage 2: Segmentation Algorithm Taxonomy

| Algorithm | Model | Strength | Fails when |
|-----------|-------|----------|------------|
| CBS (circular binary segmentation) | Recursive t-statistic breakpoint test | High precision; excellent on small focal segments | Low depth (~3x): recall drops to ~42% (worse under over-dispersed counts); ~2 orders slower; fragments across assembly gaps |
| HMM | Hidden CN states, emission + transition | Depth-robust; high recall at low coverage | Less precise on small focal segments (~5 kb: ~76% precision vs CBS ~96%, Poisson model); EM finds only local optima |
| HaarSeg | Wavelet (Haar) multiscale edge detection | Very fast; good for shallow WGS | Less precise breakpoints than CBS; threshold-sensitive |
| Fused lasso (flasso) | L1-penalized piecewise-constant fit | Smooth; tunable sparsity | Penalty hard to set; can over-smooth focal events |
| ASPCF | Allele-specific piecewise-constant fit | Joint logR+BAF segmentation (ASCAT) | Needs BAF; see allele-specific-copy-number |

**Quantitative benchmark (Zhang et al 2024, Brief Bioinform):** the cited precision/recall numbers (CBS ~42% recall at 3x; HMM ~81% recall at 3x; CBS ~96% precision vs HMM ~76% on 5 kb focal segments under a Poisson model) summarise that paper's reported direction of the trade-off. Verify the exact figures against the published tables before quoting them in print; the qualitative trade-off (depth-vs-event-size, CBS-vs-HMM) is robust across recent benchmarks but the precise percentages depend on the simulation model (Poisson vs over-dispersed negative-binomial). There is no universally correct choice.

## Decision Tree

| Scenario | Algorithm | Rationale |
|----------|-----------|-----------|
| Panel / exome, adequate depth, focal events matter | CBS | Precise on small segments |
| Shallow WGS (< ~5x), broad events | HMM or HaarSeg | CBS recall degrades at low depth |
| Heterogeneous / impure tumor | HMM (e.g. CNVkit `hmm-tumor`) | Broader state transitions absorb noise |
| Germline, near-diploid | HMM with diploid-tight priors | Priors stabilize calls near CN=2 |
| Allele-specific (need BAF) | ASPCF / FACETS joint CBS | See allele-specific-copy-number |
| Very large WGS, speed-critical | HaarSeg | Near-linear; CBS is ~100x slower |

## Bias Correction — GC Loess Normalization

**Goal:** Remove GC-content bias from a per-bin depth profile.

**Approach:** Fit a loess curve of depth versus GC content, divide each bin by its fitted value, log2-transform. This corrects GC but not replication timing — use a normal for that.

```python
import numpy as np
import pandas as pd
from statsmodels.nonparametric.smoothers_lowess import lowess

def gc_correct(bins):
    '''GC-correct a per-bin depth profile. bins: columns chrom, start, depth, gc.
    Returns log2 copy ratio relative to the GC-corrected genome median.'''
    df = bins[(bins['depth'] > 0) & bins['gc'].between(0.3, 0.7)].copy()
    fitted = lowess(df['depth'], df['gc'], frac=0.3, return_sorted=False)
    df['corrected'] = df['depth'] / fitted
    df['log2'] = np.log2(df['corrected'] / df['corrected'].median())
    return df
```

For exomes and panels, a panel of normals (per-bin median of normals, or PCA denoising) is preferred over GC-only correction because it also removes capture and replication-timing bias.

## Segmentation — CBS with DNAcopy

**Goal:** Segment a normalized log2 profile into copy-number regions.

**Approach:** Build a CNA object, smooth single-bin outliers, run CBS, then merge adjacent segments whose means differ by less than a noise-scaled threshold (`sdundo`).

```r
library(DNAcopy)

# bins: data frame with chrom, maploc (bin midpoint), log2
cna <- CNA(genomdat = bins$log2, chrom = bins$chrom, maploc = bins$maploc,
           data.type = 'logratio', sampleid = 'tumor')
cna <- smooth.CNA(cna)                          # damp single-bin outliers

# alpha = breakpoint significance; undo.splits='sdundo' merges segments whose means
# are within undo.SD noise standard deviations -- the main guard against oversegmentation.
seg <- segment(cna, alpha = 0.01, undo.splits = 'sdundo', undo.SD = 2,
               verbose = 1)
write.table(seg$output, 'tumor.segments.tsv', sep = '\t',
            quote = FALSE, row.names = FALSE)
```

## Failure Modes

### Oversegmentation / hyperfragmentation

**Trigger:** CBS `alpha` too liberal, `undo.SD` too small, or a noisy (high-MAD) profile; FACETS `cval` too low.

**Mechanism:** The breakpoint test fires on noise; the profile shatters into many tiny segments that do not correspond to real copy-number changes.

**Symptom:** Hundreds of short segments; segment count scales with noise, not biology; downstream integer CN incoherent with per-bin medians.

**Fix:** Raise `alpha` toward 0.01 or stricter, increase `undo.SD` (e.g. 2-3), or denoise the input first (better PoN, drop low-coverage bins). Three signatures (Steele 2022) had to be discarded as oversegmentation artifacts — fragmentation propagates into every downstream analysis, including copy-number signatures.

### The diploid-baseline centering trap

**Trigger:** Centering the log2 profile on its median or mode in a hyper-aneuploid or whole-genome-doubled genome.

**Mechanism:** Centering assumes the commonest log2 value is diploid. In a WGD genome the commonest state is tetraploid; centering on it shifts the whole profile so true diploid regions read as deletions and amplifications read as neutral.

**Symptom:** Genome-wide gain or loss inconsistent with biology; segmentation is fine but every call has the wrong sign.

**Fix:** Do not dept