Pangenome (HPRC) part 2: haplotypes against each other

JBrowse v5 beta

The tutorials target the JBrowse v5 beta, and the v4.3.0 release on the download page lacks some of what they show. To install the beta of JBrowse Web, run npm install -g @jbrowse/cli@next, then jbrowse create jbrowse2 --branch v5.0.0-beta.11. Desktop beta builds are coming soon.

We draw human haplotypes from the Human Pangenome Reference Consortium's release 2 side by side, one lane per haplotype drawn in the coordinates of its own assembly, and use the pangenome graph to show where neighbouring lanes match. With that view we:

at CFH, find haplotypes missing two genes
at C4, find two haplotypes sharing a copy GRCh38 lacks
at amylase, count gene copies and check them against Yilmaz et al. (2024)
at 1q21.1, tell an inversion the graph flags from an inverted paralog

Experimental

The graph view is a beta plugin, and the HPRC page lives on staging.genomes.jbrowse.org until JBrowse 5 ships. We welcome your feedback.

Prerequisites

for Reproduce it end to end: samtools built with libcurl, minimap2, Node.js for npx, python3, and the JBrowse CLI for make-pif

Where the data comes from

Show full URLs

HPRC release 2.

the graph as a gbz-base database: s3-us-west-2.amazonaws.com/…/hprc-v2.1-mc-grch38.gbz.dbhttps://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.gbz.db
our index naming its haplotypes: jbrowse.org/…/hprc-v2.1-mc-grch38.haplotype-index.f3.dbhttps://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.haplotype-index.f3.db
the assemblies: raw.githubusercontent.com/…/assemblies_release2_v1.0.index.csvhttps://raw.githubusercontent.com/human-pangenomics/hprc_intermediate_assembly/main/data_tables/assemblies_release2_v1.0.index.csv
the CAT gene annotation index, one GFF3 per haplotype: raw.githubusercontent.com/…/cat_genes_hprc_r2_v1.3.index.csvhttps://raw.githubusercontent.com/human-pangenomics/hprc_intermediate_assembly/main/data_tables/annotation/cat/cat_genes_hprc_r2_v1.3.index.csv
the release's all-vs-GRCh38 alignment, sliced for the inversion figure: s3-us-west-2.amazonaws.com/…/hprc465vsgrch38.aln.paf.gzhttps://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/release2/impg/pafs/hprc465vsgrch38.aln.paf.gz
our bubble projections of the graph, with the exact build recorded beside them: jbrowse.org/…/README.txthttps://jbrowse.org/demos/hprc/README.txt

Configuring the haplotype lanes track

One track draws the lanes, reading the release's gbz-base database (the graph stored as SQLite) directly. The HPRC page's haplotypes launch adds it, and this config is the one to adapt for your own graph, whose database and haplotype index Hosting your own graph builds. Each lane is one haplotype's walk, its route through the graph. List an assembly per haplotype in assemblyNames; the track finds each walk through the assembly's sample#haplotype alias. The HPRC config declares that alias for every release 2 haplotype (HG00097.1 is aliased HG00097#1), and one assembly entry does the same for a haplotype of your own. GRCh38 has no such alias, so assemblyNameToPanSN names its walk:

Config file

Goes in the tracks array of config.json. See Tracks.

{
  "type": "GraphTrack",
  "trackId": "hprc_v2_1_gbz_lanes",
  "name": "HPRC release 2 haplotypes vs GRCh38, read from the graph (gbz-base)",
  "assemblyNames": ["hg38", "HG00097.1", "HG00099.1"],
  "adapter": {
    "type": "GbzBaseSyntenyAdapter",
    "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.gbz.db",
    "haplotypeIndexLocation": {
      "uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.haplotype-index.f3.db"
    },
    "assemblyNames": ["hg38"],
    "assemblyNameToPanSN": { "hg38": "GRCh38#0" },
    "context": 1000
  },
  "displays": [
    { "type": "MultiWaySyntenyDisplay", "height": 600 },
    { "type": "LinearGraphDisplay" }
  ]
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "GraphTrack",
  "trackId": "hprc_v2_1_gbz_lanes",
  "name": "HPRC release 2 haplotypes vs GRCh38, read from the graph (gbz-base)",
  "assemblyNames": ["hg38", "HG00097.1", "HG00099.1"],
  "adapter": {
    "type": "GbzBaseSyntenyAdapter",
    "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.gbz.db",
    "haplotypeIndexLocation": {
      "uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.haplotype-index.f3.db"
    },
    "assemblyNames": ["hg38"],
    "assemblyNameToPanSN": { "hg38": "GRCh38#0" },
    "context": 1000
  },
  "displays": [
    { "type": "MultiWaySyntenyDisplay", "height": 600 },
    { "type": "LinearGraphDisplay" }
  ]
}'

