PCAngsd

PCAngsd is a program that estimates the covariance matrix and individual allele frequencies for low-depth next-generation sequencing (NGS) data in structured/heterogeneous populations using principal component analysis (PCA) to perform multiple population genetic analyses using genotype likelihoods. Since version 0.98, PCAngsd was re-written to be based on Cython for computational bottlenecks and parallelization.

The main method was published in 2018 and can be found here: [1]

The HWE test was published in 2019 and can be found here: [2]

Simulated low depth NGS data of 3 populations

Overview

Framework for analyzing low-depth next-generation sequencing (NGS) data in heterogeneous/structured populations using principal component analysis (PCA). Population structure is inferred by estimating individual allele frequencies in an iterative approach using a truncated SVD model. The covariance matrix is estimated using the estimated individual allele frequencies as prior information for the unobserved genotypes in low-depth NGS data.

The estimated individual allele frequencies can further be used to account for population structure in other probabilistic methods. PCAngsd can perform the following analyses:

• Covariance matrix
• Admixture estimations
• Inbreeding coefficients (both per-individual and per-site)
• HWE test
• Genome-wide selection scan
• Genotype calling
• Estimate NJ tree of samples

Older versions of PCAngsd can be found here [3].

Download and Installation

PCAngsd should work on all platforms meeting the requirements but server-side usage is highly recommended. Installation has only been tested on Linux systems.

Get PCAngsd and build

git clone https://github.com/Rosemeis/pcangsd.git
cd pcangsd/
python setup.py build_ext --inplace

Install dependencies:

The required set of Python packages are easily installed using the pip command and the 'requirements.txt file' included in the 'pcangsd' folder.

pip install --user -r requirements.txt

Quick start

PCAngsd is used by running the main caller file pcangsd.py. To see all available options use the following command:

python pcangsd.py -h

# Genotype likelihoods using 64 threads
python pcangsd.py -beagle input.beagle.gz -out output -threads 64

# PLINK files (using file-prefix, *.bed, *.bim, *.fam)
python pcangsd.py -beagle input.plink -out output -threads 64

