LD across an inversion (mosquitoes)

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

The 2La chromosomal inversion of the malaria mosquito Anopheles gambiae spans about 22 Mb and suppresses crossing over, so linkage disequilibrium (LD, the correlation between variants) runs across it as one block. We compute the LD with plink2 --r2-phased, draw it with an LDTrack, and load the inversion as a structural variant genotyped per mosquito beneath it.

Prerequisites

a JBrowse to paste the tracks into (Web or Desktop)
PLINK 2.0 (plink2)
htslib (bgzip, tabix)
samtools
curl
python3
node, for the JBrowse CLI

Where the data comes from

Show full URLs

Ag1000G phase 2 AR1 (Anopheles gambiae 1000 Genomes Consortium 2020).

The build script takes these files from their URLs, so there is nothing to download by hand.

the phased haplotypes and their sample list for chromosome arm 2L, which the commands subset to one population at a time: ngs.sanger.ac.uk/…/shapeithttps://ngs.sanger.ac.uk/production/ag1000g/phase2/AR1/haplotypes/main/shapeit/
the sample metadata the population lists come from, CMgam (Cameroon) and GAgam (Gabon): ngs.sanger.ac.uk/…/samples.meta.txthttps://ngs.sanger.ac.uk/production/ag1000g/phase2/AR1/samples/samples.meta.txt
the AgamP4 reference and its gene models, which the gene track reads: ngs.sanger.ac.uk/…/genomehttps://ngs.sanger.ac.uk/production/ag1000g/phase3/genome/
the 2La tag SNPs, the ~200 positions whose allele marks which arrangement a chromosome has, which each mosquito's karyotype is scored from (Love et al. 2019): raw.githubusercontent.com/…/2La_targets.txthttps://raw.githubusercontent.com/rrlove/compkaryo/master/compkaryo/targets/2La_targets.txt

Loading the AgamP4 assembly and genes

The LD table and the inversion calls use 2L coordinates of the AgamP4 reference, so we load that assembly and its gene models first. The gene track reads the AgamP4.12 annotation.

Config file

Goes in the assemblies array of config.json. See Assemblies.

{
  "name": "AgamP4",
  "sequence": {
    "adapter": {
      "type": "BgzipFastaAdapter",
      "uri": "https://jbrowse.org/demos/ag1000g/AgamP4.fa.bgz"
    }
  }
}

CLI (add-assembly)

jbrowse add-assembly https://jbrowse.org/demos/ag1000g/AgamP4.fa.bgz \
  --name AgamP4 \
  --type bgzipFasta

Config file

Goes in the tracks array of config.json. See Tracks.

{
  "type": "FeatureTrack",
  "trackId": "agamp4_genes",
  "name": "AgamP4.12 genes",
  "assemblyNames": ["AgamP4"],
  "adapter": {
    "type": "Gff3TabixAdapter",
    "uri": "https://jbrowse.org/demos/ag1000g/AgamP4.sorted.gff3.gz"
  }
}

CLI (add-track)

jbrowse add-track https://jbrowse.org/demos/ag1000g/AgamP4.sorted.gff3.gz \
  --trackId agamp4_genes \
  --name "AgamP4.12 genes" \
  --assemblyNames AgamP4

In the app

In 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://jbrowse.org/demos/ag1000g/AgamP4.sorted.gff3.gz

Click Next. JBrowse reads the adapter and track type off the file name. Then fill in:

Track name: AgamP4.12 genes
Assembly: AgamP4

Click Add.

Precomputing 2L LD with PLINK

JBrowse draws LD from a precomputed table: PLINK correlates the variants and PlinkLDTabixAdapter reads its output. PLINK reads a binary fileset, so we first convert the phased VCF of common variants into one. --double-id sets each family id to the sample id:

plink2 --vcf common.vcf --double-id --allow-extra-chr --make-bed --out common

Then thin the variants, correlate them, and index the table. keep.CMgam.txt lists the Cameroon samples, two tab-separated columns of the same sample id, the family/individual pair plink asks for.

# the display uploads n(n-1)/2 cells, and ~800 SNPs across an arm already
# reach screen resolution, so keep roughly one variant per 50 kb
plink2 --bfile common --allow-extra-chr --keep keep.CMgam.txt --maf 0.2 \
  --chr 2L --write-snplist --out sel
awk -F'_' -v g=50000 '{p=$2+0; if (p >= nxt) {print $0; nxt = p + g}}' \
  sel.snplist > grid.snplist

