# LD across an inversion (mosquitoes) **TL;DR:** a 22 Mb inversion reads as one block, from `plink2 --r2-phased` output through an [`LDTrack`](https://jbrowse.org/jb2-staging/docs/config/ldtrack). The same inversion also loads as a structural variant genotyped per mosquito. ## Prerequisites - [PLINK 2.0](https://www.cog-genomics.org/plink/2.0/) (`plink2`), labeled alpha for years despite being the version in general use - htslib (`bgzip`, `tabix`) - `samtools` - `curl` - `python3` - `node`, for the [JBrowse CLI](https://jbrowse.org/jb2-staging/docs/cli) ## Where the data comes from Ag1000G phase 2 AR1 ([Anopheles gambiae 1000 Genomes Consortium 2020](https://doi.org/10.1101/gr.262790.120)), whose terms of use were lifted in March 2022, so nothing here needs registration or a data-access agreement. - the phased haplotypes and their sample list for chromosome arm 2L, which the commands subset to one population at a time: https://ngs.sanger.ac.uk/production/ag1000g/phase2/AR1/haplotypes/main/shapeit/ - the sample metadata the population lists come from, `CMgam` (Cameroon) and `GAgam` (Gabon): https://ngs.sanger.ac.uk/production/ag1000g/phase2/AR1/samples/samples.meta.txt - the AgamP4 reference and its gene models, which the gene lane reads: https://ngs.sanger.ac.uk/production/ag1000g/phase3/genome/ - the 2La tag SNPs, the ~200 positions whose allele says which arrangement a chromosome carries, which each mosquito's karyotype is scored from ([Love et al. 2019](https://doi.org/10.1534/g3.119.400445)): https://raw.githubusercontent.com/rrlove/compkaryo/master/compkaryo/targets/2La_targets.txt - the finished `CMgam` LD table, rehosted so the track blocks on this page load without the build: https://jbrowse.org/demos/popgen/ag1000g_2L_CMgam.vcor.gz - the 2La genotypes per mosquito: https://jbrowse.org/demos/popgen/ag1000g_2La_CMgam.vcf.gz - the karyotype table the sample lane is grouped by: https://jbrowse.org/demos/popgen/ag1000g_2La_CMgam_samples.tsv ## The 2La inversion as one LD block Two facts set up the page: - **Inverted and standard arrangements cannot recombine in a heterozygote**, so wherever both are present the whole segment stays correlated. - **The 2La inversion in _Anopheles gambiae_ spans roughly 22 Mb** of chromosome arm 2L, past what can be computed live from a VCF, so this LD is precomputed with PLINK and read through [`PlinkLDTabixAdapter`](https://jbrowse.org/jb2-staging/docs/config/plinkldtabixadapter). The sections below build that table, load the same inversion genotyped per mosquito, and read the two together. ## Precompute the LD with PLINK The LD is a file that `plink2 --r2-phased` writes in three steps: thin the variants, correlate them, index the table. `keep.CMgam.txt` is the population, two tab-separated columns of the same sample id, the family/individual pair plink asks for. ```bash # the display uploads n(n-1)/2 cells, and ~800 SNPs across an arm is already at # screen resolution, so keep roughly one variant per 50 kb rather than every # variant the callset has 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 is the haplotype-frequency estimate rather than a correlation # between dosages, which is what the display draws; dprimeabs adds D' beside it # as a magnitude, which is how the display reads a precomputed cell. # --ld-window-r2 0 keeps the uncorrelated pairs. On PLINK 1.9 the pair is one # flag, `--r2 dprime`, and the columns come out 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 comments its own header, which is what `tabix -H` # returns. `sort-bed` is `sort -k1,1 -k2,2n` under LC_ALL=C with that `#` line # kept on top, which is what this table wants too: 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 the display reads one of its two metric columns: ```json addtrack { "type": "LDTrack", "trackId": "ag1000g_2l_cmgam", "name": "Cameroon, both arrangements segregating (r²)", "assemblyNames": ["anoGam3"], "adapter": { "type": "PlinkLDTabixAdapter", "uri": "https://jbrowse.org/demos/popgen/ag1000g_2L_CMgam.vcor.gz" }, "displays": [ { "type": "LDTrackDisplay", "ldMetric": "r2", "useGenomicPositions": true, "showLegend": true, "height": 340 } ] } ``` ## The inversion genotyped per mosquito The same inversion loads as one `` record spanning the breakpoints, genotyped across every mosquito. The [regular multi-sample variant display](https://jbrowse.org/jb2-staging/docs/user_guides/multivariant_track#regular-best-for-full-sv-detail) draws each genotype at the call's true span, so a carrier's row begins and ends at the breakpoints. Those genotypes are what the karyotype lanes in the figure below are: cells shaded by allele dosage, each lane sorted into standard, heterozygous and homozygous-inverted blocks. The `karyotype` column names the three classes by genotype, and so does the legend: `2L+a/2L+a`, `2La/2L+a`, `2La/2La`, the `+` marking the non-inverted arrangement. Load each population as a `VariantTrack` whose adapter carries the samples TSV, with a `LinearMultiSampleVariantDisplay` that orders (`groupBy`) and colors (`colorBy`) its rows by the `karyotype` column: ```json addtrack { "type": "VariantTrack", "trackId": "ag1000g_2la_karyotype_cmgam", "name": "Cameroon, one row per mosquito", "assemblyNames": ["anoGam3"], "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", "groupBy": "karyotype", "colorBy": "karyotype", "referenceDrawingMode": "skip" } ] } ``` What each setting does: - [`groupBy`](https://jbrowse.org/jb2-staging/docs/config/linearmultisamplevariantdisplay/#slot-groupby) keeps the karyotype classes contiguous, so each class reads as one block. - [`referenceDrawingMode`](https://jbrowse.org/jb2-staging/docs/config/linearmultisamplevariantdisplay/#slot-referencedrawingmode) is on its default, `skip`, which fills the lane with the reference color and paints alt cells on top: a solid grey field with the carriers' blocks on it. Rows divide the lane's height between them, so a 300-pixel lane gives each of 297 mosquitoes about a pixel. The display draws a row for every sample in the file, which makes each population its own track. ### The karyotype calls 2La is a cytologically defined arrangement whose breakpoints have been cloned and sequenced ([Sharakhov et al. 2006](https://doi.org/10.1073/pnas.0509683103)), and the call is drawn at that published extent. PCR across the junctions karyotypes single mosquitoes, checked against polytene cytology on field specimens ([White et al. 2007](https://doi.org/10.4269/ajtmh.2007.76.334)). Each mosquito's karyotype here is scored from those tag SNPs, the in-silico method MalariaGEN ships for the current Ag3 release: the mean number of alternate alleles across the tags, rounded into a genotype. That score comes out trimodal with empty space between the peaks, and the [reproduce script](#reproduce-it-end-to-end) prints the histogram and the karyotype breakdown per population. ## The block on the karyotype lanes Four lanes stack in the figure: each population's r² heatmap over its own karyotype lane, one row per mosquito, 297 from Cameroon and 69 from Gabon.
The block's edges land on the published breakpoint coordinates, and on the karyotype lane beneath the heatmap, whose cells are drawn at those same coordinates from a different file. Across the block, markers at opposite ends are about as correlated as neighbouring ones: correlation holds flat with distance over a recombination-suppressed span. Two more things stand out beyond the 2La block itself: - **The second block is _Vgsc_.** At the low-coordinate end of the arm in both panels, it is reddest along the diagonal and pales away below it: the sodium channel whose codon-995 substitutions confer pyrethroid resistance, and which this release was used to survey ([Clarkson et al. 2021](https://doi.org/10.1111/mec.15845)). Gabon shows that block too. - **Gabon's 2La span reads flat.** 64 of its 69 mosquitoes recombine freely across that span, and the MAF floor both files carry drops the variants tagging the 5 heterozygotes. ## What an LD block depends on Four things decide how strongly a block reads, and the [reproduce script](#reproduce-it-end-to-end) prints the numbers behind each: - **Common-variant density.** A sweep that went to fixation leaves few common variants to correlate. Compare density at the locus against a neutral window in the same panel. - **Whether the panel segregates the feature.** Long-range LD inside a candidate span, against an equally distant control, reads as a block only where both arrangements are present. - **Background LD.** A bottlenecked panel renders red across the whole arm, so read the absolute background alongside the inside/outside ratio. - **Number of haplotypes.** r² is a correlation between two biallelic markers, so several haplotypes at one locus fragment the block: each carries a different background, and no single pair of markers tags them all. A soft sweep leaves a patchier block than its strength suggests. ## Metric and allele-frequency floor Two metrics read the same block differently: - **D'** asks whether recombination has been seen between two markers, so it saturates near 1 wherever no recombinant haplotype has turned up. That makes it the read on where recombination stops, and the reproduce script uses it to recover the breakpoints. - **r²** asks how well one marker predicts the other, which also requires the two to be at similar frequency, so it draws the sharper boundary and reads on whether a marker can stand in for another. Switch with [`ldMetric`](https://jbrowse.org/jb2-staging/docs/config/sharedlddisplay/#slot-ldmetric); the script prints both ratios for every panel. Raising the minor allele frequency filter ([`minorAlleleFrequencyFilter`](https://jbrowse.org/jb2-staging/docs/config/sharedlddisplay/#slot-minorallelefrequencyfilter)) thins dense callsets to the common, block-tagging variants. High enough, it reaches the tagging variants themselves and the block fades. ## Reproduce it end to end [`build_ag1000g_ld.sh`](https://github.com/GMOD/jbrowse-components/blob/main/scripts/build_ag1000g_ld.sh) does the whole build for you: - downloads the phased haplotypes - runs each check above and prints the result - builds the tabix-indexed `.vcor.gz` tracks and the per-mosquito karyotype calls - writes a `config.json` opening on the inversion ```bash 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 ``` ## The same karyotype track in Drosophila [](https://jbrowse.org/jb2-staging/docs/tutorials/population_genomics) builds the same one-record karyotype track for an 11 Mb Drosophila inversion. ## See also - [](https://jbrowse.org/jb2-staging/docs/tutorials/ld_human) - [](https://jbrowse.org/jb2-staging/docs/tutorials/population_genomics) - [](https://jbrowse.org/jb2-staging/docs/user_guides/multivariant_track) - [](https://jbrowse.org/jb2-staging/docs/user_guides/variant_track) - [](https://jbrowse.org/jb2-staging/docs/config_guides/variant_track) - [Gallery: variants and populations](https://jbrowse.org/jb2-staging/gallery/#variants) ## References - 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](https://doi.org/10.1101/gr.262790.120) - Clarkson et al. (2021). [The genetic architecture of target-site resistance to pyrethroid insecticides in the African malaria vectors Anopheles gambiae and Anopheles coluzzii](https://doi.org/10.1111/mec.15845) - Love et al. (2019). [In silico karyotyping of chromosomally polymorphic malaria mosquitoes in the Anopheles gambiae complex](https://doi.org/10.1534/g3.119.400445) - Sharakhov et al. (2006). [Breakpoint structure reveals the unique origin of an interspecific chromosomal inversion (2La) in the Anopheles gambiae complex](https://doi.org/10.1073/pnas.0509683103) - White et al. (2007). [Molecular karyotyping of the 2La inversion in Anopheles gambiae](https://doi.org/10.4269/ajtmh.2007.76.334)