Bio Copy Number Copy Ratio Segmentation

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:

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

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

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.

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

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 depth-center aneuploid genomes. Anchor the diploid baseline with BAF/SNV data via an allele-specific caller, which estimates absolute ploidy. GISTIC and most callers require a correctly centered seg file as input.