JBrowse Desktop

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": "GraphTrack",
  "trackId": "hprc_v2_1_gbz_lanes",
  "name": "HPRC release 2 haplotypes vs GRCh38, read from the graph (gbz-base)",
  "assemblyNames": ["hg38", "HG00097.1", "HG00099.1"],
  "adapter": {
    "type": "GbzBaseSyntenyAdapter",
    "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.gbz.db",
    "haplotypeIndexLocation": {
      "uri": "https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.haplotype-index.f3.db"
    },
    "assemblyNames": ["hg38"],
    "assemblyNameToPanSN": { "hg38": "GRCh38#0" },
    "context": 1000
  },
  "displays": [
    { "type": "MultiWaySyntenyDisplay", "height": 600 },
    { "type": "LinearGraphDisplay" }
  ]
}

CFH: haplotypes missing two genes

Open the HPRC page and press haplotypes on the CFH / CFHR row, the complement factor H gene cluster. A lane missing CFHR3 and CFHR1 (Hughes et al. 2006), such as HG00253's second haplotype, runs from CFH straight on to CFHR4.

The CFH cluster from the haplotypes launch: RefSeq genes over one lane per structural configuration, each lane a haplotype's contig read from the graph, under its CAT genes. A lane that lacks CFHR3 and CFHR1 leaves that stretch of its neighbour unmatched in the band between them.

C4: haplotypes with a copy GRCh38 lacks

Press haplotypes on the C4A / C4B row, the complement component 4 genes. In the track menu open Lanes → Choose lanes..., press Untick shown, tick HG01978.2 and HG02004.2, and press Draw these lanes.

C4 from the HPRC page's haplotypes launch with two lanes chosen, HG01978.2 and HG02004.2, each with three copies of the C4-CYP21-TNX module, under the RefSeq genes. The band from GRCh38 leaves the third module unmatched, and the band between the two haplotypes matches it off the nodes both paths share.

Amylase: counting AMY1 copies per haplotype

Press haplotypes on the AMY1 row, the salivary amylase gene copies, and choose HG01361.1, HG00133.2, HG00133.1, NA18608.2 and HG00232.1. A lane with more copies spans more of its contig in the same width, and its label gives that span as a multiple of the window.

The amylase locus from the HPRC page's haplotypes launch with five lanes chosen, one per amylase structure, under the RefSeq genes. Each lane is drawn on the haplotype's contig under its CAT genes; a lane longer than the window gives its span as a multiple in its label, and each band draws the extra copies of the longer lane as a gap.

To read the lengths, take Display types → Graph, enter the five names in Settings → Haplotypes, then pick Layout → Walk rows and Color → Uniform to draw each haplotype's route as a bar.

The five haplotypes' walks across the amylase array in walk rows, longest first, under GRCh38's bar, each boxed with its own CAT genes so its AMY1 copies can be counted on the bar. Blue is on GRCh38's path through the graph and purple off it. Each readout gives the walk's length and its excess over GRCh38.

HG00133.1 has GRCh38's three AMY1 copies and its length, and its bar still shows a long purple stretch, because the graph routes copies of a duplication through nodes off GRCh38's path.

Right-click a bar to open that haplotype's span in a linear view on its own assembly, with its CAT genes.

Matching the amylase haplotypes to Yilmaz et al.'s structures

Yilmaz et al. (2024) name each structure by its AMY1 count. Our five haplotypes land on H1a, H2A0, H3r, H5 and H7:

Span against GRCh38's	AMY1 copies	Structure
94 kb shorter	1	H1a
72 kb shorter	2, no AMY2A	H2A0
the same	3	H3r
94 kb longer	5	H5
188 kb longer	7	H7

HG00097.2's span is closest to H5's, yet it has four AMY1 copies. Count the AMY1 boxes on each bar to read copy number.

Telling the 1q21.1 inversion from an inverted paralog

A deletion or an extra copy changes how long a lane is. An inversion keeps the same sequence and reverses it. gfatools bubble flags a bubble, a place where haplotype paths split and rejoin, as inversion when its paths disagree about orientation. Press graph on the HPRC page's CFH / CFHR row, open Filter by... → Edit filters... on the bubbles track, and enter:

jexl:feature.inversion

