# Structural variants from Hi-C
**TL;DR:** JBrowse fetches a Hi-C matrix for every _pair_ of regions on screen.
Put a chr9 window and a chr22 window in one linear view and the space between
them fills with the contacts between the two: background in a normal karyotype,
and the Philadelphia chromosome in K562.
## Prerequisites
- nothing to install to read along: every track is a public ENCODE object served
with CORS and byte ranges
- `java`, which `juicer_tools` needs for the [scan script](#run-the-scan); the
script downloads `juicer_tools` itself
- `curl`, for the same [scan script](#run-the-scan)
## Where the data comes from
Deep in situ Hi-C for GM12878 and K562 from ENCODE, plus ENCODE's own domain,
loop and compartment calls over the same two matrices.
- GM12878 in situ Hi-C (ENCSR410MDC):
https://encode-public.s3.amazonaws.com/2021/10/28/6f0cc163-86c7-4a68-baac-65af90f5a90d/ENCFF053VBX.hic
- K562 in situ Hi-C (ENCSR545YBD):
https://encode-public.s3.amazonaws.com/2021/10/28/4d332729-3463-4782-b33c-76e4fa8ff72a/ENCFF080DPJ.hic
- GM12878 contact domains (Arrowhead):
https://encode-public.s3.amazonaws.com/2021/10/28/467750ae-7aab-47b0-a304-dc5f8dff89f7/ENCFF301CUL.bedpe.gz
- GM12878 loops (HiCCUPS):
https://encode-public.s3.amazonaws.com/2021/10/28/70e6944c-1212-45f9-855c-dbc74e9a21f5/ENCFF712NKX.bedpe.gz
- GM12878 compartment eigenvector:
https://encode-public.s3.amazonaws.com/2021/10/28/5b488af0-df49-4b9b-9feb-8ad671b7eaef/ENCFF661LPK.bigWig
- K562 compartment eigenvector:
https://encode-public.s3.amazonaws.com/2021/10/28/1180b7b2-99fd-429a-bfe1-f76cc8aa751a/ENCFF699RSL.bigWig
- pseudobulk ATAC-seq by blood lineage beside the MYC domain (10x 5k-PBMC,
SnapATAC2), rehosted for the
[scATAC pseudobulk tutorial](https://jbrowse.org/jb2-staging/docs/tutorials/scatac_pseudobulk). One bigWig
per cell type, so CD14 monocytes are
https://jbrowse.org/demos/scatac_pbmc5k/CD14_Mono.bw and naive CD4 T cells
https://jbrowse.org/demos/scatac_pbmc5k/CD4_Naive.bw
## What Hi-C measures, and what it looks like
Hi-C counts how often two stretches of the genome are found touching each other
in the nucleus. JBrowse draws the result as a triangle: the diagonal runs along
the top edge, and depth below it is genomic separation, so a bin near the top is
a pair of loci close together and a bin near the bottom is a pair far apart.
A linear browser puts the matrix in the same coordinate system as genes,
annotations and read-level tracks, which is what the rest of this page reads the
contacts against.
Two features of that picture have names, and ENCODE publishes both as separate
annotation files derived from the same matrix:
- **Contact domains** (also TADs) are the square blocks sitting on the diagonal.
Inside one, everything contacts everything; across a boundary, contact drops
sharply. ENCODE calls them with
[Arrowhead](https://github.com/aidenlab/juicer/wiki/Arrowhead) and ships a
BEDPE.
- **Loops** are individual bright dots off the diagonal: two specific points
contacting each other far more than their separation predicts, usually a pair
of convergent CTCF sites. ENCODE calls them with
[HiCCUPS](https://github.com/aidenlab/juicer/wiki/HiCCUPS), also a BEDPE.
A contact domain and the loop at its corner are the same object seen two ways:
the loop is what holds the domain together, so the block in the matrix, the
Arrowhead arc and the HiCCUPS arc all end on the domain's two corners. _MYC_ is
at the left one.
All three lanes come out of one experiment. ENCODE's Arrowhead and HiCCUPS call
sets for GM12878 are derived from the in situ matrix drawn underneath them, and
that matrix is deep enough for a 600 kb block to have visible edges.
The ATAC lane comes from outside the contact map. GM12878 is a B-lymphoblastoid
line, so B-cell accessibility is the nearest public annotation of which sequence
here is regulatory. All twelve lineages are drawn, and each is more accessible
inside the domain than outside it.
The window is a domain-and-loop pair: taking every Arrowhead domain whose two
corners carry a HiCCUPS loop and ranking by that loop's contact count puts this
one near the top, with _MYC_ at its left anchor. The
[scoring script](https://github.com/GMOD/jbrowse-components/blob/main/scripts/hic_pick_loop.py)
prints that ranking and what a candidate window contains.
Which bins the matrix is drawn in decides whether any of it is visible. JBrowse
picks the largest binsize no coarser than twice the current bp-per-pixel, which
at this width is fine enough that the triangle renders as red speckle. The
figure above sets
[`resolutionBias`](https://jbrowse.org/jb2-staging/docs/config/linearhicdisplay/#slot-resolutionbias) to `2`,
stepping two levels coarser. If a Hi-C track looks like noise, this is the first
thing to change; see
[adjusting resolution](https://jbrowse.org/jb2-staging/docs/user_guides/hic_track#adjusting-resolution).
## Two chromosomes in one view
The matrix is fetched for every pair of displayed regions. With one region on
screen that is the region against itself; open a second and JBrowse also fetches
the contacts _between_ the two and draws them in the wedge between their
triangles. Type both locations into the location box, separated by a space.
That makes a genome browser a translocation detector. Contact frequency decays
with distance along whatever molecule two loci actually sit on, so two regions
on separate chromosomes only touch at the nucleus' low background rate. If they
are fused, they are neighbours, and they contact each other constantly.
K562 is derived from a chronic myeloid leukaemia patient and carries the
Philadelphia chromosome, t(9;22)(q34;q11)
([Rowley 1973](https://doi.org/10.1038/243290a0)), which joins _BCR_ on chr22 to
_ABL1_ on chr9, the fusion imatinib targets. GM12878 is a lymphoblastoid line
with a normal karyotype. Both have deep in situ Hi-C from the same lab and
pipeline in ENCODE, so the two maps are directly comparable.
The paired triangles are the same in both panels: chr9 and chr22 each fold
normally in K562.
## Depth and normalization
Two things decide what the empty wedge in the top panel means: how deeply that
file was sequenced, and how its matrix was normalized.
**Depth.** ENCODE's GM12878 "supernatant" fraction (`ENCSR730CER`, the
alternative the script carries commented out) is much shallower over this
chromosome pair than the deep in situ file the figure uses, `ENCSR410MDC`, and a
wedge empty for want of sequencing reads the same as one empty for want of a
translocation. Run the scan as it ships and GM12878 carries more contact than
K562 across the whole chr9-chr22 block; the order inverts at the junction bin.
**Normalization.** Matrix balancing exists to divide out per-bin coverage
differences, and an amplified fusion _is_ a coverage difference. Re-run the scan
below with `NORM=INTER_SCALE` and *ABL1*×*BCR* drops off the top of the table.
Both Hi-C tracks in this demo therefore set
[`selectedNormalization`](https://jbrowse.org/jb2-staging/docs/config/linearhicdisplay/#slot-selectednormalization)
to `NONE`. Balanced matrices are what to read domains and loops with;
rearrangements want the raw counts.
The control gets its own ranked list, which the scan prints below the case's.
The bin at its head is hot in GM12878, present in K562, and not a rearrangement.
## Run the scan
The figure shows the translocation; finding it is a dump and a sort. The dump is
one command per file, and it is where the normalization decision above becomes a
command-line argument:
```bash
# `observed NONE` asks for raw counts. NONE rather than a balanced vector for
# the reason above: balancing divides out per-bin coverage differences, and an
# amplified fusion is one, so a balanced dump removes what the scan looks for.
# A balanced vector is stored only at the coarser bin sizes, so asking for one
# at a fine resolution comes back as an empty file rather than an error.
# BP asks for base-pair bins rather than restriction fragments, and the number
# after it is the bin size.
# -Xmx4g because a whole chromosome pair does not fit in the default heap.
java -Xmx4g -jar juicer_tools.jar dump observed NONE \
sample.hic chr9 chr22 BP 250000 sample.chr9_chr22.txt
```
Three columns come back: bin1 start, bin2 start, contact count. Run it over the
control too, then rank the sample's bins and read the control's value for each.
[`scan_hic_translocation.sh`](https://github.com/GMOD/jbrowse-components/blob/main/scripts/scan_hic_translocation.sh)
does exactly that:
```bash
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/scan_hic_translocation.sh
bash scan_hic_translocation.sh
```
Underneath it, that is one dump per file and a sort:
```bash
# one inter-chromosomal block at one bin size. An empty output file means this
# .hic stores neither the pair nor a KR vector at that resolution; read its
# footer.
java -Xmx4g -jar juicer_tools.jar dump observed KR \
case.hic chr9 chr22 BP 25000 case.txt
# bin1, bin2, contacts. awk and not `| head`, which closes the pipe and kills
# sort with SIGPIPE mid-table under `set -o pipefail`.
sort -k3,3 -rn case.txt | awk 'NR <= 10'
```
`CASE`, `CTRL`, `CHR1`, `CHR2`, `RES` and `NORM` are all overridable, so the
same scan applies to any two `.hic` files that hold inter-chromosomal blocks.
The top row it prints pairs _ABL1_ intron 1 with the 5' end of _BCR_, the pair
of bins the canonical CML fusion joins. Drop `RES` to `10000` and the top row
lands on the junction itself, _ABL1_ intron 1 against the _BCR_ major breakpoint
cluster region.
Further down its list the scan reports a second partner for chr9 elsewhere on
chr22, well clear of the control: K562's karyotype is complex, and the ranking
is a list of candidates to open.
Purpose-built callers do this genome-wide with a trained model:
[EagleC](https://github.com/XiaoTaoWang/EagleC),
[hic_breakfinder](https://github.com/dixonlab/hic_breakfinder) and
[HiNT](https://github.com/parklab/HiNT) are the usual ones. Their output is
BEDPE, which loads here as a
[paired-arc track](https://jbrowse.org/jb2-staging/docs/config_guides/hic_track#loops-and-interactions-as-arcs)
next to the matrix it was called from.
## A and B compartments
Above the domains and loops of the first figure, the matrix separates into two
interleaved sets of regions that preferentially contact their own kind: the
gene-rich, active A compartment and the inactive B compartment. ENCODE publishes
that call for every Hi-C experiment as a
[compartment eigenvector and a set of subcompartment classes](https://jbrowse.org/jb2-staging/docs/user_guides/hic_track#compartments-and-subcompartments),
both derived from the matrix already loaded.
The band over _TCF4_ is in the B compartment in GM12878 and the A compartment in
K562, and the sequence either side of it, from the same two files and the same
pipeline, agrees.
Two settings in the figure:
- The eigenvector tracks are pinned to one shared scale, so the two lanes are
read against the same axis.
- An eigenvector names the A compartment only up to a sign, so which sign means
active is a property of the file. It is read off the gene track, since A is
the gene-rich compartment by definition.
The
[user guide section](https://jbrowse.org/jb2-staging/docs/user_guides/hic_track#compartments-and-subcompartments)
covers both, along with why the subcompartment class numbers cannot be compared
between files on their own.
## Configuring the Hi-C tracks
The `.hic` files are 20 GB and 55 GB, and only the bins on screen are requested.
```json addtrack
{
"type": "HicTrack",
"trackId": "hic_k562_insitu",
"name": "K562 in situ Hi-C (ENCODE ENCSR545YBD)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "HicAdapter",
"uri": "https://encode-public.s3.amazonaws.com/2021/10/28/4d332729-3463-4782-b33c-76e4fa8ff72a/ENCFF080DPJ.hic"
},
"displayDefaults": {
"selectedNormalization": "NONE"
}
}
```
Use ENCODE's direct S3 URLs. hic-straw's range reader cannot follow the
cross-origin redirect the portal's `@@download` links return, and the `.hic`
comes back as a 403. The S3 URL for any ENCODE file is in its metadata under
`cloud_metadata.url`.
The loop and domain BEDPEs each need one extra slot:
```json addtrack
{
"type": "FeatureTrack",
"trackId": "hic_gm12878_domains",
"name": "GM12878 contact domains (Arrowhead)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "BedpeAdapter",
"uri": "https://encode-public.s3.amazonaws.com/2021/10/28/467750ae-7aab-47b0-a304-dc5f8dff89f7/ENCFF301CUL.bedpe.gz"
}
}
```
A contact domain is a `FeatureTrack`. Arrowhead writes each domain with both
BEDPE mates set to the same interval, so read as plain features the file gives
one box per domain, with nested domains stacking into rows. Loops, whose two
mates are different places, are the paired-arc case. See the
[Hi-C track config guide](https://jbrowse.org/jb2-staging/docs/config_guides/hic_track#loops-and-interactions-as-arcs).
To color or filter either track by a column, set
[`columnNames`](https://jbrowse.org/jb2-staging/docs/config/bedpeadapter/#slot-columnnames) explicitly. The
adapter otherwise reads names off the file's own header, and juicer writes its
version banner _after_ the defline, so every column past the tenth reads back as
`undefined` and a jexl expression on one silently evaluates against nothing.
HiCCUPS writes 24 columns and Arrowhead 16; only the first ten of either are
positional.
Both callers also leave the standard BEDPE `name` and `score` columns at `.` and
put what they rank by further along: HiCCUPS' is `observed`, and Arrowhead's is
a _second_ column called `score`. Values past the tenth column arrive as
strings, so compare them with `>` and `<`, which coerce, rather than with a jexl
`==`.
## See also
- [](https://jbrowse.org/jb2-staging/docs/user_guides/hic_track)
- [](https://jbrowse.org/jb2-staging/docs/config_guides/hic_track)
- [](https://jbrowse.org/jb2-staging/docs/tutorials/chromhmm)
- [](https://jbrowse.org/jb2-staging/docs/tutorials/cancer_sv)
- [](https://jbrowse.org/jb2-staging/docs/user_guides/sv_visualization)
- [HiGlass](https://higlass.io/)