Bio applied copy number analysis

DNA Copy Number Analysis

When to Use

• Turning binned read-depth (from WGS/WES/targeted panels) into log2 ratios and segmented CNV calls
• Implementing or explaining circular binary segmentation (CBS) from first principles
• Converting a segment's mean log2 ratio into an integer copy-number state
• Correcting log2 ratios for tumor purity/ploidy in cancer samples
• Annotating CNV segments with overlapping genes for driver-gene interpretation

Version Compatibility

NumPy ≥1.24, SciPy ≥1.11, pandas ≥2.0, Python ≥3.10. No specialized CNV caller required — this is a from-scratch, tool-agnostic workflow; for production pipelines see CNVkit or GATK CNV (below).

Prerequisites

• pip install numpy scipy pandas
• Familiarity with read-depth/coverage concepts and BED-style genomic intervals
• Related skill: bio-genome-intervals-coverage-analysis for generating the raw per-bin depth

Copy Number Reference

Type	CN	Log2 ratio	Effect
Homozygous deletion	0	-inf	Gene loss
Heterozygous deletion	1	-1.0	LOH if other allele mutated
Normal diploid	2	0.0	Baseline
Single-copy gain	3	+0.585	Increased gene dosage
Amplification	4+	+1.0+	Common in oncogenes (ERBB2, MYCN)

Tumor purity: at purity p, observed CN = p*tumor_CN + (1-p)*2, which attenuates the log2 signal toward 0 as purity drops. Production tools (CNVkit, Sequenza, PURPLE) jointly estimate purity and ploidy from the data; see ## Purity Correction below for the arithmetic.

Depth Normalization

Goal: convert raw per-bin read depth into a log2 ratio centered on 0 for diploid regions. Approach: divide by a local rolling median (corrects regional GC/mappability bias), rescale so the genome-wide median is 1.0, then log2-transform.

import numpy as np
from scipy.ndimage import median_filter

def normalize_depth(depth, window=50):
    """Normalize read depth by the local rolling median to correct GC/mappability bias.

    Parameters
    ----------
    depth : np.ndarray
        Raw per-bin read depth (one value per genomic bin).
    window : int
        Size of the rolling-median window, in bins.

    Returns
    -------
    np.ndarray of normalized depth, median-rescaled to ~1.0 for diploid bins.
    """
    local_median = median_filter(depth.astype(float), size=window, mode="reflect")
    norm = depth / np.where(local_median > 0.1, local_median, 0.1)
    return norm / np.median(norm)

depth = np.random.default_rng(0).poisson(50, 500).astype(float)
log2_ratio = np.log2(np.clip(normalize_depth(depth), 0.01, None))

Circular Binary Segmentation (CBS)

Goal: split a chromosome's log2-ratio track into piecewise-constant segments at true CNV breakpoints. Approach: recursively split at the point that maximizes the absolute mean difference between the two resulting halves; stop when a segment is too small or the best split is too weak to be real signal (this is a simplified stand-in for the permutation-based test used by real CBS/DNAcopy).

def cbs_segment(log2_ratios, min_segment=10, alpha=0.01):
    """Simplified CBS: recursively split at the point of max absolute mean difference.

    Parameters
    ----------
    log2_ratios : np.ndarray
        Per-bin log2 ratios for one chromosome.
    min_segment : int
        Minimum bins per segment (stops recursion on tiny segments).
    alpha : float
        Unused placeholder for a real significance threshold (real CBS uses a
        permutation test here instead of a fixed delta cutoff).

    Returns
    -------
    list of (start_bin, end_bin, mean_log2) tuples, sorted by position.
    """
    segments = []

    def _split(lo, hi):
        if hi - lo < min_segment * 2:
            segments.append((lo, hi, log2_ratios[lo:hi].mean()))
            return
        best_delta, best_k = 0, None
        for k in range(lo + min_segment, hi - min_segment):
            delta = abs(log2_ratios[lo:k].mean() - log2_ratios[k:hi].mean())
            if delta > best_delta:
                best_delta, best_k = delta, k
        if best_k is None or best_delta < 0.3:
            segments.append((lo, hi, log2_ratios[lo:hi].mean()))
            return
        _split(lo, best_k)
        _split(best_k, hi)

    _split(0, len(log2_ratios))
    return sorted(segments)

