Dynseq binding-QTL tracks

Sequence-logo tracks compare reference and alternate SPI1 contribution scores over the same base-resolution locus.

Data use and provenance

The visualization loads the original chip_imp_ref.bw and chip_imp_alt.bw directly from the pinned dynseq-paper source revision. The score values and missing-value gaps are not modified. The Zenodo source data are distributed under CC BY 4.0.

Source attribution: dynseq tracks data, Zenodo record 6582100; Nair et al., The dynseq browser track shows context-specific features at nucleotide resolution, Nature Genetics 54, 1581–1583 (2022); and Tehranchi et al., Pooled ChIP-Seq Links Variation in Transcription Factor Binding to Complex Disease Risk, Cell 165, 730–741 (2016).

What to notice

Reference and alternate SPI1 contribution scores use the same genomic scale. At base resolution, signed sequence-logo letters show the reference bases and make the direction and magnitude of each contribution visible.

Python implementation

A named template defines the shared baseline, BigWig lookup, and indexed-FASTA sequence-logo layers. Two imports instantiate it with reference and alternate allele parameters, which select the corresponding BigWig source and title. Python authors the flatten, formula, coordinate-lookup, and filter transforms; GenomeSpy executes them in the browser as the lazy genomic window changes.

See the official GenomeSpy example for the biological context and original grammar explanation.

Code

"""Dynseq binding-QTL tracks.

Sequence-logo tracks compare reference and alternate SPI1 contribution scores
over the same base-resolution locus.
"""

import genome_spy as gs

# Choose which allele's scores and title to use for a track.
allele = gs.param("allele", value="ref")
REF_BIGWIG_URL = (
    "https://raw.githubusercontent.com/kundajelab/dynseq-paper/"
    "febc9180d72e92302d35c549002e0d56c79c536e/SPI1_bQTL/"
    "bigwigs/chip_imp_ref.bw"
)
ALT_BIGWIG_URL = (
    "https://raw.githubusercontent.com/kundajelab/dynseq-paper/"
    "febc9180d72e92302d35c549002e0d56c79c536e/SPI1_bQTL/"
    "bigwigs/chip_imp_alt.bw"
)

# Build a reusable track: each letter's height shows its contribution score.
allele_track = (
    gs.layer(
        # Draw a line at zero to separate positive and negative scores.
        gs.Chart(
            [{}],
        )
        .mark_rule()
        .encode(
            y=gs.datum(0, type="quantitative"),
            color=gs.value("gray"),
        )
        .properties(name="baseline"),
        # Stretch each colored DNA letter from zero to its score.
        gs.Chart()
        .mark_text(
            font="Source Sans Pro",
            fontWeight=700,
            size=100,
            squeeze=True,
            fitToBand=True,
            paddingX=0,
            paddingY=0,
            logoLetters=True,
        )
        .encode(
            x=gs.Locus("chrom", "pos"),
            y=gs.datum(0, type="quantitative")
            .scale(zero=True, nice=False, reverse=False)
            .axis(title="Score"),
            y2=gs.Y2("score"),
            text=gs.Text("base"),
            color=gs.Color("base:N")
            .scale(
                domain=["A", "C", "G", "T"],
                range=["green", "blue", "orange", "red"],
            )
            .legend(None),
            tooltip=[
                gs.Tooltip("base:N"),
                gs.Tooltip("score:Q"),
            ],
        )
        .properties(name="dynseq"),
    )
    .properties(
        title=gs.Title(
            text=gs.expr(
                gs.expr.if_(
                    allele == "ref",
                    "Reference allele (C)",
                    "Alternate allele (G)",
                )
            ),
            style="overlay-title",
        ),
        height=120,
    )
    .add_params(allele)
    # Match each base to its score in the chosen allele's BigWig file.
    .transform_coordinate_lookup(
        from_={
            "data": gs.lazy.bigwig(
                gs.expr(gs.expr.if_(allele == "ref", REF_BIGWIG_URL, ALT_BIGWIG_URL)),
                pixelsPerBin=1,
            ),
            "transform": [
                {
                    "type": "formula",
                    "expr": gs.datum.start,
                    "as": "pos",
                }
            ],
        },
        key=["chrom", "pos"],
        values=["score"],
    )
    # Skip bases without scores.
    .transform_filter(gs.expr.isValid(gs.datum.score))
    # Show G instead of C at rs5764238 in the alternate-allele track.
    .transform_formula(
        expr=gs.expr.if_(
            (allele == "alt") & (gs.datum.pos == 43720929),
            "G",
            gs.datum.base,
        ),
        as_="base",
        description=(
            "Show the rs5764238 alternate allele while retaining the shared "
            "reference FASTA source."
        ),
    )
)

# Use the same track twice to compare the reference and alternate alleles.
chart = (
    gs.vconcat(
        gs.import_view(template="allele-track", params={"allele": "ref"}),
        gs.import_view(template="allele-track", params={"allele": "alt"}),
    )
    .properties(
        assembly="hg38",
        data=gs.lazy.indexed_fasta("https://data.genomespy.app/genomes/hg38/hg38.fa"),
        scales=gs.scales(
            x=gs.Scale(
                domain=[
                    {"chrom": "chr22", "pos": 43720850},
                    {"chrom": "chr22", "pos": 43720960},
                ],
                zoom={
                    "extent": [
                        {"chrom": "chr22", "pos": 43719872},
                        {"chrom": "chr22", "pos": 43721985},
                    ]
                },
            )
        ),
        templates={"allele-track": allele_track},
        description=(
            "Reference and alternate SPI1 contribution logos for rs5764238. "
            "Data source: https://github.com/kundajelab/dynseq-paper"
        ),
    )
    # Split the loaded sequence into uppercase letters with genomic positions.
    .transform_flatten_sequence(field="sequence", as_=["rawPos", "base"])
    .transform_formula(expr=gs.expr.upper(gs.datum.base), as_="base")
    .transform_formula(expr=gs.datum.start + gs.datum.rawPos, as_="pos")
    .resolve_scale(y="shared")
    .resolve_axis(x="shared")
    .configure_view(fill="#FAFAFA")
)

Show specification