Complex rearrangements and derivative alleles

TL;DR: a rearrangement can take several junctions to make, and the genes it brings together say nothing about how many. Search a somatic SV callset for chains of junctions a single long read could cross, rebuild the derivative allele from the reads that span it, and show that reconstruction against the reference as a synteny view.

Prerequisites

nothing to read along. Everything below is for rebuilding the data
Command line tools (JBrowse CLI)
samtools (v1.21 or later)
minimap2
bedGraphToBigWig from the UCSC utilities
python3, for sv_multihop.py
a GRCh38 FASTA, and roughly 40 GB of free disk

On Debian/Ubuntu, apt install samtools minimap2 python3 covers three of those; bedGraphToBigWig is a single static binary from UCSC and node, for the CLI, comes from nodejs.org. sv_multihop.py is one file:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/sv_multihop.py

Where the data comes from

Every file comes out of one ONT open-data release, from its own wf-somatic-variation run (Valle-Inclán et al. 2022):

COLO829 tumor reads (ONT R10, haplotagged): https://ont-open-data.s3.amazonaws.com/colo829_2024.03/wf_somatic_variation/sup/COLO829_tumor.ht.cram
COLO829BL matched normal reads: https://ont-open-data.s3.amazonaws.com/colo829_2024.03/basecalls/colo829bl/sup/PAU59807.d052sup4305mCG_5hmCGvHg38.bam
the somatic SV calls sv_multihop.py searches: https://ont-open-data.s3.amazonaws.com/colo829_2024.03/wf_somatic_variation/sup/COLO829.wf-somatic-sv.vcf.gz
mosdepth coverage regions, tumor: https://ont-open-data.s3.amazonaws.com/colo829_2024.03/wf_somatic_variation/sup/COLO829/qc/coverage/COLO829_tumor.regions.bed.gz
mosdepth coverage regions, normal: https://ont-open-data.s3.amazonaws.com/colo829_2024.03/wf_somatic_variation/sup/COLO829/qc/coverage/COLO829_normal.regions.bed.gz
the GRCh38 build the CRAM decodes against, which derive also realigns the consensus to: https://ont-open-data.s3.amazonaws.com/colo829_2024.03/wf_somatic_variation/sup/GCA_000001405.15_GRCh38_no_alt_analysis_set.fasta

COLO829

COLO829 is a melanoma cell line with a matched normal, COLO829BL, and a community reference for somatic structural-variant calling. The tumor is sequenced deeply enough on ONT R10 that a read crosses a whole rearrangement, which is what the reconstruction below needs.

The coverage lanes beside those reads are the same run's mosdepth output in 50 kb windows, repacked as bigWig:

# awk drops the alt and decoy contigs: bedGraphToBigWig rejects a contig absent
# from the chrom.sizes outright, so one leftover row fails the conversion.
gzip -dc COLO829_tumor.regions.bed.gz | sort -k1,1 -k2,2n |
  awk 'NR==FNR{ok[$1];next} ($1 in ok)' hg38.chrom.sizes - > cov.bg
bedGraphToBigWig cov.bg hg38.chrom.sizes COLO829_tumor.coverage.bw

Multi-hop fusions

Fusion callers generally look for one junction joining two genes. Two genes can also be brought together by a series of junctions, and when the reference segments between them are short, the result is indistinguishable at the transcript level from a simple fusion.

SplitThreader made this concrete in SK-BR-3 (Nattestad et al. 2018): searching the SV graph for short paths between fusion partners found a KLHDC2-SNTB1 fusion that required three variants across three chromosomes. The same search applies to any somatic SV callset, and sv_multihop.py runs it.

Finding the chains

The search needs only the VCF. Two junctions belong to the same chain when an endpoint of one lands close enough to an endpoint of the other that a single read could carry both:

python3 sv_multihop.py chains COLO829.somatic-sv.vcf.gz --min-hops 3

100 distinct junctions in COLO829.somatic-sv.vcf.gz
4 chain(s) of >=3 junctions linked by reference segments <=20000 bp

chain 1: 3 junctions across 3 chromosome(s)
    chr3:25,359,111 <-> chr12:72,273,112
    chr3:25,359,568 <-> chr10:58,717,464
    chr10:58,717,662 <-> chr12:72,273,294
    --loci chr10:58717464,chr12:72273112,chr3:25359111

Those three junctions form a closed cycle, and the whole derivative path is under a kilobase spread across three chromosomes. The genes involved are RARB on chr3, a retinoic-acid receptor that acts as a tumor suppressor, BICC1 on chr10, and TRHDE on chr12.

