Pangenome (HPRC) part 3: repeat lengths across haplotypes

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.20. Desktop beta builds for Linux and macOS are on the releases page; the Windows build waits on code signing.

A tandem repeat is one sequence copied head to tail, and the number of copies varies from person to person. An assembled haplotype spans the whole array, so the Human Pangenome Reference Consortium's release 2 graph holds each haplotype's repeat at its full length. We read three arrays off it: the kringle IV type 2 copies in LPA and their two repeat types, a VNTR in an intron of ABCA7 that PacBio's TRGT also measured from reads, and the AMY1 gene copies in the amylase array. We:

at LPA, draw every haplotype's walk through the array as a bar with the walk rows layout, and open a tandem-repeat record's copies by type
at ABCA7, set TRGT's genotypes on the walk rows with the Repeat layout entry
choose the samples whose TRGT ticks turn red, and draw TRGT's alleles as a tandem-repeat record
at amylase, read AMY1 copies off walk rows, and stack one haplotype per structure in a linear synteny view

Experimental

The graph view is a beta plugin, and walk rows is its newest layout. The HPRC page launches it on the JBrowse 5 development build until JBrowse 5 ships. We welcome your feedback.

Prerequisites

the GraphGenomeView plugin. In your own JBrowse, install it from the plugin store; the HPRC page's launches and the hosted HPRC config the sessions below open already load it
for your own samples: TRGT and bcftools
for Which copy is which and The ABCA7 alleles across the samples: Node, to run our scripts
for Reproduce it end to end: DuckDB, and for its amylase build samtools built with libcurl, minimap2, python3 and Node.js with npm

# Debian/Ubuntu
sudo apt install samtools bcftools minimap2 nodejs
# or conda, which has every tool above
conda install -c bioconda -c conda-forge samtools bcftools minimap2 duckdb \
  nodejs trgt

The data

HPRC release 2 and PacBio's TRGT genotypes of a subset of its samples (Dolzhenko et al. 2024), most of which have a call at ABCA7.

You don't need these files to follow the page: the track configs below read the hosted results. The build script fetches them to rebuild those results. The tandem-repeat VCF step takes a cut of your own graph over your own array, and the TRGT step your own aligned HiFi reads and repeat catalogue.

What it is	File
TRGT's genotypes over the Genome in a Bottle repeat catalogue, a TRGTdb	adotto_hprc.tdb.tar
The catalogue itself	adotto_repeats.hg38.bed.gz
The assemblies, for the amylase build	assemblies_release2_v1.0.index.csv
The graph cut over the KIV-2 array, eight haplotypes' walks with their bases	hprc-v2.1-mc-grch38.kiv2.eight-haplotypes.gfa

Read by URL (no download needed)

What it is	File
Our tabix index files cut from the release's rGFA graph and from its haplotype walks, with the build recorded beside them	README.txt

Reading the LPA kringle repeat as a graph

LPA contains a tandem array of kringle IV type 2 (KIV-2) copies, tied to lipoprotein(a) levels, a heart-disease risk factor (Schmidt et al. 2016). We'll open it from the HPRC page:

Open the HPRC page, click the LPA example, and click Graph under More about this window. JBrowse opens on chr6:160,525,000-160,655,000.
Tick Show → Show bubble halos in the graph track's menu.

The strip along the top of the track draws each reference segment at its position on the ruler, in the color of its node below.

LPA over the force-directed graph track. The KIV-2 array is the haloed knot of loops on the right.

Each loop in the knot is a different number of copies. Click the array's purple label to lay the bubble's segments out alone, with a button back to the window.

Haplotypes' routes through the KIV-2 array

A walk is one haplotype's route through the graph, and the haplotypes track reads one per haplotype from tabix-indexed files. In the same window:

Type chr6:160,616,002-160,646,753, the array.
Display types → Feature display puts the rGFA graph track back to a row of segments.
Turn the haplotypes track on in the track selector. Display types → Graph in its menu draws a force-directed graph of the walks.
Check Haplotypes → The track's 8 assemblies in the same menu. Every haplotype's walks through the array have more nodes than a force-directed drawing handles, so this keeps the eight HPRC assemblies the hosted config loads, and walk rows below draw them all.

