# Parsing LD statistics

As described in the multi-population LD section, we are interested in $$\sigma_d^2$$-type statistics, which is the ratio of expectations of $$D^2$$ and $$\pi_2 = p(1-p)q(1-q)$$. Again, $$p$$ and $$q$$ are the allele frequencies at the left and right loci.

To estimate these statistics from data, we take the average of each LD statistic over all pairs of observed (biallelic) SNPs at a given recombination distance, and then divide by the observed $$\pi_2$$ in one of the populations (the “normalizing” population). As described below, we also use a block bootstrapping approach to estimate variances and covariances of observed statistics at each recombination distance, which is used in inference and computing confidence intervals.

## Binned LD decay

To estimate LD decay curves from the data, we bin all pairs of observed SNPs by recombination distance. While we can bin by physical distance (bps) separating SNPs, genetic maps are non-uniform and physical distance does not perfectly correlate with genetic distance at small scales. If we have a recombination map available, it is preferable to compute and compare statistics using that map.

Recombination rate bins are defined by bin edges, which is a list or array with length equal to the number of desired bins plus one. Bin edges should be monotonically increasing, and are thus adjacent without gaps between bins. Thus, bins are defined as semi-open intervals:

import moments.LD
import numpy as np

bin_edges = np.array([0, 1e-6, 1e-5, 1e-4])
print("Bins:")
for b_l, b_r in zip(bin_edges[:-1], bin_edges[1:]):
print(f"[{b_l}, {b_r})")

Bins:
[0.0, 1e-06)
[1e-06, 1e-05)
[1e-05, 0.0001)


There are a few considerations to keep in mind. In practice, very short distances can be problematic, because “non-standard” evolutionary processes can distort allele frequency correlations for tightly linked loci. For example, our evolutionary model does not include multi-nucleotide mutations [Harris2014] or gene conversion [Ardlie2001], both of which operate at short distances.

Thus, when working with real data we recommend omitting bins of very short recombination distances. In practice, we typically drop bins with length less than $$r=5\times 10^{-6}$$, which corresponds to roughly a few hundred bp on average in humans.

## Parsing from a VCF

The primary function of the Parsing module is computing LD statistics from an input VCF. There are a number of options available, but the primary inputs are the path to the VCF file and the bins of distances separating loci. Typically, we work in recombination distance, in which case a recombination map is also required. If we do not have a recombination map available, we can bin by base pair distances instead.

The function moments.LD.Parsing.compute_ld_statistics returns a dictionary with the bins, returned statistics, populations, and sums of each statistic over the provided bins. For example:

r_bins = np.logspace(-6, -3, 7)
ld_stats = moments.LD.Parsing.compute_ld_statistics(
vcf_path, r_bins=r_bins, rec_map_file=map_path)


### Using a recombination map

The input recombination map is specified as a text file, with the first column giving positions along the chromosome and additional column(s) defining the cumulative map(s), typically in units of cM. The header line is “Pos Map1 Map2 (etc)”, and we can use any map in the file by specifying the map_name. If no map name is given, or the specified map name does not match a genetic map in the header, we use the map in the first column.

Typically, maps are given in units of centi-Morgans, and the default behavior is to assume cM units. If the map is given in units of Morgans, we need to set cM=False.

### Populations and pop-file

We often have data from more than one population, so we need to be able to specify which samples in the VCF correspond to which populations. This is handled by passing a file that assigns each sample to a population. For example, the population file is written as

sample  pop
sid_0   pop_A
sid_1   pop_B
sid_2   pop_A
sid_3   pop_A
sid_4   pop_B
...


Then to include the population information in the function, we also pass a list of the populations to compute statistics for. Samples from omitted populations are dropped from the data.

pops = ["pop_A", "pop_B"]
ld_stats = ld_stats = moments.LD.Parsing.compute_ld_statistics(
vcf_path,
r_bins=r_bins,
rec_map_file=map_path,
pop_file=pop_file_path,
pops=pops
)


### Masking and using bed files

If there are multiple chromosomes or contigs included in the VCF, we specify which chromosome to compute statistics for by setting the chromosome flag. We can also subset a chromosome by including a bed file, which will filter out all SNPs that fall outside the region intervals given in the bed file. Bed files have the format {chrom}\t{left_pos}\t{right_pos}, which defines a semi-open interval. The path to the bed file is provided with the bed_file argument.