--max-segment is the longest reference segment one read is assumed to bridge, so set it from your own read-length distribution.

Reads at the breakpoints

At the chr3 breakpoints the tumor pileup becomes soft-clipped bases, because every read crossing the junction has its remainder aligned elsewhere. The matched normal at the same locus is clean.

Left: COLO829 tumor above COLO829BL normal at the two chr3 breakpoints, soft clipping shown. Tumor reads clip where normal reads read through. Right: the same event as a breakpoint split view over every locus the chain visits.

Soft clipping is off by default. Turn it on from the track menu with Show soft clipping. These pileups are deep enough that the track asks before downloading the window; Force load approves it for the rest of the session.

Following the chain across panels

A breakpoint split view, the right half of the figure above, stacks the loci the chain visits and draws the reads that leave one panel and arrive in another.

The reads already know which loci those are and in what order, so the view is built from them. On the tumor track, Launch → Reconstruct derivative allele... lists the routes the reads describe; pick one, set Draw as to Breakpoint split view and choose Replace current view, and the launching view is replaced by a panel per segment of that route, in the order the reads cross it, carrying the tracks that view had.

The route over the chr3 breakpoints: the tumor track menu, the routes the reads describe with the read count and segment sizes behind each one, and Breakpoint split view replacing the window with a panel per segment. The soft-clipped tails at the start return as the curves between panels. Open the session this video starts in ↗

There is one panel per segment: this chain leaves chr3 and returns to it, so it gets two chr3 panels. Every panel opens on the same span, centered on the junction its segment carries, which puts the connecting curves across the frame.

Add → Breakpoint split view builds a view whose loci you already know, one row per panel.

A single record opens the same way. Right-click it in the variant track and choose Open breakpoint split view: one dialog asks for the shape, two stacked panels or one row spanning both breakends, and for the window each panel opens at.

A BND names one partner, so the record on its own is two loci. Follow further breakends at each end reaches the rest of the rearrangement from the callset: at each end of the chain it looks for another junction leaving from the same place, and takes it when there is exactly one. On this record that is three panels, because the chr10 breakend it names has a second junction a couple of hundred bases away whose far end is on chr12. The walk then stops, since the only junction at the chr12 end returns to where the chain began.

The option assumes two junctions leaving one locus belong to one molecule. Two open continuations at a locus end the walk, since the records cannot say which molecule carries which, and so does a continuation leading back into the chain. The reads are the evidence for that assumption, which is the dialog above.

Opening the split view from the record itself: right-click the breakend, set the shape and window in the dialog, and get three panels because the chain runs chr3 to chr10 to chr12. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

For a single read, right-click it and choose Linear read vs ref. That builds a synteny view with the read as its own assembly along the bottom and every locus it touches along the top, the view Ribbon (Nattestad et al. 2021) introduced for this.

Stacked panels describe the event in reference coordinates. Laid out along the derivative, it shows the order and orientation of its pieces. The next section builds that view from the reads already on screen; the one after it rebuilds the allele's sequence, which the base-level checks need.

Reconstructing the derivative allele in the browser

A split read is already an ordered, oriented list of reference intervals, which is what a derivative allele is. With the tumor reads open at a breakpoint, the alignments track menu's Reconstruct derivative allele... groups the reads in view by the path their split alignments describe and offers each path with the number of reads that independently describe it.

A read count ranks the paths, and each row also draws its segments to scale: a rearrangement is usually a long arm carrying short inserts, and a read the aligner chopped into pieces is a row of equal blocks of the same total length.

This dataset is the reconstruction at its easiest: ONT reads tens of kilobases long, an event that moves whole arms, and 29 molecules crossing all three loci. Scored across two published somatic callsets, it recovers 129 of the 130 junctions above 10 kb or between chromosomes, 60% and 65% of those between 1 and 10 kb, and about one in ten below 1 kb — the small ones because the aligner writes them inside a read's CIGAR rather than as a split alignment, where nothing reading SA tags can see them. What the reconstruction needs, and how to weigh a route once it appears, is in SV visualization.

The result is the view type Linear read vs ref produces from a read you right-click, with the lower panel holding the path a group of reads agrees on, so a ribbon carries its whole group's support. Running Linear read vs ref on one supporting read is the single-molecule evidence behind a candidate. Whatever else was open comes along onto the reference panel, so the path is read against the genes it runs through.

Open in new view appends the reconstruction below the pileup it came from; Replace current view puts it in that view's place, which is what the figures here use.

A fold-back on one chromosome

