# Cohort-Stratified Meta-Analysis

VariantCentrifuge can combine two or more clinically distinct case-control
**strata** into one calibrated rare-variant discovery analysis using
**score-statistic meta-analysis** (issue #112). Unlike pooling with a cohort
covariate, each stratum keeps its own null model, ascertainment, and case:control
ratio; unlike combining per-gene p-values, the signed score and covariance are
preserved so burden direction and kernel (SKAT) evidence are retained.

> **Scope.** This is an analysis-engine feature for **binary traits**. Study-specific
> phenotype eligibility, covariate choice, and QC remain caller responsibilities.
> Quantitative traits, a mixed-model/GRM null, and a genotype-level saddlepoint
> (SPA) for *severe* imbalance are deferred (see *Limitations*).

## How it works

The workflow mirrors RAREMETAL / MetaSKAT / seqMeta: each stratum exports a
compact **score/covariance artifact**, then a separate step combines them.

1. **Per-stratum artifact.** For every gene, VariantCentrifuge stores the
   unweighted score vector `U = Gᵀr` (r = null residuals) and the
   projection-adjusted covariance `V = GᵀPG`, plus per-variant observed ALT-allele
   counts (`AC`) and allele numbers (`AN`). These are gene-level **aggregates** —
   no per-sample genotypes or phenotype values — so artifacts are safe to move
   between sites.
2. **Meta step (`variantcentrifuge meta`).** For each gene it aligns the strata on
   a canonical `chrom:pos:ref:alt` key (union of variants; a variant absent in a
   stratum contributes zero information there), harmonises weights from the
   **pooled** frequency, and combines:
   - **meta-burden** — common-effect score burden (`Z = 1ᵀS / √(1ᵀΣ1)`), with a
     signed direction;
   - **meta-SKAT** — kernel test via the Davies mixture chain;
   - **meta-SKAT-O** — the optimal omnibus (the reported *primary* result);
   - **heterogeneity** — Cochran's Q, I², per-stratum directions, and
     leave-one-stratum-out.
   Exome-wide multiplicity correction (BH/Bonferroni) is applied to the primary
   p-value across genes.

Weights are harmonised from the pooled ALT frequency `ΣAC/ΣAN`, using
`MAF = min(af, 1-af)` for the `Beta(MAF; 1, 25)` weight. This is essential: a
stratum's own MAF would put its `U`/`V` on an incompatible scale.

### ACAT-O vs meta-analysis

These are **different** operations and are kept distinct:

- **Within-gene ACAT-O** (see the Association Testing guide) combines *tests/masks
  within one cohort* via Cauchy combination.
- **Cross-stratum meta-analysis** (this guide) combines the *score and covariance*
  across cohorts. Per-cohort p-value combination is **never** used as a
  substitute for the score-based common-effect test.

## Creating a cohort artifact

Add `--score-artifact <path>` to a normal association run of each cohort
(binary trait, case/control samples set):

```bash
variantcentrifuge ... --perform-association --association-tests skat \
  --score-artifact /data/meta/adult.npz
```

This writes `adult.npz` and `adult.npz.meta.json` next to the usual results. It
uses the same site/sample missingness QC as the association tests but keeps genes
below the per-gene MAC filter, because a gene too rare to test in one cohort still
adds information to the meta. The artifact records the reference (`--reference`)
as build and annotation version, a hash of the variant selection (`--filters`,
`--final-filter`, `--bcftools-prefilter`) as the mask, and SHA-256 hashes of the
sample IDs for the overlap check. Run every cohort with the same reference and
filters, or the meta step will refuse to pool them.

## Study manifest

A TSV, one row per stratum:

```
stratum_id	artifact_path
adult_cohort	/data/meta/adult.npz
pediatric_cohort	/data/meta/pediatric.npz
```

Relative paths are read relative to the manifest file. The `stratum_id` names the
cohort in the output (e.g. `per_stratum_directions` = `adult_cohort:1,pediatric_cohort:-1`).

The meta step **refuses** (hard error, no silent fallback) when strata disagree on
genome build, annotation version, mask definition, or trait type, or when their
(hashed) sample sets overlap — overlap would violate the additive-covariance
independence assumption. Sample non-overlap does not by itself establish
independence for related or shared-control populations; that remains a caller
assumption.

## Running it

```bash
variantcentrifuge meta \
  --manifest study.tsv \
  --output meta_results.tsv \
  --weights beta:1,25 \
  --correction fdr
```

Only `beta:a,b` and `uniform` weights are supported (the only schemes reproducible
from the artifact). Output columns include `meta_burden_p`,
`meta_burden_direction`, `meta_skat_p`, `meta_skato_p`, `meta_skato_rho`, `het_q`,
`het_p`, `het_i2`, `per_stratum_directions`, `loo_min_p`, `primary_p`,
`primary_q`, and a `status`/`warnings` pair (e.g. a non-converged stratum is
excluded with an explicit warning).

## API

```python
from variantcentrifuge.association.meta import (
    compute_gene_score, meta_burden, meta_skat, meta_skato,
)
```

`compute_gene_score(imputed_matrix, raw_dosage, null_model, variant_keys)` builds a
`GeneScore`; the `meta_*` functions combine aligned strata. A single-stratum call
reproduces the single-cohort SKAT/SKAT-O/burden result exactly — the correctness
gate that ties the meta engine to the validated backend.

## Calibration and limitations

The common-effect score meta uses a normal/Davies approximation that is well
calibrated under **moderate** imbalance. Under **severe** case:control imbalance the
per-variant/burden normal approximation is anti-conservative in the tail; a
calibrated result there requires a genotype-level saddlepoint (SPA) backend, which
is **deferred**. The artifact stores the raw `U`/`V`, so an SPA meta path can be
added without changing the schema. Likewise, `V` is stored as a full matrix, so a
future relatedness-aware (GLMM/GRM) null drops in without a schema change.

The meta engine is cross-validated against R MetaSKAT 0.90 (`MetaSKAT_wZ`,
`combined.weight=TRUE`, `weights.beta=c(1,25)`): burden and SKAT agree to ~1e-15 and
SKAT-O to ~1e-4 (`tests/unit/test_r_parity_skato.py` stores the R values, so CI needs
no R). Below p≈1e-5 the per-rho tail probabilities use a saddlepoint approximation,
whereas R falls back to the Liu approximation there; Liu is anti-conservative in that
range, so very small p-values can be an order of magnitude or more larger (more
conservative) than R's.