Mappability annotation - flag variants where short reads mismap

🤖 Written by Claude

Follow-up to #693 (retire GeneInfo). The `Pseudogene` category we deleted there was a proxy for what curators actually want to know: **short reads mismap here, treat this call with suspicion**. Mappability measures that directly, and it is positional rather than per-gene, so it belongs on `VariantAnnotation` alongside the other VEP `--custom` tracks.

## 1. Source: UCSC Umap/Bismap mappability

UCSC hosts the Umap tracks (Karimzadeh et al., *Umap and Bismap: quantifying genome and methylome mappability*, NAR 2018, https://academic.oup.com/nar/article/46/20/e120/5086676) at:

**https://hgdownload.soe.ucsc.edu/gbdb/hg38/hoffmanMappability/**

Available for k-mer lengths 24, 36, 50 and 100, in two forms:

* `k100.Umap.MultiTrackMappability.bw` — bigWig, continuous mappability score
* `k100.Unique.Mappability.bb` — bigBed, unique-mappability regions (`bigBedToBed` for plain BED)

Pick the k-mer closest to our read lengths — k100 or k150 for typical Illumina.

Advantage over a curated gene list: it is computed genome-wide, so it cannot have the "missing the hardest genes" problem that sank the GIAB CMRG idea (see §6). PMS2, CYP21A2, HBA1/2, SMN1/2 and friends get low mappability because they genuinely are low mappability, not because someone remembered to add them.

**Two things to confirm before starting:**

1. **GRCh37.** `/gbdb/hg19/` has no `hoffmanMappability` directory — only `reMap/` and `transMap/`. Umap does publish hg19 data (via https://bismap.hoffmanlab.org and Zenodo), so the GRCh37 file likely has to come from there rather than UCSC. Needs checking, since we annotate both builds.
2. **Licence.** Not yet checked. Umap is academic (Hoffman lab) rather than NIST, so it will not have the automatic public-domain status GIAB data has. Confirm redistribution terms before vendoring anything — though note that unlike a small gene list, these are large binary tracks that likely belong in the annotation data directory rather than the repo anyway (see §3).

## 2. Where it goes: `VariantAnnotation` via VEP `--custom`

The codebase already has a first-class framework for exactly this, so this is a well-trodden path rather than new machinery:

* `VEPCustom` enum in `annotation/models/models_enums.py:130` — already carries BED (`REPEAT_MASKER`) and bigWig (`PHASTCONS_*`, `PHYLOP_*`) entries.
* `_get_custom_params_list` in `annotation/vep_annotation.py:37` — already handles both `bed` and `bw` extensions and builds the `--custom` argument.
* `ColumnVEPField` rows drive which source fields map to which `VariantAnnotation` columns.

So the work is: add a `VEPCustom` entry, add the file to the VEP config/annotation data dir, add a `ColumnVEPField`, add the field to `VariantAnnotation` (`annotation/models/models.py:1526`), and register a `VariantGridColumn`. RepeatMasker is the closest precedent for the BED route; phastCons/phyloP for the bigWig route.

Being on `VariantAnnotation` means the grid column, node filtering and VCF export all come essentially free — no separate decision needed, unlike the gene-level design this issue previously proposed.

## 3. bigWig vs BED — decide first

**bigWig** gives a continuous score. `_get_custom_params_list` already supports `summary_stats` (`vep_annotation.py:80-83`), so `summary_stats=min` over the variant span yields "worst mappability this variant touches", which is the semantically right thing for an indel or SV spanning a boundary. A float column lets curators see 0.0 vs 0.5 vs 1.0 and pick their own threshold at filter time. Precedent: the conservation tracks.

**BED/bigBed** gives a boolean overlap — simpler, smaller, no threshold argument, but throws away the gradient. Precedent: RepeatMasker.

I would lean bigWig for the richer signal, but see the SV cost below, which cuts the other way.

### The SV landmine

Worth knowing before choosing bigWig: `annotation/tasks/annotate_variants.py:533` and `annotation/sv_conservation.py` record that the four conservation `--custom` bigWig overlaps have an *"O(SV-span) cost [that] makes large SVs never finish"*. That was fixed in #1657 by dropping those `--custom` args for SVs and computing the values with pyBigWig instead.

A fifth bigWig track would hit the same wall and need the same treatment — `annotation/sv_conservation.py` is set up for it, but it is real extra work rather than a free ride. RepeatMasker shows BED overlaps are not entirely free either: it needs `num_records=10000` to stop VEP truncating with `"..."` (`vep_annotation.py:52-53`).

Given the file sizes (bigWig tracks are large binaries), these belong in the annotation data directory alongside the other VEP custom files, not vendored in the repo — so the "vendor a small file plus a data migration" plan from the earlier version of this issue no longer applies.

## 4. Gene page display

#693 asked for a badge at the top of the gene page, and that is still worth having — but with mappability it becomes a *rollup* of variant-level data rather than the primary artefact.

Open question: what the rollup is. Something like "fraction of coding bases below a mappability threshold" is defensible and gives a degree rather than a boolean, but it needs exon coordinates per build and a threshold we can justify. Simplest first cut is to ship §2 (the variant-level annotation, where the clinical value is) and treat the gene page badge as a follow-up once we can see the real distribution of values.

If we do add it: badge next to `

Gene Symbol : {{ gene_symbol }}

` (`genes/templates/genes/view_gene_symbol.html:153`), wired as a `cached_property` on `GeneSymbolViewInfo` plus an entry in the `LazyAttribute.lazy_context` list, structurally like `gene_constraint` (`genes/views/views.py:185-186`).

Note mappability is genuinely **build-specific** — GRCh37's false duplications mean a gene can be hard on one build and fine on the other. Any gene-level rollup therefore needs a genome build dimension, which is a change from the build-independent gene flag earlier drafts of this issue assumed.

## 5. Superseded: gene-list approaches

Earlier drafts proposed a vendored gene list plus a `GeneSymbol`-keyed table and a data migration. Mappability on `VariantAnnotation` is better on every axis — complete rather than curated, positional rather than per-gene, and it reuses machinery that already exists. Recording the rejected alternatives so they are not revisited:

* **A `GeneList`** — user-space objects with permissions, categories and GeneGrid columns. Using one for curated reference data is what made GeneInfo unusable.
* **`GeneVersion`** — `unique_together = ("gene", "version", "genome_build")` duplicates the value per build × consortium × version, and `import_source = FK(GeneAnnotationImport, on_delete=CASCADE)` means rows cascade away on gene annotation reimport.
* **A field on `GeneSymbol`** — pure identity table; every future reference dataset would accrete a column.
* **`GeneAnnotation`** — `variantgrid/settings/components/annotation_settings.py:7` states "GeneAnnotation is only in analyses"; it is keyed on `Gene` per `GeneAnnotationVersion` and the gene page never reads it.

## 6. Evaluated and rejected: GIAB CMRG

#693 proposed GIAB's Challenging Medically Relevant Genes and claimed "roughly 395 genes" picking up "PMS2, SMN1/2, CYP21A2, HBA1/2, STRC, NEB etc.". **Both claims are wrong**, recorded here so nobody re-treads it.

The actual file ([GRCh38_mrg_bench_gene.bed](https://raw.githubusercontent.com/usnistgov/cmrg-benchmarkset-manuscript/master/data/gene_coords/unsorted/GRCh38_mrg_bench_gene.bed), 8 KB) has **273 genes**, and of those named only **SMN1** is present — PMS2, SMN2, CYP21A2, HBA1, HBA2, STRC, NEB and CFTR are all absent.

Per `README_GIAB_medical_gene_benchmark.md` in the [FTP release](https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/CMRG_v1.00/), the set is medically relevant genes "covered <90% by v4.2.1 in HG002" but then cut down to what the hifiasm v0.11 assembly could resolve, explicitly excluding genes "not fully and accurately resolved by the hifiasm v0.11 assembly (e.g. SMN2)" and those with "very large structural or copy number variation ... (e.g. LPA, CR1, and KIR genes)".

So it is "hard genes GIAB managed to curate a benchmark for", not a difficulty list — several of the hardest are missing precisely because they are hardest. Licensing was fine (NIST, not subject to US copyright) if it is ever wanted as a cross-check.