LD at a selective sweep (human)

TL;DR: we look at linkage disequilibrium around the lactase gene, where selection for lactase persistence left one long block of correlated variants. JBrowse computes the r² triangle in the browser, straight from a phased VCF.

Prerequisites

bcftools built with libcurl, for the commands on this page and the reproduce script
htslib (tabix)
curl
python3
node, for the JBrowse CLI
bedGraphToBigWig, for the Fst lane
PLINK 2.0, for the Fst lane¹

Where the data comes from

1000 Genomes 30x high-coverage from NYGC (Byrska-Bishop et al. 2022), called natively on GRCh38, so no liftover sits between the calls and the hg38 coordinates the figures use.

phased chromosome 2, which the commands slice to a 3.4 Mb region: https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr2.filtered.SNV_INDEL_SV_phased_panel.vcf.gz
the release's own unrelated set, whose SAMPLE_NAME column is unrelated.samples. Relatives share long haplotypes for reasons that have nothing to do with a sweep: https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/1000G_2504_high_coverage.sequence.index
populations and superpopulations, narrowed to that unrelated set for panel.samples (EUR) and rest.samples (everything else): https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/20130606_g1k_3202_samples_ped_population.txt
the EUR slice the triangle is drawn from, rehosted so the figures and their live links load without the EBI round trip: https://jbrowse.org/demos/popgen/lct_1kg38_chr2_eur_wide.vcf.gz
the six-population slice the haplotype matrix reads: https://jbrowse.org/demos/popgen/lct_1kg38_chr2_6pop.vcf.gz

The gene, ClinVar and recombination lanes beside them are tracks of the hosted UCSC hg38 hub.

Reading the triangle

Start with what the picture means:

Red means two variants are almost always inherited together, white means they are independent, so the triangle shows where a stretch of chromosome travels as a unit.
It is a pairwise matrix turned on its corner, so the vertical axis is the distance between the two variants being compared.

To draw one, add an LDDisplay to an ordinary VariantTrack in an hg38 session:

Config file CLI (add-track-json)

{
  "type": "VariantTrack",
  "trackId": "kgp_lct_ld",
  "name": "LCT lactase-persistence LD, 1000G European panel (r²)",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/lct_1kg38_chr2_eur_wide.vcf.gz"
  },
  "displays": [
    {
      "type": "LDDisplay",
      "minorAlleleFrequencyFilter": 0.35,
      "useGenomicPositions": true,
      "showLegend": true,
      "forceLoad": true,
      "height": 360
    }
  ]
}

jbrowse add-track-json '{
  "type": "VariantTrack",
  "trackId": "kgp_lct_ld",
  "name": "LCT lactase-persistence LD, 1000G European panel (r²)",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/lct_1kg38_chr2_eur_wide.vcf.gz"
  },
  "displays": [
    {
      "type": "LDDisplay",
      "minorAlleleFrequencyFilter": 0.35,
      "useGenomicPositions": true,
      "showLegend": true,
      "forceLoad": true,
      "height": 360
    }
  ]
}'

What each setting does:

minorAlleleFrequencyFilter thins a dense callset to the common, block-tagging variants.
useGenomicPositions sizes each cell by genomic distance, so the block's edges land under the coordinates they are at.
forceLoad presses the FORCE LOAD button for you. r² is computed from the genotypes themselves, so a window this wide asks for more data than a track fetches unasked, and the lane arrives as a "too much data" banner instead. Set it when nobody is there to click: a figure, an embed, a notebook.
fetchSizeLimit sets that ceiling for the whole track, where forceLoad speaks for the one display that declares it.

Why there is a block here at all:

Selection favoured one variant, the allele that keeps lactase switched on into adulthood, and carried the whole run of neighbouring variants up in frequency with it, leaving them correlated (Bersaglieri et al. 2004).
That allele is rs4988235. Its dbSNP report carries the ClinVar entry and the frequency table population by population.

Cut the region out of the VCF

The triangle is drawn from what the file holds, so the slice you cut decides the picture:

Reach past both edges of the block, or those edges are not in frame.
Cut the same region twice, once over the whole release and once over the European panel the sweep happened in. r² is a correlation across every sample in the file, so the two files draw two different triangles.

# -r is a range request, so 3.4 Mb costs 3.4 Mb and not the 2.5 GB chromosome.
# -S is one sample name per line; -e drops the symbolic SV records, which are
# spans rather than the allele indicators the display correlates.
bcftools view -r chr2:133800000-137200000 -S unrelated.samples \
  -e 'ALT[0]~"<"' -Oz -o pooled.vcf.gz \
  https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.chr2.filtered.SNV_INDEL_SV_phased_panel.vcf.gz
