bio-copy-number-gatk-cnv

What this skill does

## Version Compatibility

Reference examples tested with: GATK 4.5+ (gatk4), Python 3.10+ (gcnv conda env), R 4.3+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `gatk --version` then `gatk <ToolName> --help` to confirm arguments
- gCNV requires a working `gatkcondaenv` (theano/tensorflow stack) — `gatk` will report if the Python environment is missing

GATK 4.5+ gCNV inference defaults are tuned for whole-exome data; whole-genome runs generally need parameter changes. If a tool reports an unrecognized argument, check the help for that exact GATK version rather than retrying.

# GATK CNV Workflows

**"Call CNVs the GATK way"** -> GATK has two *separate* CNV workflows that share almost no tools. Picking the wrong one is the most common mistake.

- Somatic CNV: `CollectReadCounts` -> `DenoiseReadCounts` -> `ModelSegments` -> `CallCopyRatioSegments`. Tumor copy-ratio segments, optionally allele-aware.
- Germline gCNV: `DetermineGermlineContigPloidy` -> `GermlineCNVCaller` -> `PostprocessGermlineCNVCalls`. Per-sample germline CN genotypes (VCF).

## Critical: What GATK Somatic CNV Does NOT Provide

`ModelSegments` + `CallCopyRatioSegments` produce **copy-ratio segments** and a **minor-allele fraction** per segment, and the "call" is a simple t-test emitting `+` / `-` / `0`. This is **not** integer allele-specific copy number, **not** tumor purity, and **not** ploidy. Practitioners routinely assume parity with ASCAT/FACETS and there is none. For integer allele-specific CN, purity, ploidy, LOH state, or whole-genome-doubling status, use allele-specific-copy-number (ASCAT, Sequenza, FACETS, or PureCN — PureCN can even reuse the GATK `ModelSegments` segmentation as input).

## Somatic vs Germline — Choosing the Workflow

| Question | Somatic CNV | Germline gCNV |
|----------|-------------|---------------|
| Input | One tumor (+ optional matched normal) | A cohort of constitutional samples |
| Output | Copy-ratio segments, +/-/0 call, minor-allele fraction | Integer germline CN genotype VCF per sample |
| Normalization | Tangent (projection onto PoN subspace) | PCA batching + Bayesian read-depth model |
| Cohort needed | PoN of normals for denoising | >= ~100 technically matched samples (cohort mode) |
| Use for | Tumor SCNAs, focal amplifications/deletions | Rare/de novo germline CNVs, NDD/Mendelian cohorts |

## Decision Tree by Scenario

| Scenario | Workflow | Key parameters |
|----------|----------|----------------|
| Tumor-normal WGS/WES, want SCNAs | Somatic, with matched-normal allelic counts | `PreprocessIntervals --bin-length 1000` (WGS) or `0` (WES) |
| Tumor-only somatic CNV | Somatic, no matched-normal allelic counts | Genotype hets in the case sample; expect more no-calls |
| Rare germline CNV, exome cohort >= 100 | gCNV cohort mode | Run `DetermineGermlineContigPloidy` cohort first |
| New sample vs an existing gCNV model | gCNV case mode | Must reuse identical scatter count and interval list |
| Need integer ASCN / purity / ploidy | Neither — escalate | Use allele-specific-copy-number |
| Targeted panel (< few hundred genes) | Prefer CNVkit | GATK interval models are unstable on tiny panels |

## Somatic CNV Pipeline

