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

a JBrowse to paste the tracks into (Web or Desktop); every file here is a URL, so Desktop needs nothing hosted
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

Red means two variants are almost always inherited together, white means they are independent. The triangle is a pairwise matrix turned on its corner, so the vertical axis is the distance between the two variants compared.

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

Config file CLI (add-track-json) JBrowse Desktop

{
  "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
    }
  ]
}'

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": "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 their coordinates
forceLoad presses the FORCE LOAD button for you, since r² is computed from the genotypes and a window this wide exceeds what a track fetches unasked. Set it where nobody is there to click: a figure, an embed, a notebook
fetchSizeLimit raises that ceiling for the whole track instead

The block is a selective sweep: the allele that keeps lactase switched on into adulthood, rs4988235, rose in frequency and carried its neighbouring variants with it (Bersaglieri et al. 2004). Its dbSNP report carries the ClinVar entry and the per-population frequency table.

Cut the region out of the VCF

The slice decides the picture: reach past both edges of the block, and cut the 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

The reproduce script bins r² against the causal variant by position and prints where the correlation falls away, which is where this window's width comes from. The two files also draw different variants, since minorAlleleFrequencyFilter is a frequency in whatever samples the file holds.

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 the lanes around it carry the comparison:

Fst, top. Fst scores how differently two sets of samples carry a variant. Widened well past the block, the most differentiated sites in the span are the ones inside it. It is scored per variant, since a sweep differentiates the variants on its own haplotype and leaves the rest of a bin on the background
Genetic map. The block fills the span where the deCODE map (Halldorsson et al. 2019) reads flat, with a recombination hotspot at each end. The map counts crossovers in sequenced families, so it carries no LD of its own; the HapMap and 1000 Genomes maps in the same hub are estimated from LD and cannot check a triangle independently
The two triangles. The haplotype swept in Europe. Pooling that panel with populations it never reached makes every pair of variants look less correlated than it is inside either group: the paler, patchier upper triangle

The haplotypes behind the triangle

The same VCF draws the haplotypes one lane below the triangle: a LinearMultiSampleVariantMatrixDisplay in renderingMode: 'phased' gives one row per chromosome and one column per variant, and colorBy puts population in the sidebar stripe.

Config file CLI (add-track-json) JBrowse Desktop

{
  "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
    }
  ]
}'

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": "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 is LCT and MCM6. rs4988235 is an enhancer variant in MCM6 intron 13.

Ordering is what makes a block visible. In file order the matrix is a plaid, because which allele is non-reference varies from site to site. Clustering puts near-identical chromosomes together, and a swept haplotype carries little variation of its own, so it resolves into one slab
The ClinVar lane marks the causal variant. rs4988235 falls below the frequency floor, so it is not a column, and the lane places it independently
Narrow that lane or it marks nothing. It is the hub's ClinVar track filtered with jexl:feature.phenotypeList=='LACTASE PERSISTENCE'; unfiltered it draws every ClinVar record in the window

The subsample behind the figure

Over the whole release each haplotype row falls well below a pixel and averages into 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 cuts the region out of the callset without downloading it, writes a ready-to-serve config with every lane above, and prints the two PLINK tables the window and panel choice rest on: r² against rs4988235 binned by position, and mean pairwise r² inside the block for the panel against the pooled release.

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:

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 subsampled haplotype matrix is a third file, from build_lct_haploblock.sh:

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, from plink2's .vcor or PLINK 1.9's .ld. 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#