Add PositionScore command line tool for the mutational screening figure

[vmesel]

What

The synthetic mutational screening (PositionScore) only existed inside src/training_notebooks/model_training_microexons.ipynb (cells 59 and 65). This turns it into a command line tool that takes a microexon coordinate and draws the Fig. 7A panel:

python src/positionscore.py --model src/saved_model.hdf5 \
  --genome src/data/hg38.fa --conservation src/data/hg38_cons.bw \
  --exon chrX:31126642:31126673:+ --gene DMD --panel-label A \
  --output DMD_positionscore.png --tsv DMD_positionscore.tsv

Every base of the two 100 nt flanks is replaced by each of the four nucleotides, the microexon is re-scored with the CNN, and the delta (mutant - wild type) is drawn as a 4 x 200 heatmap with the microexon diagram below it, keeping the settings of the published figure (RdBu_r, scale from -0.015 to +0.015, TGAC row order).

Closes the "Include the mutational screening in the command line tool" item of the README roadmap.

Changes

• src/deepmex.py: the script body moves into main(), so Microexon can be imported. The CLI, its arguments and its output are unchanged. New predict_batch() scores many sequence variants in a single model call, and prediction() now delegates to it, so both paths handle strands identically. Unknown bases (N) no longer raise a KeyError, and the tool no longer crashes when stdout is redirected to a file, which is how the README documents it (> result.out).
• src/positionscore.py: the screening, the figure, an optional TSV export of the per position / per base scores, and --demo to inspect the layout without the model or the reference files.
• tests/test_positionscore.py: checks against a stand in model, including the batched screening against the position by position loop of the notebook. Runs without tensorflow and without the 8 GB of reference files (make test).
• README.md / Makefile / requirements.txt: docs, make run-positionscore-example, make positionscore-demo, make test, and matplotlib.

Notes for review

• Minus strand exons are drawn in transcript orientation (column 0 is always -100 relative to the exon). The flank handling matches what deepmex.py already does when scoring a minus strand microexon: flanks are swapped and sequences are read backwards, without complementing the bases. Worth a second look at some point, but this PR deliberately keeps the existing behaviour so the figure agrees with the score the tool already reports.
• Only panel A of Fig. 7 is generated. Panel B (PSI distributions of the PTBP1 knock-down) and panel C (UCSC genome browser tracks) come from external data, not from this pipeline.
• The 800 variants are scored in batches instead of one prediction per variant, which is what made the notebook version slow.
• make download-data saves the bigwig as src/data/hg38.phastCons100way.bw while make run-example reads src/data/hg38_cons.bw; that mismatch predates this PR and is left untouched.

[vmesel]

Pushed four more commits. Two of them change behaviour and want a second opinion.

deepmex-figure: coordinate in, figure and impact ranking out

uv run deepmex-figure chrX:31126642-31126673

Strand and gene name come from RefSeq, the flanks from the UCSC API, so neither the reference files nor bedtools are needed. It writes the figure and a table of every flank base ranked by how much the model reacts to mutating it, each row carrying the genomic position, the position relative to the exon, the reference base and the delta of each substitution. For the DMD microexon the top positions are the thymines of the polypyrimidine tract, 14 to 22 nt upstream:

  #      position    rel  wt    impact   deltas
  1      31126694    -20   T     1.566   A-0.81  C-1.05  G-1.57
  2      31126690    -16   T     1.395   A-0.98  C-1.07  G-1.40

--genome/--conservation still read local files when given.

Panels A and B in one figure

Panel B (inclusion under a knock-down) is not model output, so it is read from a vast-tools INCLUSION_LEVELS_FULL table: the densities are the Beta posteriors rebuilt from the corrected read counts of the quality column, and the right hand curve is P(|dPSI| > x) with the 95% MV[dPsi] marked. Panel C is a UCSC browser screenshot and stays out of scope.

Two fixes that change what the tool outputs

The published panel is the reverse complement, not just the reverse. Extracting panel A from the Fig. 7 TIF and correlating it against the screening of the same microexon identifies the transform: r goes from 0.13 to 0.67 once the bases are complemented, and every other candidate transform scores below 0.21. Independent check: with the fix, plus and minus strand exons show the same signature upstream (purines negative, pyrimidines positive) across 38 microexons; before it, the 19 minus strand ones were inverted.

--vmax was off by a factor of 100. The 0.015 of the paper is on the 0-1 scale of the notebook, while the CLI reports 0-100. Real deltas reach ~1.0, so every figure came out fully saturated. The default is now 1.5.

Also: the DMD example passed strand +, but DMD is on the minus strand (16 RefSeq transcripts, all -).

Open question, not touched

The flanks of the published figure start one base earlier than the code computes. Correlation with the extracted panel goes 0.67 → 0.87 with the shift, and the deciding half is the downstream flank: the code reads [end+1, end+101), skipping the base immediately after the exon, where [end, end+100) looks right. I have left it alone because it changes the flanks of every prediction, not just the drawing (76.12 → 76.50 for DMD). Same question as the strand handling raised earlier — worth a look before either is changed.