PCAngsd: Difference between revisions

From software
Jump to navigation Jump to search
No edit summary
Line 1: Line 1:
This page contains information about the program PCAngsd, which estimates the covariance matrix for low depth NGS data in an iterative procedure based on genotype likelihoods and is able to perform multiple population genetic analyses in heterogeneous populations. Based on population structure inference, PCAngsd is able to estimate individual allele frequencies. These individual allele frequencies can be used in various population genetic methods for heterogeneous populations, such that PCAngsd can perform PCA (estimate covariance matrix), call genotypes, estimate individual admixture proportions, estimate inbreeding coefficients (per-individual and per-site) and perform a genome selection scan using principal components. The entire program is written in Python 2.7 and is multithreaded to take advantage of several CPUs.
This page contains information about PCAngsd, which estimates the covariance matrix for low-depth NGS data in an iterative procedure based on genotype likelihoods and is able to perform multiple population genetic analyses in structured/heterogeneous populations. Based on iterative population structure inference, PCAngsd estimates individual allele frequencies. These individual allele frequencies can be used in various analyses to account for population structure in such a way that PCAngsd can perform PCA (estimate covariance matrix), call genotypes, estimate individual admixture proportions, estimate inbreeding coefficients (per-site and per-individual) and perform genomic selection scan using principal components. The entire program is written in Python and is multithreaded to take advantage of several CPUs.
 
Since version ''0.98'', PCAngsd has been re-written in Cython for computational bottlenecks and parallelization.


[[File:Pcangsd_admix.gif|frame]]
[[File:Pcangsd_admix.gif|frame]]
Line 10: Line 12:
https://github.com/Rosemeis/pcangsd
https://github.com/Rosemeis/pcangsd


See github for Installation
<pre>
git clone https://github.com/Rosemeis/pcangsd.git
</pre>


PCAngsd should work on all platforms meeting the requirements but server-side usage is highly recommended.
See github for more information regarding 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.


==Quick start==
==Quick start==
Line 26: Line 31:


# Estimate covariance matrix and inbreeding coefficients
# Estimate covariance matrix and inbreeding coefficients
python pcangsd.py -beagle data.beagle.gz -inbreed 2 -o test -threads 10
python pcangsd.py -beagle data.beagle.gz -inbreed 1 -o test -threads 10


# Estimate covariance matrix and perform selection scan
# Estimate covariance matrix and perform selection scan
python pcangsd.py -beagle data.beagle.gz -selection 1 -o test -threads 10
python pcangsd.py -beagle data.beagle.gz -selection -o test -threads 10
</pre>
</pre>


Line 43: Line 48:
=Using PCAngsd=
=Using PCAngsd=


All the different options in PCAngsd are listed here. PCAngsd will always compute the covariance matrix, where it uses principal components to estimate individual allele frequencies in an iterative procedure. The estimated individual allele frequencies will then be used in any of the other specified options of PCAngsd.
All the different options in PCAngsd are listed here. PCAngsd will always compute and output the covariance matrix, where it uses principal components to estimate individual allele frequencies in an iterative procedure. The estimated individual allele frequencies will then be used in any of the other specified options of PCAngsd.


