Stacked genome browser tracks¶
Several lazy hg38 resources share one zoomable genomic axis: signal tracks, regulatory intervals, sequence, and RefSeq gene annotations.
Code¶
"""Stacked genome browser tracks.
Several lazy hg38 resources share one zoomable genomic axis: signal tracks,
regulatory intervals, sequence, and RefSeq gene annotations.
"""
from __future__ import annotations
import genome_spy as gs
# Start with a 20 kb region on chromosome 7.
DOMAIN = [
{"chrom": "chr7", "pos": 55100000},
{"chrom": "chr7", "pos": 55120000},
]
CCRE_WINDOW_SIZE = 1_000_000
signal_tooltip = [
gs.Tooltip("chrom:N").title("Chromosome"),
gs.Tooltip("start:Q").title("Start").format(",d"),
gs.Tooltip("end:Q").title("End").format(",d"),
gs.Tooltip("score:Q").title("Value"),
]
# Show GC content, loading just the region being viewed.
gc_track = (
gs.Chart(gs.lazy.bigwig("https://data.genomespy.app/genomes/hg38/hg38.gc5Base.bw"))
.mark_rect(color="#6c8ebf", minWidth=0.5, minOpacity=1)
.encode(
x=gs.Locus("chrom", "start"),
x2=gs.Locus("chrom", "end"),
y=gs.Y("score:Q")
.scale(domain=[0, 100])
.axis(title=None, grid=True, gridDash=[2, 2]),
tooltip=signal_tooltip,
)
.properties(
name="gc-content", title=gs.title("GC (%)", style="track-title"), height=80
)
)
# Add conservation scores to compare with the GC-content track.
conservation_track = (
gs.Chart(
gs.lazy.bigwig(
"https://hgdownload.soe.ucsc.edu/goldenPath/hg38/phyloP100way/"
"hg38.phyloP100way.bw"
)
)
.mark_rect(color="#c77c8a", opacity=0.75, minWidth=0.5)
.encode(
x=gs.Locus("chrom", "start"),
x2=gs.Locus("chrom", "end"),
y=gs.Y("score:Q")
.scale(domain=[-5, 5])
.axis(title=None, grid=True, gridDash=[2, 2]),
tooltip=signal_tooltip,
)
.properties(
name="phylop-100way",
title=gs.title("phyloP", style="track-title"),
height=80,
)
)
# Color candidate regulatory regions by type.
ccre_marks = (
gs.Chart(
gs.lazy.bigbed(
"https://data.genomespy.app/sample-data/encodeCcreCombined.hg38.bb",
windowSize=CCRE_WINDOW_SIZE,
)
)
.mark_rect(minWidth=0.5)
.encode(
x=gs.Locus("chrom", "chromStart"),
x2=gs.Locus("chrom", "chromEnd"),
color=gs.Color("ucscLabel:N")
.scale(
domain=["prom", "enhP", "enhD", "K4m3", "CTCF"],
range=["#e45756", "#f2a541", "#f6cf65", "#d99ac5", "#4f9fc4"],
)
.legend(title="cCRE type", columns=2),
tooltip=[
gs.Tooltip("name:N").title("cCRE"),
gs.Tooltip("ucscLabel:N").title("Type"),
gs.Tooltip("chrom:N").title("Chromosome"),
gs.Tooltip("chromStart:Q").title("Start").format(",d"),
gs.Tooltip("chromEnd:Q").title("End").format(",d"),
],
)
)
ccre_zoom_message = (
gs.Chart([{}])
.mark_text(text="Zoom in to see cCREs", color="#555", size=12, tooltip=None)
.encode(x=gs.value(0.5), y=gs.value(0.5))
)
ccre_track = gs.multiscale(
ccre_zoom_message,
ccre_marks,
stops=[gs.expr(CCRE_WINDOW_SIZE / gs.expr.max(gs.Expression("width"), 1))],
name="ccre",
title=gs.title("cCRE", style="track-title"),
height=32,
)
# Draw each DNA base as a colored tile with its letter on top.
sequence_rects = gs.Chart().mark_rect()
sequence_labels = (
gs.Chart()
.mark_text(
size=13,
fitToBand=True,
paddingX=1.5,
paddingY=1,
opacity=0.7,
flushX=False,
)
.encode(color=gs.value("black"), text=gs.Text("base:N"))
)
# Load the DNA sequence and reveal it gradually as the reader zooms in.
sequence_track = (
(sequence_rects + sequence_labels)
.properties(
name="sequence",
title=gs.title("Sequence", style="track-title"),
height=32,
data=gs.lazy.indexed_fasta(
"https://data.genomespy.app/genomes/hg38/hg38.fa", windowSize=30_000
),
opacity=gs.dynamic_opacity(unitsPerPixel=[100, 10], values=[0, 1]),
)
.encode(
x=gs.Locus("chrom", "pos"),
color=gs.Color("base:N")
.scale(
domain=["A", "C", "T", "G", "N"],
range=[
"#7BD56C",
"#FF9B9B",
"#86BBF1",
"#FFC56C",
"#E0E0E0",
],
)
.legend(title="DNA base", columns=2),
tooltip=[
gs.Tooltip("chrom:N").title("Chromosome"),
gs.Tooltip("pos:Q").title("Position").format(",d"),
gs.Tooltip("base:N").title("Base"),
],
)
# Give each base its own row and genomic position.
.transform_flatten_sequence(field="sequence", as_=["rawPos", "base"])
.transform_formula(expr=gs.expr.upper(gs.datum.base), as_="base")
.transform_formula(expr=gs.datum.rawPos + gs.datum.start, as_="pos")
)
# Draw exon blocks along each transcript.
exons = (
gs.Chart()
.mark_rect(minOpacity=0.2, minWidth=0.5)
.encode(
x=gs.X("exonStart:L"),
x2=gs.X2("exonEnd"),
tooltip=[
gs.Tooltip("symbol:N").title("Gene"),
gs.Tooltip("exonStart:Q").title("Exon start").format(",d"),
gs.Tooltip("exonEnd:Q").title("Exon end").format(",d"),
],
)
.transform_project(fields=["symbol", "_lane", "_start", "exons"])
.transform_flatten_compressed_exons(start="_start")
.properties(name="exons")
)
# Connect the exons with a thin line spanning the transcript.
bodies = (
gs.Chart()
.mark_rule(minLength=0.5, size=1)
.encode(
x=gs.X("_start:L"),
x2=gs.X2("_end"),
search=gs.Search("symbol"),
tooltip=[
gs.Tooltip("symbol:N").title("Gene"),
gs.Tooltip("chrom:N").title("Chromosome"),
gs.Tooltip("start:Q").title("Start").format(",d"),
gs.Tooltip("strand:N").title("Strand"),
],
)
.properties(name="bodies", title="Gene annotations")
)
# Show more gene detail as the reader zooms in.
transcripts = (
(exons + bodies)
.properties(
name="transcripts",
opacity=gs.dynamic_opacity(unitsPerPixel=[100000, 40000], values=[0, 1]),
)
.encode(color=gs.value("#909090"))
)
# Add gene names with more information available on hover.
labels = (
gs.Chart()
.mark_text(size=11, yOffset=7, tooltip=gs.HandledTooltip(handler="refseqgene"))
.encode(x=gs.X("_centroid:L"), text=gs.Text("symbol:N"))
.properties(name="labels")
)
# Put a direction arrow just beside each gene name.
arrows = (
gs.Chart()
.mark_point(yOffset=7, size=50, tooltip=None)
.encode(
x=gs.X("_centroid:L"),
dx=gs.Dx(
gs.expr(
(gs.datum._textWidth / 2 + 5)
* gs.expr.if_(gs.datum.strand == "-", -1, 1)
),
type="quantitative",
).scale(None),
color=gs.value("black"),
shape=gs.Shape("strand:N")
.scale(domain=["-", "+"], range=["triangle-left", "triangle-right"])
.legend(None),
)
.properties(
name="arrows",
opacity=gs.dynamic_opacity(unitsPerPixel=[100000, 40000], values=[0, 1]),
)
)
# Hide overlapping names, using the supplied scores to choose which to keep.
symbols = (
(labels + arrows)
.properties(name="symbols")
.transform_measure_text(field="symbol", as_="_textWidth", fontSize=11)
.transform_filter_scored_labels(
lane="_lane",
score="score",
width="_textWidth",
pos="_centroid",
padding=5,
)
)
# Load the RefSeq table for the gene shapes and labels.
refseq_track = (
gs.layer(transcripts, symbols)
.properties(
name="refseq-track",
title=gs.title("Genes", style="track-title"),
height=gs.step(23),
data=gs.Data(
url="https://data.genomespy.app/genomes/hg38/refSeqGenes-hg38-release232.tsv.gz",
format=gs.data_format(
parse=gs.parse(
symbol="string",
chrom="string",
start="integer",
length="integer",
strand="string",
score="integer",
exons="string",
)
),
),
)
.encode(
y=gs.Y("_lane:O")
.scale(
type="index",
align=0,
paddingInner=0.4,
paddingOuter=0.2,
domain=[0, 3],
reverse=True,
zoom=False,
)
.axis(None),
)
# Find each transcript's start, end, and label position.
.transform_linearize_genomic_coordinate(chrom="chrom", pos="start", as_="_start")
.transform_formula(expr=gs.datum._start + gs.datum.length, as_="_end")
.transform_formula(
expr=gs.datum._start + gs.datum.length / 2,
as_="_centroid",
)
# Put overlapping transcripts on separate rows, showing up to three rows.
.transform_collect(sort=gs.compare(["_start"]))
.transform_pileup(
start="_start",
end="_end",
as_="_lane",
preference="strand",
preferredOrder=["-", "+"],
)
.transform_filter(gs.datum._lane < 3)
)
# Stack the five tracks so they move together when zooming or panning.
chart = (
(gc_track & conservation_track & ccre_track & sequence_track & refseq_track)
.properties(
assembly="hg38",
title="Stacked hg38 genome browser tracks",
description=(
"Lazy signal, regulatory, sequence, and RefSeq tracks aligned to "
"one zoomable genomic axis."
),
axes=gs.axes(x=gs.GenomeAxis(orient="bottom", title="Genomic position")),
spacing=8,
scales=gs.scales(x=gs.Scale(domain=DOMAIN)),
)
.resolve_scale(x="shared", y="independent", color="independent")
.resolve_axis(x="shared", y="independent")
.resolve_legend(color="collected")
)