Statistical models
Mixed models allow to correct for dependence between observations due to genetic structure. Yu et al. (2006) defined a “unified mixed model”, also known as the \(Q + K\) model [4], that can accommodate both a population structure matrix (\(Q\)) and a kinship matrix (\(K\)):
$$\begin{array}{*{20}c} {y = X{\varvec{\beta}} + Q{\varvec{v}} + \underline{{Z{\varvec{u}}}} + \user2{\underline {\varepsilon } }\quad Var\left( {\varvec{u}} \right) = K\sigma_{G}^{2} \quad Var\left( {\varvec{\varepsilon}} \right) = R\sigma_{\varepsilon }^{2} } \\ \end{array}$$
(1)
where \(y\) is the vector of phenotypic trait values, \(X{\varvec{\beta}}\) represents the incidence matrix and marker effects (SNP effect in [6]); \(Q{\varvec{v}}\) are the population structure matrix and vector, respectively; \(\underline{{Z{\varvec{u}}}}\) are design matrix and vector of genetic background effects (polygene component in [6]); and \(\user2{\underline {\varepsilon } }\) is the residuals vector. The variances of the random effects, \({\varvec{u}}\) and \({\varvec{\varepsilon}}\) are also defined: \(K\) is the kinship matrix and \(\sigma_{G}^{2}\), the genetic variance; \(R\) is a matrix with off-diagonal numbers being 0 and the diagonal is the reciprocal of the number of observations underlying each genotype estimation, and \(\sigma_{\varepsilon }^{2}\) is the residual variance.
Fixed term: allele parametrization
Definition of \(X\) requires a genetic model, that is, a method to transform genetic data into an incidence matrix \(X\). Polyploid genetic models have existed for a long time [15] and have inspired more recent versions applied to SNP data [16, 17]. The simplest of them is the biallelic model (model B in [7], association mapping in [18]), which considers SNP alleles as equal to QTL alleles. In a biallelic model, the SNP dosages are used to predict genetic effects, giving the \(X\beta\) term the following form:
$$\begin{array}{*{20}c} {X_{b} \beta = \left[ {\begin{array}{*{20}l} 1 \hfill & {\delta_{1} } \hfill \\ 1 \hfill & {\delta_{2} } \hfill \\ \vdots \hfill & \vdots \hfill \\ 1 \hfill & {\delta_{n} } \hfill \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} \mu \\ \beta \\ \end{array} } \right] } \\ \end{array}$$
(2)
where \(\delta_{i}\) are the dosages (a value from 0 to ploidy) of one of the SNP alleles, \(\mu\) is the intercept and \(\beta\) the genetic effect of that SNP allele. We denote the incidence matrix as \(X_{b}\) for this modelling strategy. Note that this represents an additive model without intra or inter-locus interaction, i.e. no dominance or epistasis between alleles is modelled.
Alternatively, Identity-By-Descent (IBD) information can be used to generate an ancestral model [13], also known as a PBA model [10] or an LDLA model [9, 19]. Under the ancestral model, the dosage of each ancestral allele or haplotype in the NAM population is used to estimate genetic effects. The shape of the \(X\beta\) term then takes the form:
$$\begin{array}{*{20}c} {X_{a} \beta = \left[ {\begin{array}{*{20}l} 1 \hfill & {\delta_{11} } \hfill & {\delta_{12} } \hfill & \hfill & {\delta_{1k} } \hfill \\ 1 \hfill & {\delta_{21} } \hfill & {\delta_{22} } \hfill & \hfill & {\delta_{2k} } \hfill \\ \vdots \hfill & \vdots \hfill & \vdots \hfill & { \ldots } \hfill & \vdots \hfill \\ 1 \hfill & {\delta_{n1} } \hfill & {\delta_{n2} } \hfill & \hfill & {\delta_{nk} } \hfill \\ \end{array} } \right]\left[ {\begin{array}{*{20}l} \mu \hfill \\ {\beta_{1} } \hfill \\ {\beta_{2} } \hfill \\ \vdots \hfill \\ {\beta_{k} } \hfill \\ \end{array} } \right]} \\ \end{array}$$
(3)
In this case, the dosages of all alleles except one (the reference allele) are specified. Therefore, \(k\) is the number of alleles \(- 1\). Each \(\beta\) represents the additive genetic effect of each ancestral allele, relative to the effect of the reference ancestral.
Random term: kinship matrix calculation
In this model, a kinship matrix \(K\) is calculated using the realized relationship [4]:
$$K = \frac{{DD^{t} }}{{\Delta }}\quad {\Delta } = \overline{{diag\left( {DD^{t} } \right)}}$$
where \(D\) is a dosage matrix with markers on columns and individuals on rows, and the mean of each column is zero (column means have been subtracted for each column); and \({\Delta }\) is the mean of the diagonal of the \(DD^{t}\) matrix. If haplotypes are used instead of biallelic SNPs, \(D\) can consist of concatenated matrices similar to \(X_{a}\) (without the intercept column), so that the number of columns is equal to the total number of alleles present across all markers used. To mitigate the bias due to differences in marker density across the genome, kinship estimates are calculated on a subset of evenly distributed SNPs (one marker per cM).
Haplotyping
Haploblocks were arbitrarily defined using a sliding window of 6 consecutive SNPs with an overlap of 4 SNPs (first haplotype is SNP1-SNP2-…-SNP6, second is SNP3-SNP4…-SNP8). A haploblock of length 6 can tag a maximum of \(2^{6} = 64\) alleles if all combinations are present, although in our simulations the number of observed alleles was much lower, with the average number of observed alleles ranging from 11.23 in NAM1 to 21.8 in NAM10. To obtain a haploblock position, the average position of the 6 SNP markers was taken. Haplotypes were obtained from the simulated phased SNP genotypes generated by PedigreeSim.
Power study
Definition of QTL interval
Single marker QTL methods do not provide an estimate for the QTL interval, yet with a defined threshold and a genetic map one can interpret the p value distribution to obtain them. Since adjacent markers are not independent, and the closer to a true QTL position, the more significant the p-value becomes, one expects a chain of increasingly significant markers, pointing towards a true QTL position. Based on this, we define a QTL interval as a set of ordered markers above the significance threshold such that:
$$QTL = \left\{ {m_{1} , \ldots , m_{n} } \right\}\quad {\text{where}}\;d_{ij} < l$$
where \(d_{ij}\) is the distance between adjacent significant markers \(i\) and \(j\), and \(l\) represents a linking distance. As a result, a QTL interval is defined by a chain of significant markers, where adjacent significant markers are at a distance smaller than \(l\). Therefore, for each value of \(l\) we can define a set of detected QTL intervals. Since the choice of \(l\) is arbitrary, we performed power calculations with \(l\) from 0 to 10 cM in steps of 0.5 cM.
Significance threshold
To adjust for multiple testing, an empirical permutation threshold was calculated for each QTL analysis [20]. Thresholds were obtained with 100 permutations on a single population for each model, as threshold values did not change substantially between populations.
Power estimates
To evaluate the QTL models here presented we will use (1) QTL detection power, the probability of detecting a QTL position when present; (2) false positive rate, the probability of having a significant marker outside a QTL region; (3) QTL accuracy, the closeness of a QTL peak (position of maximum probability within an interval) to the true position and (4) QTL and marker precision, the probability that a significant QTL interval or marker is a true positive.
QTL detection power can be calculated as the proportion of true QTLs that are found by the model. While this is informative, one can easily increase detection power by increasing the number of false positives. To estimate the false positive rate, we must define the true negative markers (N). We considered as true negatives all markers outside a 10 cM interval around our true QTL positions (5 cM above and 5 cM below). We then define as false positives (FP) those markers that are above the significance threshold (they have lower p values, higher significance) and are outside the 10 cM true interval. Lastly the false positive rate is calculated as \(FP/N\).
The range of a QTL interval is defined by the positions of its leftmost and rightmost markers. QTL intervals will be considered true positives if the QTL range includes the simulated QTL position. All markers belonging to a true positive QTL interval are considered true positive markers, whereas the rest of significant markers present in other QTL intervals will be considered false positives. Isolated significant markers will be ignored.
Under this framework we can define detected QTLs, true QTLs, significant markers and true positive markers. We will use these values to calculate the precision (proportion of true positives over all positives) for both QTLs and markers.
$$QTL_{precision} = \frac{true\;positive\;QTLs}{{detected\;QTLs}}$$
$$marker_{precision} = \frac{true\;positive\; markers}{{significant\; markers}}$$
Finally, we considered the ability of a model to predict the position of QTL within an interval. We can define a QTL peak as the most significant marker within a QTL interval, as is done when applying logarithm of odds (LOD) scores. QTL accuracy can then be calculated as the average distance of a QTL peak in a true QTL to the true QTL position.
Power measures were calculated for each of the three models in 11 populations for each level of genetic diversity (total of 44 populations).
Implementation
All computations in this study were done in R [14].
Ridge regression using a restricted maximum likelihood procedure was used to obtain the mixed model estimates, which in this context are equivalent to the Best Linear Unbiased Predictions (BLUP) [21, 22]. Such calculations can be performed using the mpQTL package, where the solution algorithm, F-test approximation and p value calculation where based on the mixed.solve() function of the rrBLUP package [23].
To improve computational efficiency, the EMMAX/P3D approach was applied [24, 25], which approximates variance components once, and recycles these components at each marker position, reducing the amount of large matrix multiplications that must be performed.
Simulation
Multiparental population design and genotype simulation
Nested Association Mapping (NAM) populations were generated using PedigreeSim V2.0, a simulation software that can simulate not only diploid but also polyploid meiosis [26]. PedigreeSim generates genotypes given a genetic map, a pedigree and the genotypes of the first generation (founders) of that pedigree. Simulations were performed using Haldane’s mapping function, allowing only bivalents with random pairing and the parameter “NATURALPAIRING” set to 1.
To speed up the calculations, an adapted tetraploid potato genetic map was used [27] containing only the first five chromosomes (3509 markers representing 485 cM). The individuals used in this study were simulated in a two-stage process: firstly, ancestor individuals were generated and used to obtain ten separate populations (ancestral groups); secondly, from each ancestral group a set of NAM founders were chosen to obtain parallel NAM populations.
For each ancestral group (AG), 10 ancestor individuals were generated with random SNP scores at each marker. Each SNP position is also given an “IBD allele”, unique to each homologue of each ancestor (even if the SNP state is the same). Each ancestral group has 10 founders, and thus a total of 40 IBD alleles will segregate in each AG. These alleles we will name ancestral alleles. Each ancestor is randomly crossed (without selfing, as potato is an outbreeder) to produce a first generation of 100 individuals, which will serve as parents of the second generation. This process was repeated for 50 generations, maintaining a constant generation size of 100 individuals. Finally, 100 individuals per AG were obtained as potential parents for the creation of NAM populations.
A NAM population consisted of one central parent crossed with nine peripheral parents, without any of the subsequent inbreeding that was originally proposed for NAM crossing scheme for selfing crops [12]. Each cross produced 50 offspring, thus totalling at 460 individuals per NAM. To simulate NAMs with different degrees of genetic diversity, parents were sampled from the same or from different AGs. A NAM1 contains parents from only one AG, while a NAM5 contains parents from 5 different AGs, with the same number of parents per group when possible. When the numbers of parents per AG was not equal the central parent always originated from the AG providing the most parents. For each level of genetic diversity, 11 populations were simulated. At the end of the process, the genotypes of each individual were obtained in terms of ancestral alleles (IBD alleles) and in terms of SNP dosages.
Phenotype simulation
Phenotypes were simulated based on the simulated genotypes: genotypic values were obtained by assigning genetic effects to the ancestral alleles at pre-defined QTL positions. Each individual will then harbour four QTL alleles at each QTL position and the final phenotype is equal to the added effects of all QTL alleles plus a normally distributed noise. No interactions between alleles in one QTL or among QTL loci were simulated, and thus additive phenotypes were obtained.
We considered a situation where three unique QTL positions (at chromosome 1, 67.88 cM; chromosome 2, 61.2 cM and chromosome 4, 100.49 cM) were segregating. Each AG has a random allelic mean, and allele effects are drawn from a normal distribution around that mean. Additionally, 50 small-effect QTLs were added randomly across the genome to simulate a polygenic effect.
For further information see Additional file 1.