==Estimation of individual allele frequencies==
==Estimation of individual allele frequencies==
; -beagle [Beagle filename]
; -beagle [Beagle filename]
Input file of genotype likelihoods in Beagle format (.beagle.gz).
Input file of genotype likelihoods in Beagle format (.beagle.gz).
; -indf [Individual allele frequencies filename]
Input file of individual allele frequencies (binary). Numpy format (.indf.npy).
; -plink [Prefix for binary PLINK files]
; -plink [Prefix for binary PLINK files]
Path to PLINK files using their prefix. (.bed, .bim, .fam)
Path to PLINK files using their prefix (.bed, .bim, .fam).
; -epsilon [float]
; -plink_error [float]
Include error assumption in PLINK genotypes. (Default: 0.00)
Incorporate error model for PLINK genotypes.
; -minMaf [float]
; -minMaf [float]
Minimum minor allele frequency threshold. (Default: 0.05)
Minimum minor allele frequency threshold. (Default: 0.05)
; -threads [int]
Specify the number of thread(s) to use. (Default: 1)
; -iter [int]
; -iter [int]
Maximum number of iterations for estimation of individual allele frequencies. (Default: 100)
Maximum number of iterations for estimation of individual allele frequencies (Default: 100).
; -tole [float]
; -tole [float]
Tolerance value for update in estimation of individual allele frequencies. (Default: 1e-5)
Tolerance value for update in estimation of individual allele frequencies (Default: 1e-5).
; -maf [int]
; -maf_iter [int]
Maximum number of EM iterations for computing the population allele frequencies. (Default: 200)
Maximum number of EM iterations for computing the population allele frequencies (Default: 100).
; -maf_tole [float]
; -maf_tole [float]
Tolerance value in EM algorithm for population allele frequencies estimation. (Default: 1e-4)
Tolerance value in EM algorithm for population allele frequencies estimation (Default: 1e-4).
; -e [int]
; -e [int]
Manually select the number of eigenvalues to use in the modelling of individual allele frequencies. (Default: Automatically tested using MAP test)
Manually select the number of eigenvalues to use in the modelling of individual allele frequencies (Default: Automatically tested using MAP test).
; -o [prefix]
; -o [prefix]
Set the prefix for all output files created by PCAngsd (Default: "pcangsd").
Set the prefix for all output files created by PCAngsd (Default: "pcangsd").
; -indf_save
; -indf_save
Choose to save estimated individual allele frequencies (Binary). Numpy format (.npy).
Choose to save estimated individual allele frequencies (Binary). Numpy format (.npy).
; -expg_save
; -dosage_save
Choose to save estimated genotype dosages (Binary). Numpy format (.npy).
Choose to save estimated genotype dosages (Binary). Numpy format (.npy).
; -sites_save
; -sites_save
Choose to save the marker IDs after performing filtering using population allele frequencies. Useful for especially selection scans and per-site inbreeding coefficients.
Choose to save the marker IDs after performing filtering using population allele frequencies. Useful for especially selection scans and per-site inbreeding coefficients.
; -post_save
Choose to save the posterior genotype probabilities. Beagle format (.beagle).
; -threads [int]
Specify the number of thread(s) to use (Default: 1).


==Call genotypes==
==Call genotypes==
Line 81: Line 86:


; -geno [float]
; -geno [float]
Call genotypes with defined threshold.
Call genotypes with defined threshold (Binary). Numpy format (.npy).
; -genoInbreed [float]
; -genoInbreed [float]
Call genotypes with defined threshold also taking inbreeding into account. '''-inbreed [int]''' is required, since individual inbreeding coefficients must have been estimated prior to calling genotypes using that information.
Call genotypes with defined threshold also taking inbreeding into account. '''-inbreed [int]''' is required, since individual inbreeding coefficients must have been estimated prior to calling genotypes using that information (Binary). Numpy format (.npy).


==Admixture==
==Admixture==
Line 89: Line 94:


