# Bio Applied Copy Number Analysis

> Call CNVs from binned read-depth: GC-bias normalization, circular binary segmentation (CBS), log2-ratio-to-CN-state calling, tumor-purity correction. Use for CNV/copy-number segmentation, log2 ratio analysis, or gain/deletion/amplification calls.

- Skill: `pavel-kravchenko/bio-applied-copy-number-analysis` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-copy-number-analysis`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-copy-number-analysis/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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.

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

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

```python
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`.

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

```python
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