The smallest thing this menu produces is an allele of two segments on one chromosome, so that is where to read it first. COLO829 has one on chr9, where the tumor reads run out at 28,031,837 and resume inverted from 28,059,142: the arm turns around and continues backwards, which is a fold-back. A second call anchors 28 bp from the first, the pattern repeated breakage-fusion-bridge cycles leave behind.

The fold-back's row names the same chromosome twice, once inverted, and its strip draws it as two blocks of one color whose arrows point at each other. Rows tied on read count are ordered by segment count, so the three-segment route through the second anchor sorts above it.

Top: the candidate list at the chr9 fold-back anchors. Bottom: the two-segment fold-back drawn, two windows of chr9 above and the allele below. The ribbons cross where the arm turns around. Open this view in JBrowse ↗

More than one row here means more than one allele: reads reaching the anchors from different directions describe different routes through the same breakpoints, each offered with its own support.

The window is narrower than the event. Reconstruction reads SA tags, so the arm a read returns from can be off screen, and the whole event at this depth is more alignment than the track will fetch.

The der(3) allele, four segments across three chromosomes

The same menu at the chr3 breakpoint this page has been following returns two rows. The first is the whole event: a 52.3 kb arm of chr3, 199 bp of chr10, 183 bp of chr12 inverted, then 8.43 kb of chr3 inverted. The second is that route with the chr12 piece missing, and two reads take it; both cross the same first junction, so the disagreement is about what follows it.

The reconstruction is anchored on the window the pileup was showing, so the reference row is tens of kilobases of chr3 with the two insert loci a few pixels wide at its right-hand end, and the ribbons reaching them are hairlines. The next two sections open those junctions at base scale.

Top: the candidate list over the tumor pileup it was computed from, each row's segments drawn to scale. Bottom: the synteny view the top row draws, reference above and allele below, a ribbon per segment. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The reads describing one allele agree on its junctions: each starts and stops where its own molecule did, and each crosses the allele from whichever end it was sequenced from. Both are properties of the read, so paths are identified by their junctions alone and a chain is folded together with its reverse complement before the counting, which puts this event on one row.

The output is a proposal: reads mismapped into a repeat produce a confident-looking path, so check that the reads run through each junction and that the segments land in the genes the event is supposed to involve.

Those checks are structural. The path is assembled from where the reads' alignments start and stop, so it agrees with those alignments by construction and inherits whatever the aligner got wrong. Deciding whether the reads' own bases support a junction takes a sequence to align them to, which is the next section: derive builds the allele's consensus and realigns the spanning reads onto it, and a wrong junction shows as clipping and mismatches at that position.

Reconstructing the allele's sequence

The candidates above are structure. derive builds the allele's sequence, which is what the base-level checks below need: it pulls the reads spanning every locus, takes the longest as a backbone, polishes it into a consensus with the rest, aligns that consensus back to the reference, and realigns the reads to it.

python3 sv_multihop.py derive \
  --aln COLO829_tumor.ht.cram --ref GRCh38.fa \
  --loci chr10:58717464,chr12:72273112,chr3:25359111 \
  --out der3_RARB --name der3_RARB_BICC1_TRHDE

29 spanning reads
backbone read 8315652b-cd0f-4290-ad6b-51112f93a44a (57,134 bp)
wrote der3_RARB.derivative.fa (39,549 bp supported by >=3 reads)
wrote der3_RARB.vs_reference.paf
    derivative       0-32732   + -> chr3:25,326,821-25,359,568
    derivative   32732-32931   + -> chr10:58,717,463-58,717,662
    derivative   32932-33115   - -> chr12:72,273,111-72,273,294
    derivative   33126-39549   - -> chr3:25,352,683-25,359,111
wrote der3_RARB.derivative_segments.bed

Four contiguous segments: two chr3 arms in opposite orientations, a foldback, with short pieces of chr10 and chr12 spliced in at the turn. Those two fragments are templated insertions, short stretches of other chromosomes captured at a repair junction.

The PAF is a synteny track and the consensus is an assembly, so the reconstruction loads against the reference directly. The BED is the same segments as a feature track on the derivative, each labelled with the interval it came from. Adding --jbrowse-out config.json writes the config that wires those together (both assemblies, the synteny track, the segments and the realigned reads) and prints the URL that opens them as a synteny view.

--genes takes a tabix-indexed GFF3 and projects the reference's own gene annotation through those same segments, so each feature lands in derivative coordinates, clipped where a junction cut it and flipped where a segment is inverted:

python3 sv_multihop.py derive ... --genes ncbiRefSeq.gff.gz

