Bio Data Visualization Manhattan Qq Locuszoom

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

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).

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), ')'))

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.

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

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)

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 of the panel. Use Y-axis label "−log10(P), capped at 25"
2. Split y-axis via ggbreak::scale_y_break() (R) or axes_grid1.divider (matplotlib)
3. Two-panel plot with full y range in top, zoomed range in bottom
4. Use sqrt or asinh transform: less intuitive but preserves all data; rarely chosen for Manhattan

Miami Plot -- Two-Trait Comparison

# CMplot supports Miami natively
library(CMplot)
CMplot(list(trait1_df, trait2_df),
       plot.type = 'm',
       multraits = TRUE,
       threshold = 5e-8,
       threshold.col = 'red',
       col = list(c('#0072B2','#56B4E9'), c('#D55E00','#E69F00')),
       file = 'jpg', file.output = TRUE)

Miami plot mirrors trait 1 above the x-axis, trait 2 below. Useful for shared-locus discovery (mirrored peaks at the same locus = pleiotropy candidate).

Locuszoom-Style Regional Plot

A locuszoom plot zooms into a ~1Mb window around a lead SNP, colors SNPs by LD r² to the lead, overlays recombination rate (cM/Mb), and shows the gene track. The canonical tool is locuszoom.org (Pruim 2010 Bioinformatics 26:2336); the R package is locuszoomr (Lai 2024).

library(locuszoomr)
loc <- locus(gene = 'TCF7L2',
             flank = 5e5,
             ens_db = 'EnsDb.Hsapiens.v86',
             data = gwas_df,
             snp = 'SNP', chrom = 'CHR', pos = 'BP', p = 'P', labs = 'SNP')
# LD computed via LDlinkR / 1000G reference
loc <- link_LD(loc, pop = 'EUR', token = ldlink_token)
locus_plot(loc, labels = c('index', 'top'))

LD reference choice: ALWAYS match the GWAS population. Using a 1000G European LD reference for a Japanese GWAS produces wrong LD colorings and misleads fine-mapping.

Per-Method Failure Modes

Inflation diagnosed as polygenic signal

Trigger: λGC = 1.15; analyst concludes "polygenic," moves on.

Mechanism: Inflation can be confounding (population structure, relatedness, technical) OR polygenicity. LD-score regression separates them: intercept = confounding; slope = polygenicity.

Symptom: Top hits replicate poorly in independent cohorts.

Fix: Run LDSC (ldsc.py --h2); intercept significantly > 1 indicates confounding. Adjust with 10-20 PCs or switch to LMM.

Wrong significance threshold for the analysis

Trigger: Plotting Bonferroni-naive 5e-8 line on a TWAS or rare-variant burden plot.

Mechanism: 5e-8 is calibrated for ~1M independent common-variant tests; other analyses have different effective test counts.

Symptom: Reviewer flags "this gene doesn't pass Bonferroni" but the plot's red line is at 5e-8.

Fix: Match the threshold to the analysis (table above).

Cap with no indication

Trigger: Y-axis clipped at 25 without arrow markers for capped points.

Mechanism: Reader cannot tell how high the true peak is.

Symptom: Reviewer asks for actual p-value; the reported 1e-100 conflicts with the displayed cap at 25.

Fix: Mark capped points with ^ symbol; annotate axis "(capped at 25)"; include unclipped numerical value in caption.