bio-copy-number-cnv-annotation

What this skill does

## Version Compatibility

Reference examples tested with: bedtools 2.31+, AnnotSV 3.4+, Python 3.10+ with pybedtools 0.9+, pandas 2.2+, pysam 0.22+; R 4.3+ with clusterProfiler 4.10+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `bedtools --version`, `AnnotSV --version`
- Python: `pip show pybedtools pandas pysam`
- R: `packageVersion('clusterProfiler')`

If code throws an error, introspect the installed package and adapt the example. AnnotSV output column names change between major versions — verify against the installed version.

# CNV Annotation

**"Annotate my CNV calls with the genes they affect"** -> Overlap CNV segments with gene models, dosage-sensitivity maps, and clinical databases. The hard part is not the intersection — it is deciding *which* genes matter. A focal amplification overlapping 30 genes usually has one driver (the peak gene); a deletion's consequence depends on whether each gene is dosage-sensitive and whether the whole gene or only part is removed.

- CLI: `bedtools intersect -a cnvs.bed -b genes.bed -wa -wb`; `AnnotSV` for full annotation
- Python: `pybedtools` for interval logic; `pysam` for VCF database queries

## Annotation Strategy — Pick the Database for the Question

| Question | Resource | What it answers |
|----------|----------|-----------------|
| Which genes does this CNV span? | RefSeq/GENCODE gene BED | Raw overlap (not yet consequence) |
| Is loss of this gene damaging? | ClinGen haploinsufficiency (HI) score | Dosage sensitivity to deletion |
| Is gain of this gene damaging? | ClinGen triplosensitivity (TS) score | Dosage sensitivity to duplication |
| Is this CNV common in the population? | gnomAD-SV, DGV, 1000G CNV | Benign-frequency filtering |
| Is this a known recurrent disorder locus? | ClinGen dosage regions, DECIPHER | Genomic-disorder context |
| Is this a cancer driver? | COSMIC Cancer Gene Census, OncoKB | Oncogene vs tumor-suppressor role |
| Is there pathogenic small-variant content? | ClinVar | Coincident SNV/indel pathogenicity |
| One-shot comprehensive annotation + ranking | AnnotSV | Aggregates most of the above |

For constitutional CNV *classification* (assigning pathogenic/VUS/benign), the annotated output feeds the ACMG/ClinGen points framework — see germline-cnv-interpretation. For cohort-level recurrence and driver-peak identification, see recurrent-cnv.

## The Core Distinction: Overlap Is Not Consequence

A CNV overlapping a gene does not necessarily change that gene's dosage in a way that matters. Three refinements separate annotation from interpretation:

1. **Whole-gene vs partial overlap.** A deletion spanning an entire gene removes one copy (clean haploinsufficiency test). A deletion removing only the last two exons creates a truncated allele — a different, often more damaging, consequence. Always record the fraction of each gene covered and whether coding exons or only introns/UTRs are hit.
2. **Dosage sensitivity.** Most genes tolerate single-copy loss. ClinGen HI/TS scores (3 = sufficient evidence for dosage sensitivity, 0 = no evidence, 30 = gene associated with an autosomal-recessive phenotype, 40 = dosage sensitivity unlikely) indicate which genes' loss/gain is actually consequential.
3. **Driver vs passenger in focal events.** A focal amplification carries many genes; the driver is the one under selection, typically at the recurrence peak across a cohort (GISTIC) and a known oncogene. Annotating all 30 genes as "amplified" overstates.

## Gene Overlap with bedtools

**Goal:** Find genes overlapping each CNV segment, recording overlap extent.

**Approach:** Convert segments to BED, intersect with a gene model, keep both feature sets (`-wo` reports the overlap length) so partial vs whole-gene overlap is recoverable.

```bash
# Segments to BED (CNVkit .cns example; columns chrom/start/end/log2)
awk 'NR>1 {print $1"\t"$2"\t"$3"\t"$5}' sample.cns > sample.cnv.bed

# Intersect; -wo appends the number of overlapping bases
bedtools intersect -a sample.cnv.bed -b gencode.genes.bed -wo > cnv_gene_overlap.txt
```

## Comprehensive Annotation with AnnotSV

**Goal:** Annotate CNVs against genes, dosage maps, population frequency, and clinical databases in one pass, with a built-in pathogenicity ranking.