# --r2-phased estimates r2 from haplotype frequencies, the statistic the
# display draws
# dprimeabs adds D' as a magnitude, the form the display reads
# --ld-window-r2 0 keeps the uncorrelated pairs
# PLINK 1.9 spells the pair `--r2 dprime`, with the columns at the same offsets
plink2 --bfile common --allow-extra-chr --keep keep.CMgam.txt \
  --extract grid.snplist \
  --r2-phased cols=chrom,pos,id,dprimeabs \
  --ld-window 999999 --ld-window-kb 1000000 --ld-window-r2 0 \
  --out ag1000g_2L_CMgam

# plink2 writes tabs and a commented header, which `tabix -H` returns
# sort-bed runs `sort -k1,1 -k2,2n` under LC_ALL=C and keeps the `#` line on
# top; this table sorts on the same first two columns
jbrowse sort-bed < ag1000g_2L_CMgam.vcor |
  bgzip > ag1000g_2L_CMgam.vcor.gz
tabix -s 1 -b 2 -e 2 -f ag1000g_2L_CMgam.vcor.gz

The track over that file is an LDTrack, and color.field picks which of the two metric columns, r² or D', the display reads:

Config file

Goes in the tracks array of config.json. See Tracks.

{
  "type": "LDTrack",
  "trackId": "ag1000g_2l_cmgam",
  "name": "Cameroon, both arrangements segregating (r²)",
  "assemblyNames": ["AgamP4"],
  "adapter": {
    "type": "PlinkLDTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/ag1000g_2L_CMgam.vcor.gz"
  },
  "displays": [
    {
      "type": "LDTrackDisplay",
      "color": { "field": "r2" },
      "variantLayout": "genomic",
      "showLegend": true,
      "height": 340
    }
  ]
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "LDTrack",
  "trackId": "ag1000g_2l_cmgam",
  "name": "Cameroon, both arrangements segregating (r²)",
  "assemblyNames": ["AgamP4"],
  "adapter": {
    "type": "PlinkLDTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/ag1000g_2L_CMgam.vcor.gz"
  },
  "displays": [
    {
      "type": "LDTrackDisplay",
      "color": { "field": "r2" },
      "variantLayout": "genomic",
      "showLegend": true,
      "height": 340
    }
  ]
}'

In the app

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": "ag1000g_2l_cmgam",
  "name": "Cameroon, both arrangements segregating (r²)",
  "assemblyNames": ["AgamP4"],
  "adapter": {
    "type": "PlinkLDTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/ag1000g_2L_CMgam.vcor.gz"
  },
  "displays": [
    {
      "type": "LDTrackDisplay",
      "color": { "field": "r2" },
      "variantLayout": "genomic",
      "showLegend": true,
      "height": 340
    }
  ]
}

The inversion genotyped per mosquito

The 2La inversion also loads as one <INV> record spanning the breakpoints, genotyped across every mosquito. The regular multi-sample variant display draws each genotype at the call's true span. END is the far breakpoint and each sample column holds one GT:

#CHROM  POS       ID   REF  ALT    QUAL  FILTER  INFO                     FORMAT  AN0007-C  AN0009-C
2L      20524058  2La  N    <INV>  .     PASS    SVTYPE=INV;END=42165532  GT      0/0       0/1

The samples TSV has a name column matching the VCF sample ids and a karyotype column naming each mosquito's pair of arrangements: 2L+a/2L+a, 2La/2L+a, 2La/2La, the + marking the non-inverted arrangement.

Load each population as a VariantTrack whose adapter includes the samples TSV, with a LinearMultiSampleVariantDisplay:

Config file

Goes in the tracks array of config.json. See Tracks.

{
  "type": "VariantTrack",
  "trackId": "ag1000g_2la_karyotype_cmgam",
  "name": "Cameroon, one row per mosquito",
  "assemblyNames": ["AgamP4"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/ag1000g_2La_CMgam.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/popgen/ag1000g_2La_CMgam_samples.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantDisplay",
      "facet": {
        "field": "karyotype",
        "domain": ["2L+a/2L+a", "2La/2L+a", "2La/2La"]
      },
      "rowColor": "karyotype",
      "referenceDrawingMode": "skip"
    }
  ]
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "VariantTrack",
  "trackId": "ag1000g_2la_karyotype_cmgam",
  "name": "Cameroon, one row per mosquito",
  "assemblyNames": ["AgamP4"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/ag1000g_2La_CMgam.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/popgen/ag1000g_2La_CMgam_samples.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantDisplay",
      "facet": {
        "field": "karyotype",
        "domain": ["2L+a/2L+a", "2La/2L+a", "2La/2La"]
      },
      "rowColor": "karyotype",
      "referenceDrawingMode": "skip"
    }
  ]
}'

