TCGA ovarian cancer GISTIC landscape¶
Recurrent copy-number scores and their amplification and deletion peaks are shown over a shared hg19 genomic axis.
Data use and provenance
The results shown here are in whole or part based upon data generated by the TCGA Research Network.
The displayed GISTIC2 output files are open-access TCGA data produced by the Broad Institute TCGA Genome Data Analysis Center.
Source: Broad Institute TCGA Genome Data Analysis Center (2016), TCGA OV-TP
CopyNumber GISTIC2 Level 4, Firehose run 2016-01-28, GISTIC2.0.22, hg19.
Download the source archive.
The aligned hg19 RefSeq bodies come from the assembly-wide resource
independently prepared from the official UCSC refGene table. The track design
follows the
MutGlyph gene-annotation pattern.
What to notice¶
The upper track places amplification scores above zero and deletion scores below it. The middle track distinguishes each recurrent event’s peak, wide peak, and broader region by stroke width and opacity. The lower gene track places the strong 19p13.3 deletion peak in local context, including STK11, GPX4, PTBP1, and other RefSeq genes.
Python implementation¶
The chart loads packaged copies of the complete GISTIC score and lesion tables and opens at chr19:400,000–2,000,000. Formula, regex extraction, regex folding, filtering, and projection transforms derive the fields needed by the score and lesion tracks in the browser. The initial chr19 domain does not pre-filter the GISTIC or complete hg19 gene tables. GenomeSpy packs overlapping gene bodies, filters colliding labels by score, and keeps all three tracks on one shared zoomable locus scale.
See the official GenomeSpy example for the original input-file and transform details.
Code¶
"""TCGA ovarian cancer GISTIC landscape.
Recurrent copy-number scores and their amplification and deletion peaks are
shown over a shared hg19 genomic axis.
"""
import genome_spy as gs
from genome_spy.datasets._annotations import refseq_gene_bodies
from genome_spy.datasets._gistic import tcga_ov_gistic_data
GENOME_DOMAIN = [{"chrom": "chr1"}, {"chrom": "chrY"}]
event_colors = gs.Scale(
domain=["Amp", "Del"],
range=["#e45756", "#4c78a8"],
)
# Load the prepared GISTIC results and matching gene annotations.
data = tcga_ov_gistic_data()
genes = refseq_gene_bodies("hg19")
# Draw a zero line between amplifications and deletions.
zero_line = (
gs.Chart([{"value": 0}])
.mark_rule(color="black", opacity=0.3)
.encode(y=gs.Y("value:Q").title(None))
.properties(name="zero-line")
)
# Show amplification scores above zero and deletion scores below it.
q_values = (
gs.Chart(data["scores"])
.transform_formula(
expr=gs.datum["-log10(q-value)"] * gs.expr.if_(gs.datum.Type == "Del", -1, 1),
as_="-log10(q-value)",
)
.mark_rect(minOpacity=1)
.encode(
x=gs.Locus("Chromosome", "Start"),
x2=gs.Locus("Chromosome", "End"),
y=gs.Y("-log10(q-value):Q"),
color=gs.Color("Type:N").scale(event_colors),
)
.properties(name="q-value-rects")
)
# Mark the significance threshold on both sides of zero.
thresholds = (
gs.Chart([{"value": 0.602}, {"value": -0.602}])
.mark_rule(strokeDash=[3, 1], color="black", opacity=0.3)
.encode(y=gs.Y("value:Q").title(None))
.properties(name="q-value-thresholds")
)
# Combine the scores and reference lines in the top track.
score_track = gs.layer(zero_line, q_values, thresholds).properties(
name="gistic-q-value",
title=gs.Title(
text="GISTIC q-values from the TCGA OV-TP cohort",
style="overlay-title",
),
)
# Show where each peak lies within its wider affected region.
lesion_track = (
gs.Chart(data["lesions"])
.transform_regex_extract(
field="Unique Name",
regex=r"^(Amplification|Deletion) Peak[ ]+\d+$",
as_=["Type"],
skipInvalidInput=True,
)
.transform_filter(gs.datum.Type)
# Read the wide peak, peak, and region positions from their separate columns.
.transform_regex_fold(
columnRegex=[r"^(.*) Limits$"],
asValue=["limits"],
asKey="Segment type",
)
.transform_regex_extract(
field="limits",
regex=r"^(chr[^:]+):(\d+)-(\d+)",
as_=["Chrom", "Start", "End"],
)
.transform_project(
fields=[
"Segment type",
"Chrom",
"Start",
"End",
"Type",
"Descriptor",
"q values",
]
)
.mark_rule(minLength=2)
.encode(
x=gs.Locus("Chrom", "Start"),
x2=gs.Locus("Chrom", "End"),
y=gs.Y("Type:N")
.scale(domain=["Amplification", "Deletion"], padding=0.2)
.title(None),
color=gs.Color("Type:N").scale(
domain=["Amplification", "Deletion"],
range=["#e45756", "#4c78a8"],
),
opacity=gs.Opacity("Segment type:N").scale(
type="ordinal",
domain=["Wide Peak", "Peak", "Region"],
range=[0.3, 1, 0.3],
),
size=gs.Size("Segment type:N").scale(
type="ordinal",
domain=["Wide Peak", "Peak", "Region"],
range=[11, 15, 2],
),
)
.properties(
name="gistic-all-lesions",
title=gs.Title(text="Regions and peaks", orient="none"),
height=gs.step(20),
)
)
# 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 without removing the gene shapes beneath them.
# The supplied scores decide which names 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"),
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)
)
# Stack the scores, affected regions, and genes so they zoom together.
chart = (
gs.vconcat(score_track, lesion_track, gene_track)
.properties(
assembly="hg19",
name="gistic-track",
width="container",
scales=gs.scales(x=gs.Scale(domain=GENOME_DOMAIN)),
axes=gs.axes(x=gs.GenomeAxis(title="Genomic position")),
spacing=8,
description=(
"TCGA OV-TP GISTIC2 copy-number scores, recurrent lesions, and "
"aligned RefSeq gene bodies across the hg19 genome."
),
)
.resolve_scale(x="shared", y="independent")
.resolve_axis(x="shared", y="independent")
.configure_legend(disable=True)
.configure_view(stroke="lightgray")
)