PCAngsd accepts either genotype likelihoods in Beagle format or PLINK genotype files. Beagle files can be generated from BAM files using ANGSD. For inference of population structure in genotype data with non-random missigness, we recommend our EMU software that performs accelerated EM-PCA, however with fewer functionalities than PCAngsd (#soon).

PCAngsd will mostly output files in binary Numpy format (.npy) with a few exceptions. In order to read files in python:

import numpy as np
C = np.genfromtxt("output.cov") # Reads in estimated covariance matrix (text)
D = np.load("output.selection.npy") # Reads PC based selection statistics

R can also read Numpy matrices using the "RcppCNPy" R library:

library(RcppCNPy)
C <- as.matrix(read.table("output.cov")) # Reads in estimated covariance matrix
D <- npyLoad("output.selection.npy") # Reads PC based selection statistics

An example of generating genotype likelihoods in ANGSD and output them in the required Beagle text format.

./angsd -GL 2 -out input -nThreads 4 -doGlf 2 -doMajorMinor 1 -doMaf 2 -SNP_pval 1e-6 -bam bam.filelist

Tutorial

Please refer to the tutorial's page [4]

Options

# See all options in PCAngsd
python pcangsd.py -h

General usage

-beagle [Beagle file]

Input file of genotype likelihoods in Beagle format (.beagle.gz).

-filter [Text file]

Input file of 1's or 0's whether to keep individuals or not.

-plink [Prefix for binary PLINK files]

Path to PLINK files using their ONLY prefix (.bed, .bim, .fam).

-plink_error [float]

Incorporate errors into genotypes by specifying rate as argument.

-minMaf [float]

Minimum minor allele frequency threshold. (Default: 0.05)

-maf_iter [int]

Maximum number of EM iterations for computing the population allele frequencies (Default: 200).

-maf_tole [float]

Tolerance value in EM algorithm for population allele frequencies estimation (Default: 1e-4).

-iter [int]

Maximum number of iterations for estimation of individual allele frequencies (Default: 100).

-tole [float]

Tolerance value for update in estimation of individual allele frequencies (Default: 1e-5).

-hwe [.lrt.npy file]

Input file of LRT binary file from previous PCAngsd run to filter based on HWE.

-hwe_tole [float]

Threshold for HWE filtering of sites.

-e [int]

Manually select the number of eigenvalues to use in the modelling of individual allele frequencies (Default: Automatically tested using MAP test).

-pi [.pi.npy file]

Load previous estimation of individual allele frequencies to skip covariance estimation.

-maf_save

Choose to save estimated population allele frequencies (Binary). Numpy format (.npy).

-pi_save

Choose to save estimated individual allele frequencies (Binary). Numpy format (.npy). Can be used with the '-pi' command.

-dosage_save

Choose to save estimated genotype dosages (Binary). Numpy format (.npy).

-post_save

Choose to save the posterior genotype probabilities. Beagle format (.beagle).

-sites_save

Choose to save the kept sites after filtering which is useful for downstream analysis. Outputs a file of 1's and 0's for keeping a site or not, respectively.

-threads [int]

Specify the number of thread(s) to use (Default: 1).

-out [output prefix]

Fileprefix for all output files created by PCAngsd (Default: "pcangsd").

Selection

Perform PC-based genome-wide selection scans using posterior expectations of the genotypes (genotype dosages):

-selection

Using an extended model of FastPCA. Performs a genome-wide selection scan along all significant PCs. Outputs the selection statistics and must be converted to p-values by user. Each column reflect the selection statistics along a tested PC and they are χ²-distributed with 1 degree of freedom.

-pcadapt

Using an extended model of pcadapt. Performs a genome-wide selection scan across all significant PCs. Outputs the z-scores and must be converted to test statistics with the provided script 'pcangsd/scripts/pcadapt.R', and the test statistics are χ²-distributed with K degree of freedom.

-snp_weights

Output the SNP weights of the significant K eigenvectors.

Inbreeding

-inbreedSites

Estimate per-site inbreeding coefficients accounting for population structure and perform likehood ratio test for detecting sites deviating from HWE [5].

-inbreedSamples

Estimate per-individual inbreeding coefficients accounting for population structure which is based on an extension of ngsF for structured populations.

-inbreed_iter [int]

Maximum number of iterations for inbreeding EM algorithm. (Default: 200)

-inbreed_tole [float]

Tolerance value for inbreeding EM algorithm in estimating inbreeding coefficients. (Default: 1e-4)

Call genotypes

Genotypes can be called from posterior genotype probabilities by incorporating the individual allele frequencies as prior information.

-geno [float]

Call genotypes with defined threshold.

-genoInbreed [float]

Call genotypes with defined threshold also taking inbreeding into account. '-inbreedSamples' must also be called for using this option.

Admixture

Individual admixture proportions and ancestral allele frequencies can be estimated assuming K ancestral populations using an accelerated mini-batch NMF method.

-admix

Toggles admixture estimations. Estimates admixture proportions and ancestral allele frequencies.

-admix_K [int]

Not recommended. Override the number of ancestry components (K) to use, instead of using K=e-1.

-admix_iter [int]

Maximum number of iterations for admixture estimations using NMF. (Default: 200)

-admix_tole [float]

Tolerance value for update in admixture estimations using NMF. (Default: 1e-5)

-admix_alpha [float

Specify alpha (sparseness regularization parameter). (Default: 0)

-admix_auto [float]

Enable automatic search for optimal alpha using likelihood measure, by giving soft upper search bound of alpha.

-admix_seed [int]

Specify seed for random initializations of factor matrices in admixture estimations.

Citation

Our methods for inferring population structure have been published in GENETICS:

Inferring Population Structure and Admixture Proportions in Low Depth NGS Data

Our method for testing for HWE in structured populations has been published in Molecular Ecology Resources:

Testing for Hardy‐Weinberg Equilibrium in Structured Populations using Genotype or Low‐Depth NGS Data