LD at a selective sweep (human)

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.8. Desktop beta builds are coming soon.

We look at linkage disequilibrium around the lactase gene, where selection for lactase persistence left one long block of correlated variants. PLINK correlates the phased genotypes and JBrowse draws the triangle from its output.

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 1.9 for the r² tables and PLINK 2.0 for the Fst lane and the frequency filter¹

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 two r² tables the triangles are drawn from, one per cohort, as PLINK wrote them: https://jbrowse.org/demos/popgen/lct_1kg38_chr2_eur.ld.gz and https://jbrowse.org/demos/popgen/lct_1kg38_chr2_pooled.ld.gz
the EUR slice they were computed from, rehosted so the 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, point an LDTrack at the r² table PLINK wrote below, in an hg38 session:

Config file

{
  "type": "LDTrack",
  "trackId": "kgp_lct_ld",
  "name": "LCT lactase-persistence LD, 1000G European panel (r²)",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "PlinkLDTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/lct_1kg38_chr2_eur.ld.gz"
  },
  "displays": [
    {
      "type": "LDTrackDisplay",
      "useGenomicPositions": true,
      "showLegend": true,
      "height": 360
    }
  ]
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "LDTrack",
  "trackId": "kgp_lct_ld",
  "name": "LCT lactase-persistence LD, 1000G European panel (r²)",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "PlinkLDTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/lct_1kg38_chr2_eur.ld.gz"
  },
  "displays": [
    {
      "type": "LDTrackDisplay",
      "useGenomicPositions": true,
      "showLegend": true,
      "height": 360
    }
  ]
}'

JBrowse Desktop

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": "LDTrack",
  "trackId": "kgp_lct_ld",
  "name": "LCT lactase-persistence LD, 1000G European panel (r²)",
  "assemblyNames": ["hg38"],
  "adapter": {
    "type": "PlinkLDTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/lct_1kg38_chr2_eur.ld.gz"
  },
  "displays": [
    {
      "type": "LDTrackDisplay",
      "useGenomicPositions": true,
      "showLegend": true,
      "height": 360
    }
  ]
}

What each setting does:

useGenomicPositions sizes each cell by genomic distance, so the block's edges land under their coordinates
ldMetric picks which of the file's columns to draw. This table has both r² and D', so either reads; a file without a DP column disables the D' row rather than drawing zeros

The allele-frequency floor is not a display setting here. It is applied when the variants are picked for correlation, so it is a property of the file. That is also why the two cohorts below are a fair comparison.

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 give 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.

Correlate the variants with PLINK

Pick the common variants per cohort, then correlate every pair of them. The MAF floor keeps the table to a size a browser can draw. Every pair is a row, so n variants cost n(n-1)/2 of them.

# 0.35 is high for a MAF floor and deliberately so: it keeps the variants that
# tag the block rather than every rare one riding on it. Frequency is measured
# in THIS file's samples, so each cohort keeps its own set.
plink2 --vcf panel.snvs.vcf.gz --double-id --allow-extra-chr --output-chr chrM \
  --set-missing-var-ids @:# --maf 0.35 --chr chr2 --write-snplist --out sel

# dprime adds D' beside r², which is the display's other metric.
# --ld-window-r2 0 keeps the uncorrelated pairs, so white cells are drawn as
# white rather than left absent, and the two window flags have to be raised
# together. The defaults cut off after 10 variants or 1 Mb, whichever comes
# first, which would clip this block at both.
plink --vcf panel.snvs.vcf.gz --double-id --allow-extra-chr --output-chr chrM \
  --set-missing-var-ids @:# --extract sel.snplist \
  --r2 dprime --ld-window 999999 --ld-window-kb 4000 --ld-window-r2 0 \
  --out lct_1kg38_chr2_eur

# tabix needs real tabs and a commented header. plink pads its columns with
# spaces to align them, which is not the same thing, so squeeze the runs to
# tabs and mark the header before indexing.
awk 'NR==1{$1=$1; print "#" $0; next} {$1=$1; print}' OFS='\t' \
  lct_1kg38_chr2_eur.ld | bgzip > lct_1kg38_chr2_eur.ld.gz
tabix -s 1 -b 2 -e 2 -f lct_1kg38_chr2_eur.ld.gz

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

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

CLI (add-track-json)

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

JBrowse Desktop

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, as 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 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

A bigger span

LD across an inversion (mosquitoes) draws the same track type over a 22 Mb inversion, where the variants have to be thinned to a grid before they are correlated rather than only filtered by frequency.

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

The two are separate programs, not versions to choose between, and this page uses each where it is the simpler one. plink2 gained --r2-phased, which writes the same table as a .vcor under column names of its own, in the a6 alphas; on an earlier plink2 the flag is simply absent, which is why the r² step here is PLINK 1.9's. JBrowse's PlinkLDTabixAdapter reads either spelling: it resolves the columns from the header rather than by position. ↩

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#

Correlate the variants with PLINK#

Compute Fst per variant#

The block at two scales#

The haplotypes behind the triangle#

The subsample behind the figure#

Reproduce it end to end#

A bigger span#

See also#

References#