A node draws thicker the more walks visit it (as in Bandage), so the shared backbone is the thick line and copies on one haplotype are thin loops. Pick HG00133 under Haplotypes: its route stays dark while everything else fades, and a readout gives its length against the reference walk.

The KIV-2 cut with HG00133 picked under Haplotypes. The dark links are copies HG00133 walks and GRCh38 does not.

Every haplotype's KIV-2 array as a bar

Walk rows draw each walk as a bar, so every haplotype in the release fits in the track. In the haplotypes track's menu:

Walk rows under the Layout row and Uniform under the Color row turn each walk into a bar, longest first, blue on GRCh38's path through the graph and purple off it
Haplotypes → Every haplotype in the graph cuts the window again for the whole release

Read the bars this way:

Length counts copies. Each copy a haplotype adds lengthens its bar.
Color does not. The graph routes most copies through nodes off GRCh38's path, so most bars are largely purple.
Gold is steps the file did not load. The walk file abridges the costliest walks through the array, and a bar marks its length ≈ where the walk runs on past the window. Right-click one and take the Load the … abridged steps row, which counts them, to read it in full.
Hover a bar for its haplotype, length and excess over GRCh38, since the rows pack too thin to letter.

Every HPRC haplotype's walk across the KIV-2 array in walk rows. Purple is sequence off GRCh38's path; gold is an abridged walk's steps the file did not load.

Click the KIV-2 bubble boxed in the bubbles track. Its details give shortestAlleleLength and longestAlleleLength, the shortest and longest routes the rGFA holds there, and the bars run from about the one to about the other.

Two controls change which bars show:

Right-click a bar and pick Show only with its name to read one sample at full size.
Samples → Choose samples... under the Layout row picks names, and Samples → Every sample brings the cohort back.

Telling KIV-2A from KIV-2B

Walk rows give each haplotype's copy count. KIV-2's copies are near-identical, and the aligner's choice of which GRCh38 copy an extra one matches is arbitrary, so the walks leave open which copy is which. A record that lists each haplotype's copies settles it.

We host one as a track: a single VCF 4.5 <CNV:TR> record at the array over every haplotype whose walk reaches both flanks. Each allele lists its runs of one copy type and every copy's length, and a phased genotype puts each allele on its haplotype:

Config file

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

{
  "trackId": "hprc_kiv2_copies_all",
  "name": "LPA KIV-2 copies by unit, all HPRC haplotypes",
  "uri": "https://jbrowse.org/demos/hprc/hprc_kiv2_copies_all.vcf.gz",
  "assemblyNames": ["hg38"]
}

CLI (add-track)

jbrowse add-track https://jbrowse.org/demos/hprc/hprc_kiv2_copies_all.vcf.gz \
  --trackId hprc_kiv2_copies_all \
  --name "LPA KIV-2 copies by unit, all HPRC haplotypes" \
  --assemblyNames hg38

In the app

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track... and, in Add a track from file or URL, enter:

Main file: https://jbrowse.org/demos/hprc/hprc_kiv2_copies_all.vcf.gz

Click Next. JBrowse reads the adapter and track type off the file name. Then fill in:

Track name: LPA KIV-2 copies by unit, all HPRC haplotypes
Assembly: hg38

Click Add.

Open in JBrowse

Open this track in JBrowse Web ↗

Open in JBrowse Desktop ↗

Needs JBrowse Desktop 5.0 or newer. If nothing opens, copy the link and paste it into Desktop's Open JBrowse Web link..., on the start screen or under File → Session.

https://jbrowse.org/code/jb2/main/?config=https%3A%2F%2Fjbrowse.org%2Fucsc%2Fhg38%2Fconfig.json&session=spec-%7B%22sessionTracks%22%3A%5B%7B%22trackId%22%3A%22hprc_kiv2_copies_all%22%2C%22name%22%3A%22LPA+KIV-2+copies+by+unit%2C+all+HPRC+haplotypes%22%2C%22uri%22%3A%22https%3A%2F%2Fjbrowse.org%2Fdemos%2Fhprc%2Fhprc_kiv2_copies_all.vcf.gz%22%2C%22assemblyNames%22%3A%5B%22hg38%22%5D%7D%5D%2C%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr6%3A160%2C616%2C002-160%2C646%2C753%22%2C%22tracks%22%3A%5B%22hprc_kiv2_copies_all%22%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+LPA+KIV-2+copies+by+unit%2C+all+HPRC+haplotypes

