Phased trio analysis (1000 Genomes)
TL;DR: hap-ibd tells 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
- nothing to read along. Everything below is for building the tracks yourself
- 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 phased trio VCF: https://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: https://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/plink.GRCh38.map.zip - the hg38 reference sequence the reproduce script's own JBrowse instance opens on, rehosted: https://jbrowse.org/genomes/GRCh38/fasta/GRCh38.fa.gz
The trio VCF
A trio is a mother, father, and child sequenced together. A phased VCF tags each
variant with the haplotype it sits on (0|1 vs 1|0), so each variant can be
followed to the copy of the genome it came from.
This page uses the phased VCF above, the Kinh-Vietnamese trio HG02024, chr1 only.
Everything here is on hg38. Add the VCF with jbrowse add-track or the in-app
"Add track" workflow, both covered in the
variant track guide.
Enabling the matrix view
Switch the track to the Multi-sample variant display (matrix). Each sample becomes a row and each variant a column, with black lines tying the columns back to their genomic positions.
Enabling the phased mode
Turn on Rendering mode → Phased from the track menu:
- it splits each sample into its two haplotypes, so the three trio members become six rows
- it needs genotypes written with the
0|1separator rather than0/1; getting there from unphased calls takes a phasing program like SHAPEIT
That last move is wider than the default display will draw: it stops above its feature-density limit, where the matrix keeps going because a column is a variant rather than a position.
Reading matching haplotypes off the matrix
Every row is now a strip of colored blocks, and matching stretches between rows jump out: the child's two haplotypes match the mother's in some blocks and the father's in others. The rest of this tutorial turns that by-eye pattern into a painted track.
Finding the matching blocks programmatically
hap-ibd computes that matching, as "identical by descent" blocks. hap-ibd is built for population-scale cohorts and runs on a single trio VCF, and the run needs two things:
- a phased VCF, like the trio dataset above
- a genetic map in PLINK format (hap-ibd's README links GRCh38 ones)
Grab hap-ibd.jar from the
releases page along with
those maps.
Running hap-ibd
The trio VCF calls its chromosome 1, with no chr prefix, so reach for the
no_chr_in_chrom_field variant of the GRCh38 PLINK map:
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 roles come from the 1000 Genomes pedigree line
VN049 HG02024 HG02026 HG02025: father HG02026, mother HG02025.) Within one
child haplotype, the matching parental copy flips between the parent's copy 1
and copy 2 at each crossover. Those flips are what the track below paints.
hap-ibd's output has gaps, plus short spurious segments from the statistical phasing, so it is collapsed into clean blocks before painting.
Converting hap-ibd data into painted inheritance blocks
The goal is one row per parental haplotype (father copy 1, father copy 2, mother copy 1, mother copy 2), with the child's inherited chromosome tiled across each parent's pair of rows. A crossover then shows up as a block stepping from one row to its partner.
hapibd_to_bed.py
does the cleanup. Per child haplotype it:
- merges adjacent segments of the same parental copy into runs,
- drops short interior runs (the switch-error specks), and
- snaps each remaining crossover to the midpoint of the gap between runs so the blocks abut (real gaps, like the centromere, stay blank).
Out comes one BED9 line per block plus a parenthap label, with the father's
two copies in blues and the mother's in reds via itemRgb. Feed it
trio.ibd.gz plus the child, father, and mother sample IDs, then bgzip and
tabix -p bed so the BedTabixAdapter can read it:
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 can read the column names off
the file and the order does not shift with your locale.
Load the result as a FeatureTrack with a LinearMultiRowFeatureDisplay:
partitionFielddraws one row per distinct value it finds, soparenthapgives the four parental-haplotype rowsrowOrdersets their top-to-bottom order- a BED carrying
itemRgbis painted with it automatically, no extra color config needed
{
"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",
"partitionField": "parenthap",
"showLegend": false,
"rowOrder": ["Father hap1", "Father hap2", "Mother hap1", "Mother hap2"]
}
]
}Two things about the config above:
- The BED's
#-header line names its columns, so the adapter needs nocolumnNames, andparenthapis the one the display partitions on. showLegendis off, because the color and the row label already carry the same four categories.
Reading the painted crossovers
The four rows are each parent's two copies, blues for father HG02026 and reds for mother HG02025:
Read the rows in pairs:
- Blue rows are the child's paternal chromosome. Exactly one of them is filled at any position, and that is which of the father's two copies the child got there there; every step between the blue rows is a crossover.
- Red rows work the same way for the maternal chromosome.
That rule is the figure's own control: two filled blue rows at a position, or neither, means hap-ibd matched one child haplotype to both of the father's copies or to neither. The centromere is the blank with no markers to match on.
Relating the painting back to the genotypes
Stack the painting directly above the same VCF in the phased multi-sample variant display, which draws genotypes at their real genomic positions. Matrix mode spaces its columns evenly, on a scale of its own.
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:
The maternal chromosome does the same thing at its own boundaries. Near chr1:55.8 Mb the child's maternal haplotype steps between the mother's two copies:
The genotypes underneath switch between the two parental copies more often than real crossovers do, and the painting above summarises those switches away.
Where the boundaries come from
This 1000 Genomes VCF is statistically phased, and its haplotypes carry switch errors, which are the extra copy-switches visible in the genotype rows. hap-ibd's cM-length threshold filters most of them out, so its blocks track the real boundaries more closely and the two crossovers above are the well-supported ones; the finer blocks are approximate. hap-ibd gives paintable inheritance blocks, and crossover mapping proper uses a pedigree-aware method such as duoHMM.
Reproduce it end to end
build_khv_trio_hapibd.sh
runs the whole pipeline in one shot. It downloads the trio VCF, hap-ibd, and the
genetic map, runs hap-ibd, paints the BED with
hapibd_to_bed.py,
downloads JBrowse, and writes a config.json with the hg38 assembly plus the
VCF and hap-ibd tracks.
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_khv_trio_hapibd.sh
bash build_khv_trio_hapibd.sh # builds ./khv_trio_build/jbrowse2
npx --yes serve khv_trio_build/jbrowse2 # then open the printed URLIt needs java, python3, node, and htslib (bgzip and tabix). Opening
khv_trio_build/jbrowse2/config.json in JBrowse Desktop via File -> Session
-> Open config.json or .jbrowse file... gives the same view without serving
anything.
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
- Variant track
References
- 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.