tabix -p vcf pooled.vcf.gz

bcftools view -S panel.samples -Oz -o panel.vcf.gz pooled.vcf.gz
tabix -p vcf panel.vcf.gz

Two things to know about the window you just cut:

How wide is wide enough is a question for the file, not for the picture. The reproduce script bins r² against the causal variant by position and prints where the correlation falls away, which is where this window comes from.
The two files draw different variants. minorAlleleFrequencyFilter is a frequency in whatever samples the file holds, so a variant common in one panel and rare elsewhere clears the floor in the panel file and falls below it in the pooled one.

Compute Fst per variant

The Fst lane is plink2 --fst over panel.samples and rest.samples, written out as a bigWig for a quantitative track:

# plink2 takes the two panels as one categorical phenotype rather than as two
# sample lists, and wants FID beside IID: a #IID-only header is refused as "No
# entries correspond to loaded sample IDs" even when every ID matches
{ printf '#FID\tIID\tPOP\n'
  awk '{print $1"\t"$1"\tPANEL"}' panel.samples
  awk '{print $1"\t"$1"\tREST"}' rest.samples; } > fst_pops.txt

# method=wc is Weir and Cockerham; plink2 defaults to Hudson, which is a
# different number. report-variants is per variant rather than windowed, and
# --output-chr chrM keeps CHROM spelled chr2 rather than plink2's bare 2
plink2 --vcf pooled.vcf.gz --double-id --output-chr chrM --pheno fst_pops.txt \
  --fst POP method=wc report-variants vcols=chrom,pos,fst --out fst_site

# 1-based site to bedGraph interval, dropping the sites scored nan
awk 'NR>1 && $4!="nan" {printf "%s\t%d\t%d\t%.5f\n",$1,$2-1,$2,$4}' \
  fst_site.PANEL.REST.fst.var | sort -k1,1 -k2,2n > fst_site.bedgraph
printf 'chr2\t242193529\n' > hg38.chrom.sizes
bedGraphToBigWig fst_site.bedgraph hg38.chrom.sizes fst.bw

The block at two scales

Top, RefSeq genes and Weir and Cockerham Fst per variant across a wide span of chr2. Under the wedge, the same locus and allele-frequency floor twice, differing only in which samples went in, over that Fst lane at its own scale and the deCODE genetic map. Open in JBrowse: Wide scan ↗ (large dataset — this session takes a while to load) · The two triangles ↗ (large dataset — this session takes a while to load)

The lower frame is all block, so it cannot tell you the locus is unusual: everything in it sits on the swept haplotype. Read the lanes around it for the comparison:

Fst, top. Fst scores how differently two sets of samples carry a variant, so a variant one panel carries and the other mostly lacks scores high. Widened well past the block, the most differentiated sites in the span are the ones inside it. Read it per variant, which is what the build script scores: a sweep differentiates the variants on its own haplotype and leaves the rest of a bin on the background, so averaging a bin averages the signal away.
Genetic map. The block fills the span where the deCODE map (Halldorsson et al. 2019) reads flat, with a recombination hotspot at each end, where crossovers are frequent enough to break a haplotype up. The map counts crossovers in sequenced families, so it carries no LD of its own. If you reach for a different recombination track in that hub, check how it was made: the HapMap and 1000 Genomes maps there are estimated from LD, so they cannot check a triangle independently.
The two triangles. This haplotype swept in Europe. Pool that panel with populations the haplotype never reached and their chromosomes go into the same correlation carrying a different background, so every pair of variants looks less correlated than it is inside either group. That is the paler, patchier upper triangle.

The haplotypes behind the triangle

The same VCF draws those haplotypes in the same view, one lane below the triangle:

a LinearMultiSampleVariantMatrixDisplay in renderingMode: 'phased' gives one row per chromosome and one column per variant
colorBy puts population in the sidebar stripe, so which populations carry the block reads off the rows

Config file CLI (add-track-json)

{
  "type": "VariantTrack",
  "trackId": "kgp_lct_haplotypes",
  "name": "1000 Genomes haplotypes across LCT (one row per haplotype)",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/lct_1kg38_chr2_6pop.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/genomes/hg19/1000g.sorted.csv.gz"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantMatrixDisplay",
      "renderingMode": "phased",
      "colorBy": "population",
      "minorAlleleFrequencyFilter": 0.35,
      "forceLoad": true,
      "height": 700
    }
  ]
}