; -admix
; -admix
Toggles admixture estimations.
Toggles admixture estimations. Individual ancestry proportions are saved (Binary). Numpy format (.npy).
; -admix_alpha [int-list]
; -admix_alpha [float-list]
Specify alpha (sparseness regularization parameter). Can be specified as a sequence to try several alpha's in a single run. Fully compatible with -admix_seed. (Default: 0)
Specify alpha (sparseness regularization parameter). Can be specified as a sequence to try several alpha's in a single run (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-list]
; -admix_seed [int-list]
Specify seed for random initializations of factor matrices in admixture estimations. Can be specified as a sequence to try several different seeds in a single run. Fully compatible with -admix_alpha.
Specify seed for random initializations of factor matrices in admixture estimations. Can be specified as a sequence to try several different seeds in a single run.
; -admix_K [int]
; -admix_K [int]
Not recommended. Manually specify the number of ancestral populations to use in admixture estimations (overrides number chosen from '''-e'''). Structure explained by individual allele frequencies may therefore not reflect the manually chosen K. It is recommended to adjust '''-e''' instead of '''-admix_K'''.
Not recommended. Manually specify the number of ancestral populations to use in admixture estimations (overrides number chosen from '''-e'''). Structure explained by individual allele frequencies may therefore not reflect the manually chosen K. It is recommended to adjust '''-e''' instead of '''-admix_K'''.
; -admix_iter [int]
; -admix_iter [int]
Maximum number of iterations for admixture estimations using NMF. (Default: 100)
Maximum number of iterations for admixture estimations using NMF. (Default: 200)
; -admix_tole [float]
; -admix_tole [float]
Tolerance value for update in admixture estimations using NMF. (Default: 5e-5)
Tolerance value for update in admixture estimations using NMF. (Default: 1e-5)
; -admix_batch [int]
; -admix_batch [int]
Specify the number of batches to use in NMF method. (Default: 5)
Specify the number of batches to use in NMF method. (Default: 10)
; -admix_save
; -admix_save
Choose to save the population-specific allele frequencies (Binary). Numpy format (.npy).
Choose to save the population-specific allele frequencies (Binary). Numpy format (.npy).


==Inbreeding==
==Inbreeding==
Per-individual inbreeding coefficients incorporating population structure can be computed using three different methods. However, -inbreed 2 is recommended for low depth cases.
Per-individual inbreeding coefficients incorporating population structure can be computed using three different methods. However, -inbreed 1 is recommended for low depth cases.


