Novel genetic matching methods for handling population stratification in genome-wide association studies

Case-control pairwise matching

The Hungarian Method requires balanced and complete bipartite graphs. Therefore, in case of different cardinalities of the sets of controls O and cases A we extend the smaller set to have max(|O|,|A|) cardinality. We add \(\left ||O|-|A|\right |\) additional elements (“sinks”) to the smaller set. Then we consider the balanced bipartite graph G(V,E,w) with vertices V=O∪A, O∩A=∅, edges E⊆O×A and weights \(w:E\rightarrow \mathbb {R}\). The weights are given by the genetic similarity score of eq. (1), w(i,j)=s _ij . For the edges which are incident upon sinks the weights are set to zero. The Hungarian Method returns a matching M⊂E with max(|O|,|A|) matched pairs. From that we remove matches comprising sinks. Thus, the individuals that are matched to the sinks are also removed, ending up with min(|O|,|A|) cases and controls, respectively. By doing this, we reduce the sample size to an equal number of cases and controls. Subsequently, we perform a stratum validation as described in section ‘Stratum validation strategies’.

Case-control groupwise matching

In order not to reduce asymmetric samples (|O|≠|A|) too extensively, and therefore lose power, one can refit unmatched individuals into the sample. To achieve that, we perform an initial matching as described in section ‘Case-control pairwise matching’. Afterwards we repeat the matching process between matched cases with unmatched controls (those that were removed) and vice versa. The newly matched individuals are added to the case-control pairs of the initial matching. Thus, instead of receiving min(|O|,|A|) case-control pairs, we end up with min(|O|,|A|) small groups of at least one case and one control. The process of re-matching unmatched individuals to matched individuals is repeated until either every individual is matched or until no individual is successfully been re-matched during a single run of the re-matching process. The latter may happen because of the validation procedures which we will introduce in section ‘Stratum validation strategies’. Because of the validation and iterative re-matching it may occur that a group contains more than one case and more than one control.

We also note that one might consider multi-objective matchings similar to [47,48], in order to obtain groups already in the initial steps of the matching. However, the removal and repeated matching of invalid pairs due to the stratum validation completely compensates for the optimality-advantage of more sophisticated methods. We therefore stick to our computationally more efficient method.

Agglomerative hierarchical clustering

The Hungarian Method can be utilized as a basic building block for agglomerative hierarchical clustering. This is called the Hungarian Clustering Algorithm [37]. The algorithm is capable to cluster non-convex data sets. The widely used K–Means- or EM-based clustering approaches, for instance, have difficulties with such data sets. It is also robust to noisy data due to the hierarchical nature that prevents fast propagation of clustering errors. The number of clusters is intrinsically found as part of the process, while the performance of the algorithm is quite competitive to other clustering methods. A description of our implementation of the Hungarian Clustering Algorithm is given in Additional file 1: Appendix B. For comparison to other approaches like the spectral clustering we refer to [37].

Intra-cluster matching

Preceding to the case-control matching, cluster affiliation is obtained by a clustering algorithm (e.g. section ‘Agglomerative hierarchical clustering’). Then the matching is performed in each cluster separately. This enforces case-control pairs/groups to originate from the same stratum. The idea to perform a within-cluster matching after obtaining cluster information is also given in [19].

Vicinity check

is fulfilled the pair is valid, otherwise the pair is removed. The parameter T _cc does guarantee close relatedness in the sense that it removes pairs which are too far apart compared to the vicinities of the individuals. On the other hand it does not enforce a pair belonging to the same stratum. We expect this approach to be robust in the presence of strong asymmetric strata, in particular if strata are overlapping. In other words, the passing condition eq. (3) checks for a well-defined abundance of possible better mates (first term) or if there are no better mates at all (remaining terms).

Squared test: MCAT⁽²⁾

