# 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
- [](https://jbrowse.org/jb2-staging/docs/cli)
- [samtools](http://www.htslib.org/) (v1.21 or later)
- [minimap2](https://github.com/lh3/minimap2)
- `bedGraphToBigWig` from the
[UCSC utilities](https://hgdownload.soe.ucsc.edu/admin/exe/)
- `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](https://nodejs.org/). `sv_multihop.py` is one file:
```bash
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](https://doi.org/10.1016/j.xgen.2022.100139)):
- 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:
```bash
# 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](https://doi.org/10.1101/gr.231100.117)): 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`](https://github.com/GMOD/jbrowse-components/blob/main/scripts/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:
```bash
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.
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.
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.
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](https://doi.org/10.1093/bioinformatics/btaa1080))
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
[](https://jbrowse.org/jb2-staging/docs/user_guides/sv_visualization#what-the-reconstruction-needs).
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.
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.
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.
```bash
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:
```bash
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.
## 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.
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.
## Reproduce it end to end
[`scripts/build_cancer_sv_demo.sh`](https://github.com/GMOD/jbrowse-components/blob/main/scripts/build_cancer_sv_demo.sh)
builds everything above from public sources:
```bash
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 [](https://jbrowse.org/jb2-staging/docs/tutorials/k562_fusions) walks through.
## See also
- [](https://jbrowse.org/jb2-staging/docs/tutorials/k562_fusions)
- [](https://jbrowse.org/jb2-staging/docs/tutorials/sv_callset_review)
- [](https://jbrowse.org/jb2-staging/docs/tutorials/hic_structural_variants)
- [](https://jbrowse.org/jb2-staging/docs/user_guides/sv_visualization)
- [](https://jbrowse.org/jb2-staging/docs/user_guides/sv_inspector_view)
- [](https://jbrowse.org/jb2-staging/docs/user_guides/linear_synteny_view)
- [](https://jbrowse.org/jb2-staging/docs/tutorials/sv_visualization_cgiab)
## 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