PCAmatchR: a flexible R package for optimal case–control matching using weighted principal components

1 Introduction

Genome-wide association studies (GWAS) have discovered thousands of germline susceptibility variants associated with disease risk. The selection of cases and well-matched controls is vital for obtaining unbiased association results (Hinds et al., 2004; Luca et al., 2008). An increasingly observed trend in GWAS is the practice of maximizing study resources by only genotyping cases and then borrowing controls from large sets of pre-genotyped, disease-free individuals (Hinds et al., 2004; Luca et al., 2008). When borrowing controls from previous studies, it is essential that cases be genotyped on the same arrays and undergo the same quality control and filtering steps, to lessen the possibilities of technical artifacts in genotype calling. In some scenarios, pooling of control genotyping data with case genotyping data can be an effective approach for performing genetic association studies as long as case–control imbalances are appropriately accounted for in the statistical analysis (Hayeck et al., 2015; Ma et al., 2013; Zhou et al., 2018). However, there are instances in which close genetic matching of cases to controls may be necessary (e.g. selecting subjects for costly molecular analyses).

Combining data from different sources could lead to confounding due to population stratification, i.e. genetic differences in cases and controls based on ancestry and not disease status (Epstein et al., 2012; Lacour et al., 2015), as allele frequencies and linkage disequilibrium patterns vary substantially by ancestry (Machiela et al., 2015). Many methods exist for removing population stratification bias from pooled genotype data (Byun et al., 2017; Epstein et al., 2012; Hinds et al., 2004; Lacour et al., 2015; Luca et al., 2008). A popular strategy is the utilization of principal component analysis (PCA) (Byun et al., 2017; Price et al., 2006). Briefly, PCA involves transforming large sets of predetermined ancestry-related genetic variants to a set of linearly uncorrelated principal components (PCs), the first of which explain the highest percentage of variation in the full genetic dataset (Price et al., 2006). Thus, instead of matching cases and controls on all available genotyped variants, close matching can be efficiently achieved using only the first few PCs.

Here, we describe the use of PCAmatchR, an open-source R package which allows for matching of controls to cases using PCA results (Brown et al., 2020). While case–control matching using derived PCs is not a novel concept (Machiela et al., 2018), our methodology involves weighing each PC by the respective percent variation explained. Using an improved weighted Mahalanobis distance metric, our approach selects optimized matches which ensure PCs that explain the highest degrees of variation within the genomic data are more heavily weighted. This approach produces a set of matched controls that are more genetically similar to cases and results in reduced genomic inflation in genetic association studies.

2 Materials and methods

We developed PCAmatchR to optimally match population-based controls to cases. PCAmatchR converts PCs into a Mahalanobis distance metric for selecting well-matched controls. Mahalanobis distance is a commonly used multivariate distance metric that creates a standardized distance between two samples, where smaller Mahalanobis distance indicates greater sample similarity (Hu et al., 2013). This distance metric has been highly utilized for matching within propensity score analyses (Rubin, 1979; Zhao, 2004). The distance $\left(D_{ij}\right)$ between a vector of input PCs $\left(X\right)$ for subjects $i$ and $j$ is given by:

where $\Sigma$ is the variance covariance matrix of $X$ in the pooled case and control sample population (Stuart, 2010). In the case of weighted matching, Equation (1) is extended to allow for weights in the distance metric. The weighted Mahalanobis distance metric is defined as:

where W is given by:

and each $w_n$ represents the percent variance in data explained for each of the $n$ PCs (Hu et al., 2013; Krusińska, 1987). The percentage of variance explained by each PC is calculated by taking the PC’s corresponding eigenvalue and dividing by the sum of all eigenvalues generated in the PCA (Jolliffe et al., 2016); although other user-defined weights can also be applied. Unlike the standard Mahalanobis distance metric which underperforms in terms of matching when the number of covariates is large (greater than 8), as the metric treats each element of $X$ with equal importance (Rubin, 1979; Stuart, 2010), the weighted Mahalanobis distance metric may produce more genetically similar matches, as PCs are weighted based on their relative importance (Krusińska, 1987). PCAmatchR utilizes the optmatch package, which performs bipartite matching using the RELAX-IV minimum cost flow solver (Hansen et al., 2006). Briefly, this optimal matching algorithm jointly minimizes the total Mahalanobis distance among all possible matches, which removes the effects of possible ties (i.e. one control being the best match for multiple cases) from the final matched set (Hansen, 2007; Rosenbaum, 1989).