where M=|units| is the number of units. The risk allele frequencies of the unit i are given by \(f_{\textit {co}}^{(i)} = (2n_{11} + n_{12})/n_{\textit {co}}\) for controls and \(f_{\textit {ca}}^{(i)} = (2n_{21} + n_{22})/n_{\textit {ca}}\) for cases. Due to the signed nature of the contribution per unit in eq. (5), the test statistic \(Y_{T^{2}}\phantom {\dot {i}\!}\) cannot be expressed by a \({\chi ^{2}_{n}}\)-distribution, and we will need to employ resampling simulations in order to calculate P-values. In the following, we will call the test statistics \(Y_{T^{2}}\phantom {\dot {i}\!}\) the squared Matching Cochran-Armitage Trend (MCAT⁽²⁾) test.

Linear test: MCAT⁽¹⁾

The resulting test statistic for large M is asymptotically normal distributed with variance M, \({\lim }_{M \rightarrow \infty } Y_{U} \sim \mathcal N(0,M)\). Note, that the convergence is quite fast, therefore it is reasonable to employ this test for more than a dozen pairs/groups while the variance is finite. To be more precise, the test can be transformed to a standard normal distribution where the test statistic scales with \(1/\sqrt {M}\). Likewise, considering large clusters, the test statistic is normally distributed for small M with large numbers of individuals per unit N ⁽ⁱ⁾, \({\lim }_{N^{(i)} \rightarrow \infty, \forall i \in M} Y_{U} \sim \mathcal N(0,M)\). In the following, we will call the test statistics Y _U the linear Matching Cochran-Armitage Trend (MCAT⁽¹⁾) test compared to the squared test MCAT⁽²⁾.

(M)CAT tests with resampling simulation

We obtain P-values for the CAT and MCAT⁽²⁾ tests by utilizing resampling simulations on the basis of within-unit – i.e. within-pair, -group or -cluster – permutation of the case-control trait. The P-values are determined by the fraction of simulations, where the resulting test statistic is equal to or more extreme than the test statistic of the original set. We also allow adjustment for multiple testing by employing the minP approach [49], that has previously been used in the context of pathway association analysis [50]. One considerable strength of the minP approach is that it allows to avoid nested simulations.

Regression models with structure covariates

Regression models with structure covariates provide useful tools to perform stratified analyses without employing resampling simulations. For case-control studies we may employ logistic regression (LR) and for quantitative trait studies linear regression. Population structure covariates are obtained by calculating principal components using MDS or PCA. The P-values are calculated using the likelihood ratio test.

H ₀ simulation

While the naive CAT test with resampling within units over-compensates for stratification effects, both MCAT tests in general perform much better. The groupwise matchings show overall good results in terms of inflation (λ≈1) and consumption of the statistical level (f _p ≈0.05). We find that both validation methods (pre-cluster and vicinity) work well with the groupwise matching. In Table 3 the number of found groups for the vicinity validation is smaller than for the clustering validation. Thus, the vicinity validation is more stringent than the clustering validation. The clustering approach was able to detect all 14 populations in nearly all of the ten samples, sometimes deviating by one. The pairwise matchings reduce the sample to become equal in the number of cases and controls and remove the least fitting individuals. For that reason it is comprehensible that it tends to produce samples that are deflated (λ<1, Table 3).

The MDS approach also yields very good results if we provide the correct number of populations (LRmds14, λ=1.006). However, underestimating the number does only partially correct for population stratification (LRmds07) and also overestimating the number (LRmds28) shows a slow increase of the false-positive rate. We think the latter may happen due to cancellation of redundant components, which leads to a drop-off in the ability to correct for stratification.

Computation time of the resampling simulation (99,999 cycle, 44,900 SNPs) on the MCAT⁽²⁾ test with groupwise matching takes about 450 minutes on a single used core of an Intel®; Xeon®; E5540 CPU with 2.53GHz. Our implementation supports parallelization, which reduces the real time correspondingly.

Power simulation