Turn the track on and click the record. A Tandem repeat card opens above the details, each copy colored by its type. With this many haplotypes it starts on By allele: one bar per allele, longest first, and hovering a copy gives that allele's share of the haplotypes. By haplotype draws one bar per haplotype, and hovering a copy names its haplotype. Click a bar: the other bars fade and the Samples card narrows to the samples carrying that allele.

The KIV-2 record's Tandem repeat card, by haplotype. Each bar is one haplotype's array, each copy colored by its type.

The track selector also holds LPA KIV-2 copies by unit, eight HPRC haplotypes, few enough bars to label each by name.

To write a record like it for your own array, start from a cut of the graph over the array: a GFA whose W lines are the walks across it, from the reference nodes on either side, and whose S lines hold those nodes' bases. We host the cut for the eight haplotypes:¹

curl -fO https://jbrowse.org/demos/hprc/hprc-v2.1-mc-grch38.kiv2.eight-haplotypes.gfa

tandem-repeat-vcf.mjs splits each walk into copies wherever the reference array's first 24 bases recur, and groups copies within 1% of each other into a unit. It takes the cut and a BED row naming the array on the reference (chrom, start, end, name), and --unit-names names the units, the one with the most copies first:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/tandem-repeat-vcf.mjs

printf 'chr6\t160616002\t160646753\tKIV-2\n' > arrays.bed
node tandem-repeat-vcf.mjs hprc-v2.1-mc-grch38.kiv2.eight-haplotypes.gfa \
  --bed arrays.bed --name KIV-2 \
  --unit-names KIV-2A,KIV-2B > kiv2.vcf

Given one copy's exons as BED rows with their strand, --sites reports where the units differ inside each. At KIV-2 that is exon 1's positions 14, 41 and 86, the three synonymous sites LPA studies use to tell KIV-2B from KIV-2A:

printf 'chr6\t160617116\t160617276\texon 1\t0\t-\nchr6\t160618483\t160618665\texon 2\t0\t-\n' > kiv2_exons.bed
node tandem-repeat-vcf.mjs hprc-v2.1-mc-grch38.kiv2.eight-haplotypes.gfa \
  --bed arrays.bed --name KIV-2 \
  --sites kiv2_exons.bed > /dev/null

A repeat finder's output draws as these bars too once it is written as a <CNV:TR> record with each allele's runs and copy lengths.

Measuring the ABCA7 VNTR against TRGT's read-based genotypes

An intron of ABCA7 holds a VNTR (variable-number tandem repeat) tied to Alzheimer's disease risk (De Roeck et al. 2018). Each haplotype's walk length comes from its assembly, and PacBio's TRGT, a tandem repeat genotyper for HiFi reads, measured the same repeat in HPRC samples. We host its ABCA7 call as this track:

Config file

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

{
  "trackId": "hprc_abca7_trgt",
  "name": "TRGT repeat genotypes at ABCA7, 94 HPRC samples",
  "uri": "https://jbrowse.org/demos/hprc/hprc_abca7_trgt.vcf.gz",
  "csi": true,
  "assemblyNames": ["hg38"]
}

CLI (add-track)

jbrowse add-track https://jbrowse.org/demos/hprc/hprc_abca7_trgt.vcf.gz \
  --trackId hprc_abca7_trgt \
  --name "TRGT repeat genotypes at ABCA7, 94 HPRC samples" \
  --assemblyNames hg38 \
  --indexFile https://jbrowse.org/demos/hprc/hprc_abca7_trgt.vcf.gz.csi

In the app

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:

{
  "trackId": "hprc_abca7_trgt",
  "name": "TRGT repeat genotypes at ABCA7, 94 HPRC samples",
  "uri": "https://jbrowse.org/demos/hprc/hprc_abca7_trgt.vcf.gz",
  "csi": true,
  "assemblyNames": ["hg38"]
}

Open in JBrowse

