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
- a JBrowse to paste the tracks into (Web or Desktop); every file here is a URL, so Desktop needs nothing hosted
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 touch in the nucleus. JBrowse draws the result as a triangle: the diagonal runs along the top edge, and depth below it is genomic separation. Two features of that picture have names, and ENCODE publishes both as annotation files derived from the 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.
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. The ATAC lane comes from outside the contact map: GM12878 is a B-lymphoblastoid line, and all twelve lineages are more accessible inside the domain than outside it.
The window was chosen by taking every Arrowhead domain whose two corners carry a HiCCUPS loop and ranking by that loop's contact count; the scoring script prints that ranking.
JBrowse picks the largest binsize no coarser than twice the current
bp-per-pixel, which at this width renders the triangle as red speckle. The
figure sets
resolutionBias to 2,
stepping two levels coarser. If a Hi-C track looks like noise, change this
first; see
adjusting resolution.
Two chromosomes in one view
The matrix is fetched for every pair of displayed regions. Open a second region and JBrowse also fetches the contacts between the two, drawn in the wedge between their triangles. Type both locations into the location box, separated by a space.
Contact frequency decays with distance along whatever molecule two loci sit on, so two regions on separate chromosomes only touch at background rate. If they are fused, they contact each other constantly.
K562 carries the Philadelphia chromosome, t(9;22)(q34;q11) (Rowley 1973), joining BCR on chr22 to ABL1 on chr9. GM12878 has a normal karyotype. Both have deep in situ Hi-C from the same ENCODE lab and pipeline.
The paired triangles are the same in both panels: chr9 and chr22 each fold normally in K562.
Depth and normalization
Depth. ENCODE's GM12878 "supernatant" fraction (ENCSR730CER, which the
script carries commented out) is much shallower over this chromosome pair than
the in situ file the figure uses, ENCSR410MDC, and a wedge empty for want of
reads looks the same as one empty for want of a translocation. As shipped, the
scan finds GM12878 carrying more contact than K562 across the whole chr9-chr22
block, with the order inverting at the junction bin.
Normalization. Matrix balancing divides out per-bin coverage differences,
and an amplified fusion is one. Re-run the scan with NORM=INTER_SCALE and
ABL1×BCR drops off the top of the table. Both Hi-C tracks here set
selectedNormalization
to NONE: balanced matrices for domains and loops, raw counts for
rearrangements.
The scan prints the control's own ranked list below the case's. The bin at its head is hot in GM12878, present in K562, and not a rearrangement.
Run the scan
Finding the translocation is a dump and a sort, one dump per file:
# `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. 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 pairs ABL1 intron 1 with the 5' end of BCR. Drop RES to
10000 and it lands on the junction itself, ABL1 intron 1 against the BCR
major breakpoint cluster region. Further down, a second chr9 partner elsewhere
on chr22 sits well clear of the control; 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 domains and loops, the matrix separates into two interleaved sets of regions that contact their own kind: the gene-rich, active A compartment and the inactive B compartment. ENCODE publishes that call for every experiment as a compartment eigenvector and a set of subcompartment classes.
The band over TCF4 is in the B compartment in GM12878 and the A compartment in K562, and the sequence either side of it agrees. Two settings in the figure:
- The eigenvector tracks are pinned to one shared scale
- An eigenvector names the A compartment only up to a sign, so which sign is active is read off the gene track, A being the gene-rich compartment
The user guide section covers both, and why subcompartment class numbers cannot be compared between files.
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"}'In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:
{
"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:
{
"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"
}
}'In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:
{
"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"
}
}Arrowhead writes each domain with both BEDPE mates set to the same interval, so
as a FeatureTrack the file gives one box per domain, nested domains stacking
into rows. Loops, whose mates differ, are the paired-arc case. See the
Hi-C track config guide.
To color or filter either track by a column, set
columnNames explicitly. Juicer
writes its version banner after the defline, so names read off the header make
every column past the tenth undefined, and a jexl expression on one silently
evaluates against nothing. HiCCUPS writes 24 columns and Arrowhead 16.
Both callers leave name and score at . and put what they rank by further
along: HiCCUPS' is observed, Arrowhead's a second column called score.
Values past the tenth column arrive as strings, so compare with > and <,
which coerce, rather than ==.
See also
- User guide: Hi-C track
- Config guide: Hi-C track
- ChromHMM chromatin states
- Complex rearrangements and derivative alleles
- SV visualization
- HiGlass
Feedback on this tutorial is welcome: contact us.