Reviewing a whole SV callset

JBrowse v5 beta

The tutorials target the JBrowse v5 beta, and the v4.3.0 release on the download page lacks some of what they show. To install the beta of JBrowse Web, run npm install -g @jbrowse/cli@next, then jbrowse create jbrowse2 --branch v5.0.0-beta.13. Desktop beta builds are coming soon.

The COLO829 somatic structural-variant callset lists over a hundred junctions, the points where two distant loci are joined, and each one is worth checking against the reads. We render every junction as an image of the reads at both ends with jb2export batch, then render the matched normal the same way as the control.

Prerequisites

@jbrowse/img, which puts jb2export on your PATH
a JBrowse for the last section (Web or Desktop)

npm install -g @jbrowse/img

Where the data comes from

Show full URLs

The COLO829 somatic SV callset comes from the ONT open-data release's wf-somatic-variation run (Valle-Inclán et al. 2022), rehosted alongside the cancer SV demo.

The commands below read this file by URL, so there is nothing to download by hand.

the callset these commands fetch directly: jbrowse.org/…/COLO829.somatic-sv.vcf.gzhttps://jbrowse.org/demos/cancer_sv/COLO829.somatic-sv.vcf.gz

Read by URL (no download needed)

the tumor reads the config's COLO829_tumor_ont track streams, Oxford Nanopore R10 from the ONT open-data release: ont-open-data.s3.amazonaws.com/…/COLO829_tumor.ht.cramhttps://ont-open-data.s3.amazonaws.com/colo829_2024.03/wf_somatic_variation/sup/COLO829_tumor.ht.cram
the matched normal the config's COLO829BL_normal_ont track streams: ont-open-data.s3.amazonaws.com/…/PAU59807.d052sup4305mCG_5hmCGvHg38.bamhttps://ont-open-data.s3.amazonaws.com/colo829_2024.03/basecalls/colo829bl/sup/PAU59807.d052sup4305mCG_5hmCGvHg38.bam

Rendering every tumor junction with jb2export batch

jb2export batch renders one image per record, a chromosome a row. A deletion, duplication or inversion has both ends on one chromosome, so it is one row: the two ends side by side under the arc of the reads joining them. A junction between two chromosomes is a breakpoint split view, one panel a chromosome, with the reads joining them drawn as curves:

curl -fO https://jbrowse.org/demos/cancer_sv/COLO829.somatic-sv.vcf.gz
# --track: the track id, then display settings as key:value pairs
jb2export batch --vcf COLO829.somatic-sv.vcf.gz \
  --config https://jbrowse.org/demos/cancer_sv/config.json --assembly hg38 \
  --track COLO829_tumor_ont height:240 \
  --outDir tumor --flank 600 --width 1100

