# Linear / mean baselines on the VCC-2025 H1 data — GlioState own analysis (draft, pseudobulk only)

Run date 2026-10-16 (Pacific Time). Purpose: check, on one public dataset, how three baseline families (control mean, training mean, a linear model with a PCA-derived perturbation embedding) behave under different metric and reference choices, with a technical-duplicate positive control and a shuffled-embedding null. It is **not** a peer-reviewed result and evaluates **no** deep-learning model.

## Data
The same public 2025 Virtual Cell Challenge files as `../vcc_baseline/` (H1 hESC CRISPRi Perturb-seq; train 150 / validation 50 / test 100 perturbations; 18,080 genes; control `non-targeting`, 38,176 control cells in the training file). Not redistributed. See `../vcc_baseline/README.md` for the download and for `01_pseudobulk.py`, which writes `train.npz`, `val.npz`, `test.npz` (means of all cells, half 0 and half 1, per target).

## Reproduce
```
python3 -m venv venv && venv/bin/pip install numpy scipy
venv/bin/python 03_linear_baselines.py DATADIR .      # DATADIR holds train.npz val.npz test.npz; ~2 min CPU; writes results.json
```

## What is computed
- Truth for every predictor = half-0 pseudobulk of each held-out perturbation (validation 50 + test 100 = 150 perturbations; "pooled" = all 150). Control reference = mean of all control cells.
- Predictors: control mean; mean of the 150 training perturbation means (perturbation-weighted; VC0048 'mean', cell-eval baseline construction); cell-weighted training mean (VC0052 'perturbed mean'); linear model of VC0048 (Y ≈ G W Pᵀ + b, ridge λ = 0.1, W by the closed form in its Methods; G = PCA scores of the centred training matrix with genes as observations, P = rows of G for the training target genes) for K = 1, 2, 5, 10, 20, 50; the same with the target gene set to the median training knockdown ratio (0.095 × control); technical duplicate (half 1 as prediction; positive control, VC0146 definition); and a **shuffled-embedding null** (K = 10 model, but each held-out perturbation is given another held-out perturbation's target-gene embedding; 20 shuffles averaged).
- Metrics per perturbation: MSE over genes; L2 distance on the 1,000 genes most highly expressed in control (VC0048's error readout); Pearson delta with control reference and with the training-mean reference (VC0052 uses a perturbed-centroid reference; ours is the perturbation-weighted training mean); PDS with L1 and L2 (1 − rank/N within validation or test; target gene excluded). Each metric is computed with the target gene excluded and with it included.
- Paired bootstrap (4,000 resamples over the 150 perturbations) for linear K = 10 minus training mean, and for linear K = 10 minus its shuffled null.
- Dynamic-range fraction (VC0146 definition) with the technical duplicate as positive and the training mean as negative control (control mean as negative control for the training-mean-referenced Pearson delta, where the training mean itself has no defined delta).

## Caveats (read before citing)
- Our re-implementation. The PCA convention for G is our reading of "PCA on Y_train" in VC0048's Methods; the authors' code was not run. A different convention could change the linear-model results.
- Pseudobulk only; not the official cell-eval; no differential-expression metrics (so no weighted MSE, no interpolated duplicate, no DEG-restricted Pearson); no deep-learning model evaluated.
- One cell line (H1), 2025 in-context task with 150 training perturbations; nothing here transfers automatically to K562/RPE1 data or to the 2026 zero-shot task.
- Truth is a half-cell pseudobulk, so all scores include sampling noise; the technical duplicate is a noisy positive control (two half-cell estimates).
- One random split of cells (seed 0); confidence intervals resample perturbations only and ignore cell-split variance.
- Perturbation-weighted training mean is the reference for "Pearson delta (training-mean reference)"; this differs from VC0052's perturbed-centroid reference computed on all perturbed cells.

---
# Update 2026-10-17 (Pacific Time): multi-split repeat, K562/RPE1, official cell-eval

## 1. 20 random perturbation-level splits (`04_multi_split.py`, `06_summarize.py`)
`python 04_multi_split.py DATASET N_TRAIN N_SPLITS OUT.json [PROCS]` with DATASET h1 | k562 | rpe1; `python 06_summarize.py` writes `multi_split_summary.json` (mean, SD, min, max; paired differences with counts of splits in each direction).
- **What is split:** which perturbations are training (N_TRAIN) and which are held out (150). Split seeds 1000..1019. The cell-level split is NOT repeated (H1: the seed-0 halves of `01_pseudobulk.py`; truth for held-out perturbations = half 0, technical duplicate = half 1; training matrix = all-cell means).
- **H1 has only 300 perturbations** (train 150 + validation 50 + test 100), pooled; each split re-partitions the same 300 into 150/150, so splits overlap heavily and the SD reflects sensitivity to which perturbations are held out, not sampling of new perturbations. Control = non-targeting mean of the training file.
- **K562 and RPE1:** Replogle et al. 2022 `K562_essential_raw_bulk_01.h5ad` (79.8 MB) and `rpe1_raw_bulk_01.h5ad` (95.4 MB) from figshare article 20029387 (public; not redistributed here). Rows = per-guide-set mean UMI counts per cell; we normalise each row to 10,000 counts and log1p, average rows with the same target gene, average non-targeting rows (>= 25 filtered cells) as control, keep perturbations whose target gene is a measured gene and with >= 25 cells: 1,821 (K562, 8,563 genes) and 1,851 (RPE1, 8,749 genes). Truth = that pseudobulk (no half split, so **no technical duplicate**). Splits: 150 train / 150 held-out, and 900 train / 150 held-out (`*_ntrain900.json`).
- Predictors: control mean, training mean, linear model K=10 and K=50 (same closed form/PCA convention as `03_linear_baselines.py`), K=10 shuffled-embedding null (5 shuffles per split, averaged), target-gene overwrite (median training knockdown ratio x control) applied to the mean and to K=10. Metrics: MSE, Pearson delta (control and training-mean reference), PDS-L1/L2 (N = 150, 1 - rank/N, target excluded); MSE and Pearson delta also with the target included.
- Runtime: ~1 min per H1 split on 1 core; full set including 900-train runs about 1.5 h with 4 processes.

### Key results (mean ± SD over 20 splits; target excluded; MSE x1e-3)
| Dataset | Predictor | MSE | Pearson delta (ctrl) | PDS-L1 | PDS-L2 |
|---|---|---|---|---|---|
| H1 150/150 | training mean | 3.74±0.44 | 0.161±0.012 | 0.506 | 0.553 |
| | linear K=10 | 4.04±0.39 | 0.111±0.015 | 0.508±0.005 | 0.552±0.009 |
| | linear K=50 | 4.61±0.31 | 0.089±0.012 | 0.513±0.006 | 0.556±0.010 |
| | shuffled null | 4.10±0.39 | 0.100±0.012 | 0.506±0.002 | 0.552±0.006 |
| | technical duplicate | 1.13±0.11 | 0.649±0.017 | 0.953 | 0.990 |
| K562 150/150 | training mean | 3.36±0.20 | 0.232±0.012 | 0.505 | 0.511 |
| | linear K=10 | 3.39±0.20 | 0.222±0.016 | 0.521±0.008 | 0.530±0.010 |
| | shuffled null | 3.59±0.22 | 0.196±0.010 | 0.505 | 0.511 |
| RPE1 150/150 | training mean | 7.86±0.53 | 0.476±0.032 | 0.504 | 0.507 |
| | linear K=10 | 7.98±0.44 | 0.466±0.032 | 0.524±0.007 | 0.530±0.009 |
| | shuffled null | 8.44±0.52 | 0.459±0.031 | 0.504 | 0.507 |
| K562 900/150 | training mean / linear K=10 | 3.4 / 3.3 (K=10 lower in 20/20 splits) | 0.243 / 0.266 | 0.504 / 0.525 | 0.510 / 0.535 |
| RPE1 900/150 | training mean / linear K=10 | 7.8 / 7.5 (lower in 19/20) | 0.486 / 0.490 | 0.504 / 0.532 | 0.507 / 0.541 |

Target-gene overwrite, Pearson delta with the target included: H1 training mean 0.156 -> 0.344; K562 0.231 -> 0.253; RPE1 0.475 -> 0.481. Full tables: `multi_split_summary.json`.

### What changed relative to the single-split conclusions above
- H1: reproduced — linear K=10/K=50 worse than the training mean on MSE (20/20 splits) and control-referenced Pearson delta (20/20); K=10 PDS equal to the mean. Softened — "shuffled-embedding null indistinguishable" becomes "K=10 slightly better than its null on Pearson delta (+0.011±0.010; 18/20 splits) and MSE (-0.07e-3; 19/20)", i.e. very weak embedding information.
- K562/RPE1 (different technology, normalisation, gene set, panel): linear K=10 ~ mean on MSE and ~0.01 lower on Pearson delta at 150 training perturbations, but higher than both mean and null on PDS in 20/20 splits (+0.02); at 900 training perturbations it also beats the mean on MSE. The general statement "a linear model built from training data does not beat the mean on any metric" is therefore **not** supported beyond H1 / these settings.
- Target-gene overwrite effect is large on H1, small on K562/RPE1.
- Not tested: whether the K562/RPE1 vs H1 difference is due to cell line, data processing or panel.

## 2. Official cell-eval (`05_cell_eval_prep.py`, `08_run_cell_eval.sh`, `07_cell_eval_xcheck.py`, `cell_eval_results.json`)
cell-eval 0.8.2 (`pip install cell-eval`, Python 3.13 venv; installs without problems) was run with `--profile full` on single-cell input built from the H1 *training* h5ad (the 15.5 GB file was downloaded to /workspace/replogle for this and is not in the repo; only 150 perturbations have cells there because the validation/test files were not used). One split only: 100 perturbations train the predictors, 50 are held out with <= 100 real cells each, 1,000 real control cells; predicted profiles are replicated per cell (cell-eval baseline construction); pseudobulk predictors as above, linear predictions clipped at 0 because cell-eval rejects negative log1p values. Rerun note: a first run had a row/label misalignment in the prep script (rows sorted in file order, labels in pick order); it was found because the median real knockdown ratio was 0.99, fixed (`np.sort`), and all numbers below come from the fixed rerun (real target ratio median 0.089).
Our own pseudobulk MSE and Pearson delta computed from the same h5ad files agree with cell-eval's to 4 digits (train mean: Pearson delta 0.2477, MSE 4.31e-3 in both), which validates the re-implementation of these two metrics.

| Predictor (50 held-out perts) | MSE x1e-3 | Pearson delta | PDS-L1 | PDS-L2 | DE direction match | DE sig-gene recall |
|---|---|---|---|---|---|---|
| control mean | 4.38 | 0.226 | 0.515 | 0.556 | 0.716 | 0.418 |
| training mean | 4.31 | 0.248 | 0.515 | 0.564 | 0.768 | 0.555 |
| linear K=10 | 5.25 | 0.159 | 0.515 | 0.559 | 0.684 | 0.582 |
| linear K=50 | 6.54 | 0.132 | 0.508 | 0.547 | 0.673 | 0.628 |
| K=10 shuffled null | 5.23 | 0.164 | 0.540 | 0.578 | 0.683 | 0.559 |
| training mean + target overwrite | 4.18 | 0.319 | 0.515 | 0.564 | 0.795 | 0.595 |

Reading: ordering on MSE and Pearson delta matches the pseudobulk analysis (training mean best; linear worse; K=10 and its null equal). On cell-eval's PDS the shuffled null is even above the real K=10 (single split, 50 perts: noise-level). Linear models score higher on DE-recall metrics (0.58–0.63 vs 0.56) but that metric is dominated by the number of called DE genes in the prediction and we did not analyse it further. One split of 50 perturbations from the training file is not comparable to the 150-held-out tables (different perturbations, fewer cells in the truth, different split); no CI. cell-eval's `pearson_delta` includes the target gene, which is why overwrite helps. We did not run `--ceiling`, nor the official leaderboard configuration (`cell-eval prep`, vcc profile).