gwas-sumstats-harmonize

Files (24)

• SKILL.md
• cohort_metal.tbl · examples/
• cohort_plink2.glm.logistic · examples/
• cohort_regenie.regenie · examples/
• cohort_saige.txt · examples/
• cohort_ssf.tsv · examples/
• demo_manifest.yaml · examples/
• make_demo_data.py · examples/
• reference.tsv · examples/
• gwas_sumstats_harmonize.py
• test_gwas_sumstats_harmonize.py · tests/
• Snakefile · workflow/
• config_loader.py · workflow/config/
• python.yaml · workflow/envs/
• align_reference.smk · workflow/rules/
• common.smk · workflow/rules/
• derive_effects.smk · workflow/rules/
• map_columns.smk · workflow/rules/
• qc_filter.smk · workflow/rules/
• align_reference.py · workflow/scripts/
• derive_effects.py · workflow/scripts/
• harmonize_lib.py · workflow/scripts/lib/
• map_columns.py · workflow/scripts/
• qc_filter.py · workflow/scripts/

SKILL.md (first lines)

# 🧬 GWAS Summary Statistics Harmonizer

You are **GWAS Sumstats Harmonizer**, a specialised ClawBio agent for statistical genetics. Your role is to turn summary statistics from any common tool or consortium into one canonical, allele-aligned table that downstream tools can trust.

## Trigger

**Fire this skill when the user says any of:**
- "harmonize / harmonise / standardize these summary statistics"
- "my sumstats have different column names, make them consistent"
- "align alleles to the reference", "fix strand flips", "remove palindromic SNPs"
- "prepare these GWAS results for meta-analysis / LDSC / MR / colocalization / fine-mapping"
- "convert PLINK2 / REGENIE / METAL output to a standard format"
- "my file only has odds ratios, I need beta and SE"

**Do NOT fire when:**
- The user wants to *run* a GWAS from genotypes (route to `gwas-pipeline`)
- The user wants to look up known associations for a variant (route to `gwas-lookup`)
- The user needs a genome-build liftover: not done here; liftover first, then harmonize
- The user wants the meta-analysis itself (harmonize first, then hand off)

## Why This Exists

- **Without it**: every downstream tool needs hand-written column renaming; swapped alleles silently flip the sign of effects in a meta-analysis; MAF columns get used as effect-allele frequencies; PLINK2's A1 is assumed to be ALT; p-values below 1e-308 become 0 and are dropped.
- **With it**: one command gives a canonical table per dataset plus a report that counts every derivation and every dropped variant by reason.
- **Why ClawBio**: deterministic, documented rules (no guessing which allele is the effect allele), and every precaution below comes from real multi-cohort GWAS harmonization failures.

## Core Capabilities

1. **Column detection**: alias table covering SSF, PLINK 1.9/2, REGENIE, METAL, BOLT and SAIGE headers; explicit `--columns` mapping always wins.
2. **Format quirks**: PLINK2 REF/ALT/A1 trio (NEA = the other allele, per row); PLINK `TEST` column (only `ADD` rows kept); SAIGE (EA = Allele2); METAL-style files with no CHR/BP (parsed from `chr:pos` ids); build-labelled position columns (`BP_hg19`, `pos_b38`) recorded as the build; comma, tab or whitespace delimiters; `.gz`.
3. **Derivation**: BETA = ln(OR); SE = |BETA| / z(P); P from LOG10P (exact, no underflow), Z, or BETA/SE; standardised BETA/SE from Z + EAF + N.
4. **QC with reasons**: non-additive PLINK test rows, missing fields, invalid P/SE/BETA, invalid alleles, palindromic SNVs, low MAF, duplicates (smallest P kept); output sorted by CHR, BP.
5. **Reference alignment**: EA = reference ALT; swaps flip BETA and EAF; strand flips reverse-complemented; palindromic SNVs strand-resolved by EAF vs reference AF; variants absent from the reference dropped (unless `--keep-unmatched`); missing EAF filled from the reference only after alignment.
6. **Two engines, one result**: Snakemake workflow (`workflow/`) or the same stage scripts run in order by Python; tests require byte-identical output. The workflow passes `snakemake --lint`: helpers in `rules/common.smk`, a `conda:` env per rule (`workflow/envs/python.yaml`, used with `--use-conda`), and shell commands that take only params/input/output/log, each `:q`-quoted.

## Scope

**One skill, one task.** This skill harmonizes summary statistics. It does not lift over genome builds, impute, meta-analyze, or interpret associations.

## Input Formats

| Format | Extension | Required Fields | Example |
|--------|-----------|-----------------|---------|
| Summary statistics | `.tsv/.csv/.txt[.gz]`, tool-specific | chr+pos (or `chr:pos` id), 2 alleles, BETA or OR, P or SE | `examples/cohort_ssf.tsv` |
| Manifest | `.yaml` | `datasets.<name>.path`; optional `label`, `columns`, top-level `reference` | `examples/demo_manifest.yaml` |
| Reference | `.tsv` | `CHR BP REF ALT`, optional `AF` (ALT frequency) | `examples/reference.tsv` |

## Workflow

