Rainfall plot

Inter-mutation distance across the genome for a single tumor sample. Point color encodes substitution class and arrow annotations call out compact kataegis-like clusters.

Data use and provenance

The mutation calls are the TCGA BRCA sample distributed with the MIT-licensed maftools rainfall example. During data loading, the package selects the most mutated sample, calculates inter-event distances, and detects compact six-mutation windows. The trimmed MAF omits NCBI_Build; its coordinates and gene-symbol overlaps identify it as hg19. The displayed hg19 RefSeq gene bodies come from the assembly-wide resource independently prepared from the official UCSC refGene table. The track design follows the MutGlyph gene-annotation pattern. GenomeSpy performs the lane packing and scored-label filtering while rendering.

What to notice

The focal chr8 view contains three nearby compact mutation clusters. The gene track places CPQ, TSPYL5, MTDH, LAPTM4B, and neighboring RefSeq genes on the same zoomable locus scale, providing context without implying that a nearby gene caused a cluster.

Python implementation

The rainfall points, cluster arrows, and full hg19 gene table are vertically concatenated. Their x scale and genome axis are shared, while their unrelated y scales remain independent. The domain opens on chr8 but does not pre-filter the gene table, so annotation bodies remain available when panning. GenomeSpy packs overlapping genes into three lanes and retains non-overlapping labels according to their display score.

Code

"""Rainfall plot.

Inter-mutation distance across the genome for a single tumor sample. Point
color encodes substitution class and arrow annotations call out compact
kataegis-like clusters.
"""

from __future__ import annotations

import genome_spy as gs
from genome_spy.datasets._annotations import refseq_gene_bodies
from genome_spy.datasets._mutation import brca_rainfall_data
from genome_spy.schema import GenomeAxis, Legend, Scale


CONVERSION_ORDER = ["C>T", "C>G", "C>A", "T>C", "T>A", "T>G"]
CONVERSION_COLORS = ["#f64b3c", "#4f63c9", "#2891e8", "#f6b617", "#4caf50", "#f7931a"]
GENOME_DOMAIN = [{"chrom": "chr1"}, {"chrom": "chrY"}]


# Load mutation distances and cluster annotations for one tumor sample.
data = brca_rainfall_data()
points = data["points"]
change_points = data["change_points"]
y_domain = [0.0, data["y_max"]]
sample_name = data["sample"]
assembly = data["reference_build"]
# Load gene annotations for the same genome version.
genes = refseq_gene_bodies(assembly)

# Label chromosomes and give neighboring chromosomes alternating backgrounds.
axis = (
    GenomeAxis()
    .title("Genomic position")
    .chromGrid(True)
    .chromGridOpacity(0.16)
    .chromGridFillEven("#f5f7fb")
    .chromGridFillOdd("#ffffff")
    .chromLabels(True)
    .chromLabelFontSize(11)
    .chromTicks(True)
    .grid(False)
)

mutation_scale = Scale().domain(CONVERSION_ORDER).range(CONVERSION_COLORS)
mutation_legend = (
    Legend()
    .orient("top-right")
    .direction("horizontal")
    .columns(3)
    .symbolSize(80)
    .title("Substitution class")
)

# Plot the distance from each mutation to the previous one; low points are close
# together. offset=1 accounts for mutation positions being numbered from one.
rainfall_points = (
    gs.Chart(points)
    .mark_point(size=18, filled=True, opacity=0.95)
    .encode(
        x=gs.Locus("chrom", "pos", offset=1),
        y=gs.Y("log10_distance:Q")
        .scale(reverse=False, domain=y_domain)
        .title("log10 inter-event distance"),
        color=gs.Color("con_class:N").scale(mutation_scale).legend(mutation_legend),
    )
)

# Draw pointers to the annotated mutation clusters.
change_point_stems = (
    gs.Chart(change_points)
    .mark_rule(color="#111111", size=1.2)
    .encode(
        x=gs.Locus("chrom", "start", offset=1),
        y=gs.Y(gs.datum(0), type="quantitative")
        .scale(reverse=False, domain=y_domain)
        .title("log10 inter-event distance"),
        y2=gs.Y2("arrow_y"),
    )
)

