Pangenome (HPRC) part 2, who carries what
2 videos in this tutorial
Part 1 ends with one allele and the single haplotype the graph credits it to. Here we use the rest of HPRC release 2 to find everyone who carries that sequence, on the same GRCh38 axis. The tracks show where the graph varies and by how much, what each alternative is, and which of the 464 haplotypes walk it. Every track is a URL that JBrowse reads a window at a time.
Prerequisites
- the session from part 1: hg38 with its genes, and the rGFA segments track loaded on it
- the GraphGenomeView plugin,
for the tracks that use
MinigraphBubbleAdapterandRgfaTabixAdapter; both callsets are a URL you can paste and need no plugin
Where the data comes from
The data is HPRC release 2. JBrowse reads its wave callset and snarl-level carriage file directly from S3, and reads the bubble and allele projections from our host.
- the decomposed variant callset, 464 haplotypes, read straight off S3: https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/minigraph-cactus/v2.1/hprc-v2.1-mc-grch38/hprc-v2.1-mc-grch38.wave.vcf.gz
- the undecomposed, snarl-level carriage file, 462 haplotypes: https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/minigraph-cactus/v2.1/hprc-v2.1-mc-grch38/hprc-v2.1-mc-grch38.pgbi.vcf.gz
- the release's all-vs-GRCh38 alignment, sliced for the inversion synteny figure: https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/impg/pafs/hprc465vsgrch38.aln.paf.gz
- the CAT gene annotation index, one GFF3 per haplotype, which the inversion figure's lanes carry: https://raw.githubusercontent.com/human-pangenomics/hprc_intermediate_assembly/main/data_tables/annotation/cat/cat_genes_hprc_r2_v1.3.index.csv
- our bubble and allele projections, with the exact build recorded beside them: https://jbrowse.org/demos/hprc/README.txt
The bubble track
A bubble is a place where haplotypes diverge and rejoin. The bubble track reads one file and reports where the graph varies and by how much:
{
"type": "FeatureTrack",
"trackId": "hprc_minigraph_bubbles",
"name": "HPRC release 2 bubbles",
"assemblyNames": ["hg38"],
"adapter": {
"type": "MinigraphBubbleAdapter",
"uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.bubbles.bed.gz",
"assemblyNameToPanSN": { "hg38": "GRCh38" }
}
}jbrowse add-track-json '{
"type": "FeatureTrack",
"trackId": "hprc_minigraph_bubbles",
"name": "HPRC release 2 bubbles",
"assemblyNames": ["hg38"],
"adapter": {
"type": "MinigraphBubbleAdapter",
"uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.bubbles.bed.gz",
"assemblyNameToPanSN": { "hg38": "GRCh38" }
}
}'In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add pangenome graph track, and fill in:
- File type:
Minigraph bubbles (tabix BED) - Path:
https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.bubbles.bed.gz - Sample name in the graph:
GRCh38 - Track name:
HPRC release 2 bubbles - Assembly:
hg38
The MinigraphBubbleAdapter labels each bubble with its shortest and longest
allele, and one bubble in the HLA class II window spans tens of kilobases
depending on the haplotype.1 HPRC publishes no bubble file, so we
built this one too, with gfatools bubble.
A whole chromosome, one node per bubble
The segments track part 1 loads draws one node per segment. A window past a few hundred kilobases holds more nodes than any layout can place.
The bubble file also gives a coarser level of detail. We collapse each bubble to
a single node and keep the invariant reference between bubbles as backbone, and
the same graph then fits on a screen. We host that tier beside the fine index.
build_bubble_tier.sh
builds it from the bubble file in one pass, at a threshold of 10,000. A tier is
a prefix, so you choose a level of detail by choosing a file:
{
"type": "FeatureTrack",
"trackId": "hprc_tier",
"name": "HPRC release 2 graph: bubble tier (one node per bubble)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "RgfaTabixAdapter",
"uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.tier10000",
"assemblyNameToPanSN": { "hg38": "GRCh38" }
}
}jbrowse add-track-json '{
"type": "FeatureTrack",
"trackId": "hprc_tier",
"name": "HPRC release 2 graph: bubble tier (one node per bubble)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "RgfaTabixAdapter",
"uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.tier10000",
"assemblyNameToPanSN": { "hg38": "GRCh38" }
}
}'In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add pangenome graph track, and fill in:
- File type:
rGFA segments (tabix BED pair) - Path:
https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.tier10000.segs.bed.gz - Sample name in the graph:
GRCh38 - Track name:
HPRC release 2 graph: bubble tier (one node per bubble) - Assembly:
hg38
The same bubble file also plots as a curve of where the graph varies and by how
much. MinigraphBubbleAdapter already reports each bubble's segment count as
its score, so only the track type changes:
{
"type": "QuantitativeTrack",
"trackId": "hprc_bubble_score",
"name": "HPRC release 2 graph: variability (segments per bubble)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "MinigraphBubbleAdapter",
"uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.bubbles.bed.gz",
"assemblyNameToPanSN": { "hg38": "GRCh38" }
}
}jbrowse add-track-json '{
"type": "QuantitativeTrack",
"trackId": "hprc_bubble_score",
"name": "HPRC release 2 graph: variability (segments per bubble)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "MinigraphBubbleAdapter",
"uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.bubbles.bed.gz",
"assemblyNameToPanSN": { "hg38": "GRCh38" }
}
}'In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add pangenome graph track, and fill in:
- File type:
Minigraph bubbles (tabix BED) - Path:
https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.bubbles.bed.gz - Sample name in the graph:
GRCh38 - Track name:
HPRC release 2 graph: variability (segments per bubble) - Assembly:
hg38
Use the two granularities together: the tier to find an event, and the fine index to open it. The node's menu moves you from one to the other. On a tier node, Open in hg38 puts the linear view on the bubble's span. A drag across that span offers Graph genome view (this selection) as a submenu naming both tracks, and the fine cut opens as a second pane under the tier's. On a whole chromosome a bubble is narrower than a pixel, so find it in the tier lane or the variability curve first and cut the tier around it.2
The allele inventory
The bubbles report where the graph varies. A third hosted file reports what the variation is. It has one row per allele in the graph, anchored on GRCh38, and we derived it from the two indexes part 1 loads.
{
"type": "AlignmentsTrack",
"trackId": "hprc_minigraph_alleles",
"name": "HPRC release 2 graph: allele inventory",
"assemblyNames": ["hg38"],
"adapter": {
"type": "BedTabixAdapter",
"uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.alleles.bed.gz"
}
}jbrowse add-track https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.alleles.bed.gz \
--trackId hprc_minigraph_alleles \
--name "HPRC release 2 graph: allele inventory" \
--assemblyNames hg38 \
--trackType AlignmentsTrackIn 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": "AlignmentsTrack",
"trackId": "hprc_minigraph_alleles",
"name": "HPRC release 2 graph: allele inventory",
"assemblyNames": ["hg38"],
"adapter": {
"type": "BedTabixAdapter",
"uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.alleles.bed.gz"
}
}An AlignmentsTrack over a BED draws the allele sizes. Each row carries a
CIGAR against the reference span it replaces (2062M63348I). The alignments
display draws any row that has a CIGAR, so the alleles pack into rows and each
insertion draws at its real magnitude. At the CFH cluster on chr1, one of those
rows is the 84,684 bp deletion between CFHR3 and CFHR1, which part 1
draws on the same coordinates as a graph edge.
The whole graph holds a few hundred thousand alleles, so a wide window is dense.
The
graph genome view guide
explains the columns and how the walk derives them. It also gives two filters
that make a lane this size readable: jexl:abs(feature.delta)>10000 for size,
and jexl:feature.nested==0 before reading lengths in bulk.
The variant callset
The wave.vcf.gz ships its index beside it, so JBrowse reads only the slice you
are viewing out of the 2.3 GB file. The adapter needs only the S3 URL, on a
VariantTrack with the multi-sample display:
{
"type": "VariantTrack",
"trackId": "hprc2_wave_grch38",
"name": "HPRC2 pangenome",
"assemblyNames": ["hg38"],
"adapter": {
"type": "VcfTabixAdapter",
"uri": "https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/minigraph-cactus/v2.1/hprc-v2.1-mc-grch38/hprc-v2.1-mc-grch38.wave.vcf.gz"
},
"displays": [
{
"type": "LinearMultiSampleVariantDisplay",
"renderingMode": "phased"
}
]
}jbrowse add-track-json '{
"type": "VariantTrack",
"trackId": "hprc2_wave_grch38",
"name": "HPRC2 pangenome",
"assemblyNames": ["hg38"],
"adapter": {
"type": "VcfTabixAdapter",
"uri": "https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/minigraph-cactus/v2.1/hprc-v2.1-mc-grch38/hprc-v2.1-mc-grch38.wave.vcf.gz"
},
"displays": [
{
"type": "LinearMultiSampleVariantDisplay",
"renderingMode": "phased"
}
]
}'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": "VariantTrack",
"trackId": "hprc2_wave_grch38",
"name": "HPRC2 pangenome",
"assemblyNames": ["hg38"],
"adapter": {
"type": "VcfTabixAdapter",
"uri": "https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/minigraph-cactus/v2.1/hprc-v2.1-mc-grch38/hprc-v2.1-mc-grch38.wave.vcf.gz"
},
"displays": [
{
"type": "LinearMultiSampleVariantDisplay",
"renderingMode": "phased"
}
]
}renderingMode: "phased" splits each phased sample column into its two
haplotypes, drawing one row per haplotype. Co-inherited blocks are visible only
in that form.
The VCF is fully decomposed, so chr6:32,450,000-32,650,000 (the window in the
figure below) holds over fourteen thousand records. Most are SNPs and the rest
are small indels. This file already contains the structural tier. Add this
filter from Edit filters to select it:
jexl:feature.INFO.LV[0]==0 && alleleLength(feature)>=50The same window then drops to a couple of hundred sites.
The filter needs both halves:
alleleLengthis the longest allele the record describes. A filter onend - startwould keep only deletions, since an insertion consumes no reference.LVis the record's level in vg's snarl tree, andLV==0keeps the top-level sites. This file writes a nested child as a separate record beside its parent, withPSnaming that parent. WithoutLV==0the panel paints some events twice, at two positions.
With the filter on, few enough alleles remain to draw each at its genomic position, lined up with the genes above. Clustering → Cluster rows by genotype... → Run clustering in the track menu reorders the rows by genotype similarity and draws a dendrogram beside them. Comparing the graph with the callset puts that matrix beside the graph that the same alleles came from.
Frequency is in the file. AC, AF, AN and NS are on every record, so
jexl:feature.INFO.AF[0]>0.05 selects the common alleles without clustering
anything.
The display widens each insertion cell to a marker sized by the inserted bp, in
that haplotype's genotype color
(showInsertionGlyphs).
Only haplotypes carrying the allele widen.
Carriage at the graph's own granularity
The callset above is decomposed, so one graph bubble becomes many records, and a
column is a primitive variant. Release 2 also publishes the undecomposed form,
one record per snarl, as the pgbi.vcf.gz beside each build. Read that file
to find who carries a given bubble, because its rows are the graph's alleles.
{
"type": "VariantTrack",
"trackId": "hprc2_pgbi_grch38",
"name": "HPRC2 pangenome carriage (snarl-level, 462 haplotypes)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "VcfTabixAdapter",
"uri": "https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/minigraph-cactus/v2.1/hprc-v2.1-mc-grch38/hprc-v2.1-mc-grch38.pgbi.vcf.gz"
},
"displays": [
{ "type": "LinearMultiSampleVariantDisplay", "renderingMode": "phased" }
]
}jbrowse add-track-json '{
"type": "VariantTrack",
"trackId": "hprc2_pgbi_grch38",
"name": "HPRC2 pangenome carriage (snarl-level, 462 haplotypes)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "VcfTabixAdapter",
"uri": "https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/minigraph-cactus/v2.1/hprc-v2.1-mc-grch38/hprc-v2.1-mc-grch38.pgbi.vcf.gz"
},
"displays": [
{ "type": "LinearMultiSampleVariantDisplay", "renderingMode": "phased" }
]
}'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": "VariantTrack",
"trackId": "hprc2_pgbi_grch38",
"name": "HPRC2 pangenome carriage (snarl-level, 462 haplotypes)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "VcfTabixAdapter",
"uri": "https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/minigraph-cactus/v2.1/hprc-v2.1-mc-grch38/hprc-v2.1-mc-grch38.pgbi.vcf.gz"
},
"displays": [
{ "type": "LinearMultiSampleVariantDisplay", "renderingMode": "phased" }
]
}The snarl file ships with an index, so JBrowse downloads only the slice in view.
The C4 window loads in a couple of seconds over HTTP. The file carries 231
sample columns, without CHM13, so phased mode draws 462 rows and AN tops out
at 462 too.
This file adds the AT field, which records each allele as the traversal it
takes through the graph. A pggb VCF carries the same AT field. The wave file
drops it, and its header shows the bcftools annotate -x INFO/AT command that
did so. The same LV==0 filter cuts this lane to
top-level sites, the tier that the bubble track holds. A
record matches a bubble by interval.3
With this lane loaded, two compatible readings sit over one coordinate. The allele inventory gives the haplotype the graph credits an allele to. A genotype column here gives the haplotypes that walk the allele.
Comparing the graph with the callset
The graph and the callset describe the same variation at two resolutions. minigraph records structural variation (roughly >50 bp) and collapses everything smaller, so SNPs are absent from the graph even though every one is in the VCF. Filter the callset to that same tier, and the two describe the same events from opposite ends. The graph records an allele and its length, and collapsing the smaller variation makes the allele easy to find. The callset kept the samples, so it records who carries each allele.
The events line up across the two. Mark an interval in the linear view, and the mark crosses the genes, the segments lane and the genotype matrix in one column. The reference-position ramp gives the graph's backbone at that position the same hue as the segments above it.
Inversions
Insertions are nodes and deletions are edges. An inversion is the same reference
sequence, walked backwards. gfatools bubble sets a column when a bubble's
paths disagree about orientation, and the adapter exposes that column as an
inversion boolean. Edit filters on the bubble track cuts the lane to
inversions:
jexl:feature.inversionbuild_hprc_inversion_synteny.sh
reduces the 1q21.1 bubble to two haplotype rows from HPRC's published
all-vs-GRCh38 PAF, one carrying the inversion and one not. Each row shows the
CAT gene annotation for that haplotype. The boxed pair on each row is the same
two genes, PPIAL4F and PPIAL4E. On the carrier PPIAL4F comes first, and on
the non-carrier PPIAL4E does. The hg38 row between them agrees with the
non-carrier.
Reproduce it end to end
Preparing your own graph builds the
bubble file, its coarse tier and the
allele inventory, and its commands run on any rGFA.
Pointed at hprc-v2.1-mc-grch38.sv.gfa.gz, it writes the files we host.
README.txt beside those files
records their provenance: source, size, exact commands and build date.
We do not rebuild carriage here, because HPRC publishes it as a file, tabix-indexed like the callset. Rebuilding it takes a 464-assembly download and a mapping run. The run makes this one call per sample:
# --call asks, at every bubble, which path this assembly takes through the graph
minigraph -cxasm --call -t"$(nproc)" graph.gfa assembly.fa > sample.call.bedbuild_minigraph_paths.sh
wraps that call in the per-sample loop and the join, and
the same page walks
through the script.
A separate script builds the inversion figure from release 2's published all-vs-GRCh38 PAF:
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_hprc_inversion_synteny.sh
bash build_hprc_inversion_synteny.sh # writes ./hprc_inversion_synteny_build/build_hprc_inversion_synteny.sh
classifies each haplotype at that bubble. It keeps the haplotypes whose
alignments reverse the block while the sequence on either side stays forward. It
prints the split it finds, then slices out one haplotype of each kind.
See also
- Pangenome (HPRC) part 1, reading the graph
- Pangenome (HPRC) part 3, every haplotype in its own coordinates
- Pangenome (preparing your own graph)
- Graph genome view
- Multi-sample variant display
References
- HPRC release 2, the release this page opens: the Minigraph-Cactus graph and the two callsets deconstructed from it.
- Li H.
The rGFA format and
gfatools, whose
bubblesubcommand calls the bubbles this page reads.
Notes
-
Each bubble's description also carries a path count, which counts the routes through the bubble and not the haplotypes observed. gfatools saturates the count at
2147483647, and the track describes those bubbles as having more paths than gfatools counts. ↩ -
The view refuses a cut over 5 Mb. That width limit stands in for node count, which works at segment granularity and fails on a tier. A
GraphGenomeViewpointed at a tier therefore raisesmaxRegionBpto the span it draws. The node limit still applies:maxGraphNodescounts the nodes that actually came back. ↩ -
IDandATname base-level integer nodes (>161001867>161004536), wheresv.gfausessNNNNNsegment ids. ↩
Feedback on this tutorial is welcome: contact us.