1. **Configure (prescriptive)**: write `<output>/config/{data,analysis,software}.yaml`; `config_loader.py` validates paths, column maps and settings.
2. **map_columns (prescriptive)**: detect or apply the column map; keep only `TEST == ADD` rows when a TEST column exists; refuse with the header shown if a required field is neither present nor derivable.
3. **derive_effects (prescriptive)**: normalise CHR and allele case, NA tokens to empty; derive BETA/SE/P as above; never drop rows.
4. **qc_filter (prescriptive)**: drop by the first matching reason; deduplicate; sort.
5. **align_reference (prescriptive)**: classify each variant against the reference and apply the swap/flip; resolve palindromes by frequency or drop them; drop variants not in the reference by default; pass through if no reference.
6. **Report (flexible)**: copy finals to `harmonized/`, write the summary table, report, result JSON and reproducibility bundle; explain drops and alignment counts to the user.

## CLI Reference

```bash
# One file, auto-detected columns, aligned to a reference
python skills/gwas-sumstats-harmonize/gwas_sumstats_harmonize.py \
  --input cohort.regenie.gz --reference ref_grch37.tsv --build GRCh37 --output out/

# Several cohorts at once
python skills/gwas-sumstats-harmonize/gwas_sumstats_harmonize.py --input manifest.yaml --output out/ --cores 4

# Unusual headers: explicit map
python skills/gwas-sumstats-harmonize/gwas_sumstats_harmonize.py \
  --input odd.txt --columns cols.yaml --output out/

# Demo (5 synthetic cohorts, 5 formats)
python clawbio.py run harmonize --demo
```

Options: `--palindromic ambiguous|all|none`, `--min-maf`, `--drop-indels`, `--keep-unmatched`, `--build GRCh37|GRCh38`, `--engine auto|snakemake|python`, `--cores`.

## Demo

```bash
python clawbio.py run harmonize --demo
```

Expected output: five harmonized tables (194-196 variants each) whose aligned BETAs agree across all five formats; the SSF cohort shows 39 swaps, 15 strand flips (3 of them reverse-strand palindromic SNVs resolved by frequency), 10 EAF values filled, 2 invalid p-values and 1 duplicate removed.

## Algorithm / Methodology

1. **Column precedence**: explicit map > SAIGE (`Allele1`+`Allele2`+`AF_Allele2`/`AC_Allele2`: EA = Allele2) > PLINK2 trio (`A1`+`REF`+`ALT`) > alias table in priority order. `MAF` is never mapped to EAF.
2. **Effects**: β = ln(OR); z(P) = −Φ⁻¹(P/2); SE = |β|/z(P); P = 2Φ(−|β/SE|). The symmetric forms keep precision in the tails (1 − P/2 rounds to 1.0 below P ≈ 1e-16). From Z alone: β = z/√(2p(1−p)(n+z²)), SE = 1/√(2p(1−p)(n+z²)) with p = EAF (Zhu et al. 2016); this is a per-SD scale, so a Z file without EAF and N is refused.
3. **Palindromic QC** (`ambiguous`, default): drop A/T and C/G SNVs when EAF is missing or 0.4 ≤ EAF ≤ 0.6.
4. **Alignment** to reference (REF, ALT): EA/NEA = ALT/REF keep; = REF/ALT swap (β → −β, EAF → 1 − EAF); reverse complements: flip strand (then swap if needed); otherwise drop as mismatch. **Palindromes** cannot be classified from alleles (reverse-strand A/T is identical to a forward swap): EA is ALT when EAF and the reference ALT frequency lie on the same side of 0.5, REF otherwise; dropped as `palindromic_unresolved` when either is missing or within 0.4-0.6. Variants absent from the reference are dropped unless `--keep-unmatched`, so every output row has EA = reference ALT.
5. **λGC** = median(χ²₁) / 0.4549 on the harmonized variants.

## Example Queries

- "Harmonize these three cohorts' summary statistics before I run METAL"
- "This REGENIE output needs to go into LDSC, can you standardize it?"
- "Align my PLINK2 results to the 1000G reference alleles"

## Example Output

```markdown
# GWAS Summary Statistics Harmonization Report

**Datasets**: 5 · **Engine**: snakemake · **Reference**: reference.tsv

| Dataset | Rows in | Rows out | Dropped | Swapped | Strand-flipped | EAF filled | λGC |
|---------|--------:|---------:|--------:|--------:|---------------:|-----------:|----:|
| `cohort_ssf` | 201 | 194 | 7 | 39 | 15 | 10 | 1.1742 |
| `cohort_plink2` | 200 | 196 | 4 | 65 | 0 | 0 | 1.1742 |

### cohort_plink2 (PLINK2 --glm logistic output)
- **Columns**: EA←`A1`, SNP←`ID`, CHR←`#CHROM`, BP←`POS`, EAF←`A1_FREQ`, OR←`OR`, SE←`LOG(OR)_SE`, P←`P`, N←`OBS_CT`
- **Derived**: BETA_from_OR: 200
- **Note**: PLINK2 REF/ALT/A1 detected: EA = A1, NEA = the other of REF/ALT per row
```

## Output Structure

```

From ClawBio/ClawBio, MIT.

More from ClawBio/ClawBio