Open this track in JBrowse Web ↗

Open in JBrowse Desktop ↗

Needs JBrowse Desktop 5.0 or newer. If nothing opens, copy the link and paste it into Desktop's Open JBrowse Web link..., on the start screen or under File → Session.

https://jbrowse.org/code/jb2/main/?config=https%3A%2F%2Fjbrowse.org%2Fucsc%2Fhg38%2Fconfig.json&session=spec-%7B%22sessionTracks%22%3A%5B%7B%22trackId%22%3A%22hprc_abca7_trgt%22%2C%22name%22%3A%22TRGT+repeat+genotypes+at+ABCA7%2C+94+HPRC+samples%22%2C%22uri%22%3A%22https%3A%2F%2Fjbrowse.org%2Fdemos%2Fhprc%2Fhprc_abca7_trgt.vcf.gz%22%2C%22csi%22%3Atrue%2C%22assemblyNames%22%3A%5B%22hg38%22%5D%7D%5D%2C%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr19%3A1%2C049%2C000-1%2C050%2C500%22%2C%22tracks%22%3A%5B%22hprc_abca7_trgt%22%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+TRGT+repeat+genotypes+at+ABCA7%2C+94+HPRC+samples

TRGT writes the same fields for your own samples, one VCF each, joined by trgt merge:

# one run per sample: HiFi reads aligned to GRCh38, and the repeat catalogue
trgt genotype --genome GRCh38.fa --reads sample.bam \
  --repeats adotto_repeats.hg38.bed --output-prefix sample
# genotype writes its VCF unsorted, and merge reads sorted, indexed input
bcftools sort -Oz -o sample.sorted.vcf.gz sample.vcf.gz
bcftools index -t sample.sorted.vcf.gz
trgt merge --vcf *.sorted.vcf.gz --genome GRCh38.fa --output-type z --output merged.vcf.gz

The session below lays the walks out as bars beside the genes, the catalogue's VNTR row and the TRGT genotypes:

Config file

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

{
  "defaultSession": {
    "name": "ABCA7 VNTR, walk rows",
    "sessionTracks": [
      {
        "type": "FeatureTrack",
        "trackId": "abca7_vntr",
        "name": "Tandem repeat catalogue (adotto)",
        "assemblyNames": ["hg38"],
        "adapter": {
          "type": "FromConfigAdapter",
          "features": [
            {
              "uniqueId": "abca7_vntr",
              "refName": "chr19",
              "start": 1049407,
              "end": 1050096,
              "name": "ABCA7 VNTR"
            }
          ]
        }
      }
    ],
    "views": [
      {
        "type": "LinearGenomeView",
        "assembly": "hg38",
        "loc": "chr19:1,049,000-1,050,500",
        "tracks": [
          {
            "trackId": "hg38_ncbiRefSeq_ucsc",
            "type": "LinearBasicDisplay",
            "filter": ["geneLike"],
            "displayMode": "compact",
            "height": 40
          },
          {
            "trackId": "abca7_vntr",
            "type": "LinearBasicDisplay",
            "height": 40
          },
          {
            "trackId": "hprc_abca7_trgt",
            "type": "LinearVariantDisplay",
            "height": 40
          },
          {
            "trackId": "hprc_v2_1_walk_lanes",
            "type": "LinearGraphDisplay",
            "layoutMode": "walkrows",
            "color": "uniform",
            "subgraphHaplotypes": [],
            "repeatTrackId": "hprc_abca7_trgt",
            "height": 600
          }
        ]
      }
    ]
  }
}

CLI (set-default-session)