```bash
# 1. Preprocess and annotate intervals (WES: bin-length 0 = use exome targets as-is)
gatk PreprocessIntervals -R ref.fa -L targets.interval_list \
    --bin-length 0 --interval-merging-rule OVERLAPPING_ONLY -O preprocessed.interval_list
gatk AnnotateIntervals -R ref.fa -L preprocessed.interval_list \
    --interval-merging-rule OVERLAPPING_ONLY -O annotated.tsv     # GC content for FilterIntervals

# 2. Collect read counts (each BAM)
gatk CollectReadCounts -R ref.fa -I sample.bam -L preprocessed.interval_list \
    --interval-merging-rule OVERLAPPING_ONLY -O sample.counts.hdf5

# 3. Build the panel of normals (tangent-normalization basis)
# --minimum-interval-median-percentile 5.0 is the GATK CNV tutorial value (tool default 10.0)
gatk CreateReadCountPanelOfNormals \
    -I normal1.counts.hdf5 -I normal2.counts.hdf5 -I normalN.counts.hdf5 \
    --annotated-intervals annotated.tsv \
    --minimum-interval-median-percentile 5.0 -O cnv.pon.hdf5

# 4. Denoise tumor against the PoN (tangent normalization)
gatk DenoiseReadCounts -I tumor.counts.hdf5 --count-panel-of-normals cnv.pon.hdf5 \
    --standardized-copy-ratios tumor.standardizedCR.tsv \
    --denoised-copy-ratios tumor.denoisedCR.tsv

# 5. Allelic counts at common biallelic SNPs (tumor and matched normal)
gatk CollectAllelicCounts -R ref.fa -I tumor.bam -L common_snps.interval_list \
    -O tumor.allelicCounts.tsv
gatk CollectAllelicCounts -R ref.fa -I normal.bam -L common_snps.interval_list \
    -O normal.allelicCounts.tsv

# 6. Joint segmentation of copy ratio and allele fraction
gatk ModelSegments --denoised-copy-ratios tumor.denoisedCR.tsv \
    --allelic-counts tumor.allelicCounts.tsv \
    --normal-allelic-counts normal.allelicCounts.tsv \
    --output-prefix tumor -O segments/

# 7. Call each segment +/-/0 (simple t-test against the copy-ratio baseline)
gatk CallCopyRatioSegments -I segments/tumor.cr.seg -O segments/tumor.called.seg
```

`AnnotateIntervals` (step 1) and supplying `--annotated-intervals` to the PoN are frequently skipped — they enable explicit GC-bias correction and are recommended.

## Germline gCNV Pipeline

```bash
# 1. Determine contig ploidy across the cohort (karyotype + global depth)
gatk DetermineGermlineContigPloidy -L preprocessed.interval_list \
    --interval-merging-rule OVERLAPPING_ONLY \
    -I sample1.counts.hdf5 -I sampleN.counts.hdf5 \
    --contig-ploidy-priors ploidy_priors.tsv --output-prefix cohort -O ploidy-calls/

# 2. FilterIntervals — remove low-mappability / extreme-GC / low-count intervals
gatk FilterIntervals -L preprocessed.interval_list --annotated-intervals annotated.tsv \
    -I sample1.counts.hdf5 -I sampleN.counts.hdf5 \
    --interval-merging-rule OVERLAPPING_ONLY -O filtered.interval_list

# 3. GermlineCNVCaller, cohort mode (builds the model AND calls the cohort)
gatk GermlineCNVCaller --run-mode COHORT -L filtered.interval_list \
    --interval-merging-rule OVERLAPPING_ONLY \
    --contig-ploidy-calls ploidy-calls/cohort-calls \
    -I sample1.counts.hdf5 -I sampleN.counts.hdf5 \
    --output-prefix cohort -O gcnv-calls/

# 4. Post-process per sample into a genotyped VCF
gatk PostprocessGermlineCNVCalls \
    --calls-shard-path gcnv-calls/cohort-calls \
    --model-shard-path gcnv-calls/cohort-model \
    --contig-ploidy-calls ploidy-calls/cohort-calls \
    --sample-index 0 \
    --output-genotyped-intervals sample0.intervals.vcf.gz \
    --output-genotyped-segments sample0.segments.vcf.gz \
    --output-denoised-copy-ratios sample0.denoisedCR.tsv
```

Case mode (`--run-mode CASE`) scores a new sample against the cohort `*-model` shards; it must use the **identical** `filtered.interval_list` and the **same scatter count** as the cohort run, or it fails or produces incomparable calls.

## Failure Modes

### Tangent normalization removes a real CNV

**Trigger:** PoN is small (< ~20 normals) or contains samples that share a recurrent CNV (e.g. a common germline CNV, or a PoN accidentally built from tumors).

**Mechanism:** `DenoiseReadCounts` projects the tumor coverage profile onto the subspace spanned by the PoN's principal components. Any copy-number pattern present in that subspace is treated as "systematic noise" and subtracted. A CNV shared by PoN members is therefore normalized out of the tumor.

**Symptom:** A known event (recurrent amplification/deletion, or a common germline CNV) is absent from `denoisedCR.tsv`; denoised profile is suspiciously flat at that locus.

**Fix:** Build the PoN from >= 20-40 unrelated, tumor-free, process-matched normals. Never put tumors in the PoN. Cross-check against the `standardizedCR.tsv` (pre-tangent) profile — if the event is there but gone after denoising, the PoN ate it.

### Mistaking ModelSegments output for allele-