Type chr1:144,260,000-144,610,000, the 1q21.1 locus, where the track flags one bubble. The flag fits a polymorphic inversion (present in some people) and an inverted paralog (a reversed duplicate copy) equally, and the graph draws the breakpoints as two deletion arcs because its edges record no orientation. Haplotype alignments tell the two apart. The figure below slices one haplotype with the inversion and one without out of HPRC's all-vs-GRCh38 PAF, each with its CAT annotation (Reproduce it end to end builds it), and the hg38 row between them agrees with the one without.

The 1q21.1 bubble the graph flags as an inversion, drawn as alignments. The pink ribbons are each haplotype's alignment to hg38, and a ribbon that crosses itself is an inversion. Between the two haplotype rows: the RefSeq genes, the bubbles track filtered to inversions, and the rGFA segments. The boxed pair on each row, PPIAL4F and PPIAL4E, is in opposite order on the two haplotypes.

Reproduce it end to end

The amylase build works in four steps, and the commands below run each one on any gbz-base database and any bgzipped, indexed assembly:

Ask the graph where each haplotype crosses two single-copy windows, one either side of the array. The distance between the two is that haplotype's span across the locus, and the spans fall into the published structures.
Pick one haplotype per structure and fetch only the locus from its assembly.
Count each haplotype's gene copies, and align it to its neighbour in the stack, so each band is an alignment between two haplotypes and can match copies GRCh38 lacks.
Shift each alignment from the fetched piece's coordinates back onto the whole contig, so every row draws in its own assembly's coordinates under its genes.

Get every haplotype's path through a window:

npx --yes -p @gmod/gbz-base gbz-base-query graph.gbz.db \
  --haplotype-index haplotypes.db \
  --sample GRCh38 --contig chr1 --interval 103540000..103541000 \
  --context 0 --alignments > window.json

Fetch one haplotype's copy of the locus:

samtools faidx HG00232_hap1.fa.gz 'HG00232#1#CM089991.1:103491008-103991760' > HG00232.1.fa

Count gene copies, keeping hits over 90% of the gene at 97% identity:

minimap2 -c -x asm20 -N 50 -p 0.5 HG00232.1.fa genes.fa |
  awk '($4-$3)/$2>=0.9 && $10/$11>=0.97 { c[$1]++ } END { for (g in c) print g, c[g] }'

Align two haplotypes to draw them as synteny:

# --secondary=no keeps the primary chain; the secondary ones are paralogous
# copies aligning to each other
minimap2 -c --eqx -x asm20 --secondary=no HG00232.1.fa NA18608.2.fa > adjacent.paf

The whole build, for four of the amylase haplotypes above and hg38:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_amylase_haplotypes.sh
bash build_amylase_haplotypes.sh

The build script writes the config.json that this session opens, with the five assemblies, their gene tracks and the amylase_adjacent alignments. The session below opens that config; point config= at your copy:

Config file

Goes at the top level of config.json, replacing any defaultSession there. See Default session.

{
  "defaultSession": {
    "name": "Amylase haplotypes from one AMY1 copy to seven, each aligned to the next",
    "views": [
      {
        "type": "LinearSyntenyView",
        "views": [
          {
            "assembly": "HG01361.1",
            "loc": "CM089019.1:103,831,655-104,050,048",
            "tracks": ["hprc_genes_HG01361_1"]
          },
          {
            "assembly": "hg38",
            "loc": "chr1:103,520,894-103,832,637",
            "tracks": ["hg38_ncbiRefSeq_ucsc"]
          },
          {
            "assembly": "HG00133.1",
            "loc": "CM090045.1:103,669,666-103,981,330",
            "tracks": ["hprc_genes_HG00133_1"]
          },
          {
            "assembly": "NA18608.2",
            "loc": "CM089849.1:103,796,766-104,203,421",
            "tracks": ["hprc_genes_NA18608_2"]
          },
          {
            "assembly": "HG00232.1",
            "loc": "CM089991.1:103,491,008-103,991,760",
            "tracks": ["hprc_genes_HG00232_1"]
          }
        ],
        "tracks": [
          ["amylase_adjacent"],
          ["amylase_adjacent"],
          ["amylase_adjacent"],
          ["amylase_adjacent"]
        ],
        "color": { "field": "strand" },
        "drawCurves": true,
        "levelHeights": [110, 110, 110, 110]
      }
    ]
  }
}

CLI (set-default-session)

