Run SpliceAI on demand for individual variants the precalculated scores miss

🤖 Written by Claude

Split out of #720, which has become the general "annotation pipeline types" issue. The framework part of that landed as the AnnotSV split (`ANNOTSV` pipeline type + `annotation/pipelines/` registry); this is the original ask that started #720 — running SpliceAI ourselves to cover what the precalculated scores miss.

## The constraint that drives the whole design

**SpliceAI runs are sparse, manual, and repeated over time.** A run is a handful of specially chosen variants, not a range. One variant in a range lock gets run today; a second variant in the *same* range lock gets run months later, possibly against a different SpliceAI version.

That breaks the provenance model everything else in the annotation pipeline uses. For VEP and AnnotSV, a FINISHED `AnnotationRun` over a range lock means *every* variant in that range has been handled — which is what range locks are for, and why #720 explicitly decided against per-variant provenance. SpliceAI cannot work that way: its coverage of a range lock is partial and grows in unpredictable steps.

Two things follow, before any pipeline design:

1. **Coverage is recorded per variant**, not per range lock.
2. **`AnnotationRun` is the wrong container for a manual run.** It is coupled to a range lock throughout — `variant_annotation_version`, `genome_build` and `annotation_consortium` all read through `annotation_range_lock`, `_trigger_dispatch` takes `annotation_run.annotation_range_lock.version_id`, and #1654 deliberately arranged that no rangeless `AnnotationRun` is ever committed. A manual run of three variants should not go through the lease/dispatch/reclaim machinery built for multi-hour VEP batches.

## 1. Where the record goes: a run collection against `Variant`

**Against `Variant`, via a collection that carries the SpliceAI details.**

Not `Allele`: SpliceAI scores a variant against a specific reference sequence, so a GRCh37 score can't be assumed valid for the GRCh38 representation of the same allele — liftover isn't guaranteed right.

Not `VariantAnnotation` either. The manual-run fact is version-independent — "we spent ~10s of CPU on this variant in this build" stays true across annotation rebuilds — and if it lived on the annotation row, an archived prior partition (`DataArchiveMixin`) would take both the score *and* the record of ever having computed one. You'd then be unable to tell that a re-run was needed. On `Variant` the knowledge survives.

`Variant` has no version context of its own, so the collection carries the SpliceAI details:

```python
class SpliceAIRun(TimeStampedModel):
""" #1755: one manual SpliceAI invocation - a batch of specially chosen variants.

SpliceAI is not run in bulk (~9.5s CPU/variant), so coverage is sparse and grows in unpredictable
steps: a range may be partly covered now and covered further months later. That rules out the
range-lock provenance the sweeping pipelines use (#720), so runs are recorded as collections. """
# AnnotationPipelineVersion (added for AnnotSV in #720) already records "what version of a non-VEP
# tool ran", keyed by genome_build. For SpliceAI: code_version is the `spliceai` package version
# actually invoked, data_version is the masked/raw flavour - there is no separate data bundle, the
# model weights ship with the package, so flavour is the only remaining axis.
pipeline_version = models.ForeignKey(AnnotationPipelineVersion, on_delete=PROTECT)
user = models.ForeignKey(User, null=True, on_delete=SET_NULL)

class SpliceAIRunVariant(models.Model):
""" A variant submitted to a run.

Written even when SpliceAI returned nothing, which is the point: `spliceai_max_ds IS NULL` cannot
otherwise distinguish "never ran here" from "ran, no splice effect", and SpliceAI legitimately
returns nothing for variants outside genes or longer than 2 * -D. Without this the expensive tool
is re-run forever on variants already known to score nothing, and a curator can't tell a computed
no-effect from an absence - which is most of the value of running it. """
run = models.ForeignKey(SpliceAIRun, related_name="variants", on_delete=CASCADE)
variant = models.ForeignKey(Variant, on_delete=CASCADE)

class Meta:
unique_together = ("run", "variant")
```

*"Has this variant been manually run in this build?"* is then
`SpliceAIRunVariant.objects.filter(variant=v, run__pipeline_version__genome_build=build)` — index `variant`, since that's the hot lookup on the variant page.

This tracks both things #720 asked for: **which SpliceAI version**, if they ever change it, and **masked or raw**. Flavour matters — the default settings ship the *masked* precomputed files (`spliceai_scores.masked.*`) while other configs use raw, and the two aren't comparable, so a manual run should match its deployment's flavour and record which it used. The existing `" "` convention (`_spliceai_label`, `vep_annotation.py:347`) that populates `VariantAnnotationVersion.spliceai` is the string to compare against.

### Selection

**Only variants with `spliceai_max_ds IS NULL`.** The dbNSFP-derived precomputed scores are effectively free and always win; manual runs exist purely to fill what that cache can't cover:

```python
variantannotation__spliceai_max_ds__isnull=True # no precomputed score
# and no existing SpliceAIRunVariant for this variant in this build
```

### Scores go in the normal `VariantAnnotation` columns, carried forward on version roll