; -inbreed 1
; -inbreed 1
A maximum likelihood estimator computed by an EM algorithm. Only allows for F-values between 0 and 1. Based on [https://www.cambridge.org/core/journals/genetics-research/article/maximum-likelihood-estimation-of-individual-inbreeding-coefficients-and-null-allele-frequencies/2DEBA0C0C2B92DF0EE89BD27DFCAD3FB].
Simple estimator computed by an EM algorithm. Allows for F-values between -1 and 1. Based on [http://genome.cshlp.org/content/23/11/1852.full ngsF].
; -inbreed 2
; -inbreed 2
Simple estimator also computed by an EM algorithm. Allows for F-values between -1 and 1. Based on [http://genome.cshlp.org/content/23/11/1852.full ngsF].
A maximum likelihood estimator also computed by an EM algorithm. Only allows for F-values between 0 and 1. Based on [https://www.cambridge.org/core/journals/genetics-research/article/maximum-likelihood-estimation-of-individual-inbreeding-coefficients-and-null-allele-frequencies/2DEBA0C0C2B92DF0EE89BD27DFCAD3FB].
; -inbreed 3
; -inbreed 3
Estimator using an estimated kinship matrix. Allows for F-values between -1 and 1. Based on [http://www.cell.com/ajhg/abstract/S0002-9297(15)00493-0 PC-Relate].   
Estimator using an estimated kinship matrix. Allows for F-values between -1 and 1. Based on [http://www.cell.com/ajhg/abstract/S0002-9297(15)00493-0 PC-Relate].   
Line 126: Line 133:
Use likelihood ratio tests (.lrt.sites.gz) generated from '''-inbreedSites''' to filter out variable sites using a given threshold for HWE test p-value:
Use likelihood ratio tests (.lrt.sites.gz) generated from '''-inbreedSites''' to filter out variable sites using a given threshold for HWE test p-value:


; -HWE_filter [LRT filename]
; -hwe [LRT filename]
; -HWE_tole [float]
; -hwe_tole [float]
Tolerance value for HWE test. (Default: 1e-6)
Tolerance value for HWE test. (Default: 1e-6)


==Selection==
==Selection==
A genome selection scan can be computed using two different methods based on posterior expectations of the genotypes (genotype dosages):
A genome selection scan can be computed based on posterior expectations of the genotypes (genotype dosages):


; -selection 1
; -selection
Using an extended model of [http://www.cell.com/ajhg/abstract/S0002-9297(16)00003-3 FastPCA]. Performs a genome 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.
Using an extended model of [http://www.cell.com/ajhg/abstract/S0002-9297(16)00003-3 FastPCA]. Performs a genome 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.
; -selection 2
Using an extended model of [http://onlinelibrary.wiley.com/doi/10.1111/1755-0998.12592/abstract PCAdapt]. Outputs the selection statistics and must be converted to p-values by user. Selection statistics are χ²-distributed with '''-e''' degrees of freedom (number of significant eigenvectors).


==Relatedness==
==Relatedness==
'''Work in progress...'''
Estimate kinship matrix based on method Based on [http://www.cell.com/ajhg/abstract/S0002-9297(15)00493-0 PC-Relate]:
Estimate kinship matrix based on method Based on [http://www.cell.com/ajhg/abstract/S0002-9297(15)00493-0 PC-Relate]:


; -kinship
; -kinship
Automatically estimated if '''-inbreed 3''' has been selected.
Automatically estimated if '''-inbreed 3''' has been selected.
Remove related individuals based on kinhsip matrix of previous run:
; -relate [Kinship filename]
; -relate_tole [float]
Threshold for kinship coefficients for removing individuals (Default: 0.0625).


=Citation=
=Citation=
Our methods for inferring population structure have been published in GENETICS:
Our methods for inferring population structure have been published in GENETICS:
[http://www.genetics.org/content/early/2018/08/21/genetics.118.301336 Inferring Population Structure and Admixture Proportions in Low Depth NGS Data]


[http://www.genetics.org/content/early/2018/08/21/genetics.118.301336 Inferring Population Structure and Admixture Proportions in Low Depth NGS Data]
Our method for estimating per-site inbreeding sites and testing for HWE has been published in Molecular Ecology Resources:
[https://onlinelibrary.wiley.com/doi/abs/10.1111/1755-0998.13019 Testing for Hardy‐Weinberg Equilibrium in Structured Populations using Genotype or Low‐Depth NGS Data]

Revision as of 09:07, 26 July 2019

This page contains information about PCAngsd, which estimates the covariance matrix for low-depth NGS data in an iterative procedure based on genotype likelihoods and is able to perform multiple population genetic analyses in structured/heterogeneous populations. Based on iterative population structure inference, PCAngsd estimates individual allele frequencies. These individual allele frequencies can be used in various analyses to account for population structure in such a way that PCAngsd can perform PCA (estimate covariance matrix), call genotypes, estimate individual admixture proportions, estimate inbreeding coefficients (per-site and per-individual) and perform genomic selection scan using principal components. The entire program is written in Python and is multithreaded to take advantage of several CPUs.

Since version 0.98, PCAngsd has been re-written in Cython for computational bottlenecks and parallelization.

Simulated low depth NGS data of 3 populations


Download

The program can be downloaded from Github: https://github.com/Rosemeis/pcangsd

git clone https://github.com/Rosemeis/pcangsd.git

See github for more information regarding 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.

Quick start

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

# Only estimate covariance matrix using 10 threads
python pcangsd.py -beagle data.beagle.gz -o test -threads 10

# Estimate covariance matrix and individual admixture proportions
python pcangsd.py -beagle data.beagle.gz -admix -o test -threads 10

# Estimate covariance matrix and inbreeding coefficients
python pcangsd.py -beagle data.beagle.gz -inbreed 1 -o test -threads 10

# Estimate covariance matrix and perform selection scan
python pcangsd.py -beagle data.beagle.gz -selection -o test -threads 10

Input

The only input PCAngsd needs and accepts are genotype likelihoods in Beagle format. ANGSD can be easily be used to compute genotype likelihoods and output them in the required Beagle format.

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

See ANGSD for more information on how to compute the genotype likelihoods and call SNPs.

Using PCAngsd

All the different options in PCAngsd are listed here. PCAngsd will always compute and output the covariance matrix, where it uses principal components to estimate individual allele frequencies in an iterative procedure. The estimated individual allele frequencies will then be used in any of the other specified options of PCAngsd.

Estimation of individual allele frequencies

-beagle [Beagle filename]

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

-plink [Prefix for binary PLINK files]

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

-plink_error [float]

Incorporate error model for PLINK genotypes.

-minMaf [float]

Minimum minor allele frequency threshold. (Default: 0.05)

-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).

-maf_iter [int]

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

-maf_tole [float]

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

-e [int]

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

-o [prefix]

Set the prefix for all output files created by PCAngsd (Default: "pcangsd").

-indf_save

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

-dosage_save

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

-sites_save

Choose to save the marker IDs after performing filtering using population allele frequencies. Useful for especially selection scans and per-site inbreeding coefficients.

-post_save

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

-threads [int]

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

Call genotypes

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

-geno [float]

Call genotypes with defined threshold (Binary). Numpy format (.npy).

-genoInbreed [float]

Call genotypes with defined threshold also taking inbreeding into account. -inbreed [int] is required, since individual inbreeding coefficients must have been estimated prior to calling genotypes using that information (Binary). Numpy format (.npy).

Admixture

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

-admix

Toggles admixture estimations. Individual ancestry proportions are saved (Binary). Numpy format (.npy).

-admix_alpha [float-list]

Specify alpha (sparseness regularization parameter). Can be specified as a sequence to try several alpha's in a single run (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-list]

Specify seed for random initializations of factor matrices in admixture estimations. Can be specified as a sequence to try several different seeds in a single run.

-admix_K [int]

Not recommended. Manually specify the number of ancestral populations to use in admixture estimations (overrides number chosen from -e). Structure explained by individual allele frequencies may therefore not reflect the manually chosen K. It is recommended to adjust -e instead of -admix_K.

-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_batch [int]

Specify the number of batches to use in NMF method. (Default: 10)

-admix_save

Choose to save the population-specific allele frequencies (Binary). Numpy format (.npy).

Inbreeding

Per-individual inbreeding coefficients incorporating population structure can be computed using three different methods. However, -inbreed 1 is recommended for low depth cases.

-inbreed 1

Simple estimator computed by an EM algorithm. Allows for F-values between -1 and 1. Based on ngsF.

-inbreed 2

A maximum likelihood estimator also computed by an EM algorithm. Only allows for F-values between 0 and 1. Based on [1].

-inbreed 3

Estimator using an estimated kinship matrix. Allows for F-values between -1 and 1. Based on PC-Relate.

-inbreed_iter [int]

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

-inbreed_tole [float]

Tolerance value for the EM algorithms for inbreeding coefficients estimation. (Default: 1e-4)


Per-site inbreeding coefficients incorporating population structure alongside likehood ratio tests for HWE can be computed as follows:

-inbreedSites

Use likelihood ratio tests (.lrt.sites.gz) generated from -inbreedSites to filter out variable sites using a given threshold for HWE test p-value:

-hwe [LRT filename]
-hwe_tole [float]

Tolerance value for HWE test. (Default: 1e-6)

Selection

A genome selection scan can be computed based on posterior expectations of the genotypes (genotype dosages):

-selection

Using an extended model of FastPCA. Performs a genome 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.

Relatedness

Estimate kinship matrix based on method Based on PC-Relate:

-kinship

Automatically estimated if -inbreed 3 has been selected.

Remove related individuals based on kinhsip matrix of previous run:

-relate [Kinship filename]
-relate_tole [float]

Threshold for kinship coefficients for removing individuals (Default: 0.0625).

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 estimating per-site inbreeding sites and testing for HWE has been published in Molecular Ecology Resources: Testing for Hardy‐Weinberg Equilibrium in Structured Populations using Genotype or Low‐Depth NGS Data