In the app

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": "ag1000g_2la_karyotype_cmgam",
  "name": "Cameroon, one row per mosquito",
  "assemblyNames": ["AgamP4"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/popgen/ag1000g_2La_CMgam.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/popgen/ag1000g_2La_CMgam_samples.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantDisplay",
      "facet": {
        "field": "karyotype",
        "domain": ["2L+a/2L+a", "2La/2L+a", "2La/2La"]
      },
      "rowColor": "karyotype",
      "referenceDrawingMode": "skip"
    }
  ]
}

facet gives each karyotype class its own labelled band, and its domain orders the bands by dosage. referenceDrawingMode skip fills the track with the reference color and paints alt cells on top. The display draws a row for every sample in the file and divides the track height among them, so each population gets its own track. Gabon's two tracks are the same configs with CMgam replaced by GAgam in the trackIds and file names:

https://jbrowse.org/demos/popgen/ag1000g_2L_GAgam.vcor.gz
https://jbrowse.org/demos/popgen/ag1000g_2La_GAgam.vcf.gz
https://jbrowse.org/demos/popgen/ag1000g_2La_GAgam_samples.tsv

Scoring each mosquito's 2La karyotype from tag SNPs

The build script draws the <INV> call at the published 2La extent (Sharakhov et al. 2006, White et al. 2007) and scores each mosquito's karyotype as its mean alternate-allele count across the tag SNPs, rounded to a genotype, as MalariaGEN does for its Ag3 release. The score is trimodal, which the reproduce script checks.

Comparing the 2La LD block with karyotypes in Cameroon and Gabon

Stack the r² track of each population over the karyotype track of the same population, one row per mosquito.

Ag1000G chromosome arm 2L, the same window and settings throughout. Top: the published extents of 2La and of Vgsc, the two loci the blocks below sit on. r² fills the 2La extent in Cameroon, which segregates both arrangements, and is empty over that span in Gabon, which is near-fixed for the standard arrangement.

The block's edges line up with the published breakpoint coordinates.

A second block, at the low-coordinate end of the arm in both populations, is Vgsc, a sodium channel gene whose codon-995 substitutions confer pyrethroid insecticide resistance (Clarkson et al. 2021).
Gabon is near-fixed for the standard arrangement, so 2La is rare enough that the variants tagging it fall below the --maf cutoff.

Reproduce it end to end

build_ag1000g_ld.sh uses the published 2La span only as a test window, and prints the evidence for each choice it makes from there:

Each panel is one population, since correlation pooled across populations invents linkage none of them has. The script prints, per population, mean D' between variants more than 5 Mb apart inside the test window and outside it, because only a population with both arrangements can show the block.
It keeps common variants, thins them to a grid, and writes each panel's r² and D' table.
It bins D' to distant partners along the whole arm. The steps up and down are the inversion's breakpoints, recovered from the data alone.
It scores each mosquito's karyotype from the tag SNPs, prints the score histogram, which has to come out with three peaks, and writes the <INV> calls and a config.json opening on the inversion.

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_ag1000g_ld.sh
bash build_ag1000g_ld.sh              # writes ./ag1000g_ld_build/jbrowse2
npx --yes serve ag1000g_ld_build/jbrowse2

Citations

Anopheles gambiae 1000 Genomes Consortium (2020). Genome variation and population structure among 1142 mosquitoes of the African malaria vector species Anopheles gambiae and Anopheles coluzzii
Clarkson et al. (2021). The genetic architecture of target-site resistance to pyrethroid insecticides in the African malaria vectors Anopheles gambiae and Anopheles coluzzii
Love et al. (2019). In silico karyotyping of chromosomally polymorphic malaria mosquitoes in the Anopheles gambiae complex
Sharakhov et al. (2006). Breakpoint structure reveals the unique origin of an interspecific chromosomal inversion (2La) in the Anopheles gambiae complex
White et al. (2007). Molecular karyotyping of the 2La inversion in Anopheles gambiae

Feedback on this tutorial is welcome: contact us.

LD across an inversion (mosquitoes)

Prerequisites#

Where the data comes from#

Loading the AgamP4 assembly and genes#

Precomputing 2L LD with PLINK#

The inversion genotyped per mosquito#

Scoring each mosquito's 2La karyotype from tag SNPs#

Comparing the 2La LD block with karyotypes in Cameroon and Gabon#

Reproduce it end to end#

See also#

Citations#

Prerequisites

Where the data comes from

Loading the AgamP4 assembly and genes

Precomputing 2L LD with PLINK

The inversion genotyped per mosquito

Scoring each mosquito's 2La karyotype from tag SNPs

Comparing the 2La LD block with karyotypes in Cameroon and Gabon

Reproduce it end to end

See also

Citations