jbrowse set-default-session --session - << 'EOF'
{
  "name": "ABCA7 VNTR, walk rows",
  "sessionTracks": [
    {
      "type": "FeatureTrack",
      "trackId": "abca7_vntr",
      "name": "Tandem repeat catalogue (adotto)",
      "assemblyNames": ["hg38"],
      "adapter": {
        "type": "FromConfigAdapter",
        "features": [
          {
            "uniqueId": "abca7_vntr",
            "refName": "chr19",
            "start": 1049407,
            "end": 1050096,
            "name": "ABCA7 VNTR"
          }
        ]
      }
    }
  ],
  "views": [
    {
      "type": "LinearGenomeView",
      "assembly": "hg38",
      "loc": "chr19:1,049,000-1,050,500",
      "tracks": [
        {
          "trackId": "hg38_ncbiRefSeq_ucsc",
          "type": "LinearBasicDisplay",
          "filter": ["geneLike"],
          "displayMode": "compact",
          "height": 40
        },
        {
          "trackId": "abca7_vntr",
          "type": "LinearBasicDisplay",
          "height": 40
        },
        {
          "trackId": "hprc_abca7_trgt",
          "type": "LinearVariantDisplay",
          "height": 40
        },
        {
          "trackId": "hprc_v2_1_walk_lanes",
          "type": "LinearGraphDisplay",
          "layoutMode": "walkrows",
          "color": "uniform",
          "subgraphHaplotypes": [],
          "repeatTrackId": "hprc_abca7_trgt",
          "height": 600
        }
      ]
    }
  ]
}
EOF

Open in JBrowse

Open this session in JBrowse Web ↗

Each row is one haplotype's walk between the flanking reference nodes. Hover a row for its haplotype and length.

The session names the TRGT track as its repeatTrackId, so the Layout row of the track menu gains a Repeat entry. Pick the ABCA7 record: each walk gets a black tick at the allele TRGT called for it (paired by length), red where the two lie more than 10% apart. At full size, as in the next section, the bars also divide into motif-length units.

The ABCA7 VNTR in walk rows with the TRGT record picked under Repeat. Ticks are TRGT's calls, red where a call is far from its walk.

Samples where reads and assemblies disagree

To compare samples, keep the ABCA7 record picked under Repeat. Samples → Choose samples..., under the Layout row of the graph track's menu, takes HG00099, HG03688, HG00741, HG02647, HG01943, HG02559 and HG04199 in that order. Their walks show in pairs in that order and fall into groups:

HG00099, HG03688 and HG00741 tick at each bar's end: reads and assemblies agree.
HG02647 and HG01943 turn red: TRGT calls each near-homozygous, but a walk sits far from both calls (both walks, in HG01943).
HG02559 and HG04199 each have a walk with no verdict. No read spans one of HG02559's alleles, so the tick on HG02559#1 is grey, and HG04199's assembly does not span the repeat: its second walk's readout shows that the contig ends before it.

Click the TRGT record to see why. Its sample table gives AL, the allele lengths behind each tick, and SD, the number of reads spanning each allele; one of HG02559's alleles has an SD of 0.

Samples' walks through the ABCA7 VNTR, each marked with TRGT's allele as a tick.

The ABCA7 alleles across the samples

TRGT writes each allele as its full sequence, so the details are a wall of CCCCGTGAGC. Our script rewrites the record as <CNV:TR> alleles, each a run of the locus's 51 bp motif:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/trgt-to-cnv-tr.mjs

# the hosted record has no MS field, TRGT's per-run spans, so each allele
# becomes one run of the motif, its copies the allele's length over 51 bp
node trgt-to-cnv-tr.mjs hprc_abca7_trgt.vcf.gz > hprc_abca7_cnvtr.vcf

Sort, bgzip and index the output, then add it as a track:

Config file

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

{
  "trackId": "hprc_abca7_cnvtr",
  "name": "TRGT alleles at ABCA7 as repeat records, 94 HPRC samples",
  "uri": "https://jbrowse.org/demos/hprc/hprc_abca7_cnvtr.vcf.gz",
  "assemblyNames": ["hg38"]
}

CLI (add-track)

jbrowse add-track https://jbrowse.org/demos/hprc/hprc_abca7_cnvtr.vcf.gz \
  --trackId hprc_abca7_cnvtr \
  --name "TRGT alleles at ABCA7 as repeat records, 94 HPRC samples" \
  --assemblyNames hg38

In the app

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track... and, in Add a track from file or URL, enter:

Main file: https://jbrowse.org/demos/hprc/hprc_abca7_cnvtr.vcf.gz

Click Next. JBrowse reads the adapter and track type off the file name. Then fill in:

Track name: TRGT alleles at ABCA7 as repeat records, 94 HPRC samples
Assembly: hg38

Click Add.