PCAmatchR takes as input PCs and eigenvalues and directly outputs optimal case and control matches. The first step is to choose a subset of genotyped single nucleotide polymorphisms (SNPs) for use as input when calculating the PCs (Price et al., 2008; Tian et al., 2008; Yu et al., 2008). PCA can then be conducted in the normal way using any standard tool, e.g. PLINK -- pca, on a combined dataset including both the case and control data. The user should retain the output matrix of derived PCs and the corresponding eigenvalues to perform the matching procedure. PCAmatchR has the flexibility to both match any number of controls to each case (e.g. 1:n matching) and further create exact matches based on user-defined covariates (e.g. sex), assuming sufficient control population size.

3 Implementation

To demonstrate the utility of PCAmatchR, we performed a hypothetical example of case–control matching using data from the 1000 Genomes Project (Auton et al., 2015; Sudmant et al., 2015) Phase 3 data release, which contains genotype data from 2504 individuals from 26 distinct populations (available at https://www.cog-genomics.org/plink/2.0/resources). Using a set of ancestry informative SNPs (Yu et al., 2008), we performed a PCA on all available genotyped participants using PLINK. The first 20 PCs explained 16% of the variability of the genomic data with the first PC explaining 9% (Supplementary Fig. S1a).

For a sample analysis, we selected all individuals from the CEU [Utah Residents (CEPH) with Northern and Western European Ancestry] population as cases (N = 99), whereas all remaining samples were used as our control population (N = 2405) (Supplementary Figs S1b and S2a). Cases were 1:1 matched to controls based on the weighted Mahalanobis distance metric using the first 20 PCs (Supplementary Fig. S2b and Supplementary Code). Using the plotting functionality within PCAmatchR, the connections between cases and matched controls are depicted in Figure 1,further providing a visual representation of the multivariate distance between matches. To show that the weighted matching procedure correctly identified individuals with higher genetic similarity, we calculated a relationship matrix, which estimates the SNP-based genetic relationship between samples, using PLINK -- make-rel (Yang et al., 2011). The genetic relationships between matches were extracted and averaged across all matches to give an overall estimate of genetic similarity. For comparison, we also performed the matching procedure using the unweighted Mahalanobis distance metric. Overall, controls selected using the weighted Mahalanobis distance metric had higher average genomic similarity to the CEU cases than the controls selected using the unweighted metric (0.1216 versus 0.1172, Supplementary Fig. S3 and Supplementary Table S1). Finally, we performed GWAS analyses, based on the ancestry SNP panel variants, and observed the weighted Mahalanobis distance metric produced a smaller lambda value than the matches selected using the unweighted metric (Supplementary Fig. S4 and Supplementary Table S1), indicating the weighted matching procedure within PCAmatchR effectively removed more potential population stratification bias in the pooled sample, compared to the unweighted distance-based matches.

For completeness, we replicated these analyses for each of the 26 available 1000 Genomes Project populations, selecting each individual population as a separate case sample. Across the populations, matches derived using the proposed weighted Mahalanobis distance metric outperformed the unweighted metric matches in terms of genomic similarity within 21 (81%) of the populations and in terms of lambda value within 23 (88%) of the populations (Supplementary Fig. S5 and Supplementary Table S1). The importance of matching was further demonstrated using comparisons to unmatched samples with randomly selected controls (size equal to the number of cases). Weighted Mahalanobis distance matching greatly reduced the effects of population stratification compared to both baseline levels as well as those achieved by current unweighted matching techniques (Supplementary Table S1).

Sensitivity analyses were performed to assess the effects of the number of PCs included when matching. Weighted matching was less sensitive to the number of included PCs compared to unweighted matching (Supplementary Table S2). When the number of PCs increased (>20), unweighted matching poorly removed the effects of population stratification from analyses. In practice, the number of PCs included for matching should be based on the percent variance explained by the eigenvalues.

4 Discussion

PCAmatchR is an efficient and easy-to-use open access tool to match cases and controls based on user defined PCs, which can aid in removing the effects of population stratification for association analyses as well as facilitate the selection of biological samples for subsequent costly assays (e.g. methylation arrays and whole-genome sequencing). When applied to available data from the 1000 Genomes Project, PCAmatchR demonstrated the ability of weighted PCs to increase genomic similarity between cases and controls and reduce genomic inflation.

Supplementary Material