bio-copy-number-gatk-cnv
$
npx mdskill add GPTomics/bioSkills/bio-copy-number-gatk-cnvCall copy number variants using GATK best practices
- Detects somatic and germline CNVs from WGS or WES data.
- Depends on GATK CLI tools for read counting, denoising, and modeling.
- Decides workflow steps based on input type like tumor-normal pairs.
- Delivers segmented copy ratio states via standard output files.
SKILL.md
.github/skills/bio-copy-number-gatk-cnvView on GitHub ↗
---
name: bio-copy-number-gatk-cnv
description: Call copy number variants using GATK best practices workflow. Supports both somatic (tumor-normal) and germline CNV detection from WGS or WES data. Use when following GATK best practices or integrating CNV calling with other GATK variant pipelines.
tool_type: cli
primary_tool: gatk
---
## Version Compatibility
Reference examples tested with: GATK 4.5+
Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
# GATK CNV Workflow
**"Call CNVs using GATK best practices"** → Collect read counts, build a panel of normals, denoise tumor coverage, model segments with allelic counts, and call copy ratio states.
- CLI: `gatk CollectReadCounts` → `gatk DenoiseReadCounts` → `gatk ModelSegments` → `gatk CallCopyRatioSegments`
## Somatic CNV Workflow Overview
```
1. PreprocessIntervals → intervals.interval_list
2. CollectReadCounts → sample.counts.hdf5
3. CreateReadCountPanelOfNormals → pon.hdf5
4. DenoiseReadCounts → sample.denoised.tsv
5. CollectAllelicCounts → sample.allelicCounts.tsv
6. ModelSegments → sample.modelFinal.seg
7. CallCopyRatioSegments → sample.called.seg
```
## Step 1: Preprocess Intervals
**Goal:** Prepare genomic intervals for read counting, handling both WES and WGS modes.
**Approach:** Use PreprocessIntervals to bin or merge target intervals with appropriate padding.
```bash
# For WES/targeted
gatk PreprocessIntervals \
-R reference.fa \
-L targets.interval_list \
--bin-length 0 \
--interval-merging-rule OVERLAPPING_ONLY \
-O preprocessed.interval_list
# For WGS
gatk PreprocessIntervals \
-R reference.fa \
--bin-length 1000 \
--padding 0 \
-O wgs.interval_list
```
## Step 2: Collect Read Counts
**Goal:** Count reads per interval for each sample.
**Approach:** Run CollectReadCounts on each BAM against the preprocessed intervals.
```bash
# For each sample
gatk CollectReadCounts \
-R reference.fa \
-I sample.bam \
-L preprocessed.interval_list \
--interval-merging-rule OVERLAPPING_ONLY \
-O sample.counts.hdf5
```
## Step 3: Create Panel of Normals
**Goal:** Build a reference panel from multiple normal samples for denoising.
**Approach:** Combine normal sample count HDF5 files into a single panel-of-normals using PCA-based denoising.
```bash
# Combine multiple normal samples
gatk CreateReadCountPanelOfNormals \
-I normal1.counts.hdf5 \
-I normal2.counts.hdf5 \
-I normal3.counts.hdf5 \
--minimum-interval-median-percentile 5.0 \
-O cnv_pon.hdf5
```
## Step 4: Denoise Read Counts
**Goal:** Remove systematic noise from tumor read counts using the panel of normals.
**Approach:** Apply DenoiseReadCounts with the PoN to produce standardized and denoised copy ratio profiles.
```bash
# Using panel of normals
gatk DenoiseReadCounts \
-I tumor.counts.hdf5 \
--count-panel-of-normals cnv_pon.hdf5 \
--standardized-copy-ratios tumor.standardized.tsv \
--denoised-copy-ratios tumor.denoised.tsv
```
## Step 5: Collect Allelic Counts
**Goal:** Capture allele-specific information at known heterozygous SNP sites for LOH detection.
**Approach:** Run CollectAllelicCounts against common SNP sites to generate allelic count profiles.
```bash
# From known SNP sites (for LOH detection)
gatk CollectAllelicCounts \
-R reference.fa \
-I tumor.bam \
-L common_snps.vcf \
-O tumor.allelicCounts.tsv
```
## Step 6: Model Segments
**Goal:** Jointly segment copy ratio and allelic data to identify CNV regions.
**Approach:** Run ModelSegments with denoised ratios and allelic counts from both tumor and matched normal.
```bash
# Somatic with matched normal allelic counts
gatk ModelSegments \
--denoised-copy-ratios tumor.denoised.tsv \
--allelic-counts tumor.allelicCounts.tsv \
--normal-allelic-counts normal.allelicCounts.tsv \
--output-prefix tumor \
-O results/
# Output files: tumor.cr.seg, tumor.modelFinal.seg, tumor.hets.tsv
```
## Step 7: Call Copy Ratio Segments
**Goal:** Assign amplification, deletion, or neutral calls to each segment.
**Approach:** Apply CallCopyRatioSegments to convert continuous log2 ratios into discrete CN states.
```bash
gatk CallCopyRatioSegments \
-I results/tumor.cr.seg \
-O results/tumor.called.seg
```
## Plotting
**Goal:** Visualize denoised copy ratios and modeled segments with allelic information.
**Approach:** Use GATK PlotDenoisedCopyRatios and PlotModeledSegments to generate standardized plots.
```bash
# Plot copy ratios and segments
gatk PlotDenoisedCopyRatios \
--standardized-copy-ratios tumor.standardized.tsv \
--denoised-copy-ratios tumor.denoised.tsv \
--sequence-dictionary reference.dict \
--minimum-contig-length 46709983 \
--output-prefix tumor \
-O plots/
# Plot segments with allelic information
gatk PlotModeledSegments \
--denoised-copy-ratios tumor.denoised.tsv \
--allelic-counts results/tumor.hets.tsv \
--segments results/tumor.modelFinal.seg \
--sequence-dictionary reference.dict \
--minimum-contig-length 46709983 \
--output-prefix tumor \
-O plots/
```
## Germline CNV Workflow
```bash
# For germline: use cohort mode
# 1. Collect counts (same as above)
# 2. Determine contig ploidy
gatk DetermineGermlineContigPloidy \
-I sample1.counts.hdf5 \
-I sample2.counts.hdf5 \
--model cohort_ploidy_model \
--contig-ploidy-priors ploidy_priors.tsv \
-O ploidy-calls/
# 3. Call germline CNVs
gatk GermlineCNVCaller \
--run-mode COHORT \
-I sample1.counts.hdf5 \
-I sample2.counts.hdf5 \
--contig-ploidy-calls ploidy-calls/ploidy_calls \
--annotated-intervals annotated_intervals.tsv \
--output-prefix cohort \
-O germline_cnv_calls/
# 4. Post-process calls per sample
gatk PostprocessGermlineCNVCalls \
--calls-shard-path germline_cnv_calls/cohort-calls \
--model-shard-path germline_cnv_calls/cohort-model \
--sample-index 0 \
--contig-ploidy-calls ploidy-calls/ploidy_calls \
--sequence-dictionary reference.dict \
--output-genotyped-intervals sample1.genotyped.tsv \
--output-denoised-copy-ratios sample1.denoised.tsv \
-O sample1_segments.vcf
```
## Complete Somatic Pipeline Script
```bash
#!/bin/bash
REFERENCE=reference.fa
INTERVALS=targets.interval_list
PON=cnv_pon.hdf5
SNP_SITES=common_snps.vcf
TUMOR=$1
NORMAL=$2
OUTDIR=$3
mkdir -p $OUTDIR
# Collect read counts
gatk CollectReadCounts -R $REFERENCE -I $TUMOR -L $INTERVALS \
-O $OUTDIR/tumor.counts.hdf5
gatk CollectReadCounts -R $REFERENCE -I $NORMAL -L $INTERVALS \
-O $OUTDIR/normal.counts.hdf5
# Denoise
gatk DenoiseReadCounts -I $OUTDIR/tumor.counts.hdf5 \
--count-panel-of-normals $PON \
--standardized-copy-ratios $OUTDIR/tumor.standardized.tsv \
--denoised-copy-ratios $OUTDIR/tumor.denoised.tsv
# Allelic counts
gatk CollectAllelicCounts -R $REFERENCE -I $TUMOR -L $SNP_SITES \
-O $OUTDIR/tumor.allelicCounts.tsv
gatk CollectAllelicCounts -R $REFERENCE -I $NORMAL -L $SNP_SITES \
-O $OUTDIR/normal.allelicCounts.tsv
# Model and call
gatk ModelSegments \
--denoised-copy-ratios $OUTDIR/tumor.denoised.tsv \
--allelic-counts $OUTDIR/tumor.allelicCounts.tsv \
--normal-allelic-counts $OUTDIR/normal.allelicCounts.tsv \
--output-prefix tumor -O $OUTDIR/
gatk CallCopyRatioSegments -I $OUTDIR/tumor.cr.seg -O $OUTDIR/tumor.called.seg
```
## Key Output Files
| File | Description |
|------|-------------|
| .counts.hdf5 | Raw read counts per interval |
| .denoised.tsv | Denoised log2 copy ratios |
| .modelFinal.seg | Segmented copy ratios with confidence |
| .called.seg | Final called segments with CN state |
| .hets.tsv | Heterozygous SNP allelic counts |
## Related Skills
- copy-number/cnvkit-analysis - Alternative CNV caller
- copy-number/cnv-visualization - Plotting results
- alignment-files/bam-statistics - Input BAM QC
- variant-calling/variant-calling - SNP calling for allelic counts
More from GPTomics/bioSkills
- bio-admet-predictionPredicts ADMET properties using ADMETlab 3.0 API or DeepChem models. Estimates bioavailability, CYP inhibition, hERG liability, and 119 toxicity endpoints with uncertainty quantification. Filters for PAINS and other structural alerts. Use when filtering compounds for drug-likeness or prioritizing leads by predicted safety.
- bio-alignment-amplicon-clippingTrim PCR primers from aligned reads in amplicon-panel BAMs using samtools ampliconclip. Use when processing SARS-CoV-2 ARTIC, hereditary cancer panels, ctDNA hot-spot panels, or any amplicon assay where primer-derived bases would falsely confirm reference at primer footprints.
- bio-alignment-filteringFilter alignments by flags, mapping quality, and regions using samtools view and pysam. Use when extracting specific reads, removing low-quality alignments, or subsetting to target regions.
- bio-alignment-indexingCreate and use BAI/CSI indices for BAM/CRAM files using samtools and pysam. Use when enabling random access to alignment files or fetching specific genomic regions.
- bio-alignment-ioRead, write, and convert multiple sequence alignment files using Biopython Bio.AlignIO. Supports Clustal, PHYLIP, Stockholm, FASTA, Nexus, and other alignment formats for phylogenetics and conservation analysis. Use when reading, writing, or converting alignment file formats.
- bio-alignment-msa-parsingParse and analyze multiple sequence alignments using Biopython. Extract sequences, identify conserved regions, analyze gaps, work with annotations, and manipulate alignment data for downstream analysis. Use when parsing or manipulating multiple sequence alignments.
- bio-alignment-msa-statisticsCalculate alignment statistics including sequence identity, conservation scores, substitution matrices, and similarity metrics. Use when comparing alignment quality, measuring sequence divergence, and analyzing evolutionary patterns.
- bio-alignment-multiplePerform multiple sequence alignment using MAFFT, MUSCLE5, ClustalOmega, or T-Coffee. Guides tool and algorithm selection based on dataset size, sequence divergence, and downstream application. Use when aligning three or more homologous sequences for phylogenetics, conservation analysis, or evolutionary studies.
- bio-alignment-pairwisePerform pairwise sequence alignment using Biopython Bio.Align.PairwiseAligner. Use when comparing two sequences, finding optimal alignments, scoring similarity, and identifying local or global matches between DNA, RNA, or protein sequences.
- bio-alignment-sortingSort alignment files by coordinate or read name using samtools and pysam. Use when preparing BAM files for indexing, variant calling, or paired-end analysis.