The scores themselves go straight into the existing `VariantAnnotation.spliceai_*` columns — **not** a SpliceAI-only annotation store. Far simpler, and every existing consumer (`DamageNode.spliceai` filtering, grid columns, VCF export) then works with no change at all, which was the point back in #720.

Those columns are versioned, so a new `VariantAnnotationVersion` rebuilds them from VEP + dbNSFP and would wipe a manual score. At ~9.5s CPU each that isn't acceptable, so **on an annotation version upgrade, carry the previous version's manually-run scores forward** — the older rows still exist in their own partition, so this is a copy, not a re-run. `SpliceAIRunVariant` is what says which variants to copy.

That flag is also what makes the copy *safe*, and is why a bare "copy any previous SpliceAI score" would be wrong: a score that came from dbNSFP and has since been dropped from the cache should not be resurrected — only ones we computed ourselves.

Three details for whoever builds it:

* **Only carry forward a matching flavour.** Compare the run's `data_version` against the new version's `VariantAnnotationVersion.spliceai`. Carrying a *masked* manual score into a version annotated with *raw* precomputed data mixes two incomparable scales.
* **Hook it where ordering is guaranteed.** The new row only exists once the `STANDARD` VEP run for that range has imported, so the natural hook is the import lane for `STANDARD` runs — bounded by the range lock's variant ids against a small table, and skipped entirely when nothing in range is flagged. A post-hoc management command is simpler but easy to forget, and forgetting it loses scores silently.
* If a prior partition has been archived, fall back to re-running rather than silently leaving the columns null — which is possible precisely because the run record isn't in the archived partition.

## 2. The manual path (build this first)

The primary mode, and the cheap one.

* An action on the variant page — "Run SpliceAI" — for a curator looking at a variant with no splice prediction, queued to a worker.
* Accepts a small set of variants: one variant, or a list (a management command taking variant IDs / a VCF, for the "run these 40" case).
* Writes a VCF, runs `spliceai`, imports the result into the existing `spliceai_*` columns, and records a `SpliceAIRunVariant` for **every** variant submitted — including those that came back with nothing.
* Surfaces state on the variant page: not run / running / run at version X on date Y / no splice effect. A curator needs to distinguish "no score" from "not computed" — that distinction is most of the value here.

At the benchmarked ~9.5s CPU per variant this is entirely tolerable for a user-initiated lookup, and it needs none of the dispatcher machinery.

## 3. Batch mode: measure before building

The precalculated scores (via dbNSFP) cover *"all possible substitutions (snv), 1 base insertions and 1-4 base deletions (indel) within genes"*. The gap is roughly **indels of 2–50bp inside genes** — the case that started #720:

```
NC_000019.10:g.11113506_11113523del
19:11113504 GCGCTGATGCCCTTCTCTC>G LDLR, DS_AL 0.98
```

A real splice hit, invisible to us. But sweeping that gap is a different proposition from manual runs, and the benchmarks on #720 are brutal:

| | throughput |
|---|---|
| vg test, 38,889 records | 6,162 CPU-minutes ≈ **9.5s CPU/variant** |
| desktop, 4 cores / 8 threads @ 3.6GHz, 355 records | 6m31 wall (21m57 CPU) ≈ **0.9 variants/sec** |
| Quadro P600 (384 CUDA cores, 2GB), 355 records | 6m50 wall (7m01 CPU) — **no faster than 8 CPU threads**, but ~3x less CPU |

At ~0.9 variants/sec, 1M gap variants is about **two weeks of wall-clock on one box**. So this only gets built if the gap turns out small, and that has to be counted first.

The canary from #720: every dbNSFP record has `cadd_raw_rankscore`, so its absence means dbNSFP — and therefore the SpliceAI cache — didn't cover that variant.

```python
from annotation.models import AnnotationVersion, VariantAnnotationPipelineType
from annotation.annotation_version_querysets import get_variants_qs_for_annotation
from snpdb.models import GenomeBuild

av = AnnotationVersion.latest(GenomeBuild.grch38())
qs = get_variants_qs_for_annotation(av, pipeline_type=VariantAnnotationPipelineType.STANDARD,
annotated=True)
qs.filter(variantannotation__cadd_raw_rankscore__isnull=True,
variantannotation__spliceai_max_ds__isnull=True).count()
```

Break it down by variant length and in-gene/out-of-gene: SpliceAI silently skips anything longer than `2 * -D` (default 50) and returns nothing outside genes, so those are gap-but-unreachable and aren't work.

If the number is small enough to be worth it, batch mode is a `SPLICEAI` pipeline type on the existing registry, and small — the AnnotSV split did the load-bearing work:

```python
PipelineDef(SpliceAIRunner(),
depends_on=VariantAnnotationPipelineType.STANDARD,
blocks_vcf_import=False,
enabled_setting="ANNOTATION_SPLICEAI_ENABLED"),
```

with `get_variants_qs` filtering to the gap set **and excluding variants that already have a `SpliceAIRunVariant`** — which is what lets batch and manual coexist without re-running each other's work. `blocks_vcf_import=False` matters far more here than it did for AnnotSV, given the runtimes.