jbrowse set-default-session --session - << 'EOF'
{
  "name": "Amylase haplotypes from one AMY1 copy to seven, each aligned to the next",
  "views": [
    {
      "type": "LinearSyntenyView",
      "views": [
        {
          "assembly": "HG01361.1",
          "loc": "CM089019.1:103,831,655-104,050,048",
          "tracks": ["hprc_genes_HG01361_1"]
        },
        {
          "assembly": "hg38",
          "loc": "chr1:103,520,894-103,832,637",
          "tracks": ["hg38_ncbiRefSeq_ucsc"]
        },
        {
          "assembly": "HG00133.1",
          "loc": "CM090045.1:103,669,666-103,981,330",
          "tracks": ["hprc_genes_HG00133_1"]
        },
        {
          "assembly": "NA18608.2",
          "loc": "CM089849.1:103,796,766-104,203,421",
          "tracks": ["hprc_genes_NA18608_2"]
        },
        {
          "assembly": "HG00232.1",
          "loc": "CM089991.1:103,491,008-103,991,760",
          "tracks": ["hprc_genes_HG00232_1"]
        }
      ],
      "tracks": [
        ["amylase_adjacent"],
        ["amylase_adjacent"],
        ["amylase_adjacent"],
        ["amylase_adjacent"]
      ],
      "color": { "field": "strand" },
      "drawCurves": true,
      "levelHeights": [110, 110, 110, 110]
    }
  ]
}
EOF

Open in JBrowse

Open this session in JBrowse Web ↗

HG01361.1, hg38, HG00133.1, NA18608.2 and HG00232.1, one AMY1 copy at the top to seven at the bottom, each aligned to the row under it by minimap2 and colored by strand. Each step up in copies opens a wedge over the genes only the longer row has.

A separate script builds the inversion figure from release 2's published all-vs-GRCh38 PAF. An assembler can deposit a contig in either orientation, so a reversed alignment alone is ambiguous; the test for an inversion is a block that reverses while the sequence either side of it stays forward. The script:

streams the PAF once, keeping every haplotype's alignments over a window that reaches well past the bubble on both sides
sorts each haplotype into inverted (bubble reversed, flanks forward), forward (forward throughout) or neither, and prints how many fall in each
keeps one inverted and one forward haplotype, with each one's contig length and CAT genes for the drawn window

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_hprc_inversion_synteny.sh
bash build_hprc_inversion_synteny.sh

Whole-genome synteny from a GFA's walks

A GFA that records each haplotype's walk converts to PAF against the reference walk, one record per stretch of shared nodes, for a whole-genome synteny view.

Convert the walks you want:

curl -fO https://raw.githubusercontent.com/cmdcolin/gfa-to-pairwise-paf/v1.0.0/gfa_to_pairwise_paf.py
# --contig-lengths: each assembly's .fai, since walks omit contig lengths
gzip -dc graph.gfa.gz | python3 gfa_to_pairwise_paf.py --reference GRCh38#0 \
  --queries HG01109#1,HG00099#1 --contig-lengths contigs.fai > graph.paf

Index the PAF for JBrowse:

jbrowse make-pif graph.paf --csi --out graph.pif.gz

The HPRC build downloads the 63 GB release graph, converts eight haplotypes' walks against GRCh38, indexes the PAF, and fetches each haplotype's CAT genes:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_hprc_multiway_synteny.sh
bash build_hprc_multiway_synteny.sh

Its output is hosted and opens at the CFH deletion, one lane per haplotype under GRCh38.

Hosting your own graph builds the gbz-base database for a graph of your own.

External links

Li H. gfatools, whose bubble subcommand flags the inversion.
gbz-base, which stores a GBZ as the SQLite database a window is range-requested out of.

Citations

HPRC release 2, the release whose graph, assemblies and CAT annotations this page reads.
Hughes AE et al. A common CFH haplotype, with deletion of CFHR1 and CFHR3, is associated with lower risk of age-related macular degeneration. Nature Genetics, 2006.
Yilmaz F, et al. Reconstruction of the human amylase locus reveals ancient duplications seeding modern-day variation. Science (2024). https://doi.org/10.1126/science.adn0609
Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics (2018). https://doi.org/10.1093/bioinformatics/bty191

Feedback on this tutorial is welcome: contact us.

Pangenome (HPRC) part 2: haplotypes against each other

Prerequisites#

Where the data comes from#

Configuring the haplotype lanes track#

CFH: haplotypes missing two genes#

C4: haplotypes with a copy GRCh38 lacks#

Amylase: counting AMY1 copies per haplotype#

Matching the amylase haplotypes to Yilmaz et al.'s structures#

Telling the 1q21.1 inversion from an inverted paralog#

Reproduce it end to end#

Whole-genome synteny from a GFA's walks#

See also#

External links#

Citations#