# Reproduce the BRCA1 baseline comparison on a CPU

This exercise recalculates a published comparison on real experimental labels. It uses stored Evo 2 predictions, so it needs no model download, GPU or cloud account. It does not train a classifier or run Evo 2 on your laptop.

The question is narrow: does Evo 2 7B separate experimental BRCA1 loss-of-function variants better than CADD, and does the answer change when the CADD version changes?

## Run

You need Python 3.11 or newer, [uv](https://docs.astral.sh/uv/), and the `tabix` executable from [HTSlib](https://www.htslib.org/). Read the data terms below before fetching the inputs.

Save `reproduce_brca1.py`, then run:

```sh
uv run --script reproduce_brca1.py --fetch
```

The script declares its Python dependencies. It downloads a small label spreadsheet, the study's stored predictions and result summary, and an indexed region from CADD v1.7. **It does not download the 79 GB genome-wide CADD table.** Inputs go into `.brca1-cache`; all are checked against recorded hashes. The output is `brca1-reproduced.json`.

Once the inputs are cached, rerun without network access:

```sh
uv run --script reproduce_brca1.py
```

If you prefer an existing Python environment, install the versions declared at the top of the script and run it with `python reproduce_brca1.py --fetch`.

## Expected output

```text
Region     n     Evo 2   CADD1.3 CADD1.7  Evo2-CADD1.7 [95% CI]
all        3893  0.875   0.804   0.881   -0.006 [-0.025, +0.014]
coding     2768  0.840   0.841   0.892   -0.052 [-0.073, -0.030]
noncoding  1125  0.953   0.875   0.932   +0.021 [-0.003, +0.045]
```

The older CADD version makes Evo 2's noncoding result look stronger. Against CADD v1.7, the paired confidence interval includes zero. This comparison does not establish an advantage for Evo 2, or prove that the models are equivalent. CADD v1.7 separates the coding variants better in this test.

The JSON retains unrounded values, intervals, source URLs, checksums and the exact study revision. The script also checks its results against the published summary, failing if an AUROC or interval differs by more than `1e-12`.

## What the code does

1. Loads Findlay et al.'s 3,893 BRCA1 single-nucleotide variants and experimental functional classes. The positive label is `LOF`; `FUNC` and `INT` form the negative class, matching the existing study. Intermediate variants are therefore included, not silently discarded.
2. Matches labels, stored Evo 2 scores and CADD v1.7 on GRCh37 chromosome, position, reference allele and alternate allele. It stops on duplicate keys, missing scores or a changed cohort.
3. Uses `-delta` for Evo 2, where `delta` is the alternate-window mean log-likelihood minus the reference-window mean log-likelihood. Larger values therefore mean more likely loss of function. CADD uses its raw score in the same direction.
4. Calculates AUROC on all, coding and noncoding variants. Coding means Missense, Nonsense or Synonymous in the supplied consequence annotation; all other annotations form the noncoding group.
5. Resamples **positions**, keeping different alternate alleles at the same position together. It uses the same 2,000 resamples for every predictor, with seed 7. Subtracting scores from these paired resamples gives the uncertainty on the difference, rather than comparing separate confidence intervals by eye.

This is a reconstruction of [amendment 7 of the cloud study](https://github.com/rewire-bio/evo2-cloud-operations/blob/f8ae11c9a1c45505487c978b4cb7e064c4e59905/followup/CURRENT-PREDICTORS.md), restricted to CADD so the exercise avoids the much larger annotated and REVEL downloads. The original analysis also compares SpliceAI and REVEL.

## What it establishes

It checks that the stored predictions and published input labels reproduce the reported paired comparison. It shows why the baseline version, coverage, variant class and uncertainty can change a model-selection conclusion.

It is **one gene and one assay**, with no patient-level clinical validation. Most noncoding loss-of-function variants in this set are splice-related. Grouping the bootstrap by position does not make the set an independent sample of genes, remove possible training-data overlap, or establish performance on distal regulatory variants. No prediction threshold is fitted or validated here. New data and a prespecified calibration/held-out evaluation would be needed for that claim.

## Sources and data terms

- **Experimental labels:** [Findlay et al., Nature 2018](https://pmc.ncbi.nlm.nih.gov/articles/PMC6181777/). The authors make function scores available for nonprofit uses; commercial use requires a licence. The script fetches Supplementary Table 1 through the exact Arc repository copy used in the study. That repository's software licence does not relicense the experimental labels.
- **Evo 2 predictions:** the public [recorded cloud-study run](https://github.com/rewire-bio/evo2-cloud-operations/blob/f8ae11c9a1c45505487c978b4cb7e064c4e59905/evidence/runs/recorded-b75c67aa/platforms/P2/attempt-1/w1_brca1_scores.tsv), pinned to a commit and SHA-256. These are research predictions from experimental variants, not patient records.
- **CADD:** [official download and terms](https://cadd.bihealth.org/download). Scores are free for non-commercial applications; commercial applications require a licence. This exercise retrieves the GRCh37 v1.7 region `17:41197000-41277000` using tabix. The extracted rows must match the published study's SHA-256.

The supplied bundle contains code and aggregate results only. It does not redistribute the experimental label table or per-variant CADD values. Keep the downloaded cache local. Availability for download is not a licence for every use.

## Verification

Tested on 15 October 2026 with Python 3.11.13 on macOS arm64, numpy 2.3.3, pandas 2.3.3, openpyxl 3.1.5 and tabix 1.24, using a fresh isolated environment and fresh network downloads. The script took 14.7 seconds after dependency installation on this machine. A second run with cached inputs and uv's offline mode reproduced the same results in 3.5 seconds. Network and runtime will differ elsewhere. `tested-results.json` records the successful reproduction; it contains aggregates only.