### Computing a subset of statistics

Sometimes we may wish to only compute a subset of possible LD statistics. By default, the parsing function computes all statistics possible for the number of populations provided. Instead, we can specify the stats_to_compute, which is a list (of length 2) of lists. The first list are the LD statistics to return, and the second list has the heterozygosity statistics to return. Statistic names follow the convention in moments.LD.Util.moment_names(num_pops), and should be formatted accordingly.

### Phased vs unphased data

We can compute LD statistics from either phased or unphased data. The default behavior is to assume that phasing is unknown, and use_genotypes is True by default. If we want to compute LD using phased data, we set use_genotypes=False, and parsing uses phased haplotypes instead. In general, phasing errors can bias LD statistics, sometimes significantly, and using genotypes instead of haplotypes only slightly increases uncertainty in most cases. Therefore, we usually recommend leaving use_genotypes=True.

## Computing averages and covariances over regions

From moments.LD.Parsing.compute_ld_statistics(), we get LD statistic sums from the regions in a VCF, perhaps constrained by a bed file. Our strategy is to divide our data into some large number of roughly equally sized chunks, for example 500 regions across all 22 autosomes in human data. We then compute LD statistics independently for each region (it helps to parallelize that step, using a compute cluster). From those outputs, we can then compute average statistics genome-wide, as well as covariances of statistics within each bin. Those covariances are needed to be able to compute likelihoods and run optimization.

The outputs of compute_ld_statistics are compiled in a dictionary, where the keys are unique region identifiers, and items the outputs of that function. For example:

region_stats = {
0: moments.LD.Parsing.compute_ld_statistics(VCF, bed_file="region_0.bed", ...),
1: moments.LD.Parsing.compute_ld_statistics(VCF, bed_file="region_1.bed", ...),
2: moments.LD.Parsing.compute_ld_statistics(VCF, bed_file="region_2.bed", ...),
...
}


Mean and variance-covariance matrices are computed by calling bootstrap_data, passing the region statistics dictionary, and optionally the index of the population to normalize $$\sigma_d^2$$ statistics by. By default, the normalizing population is the first (index 0).

mv = moments.LD.Parsing.bootstrap_data(region_stats)


mv contains the bins, statistics, and populations, as well as lists of mean statistics and variance-covariance matrices. This data can then be directly compared to model expectations and used in inference.

## Example

Using msprime [Kelleher2016], we’ll simulate some data under an isolation-with-migration (IM) model and then compute LD and heterozygosity statistics using the LD.Parsing methods. First, the simulation will use the demes-msprime interface, which are then written as a VCF.

The YAML-file specifying the model is

description: A simple isolation-with-migration model
time_units: generations
demes:
- name: anc
epochs: [{start_size: 10000, end_time: 1500}]
- name: deme0
ancestors: [anc]
epochs: [{start_size: 2000}]
- name: deme1
ancestors: [anc]
epochs: [{start_size: 20000}]
migrations:
- demes: [deme0, deme1]
rate: 1e-4


And we use msprime to simulate 1Mb of data, using a constant recombination and mutation rate.

import msprime
import demes
import os

# set up simulation parameters
L = 1e6
u = r = 1.5e-8
n = 10

demog = msprime.Demography.from_demes(g)

trees = msprime.sim_ancestry(
{"deme0": n, "deme1": n},
demography=demog,
sequence_length=L,
recombination_rate=r,
random_seed=321,
)

trees = msprime.sim_mutations(trees, rate=u, random_seed=123)

with open("data/im-parsing-example.vcf", "w+") as fout:
trees.write_vcf(fout)


This simulation had 10 diploid individuals per population, and msprime/tskit writes their IDs as tsk_0, tsk_1, etc:

##fileformat=VCFv4.2
##source=tskit 0.3.5
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=1,length=1000000>
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	tsk_0	tsk_1	tsk_2	tsk_3	tsk_4	tsk_5	tsk_6	tsk_7	tsk_8	tsk_9	tsk_10	tsk_11	tsk_12	tsk_13	tsk_14	tsk_15	tsk_16	tsk_17	tsk_18	tsk_19
1	221	.	A	T	.	PASS	.	GT	1|1	1|0	1|1	1|1	0|1	0|1	0|0	1|1	0|1	1|1	1|1	1|1	1|1	1|1	1|1	1|1	1|1	1|1	1|1	1|1
1	966	.	A	C	.	PASS	.	GT	0|0	0|0	0|0	0|0	0|0	0|0	0|0	1|1	0|0	0|0	0|0	0|0	0|1	0|0	0|0	0|0	0|0	0|0	0|1	1|0
1	1082	.	G	A	.	PASS	.	GT	0|0	0|1	0|0	0|0	1|0	1|0	1|1	0|0	1|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0
1	1133	.	G	T	.	PASS	.	GT	0|0	0|1	0|0	0|0	1|0	1|0	1|1	0|0	1|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0


To parse this data, we need the file that maps samples to populations and the recombination map file (the total map length is found by $$1 \times 10^{6} \text{ bp} \times 1.5 \times 10^{-8} \text{ M/bp} \times 100 \text{ cM/M}$$):

sample  pop
tsk_0   deme0
tsk_1   deme0
tsk_2   deme0
tsk_3   deme0
tsk_4   deme0
tsk_5   deme0
tsk_6   deme0
tsk_7   deme0
tsk_8   deme0
tsk_9   deme0
tsk_10  deme1
tsk_11  deme1
tsk_12  deme1
tsk_13  deme1
tsk_14  deme1
tsk_15  deme1
tsk_16  deme1
tsk_17  deme1
tsk_18  deme1
tsk_19  deme1

Pos Map(cM)
0   0
1000000 1.5


With all this, we can now compute LD based on recombination distance bins:

r_bins = np.array(
[0, 1e-6, 2e-6, 5e-6, 1e-5, 2e-5, 5e-5, 1e-4, 2e-4, 5e-4, 1e-3]
)

vcf_file = "data/im-parsing-example.vcf"
map_file = "data/im-parsing-example.map.txt"
pop_file = "data/im-parsing-example.samples.pops.txt"
pops = ["deme0", "deme1"]
ld_stats = moments.LD.Parsing.compute_ld_statistics(
vcf_file,
rec_map_file=map_file,
pop_file=pop_file,
pops=["deme0", "deme1"],
r_bins=r_bins,
report=False,
)


The output, ld_stats, is a dictionary with the keys bins, stats, pops, and sums. To get the average statistics over multiple regions (here, we only have a single region that we simulated), we use means_from_region_data:

means = moments.LD.Parsing.means_from_region_data(
{0: ld_stats}, ld_stats["stats"], norm_idx=0
)


This provides $$\sigma_d^2$$-type statistics relative to pi_2 in deme0, and relative heterozygosities (also relative to deme0). These statistics were computed from only a single relatively small region, so they will be quite noisy. But we can still compare to expectations under the input IM demographic model.

import demes

y = moments.Demes.LD(
g,
sampled_demes=["deme0", "deme1"],
rho=4 * g["anc"].epochs[0].start_size * r_bins,
)

# stats are computed at the bin edges - average to get midpoint estimates
y = moments.LD.LDstats(
[(y_l + y_r) / 2 for y_l, y_r in zip(y[:-2], y[1:-1])] + [y[-1]],
num_pops=y.num_pops,
pop_ids=y.pop_ids,
)

# plot LD decay curves for some statistics
moments.LD.Plotting.plot_ld_curves_comp(
y,
means[:-1],
[],
rs=r_bins,
stats_to_plot=[
["DD_0_0", "DD_0_1", "DD_1_1"],
["Dz_0_0_0", "Dz_0_1_1", "Dz_1_1_1"],
["pi2_0_0_1_1", "pi2_0_1_0_1", "pi2_1_1_1_1"]
],
labels=[[r"$D_0^2$", r"$D_0 D_1$", r"$D_1^2$"],
[r"$Dz_{0,0,0}$", r"$Dz_{0,1,1}$", r"$Dz_{1,1,1}$"],
[r"$\pi_{2;0,0,1,1}$", r"$\pi_{2;0,1,0,1}$", r"$\pi_{2;1,1,1,1}$"]
],
plot_vcs=False,
fig_size=(8, 3),
show=True,
)


### Bootstrapping over multiple regions

