Manhattan plot

Genome-wide association hits on a locus-aware chromosome axis. Alternating chromosome colors separate the blocks, dashed rules mark significance thresholds, and the strongest peaks are outlined for emphasis.

Move the −log10 p cutoff slider to change the red significance line. Points above the line turn red immediately; a higher cutoff is stricter. GenomeSpy updates the rule and highlighted set in the browser.

Data use and provenance

The packaged table is a subset of HapMap variants from the MIT-licensed manhattanly and Plotly datasets projects. Genomic coordinates and annotations are real, but the association statistics are simulated and are not biological findings. During data loading, the package removes invalid p-values and calculates -log10(p). GenomeSpy then renders the prepared table and applies the interactive threshold in the browser.

Code

"""Manhattan plot.

Genome-wide association hits on a locus-aware chromosome axis. Alternating
chromosome colors separate the blocks, dashed rules mark significance
thresholds, and the strongest peaks are outlined for emphasis.
"""

import genome_spy as gs
from genome_spy.datasets._hapmap import hapmap_manhattan_data
from genome_spy.schema import GenomeAxis, Scale


GENOME_WIDE_P = 5e-8
SUGGESTIVE_P = 1e-5


# Load association results and the starting significance thresholds.
data, _top_hits, domains = hapmap_manhattan_data(
    genome_wide_p=GENOME_WIDE_P,
    suggestive_p=SUGGESTIVE_P,
)
# Moving this slider redraws the red guide and outlines the hits above it.
significance_cutoff = gs.param(
    "manhattanSignificanceCutoff",
    value=round(domains["genome_wide_y"], 1),
    bind=gs.binding_range(
        min=3,
        max=domains["y_domain"][1],
        step=0.1,
        name="−log10 p cutoff: ",
    ),
)

# Alternate colors to make neighboring chromosomes easier to distinguish.
chrom_colors = Scale().range(["#5b8fd6", "#8f98a3"])

# Label chromosomes and add faint backgrounds to separate them.
axis = (
    GenomeAxis()
    .title("Genomic position")
    .chromGrid(True)
    .chromGridOpacity(0.14)
    .chromGridFillEven("#f4f6fb")
    .chromGridFillOdd("#ffffff")
    .chromLabels(True)
    .chromLabelFontSize(11)
    .chromTicks(True)
    .chromTickSize(10)
    .labelFontSize(10)
    .grid(False)
)

# Draw one point per variant at its genomic position.
points = (
    gs.Chart()
    .mark_point(size=20, filled=True, opacity=0.82)
    .encode(
        x=gs.Locus("chrom", "BP").scale(assembly="hg18").axis(axis),
        y=gs.Y("neglog:Q")
        .scale(reverse=False, domain=domains["y_domain"])
        .title("−log10 p"),
        color=gs.Color("chrom_group:N")
        .scale(chrom_colors)
        .legend(title="Chromosome group"),
    )
)

# Move the red significance line with the slider.
genome_wide_rule = (
    gs.Chart([{}])
    .mark_rule(strokeDash=[6, 4], size=1.4, color="#c53b2c")
    .encode(
        y=gs.Y(gs.datum(significance_cutoff), type="quantitative")
        .scale(reverse=False, domain=domains["y_domain"])
        .title("−log10 p")
    )
)

# Keep a second, less strict threshold visible for comparison.
suggestive_rule = (
    gs.Chart([{"threshold": domains["suggestive_y"]}])
    .mark_rule(strokeDash=[2, 4], size=1.2, color="#d48b31")
    .encode(
        y=gs.Y("threshold:Q")
        .scale(reverse=False, domain=domains["y_domain"])
        .title("−log10 p")
    )
)

# Highlight variants above the slider's current threshold.
highlight_points = (
    gs.Chart()
    .transform_collect()
    .transform_filter(gs.datum.neglog >= significance_cutoff)
    .mark_point(
        size=48,
        filled=True,
        color="#c53b2c",
        stroke="black",
        strokeWidth=0.5,
    )
    .encode(
        x=gs.Locus("chrom", "BP").scale(assembly="hg18"),
        y=gs.Y("neglog:Q")
        .scale(reverse=False, domain=domains["y_domain"])
        .title("−log10 p"),
    )
)

# Combine the points and guide lines, then attach the slider.
association_track = genome_wide_rule + suggestive_rule + points + highlight_points

chart = association_track.properties(
    assembly="hg18",
    title="HapMap genome-wide association scan",
    data=data,
    description=(
        "A Manhattan plot with a locus-aware chromosome axis and an "
        "interactive significance threshold."
    ),
).add_params(significance_cutoff)

Show specification