We create a series of data sets by moving along the simulated chromosome. A particular encountered SNP is assigned to be associated and we assume a relative disease risk of 1.5 under a multiplicative model. Based on this model assumption and on strata-specific baseline allele frequency we simulate cases and controls from each of the 14 population samples of section ‘Simulation of stratified population samples’ with the distribution given in Table 2. In this way, we mimic over/under-sampling of cases from different strata, thereby generating stratified data sets that each contain one true SNP association. By construction, the biased sampling from different strata may blur the true association effect in the joint sample. In total, our procedure yields 11,010 data sets with the case-control status under this model. We use 10,000 SNPs for the matching and keep a gap of at least 1,000 SNPs to the analyzed SNP to guarantee that there is no LD between SNPs that are used for the matching and the associated SNP. We determine the P-values as described in section ‘ H ₀ simulation’. On the results of the power study we perform genomic control: from the P-values we calculate the test statistics by an inverse χ ²-distribution, then correct those test statistics by the inflation factors of Table 3, and finally calculate again the P-values according to a χ ²-distribution.

In Figure 1 we illustrate the power vs. the nominal error level for both MCAT⁽²⁾ and MCAT⁽¹⁾ test and the logistic regression test with covariates. In Table 4 we list the power for three selected nominal error levels (0.01, 0.001, 0.0001) for all employed tests. We observe that the unstructured association testing (CAT - red line) is fully outperformed by all structured association testing methods. The logistic regression model with an over-estimated count of MDS covariates (LRmds28 - dashed gray line) performs considerably weaker than the model with the optimal count (LRmds14 - solid gray line). The naive CAT tests, with resampling simulations performed within units, drop off rapidly for already very large nominal levels (Table 4). Both MCAT tests (colored lines) for the groupwise matching and the cluster are competitive with LRmds14. Reducing the sample size to matched pairs (blue and turquoise lines) reduces power considerably. We conclude, that groupwise matching should definitely be favored over pairwise matchings.

			Power 1− β
Test		Units	Validation	α =0 . 01	α =0 . 001	α =0 . 0001
CAT	AT	–	–	0.665	0.450	0.271
CAT	RSU	Pairs	Cluster	0.881	0.782	0.662
CAT	RSU	Pairs	Vicinity	0.846	0.710	0.558
CAT	RSU	Groups	Cluster	0.809	0.705	0.594
CAT	RSU	Groups	Vicinity	0.833	0.729	0.607
CAT	RSU	Clusters	–	0.808	0.704	0.597
MCAT(2)	RSU	Pairs	Cluster	0.884	0.783	0.658
MCAT(2)	RSU	Pairs	Vicinity	0.823	0.705	0.549
MCAT(2)	RSU	Groups	Cluster	0.907	0.808	0.685
MCAT(2)	RSU	Groups	Vicinity	0.872	0.745	0.591
MCAT(2)	RSU	Clusters	–	0.894	0.782	0.641
MCAT(1)	AT	Pairs	Cluster	0.862	0.733	0.579
MCAT(1)	AT	Pairs	Vicinity	0.813	0.641	0.454
MCAT(1)	AT	Groups	Cluster	0.903	0.797	0.671
MCAT(1)	AT	Groups	Vicinity	0.878	0.751	0.603
MCAT(1)	AT	Clusters	–	0.898	0.786	0.650
LRmds	LRT	7 PCs	–	0.871	0.741	0.588
LRmds	LRT	14 PCs	–	0.908	0.805	0.678
LRmds	LRT	28 PCs	–	0.901	0.798	0.664

11,010 SNPs, 1845 individuals, nominal error level α and power 1−β. 10,000 independent SNPs were used to obtain structure information. The abbreviations in the second column are: AT asymptotic test, RSU resampling simulation within units (99,995 cycles), LRT likelihood ratio test. Column N shows the number of individuals included and in brackets the number of pairs p, groups g and clusters c.