wrote der3_RARB.derivative_genes.gff3 (44 features from 41 reference rows)

This allele carries RARB's first coding exon and its start codon, then the 183 bp of chr12 that the second junction splices in, which is TRHDE coding sequence in reverse, then RARB again inverted.

Ribbons below are colored by the reference chromosome they come from, so the wide green one is the chr3 arm and the crossing ribbons at right are the chr10 and chr12 inserts with chr3 returning inverted.

That last segment is what names the event: an interval the allele has already carried, read back on the other strand, so the derivative turns around on itself, a fold-back. The turn leaves that stretch in the allele twice in opposite orientations, which is the same thing as an inverted duplication, and the two templated inserts are what sits at the turn. Fold-backs are the canonical opening move of a breakage-fusion-bridge cycle.

A read lane sits under each row. Against hg38 it draws split alignments only, so the band over it counts the molecules carrying a junction; on the allele it draws every read realigned there.

The allele's lane thins partway along. The tumor has two chromosome 3s: the allele begins as sequence they share, so reads off the intact homolog align down it and stop where the derivative leaves chr3, and reads off the rearranged copy carry on. The shaded band marks what only the rearranged copy reaches.

Between the two, each junction is drawn once as an arc joining its two ends, with a short tick at each foot lying over the sequence that end keeps: ticks pointing away from each other are a deletion-type join, toward each other a duplication-type, and parallel an inversion.

The reconstructed derivative against its three source loci: RefSeq genes above, the same annotation projected onto the allele below, each segment labelled with the interval it came from, a read lane under each row, and the junctions drawn once each as arcs over the hg38 lane. Shaded on the allele: the stretch only the rearranged chromosome reaches, where the read lane thins to one allele's worth. Open this view in JBrowse ↗

Checking the reconstruction

Zoomed to the kilobase holding the junctions, the two inserts are the same width as the arms either side of them. Realigned against the derivative, reads the reference tore into pieces run straight through: none clips at a junction, and depth holds flat across them.

Each hg38 window runs past the segment the allele takes, so the bare reference either side of the reads is what this allele leaves behind. That lane draws split alignments only, and its coverage band counts the same subset: the reads carrying a junction, stepping down as each arm runs out.

The stitching at base scale: chr3 runs out, chr10 follows, then chr12 inverted, then chr3 resumes backwards. Above, the same molecules against hg38, split alignments only, each row stopping at a junction with a connector to the piece it continues on; below, the allele's segments over the reads realigned to them at flat depth. Open this view in JBrowse ↗

Following one read across the two alignments is what a breakpoint split view does: soft clipping is shown on both sides, and a curve joins each molecule's pieces. The hg38 side carries a panel per locus the allele visits, so every connector runs between two segments that are both on screen; a dashed connector means the read passes through a segment no panel is showing.

COLO829 tumor ONT reads over one junction, twice. Against hg38 (left, split alignments only) they stop at chr3:25,359,568 with their tails clipped; realigned to the derivative (right) they cross at flat depth. The panes are at different zooms. Open in JBrowse: hg38 ↗ (large dataset — this session takes a while to load) · derivative ↗

Reproduce it end to end

scripts/build_cancer_sv_demo.sh builds everything above from public sources:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_cancer_sv_demo.sh
bash build_cancer_sv_demo.sh    # builds ./cancer_sv_build/jbrowse2
npx --yes serve cancer_sv_build/jbrowse2

It fetches the ONT COLO829 somatic SV calls and coverage and runs both sv_multihop.py steps against the tumor CRAM over HTTP. The same script builds the K562 half of the demo, which Gene fusion calls and the DNA behind them walks through.

References

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
Nattestad M, et al. Complex rearrangements and oncogene amplifications revealed by long-read DNA and RNA sequencing of a breast cancer cell line. Genome Research (2018). https://doi.org/10.1101/gr.231100.117
Nattestad M, Aboukhalil R, Chin CS, Schatz MC. Ribbon: intuitive visualization for complex genomic variation. Bioinformatics (2021). https://doi.org/10.1093/bioinformatics/btaa1080

Feedback on this tutorial is welcome: contact us.

Complex rearrangements and derivative alleles

Prerequisites#

Where the data comes from#

COLO829#

Multi-hop fusions#

Finding the chains#

Reads at the breakpoints#

Following the chain across panels#

Reconstructing the derivative allele in the browser#

A fold-back on one chromosome#

The der(3) allele, four segments across three chromosomes#

Reconstructing the allele's sequence#

Checking the reconstruction#

Reproduce it end to end#

See also#

References#