Normally, we’ll want more data than from a single 1Mb region to compute averages and variances of statistics. Using the same approach as the above example, ld_stats for 100 replicates we computed (see example in the moments repository here). From this, each replicate set of statistics were placed in a dictionary, as rep_stats = {0: ld_stats_0, 1: ld_stats_1, ..., 99: ld_stats_99}. This dictionary can then be used to compute means and covariances of statistics.

mv = moments.LD.Parsing.bootstrap_data(ld_stats)


By simulating more data, the LD decay curves are much less noisy, and by simulating multiple replicates, we also compute the variance-covariance matrices for each bin and can include standard errors in the plots.

# plot LD decay curves for some statistics
moments.LD.Plotting.plot_ld_curves_comp(
y,
mv["means"][:-1],
mv["varcovs"][:-1],
rs=r_bins,
stats_to_plot=[
["DD_0_0", "DD_0_1", "DD_1_1"],
["Dz_0_0_0", "Dz_0_1_1", "Dz_1_1_1"],
["pi2_0_0_1_1", "pi2_0_1_0_1", "pi2_1_1_1_1"]
],
labels=[[r"$D_0^2$", r"$D_0 D_1$", r"$D_1^2$"],
[r"$Dz_{0,0,0}$", r"$Dz_{0,1,1}$", r"$Dz_{1,1,1}$"],
[r"$\pi_{2;0,0,1,1}$", r"$\pi_{2;0,1,0,1}$", r"$\pi_{2;1,1,1,1}$"]
],
plot_vcs=True,
fig_size=(8, 3),
show=True,
)


Note

The means-covariances data is required for inference using LD statistics. In Inferring demography with LD, we’ll use the same mv data dictionary to refit the IM model as an example.

## LD statistics in genotype blocks

moments.LD.Parsing also includes some functions for computing LD from genotype (or haplotype) blocks. Genotype blocks are arrays of shape $$L\times n$$, where L is the number of loci and n is the sample size. We assume a single population, and so we compute $$D^2$$, $$Dz$$, $$\pi_2$$, and $$D$$, either pairwise or averaged over all pairwise comparisons.

If we have a genotype matrix containing n diploid samples, genotypes are coded as 0, 1, and 2, and we set genotypes=True. If we have a haplotype matrix with data from n haploid copies, genotypes are coded as 0 and 1 only, and we set genotypes=False.

For example, given a single genotype matrix, we compute all pairwise statistics and average statistics as shown below:

L = 10
n = 5
G = np.random.randint(3, size=L * n).reshape(L, n)

# all pairwise comparisons:
D2_pw, Dz_pw, pi2_pw, D_pw = moments.LD.Parsing.compute_pairwise_stats(G)

# averages:
D2_ave, Dz_ave, pi2_ave, D_ave = moments.LD.Parsing.compute_average_stats(G)


Similarly, we can compute the pairwise or average statistics between two genotype matrices. The matrices can have differing number of loci, but they must have the same number of samples, as the genotype matrices are assumed to come from different regions within the same samples.

L2 = 12
n = 5

G2 = np.random.randint(3, size=L2 * n).reshape(L2, n)

# all pairwise comparisons:
D2_pw, Dz_pw, pi2_pw, D_pw = moments.LD.Parsing.compute_pairwise_stats_between(G, G2)

# averages:
D2_ave, Dz_ave, pi2_ave, D_ave = moments.LD.Parsing.compute_average_stats_between(G, G2)


Note

Computing LD in genotype blocks uses C-extensions that are not built by default, so are only available if these are built when compiling the C-extensions. In order to use these methods, we need to build these extensions using the --ld_extensions flag, as python setup.py build_ext --ld_extensions -i.

## References

Ardlie2001

Ardlie, Kristin, et al. “Lower-than-expected linkage disequilibrium between tightly linked markers in humans suggests a role for gene conversion.” The American Journal of Human Genetics 69.3 (2001): 582-589.

Harris2014

Harris, Kelley, and Rasmus Nielsen. “Error-prone polymerase activity causes multinucleotide mutations in humans.” Genome research 24.9 (2014): 1445-1454.

Kelleher2016

Kelleher, Jerome, Alison M. Etheridge, and Gilean McVean. “Efficient coalescent simulation and genealogical analysis for large sample sizes.” PLoS computational biology 12.5 (2016): e1004842.