LD_decay

Streamlining the calculation of LD, plotting of LD decay, and decision of LD thresholds for downstream analyses

Setup environment

To setup the conda environment with plink2 and python v3, run the below command: conda create --name LD_decay --file explicit.env.txt

This environment contains plink2, python, numpy/pandas/matplotlib, and dependencies.

Then simply activate the environment with conda activate LD_decay and you should be able to run all scripts within this repo.

In case you can't create the environment with the explicit file, create it with: conda create -n LD_decay -c conda-forge -c bioconda plink2 python pandas numpy matplotlib scipy

Example work flow

Run plinkv2 to calculate LD

Following methods from the speciation genomics tutorial link

$ conda activate LD_decay
$ plink2 --vcf ${vcf_file} --double-id --allow-extra-chr \
	--set-missing-var-ids @:# \
	--mind 0.5 \
	--r2-unphased --ld-window 100 --ld-window-kb 1000 \
	--ld-window-r2 0 \
	--threads ${n_threads} \
	--out ${outname}

This command uses an input vcf file, it allows for extra chromosome names (outside those used in humans), it sets the name of variants to @:# (scaffold:position), it removes individuals with more than 50% missing data, it calculates r2 from unphased genotype data, it sets a minimum window of 100 bp and a maximum window of 1 Mb, it does not filter the output based on r2, it uses n_threads, and produces an output with ${outname}.

Any value presented as ${value} refers to a variable set within the Unix environment. For instance, I might have these three variables set in my script:

vcf_file=/path/to/input.vcf.gz
outname=/path/to/ouput.vcor
threads=8

If needed (e.g., you didn't filter your vcf file yet), add additional filters (e.g., --geno to filter by missing data per variant, --maf to filter variants by minor allele frequency, etc).

This step will produce a .vcor file, which contains a table of pairwise LD values between SNPs. This file will be the input for the LD_plotfit.py script below.

Run plot and fit script

Information on script:

$ python LD_plotfit.py -h
usage: LD_plotfit.py [-h] [--output_file OUTPUT_FILE] [--bin_size BIN_SIZE]
                     [--max_distance MAX_DISTANCE] [--no_fit]
                     [--models {exponential,quadratic_exp} [{exponential,quadratic_exp} ...]]
                     [--slope_threshold SLOPE_THRESHOLD]
                     input_file

Plot distribution of LD distances with binning and curve fitting (exponential
and quadratic exponential).

positional arguments:
  input_file            Path to the input file

options:
  -h, --help            show this help message and exit
  --output_file OUTPUT_FILE
                        Path to the output plot file
  --bin_size BIN_SIZE   Size of distance bins in bp (default: 1000)
  --max_distance MAX_DISTANCE
                        Maximum distance to consider (default: use all data)
  --no_fit              Skip curve fitting
  --models {exponential,quadratic_exp} [{exponential,quadratic_exp} ...]
                        Models to try for curve fitting (default: both
                        exponential and quadratic exponential)
  --slope_threshold SLOPE_THRESHOLD
                        Threshold for slope to consider as negligible (default:
                        0.01 = 1 percent)

$ python scripts/LD_plotfit.py data/matz.vcor --max_distance 100000 --output_file results/matz_LDdecay