Companions

The overarching objective was to allow for broad participation in the biomedical assessments. Thus, all living Project 1 (national survey) respondents were considered eligible for participation if their existing health information indicated an ability to travel to the clinic without excessive risk to the respondent or project staff. Siblings of main sample respondents were not part of the recruitment pool (primarily because of cost), but members of the twin sample were included. Members of the Milwaukee sample of African Americans, newly recruited at MIDUS II, were also part of the recruitment pool. Eligible respondents were first sent a letter explaining what the biological project was about. A brochure sent with the letter sketched the key objectives of the biomedical assessments, outlined what would be included in the clinic visit, and explained how financial matters related to respondents’ time and travel would be handled. Follow-up phone calls were then made to provide additional details and answer any questions the respondent might have. All travel expenses to and from the clinics were covered, and project staff also helped arrange travel itineraries. For aged individuals, or those concerned about traveling alone, an option was provided to travel to the clinic with a companion. Respondents were given $200 in consideration of their two-day visit to the medical clinic. For some, childcare costs were also provided. The study was approved by the Institutional Review Board at each participating center, and informed written consent was obtained for all participants.

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.

To count alignment breakpoints, we mapped all assemblies to the corresponding reference genomes with minimap2 under the option “--paf-no-hit -cxasm20 -r2k -z1000,500”. We used the companion script paftools.js to collect various metrics (command line: “paftools.js asmstat -q50000 -d.1”, where “-q” sets the minimum contig length and “-d” sets the max sequence divergence). To count substitutions and gaps, we applied a different minimap2 setting “-cxasm5 --cs -r2k”. This setting introduces more alignment breakpoints but avoids poorly aligned regions harboring spuriously high number of differences that are likely caused by large-scale variations and skew the counts. We used “paftools.js call” to call variations.