Phased trio analysis (1000 Genomes)
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.13. Desktop beta builds are coming soon.
hap-ibd identifies which stretches of a phased child's genome came down from the mother and which from the father. We paint those as one colored row per parental haplotype, so a meiotic crossover reads as a color change along the row.
Prerequisites
- a JBrowse to open them in: Desktop takes a local file by path, Web through Add track
- the
hg38assembly set up in JBrowse (assemblies guide) - Java 8+, for hap-ibd
python3node- htslib (
bgzip,tabix)
On Debian/Ubuntu, apt install tabix python3 default-jre covers most of it;
node comes from nodejs.org, and hap-ibd.jar is a
single download from its
releases page.
Where the data comes from
1000 Genomes Project phased low-coverage calls (1000 Genomes Project Consortium 2015), the Kinh-Vietnamese trio HG02024 (child), HG02026 (father) and HG02025 (mother), chr1 only.
The build script takes these files from their URLs, so there is nothing to download by hand.
- the phased trio VCF: hgdownload.soe.ucsc.edu/…/HG02024_VN049_KHVTrio.chr1.vcf.gzhttps://hgdownload.soe.ucsc.edu/gbdb/hg38/1000Genomes/trio/HG02024_VN049_KHV/HG02024_VN049_KHVTrio.chr1.vcf.gz
- the GRCh38 PLINK genetic map hap-ibd needs, the
no_chr_in_chrom_fieldvariant, since the trio VCF calls its chromosome1rather thanchr1: bochet.gcc.biostat.washington.edu/…/plink.GRCh38.map.ziphttps://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/plink.GRCh38.map.zip
Loading the hg38 assembly
The trio calls are on GRCh38, and the hap-ibd blocks below use its chr1 coordinates, so we load that assembly first.
Goes in the assemblies array of config.json. See Assemblies.
{
"name": "hg38",
"uri": "https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz",
"refNameAliases": {
"uri": "https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/hg38_aliases.txt"
},
"cytobands": "https://jbrowse.org/genomes/GRCh38/cytoBand.txt"
}jbrowse add-assembly https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz \
--name hg38 \
--refNameAliases https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/hg38_aliases.txt \
--config '{"cytobands":"https://jbrowse.org/genomes/GRCh38/cytoBand.txt"}'hg38 is one of the genomes JBrowse Desktop hosts, with gene tracks already set up: on the start screen click Show all available genomes and pick it. To load the files of this config instead:
In JBrowse Desktop, Open new genome on the start screen (or File → Open genome... in a session), then Open from a URL and paste, one per line:
https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz
https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz.fai
https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz.gziJBrowse reads the format off the file name. Then fill in:
- Genome name:
hg38 - refName aliases (under More options):
https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/hg38_aliases.txt - cytobands (under More options):
https://jbrowse.org/genomes/GRCh38/cytoBand.txt
Loading the trio's phased VCF
A phased VCF tags each variant with the haplotype it sits on (0|1 vs 1|0),
so you can follow each variant to the copy of the genome it came from.
The VCF loads on hg38 as an ordinary VariantTrack
(variant track guide). For your own trio,
swap uri for a bgzipped VCF with its .tbi beside it and the same chromosome
naming as the assembly. In JBrowse Web you can instead paste the URL into File
→ Open track..., which infers the adapter and the index.
Goes in the tracks array of config.json. See Tracks.
{
"type": "VariantTrack",
"trackId": "khv_trio_vcf",
"name": "KHV trio phased calls (chr1)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "VcfTabixAdapter",
"uri": "https://hgdownload.soe.ucsc.edu/gbdb/hg38/1000Genomes/trio/HG02024_VN049_KHV/HG02024_VN049_KHVTrio.chr1.vcf.gz"
}
}jbrowse add-track https://hgdownload.soe.ucsc.edu/gbdb/hg38/1000Genomes/trio/HG02024_VN049_KHV/HG02024_VN049_KHVTrio.chr1.vcf.gz \
--trackId khv_trio_vcf \
--name "KHV trio phased calls (chr1)" \
--assemblyNames hg38In 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://hgdownload.soe.ucsc.edu/gbdb/hg38/1000Genomes/trio/HG02024_VN049_KHV/HG02024_VN049_KHVTrio.chr1.vcf.gz
Click Next. JBrowse reads the adapter and track type off the file name. Then fill in:
- Track name:
KHV trio phased calls (chr1) - Assembly:
hg38
Click Add.
Showing the trio VCF as a genotype matrix
In the track menu, choose Display types → Multi-sample variant display (the multi-sample variant display), then check Show... → Show as genotype matrix. Each sample becomes a row and each variant a column, with black lines tying the columns back to their genomic positions.
Splitting each sample into two haplotype rows
Choose Rows → Per haplotype from the track menu:
- each sample splits into its two haplotypes, so the three trio members become six rows
- per-haplotype rows need phased genotypes, written
0|1; unphased calls (0/1) need a phasing program such as SHAPEIT first
Running hap-ibd to find segments shared with each parent
hap-ibd computes the matching stretches as "identical by descent" (IBD) segments. It takes a phased VCF and a genetic map in PLINK format.
The trio VCF calls its chromosome 1, with no chr prefix, so the run uses the
no_chr_in_chrom_field variant of the GRCh38 PLINK map:
# min-seed: shortest shared stretch (cM) hap-ibd starts a segment from
# min-output: shortest segment (cM) it writes
# both default to 2.0, so 1.0 also reports segments between 1 and 2 cM
java -jar hap-ibd.jar \
gt=HG02024_VN049_KHVTrio.chr1.vcf.gz \
map=plink.chr1.GRCh38.map \
out=trio min-seed=1.0 min-output=1.0The output is trio.ibd.gz, one row per shared segment, with columns sample1,
hap1, sample2, hap2, chrom, start, end, cM-length. In a trio every segment pairs
the child with one parent, and the child's two haplotypes split cleanly between
them:
| child haplotype | matches parent | inherited copy |
|---|---|---|
| HG02024:1 | HG02026 (father) | paternal |
| HG02024:2 | HG02025 (mother) | maternal |
The 1000 Genomes pedigree line VN049 HG02024 HG02026 HG02025 gives the roles:
father HG02026, mother HG02025.
Converting hap-ibd data into painted inheritance blocks
hap-ibd's output has gaps and short spurious segments from the statistical
phasing, so
hapibd_to_bed.py
merges them into clean blocks. Per child haplotype it:
- merges adjacent segments of the same parental copy into runs
- drops short interior runs, which are switch errors
- snaps each remaining crossover to the midpoint of the gap between runs so the blocks abut; a gap too wide to bridge stays blank
The script writes one BED9 line per block plus a parenthap label for the
painted track's four rows (father copy 1, father copy 2, mother copy 1, mother
copy 2), and its itemRgb colors the father's two copies blue and the mother's
red. It takes trio.ibd.gz and the child, father and mother sample IDs; bgzip
and tabix -p bed then index the output:
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/hapibd_to_bed.py
python3 hapibd_to_bed.py trio.ibd.gz HG02024 HG02026 HG02025 trio.hapibd.bed
jbrowse sort-bed trio.hapibd.bed | bgzip > trio.hapibd.bed.gz
tabix -p bed trio.hapibd.bed.gzsort-bed keeps the #-header line on top and
sorts the rest under LC_ALL=C, so the adapter reads the column names from the
header, needs no columnNames, and sees the same order in every locale.
Load trio.hapibd.bed.gz as a FeatureTrack with a
LinearMultiRowFeatureDisplay that draws one row per parenthap value, in
rows.domain order, painting each block with its BED itemRgb.
showLegend is
off because the row labels already name the four categories:
Goes in the tracks array of config.json. See Tracks.
{
"type": "FeatureTrack",
"trackId": "khv_trio_hapibd",
"name": "KHV trio hap-ibd haplotype blocks (chr1)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "BedTabixAdapter",
"disableGeneHeuristic": true,
"uri": "trio.hapibd.bed.gz"
},
"displays": [
{
"type": "LinearMultiRowFeatureDisplay",
"rows": {
"field": "parenthap",
"domain": ["Father hap1", "Father hap2", "Mother hap1", "Mother hap2"]
},
"showLegend": false
}
]
}jbrowse add-track-json '{
"type": "FeatureTrack",
"trackId": "khv_trio_hapibd",
"name": "KHV trio hap-ibd haplotype blocks (chr1)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "BedTabixAdapter",
"disableGeneHeuristic": true,
"uri": "trio.hapibd.bed.gz"
},
"displays": [
{
"type": "LinearMultiRowFeatureDisplay",
"rows": {
"field": "parenthap",
"domain": ["Father hap1", "Father hap2", "Mother hap1", "Mother hap2"]
},
"showLegend": false
}
]
}'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": "khv_trio_hapibd",
"name": "KHV trio hap-ibd haplotype blocks (chr1)",
"assemblyNames": ["hg38"],
"adapter": {
"type": "BedTabixAdapter",
"disableGeneHeuristic": true,
"uri": "trio.hapibd.bed.gz"
},
"displays": [
{
"type": "LinearMultiRowFeatureDisplay",
"rows": {
"field": "parenthap",
"domain": ["Father hap1", "Father hap2", "Mother hap1", "Mother hap2"]
},
"showLegend": false
}
]
}trio.hapibd.bed.gz is relative to a config.json. Replace it with its URL or its path on this computer.
Reading crossovers off the painted blocks
The four rows of the painted track are each parent's two copies, blue for father HG02026 and red for mother HG02025:
The blue rows together are the child's paternal chromosome. Where one of them is filled, it is the father's copy the child inherited there, so hap-ibd places a paternal crossover at each step between the blue rows. The red rows are the maternal chromosome in the same way.
The control in the figure is that no position has both blue rows filled, or both red rows, which would mean hap-ibd matched one child haplotype to both of a parent's copies. That holds along the whole chromosome. A position with neither row filled, such as the gap all four rows share around the 50M tick, is one where hap-ibd found no segment long enough to report. The blocks run straight through the centromere, because hap-ibd joins the markers on either side of it.
Comparing the painted blocks with the raw genotypes
To line the painted blocks up with the genotypes underneath:
- Drag the painting's track label above the VCF's.
- Uncheck Show... → Show as genotype matrix on the VCF track, so the phased multi-sample variant display draws each genotype at its genomic position.
- To show one parent's two rows only, as the figures below do, set the
painting's
rows.keptto those names, for example["Father hap1", "Father hap2"]. - To name the VCF's rows after the painting's rows, set the variant display's
rows.labelsto{ "HG02024 HP0": "Child hap1", "HG02024 HP1": "Child hap2", "HG02025 HP0": "Mother hap1", "HG02025 HP1": "Mother hap2", "HG02026 HP0": "Father hap1", "HG02026 HP1": "Father hap2" }, or rename the rows in Edit colors/arrangement... in its track menu.
Zoom to a few hundred kb around one boundary, where the block-step is obvious and the genotype columns resolve into individual variants. Start with the paternal crossover near chr1:29.7 Mb:
Near chr1:55.8 Mb the child's maternal haplotype steps between the mother's two copies:
The 1000 Genomes VCF is statistically phased, so its genotypes switch between the two parental copies more often than real crossovers do. hap-ibd's cM-length threshold filters most of those switches out of the painting, so the two crossovers above are well supported and the finer blocks are approximate. For crossover mapping, use a pedigree-aware method such as duoHMM.
Reproduce it end to end
build_khv_trio_hapibd.sh
runs the whole pipeline. It:
- downloads the trio VCF, hap-ibd and the genetic map
- runs hap-ibd and paints the BED with
hapibd_to_bed.py - downloads JBrowse
- writes
khv_trio_build/jbrowse2with aconfig.jsonholding the hg38 assembly plus the VCF and hap-ibd tracks
Serve the folder and open the URL serve prints:
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_khv_trio_hapibd.sh
bash build_khv_trio_hapibd.sh
npx --yes serve khv_trio_build/jbrowse2In JBrowse Desktop, File → Session → Open config.json or .jbrowse file...
opens khv_trio_build/jbrowse2/config.json directly.
See also
- Local ancestry (Dog10K)
- QTL mapping (BXD mice)
- Structural variants (1000 Genomes)
- LD at a selective sweep (human)
- Multi-row feature track
- Multi-sample variant display
- Config guide: Variant track
Citations
- 1000 Genomes Project Consortium (2015). A global reference for human genetic variation
- Zhou et al. (2020). A fast and simple method for detecting identity-by-descent segments in large-scale data, hap-ibd
- O'Connell et al. (2014). A general approach for haplotype phasing across the full spectrum of relatedness, duoHMM
Feedback on this tutorial is welcome: contact us.