---
name: gwasengine
description: >
  GWASEngine: Complete pure-Python GWAS analysis pipeline. Use for: GWAS analysis,
  polygenic risk scores (PRS), LD clumping, fine-mapping, LD Score Regression (LDSC),
  SNP heritability estimation, quality control of genotype data. Triggers on:
  "GWAS", "genome-wide association", "PRS", "polygenic risk score", "fine-mapping",
  "LDSC", "heritability", "Manhattan plot", "QQ plot", "association study",
  "clumping", "causal variant", "SNP", "genotype".
---

# GWASEngine — Pure Python GWAS Analysis

> **Python**: Use `/torch/venv3/pytorch/bin/python3` — numpy, scipy, pandas, scikit-learn, plotly installed.

## Environment Check

```bash
/torch/venv3/pytorch/bin/python3 -c "import numpy, scipy, pandas, sklearn, plotly; print('all deps ok')"
```

## Core API

```python
from gwasengine import run_gwas_engine

summary = run_gwas_engine(
    trait="quantitative",      # or "binary"
    out_dir="gwas_output",
    n_samples=2000,
    n_snps=10000,
    h2_snp=0.3,               # SNP heritability
    n_causal=20,              # causal variants
    n_genetic_pcs=4,          # PCs for stratification correction
    run_prs=True,
    run_finemapping=True,
    run_ldsc=True,
)
```

## Module Reference

### `generate_synthetic_gwas(n_samples, n_snps, n_causal, h2_snp, rng_seed)`

Generate synthetic genotype-phenotype data with 3 ancestry clusters, LD block structure (Cholesky), and causal SNPs with realistic effect sizes.

**Returns** `GWASData(n, m)` namedtuple: `G` (n×m genotype matrix, additive coded 0/1/2), `phenotype` (vector), `pos` (positions), `snp_ids`, `chrom`, `causal_mask`.

### `quality_control(data, min_maf=0.01, min_hwe=1e-6, min_call_rate=0.95, remove_relateds=True)`

Sample + variant QC. Returns `(data_qc, snp_mask, sample_mask, qc_report)`.

### `run_gwas_linear(data, add_genetic_pcs=4)`

Univariate linear regression per SNP. Handles n < m via residualization. Returns `(results_df, lambda_gc)`.

### `ld_clumping(gwas_results, G, pos, r2_threshold=0.1, p_threshold=5e-5)`

Gabriel LD clumping. Returns lead SNPs DataFrame.

### `compute_prs(G, snp_ids, gwas_results, p_threshold=1e-3)`

Clumping+Thresholding PRS. Returns per-sample polygenic scores.

### `fine_map_locus(gwas_results, locus_snps, n_samples)`

Wakefield Approximate Bayes Factors → PIPs + 95% credible set.

### `ld_score_regression(gwas_results, G, pos, n_samples)`

LDSC: regress χ² on LD scores → h²_SNP estimate + intercept.

### `visualize_gwas(gwas_results, lambda_gc, lead_snps, prs, phenotype, ldsc_result, fine_map_df, output_path)`

Plotly 6-panel HTML dashboard: Manhattan, QQ, PRS histogram, PIPs, LDSC, Top Hits.

## Quick Demo

```bash
cd /root/gwasengine
/torch/venv3/pytorch/bin/python3 -m gwasengine
# → generates gwas_output/gwas_analysis.html
```

## Output Files

```
gwas_output/
├── gwas_analysis.html   # 6-panel interactive dashboard
├── gwas_results.csv     # all SNP statistics (beta, se, p_value, r2)
├── lead_snps.csv        # independent genome-wide significant loci
├── prs.csv             # per-sample polygenic risk scores
├── finemapping.csv     # PIPs and 95% credible set
└── summary.json        # pipeline metrics
```

## Expected Results (2000 samples, 10k SNPs, h²=0.30, 20 causal)

| Metric | Expected |
|--------|----------|
| Top p-value | < 5×10⁻⁸ |
| λ_GC | 0.96–1.10 |
| PRS R² | 0.05–0.15 |
| LDSC h²_SNP | 0.25–0.35 |
| Pipeline time | ~30s CPU |

## Key Implementation Details

- **QC**: MAF ≥ 0.01, HWE p ≥ 1×10⁻⁶, call rate ≥ 95%. PCA on genome-wide SNPs for ancestry correction.
- **GWAS**: O(n×m) vectorized linear regression. Covariate residualization handles n < m.
- **LD**: Pairwise r² computed in 500kb windows. Gabriel block detection.
- **PRS**: C+T at p < 10⁻³ threshold. Effect sizes from GWAS betas.
- **Fine-map**: Wakefield ABF: `BF ≈ sqrt(N)/|z| × exp(χ²/2)`. Prior on causality = 1e-4.
- **LDSC**: `χ² = 1 + N×h²/M×LD + a`. Intercept > 1 indicates confounding.