Even in batch mode the per-variant records still get written (one `SpliceAIRun` per batch). Range-lock provenance is not available to this pipeline in either mode, because a manual run may already have covered part of the range.

## 4. Import

SpliceAI emits VCF natively with a single `SpliceAI=ALLELE|SYMBOL|DS_AG|DS_AL|DS_DG|DS_DL|DP_AG|DP_AL|DP_DG|DP_DL` INFO field — a different shape from the per-field `SpliceAI_pred_DS_AG` etc. that VEP's plugin produces and `annotation/vep_columns.py` maps.

A small dedicated inserter modelled on `annotation/vcf_files/bulk_annotsv_tsv_inserter.py` (parse → `bulk_update` scoped by `version`) rather than extending `BulkVEPVCFAnnotationInserter`. The AnnotSV split showed that's ~200 lines and leaves the VEP insert path untouched. All the target columns are on `VariantAnnotation` (not `AbstractVariantAnnotation`), so `VariantTranscriptAnnotation` needs nothing:

`spliceai_pred_ds_{ag,al,dg,dl}`, `spliceai_pred_dp_{ag,al,dg,dl}`, `spliceai_gene_symbol`, and `spliceai_max_ds` (max of the four DS values — see `_add_spliceai_max_ds` in `bulk_vep_vcf_annotation_inserter.py`, and `VariantAnnotation.backfill_spliceai_max_ds` for the chunked-update precedent).

Writing into the same columns means existing SpliceAI filtering (`DamageNode.spliceai`, grid columns, VCF export) works with no further change — which was the point in #720.

## 5. Install and config

From #720:

```
pip install tensorflow
pip install spliceai
```

```
cd /data/annotation/fasta
wget --quiet -O - http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz | gzip -d | bgzip > hg19.fa.gz
samtools faidx hg19.fa.gz
wget --quiet -O - http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz | gzip -d | bgzip > hg38.fa.gz
samtools faidx hg38.fa.gz
```

```
spliceai -I in.vcf -O out.vcf -R /data/annotation/fasta/hg38.fa.gz -A grch38
```

CUDA needed a purge + reinstall of the nvidia drivers, then the [TensorFlow pip instructions](https://www.tensorflow.org/install/pip#linux_1), plus:

```
CUDNN_PATH=$(dirname $(python -c "import nvidia.cudnn;print(nvidia.cudnn.__file__)"))
export LD_LIBRARY_PATH=$CUDNN_PATH/lib:$CONDA_PREFIX/lib/:$LD_LIBRARY_PATH
```

Check with `python3 -c "import tensorflow as tf; print(tf.config.list_physical_devices('GPU'))"`, and `nvidia-smi` / `rocm-smi` / `nvtop` to confirm the GPU is actually being used.

Test data: [splice_ai_test.vcf.gz](https://github.com/SACGF/variantgrid/files/14900668/splice_ai_test.vcf.gz)

The public REST service (`https://spliceai-38-xwkwwwxdwq-uc.a.run.app/spliceai/?hg=38&variant=chr19-11113504-GCGCTGATGCCCTTCTCTC-G`, 15.3s for that variant) is fine for spot-checking a result by hand, but not as an implementation: #720 already notes *"they probably wouldn't like us hammering on it too much"*, and it would send variant coordinates off-site.

## 6. Known CLI behaviour to handle

* **Symbolic variants are silently skipped** — exclude them from the dump, and don't treat a missing output record as an error.
* **Variants longer than `2 * -D`** (default 50) are skipped, silently. Either bound the input by length or raise `-D`; either way the skips need accounting for, the way #1701 does for VEP, or we can't tell "scored as no-effect" from "never scored" — which is the same distinction §1 exists to preserve.
* **Variants outside genes** get no annotation. Same accounting problem.
* Whether a GPU is worth provisioning is open — the P600 result says a small card buys nothing in wall-clock. Benchmark a current card before buying one; a GPU-pinned pipeline would want its own Celery queue (related: #1653).

## Definition of done

1. `SpliceAIRun` / `SpliceAIRunVariant` (§1), plus the carry-forward on annotation version roll.
2. Manual path (§2): variant-page action + management command for a list of variants, selecting only `spliceai_max_ds IS NULL`, recording every submitted variant including no-score results, and showing not-run vs no-effect on the variant page.
3. Gap set counted per build, broken down by length and in-gene/out-of-gene — **post the numbers here before building batch mode**.
4. Batch mode (§3) only if those numbers justify it.

Start with the gap-count query using get_variants_qs_for_annotation and inspect annotation/annotation_version_querysets.py to measure the feasible scope. Then read annotation/vcf_files/bulk_annotsv_tsv_inserter.py and the existing annotation pipeline registry and version models. Done means the manual path records every submitted variant, imports SpliceAI results into the existing columns, preserves scores across compatible versions, and exposes run state to curators.

Written by the indexing model from the issue text.