Association mapping has permitted the discovery of major QTL in many species. It can be applied to existing populations and, as a consequence, it is generally necessary to take into account structure and relatedness among individuals in the statistical model to control false positives. We analytically studied power in association studies by computing noncentrality parameter of the tests and its relationship with parameters characterizing diversity (genetic differentiation between groups and allele frequencies) and kinship between individuals. Investigation of three different maize diversity panels genotyped with the 50k SNPs array highlighted contrasted average power among panels and revealed gaps of power of classical mixed models in regions with high linkage disequilibrium (LD). These gaps could be related to the fact that markers are used for both testing association and estimating relatedness. We thus considered two alternative approaches to estimating the kinship matrix to recover power in regions of high LD. In the first one, we estimated the kinship with all the markers that are not located on the same chromosome than the tested SNP. In the second one, correlation between markers was taken into account to weight the contribution of each marker to the kinship. Simulations revealed that these two approaches were efficient to control false positives and were more powerful than classical models.
QUANTITATIVE traits are determined by the polymorphism of many genes or genomic regions with small effects, i.e., quantitative trait loci (QTL). Understanding the genetic architecture of such traits, which supposes the identification of these causal loci, is now facilitated by a dramatic increase in the number of molecular markers available. This makes it possible to conduct genome-wide association studies (GWAS), in which phenotypes and genotypes of individuals in highly diverse panels are used to detect QTL (Lynch and Walsh 1998). Such panels have accumulated numerous historical recombinations, leading to a low extent of linkage disequilibrium (LD). Compared to linkage mapping, more markers are therefore needed to capture causal signals but with a much higher mapping resolution (Rafalski and Morgante 2004). Major genes were identified by this approach in human, animal, and plant genetics (Ozaki et al. 2002; Beló et al. 2007; Jones et al. 2008). However, contrary to linkage mapping populations, LD in association mapping panels is not only due to genetic linkage, but can also be caused by population structure, relatedness, drift, and selection (Jannink and Walsh 2002; Flint-Garcia et al. 2003). The contribution of these factors relative to linkage can be evaluated statistically (Mangin et al. 2012) and proved for instance to be substantial in grapevine and maize (Mangin et al. 2012; Bouchet et al. 2013). This component of LD due to population structure and relatedness can generate false positives and thus must be taken into account in association mapping models to control false positives (Ewens and Spielman 1995; Thornsberry et al. 2001). Once these effects are correctly modeled, only marker-trait associations due to linkage should be detected.
Population structure can be estimated with softwares such as STRUCTURE (Pritchard et al. 2000) and ADMIXTURE (Alexander et al. 2009) or by principal component analysis on the genotypic data (Price et al. 2006). These methods permit the estimation of a structure matrix (Q) attributing the admixture coefficient of each individual in each group. Relatedness (K matrix) can be estimated in different ways including identity by state (IBS), or estimators of identity by descent (IBD) considering marker allelic frequencies (Vanraden 2008; Astle and Balding 2009). Yu et al. (2006) proposed a mixed-model approach (Q + K) to detect QTL in the context of association mapping. This model has the advantage of controlling false-positive rate by including a fixed structure effect (through Q) and/or a random polygenic effect (through K). It was used in many association mapping studies and permitted the detection of QTL in humans, animals, and plants (K. Zhao et al. 2007; Huang et al. 2010; Kang et al. 2010; Price et al. 2010; Zhang et al. 2010; Bouchet et al. 2013; Romay et al. 2013). However, one of the main drawbacks of these structure and relatedness corrections is that it also reduces the number of detectable true positives, particularly if the trait is correlated to the population structure (Larsson et al. 2013). Also, including the tested SNP in the computation of K is expected to decrease power at this SNP (Listgarten et al. 2012). To increase the power of GWAS, some authors therefore proposed using only a subset of SNPs as covariates or to estimate genetic similarity (Listgarten et al. 2012; Bernardo 2013). Speed et al. (2012) proposed weighting the contribution of the SNPs in the kinship estimation to increase the accuracy of heritability estimates.
It is particularly important to evaluate the power of panels and statistical approaches to discover QTL. Power may be analytically investigated using the noncentrality parameter of the test statistics. This strategy has first been applied in linkage mapping, where several authors showed how power is influenced by the size of the population, heritability, the effect captured by the marker, and the allelic frequencies (Soller et al. 1976; Knapp and Bridges 1990; Rebai and Goffinet 1993; Charcosset and Gallais 1996). Such an analytical approach has also been applied in association studies in human and animal genetics (Sham et al. 2000; Purcell et al. 2003; Wang 2008; Teyssèdre et al. 2012). Alternatively, the estimation of power has also been addressed through simulation studies (see, for instance, Yu et al. 2006; H. Zhao et al. 2007; Stich and Melchinger 2009; Erbe et al. 2010; MacLeod et al. 2010; Bradbury et al. 2011). We can retain from these studies that power of association mapping diminishes with structure and relatedness in addition to the parameters identified in linkage analysis and that the way of estimating K has an effect on power (Stich et al. 2008). To our knowledge no study was conducted to compare analytically the power along the genome in different association mapping designs.
In this study we analytically derived the power at each marker for the classical mixed model involving relatedness between individuals (Yu et al. 2006). This analytical expression of power makes it possible to study the effect of different parameters on local power along the genome. We first used it to compare three diversity panels with different diversity patterns. We highlighted a loss of power due to the use of the genotypic information both to test marker effect and to estimate K, and this was particularly strong in regions of high LD. We therefore evaluated two alternative estimation strategies of the kinship matrix to increase power in GWAS. In the first one, we used an estimated K matrix specific to each chromosome: only the markers that are physically unlinked to the tested SNP are used to estimate K. In the second one, we weighted the contribution of each marker in the estimation of K by taking into account intrachromosomic LD. We compared in simulations based on true genotypes of maize inbreds the efficiency of the different strategies to detect QTL and to control false positives.
Materials and Methods
Statistical models for association mapping and power evaluation
Mixed models are now routinely used to control type I error in GWAS (Yu et al. 2006). Relatedness among individuals is taken into account by considering that the random polygenic effects are not independent, with a covariance matrix determined by kinship (K, with as many rows and columns as individuals, N). As K includes information on both population structure and relatedness, it is in general not useful to consider admixture information as fixed effects covariates (Astle and Balding 2009). We therefore considered the following statistical model (denoted by MK),where Y is the vector of N phenotypes, μ is the intercept, 1 is a vector of N 1, is the vector of N genotypes at the tested locus (0 and 1 corresponding to homozygotes and 0.5 to heterozygotes), is the additive effect of locus l to be estimated, U ∼ N(0, ) is the vector of random polygenic effects, being the residual polygenic variance, E ∼ N(0, ) is the vector of remaining residual effects with variance , I is an identity matrix of size equal to the number of individuals (N), and U and E are independent.
Locus effects in this mixed model can be tested using Wald statistics (Wald 1943). In the general case, a given linear combination of fixed effects (H0 hypothesis) can be tested against (the alternative hypothesis H1) usingwhere is a vector of fixed effect estimates, L is a linear combination, and and are the REML estimates of and .
In GWAS we test the particular linear combination: against , with if the only fixed effects are the intercept and the marker additive effect. Note that the approach could be extended to more complex effects such as dominance by adding extra term(s) in fixed effects. When the variances are known, W follows a χ2 distribution: , where and λ is the noncentrality parameter (NCP). The noncentrality parameter is equal toUnder H0, , whereas under H1, λ is positive. Power can thus be determined as the probability P(), λ being the NCP and the value of the central (1 − α) quantile, where α corresponds to the chosen type I error level. The power of the test increases as the NCP increases. λ depends on the QTL effect (the magnitude of departure from H0), the marker genotypes, and the variance and covariance components. Hence in addition to the number of individuals, power can be influenced by the marker genotypes, the marker effect (), the heritability (through and ), and the relatedness between individuals (K).
Analytical evaluation of the impact of panel characteristics on power
When genotypic data are available in a given association mapping panel, it is possible to analytically evaluate power at each marker thanks to the above formula. Consider a panel in which N individuals were genotyped at M markers (SNPs). The potential power at a given marker can be investigated by setting a QTL effect βl, a background genetic variance , and a residual variance to reach a given heritability h2. Power at a given marker can then be related to parameters characterizing the marker in the panel of interest. It is expected first to depend on allele frequencies, which can be characterized by the minor allele frequency (MAF). Also, according to the analytical expression of the NCP, power at a marker in MK can be influenced by its correlation with the kinship that reflects both the structure of the panel and the relationships between individuals. It is thus interesting to relate power at a given marker to its Nei index of differentiation (Fst) among genetic groups (Nei 1973) and to its correlation with the kinship matrix. Let us denote by K_Ml the kinship matrix evaluated from the considered marker l only. To define how power at a given marker is affected by its correlation to K, one can calculate the correlation between K_Ml and K at each marker. This correlation between local and global kinship is further referred to as CorK. These statistics (Fst, MAF, CorK, and analytical power) can be calculated for each marker in any association mapping panel.
In this article, we applied this strategy to three maize panels (see below). We represented the relationship between MAF, Fst, CorK, and local power with the two following approaches. In the first one, analytical power was represented as level plots considering MAF and Fst as x- and y-axes, with the R function level.plot. The same procedure was applied to MAF and CorK. In the second approach, cubic smoothing splines were adjusted along the genome to the Fst, CorK, and power for the markers with a MAF above 0.4, using the R function smooth.spline (Hastie and Tibshirani 1990).
In practice the kinship matrix K is unknown and must be estimated. One classically used estimator was proposed by Astle and Balding (2009) and is defined aswhere Gi,l and Gj,l are the genotypes of individuals i and j at marker l (Gi,l = 0 or 1 for homozygotes, 0.5 for heterozygotes), 1, and is the variance of Gi,l, respectively. One problem that might arise from this formula and other classical estimators as the identity by state, or the formula of Vanraden (2008), is that LD between SNPs is not taken into account. As a result more weight is given in the kinship estimation to the regions of the genome that carry several markers in strong LD and power may be lower in these regions.
We therefore considered two alternative approaches to limit this effect. In the first one, the kinship matrix (K_Chr) was estimated with all the markers other than those located on the same chromosome as the marker being tested. If the markers located on the other chromosomes are sufficient to reliably estimate relatedness, this method is expected to reasonably control the risk of detecting false positives and avoids considering in the kinship matrix markers linked with the tested markerwhere c is the considered chromosome and is the number of markers not located on chromosome c.
In the second approach we used all the markers as estimators of relatedness but we weighted the contribution of each marker. The kinship estimator K_Freqi,j can be understood as follows: each marker l yields an estimator of the true kinship coefficient kij between individuals i and j, which are then averaged over all markers to obtain . This average would be optimal if all estimators had the same variance and were independent. In practice none of these conditions is satisfied: the error variance of each estimator depends on the MAF of the marker, and LD between markers generates correlations between markers. As a consequence, estimators with poor precision (high error variance) will have the same weight as estimators with high precision. Moreover, m highly correlated estimators will accumulate a weight of m/L without providing m independent information; i.e., too much weight is attributed to highly correlated estimators. Alternatively, one may look for the weighted combination , which is the best linear combination of coefficient , to estimate without bias. Define and as the mean and variance of estimator over all couples of individuals (i,j) having the same kinship . Note Δ the covariance matrix between estimators , i.e., , the vector of weights, and the vector of marker estimators. Then satisfies
In this formulation the optimal weights may be negative; we added extra constraints to ensure the positivity of the weights, leading to the following optimization program:
, for all l. (1)
In practice, obtaining the optimal weights requires (i) the knowledge of matrix Δ and (ii) solving the optimization problem (1). The exact expression of matrix Δ is unknown, but one can estimate this matrix from the panel data using the classical moment estimator:The resulting estimated matrix is then plugged into the optimization program (1). Then to solve the optimization program, one should note that (1) is a quadratic problem with linear constraints and therefore can be solved using classical optimization techniques (in this article we used the R package solve.QP, which implements the dual method of Goldfarb and Idnani 1983).
The main limitation of this strategy lies in step (i): when estimating the covariance, one actually replaces the expectation over all couples having the same kinship by averaging over all couples in the panel—assuming implicitly that they all have the same kinship. Even if the kinship differs between couples, this weighting increases the contribution of markers with a high diversity (leading to a high precision) and not highly correlated with other markers. It therefore corrects the two drawbacks of the naive averaged estimator mentioned earlier.
Let us denote the statistical model for association mapping described above by MK_Freq, MK_Chr and MK_LD with K estimated as K_Freq, K_Chr, and K_LD, respectively.
Simulation-based evaluation of the impact of the estimation of K on false-positive control and power
The closed-form expression of the noncentrality parameter already revealed that kinship affects power. Comparing the impact of different kinship estimators on power implies evaluation of their ability to guarantee the expected nominal control of false positives under different hypotheses on trait genetic determinism. To this end, we simulated traits influenced by L biallelic QTL (SNPs). In a first step, QTL were sampled randomly among the SNPs located on all the chromosomes except one. The chromosome without QTL (further referred to as H0 chromosome) was used to estimate the false-positive rate. All the H0 markers (the markers on the H0 chromosome) were tested with the above-mentioned statistical models for each run of simulation. The efficiency of the different estimations of K to control false positives was evaluated by comparing expected and observed quantiles of H0 P-values and histograms of H0 P-values. In a second step we applied the same procedure, but sampling the QTL among the M SNPs (on all chromosomes). A QTL was declared detected when the P-value of the corresponding SNP in the genetic model was below the significance threshold. Power of a given model was computed as the number of QTL that were detected. We also applied a less restrictive definition of QTL detection, considering that a QTL could be detected by SNPs located near it. To do so, another analysis was conducted in which markers within a given genetic distance of a QTL were considered H1 markers and the others H0 markers. The realized false discovery rate (FDR) found was defined as the proportion of H0 markers among the markers declared significant. Power of QTL detection was estimated by considering that a QTL was detected when at least one of the corresponding H1 markers had a significant P-value. This general method will be exemplified with parameters specific to three maize panels, described below.
Genetic material and genotyping data
The above-mentioned power analyses (analytical evaluation of power and simulation based evaluation of alternative methods) were applied to three diversity panels of maize. (See File S1, File S2, and File S3.) The first panel (called C-K) was described in Camus-Kulandaivelu et al. (2006). It is composed of 375 inbred lines covering American and European diversity. It includes Tropical, Dent, and Flint lines. The second and third panels are the Dent and Flint panels of the “Cornfed” project (CF-Dent and CF-Flint), described in Rincent et al.(2012). They include lines of the C-K panel and lines derived from recent breeding schemes. Both are composed of 300 lines. These panels were genotyped with the 50k SNPs array described in Ganal et al.(2011), as presented in Bouchet et al. (2013) and Rincent et al. (2012). Individuals with marker missing rate higher than 0.1 and/or heterozygosity rate higher than 0.05 were eliminated. Markers with missing rate higher than 0.2 and/or average heterozygosity rate higher than 0.15 were eliminated. Markers missing the rate and/or with average heterozygosity >0.2 and 0.15, respectively, were eliminated. In each panel, few individuals were highly related. One individual was removed for pairs identical for >98% of the loci. In total 315, 277, and 267 individuals and 44487, 45434, and 44255 markers passed the genotyping filter criteria for the C-K, CF-Dent, and CF-Flint designs, respectively. Missing genotypes (<2% in all panels) were imputed with the software BEAGLE (Browning and Browning 2009). Panels were all adjusted to 267 individuals to compare power for a same population size. Individuals removed were chosen at random. To avoid the ascertainment bias noted by Ganal et al. (2011), we used only the markers that were developed by comparing the sequences of nested association mapping founder lines (PANZEA SNPs; Gore et al. 2009) in the estimation of admixture and relationship coefficients (29996, 30119, and 29132 markers passed the filter criteria for the C-K, CF-Dent, and CF-Flint lines respectively).
Admixture in the CF-Dent and CF-Flint panels was investigated using the SNP data with the software ADMIXTURE (Alexander et al. 2009), with a number of groups equal to four, determined according to the cross-validation procedure presented in ADMIXTURE. For the C-K panel we used the admixture in five groups estimated by Camus-Kulandaivelu et al. (2006) using 55 SSRs chosen for their broad genome coverage and reproducibility. We estimated the differentiation index among genetic groups (Fst, Nei 1973) at each marker using the R package r-hierfstat (Goudet 2005).
Finally, the relationship between LD and power along the genome can be empirically investigated using two different measures of LD. Raw LD can be estimated as the squared correlation between allelic doses at two loci (r2). Linkage related LD (denoted by r2K) can be estimated using the algorithm proposed by Mangin et al. (2012), which corrects r2 by K_Freq. LD within these panels (r2), possibly corrected by K_Freq (r2K), was estimated within a sample of 4000 markers regularly spaced on the physical map.
For analytical investigation of power in the three maize panels, the total additive genetic variance was set to 1000, was set to 17.9, which corresponds to a QTL explaining 8% of the total genetic variance if it had a MAF of 0.5, and was chosen to obtain a heritability of 0.8. Under these hypotheses, analytical power was investigated for an α type I risk equal to 1.25 × 10−6, which led to a risk of 0.05 with a Bonferroni correction on 40,000 tests. We also considered less-stringent thresholds corresponding to Bonferroni corrections on 4000 and 400 tests, although the number of tests was always the same. Power under these hypotheses was calculated in R 3.0.0 (R Development Core Team 2006) for each marker.
To estimate kinship with the different formulas presented above, we considered all the individuals to be inbred and we estimated as . For comparing the different methods for kinship estimation, we simulated traits influenced by 50 or 100 biallelic QTL (QTL effects follow a geometric series as in Lande and Thompson 1990, with parameter a set to 0.96 and 0.98 when 50 or 100 QTL were simulated, respectively). The sign of allelic effect at a given locus was assigned randomly. Genotypic values of the individuals were calculated as the sum of the allelic effects at these QTL. Phenotypes were obtained by adding a residual noise following a normal distribution with mean 0 and variance equal to , where the heritability is set to 0.8. We performed 100 runs of simulations for each scenario using the R 3.0.0 software (R Development Core Team 2013). Each chromosome was used 10 times as the H0 chromosome. For all simulations, the statistical tests were performed with EMMAX (Kang et al. 2010) to reduce computational time and then with ASREML-R (Gilmour et al. 2006) on the markers that had a P-value <0.001 with EMMAX. For P-values >0.001, P-values obtained with EMMAX and ASREML-R were very close and highly correlated. As investigations of the two criteria for QTL detection (causal factor only or window around it) led to very comparable results with respect to the main focus of our study, results considering a window around causal factor are presented as supporting information, Table S1.
Diversity and linkage disequilibrium in maize panels
Diversity and LD were investigated within the different panels to provide elements on their ability to detect QTL (i.e., their power) along the genome. On average, the MAF was lower in the CF-Flint than in the other panels. Differentiation among genetic groups (Fst) was higher for CF-Dent (0.15) than for C-K (0.11) and CF-Flint (0.08) (Table 1). The raw LD (r2) and its correction by kinship (r2K) were variable between and within panels (Figure 1). LD was on average higher in the dent panel. Within each panel, it was higher for centromeric than for telomeric regions. High r2 values were observed between physically linked markers but also unlinked markers. This last situation occurred mainly between centromeric regions (Figure 1A, chromosomes 5, 7, and 8 and Figure 1B, chromosome 7). Interchromosomic LD was reduced to a large extent when considering r2K rather than r2. Taking into account covariance between individuals (r2K) also reduced intrachromosomic LD, in particular between distant blocks with high LD (Figure 1B chromosome 10). Considering r2K instead of r2 globally had the strongest impact in the CF-Dent panel.
Relationship between MAF, Fst, CorK, and power
Above-described parametrization of QTL effects was used to investigate the influence of MAF, Fst, and the correlation between local and global covariance matrices (estimated as CorK_Freq) on power in the three maize panels. Level plots (Figure 2) showed that the MAF, the Fst, and CorK_Freq had important effects on power, with very similar graphs in all the panels. The highest power was achieved when MAF was high and Fst or CorK_Freq was low. When the MAF was <0.1, power was close to 0 even if the marker had a low Fst or low CorK_Freq. Some regions of the level plots were not covered by the available markers (regions in white on Figure 2); in particular there was no marker with a CorK_Freq <0.03. Note that the graphs obtained using K_Chr (or the IBS) were similar to those obtained with K_Freq and led to the same general conclusions (results not shown).
The parameters related to power (MAF, Fst, CorK_Freq) varied between panels (Table 1, see above). As a consequence of the above-described relationships, the mean analytical power of statistical model MK_Freq varied between the three panels (Table 1) and was higher in the C-K panel (11.3%) than in the CF-Dent and CF-Flint panels (<9.0%).
Variation of analytical power and CorK along chromosomes
Power scans (analytical power at each marker plotted against its physical position) of model MK_Freq revealed an extreme variability along the genome in the three panels (Figure 3). In all panels, power at a given location ranged from zero to a maximal value, which depended on the position according to a V-shaped curve (Figure 3 and Figure 4). This maximal value was the lowest near centromeres and the highest near telomeres. This global trend was particularly strong in the CF-Dent panel and less pronounced in the C-K panel, for which the maximum power was stable for larger segments. The V-shaped curve also had different local trends for the different chromosomes for a given panel. For instance in the CF-Flint panel, depletion in power in centromeric region was longer for chromosome 7 than for chromosome 6 (Figure 3C).
Power of model MK_Freq was in accordance with trends of CorK_Freq along the genome. Correlation between the covariance matrix at the marker and the global covariance matrix (K_Freq and K_Chr) was significantly lower for K_Chr than for K_Freq, and particularly in the pericentromeric regions (Figure 4). We observed that peaks of Fst corresponded generally to peaks of both correlations (CorK_Freq and CorK_Chr) (Figure 4B, chromosome 7, and Figure 4, A and C, chromosome 8). Conversely, pericentromeric regions with low Fst corresponded to a peak of CorK_Freq and a drop of CorK_Chr (Figure 4B, chromosomes 8 and 10, and Figure 4C chromosome 7). CorK_Freq, CorK_Chr, and the difference between these two parameters were higher in the CF-Dent panel than in the two others.
Simulation-based assessment of kinship estimation on false-positive control and power
Simulating different genetic models using the genotypes of the three panels allowed the comparison of the efficiency of the three statistical models to control false positives and to detect QTL. The efficiency to control false positives depended on the genetic model (number of QTL), the panel, and the estimation procedure for K (Table 2). The distribution of the P-values under H0 revealed that MK_Freq was conservative (Figure 5A) whereas the alternative models MK_Chr and MK_LD gave distributions closer to the expected one (Figure 5, B and C). The observed P-value quantiles were closer to the expected P-value quantiles with MK_Chr and MK_LD than with MK_Freq (Table 2). MK_Freq resulted in fewer than expected small P-values under H0; for example, in the CF-Dent panel we observed only half of the P-values that were expected to be <0.001. Observed P-value quantiles with MK_Chr and MK_LD were very close to the expected P-value quantiles, although also most of the time below it.
The second step of the simulations revealed the ability of the different statistical models to detect QTL in the different panels. With the usual Bonferroni correction, only few QTL were detected (Table 3). In each scenario MK_Chr and MK_LD were more powerful than MK_Freq. For example, they respectively permitted the detection of 2.1, 1.3, and 1.2 QTL (SNP considered as QTL) on average in the CF-Dent panel when 50 QTL were segregating. The difference of power (proportion of SNP considered as QTL detected) between the different models was more important for less stringent significance threshold. The difference of power between MK_Chr and MK_Freq was the highest in the CF-Dent panel. More QTL were found in the scenario with 50 QTL than in the scenario with 100 QTL. This was expected, QTL having a lower effect on the trait in the 100 than in the 50-QTL scenario.
Discussion and Conclusions
Analytical investigation of potential power along the genome with usual model (MK_Freq)
Power is a key parameter in association mapping, because it indicates how likely the discovery of a QTL is. We presented a general method based on noncentrality parameter to derive analytically theoretical power at each marker locus in a given panel of individuals. It was applied to three different association mapping panels. While being adjusted to the same population size, these different panels had different average power. They also displayed different local patterns of power along the genome.
Power could be related to three parameters characterizing each marker: its MAF, its differentiation index among genetic groups (Fst), and the correlation between its individual kinship matrix with that estimated with all the markers (CorK_Freq when K_Freq is considered). Power at a marker with a low MAF is limited, even if this marker is orthogonal to structure and kinship (Figure 2 and Figure 3). This effect was highlighted already for linkage studies (Soller et al. 1976; Charcosset and Gallais 1996) and GWAS (Lonsdale et al. 2013) and can be explained by the fact that when one of the two alleles is rare, the marker cannot contribute much to the genetic variation. The correlation between kinship at single markers and the global kinship had a strong negative effect on power (Figure 2). The Fst among the admixture group also had an important effect on local power (Figure 2 and Figure 3). This confirmed that admixture is reflected by the kinship matrix, because differentiated regions had a low power although we used a model with relatedness but no admixture (MK_Freq). The level plots showing analytical power at different MAF and CorK_Freq were very similar in the three panels (Figure 2, A2, B2, C2), but those showing power at different MAF and Fst differed (Figure 2, A1, B1 and C1). This suggests that group differentiation has different relative contribution to local kinship variation in the different panels. At a given pair of MAF and Fst values, power was lower in the CF-Dent and CF-Flint panels than in the C-K panel, whereas five groups were used in this panel instead of four in the two others. The C-K panel is composed of highly diverse groups (Tropical, Dent, and Flint lines) and so the admixture matrix captured ancestral population structure but only a small part of kinship. On the opposite, the CF-Dent and CF-Flint panels are composed of less heterogeneous material and so the admixture matrix captured more relatedness. Finally, the shape of the level plots (Figure 2) also suggested that the effect of the different parameters affecting power was not additive. For example Fst and CorK_Freq had a stronger effect on power for markers with higher MAF, and MAF had a stronger effect on power for less differentiated markers. These results show that controlling false positives using the K_Freq model also implies reducing power at differentiated markers (Larsson et al. 2013). It is interesting to note that no marker had a CorK_Freq <0.03 (Figure 2). To investigate the maximum power that could be reached theoretically, we generated for each panel a vector of zeros and ones simulating a marker genotype and applied a simple exchange algorithm until analytical power reached a maximum. These virtual markers (one for each panel) had analytical power much higher (>0.8) than the maximal analytical power of the existing SNPs (<0.44 in each panel). They had a MAF of 0.5 and a CorK_Freq value <0.017. This difference illustrates that the maximum power is strongly constrained by the evolution process that led to the panels.
Both Fst and CorK_Freq appeared highly variable along the genome in each panel. High differentiation (Fst) was observed in particular in pericentromeric regions (Figure 3, A and C, chromosome 8, and Figure 3B, chromosome 7). Pericentromeric regions are known to be more structured than telomeric regions (Carneiro et al. 2009; Franchini et al. 2010) because of lower recombination rates. CorK_Freq was also higher in regions of high LD (mostly pericentromeric regions; see Figure 1 and Figure 4). Beyond the effect of group differentiation, markers in regions of high LD are indeed correlated to many other SNPs that all contribute to the estimation of K_Freq. These LD and Fst features led to the observed V-shape analytical power curve along the chromosome, particularly in the CF-Dent panel in which LD was more extended (Figure 1 and Figure 3). This is in good agreement with published manhattan plots of GWAS results, which showed a reduced number of low P-values in the centromeric regions (Bouchet et al. 2013; Larsson et al. 2013). In our three panels, we observed that this problem also arose with other classical estimators of relatedness (results not shown) such as the IBS estimator or the first estimator provided in Vanraden (2008, p. 4416)
As MAF, Fst, LD extent, and consequently CorK_Freq were different in the three panels (Table 1), average power was highly variable among the three panels (adjusted for the same population size). Among the three diversity panels, the C-K panel appeared to be the most powerful on average due to its higher MAF, lesser LD extent, and its lower relatedness. It should be noted that this analytical study assumed that the variance components were known. It was therefore necessary to confirm these results with simulations.
Simulation-based comparison of type I risk and power of statistical models associated with different estimations of K
Removing the markers on the same chromosome than the one tested (MK_Chr) permitted the decrease of the correlation between the kinship at the tested SNP and the global covariance (CorK_Chr in Figure 4). CorK_Chr remained nevertheless high in structured regions (high Fst), i.e., regions with important differentiation between genetic groups (Figure 4, A and C, chromosome 8), which suggests that K_Chr was efficient to estimate covariance between individuals.
To evaluate models involving different kinship estimators for their ability (i) to control false positives at nominal levels and (ii) to detect QTL, we conducted simulations based on the genotypes of the diversity panels. Using all the markers to estimate the kinship matrix (MK_Freq) led to an overcorrection of the H0–P-values (Table 2 and Figure 5). This was particularly the case in the panel with the highest level of LD (CF-Dent). Under H0, the P-value distributions of the two alternative models were much closer to the expected distribution and revealed that these approaches were also efficient to control false positives (Figure 5). Results obtained with MK_Chr showed that molecular information carried by 9 of the 10 chromosomes was sufficient to reliably estimate covariance between individuals to control for false positives.
Knowing that the three estimations of the kinship matrix (K_Freq, K_Chr, and K_LD) were efficient to control false positives, we could compare their power in a second step of simulations. QTL were sampled from the 10 chromosomes, and the power of MK_Freq, MK_Chr, and MK_LD at different thresholds was evaluated at the SNPs/QTL. The alternative models were more powerful than the usual model MK_Freq (Table 3). In particular, estimating the covariance matrix using the markers on the nontested chromosome (MK_Chr) resulted in higher power in each scenario in each panel. As expected the gain of power was higher in the panel with more extended LD (CF-Dent). The gain of power was lower with MK_LD, but we suppose that this approach could be improved by taking into account gene density along the genome, or a priori information on genetic architecture, and by using a better estimate of the covariance between the marker-based estimators when computing optimal marker weights. Note that further research on the K_LD estimator should also consider its scalability when applied to very-high-dimensional data sets.
To check the stability of these results, when considering that a QTL could be detected by SNPs located near it, we used another simulation approach, in which SNPs within a genetic window around the QTL positions were considered as H1 markers and the others as H0 markers. The results (Table S1) confirmed that at a given realized FDR, the alternative models and in particular MK_Chr were more powerful than the traditional model (MK_Freq). Considering that true discoveries were within 5 cM of the QTL, MK_Freq had a power to detect QTL of 11%, MK_Chr of 26% and MK_LD of 19% at a realized FDR of 10%, when 100 QTL were simulated in the CF-Dent panel.
In conclusion, the derivation of analytical power permitted highlighting which parameters are linked to power in Association Mapping. In particular the kinship between individuals (K) clearly influenced the noncentrality parameter. Analytical power scan in three diversity panels also confirmed that the way of estimating K can affect power. In particular, the usual model (MK_Freq) has a low power in regions of high LD. We considered two alternative approaches to recover this gap of power, and we could show with simulations based on true genotypes that they were more powerful than the usual models at given type I risks.
L. Moreau, T. Mary-Huard and A. Charcosset conducted this research in the framework of Amaizing Investissement d’Avenir program. The authors thank the reviewers and the editor for their comments which improved the manuscript. This research was jointly supported as “Cornfed project” by the French National Agency for Research (ANR), the German Federal Ministry of Education and Research (BMBF), and the Spanish ministry of Science and Innovation (MICINN, research project EUI2008-3635). R. Rincent is jointly funded by Limagrain, Biogemma, Kleinwanzlebener Saatzucht, and the French ANRt.
Communicating editor: N. A Rosenberg
- Received November 15, 2013.
- Accepted February 9, 2014.
- Copyright © 2014 by the Genetics Society of America