PISA squid plot

Brush input and output positions to explore base-to-base effects around the Drosophila sog enhancer, with accessibility, contribution, and motif tracks.

Each diagonal link joins an influential input base to an affected ATAC-seq output position around the Drosophila sog enhancer, following Figure 2c of McAnany et al. Accessibility, contribution scores, and motif intervals provide context on the same genomic scale.

What to notice

Positive effects are red and negative effects are blue; opacity shows effect magnitude. Drag on Accessibility to select output positions, or on Contribution score to select input positions. Links matching both active brushes retain their effect encoding and draw above muted links. An empty brush leaves its endpoint unconstrained. Bars and sequence letters outside an active brush turn gray.

Hold Shift and hover to highlight a link, including one outside the brushes. With no active brush or Shift-hover, all links look normal. Double-click a track to clear its brush. Pan and zoom to inspect the locus; zooming in reveals base-colored letters in the contribution track.

How the chart is built

Python authors the chart and references the prepared remote Parquet tables; it does not load or recalculate PISA effects. GenomeSpy loads the tables and executes the declarative filters and formulas in the browser. These select the margin tracks, calculate absolute effect magnitude, and format motif labels.

The link view defines one highlightedLink predicate combining Shift-hover with interval brushes projected onto x (input) and x2 (output). gs.when({"ref": "highlightedLink"}) reuses it for color, opacity, and draw order. The + and & operators layer the motif annotations and stack the tracks.

See the PISA interaction matrix for the full effect field and the official GenomeSpy squid example for the original JSON specification.

Data use and provenance

This example uses a Parquet extract of the dm6 sog locus from the supporting data for McAnany et al., Positional interpretation of cis-regulatory code and nucleosome organization with deep learning models, prepared with the GenomeSpy recipe. The extract is distributed under GPL-2.0-or-later; the accessibility model’s training data are GEO accession GSE218852.

Code

"""PISA squid plot.

Brush input and output positions to explore base-to-base effects around the
Drosophila sog enhancer, with accessibility, contribution, and motif tracks.
"""

import genome_spy as gs

DATA = "https://data.genomespy.app/datasets/bpreveal-pisa/v3/"
tracks = gs.Data(url=DATA + "fig2cd-atac-tracks.parquet", format={"type": "parquet"})
MUTED = "#d8dbe0"

# Each track updates a selection stored in the parent view's parameter scope.
output_brush = gs.selection_interval(
    "accessibilityRegion",
    push="outer",
    persist=False,
    encodings=["x"],
    extent="view",
    on="mousedown",
    clear="dblclick",
    zoom=False,
)
input_brush = gs.selection_interval(
    "contributionRegion",
    push="outer",
    persist=False,
    encodings=["x"],
    extent="view",
    on="mousedown",
    clear="dblclick",
    zoom=False,
)
hover = gs.selection_point(
    "pisaLinkHover",
    push="outer",
    persist=False,
    on="pointerover[event.shiftKey]",
    clear="mousemove[!event.shiftKey]",
)

accessibility = (
    gs.Chart(tracks)
    .transform_filter(gs.datum.track == "prediction")
    .mark_rect(minOpacity=1)
    .encode(
        x=gs.X("position:I").title("dm6 chrX"),
        y=gs.Y("value:Q").scale(zero=True).axis(None),
        y2=gs.datum(0),
        color=gs.when(output_brush)
        .then(gs.value("#332288"))
        .otherwise(gs.value(MUTED)),
        tooltip=[
            gs.Tooltip("position").title("Output position"),
            gs.Tooltip("value").title("Prediction").format(".4g"),
        ],
    )
    .properties(
        name="prediction",
        title=gs.Title(text="Accessibility", style="overlay-title"),
        height={"grow": 0.16},
        cursor="text",
    )
    .add_params(output_brush)
)

# A projected brush tests the input (x) or output (x2) endpoint of a link.
# The final OR makes Shift-hover isolate links when both brushes are empty.
highlighted_link = {
    "or": [
        {"param": "pisaLinkHover", "empty": False},
        {
            "and": [
                {"param": "accessibilityRegion", "project": {"x": "x2"}},
                {"param": "contributionRegion", "project": {"x": "x"}},
                {
                    "or": [
                        {"param": "pisaLinkHover"},
                        {
                            "param": "accessibilityRegion",
                            "project": {"x": "x2"},
                            "empty": False,
                        },
                        {
                            "param": "contributionRegion",
                            "project": {"x": "x"},
                            "empty": False,
                        },
                    ]
                },
            ]
        },
    ]
}
highlight = gs.when({"ref": "highlightedLink"})
effect_color = (
    gs.Color("effect:Q")
    .scale(
        domain=[
            -0.216404,
            -0.173123,
            -0.129843,
            -0.086562,
            -0.043281,
            0,
            0.043281,
            0.086562,
            0.129843,
            0.173123,
            0.216404,
        ],
        range=[
            "#053061",
            "#2166ac",
            "#4393c3",
            "#92c5de",
            "#d1e5f0",
            "#ffffff",
            "#fddbc7",
            "#f4a582",
            "#d6604d",
            "#b2182b",
            "#67001f",
        ],
        clamp=True,
    )
    .legend(
        title="PISA (log2(fc))",
        orient="left",
        direction="vertical",
        gradientLength=84,
        gradientThickness=18,
        gradientStrokeColor="#777777",
        gradientStrokeWidth=0.5,
        values=[-0.2, 0, 0.2],
    )
)
links = (
    gs.Chart(gs.Data(url=DATA + "fig2c-atac-links.parquet", format={"type": "parquet"}))
    .transform_calculate(absEffect=gs.expr.abs(gs.datum.effect))
    .mark_link(linkShape="diagonal", orient="vertical", size=1.5, minPickingSize=2)
    .encode(
        x=gs.X("source:I").title("dm6 chrX").buildIndex(False),
        x2=gs.X2("target"),
        y=gs.datum(0, type="quantitative").scale(domain=[0, 1]).axis(None),
        y2=gs.datum(1),
        order=highlight.then(gs.value(1)).otherwise(gs.value(0)),
        color=highlight.then(effect_color).otherwise(gs.value(MUTED)),
        opacity=highlight.then(
            gs.Opacity("absEffect:Q")
            .scale(
                domain=[0, 0.216404],
                range=[0, 1],
                clamp=True,
            )
            .legend(None)
        ).otherwise(gs.value(0.1)),
        tooltip=[
            gs.Tooltip("source").title("Input position"),
            gs.Tooltip("target").title("Output position"),
            gs.Tooltip("effect").title("PISA effect (log2 fold change)").format(".4f"),
        ],
    )
    .properties(name="pisa-links", predicates={"highlightedLink": highlighted_link})
    .add_params(hover)
)

