Annotate genomic intervals

Download the notebook

Brush the chart, enter a name, and click Save annotation. New marks appear in Your annotations. Click Download BED to keep them.

Loading chart…

This demo keeps annotations in your browser; reloading clears them. Use the notebook above to keep the records in a Python list called annotations.

The assembly is hg38. BED uses 0-based, half-open intervals and contains chromosome, start, end, and name. Notes are not included in BED. Brushes spanning chromosomes are rejected. Gene bodies cover the initial region only.

Code

The chart and Python hooks below come directly from the downloadable notebook. The web demo uses JavaScript for the same interactions; Python hooks run in a notebook.

Chart specification (Python)
import genome_spy as gs
from genome_spy.datasets._annotations import refseq_gene_bodies

ASSEMBLY = "hg38"
DOMAIN = [{"chrom": "chr5", "pos": 177500000}, {"chrom": "chr5", "pos": 177700000}]
genes = refseq_gene_bodies(ASSEMBLY)
genes = genes.loc[
    (genes.chrom == "chr5")
    & (genes.start < DOMAIN[1]["pos"])
    & (genes.end > DOMAIN[0]["pos"])
]

signal = (
    gs.Chart(
        data=gs.lazy.bigwig("https://data.genomespy.app/genomes/hg38/hg38.gc5Base.bw")
    )
    .mark_rect(color="#547aa5", minWidth=0.5)
    .encode(
        x=gs.Locus("chrom", "start"),
        x2=gs.Locus("chrom", "end"),
        y=gs.Y("score:Q").scale(domain=[0, 100]).title("GC (%)"),
    )
    .properties(height=130)
)
gene_track = (
    gs.Chart(genes)
    .mark_rect(color="#8c929c", cornerRadius=2)
    .encode(
        x=gs.Locus("chrom", "start"),
        x2=gs.Locus("chrom", "end"),
        y=gs.Y("symbol:N").axis(title=None),
        tooltip=["symbol:N", "chrom:N", "start:Q", "end:Q"],
    )
    .properties(height=120, title="RefSeq gene bodies (shown region)")
)
saved_marks = (
    gs.Chart(data={"name": "annotations"})
    .mark_rect(color="#d89632", opacity=0.75)
    .encode(
        x=gs.Locus("chrom", "start"),
        x2=gs.Locus("chrom", "end"),
        tooltip=["name:N", "note:N", "chrom:N", "start:Q", "end:Q"],
    )
)
saved_labels = saved_marks.mark_text(align="left", dx=4, color="#253047").encode(
    text="name:N"
)
saved_track = (saved_marks + saved_labels).properties(
    height=26, title="Your annotations"
)
chart = (
    (signal & gene_track & saved_track)
    .properties(
        assembly=ASSEMBLY,
        datasets={"annotations": []},
        width=760,
        spacing=12,
        padding=gs.Paddings(top=10, right=50, bottom=10, left=10),
        scales=gs.scales(x=gs.Scale(domain=DOMAIN)),
    )
    .resolve_scale(y="independent")
    .add_params(
        gs.selection_interval(
            "brush", encodings=["x"], extent="container", on={"type": "mousedown"}
        )
    )
)
Python hooks and controls
import asyncio
import html
from datetime import datetime
from pathlib import Path

import ipywidgets as widgets
import pandas as pd
from IPython.display import FileLink, display

# Rerunning starts a fresh session.
for task_name in ("connection_task", "save_task"):
    if old_task := globals().get(task_name):
        old_task.cancel()
if old_widget := globals().get("widget"):
    old_widget.close()

widget = chart.widget(inline=True, controls=False)
annotations = []
selection = None
save_task = None
display(widget)
name_input = widgets.Text(description="Name:", placeholder="e.g. GC-rich-region")
note_input = widgets.Textarea(description="Note:", placeholder="Optional")
save_button = widgets.Button(description="Save annotation", disabled=True)
export_button = widgets.Button(description="Export BED", disabled=True)
status = widgets.HTML("Connecting…")
table = widgets.HTML("Saved regions appear in the annotation track above.")
export_output = widgets.Output()


def bed_interval(snapshot):
    # Locus snapshots contain zero-based positions, already rounded down by GenomeSpy.
    endpoints = snapshot["complexIntervals"].get("x")
    if not snapshot["active"] or not endpoints:
        raise ValueError("Brush a region first.")
    left, right = endpoints
    if left["chrom"] != right["chrom"]:
        raise ValueError("Choose a region within one chromosome.")
    start, end = left["pos"], right["pos"]
    if not isinstance(start, int) or not isinstance(end, int) or not 0 <= start < end:
        raise ValueError("Select at least one whole base.")
    # Keep a BED half-open interval: do not add 1 to either endpoint.
    return {"chrom": left["chrom"], "start": start, "end": end}


async def save_annotation():
    save_button.disabled = True
    try:
        name = name_input.value.strip()
        note = note_input.value.strip()
        if not name or any(character.isspace() for character in name):
            raise ValueError("Use a name without spaces (e.g. GC-rich-region).")
        async with asyncio.timeout(30):
            region = bed_interval(await selection.get_value())
            record = {**region, "name": name, "note": note}
            await api.datasets.set("annotations", [*annotations, record])
            annotations.append(record)
            table.value = pd.DataFrame(annotations).to_html(index=False, escape=True)
            export_button.disabled = False
            name_input.value = ""
            note_input.value = ""
            status.value = f"Saved {len(annotations)} annotation(s)."
            await selection.clear()
    except Exception as error:
        status.value = html.escape(f"{type(error).__name__}: {error}")
    finally:
        save_button.disabled = False


def start_save():
    global save_task
    if save_task is None or save_task.done():
        save_task = asyncio.create_task(save_annotation())


def on_save(button):
    connection_task.get_loop().call_soon_threadsafe(start_save)


def export_bed(button):
    if not annotations:
        return
    path = Path(f"annotations.hg38.{datetime.now():%Y%m%d-%H%M%S-%f}.bed")
    pd.DataFrame(annotations)[["chrom", "start", "end", "name"]].to_csv(
        path, sep="\t", header=False, index=False, mode="x"
    )
    with export_output:
        export_output.clear_output()
        display(FileLink(str(path)))


async def connect():
    global api, selection
    try:
        async with asyncio.timeout(30):
            api = await widget.get_embed_api()
            selection = await api.params.get_selection("brush")
        status.value = "Ready — brush the chart, enter a name, then save."
        save_button.disabled = False
    except Exception as error:
        status.value = html.escape(
            f"Connection failed: {type(error).__name__}: {error}"
        )


save_button.on_click(on_save)
export_button.on_click(export_bed)
display(
    widgets.VBox(
        [
            name_input,
            note_input,
            widgets.HBox([save_button, export_button]),
            status,
            table,
            export_output,
        ]
    )
)
# Let this cell finish so the notebook can receive browser replies.
connection_task = asyncio.create_task(connect())

Data use and provenance

GC content comes from UCSC hg38, served by GenomeSpy. Gene bodies are packaged UCSC/NCBI RefSeq annotations with overlapping transcripts collapsed—not exon models. Python selects the gene subset; GenomeSpy loads the signal and renders the tracks. See data notices.