Skip to main content

Reviewing a whole SV callset

TL;DR: jb2export batch renders every junction in a somatic callset as a breakpoint split view, so triage becomes a directory of images. The matched normal, rendered the same way, is the control.

Prerequisites

  • @jbrowse/img, which puts jb2export on your PATH
  • nothing to download: the callset, the tumor reads and the matched normal are all hosted
npm install -g @jbrowse/img

Where the data comes from

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

COLO829 and its matched normal

COLO829 is a melanoma cell line with a matched normal, COLO829BL, and a community reference for somatic structural-variant calling. The multi-hop tutorial follows one event in this callset all the way down; this one renders every junction at a glance.

The contact sheet

A junction joins two loci, and those two loci are exactly the two panels of a breakpoint split view. So a callset renders straight into a review queue:

curl -fO https://jbrowse.org/demos/cancer_sv/COLO829.somatic-sv.vcf.gz
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
Warning: skipped 35 record(s), e.g. line 271: names no second locus (INS)
[########################] 100% 100/100
wrote 100/100 images to tumor

100 junctions. Insertions name one locus, so there is no second panel to stack and the run counts them out in its warning; what remains collapses because a caller writes each breakend pair twice. That is the same 100 that sv_multihop.py chains reports on this file in the multi-hop tutorial, junction for junction and in the same order.

They agree because neither parses the ALT bracket by hand, which goes wrong four ways, all of them silent:

  • the replacement string may carry inserted sequence either side of the bracket (GTGATGGATTCA[CHR12:72273112[)
  • callers upper-case the mate contig, and CHR12 is not a region hg38 has
  • END= matches inside CIEND=, and the first hit wins
  • the two records of one breakend pair name the same translocation twice

One image per row, written as 1_chr1_33053494-chr6_2919922_gridss12o.png: index first so the directory sorts in callset order, coordinates next so you can find the one you are looking at, and the caller's own ID last so you can go back to the VCF row it came from. A file with no ID column falls back to junction_<n>.

--flank is the setting that decides the picture: a caller's breakend is one base, so the flank is what frames the panel. --dryRun prints the file and loci of every row and renders nothing, and --limit 20 renders the first few, so you can check the framing before committing to the whole callset.

Two flags for a long run:

  • --resume skips a row whose image is already in --outDir, so an interrupted callset continues from where it stopped.
  • --manifest writes manifest.tsv beside the images: one row per junction with its file, both loci, its name, and whether it rendered. The status column is where the failed rows stay readable after the run's output has scrolled past.

--passOnly drops the records the caller has already filtered out. --limit takes the first N in file order, so on an unfiltered callset the two go together.

The reads stream from the hosted CRAM and each image is rendered server-side. The run is a single process, so the module graph loads once for the whole callset, and a --config URL or a --hub is fetched once.

A row that cannot be rendered is reported and the run continues, so a translocation into a contig the assembly does not have costs you that row alone.

A junction is two loci, so that is what batch draws. A connector drawn dashed means the read carrying it has a segment at a locus the frame does not show, and these reads also visit chr10, so this junction wants a third panel. The control belongs beside it: one render per sample, the same --loc list and the same --width.

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

Reads at 1 px apiece (featureHeight:super-compact) is what keeps six pileups on one screen.

A curve says two loci are joined in this sample; a contig says in what order and in which orientation, which takes a reconstruction step. The multi-hop tutorial builds that contig from these same reads, and rendering it is another jb2export run with a different --assembly, since a derivative allele is an assembly like any other.

The same export over the normal

The whole callset gets the same treatment, one directory per track:

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

Put the two directories side by side and the somatic calls are the ones with curves in tumor/ and none in normal/, rendered at the same flank and width.

Reading the sheet

What each picture says:

  • a fan of curves at both breakends is the junction as the reads describe it
  • nothing connecting the panels means the reads do not support the caller's coordinates, which is either a false call or a breakpoint placed far enough off that --flank missed it. Re-render that row wider
  • curves in the normal too means germline, not somatic
  • a dense fan in a region of ragged coverage is usually a repeat. The connectors are drawn from what the aligner said, so a read mismapped into a repeat contributes a confident-looking curve

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 for the interactive version of the picture you just looked at, with the gene track and the read details attached.

For a junction that turns out to be one hop of something larger, the alignments track menu's Reconstruct derivative allele... groups the reads in view by the route their split alignments describe. COLO829's der(3), the junction in the figures above, is three junctions across three chromosomes, and the multi-hop tutorial follows it the rest of the way.

The junction the figures above are of, taken the rest of the way in the browser: the tumor track menu, the routes the reads describe with the count behind each one, and Breakpoint split view replacing the window with one panel per segment of the route. Open the session this video starts in ↗

Other callers

The recipe follows the format. Anything that writes breakends or symbolic SVs to a VCF goes through the same two commands:

  • cuteSV, Sniffles, pbsv, Delly, Manta, GRIDSS all write a VCF that sv_multihop.py bedpe reads directly
  • LINX publishes its clusters and chained links as TSVs rather than VCF. Convert the junction columns to the six BEDPE columns with awk and the rest of this page is unchanged; the cluster and chain ids make good --outDir names, so one directory per cluster gives you a chromothripsis event as a contact sheet
  • PURPLE copy-number segments are not junctions. Convert the segment TSV to a bedGraph, bedGraphToBigWig it, and add it as a --bigwig so every image carries the copy number under the reads

Ordering breakends into a derivative chromosome needs allele-specific copy number and a centromere constraint, which is what LINX does with PURPLE's purity and ploidy; run it, and load its output here.

Reproduce it end to end

Everything on this page is the commands above against hosted files. The figure is three jb2export invocations: two breakpoint renders with one --loc per panel, one per sample, and a plain render of the derivative assembly, which are the sv_review_tumor, sv_review_normal and sv_review_derivative specs in website/scripts/specs/jbrowse-img.ts.

See also

References