CN State Calling

Goal: map each segment's mean log2 ratio to an integer copy-number state. Approach: walk from the highest CN state downward, returning the first state whose lower log2 boundary is met — boundaries are midpoints between the theoretical log2 values for adjacent CN states (log2(CN/2) for a diploid sample).

def call_cn_state(mean_log2, thresholds=None):
    """Map a segment's mean log2 ratio to an integer copy-number state (diploid baseline)."""
    if thresholds is None:
        # lower log2 boundary for each CN state; midpoints between log2(cn/2) values
        thresholds = {0: float("-inf"), 1: -1.5, 2: -0.4, 3: 0.35, 4: 0.80, 5: 1.25}
    for cn in sorted(thresholds.keys(), reverse=True):
        if mean_log2 >= thresholds[cn]:
            return cn
    return 0

segs = cbs_segment(log2_ratio)
calls = [{"start": s, "end": e, "log2": round(m, 3), "CN": call_cn_state(m)} for s, e, m in segs]

Purity Correction

Goal: back out the true tumor-cell copy number from an observed log2 ratio when tumor purity is known (or estimated). Approach: invert observed_CN = purity * tumor_CN + (1 - purity) * 2.

def purity_correct(observed_log2, purity, ploidy=2):
    """Recover the tumor-cell copy number from an observed log2 ratio and known purity.

    observed = purity * tumor_CN + (1 - purity) * ploidy
    => tumor_CN = (observed - (1 - purity) * ploidy) / purity
    """
    observed_cn = ploidy * (2 ** observed_log2)
    tumor_cn = (observed_cn - (1 - purity) * ploidy) / purity
    return max(0.0, tumor_cn)

Gene-Level Annotation

Goal: attach gene symbols to CNV segments so amplified oncogenes / deleted tumor suppressors can be flagged. Approach: a simple interval-overlap join between a segment table and a gene BED table (swap in pyranges/bedtools for large-scale annotation; see bio-genome-intervals-interval-arithmetic).

import pandas as pd

def annotate_segments(segments_df, genes_df):
    """Find genes overlapping each CNV segment (both DataFrames use chrom/start/end).

    segments_df: columns chrom, start, end, CN
    genes_df: columns chrom, start, end, gene
    """
    hits = []
    for chrom, sub_genes in genes_df.groupby("chrom"):
        sub_segs = segments_df[segments_df["chrom"] == chrom]
        for _, seg in sub_segs.iterrows():
            overlapping = sub_genes[(sub_genes["start"] < seg["end"]) & (sub_genes["end"] > seg["start"])]
            for gene in overlapping["gene"]:
                hits.append({"chrom": chrom, "start": seg["start"], "end": seg["end"],
                             "CN": seg["CN"], "gene": gene})
    return pd.DataFrame(hits)

Downstream: amplified oncogenes (COSMIC Cancer Gene Census), deleted tumor suppressors (TP53, BRCA1/2, RB1). Focal events (<3 Mb) are more often drivers than arm-level events. TCGA data: gdc.cancer.gov portal, TCGAbiolinks (R), or xenahubs.net; GISTIC2 identifies recurrent amplifications/deletions across a cohort.

Pitfalls

• Coordinate systems: BED uses 0-based half-open; VCF/GFF use 1-based inclusive — mixing them causes off-by-one errors
• CN=0 gives log2 = -inf: clip depth/ratios before log-transforming, or homozygous deletions will produce -inf/NaN that break segmentation and plotting
• Purity/ploidy confounding: without purity correction, low-purity samples show attenuated log2 signal that can be mistaken for a smaller CNV or missed entirely
• Batch effects: always check for batch confounding (sequencing run, capture kit) before interpreting biological signal
• Multiple testing: apply FDR correction (Benjamini-Hochberg) when testing thousands of bins/segments simultaneously

See Also

• bio-copy-number-cnvkit-analysis — production CNV calling from BAM files
• bio-copy-number-gatk-cnv — GATK's CNV calling workflow
• bio-copy-number-cnv-annotation — richer gene/pathway annotation of CNV calls
• bio-copy-number-cnv-visualization — genome-wide CNV plotting