jbrowse add-track-json '{
  "type": "VariantTrack",
  "trackId": "kgp_lct_haplotypes",
  "name": "1000 Genomes haplotypes across LCT (one row per haplotype)",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/lct_1kg38_chr2_6pop.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/genomes/hg19/1000g.sorted.csv.gz"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantMatrixDisplay",
      "renderingMode": "phased",
      "colorBy": "population",
      "minorAlleleFrequencyFilter": 0.35,
      "forceLoad": true,
      "height": 700
    }
  ]
}'

Run the clustering two ways:

from the track menu, Clustering → Cluster rows by genotype...
baked into a session with the runClustering and clusterRegion model properties, which is what the figure below does

An LD triangle over the haplotypes it summarises: 1000 Genomes chromosomes at LCT/MCM6, one row each, clustered by genotype. The pale slab is one cluster of near-identical chromosomes, uniform across the block that fills the triangle above. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The highlight across both lanes is the 89 kb of LCT and MCM6 that selection acted on.

Ordering is what makes a block visible. In file order the same matrix is a plaid at any size, because a block is a set of alleles travelling together and which of them is the non-reference allele varies from site to site. Clustering puts near-identical chromosomes next to each other, and a swept haplotype carries little variation of its own, so it resolves into one slab.
The ClinVar lane marks the causal variant, not the matrix. rs4988235 falls below the figure's own frequency floor, so it is not one of the columns, and that lane places it independently of the rows the clustering ran on.
Narrow that lane or it marks nothing. It is the hub's ClinVar track filtered with jexl:get(feature,'phenotypeList')=='LACTASE PERSISTENCE'; unfiltered it draws every ClinVar record in the window, congenital lactase deficiency in the same gene included.

The subsample behind the figure

Over the whole release each haplotype row falls well below a pixel in a lane this tall and averages into its neighbours, leaving a flat wash whatever the ordering. This figure reads a subsample of six populations instead, built by the third script under Reproduce it end to end.

Reproduce it end to end

build_lct_ld.sh does the whole build for you:

cuts the region out of the 1000 Genomes 30x callset without downloading it
writes a ready-to-serve config carrying both LD lanes, the Fst lane, the genetic map and the haplotype matrix
prints the two PLINK tables the window and the panel choice rest on: r² against rs4988235 binned by position, and mean pairwise r² inside the block for the panel against the pooled release

The VCFs it writes hold genotypes only, and JBrowse computes the r² from them as it draws the triangle.

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_lct_ld.sh
bash build_lct_ld.sh                  # builds ./lct_ld_build/jbrowse2
npx --yes serve lct_ld_build/jbrowse2 # then open the printed URL

The wide Fst lane is a second file, from build_lct_fst_scan.sh, because the slice above it stops a little either side of the block:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_lct_fst_scan.sh
bash build_lct_fst_scan.sh            # builds ./lct_fst_scan_build

The million-site slice it reads takes a few minutes to come down.

The subsampled haplotype matrix is a third file, from build_lct_haploblock.sh, since the whole release draws a flat wash at that lane height:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_lct_haploblock.sh
bash build_lct_haploblock.sh          # builds ./lct_haploblock_build

When the cohort is too large to correlate live

An LDTrack reads r² PLINK has already computed, so the browser fetches a table of pairs rather than the genotypes behind them. The adapter takes plink2's .vcor or PLINK 1.9's .ld, and LD across an inversion (mosquitoes) goes that way over a 22 Mb inversion.

References

1000 Genomes Project Consortium (2015). A global reference for human genetic variation
Bersaglieri et al. (2004). Genetic signatures of strong recent positive selection at the lactase gene
Byrska-Bishop et al. (2022). High-coverage whole-genome sequencing of the expanded 1000 Genomes Project cohort including 602 trios
Halldorsson et al. (2019). Characterizing mutagenic effects of recombination through a sequence-level genetic map

Notes

PLINK 1.9 does the same work under different spellings, and writes .ld where plink2 writes .vcor. JBrowse reads either. ↩

Feedback on this tutorial is welcome: contact us.

LD at a selective sweep (human)

Prerequisites#

Where the data comes from#

Reading the triangle#

Cut the region out of the VCF#

Compute Fst per variant#

The block at two scales#

The haplotypes behind the triangle#

The subsample behind the figure#

Reproduce it end to end#

When the cohort is too large to correlate live#

See also#

References#