bio-data-visualization-manhattan-qq-locuszoom

What this skill does

## Version Compatibility

Reference examples tested with: qqman 0.1.9 (R), CMplot 4.5+ (R), matplotlib 3.8+, pandas 2.2+, scipy 1.12+, plinkQC 0.3+. For locuszoom-style: locuszoomr 0.3+ (R) or pyranges + matplotlib.

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

# Manhattan, QQ, and Locuszoom Plots

**"Plot my GWAS results"** -> Render per-variant -log10(p) across the genome (Manhattan), compare expected vs observed p quantiles (QQ + λGC), overlay two traits with mirrored axes (Miami), and zoom into a locus with LD-colored points + recombination rate + gene track (locuszoom). The choices that matter: significance thresholds, axis truncation for ultra-significant peaks, lead-SNP labeling, and LD reference selection for regional plots.

- R: `qqman::manhattan` / `qqman::qq` (Turner 2018), `CMplot::CMplot`, `locuszoomr::locus_plot`
- Python: `matplotlib` + `pandas` for custom; `assocplots` for ready-made

## The Single Most Important Modern Insight -- The Threshold Is Always Conditional

The "genome-wide significant" line at `p < 5e-8` (Pe'er 2008 *Genet Epidemiol* 32:381) is calibrated for **European-ancestry common-variant GWAS** assuming ~1M effectively independent tests. It is the wrong threshold for:

- **Whole-genome sequencing** including rare variants (~5e-9 EUR, ~1e-9 AFR; Pulit 2017 *Genet Epidemiol* 41:145; Xu 2014 *Genet Epidemiol* 38:281)
- **Non-European ancestry** with different LD structure (typically more stringent)
- **TWAS / PWAS** with ~20,000 tested genes (Bonferroni 2.5e-6)
- **Multi-ethnic meta-analysis** (5e-9 by convention for trans-ancestry)
- **Burden / SKAT rare-variant tests** (per-gene; ~2.5e-6)
- **Locus-wise fine-mapping** (within-locus testing, no genome-wide correction needed)

A Manhattan plot's significance line is a contract with the reader about which multiple-testing regime applies. Mismatched thresholds over- or under-report hits.

## Decision Tree by Analysis

| Analysis | Genome-wide threshold | Suggestive threshold | Reference |
|----------|----------------------|---------------------|-----------|
| Common-variant GWAS (Eur) | 5e-8 | 1e-5 | Pe'er 2008 *Genet Epidemiol* 32:381 |
| Whole-genome sequencing (all variants, EUR) | 5e-9 | 5e-8 | Pulit 2017 *Genet Epidemiol* 41:145; Xu 2014 *Genet Epidemiol* 38:281 |
| Non-European ancestry (empirical per pop) | ~3.24e-8 AFR; ~9.26e-8 EAS | – | Kanai 2016 *J Hum Genet* 61:861 |
| TWAS (~20k genes) | 2.5e-6 (Bonferroni) | 1e-4 | Standard practice |
| PWAS (~5k proteins) | 1e-5 | 1e-4 | Standard practice |
| eQTL trans (genome-wide per probe) | Bonferroni over genes × variants | Per-tissue | GTEx convention |
| eQTL cis (within 1Mb) | nominal p < 1e-5 with permutation | – | GTEx FastQTL |
| Rare-variant gene burden | 2.5e-6 | 1e-4 | Bonferroni 20k genes |
| Trans-ancestry meta-analysis | 5e-9 | – | Convention |

## Genomic Inflation (λGC) -- The Mandatory QC Step

**Goal:** Quantify whether observed p-values are inflated relative to the chi-square null, indicating cryptic population structure, relatedness, or technical artifacts.

**Approach:** Convert observed p to chi-square; compute median chi-square divided by 0.4549 (the median of chi-square_1; Devlin-Roeder 1999); plot expected vs observed quantiles (QQ plot).

```r
library(qqman)
chisq <- qchisq(1 - df$P, df = 1)
lambda <- median(chisq) / 0.4549
# lambda = 1.0 -> no inflation
# lambda > 1.1 -> investigate; could indicate confounding
# lambda > 1.2 -> almost certainly confounded; principal components or LMM needed

qq(df$P, main = paste('QQ plot (lambda =', round(lambda, 3), ')'))
```

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

chisq = stats.chi2.isf(df['P'], df=1)
lambda_gc = np.median(chisq) / stats.chi2.ppf(0.5, df=1)

# QQ plot
expected = -np.log10(np.arange(1, len(df)+1) / (len(df) + 1))
observed = -np.log10(np.sort(df['P']))
fig, ax = plt.subplots(figsize=(4, 4))
ax.scatter(expected, observed, s=2)
ax.plot([0, max(expected)], [0, max(expected)], 'r--')
ax.set_xlabel(r'Expected $-\log_{10}(p)$')
ax.set_ylabel(r'Observed $-\log_{10}(p)$')
ax.set_title(f'QQ ($\\lambda_{{GC}} = {lambda_gc:.3f}$)')
```

**Interpretation**:
- λ = 1.00 ± 0.02 — well-calibrated
- λ > 1.05 — possible inflation; consider sample-size adjustment (`λ_1000 = 1 + (λ - 1) * 1000/n`)
- λ > 1.10 — confounded; population structure not removed; rerun with PC adjustment or LMM (BOLT-LMM, GEMMA, SAIGE)
- λ < 1.00 — deflation; usually a bug (wrong test statistic, conservative p-values)

Inflation can also be **legitimate polygenic signal** (Yang 2011 *Eur J Hum Genet* 19:807). Distinguish via LD-score regression intercept: confounding inflates intercept; polygenic signal inflates slope.

## Small-N and Rare-Variant Regimes -- When Standard Asymptotics Break

Standard logistic regression / Wald test p-values are anti-conservative when (a) case count < 200, (b) per-variant minor-allele count < 20, (c) case-control ratio is severely unbalanced (typical in EHR-derived cohorts). λGC may look normal but per-variant p-values are inflated independently — a Manhattan plot of these p-values is misleading regardless of inflation diagnostics.

| Regime | Test choice | Tool |
|--------|-------------|------|
| Balanced case-control, N>5000, MAF>0.01 | Standard logistic / linear regression | PLINK, REGENIE |
| Unbalanced (case fraction <10%), large N | SPA-corrected logistic regression | SAIGE, REGENIE Firth/SPA |
| Small N (<5000) | Penalized regression with bias correction | SAIGE Firth, REGENIE |
| Rare variants (MAC <20) | Gene-burden or SKAT-O | STAAR, REGENIE burden |

Specifically: **SAIGE** (Zhou 2018 *Nat Genet* 50:1335) and **REGENIE** (Mbatchou 2021 *Nat Genet* 53:1097) implement saddlepoint-approximation (SPA) and Firth-bias correction. Use them for any cohort with severe case-control imbalance; the Manhattan / QQ output is then defensibly calibrated.

## Manhattan Plot -- Canonical Layout

```r
library(qqman)
manhattan(df,
          chr = 'CHR', bp = 'BP', p = 'P', snp = 'SNP',
          col = c('#0072B2', '#56B4E9'),
          genomewideline = -log10(5e-8),
          suggestiveline = -log10(1e-5),
          ylim = c(0, max(-log10(df$P)) * 1.1),
          highlight = lead_snps,
          annotatePval = 5e-8,
          annotateTop = TRUE)
```

```python
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

def manhattan_plot(df, p_threshold=5e-8, suggestive=1e-5, y_cap=None):
    df = df.sort_values(['CHR', 'BP']).copy()
    df['neg_log10_p'] = -np.log10(df['P'])
    if y_cap:
        df['neg_log10_p'] = df['neg_log10_p'].clip(upper=y_cap)

    df['x'] = np.arange(len(df))
    chr_ticks = df.groupby('CHR')['x'].median()
    chr_colors = ['#0072B2', '#56B4E9']

    fig, ax = plt.subplots(figsize=(10, 4))
    for i, (chrom, group) in enumerate(df.groupby('CHR')):
        ax.scatter(group['x'], group['neg_log10_p'],
                   c=chr_colors[i % 2], s=3, rasterized=True)
    ax.axhline(-np.log10(p_threshold), color='red', linestyle='--', lw=0.5)
    ax.axhline(-np.log10(suggestive), color='grey', linestyle='--', lw=0.5)
    ax.set_xticks(chr_ticks)
    ax.set_xticklabels(chr_ticks.index, rotation=0)
    ax.set_xlabel('Chromosome')
    ax.set_ylabel(r'$-\log_{10}(p)$')
    return fig
```

## Extreme Tail Handling -- The Y-Axis Cap

Genome-wide significant peaks routinely reach -log10(p) = 100+ (e.g., GWAS of BMI at FTO). The visual effect: one peak fills the y-axis, all other signal is invisible.

**Fixes (ordered by preference):**

1. **Cap and indicate:** `y_cap = 25`; clip points above the cap; mark capped points with `^` arrow at the top o