# Overlay motif intervals and labels near the bottom of the link panel.
motif_blocks = (
    gs.Chart()
    .mark_rect(minOpacity=1)
    .encode(
        color=gs.Color("motifLabel:N")
        .scale(
            domain=["M1bp", "Gaga", "Zelda"],
            range=["#bbcc33", "#44bb99", "#99ddff"],
        )
        .legend(orient="left", direction="vertical", symbolOpacity=1),
        tooltip=[
            gs.Tooltip("motifLabel").title("Motif"),
            gs.Tooltip("start").title("Start"),
            gs.Tooltip("end").title("End"),
            gs.Tooltip("strand").title("Strand"),
            gs.Tooltip("score").title("Score"),
        ],
    )
)
motif_labels = (
    gs.Chart()
    .mark_text(
        align="center",
        baseline="middle",
        paddingX=3,
        tooltip=None,
    )
    .encode(text="motifLabel", color=gs.value("black"))
)
motifs = (
    (motif_blocks + motif_labels)
    .properties(
        name="motifs",
        data=gs.Data(
            url=DATA + "fig2cd-atac-motifs.parquet",
            format={"type": "parquet"},
        ),
    )
    .transform_calculate(
        motifLabel="datum.name == 'm1bp' ? 'M1bp' : datum.name == 'gaga' ? 'Gaga' : 'Zelda'"
    )
    .encode(
        x=gs.X("start:I").title("dm6 chrX"),
        x2="end",
        y=gs.value(gs.expr("4 / height")),
        y2=gs.value(gs.expr("20 / height")),
    )
)

# Switch from compact bars to base-colored contribution letters as we zoom in.
bars = (
    gs.Chart()
    .mark_rect(minOpacity=1)
    .encode(
        color=gs.when(input_brush).then(gs.value("#332288")).otherwise(gs.value(MUTED)),
    )
)
logo = (
    gs.Chart()
    .mark_text(
        font="Source Sans Pro",
        fontWeight=700,
        size=100,
        squeeze=True,
        fitToBand=True,
        paddingX=0,
        paddingY=0,
        logoLetters=True,
    )
    .encode(
        text="base",
        color=gs.when(input_brush)
        .then(
            gs.Color("base:N")
            .scale(
                type="ordinal",
                domain=["A", "C", "G", "T", "N"],
                range=["#009E73", "#0072B2", "#F0E442", "#D55E00", "#BDBDBD"],
            )
            .legend(None)
        )
        .otherwise(gs.value(MUTED)),
    )
)
contribution = (
    gs.multiscale(
        bars,
        logo,
        stops={
            "channel": "x",
            "values": [0.15],
            "transition": {"type": "lerp", "halfLife": 60},
        },
    )
    .properties(
        name="importance",
        data=tracks,
        height={"grow": 0.16},
        cursor="text",
        title=gs.Title(text="Contribution score", style="overlay-title"),
    )
    .transform_filter(gs.datum.track == "importance")
    .encode(
        x=gs.X("position:I").title("dm6 chrX"),
        y=gs.datum(0, type="quantitative").scale(zero=True).axis(None),
        y2="value",
        tooltip=[
            gs.Tooltip("position").title("Input position"),
            gs.Tooltip("base").title("Base"),
            gs.Tooltip("value").title("Contribution").format(".4g"),
        ],
    )
    .add_params(input_brush)
)

chart = (
    (accessibility & (links + motifs).resolve_scale(color="independent") & contribution)
    .properties(
        spacing=0,
        padding={"top": 8, "right": 30, "bottom": 8, "left": 10},
        scales=gs.scales(x=gs.Scale(domain=[15646649, 15647250], zoom=True)),
    )
    .add_params(
        gs.param("accessibilityRegion", value=None),
        gs.param("contributionRegion", value=None),
        gs.param("pisaLinkHover", value=None),
    )
    .resolve_scale(x="shared")
    .resolve_axis(x="shared")
    .resolve_legend(color="collected")
    .configure_axis(domain=False)
    .configure_legend(
        labelFontSize=11,
        titleFontSize=11,
        layout={"left": {"anchor": "middle", "wrap": False}},
    )
    .configure_legend_track(style=None)
    .configure_title(fontSize=12, fontWeight="normal", offset=2)
    .configure_view(stroke="transparent", strokeWidth=0)
)

Show generated specification