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 distributed through the 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 generated specification