[########################] 100% 135/135
wrote 135/135 images to tumor

For your own callset, point --config at a JBrowse config that holds the assembly and an alignments track:

--assembly is the assembly's name in that config and --track the track's trackId.
The track's reads need their index beside them (.bai or .crai).
The VCF's chromosome names must match the assembly's.

jb2export batch --vcf calls.vcf.gz \
  --config your/config.json --assembly <assemblyName> --track <trackId> \
  --outDir out

A record that fits one window is drawn as a single panel: an insertion names one locus, and a deletion shorter than --flank has both ends in one frame.

batch reads the breakend notation in the ALT column with @gmod/vcf, and renders once a breakend pair that a caller writes as two records naming the same translocation.¹

batch writes one image per record, named 002_chr1_33053494-chr6_2919922_r_0_0.png: the record's index, so the directory sorts in callset order, then the coordinates and the caller's ID where the record has one.

A breakend is one base, so --flank sets the window drawn around it. Further options check the framing and manage a long run:

--dryRun prints the file and loci of every row and renders nothing, and --limit 20 renders the first few, to check the framing before the whole callset
--resume skips a row whose image is already in --outDir. A --limit run names its images as the whole run will, so the whole run picks them up
--manifest writes manifest.tsv beside the images, one row per image: file, loci, name, EVENT, whether it rendered, and line (the record's line in the VCF, which joins a row back to any column of the callset). Three more columns describe the image:
- links is the count of split reads (reads aligned in pieces to different places) with pieces in more than one window
- alt is the reads with the record's ALT over the reads covering it, 17/41, one pair per alignments track, for a record on one chromosome whose ALT says what a read carries: an SNV, an indel, a <DEL> or an <INS>
- spec is the view as a session spec, which opens the same windows and tracks in JBrowse
--passOnly drops records the caller filtered out. --limit takes the first N in file order, so on an unfiltered callset the two go together
--jobs sets how many processes render, each about a gigabyte. The default is half the cores, up to four

During a run, batch:

streams the reads from the hosted CRAM and fetches a --config URL or --hub once
reports a row it cannot render and continues, and --resume retries that row
loads every panel as if you had pressed Force load, because deep long reads can put even a --flank window over a track's size limit

COLO829's der(3), a derivative chromosome 3 joined from pieces of chr3, chr10 and chr12, has reads that visit all three loci. We render that event with jb2export breakpoint and one --loc per panel. The matched normal gets the same --loc list and --width, and sits beside it as the control:

jb2export breakpoint \
  --config https://jbrowse.org/demos/cancer_sv/config.json --assembly hg38 \
  --track COLO829_tumor_ont height:130 force:true featureHeight:super-compact \
  --loc chr3:25,358,511-25,359,711 \
  --loc chr10:58,716,962-58,718,162 \
  --loc chr12:72,272,512-72,273,712 \
  --width 1000 --out der3_tumor.png

The three loci of COLO829's der(3), chr3 then chr10 then chr12, at the same width in every panel. The tumor nanopore reads have a solid curve at every breakend and the matched normal has none at them. On the right, the same three loci as one reconstructed contig.

featureHeight:super-compact in the jb2export breakpoint command draws each read 1 px tall, which fits six pileups on one screen.

A dashed connector marks a read with a segment at a locus outside the frame. The complex rearrangements tutorial builds the reconstructed contig from these reads, and rendering it is another jb2export run with the contig as --assembly.

Rendering the matched normal as the control

Render the matched normal into a second directory:

jb2export batch --vcf COLO829.somatic-sv.vcf.gz \
  --config https://jbrowse.org/demos/cancer_sv/config.json --assembly hg38 \
  --track COLO829BL_normal_ont height:240 \
  --outDir normal --flank 600 --width 1100

A somatic call has curves in tumor/ and none in normal/. The same file name in both directories puts each call beside its control:

Three rows of the two batch directories, tumor on the left and the matched normal on the right, each labelled with its file name. The chr7 junction has a fan of curves in the tumor and none in the normal. The chr1 to chr19 junction has curves in both. The chr2 deletion has no curve in either; one tumor read has it as a gap through both panels.

The caller filed the chr1 to chr19 junction as somatic, but the curves in the normal mark it as germline.

Interpreting the curves: somatic, germline or unsupported

A fan of curves at both breakends is the junction as the reads describe it.
With no curve between the panels, the reads give no support for the caller's coordinates: either the call is false, or the breakpoint is far enough off that --flank missed it. Re-render that row wider.
Curves in the normal as well mean the variant is germline.
A dense fan in a region of ragged coverage is usually a repeat. JBrowse draws the connectors from the aligner's output, so a read mismapped into a repeat adds a confident-looking curve.

Add --manifest to both batch runs. Sort tumor/manifest.tsv on links to put the calls no split read joins at the top; the same column in normal/manifest.tsv shows which calls the normal has too. A deletion short enough for one alignment to hold draws a gap through both windows and no arc, so its reads are in alt and its links is zero.

Opening a call in the browser

Take the coordinates from an image's filename, open the SV inspector on the same VCF, and click through to the breakpoint split view, with the gene track and read details attached.

Rendering calls from other SV callers

Anything that writes breakends or symbolic SVs (such as <DEL>) to a VCF goes through jb2export batch and jb2export breakpoint:

cuteSV, Sniffles, pbsv, Delly, Manta and GRIDSS write a VCF that --vcf reads directly.
LINX writes clusters and chained links as TSVs. Convert the junction columns to the six BEDPE columns with awk and pass the file as --bedpe in place of --vcf; one --outDir per cluster renders a chromothripsis event (a chromosome shattered and rejoined at once) as one directory.
PURPLE writes copy-number segments. Convert the segment TSV to a bedGraph, run bedGraphToBigWig on it, and add it as a --bigwig to draw the copy number under the reads in every image.

The COLO829 callset does not group its junctions, so the der(3) figure above needed a hand-written --loc list. A caller that files the junctions of one rearrangement under VCF 4.4's EVENT key, as DRAGEN and the C-GIAB benchmark do, gets that image from batch itself: every event visiting more than two loci is drawn once more as event_<n>_<label>, one panel per locus in contig order. Severus writes the same grouping as CLUSTERID, and the SV inspector guide shows how to rename it to EVENT.

Ordering breakends into a derivative chromosome needs allele-specific copy number and a centromere constraint, which LINX derives from PURPLE's purity and ploidy.

Citations

Valle-Inclán JE, et al. A multi-platform reference for somatic structural variation detection. Cell Genomics (2022). https://doi.org/10.1016/j.xgen.2022.100139
Shale C, et al. Unscrambling cancer genomes via integrated analysis of structural variation and copy number. Cell Genomics (2022). https://doi.org/10.1016/j.xgen.2022.100112

Notes

The parser also handles cases a hand-written one gets wrong without raising an error: inserted sequence either side of the bracket (GTGATGGATTCA[CHR12:72273112[), mate contigs that callers upper-case (hg38 has no contig CHR12), and an END= search that would otherwise match inside CIEND=, where the first hit wins. ↩

Feedback on this tutorial is welcome: contact us.

Reviewing a whole SV callset

Prerequisites#

Where the data comes from#

Rendering every tumor junction with jb2export batch#

Rendering the matched normal as the control#

Interpreting the curves: somatic, germline or unsupported#

Opening a call in the browser#

Rendering calls from other SV callers#

See also#

Citations#