XP-CLR

Code to compute xp-clr values as per Chen, Patterson & Reich 2010. This implementation was written due to found bugs in the source of the original tool.

Installation

Clone this git repository into your working directory and cd.

python setup.py install

Also available via conda via the bioconda channel:

conda install xpclr -c bioconda

Use

This is a python module with a convenience script attached. It is designed to run on hdf5 files representing genetic data as generated using scikit-allel. Support is available for VCF and text files in same format as original XPCLR tool, but has not been optimised for this. Support to be added for zarr format soon.

Interface is under development and may change/break in future.

Documentation

usage: xpclr [-h] --out OUT [--format FORMAT] [--input INPUT]
             [--gdistkey GDISTKEY] [--samplesA SAMPLESA] [--samplesB SAMPLESB]
             [--rrate RRATE] [--map MAP] [--popA POPA] [--popB POPB] --chr
             CHROM [--ld LDCUTOFF] [--phased] [--verbose VERBOSE]
             [--maxsnps MAXSNPS] [--minsnps MINSNPS] [--size SIZE]
             [--start START] [--stop STOP] [--step STEP]

Tool to calculate XP-CLR as per Chen, Patterson, Reich 2010

optional arguments:
  -h, --help            show this help message and exit
  --out OUT, -O OUT     output file
  --format FORMAT, -F FORMAT
                        input expected. One of "vcf" (default), "hdf5", or
                        "txt"
  --input INPUT, -I INPUT
                        input file vcf or hdf5
  --gdistkey GDISTKEY   key for genetic position in variants table of hdf5/VCF
  --samplesA SAMPLESA, -Sa SAMPLESA
                        Samples comprising population A. Comma separated list
                        or path to file with each ID on a line
  --samplesB SAMPLESB, -Sb SAMPLESB
                        Samples comprising population B. Comma separated list
                        or path to file with each ID on a line
  --rrate RRATE, -R RRATE
                        recombination rate per base
  --map MAP             input map file as per XPCLR specs
  --popA POPA           filepath to population A genotypes
  --popB POPB           filepath to population A genotypes
  --chr CHROM, -C CHROM
                        Which contig analysis is based on
  --ld LDCUTOFF, -L LDCUTOFF
                        LD cutoff to apply for weighting
  --phased, -P          whether data is phased for more precise r2 calculation
  --verbose VERBOSE, -V VERBOSE
                        How verbose to be in logging. 10=DEBUG, 20=INFO,
                        30=WARN, 40=ERROR, 50=CRITICAL
  --maxsnps MAXSNPS, -M MAXSNPS
                        max SNPs in a window
  --minsnps MINSNPS, -N MINSNPS
                        min SNPs in a window
  --size SIZE           window size in base pairs
  --start START         start base position for windows
  --stop STOP           stop base position for windows
  --step STEP           step size for sliding windows

File formats

hdf5 as generated from vcf by scikit-allel.

Or

vcf via scikit-allel

Or

.geno Space delimited file, containing 0/1s with sample haplotypes as columns, and rows as SNPs.

.map Space delimited file. 6 columns: ID, chromosome, Genetic Distance, Position, REF, ALT.

For examples of these files see the fixture folder used for testing.

Outputs

modelL: The likelihood of the best fitting selection coefficient

nullL: The likelihood of the null model (selection coefficient = 0.0).

sel_coef: The best fitting selection coefficient

nSNPs: Number of SNPs in a window. Suggest to ignore windows where this is small

nSNPs avail: When number of SNPs in a window is greater than the maximum specified. If this is consistently much larger than nSNPs, it suggests your window size is too large.

start, stop: start and stop positions of the windows.

pos_start, pos_stop: actual limits of the used SNPs.

xpclr: The log likelihood ratio of the best fitting model vs the null: 2 * (modelL - nullL)

xpclr_norm: Normalized xpclr ie (xpclr - np.nanmean(xpclr))/np.nanstd(xpclr)