Manhattan plot
Genomics — GWAS summary statistics
Example from the compendium of canonical charts
Python Code
"""Genomics — Manhattan plot (GWAS summary statistics)."""
import pandas as pd
import numpy as np
import plotly.graph_objects as go
URL = "https://raw.githubusercontent.com/Cloufield/GWASTutorial/main/02_Linux_basics/sumstats.txt"
# Alternating palette for chromosomes
CHR_COLORS = [VIOLET, TEAL, PINK, GREEN, PINK, VIOLET, TEAL, PINK, GREEN, PINK,
VIOLET, TEAL, PINK, GREEN, PINK, VIOLET, TEAL, PINK, GREEN, PINK,
VIOLET, TEAL, PINK]
# GRCh38 approximate chromosome sizes (bp) — used for consistent x-axis spacing
CHR_SIZES = {
"1": 248956422, "2": 242193529, "3": 198295559, "4": 190214555,
"5": 181538259, "6": 170805979, "7": 159345973, "8": 145138636,
"9": 138394717, "10": 133797422, "11": 135086622, "12": 133275309,
"13": 114364328, "14": 107043718, "15": 101991189, "16": 90338345,
"17": 83257441, "18": 80373285, "19": 58617616, "20": 64444167,
"21": 46709983, "22": 50818468,
}
def generate():
print("fetching GWAS summary statistics …")
try:
df = fetch_csv(URL, sep="\t")
print(f" {len(df):,} rows, cols: {list(df.columns)}")
# Expected: #CHROM, POS, P
chrom_col = next((c for c in df.columns if "chrom" in c.upper()), df.columns[0])
pos_col = next((c for c in df.columns if c.upper() in ("POS", "BP", "POSITION")), df.columns[1])
p_col = next((c for c in df.columns if c.upper() in ("P", "PVAL", "P_VALUE", "PVALUE")), df.columns[2])
df = df[[chrom_col, pos_col, p_col]].copy()
df.columns = ["chrom", "pos", "p"]
df["chrom"] = (df["chrom"].astype(str)
.str.replace("#", "")
.str.replace("chr", "", case=False)
.str.strip())
df["pos"] = pd.to_numeric(df["pos"], errors="coerce")
df["p"] = pd.to_numeric(df["p"], errors="coerce")
df = df.dropna()
# keep standard chromosomes
valid_chroms = [str(i) for i in range(1, 23)] + ["X", "Y"]
df = df[df["chrom"].isin(valid_chroms)].copy()
# sort by chrom numerically
df["chrom_n"] = pd.Categorical(
df["chrom"],
categories=valid_chroms,
ordered=True
)
df = df.sort_values(["chrom_n", "pos"]).reset_index(drop=True)
# Require at least 11 chromosomes for a useful plot; fall back to synthetic otherwise
n_chroms = df["chrom"].nunique()
print(f" {n_chroms} chromosomes present")
if n_chroms < 11:
print(f" only {n_chroms} chromosomes — insufficient coverage, using synthetic")
raise ValueError(f"only {n_chroms} chromosomes in real data")
used_synthetic = False
except Exception as e:
print(f" fetch failed ({e}); generating synthetic GWAS data")
used_synthetic = True
rng = np.random.default_rng(42)
rows = []
for ch in range(1, 23):
n = rng.integers(3000, 6000)
pos = np.sort(rng.integers(1, 250_000_000, n))
p = rng.uniform(0, 1, n)
# add some significant hits on chrs 3, 7, 11
if ch in (3, 7, 11):
peak_pos = rng.integers(50_000_000, 200_000_000, 3)
for pp in peak_pos:
nearby = np.abs(pos - pp) < 500_000
p[nearby] = rng.uniform(1e-10, 5e-8, nearby.sum())
rows.append(pd.DataFrame({"chrom": str(ch), "pos": pos, "p": p}))
df = pd.concat(rows, ignore_index=True)
valid_chroms = [str(i) for i in range(1, 23)]
df["chrom_n"] = pd.Categorical(df["chrom"], categories=valid_chroms, ordered=True)
df = df.sort_values(["chrom_n", "pos"]).reset_index(drop=True)
# Build cumulative position using standard chromosome sizes for even spacing
valid_chroms_auto = [str(c) for c in df["chrom_n"].cat.categories.tolist()]
offset_map = {}
cumulative = 0
for ch in valid_chroms_auto:
offset_map[ch] = int(cumulative)
size = int(CHR_SIZES.get(ch, 100_000_000))
cumulative = int(cumulative) + size + 5_000_000
df["cum_pos"] = df["pos"].astype(int) + df["chrom"].map(offset_map)
# -log10(p)
df["neg_log_p"] = -np.log10(np.clip(df["p"], 1e-300, 1))
# Chromosome midpoints for x-axis ticks
tick_pos = df.groupby("chrom_n", observed=True)["cum_pos"].mean()
tick_labels = tick_pos.index.tolist()
GWAS_THRESHOLD = -np.log10(5e-8) # ~7.3
# Subsample non-significant points to keep JSON manageable
sig_mask = df["neg_log_p"] >= GWAS_THRESHOLD
ns_max = 8000
ns_df = df[~sig_mask]
if len(ns_df) > ns_max:
ns_df = ns_df.sample(ns_max, random_state=42)
df_plot = pd.concat([ns_df, df[sig_mask]]).sort_values(["chrom_n", "pos"])
# Build traces per chromosome (use subsampled df_plot for non-significant)
traces = []
chroms_sorted = df["chrom_n"].cat.categories.tolist()
for i, ch in enumerate(chroms_sorted):
sub = df_plot[df_plot["chrom_n"] == ch]
if len(sub) == 0:
continue
color = CHR_COLORS[i % len(CHR_COLORS)]
traces.append(go.Scattergl(
x=sub["cum_pos"],
y=sub["neg_log_p"],
mode="markers",
marker=dict(size=3, opacity=0.7, color=color),
name=f"Chr {ch}",
showlegend=False,
hovertemplate=f"Chr {ch}<br>pos=%{{customdata}}<br>-log10(p)=%{{y:.2f}}<extra></extra>",
customdata=sub["pos"],
))
# Significance threshold line
x_max = df["cum_pos"].max()
traces.append(go.Scatter(
x=[0, x_max],
y=[GWAS_THRESHOLD, GWAS_THRESHOLD],
mode="lines",
line=dict(color="red", dash="dash", width=1.2),
name="p = 5×10⁻⁸",
hoverinfo="skip",
))
# Annotate peaks above threshold
sig = df[df["neg_log_p"] > GWAS_THRESHOLD]
annotations = []
if len(sig) > 0:
# cluster peaks by chromosome, take top per chrom
peak_per_chr = sig.loc[sig.groupby("chrom_n", observed=True)["neg_log_p"].idxmax()]
for _, row in peak_per_chr.head(6).iterrows():
annotations.append(dict(
x=row["cum_pos"],
y=row["neg_log_p"],
text=f"Chr {row['chrom_n']}",
showarrow=True,
arrowhead=2,
arrowsize=0.8,
ax=0, ay=-25,
font=dict(size=10),
))
fig = go.Figure(traces)
fig.update_layout(
xaxis=dict(
title="Chromosome",
tickvals=[tick_pos[ch] for ch in tick_labels if ch in tick_pos],
ticktext=[str(ch) for ch in tick_labels if ch in tick_pos],
tickangle=0,
showgrid=False,
),
yaxis=dict(title="−log₁₀(p)", zeroline=False),
annotations=annotations,
legend=dict(orientation="h", y=1.05),
margin=dict(t=40, b=60, l=70, r=40),
height=520,
)
return fig
if __name__ == "__main__":
generate()Made with Plotly