**Approach:** Export CNVs to VCF or BED and run AnnotSV; it returns a "full" line per gene plus a "split" summary, with an ACMG-aligned rank (1-5).

```bash
AnnotSV \
    -SVinputFile sample.cnv.vcf \
    -genomeBuild GRCh38 \
    -annotationMode both \
    -outputFile sample_annotated.tsv

# Output includes: overlapped genes, ClinGen HI/TS, gnomAD-SV/DGV frequency, OMIM,
# ClinVar, DECIPHER, and an ACMG-class rank per SV.
```

AnnotSV's rank is a useful triage signal, not a final classification — confirm against the ClinGen points framework for clinical reporting.

## Dosage-Sensitivity and Driver Annotation

**Goal:** Tag each affected gene with its dosage sensitivity and, for tumors, its driver role, so passengers can be separated from drivers.

**Approach:** Join the gene-overlap table to the ClinGen dosage map (HI/TS scores) and to the COSMIC Cancer Gene Census; flag CNVs whose direction matches a known mechanism (oncogene amplified, tumor suppressor deleted).

```python
import pandas as pd

def annotate_dosage_and_drivers(overlap_tsv, clingen_dosage, cgc_file):
    '''Tag overlapped genes with ClinGen HI/TS and COSMIC driver role.'''
    cols = ['cnv_chrom', 'cnv_start', 'cnv_end', 'log2',
            'gene_chrom', 'gene_start', 'gene_end', 'gene', 'overlap_bp']
    df = pd.read_csv(overlap_tsv, sep='\t', names=cols)
    df['gene_len'] = df['gene_end'] - df['gene_start']
    df['gene_frac_covered'] = (df['overlap_bp'] / df['gene_len']).clip(upper=1.0)
    df['whole_gene'] = df['gene_frac_covered'] >= 0.99

    dosage = pd.read_csv(clingen_dosage, sep='\t')  # gene, HI_score, TS_score
    df = df.merge(dosage, on='gene', how='left')

    cgc = pd.read_csv(cgc_file, sep='\t')
    role = dict(zip(cgc['Gene Symbol'], cgc['Role in Cancer']))
    df['driver_role'] = df['gene'].map(role)

    # Direction-consistent driver hits: oncogene amplified or TSG deleted.
    df['driver_hit'] = (
        ((df['log2'] > 0.3) & df['driver_role'].fillna('').str.contains('oncogene')) |
        ((df['log2'] < -0.3) & df['driver_role'].fillna('').str.contains('TSG')))
    return df
```

## Population-Frequency Filtering

**Goal:** Remove common, presumed-benign CNVs before clinical interpretation.

**Approach:** Reciprocal-overlap match each CNV against a population SV catalog (gnomAD-SV, DGV); a CNV with high reciprocal overlap to a common population CNV of the same type is likely benign.

```bash
# 50% reciprocal overlap (-f 0.5 -r): same-type, similar-extent population match.
# Reciprocal overlap, not one-sided, prevents a tiny CNV inside a huge population CNV
# (or vice versa) from being wrongly matched.
bedtools intersect -a sample.cnv.bed -b gnomad_sv.bed -f 0.5 -r -wa -wb \
    > cnv_population_match.txt
```

## Pathway Enrichment of Affected Genes

**Goal:** Test whether genes in amplified (or deleted) regions are enriched for pathways.

**Approach:** Extract genes by CNV direction, map to Entrez IDs, run GO/KEGG enrichment. Caveat: CNVs are large and gene-dense, so enrichment is biased toward whatever pathways cluster in CNV-prone genomic regions — interpret as hypothesis-generating.

```r
library(clusterProfiler)
library(org.Hs.eg.db)

amp_genes <- unique(cnv_annot$gene[cnv_annot$log2 > 0.3])
entrez <- na.omit(mapIds(org.Hs.eg.db, keys = amp_genes,
                         keytype = 'SYMBOL', column = 'ENTREZID'))
go_bp <- enrichGO(gene = entrez, OrgDb = org.Hs.eg.db, ont = 'BP',
                  pAdjustMethod = 'BH', qvalueCutoff = 0.05)
```

## Failure Modes

### Genome-build mismatch between CNVs and annotation

**Trigger:** CNV coordinates on GRCh37 intersected with a GRCh38 gene model (or vice versa).

**Mechanism:** Coordinates silently s