Open in JBrowse

Open this track in JBrowse Web ↗

Open in JBrowse Desktop ↗

Needs JBrowse Desktop 5.0 or newer. If nothing opens, copy the link and paste it into Desktop's Open JBrowse Web link..., on the start screen or under File → Session.

https://jbrowse.org/code/jb2/main/?config=https%3A%2F%2Fjbrowse.org%2Fucsc%2Fhg38%2Fconfig.json&session=spec-%7B%22sessionTracks%22%3A%5B%7B%22trackId%22%3A%22hprc_abca7_cnvtr%22%2C%22name%22%3A%22TRGT+alleles+at+ABCA7+as+repeat+records%2C+94+HPRC+samples%22%2C%22uri%22%3A%22https%3A%2F%2Fjbrowse.org%2Fdemos%2Fhprc%2Fhprc_abca7_cnvtr.vcf.gz%22%2C%22assemblyNames%22%3A%5B%22hg38%22%5D%7D%5D%2C%22views%22%3A%5B%7B%22type%22%3A%22LinearGenomeView%22%2C%22assembly%22%3A%22hg38%22%2C%22loc%22%3A%22chr19%3A1%2C049%2C000-1%2C050%2C500%22%2C%22tracks%22%3A%5B%22hprc_abca7_cnvtr%22%5D%7D%5D%7D&sessionName=JBrowse+docs%3A+TRGT+alleles+at+ABCA7+as+repeat+records%2C+94+HPRC+samples

Open chr19:1,049,000-1,050,500 and click the record. The Tandem repeat card starts on By allele here too. Hover a copy for its allele, that allele's share and its length.

The ABCA7 VNTR record's Tandem repeat card, by allele. The dotted line is GRCh38's allele.

Amylase: counting AMY1 copies per haplotype

A gene array varies the same way a tandem repeat does, in whole gene copies. We draw haplotypes' routes through the salivary amylase array as walk rows, one per structure:

On the HPRC page, click the AMY1 example, then Haplotypes under More about this window.
Display types → Graph in the track menu.
Enter HG01361.1, HG00133.2, HG00133.1, NA18608.2 and HG00232.1 in Graph settings... → Haplotypes.
Walk rows under the Layout row and Uniform under the Color row draw each haplotype's route as a bar.

Haplotypes' walks across the amylase array, each over its own CAT genes. Purple is sequence off GRCh38's path.

Count the AMY1 boxes on each bar to read copy number. As at KIV-2, purple does not count copies: HG00133.1 has GRCh38's three AMY1 copies and its length, and its bar still shows a long purple stretch. Each bar boxes the CAT genes of that haplotype's assembly, and GRCh38's row the view's GRCh38 genes. CAT names each copy after the GRCh38 copy it projected, so one name can label several boxes.

Yilmaz et al. (2024) name each structure by its AMY1 count:

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

Each haplotype aligned to the next

Walk rows count copies; a synteny stack shows which copy of one haplotype matches which copy of its neighbour. The commands below run on any gfa-to-tabix --walks file and any bgzipped, indexed assembly.

Fetch the script that reads a window's walks, then the walk reader it imports:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/walk_spans.mjs

# the majors walk_spans.mjs is written against; a new major can change the API
npm install @jbrowse/graph-walks@4 @gmod/tabix@3 generic-filehandle2@2

Where a haplotype's walk crosses a single-copy window on either side of the array gives its span across the locus, and the span sorts it into one of the structures above. Get where every walk crosses one window, one line per walk with its contig coordinates:

# the walk file is a path or a URL, with its .tbi beside it
node walk_spans.mjs graph.GRCh38.walks.bed.gz 'GRCh38#0#chr1' 103540000 103541000 > left.tsv

Fetch one haplotype's copy of the locus:

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

Cut the three genes out of GRCh38 as genes.fa, one record per gene:

# AMY2B, AMY2A and AMY1A as GRCh38 holds them, each record renamed to its gene
samtools faidx GRCh38.fa chr1:103554644-103579534 | sed '1s/.*/>AMY2B/' > genes.fa
samtools faidx GRCh38.fa chr1:103616651-103625780 | sed '1s/.*/>AMY2A/' >> genes.fa
samtools faidx GRCh38.fa chr1:103655519-103664554 | sed '1s/.*/>AMY1/' >> genes.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 session below opens the script's config.json; 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"],
        "encoding": { "color": "strand" },
        "mark": { "interpolate": "bump" },
        "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"],
      "encoding": { "color": "strand" },
      "mark": { "interpolate": "bump" },
      "levelHeights": [110, 110, 110, 110]
    }
  ]
}
EOF

Open in JBrowse

Open this session in JBrowse Web ↗

Haplotypes from the fewest AMY1 copies at the top to the most at the bottom, each aligned to the row below by minimap2.

Reproduce it end to end

PacBio publishes TRGT's calls for a subset of HPRC samples as a TRGTdb, most of which have a call at ABCA7, and the build turns that locus into the VCF this page loads, without re-genotyping anything. It:

downloads the TRGTdb and the repeat catalogue
reads the locus's alleles and every sample's calls out of the database
writes one VCF record: each allele's sequence as an ALT, and per sample the genotype, each allele's length (AL) and the reads spanning it (SD)

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

The TRGTdb is a directory of Parquet tables, and DuckDB reads a locus's calls straight out of it:

# every sample's allele numbers and spanning-read counts at one locus
duckdb -json -c "
  select regexp_extract(filename, 'sample\.(.*)\.pq', 1) as sample,
    allele_number, spanning_reads
  from read_parquet('hprc_100.tdb/sample.*.pq', filename = true)
  where LocusID = (select LocusID from read_parquet('hprc_100.tdb/locus.pq')
    where chrom = 'chr19' and start = 1049407)"

The amylase build runs the commands under Each haplotype aligned to the next for the haplotypes in that session and hg38. It:

reads where each haplotype's walk crosses a single-copy window on either side of the array, which gives its span and so its structure
picks one haplotype per structure and fetches only the locus from its assembly
counts each haplotype's gene copies, and aligns it to its neighbour in the stack
shifts each alignment back onto the whole contig's coordinates

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

External links

ggbandage, which draws KIV-2 and amylase copy number as walk rows in ggplot2.
ggtubemap, which draws a tandem repeat near LPA as a tube map in ggplot2.

Citations

Schmidt et al. (2016). Structure, function, and genetics of lipoprotein(a)
Wick et al. (2015). Bandage: interactive visualization of de novo genome assemblies
De Roeck et al. (2018). An intronic VNTR affects splicing of ABCA7 and increases risk of Alzheimer's disease
Dolzhenko et al. (2024). Characterization and visualization of tandem repeats at genome scale
TRGT repeat catalogues and HPRC genotypes
HPRC release 2
Yilmaz et al. (2024). Reconstruction of the human amylase locus reveals ancient duplications seeding modern-day variation
Li (2018). Minimap2: pairwise alignment for nucleotide sequences

Notes

gbz-base-query --keep cut it from HPRC's gbz-base database, and README.txt gives the command. A walk file holds node bases only when gfa-to-tabix --walks --sequences writes it. ↩

Feedback on this tutorial is welcome: contact us.

Pangenome (HPRC) part 3: repeat lengths across haplotypes

Prerequisites#

The data#

Reading the LPA kringle repeat as a graph#

Haplotypes' routes through the KIV-2 array#

Every haplotype's KIV-2 array as a bar#

Telling KIV-2A from KIV-2B#

Measuring the ABCA7 VNTR against TRGT's read-based genotypes#

Samples where reads and assemblies disagree#

The ABCA7 alleles across the samples#

Amylase: counting AMY1 copies per haplotype#

Each haplotype aligned to the next#

Reproduce it end to end#

See also#

External links#

Citations#

Notes

Prerequisites

The data

Reading the LPA kringle repeat as a graph

Haplotypes' routes through the KIV-2 array

Every haplotype's KIV-2 array as a bar

Telling KIV-2A from KIV-2B

Measuring the ABCA7 VNTR against TRGT's read-based genotypes

Samples where reads and assemblies disagree

The ABCA7 alleles across the samples

Amylase: counting AMY1 copies per haplotype

Each haplotype aligned to the next

Reproduce it end to end

See also

External links

Citations