# Add an arrowhead at the tip of each pointer.
change_point_heads = (
    gs.Chart(change_points)
    .mark_point(shape="triangle-up", size=60, filled=True, color="#111111")
    .encode(
        x=gs.Locus("chrom", "start", offset=1),
        y=gs.Y("arrow_y:Q").scale(reverse=False, domain=y_domain),
    )
)

# Combine the mutation points and cluster pointers.
rainfall_track = (rainfall_points + change_point_stems + change_point_heads).properties(
    name="rainfall-track",
    title=sample_name,
    height=300,
)

# Choose the gene details to show on hover.
gene_tooltip = [
    gs.Tooltip("symbol:N").title("Gene"),
    gs.Tooltip("identifier:N").title("RefSeq locus"),
    gs.Tooltip("chrom:N").title("Chromosome"),
    gs.Tooltip("start:Q").title("Start").format(",d"),
    gs.Tooltip("end:Q").title("End").format(",d"),
    gs.Tooltip("strand:N").title("Strand"),
]

# Draw genes as arrows showing their reading direction; reveal them on zoom.
gene_bodies = (
    gs.Chart()
    .mark_arrow(
        style="arrow-block",
        fill="#d5d9de",
        stroke="#59636e",
        strokeWidth=1,
        yOffset=8,
        size=7,
        tooltip=gs.HandledTooltip(handler="default"),
    )
    .encode(
        x=gs.Locus("chrom", "start"),
        x2=gs.Locus("chrom", "end"),
        direction=gs.Direction("strand:N").scale(
            domain=["+", "-"], range=["forward", "reverse"]
        ),
        tooltip=gene_tooltip,
    )
    .properties(
        opacity=gs.dynamic_opacity(unitsPerPixel=[100000, 40000], values=[0, 1])
    )
)

# Hide overlapping names, using the supplied scores to choose which to keep.
gene_labels = (
    gs.Chart()
    .transform_measure_text(field="symbol", as_="label_width", fontSize=11)
    .transform_filter_scored_labels(
        pos="linear_start",
        pos2="linear_end",
        asMidpoint="label_position",
        score="score",
        width="label_width",
        lane="lane",
        padding=5,
    )
    .mark_text(
        baseline="middle",
        align="center",
        # Keep names inside the left and right edges of the track.
        clip="x",
        yOffset=-2,
        size=11,
        color="#20262d",
        tooltip=gs.HandledTooltip(handler="default"),
    )
    .encode(
        x=gs.X("label_position:L"),
        text="symbol:N",
        tooltip=gene_tooltip,
    )
)

# Combine gene shapes and names, leaving room above each row for the text.
gene_track = (
    (gene_bodies + gene_labels)
    .properties(
        name="refseq-genes",
        data=genes,
        title=gs.title("RefSeq genes", orient="left", offset=8),
        height=gs.step(24),
        padding=gs.Paddings(top=10),
    )
    .encode(
        y=gs.Y("lane:O")
        .scale(
            type="index",
            domain=[0, 3],
            reverse=True,
            align=0,
            paddingInner=0.4,
            paddingOuter=0.5,
            zoom=False,
        )
        .axis(None)
    )
    .transform_linearize_genomic_coordinate(
        chrom="chrom",
        pos=["start", "end"],
        as_=["linear_start", "linear_end"],
    )
    # Put overlapping genes on separate rows, showing up to three rows.
    .transform_collect(sort=gs.compare(field=["linear_start", "linear_end"]))
    .transform_pileup(
        start="linear_start",
        end="linear_end",
        as_="lane",
        preference="strand",
        preferredOrder=["-", "+"],
    )
    .transform_filter(gs.datum.lane < 3)
)

# Put genes below the mutation plot so both tracks move together when zooming.
chart = (
    (rainfall_track & gene_track)
    .properties(
        assembly=assembly,
        width="container",
        scales=gs.scales(x=gs.Scale(domain=GENOME_DOMAIN)),
        axes=gs.axes(x=axis),
        spacing=8,
        description=(
            "A genome-wide rainfall plot with substitution classes, arrowed "
            "change-point annotations, and aligned RefSeq gene bodies."
        ),
    )
    .resolve_scale(x="shared", y="independent")
    .resolve_axis(x="shared", y="independent")
    .configure_view(stroke="lightgray")
)

Show specification