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, whichjuicer_toolsneeds for the scan script; the script downloadsjuicer_toolsitselfcurl, for the same scan script
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. 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 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, 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 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 to 2,
stepping two levels coarser. If a Hi-C track looks like noise, this is the first
thing to change; see
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), 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
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:
# `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.txtThree 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
does exactly that:
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/scan_hic_translocation.sh
bash scan_hic_translocation.shUnderneath it, that is one dump per file and a sort:
# 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, hic_breakfinder and HiNT are the usual ones. Their output is BEDPE, which loads here as a paired-arc track 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, 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 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.
{
"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"
}
}jbrowse add-track https://encode-public.s3.amazonaws.com/2021/10/28/4d332729-3463-4782-b33c-76e4fa8ff72a/ENCFF080DPJ.hic \
--trackId hic_k562_insitu \
--name "K562 in situ Hi-C (ENCODE ENCSR545YBD)" \
--assemblyNames hg38 \
--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:
{
"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"
}
}jbrowse add-track-json '{
"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.
To color or filter either track by a column, set
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
- Hi-C track
- Hi-C track
- ChromHMM chromatin states
- Complex rearrangements and derivative alleles
- SV visualization
- HiGlass
Feedback on this tutorial is welcome: contact us.