Volcano plot¶
Effect size against −log10 p-value, with color highlighting points that clear both an effect-size and a significance cutoff.
Use the two sliders to define how large an effect must be and how strong its p-value evidence must be. The guide lines move with the sliders, and points that pass both cutoffs are recolored 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. The association p-values and effect sizes are simulated, so the plot
is for visualization demonstration only. During data loading, the package
calculates -log10(p). GenomeSpy then renders the prepared table and classifies
the points again whenever a threshold slider moves.
Code¶
"""Volcano plot.
Effect size against −log10 p-value, with color highlighting points that clear
both an effect-size and a significance cutoff.
"""
import genome_spy as gs
from genome_spy.datasets._hapmap import hapmap_volcano_data
from genome_spy.schema import Scale
# Load association results and starting axis limits and cutoffs.
data, domains = hapmap_volcano_data()
X_DOMAIN = domains["x_domain"]
Y_DOMAIN = domains["y_domain"]
EFFECT_CUTOFF = domains["effect_cutoff"]
NEGLOG_P_CUTOFF = domains["neglog_pvalue_cutoff"]
# Grow points gently when zooming in, without letting them become too large.
ZOOM_LEVEL = gs.Expression("zoomLevel")
POINT_SIZE = gs.expr(gs.expr.min(16 * gs.expr.pow(ZOOM_LEVEL, 0.75), 64))
# Add sliders that update the guide lines and point colors together.
effect_cutoff = gs.param(
"hapmapEffectCutoff",
value=EFFECT_CUTOFF,
bind=gs.binding_range(
min=0,
max=X_DOMAIN[1],
step=0.1,
name="Absolute effect cutoff: ",
),
)
significance_cutoff = gs.param(
"hapmapSignificanceCutoff",
value=NEGLOG_P_CUTOFF,
bind=gs.binding_range(
min=0,
max=Y_DOMAIN[1],
step=0.25,
name="−log10 p cutoff: ",
),
)
# Classify points as protective or risk only when they pass both cutoffs.
ASSOCIATION_EXPRESSION = gs.expr.if_(
(gs.datum.neglog >= significance_cutoff)
& (gs.expr.abs(gs.datum.EFFECTSIZE) >= effect_cutoff),
gs.expr.if_(gs.datum.EFFECTSIZE < 0, "protective", "risk"),
"n.s.",
)
# Use blue for protective effects, red for risk, and grey for the rest.
association_colors = Scale(
domain=["protective", "n.s.", "risk"],
range=["#3e8cb6", "#c9d1d9", "#c53b2c"],
)
# Plot effect size against significance, coloring points by the current cutoffs.
points = (
gs.Chart()
.transform_collect()
.transform_formula(expr=ASSOCIATION_EXPRESSION, as_="association")
.mark_point(size=POINT_SIZE, filled=True, opacity=0.6)
.encode(
x=gs.X("EFFECTSIZE:Q")
.scale(domain=X_DOMAIN, zoom=True)
.title("Effect size (beta)"),
y=gs.Y("neglog:Q")
.scale(reverse=False, domain=Y_DOMAIN, zoom=True)
.title("−log10 p"),
color=gs.Color("association:N")
.scale(association_colors)
.legend(title="Association"),
)
)
# Mark the effect-size cutoff on both sides of zero.
effect_cutoffs = (
gs.Chart([{"side": -1}, {"side": 1}])
.transform_collect()
.transform_formula(expr=gs.datum.side * effect_cutoff, as_="x")
.mark_rule(strokeDash=[4, 4], size=1, color="#8f98a3")
.encode(x=gs.X("x:Q").scale(domain=X_DOMAIN, zoom=True).title("Effect size (beta)"))
)
# Move the horizontal significance line with its slider.
significance_rule = (
gs.Chart([{}])
.transform_collect()
.transform_formula(expr=significance_cutoff, as_="y")
.mark_rule(strokeDash=[4, 4], size=1, color="#8f98a3")
.encode(
y=gs.Y("y:Q").scale(reverse=False, domain=Y_DOMAIN, zoom=True).title("−log10 p")
)
)
# Combine the guide lines and points, then attach the sliders.
chart = (
(effect_cutoffs + significance_rule + points)
.properties(
title="HapMap association volcano",
data=data,
description=(
"Effect size versus significance, with interactive effect-size and "
"p-value cutoffs controlling the guide lines and point colors."
),
)
.add_params(effect_cutoff, significance_cutoff)
)