Anopheles gambiae

We compared the accuracy of pixy’s estimates of π and d_XY with several popular existing tools: VCFTOOLS, ANGSD, POPGENOME, and SCIKIT-ALLEL (Danecek et al., 2011 (link); Korneliussen et al., 2014 (link); Miles et al., 2019 ; Pfeifer et al., 2014 (link)). We computed π using PIXY, VCFTOOLS, POPGENOME, SCIKIT-ALLEL, and ANGSD. Note that ANGSD was only applied to the empirical data, since its diversity functions are not designed to work with VCFs. We computed d_XY using PIXY, POPGENOME, SCIKIT-ALLEL, and the ANGSD companion script “calcDxy” (https://github.com/mfumagalli/ngsPopGen/blob/master/scripts/calcDxy.R). For VCFtools, we used the “--window-pi” method to estimate windowed π. For scikit-allel, we used the allel.sequence_diversity and allel.sequence_divergence functions to estimate windowed π and d_XY, respectively. For PopGenome, we used the nuc.diversity.within and nuc.diversity. between functions, following the recommendations in the manual. We stress that the PopGenome manual explicitly warns that computing π and d_XY in the presence of missing data will result in biased estimates (Pfeifer et al., 2014 (link)). We have chosen to include it here because PopGenome is commonly used to estimate π and d_XY in spite of this warning.
We first used our simulated data sets to examine pixy’s accuracy in comparison to these existing methods. To obtain two simulated populations for evaluating d_XY, we split the 100 simulated samples into two random groups. To standardize sample sizes between π and d_XY estimates, we computed π using the first half of the simulated individuals (n = 50), and d_XY by designating the first half of the individuals as drawn from “Population 1” and the second half as drawn from “Population 2” (each with n = 50 individuals). We computed π and d_XY in 10 kb windows in each of the VCFs with variable missing data. pixy was run using default settings, and each pre-existing method was applied using the functions described above (see code supplement).
We then compared the accuracy of each method using the empirical Anopheles gambiae data. To do this, we first applied a basic genotype-level hard filter (DP > =10, GQ >= 40|RGQ >= 40) to the invariant sites VCF produced by GATK. We also removed all variants apart from biallelic SNPs – like other existing methods, pixy does not support the calculation of summary statistics for INDELs or other structural variants. The filtered VCF was used as the input file for all methods apart from ANGSD (see below). We then computed π (PIXY, VCFTOOLS, ANGSD, POPGENOME, SCIKIT-ALLEL) and d_XY (PIXY, POPGENOME, SCIKIT-ALLEL, ANGSD) in 10 kb windows. We computed π separately for the BFS and KES populations. For ANGSD, the BAM files generated from the Anopheles BFS and KES populations were used as input, resulting in estimates of both π (ANGSD’s “pairwise theta”) and d_XY (obtained via a companion script: calcDxy – by Joshua Penalba, https://github.com/mfumagalli/ngsPopGen/blob/master/scripts/calcDxy.R). In the case of π, we explicitly divided the raw estimates of pairwise theta by the number of sites (nSites) provided by ANGSD, and not the window size (10,000).
Full details and scripts for all of the above procedures are available at https://github.com/ksamuk/pixy_analysis.