Introduction

The R package qtlpoly (v. 0.2.1) (Pereira et al. 2020) is an under development software to map multiple quantitative trait loci (QTL) in full-sib families of outcrossing autopolyploid species. It is based on the following random-effect model:

\[\boldsymbol{y} = \boldsymbol{1}\mu + \sum_{q=1}^Q \boldsymbol{g}_q + \boldsymbol{e}\]

where the vector of phenotypic values from a specific trait \(\boldsymbol{y}\) is a function of the fixed intercept \(\mu\), the \(q = 1, \dots, Q\) random QTL effects \(\boldsymbol{g}_q \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{G}_q\sigma^2_q)\), and the random environmental error \(\boldsymbol{e} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I}\sigma^2)\). \(\boldsymbol{G}_q\) is the additive relationship matrix between full-sibs derived from the genotype conditional probabilities of QTL \(q\). See model details in Pereira et al. (2020).

Variance components associated with putative QTL (\(\sigma^2_q\)) are tested using score statistics from the R package varComp (v. 0.2-0) (Qu, Guennel, and Marshall 2013). Final models are fitted using residual maximum likelihood (REML) from the R package sommer (v. 3.6) (Covarrubias-Pazaran 2016). Plots for visualizing the results are based on ggplot2 (v. 3.1.0) (Wickham 2016).

This tutorial used R version 4.0.3 (2020-10-10) running on either Ubuntu 18.04.1 LTS (64-bit) or Windows 10.

Install and load the qtlpoly package and data

As mentioned, the package qtlpoly depends on a couple of functions from sommer (v. 3.6) and varComp (v. 0.2-0). varComp has been archived from CRAN, while sommer has been constantly updated. In order to avoid conflict with updates in functions and object structures, we decided to stick with an earlier sommer version (v. 3.6), so you will need to downgrade if you have installed the most recent version of sommer.

qtlpoly package is available at GitHub. You can install all needed packages within R using the functions from the R package devtools:

Then, use the function library() – or require() – to load the package, and the function data() to preload the simulated linkage map and phenotypic data that comes with the package:

These simulated data were made simpler, so that one could run the functions along with this tutorial in a regular personal computer (with 4 cores and 6 GB of RAM, minimum). For real data analyses, you may need to run an R script in a cluster with more cores and RAM, though. In general, the computational needs depend on ploidy level, population size and number of markers.

The data was based on a full-sib family (\(N = 300\)) simulated from a cross between two parents from a putative autohexaploid species with three chromosomes. The simulated linkage map consists of three linkage groups (LGs) with 1,310 markers, spanning 274.76 centiMorgans (cM) in total. Please see details on the map construction for this simulated cross using mappoly here.

The phenotypic data contains three traits, and it is a simple table with traits in columns and individuals in rows. For example:

Trait means were simulated as zero, i.e. \(\mu = 0\) . Environmental variances were simulated from a standard normal distribution, i.e. \(\sigma^2_e = 1\). QTL were simulated as described below:

  • ‘T32’: 3 QTL, with heritabilities of 0.20 (LG 1 at 32.03 cM), 0.15 (LG 1 at 95.02 cM) and 0.30 (LG 2 at 40.01 cM), i.e. \(\sigma^2_q = \{0.60, 0.45, 0.80\}\);
  • ‘T17’: 1 QTL, with heritability of 0.17 (LG 3 at 34.51 cM), i.e. \(\sigma^2_q = 0.20\);
  • ‘T45’: no QTL.

Conditional probabilities

Our qtlpoly package is fully integrated with mappoly (Mollinari and Garcia 2019), so that if you have a map built using another software, you will need to re-estimate it using functions from mappoly – or, at least, mimic the ‘mappoly.map’ object structure (good luck with that!). Please refer to its tutorial in order to learn more about the software.

To install mappoly, simply use the following command:

mappoly uses hidden Markov model (HMM) adapted from Lander and Green (1987) to autopolyploids to estimate the genotype conditional probabilities of putative QTL given an ordered set of markers comprising a linkage map (Mollinari and Garcia 2019). These conditional probabilities are ultimately used for QTL mapping by qtlpoly.

To compute the conditional probabilities for each linkage group, one should use the function calc_genoprob() from mappoly:

Prepare data

Once the conditional probabilities have been calculated, you can use the function read_data() to read both ‘genoprob’ and ‘pheno’ objects. Row names (individuals) on ‘pheno’ must match the individual names from the ‘genoprob’ object. The argument ploidy = 6 informs the ploidy level of the cross, i.e. a hexaploid cross in this case. Finally, step = 1 provides the step size, so that only positions at every 1 cM (approximately) will be tested:

If you want to see detailed information on the map (number of positions per linkage group) and the phenotypes (number of individuals with non-missing values per trait), use a print() function with detailed = TRUE:

Notice that regardless of the number of markers per linkage group, we will rather evaluate only a limited number of positions, given the step argument (in our example, every 1 cM). Because of the quite large linkage disequilibrium blocks in mapping populations, specially in full-sib families with relatively small sample sizes (let us say, \(N < 300\)), there is no need to test every single marker in the map. Moreover, with such marker density (e.g., 5+ markers/cM), a lot of markers will contain redundant information for QTL mapping purposes anyway. step = NULL will allow to test every single marker in a high density map if one wants to do so, but be sure that computational time (and needed RAM for loading ‘qtlpoly.data’ object) will be hugely increased.

Perform QTL detection

Building a multiple QTL model is considered a model selection problem, and there are several ways to approach it. Here, we adapted the algorithm proposed by Kao, Zeng, and Teasdale (1999) for fixed-effect multiple interval mapping (MIM) for diploids to our random-effect MIM (REMIM) for polyploids, which is summarized as follows:

  1. Null model: for each trait, a model starts with no QTL \[\boldsymbol{y} = \boldsymbol{1}\mu + \boldsymbol{e}\]
  2. Forward search: QTL (\(q = 1, \dots, Q\)) are added one at a time, conditional to the one(s) (if any) already in the model, under a less stringent threshold (e.g., \(P < 0.01\)) \[\boldsymbol{y} = \boldsymbol{1}\mu + \sum_{q=1}^Q \boldsymbol{g}_q + \boldsymbol{e}\]
  3. Model optimization: each QTL \(r\) is tested again conditional to the remaining one(s) in the model under a more stringent threshold (e.g., \(P < 10^{-4}\)) \[\boldsymbol{y} = \boldsymbol{1}\mu + \boldsymbol{g}_r + \sum_{q \neq r} \boldsymbol{g}_q + \boldsymbol{e}\] Steps 1 and 2 are repeated until no more QTL can be added to or dropped from the model, and positions of the remaining QTL do not change. After the first model optimization, the following forward searches use the more stringent threshold (e.g., \(P < 10^{-4}\)) as the detection power is expected to increase once QTL have already been added to the model.
  4. QTL profiling: score statistics for the whole genome are updated conditional to the final set of selected QTL.

It is worth to mention that there is no such definitive way to perform model selection. One may try to come up with alternative ways of building a multiple QTL model, some of which will be mentioned later. qtlpoly tries to provide functions flexible enough, so that the users can build a multiple QTL model on their own, manually. One strategy we have used the most and tested through simulations implements the above-mentioned algorithm of building a multiple QTL model, and will be demonstrated here.

Add/drop QTL to/from the model

Suppose you have information (from previous work, for instance) that can help you to decide whether a QTL should be in the model or not. You can take advantage of the function modify_qtl() to either include or exclude a certain putative QTL from the model. For example, let us say we want to add a QTL to the trait 3 ‘T17’ model on LG 3 at 13.03 cM (position number 184), which has not reached the most stringent threshold when a second search was performed after (see the results from the step 2. Model optimization). We simply choose the pheno.col and the position where we want a QTL to be added, i.e. add.qtl = 184:

The function print() shows how the model has been modified:

Notice that since the model has not been profiled yet, statistics may change.

Now, suppose you want to drop the least significant QTL from the trait 1 ‘T32’ model. You would have to use drop.qtl = 100 instead:

Once modifications have been suggested, run the function profile() in order to obtain updated statistics:

Using print(), notice how positions and statistics may have changed:

It is worth to mention that a modification should be considered carefully, as we may not have enough statistical support to hold it. Nevertheless, modification might be interesting to study how much of the phenotypic variance an added QTL can explain, regardless of its \(P\)-value, for example.

Plot profiles

Given the final profiled models, you can plot either individual or joint \(LOP\) profiles using the function plot_profile(). For individual plots, one need to specify the pheno.col argument:

Triangles show where the mapped QTL peaks are located. By providing values to the y-axis limits using ylim, one makes sure that all plots will be at the same scales.

If one wants to plot some or all QTL profiles at once for comparison, pheno.col needs to be specified as a vector of numbers representing the phenotype columns. As default, pheno.col = NULL includes all evaluated traits:

The argument grid organizes the multiple plots as a grid if TRUE, or superimposed profiles if FALSE.

Support intervals

You can visualize the QTL distributed along the linkage map, together with their support intervals using the function plot_sint():

Notice that for highly significant QTL (\(P < 2.22 \times 10^{-16}\)), computing support intervals may not be very reliable, which is expected to be narrower than the support interval for the QTL on LG 2 for ‘T32’. We used the same colors as in plot profiles to keep them consistent through plots. Obviously, no QTL are shown for the trait 3 ‘T45’.

Fit multiple QTL models

Once final models have been defined, one may use REML to estimate their parameters by using the function fit_model():

Here, we used the joint probabilities of the putative QTL genotypes (probs= "joint"), instead of the parent marginal ones (probs = "marginal"). The latter may be interesting if one of the parents seems to not contribute much for a specific QTL (because all QTL allele effects are about the same, for example). Also, if polygenes = "none", all individual QTL have their variance component estimated. One may choose to do polygenes = "most" if all but one QTL have to be combined as a polygenic effect, or polygenes = "all" if all QTL have to be combined as a polygenic effect. Finally, keep = TRUE (the default) saves matrices used on the calculations and additional results from sommer. One may want to assess standard errors from parameter estimates there, for example.

A summary() function shows parameter estimates together with the QTL heritabilities computed from the variance component estimates:

At this point, we are able to compare our results with the simulated values presented at install and load the qtlpoly package and data section:

  • ‘T32’: 3 QTL were detected, with heritabilities of 0.19 (LG 1 at 33.07 cM), 0.13 (LG 1 at 104.28 cM) and 0.28 (LG 2 at 41.20 cM), when we expected 0.20 (LG 1 at 32.03 cM), 0.15 (LG 1 at 95.02 cM) and 0.30 (LG 2 at 40.01 cM), respectively;
  • ‘T17’: 1 QTL was detected, with heritability of 0.19 (LG 3 at 29.19 cM), when we expected 0.17 (LG 3 at 34.51 cM);
  • ‘T45’: no QTL were detected, as expected.

Plot QTL

After models have been fitted, one way of visualization you can use is the one provided by the function plot_qtl():

Dots are located on the respective chromosome positions of the QTL peaks. Size of the dots corresponds to the specific QTL heritability. Color of the dots corresponds to the \(P\)-values, and helps to identify the most or the less significant QTL. drop = FALSE makes sure all traits are displayed, even if they do not had any QTL detected.

Estimate allele effects

Additive effects contributing to the overall mean by each allele individually as well as their combinations within each parent may be computed using the function qtl_effects() as follows:

A plot() function allows the user to visualize these contributions graphically:

For each QTL, one will be able to see which alleles contribute the most to increase or decrease the phenotypic mean. This is interesting to plan marker-assisted selection strategies, for example, as specific alleles can be selected based on the marker dosage at or around the QTL peak and with assistance of the population linkage map from mappoly. For example, if one is interested in decreasing the value of the trait ‘T17’, the allele i, combined with j and l from the parent ‘P2’ must be targeted. On the other hand, the parent ‘P1’ does not seem to contribute much to the trait variation, although selecting for the alleles b through e must be the best way to achieve one’s goal of decreasing the phenotypic value of ‘T17’.

Predict breeding values

Finally, with the estimates of the final models in hands, one can use them to perform predictions of the breeding values as follows:

A plot() function shows the distribution of the genotypic values for the population:

This may also be interesting for those populations whose individuals have been genotyped, but not phenotyped, and you still want to consider them when selecting the best genotypes.

Compare with a fixed-effect approach

A previous fixed-effect interval mapping (here named FEIM) model has been proposed as a first approach to map QTL in autopolyploid species (Hackett, Bradshaw, and Bryan 2014). It consists of a single-QTL model, where every position is tested according to the model:

\[Y = \mu_C + \sum_{i=2}^{m} \alpha_i X_i + \sum_{i=m+2}^{2m} \alpha_i X_i\]

where \(\mu_C\) is the intercept, and \(\alpha_i\) and \(X_i\) are the main effects and indicator variables for allele \(i\), respectively, and \(n\) is the ploidy level. The constraints \(\alpha_{1} = 0\) and \(\alpha_{m+1} = 0\) are imposed to satisfy the condition \(\sum_{i=1}^m X_i = m/2\) and \(\sum_{i=m+1}^{2m} X_i = m/2\), so that \(\mu_C\) is a constant hard to interpret due to these constraints. Notice that the higher the ploidy level, the more effects have to be estimated, i.e. tetraploid models have six main effects, hexaploid models have 10 effects, octoploid models will have 14 effects (i.e. \(2m-2\)).

Different from REMIM (where we test the variances), here the interest is to know if the average allele effects are different from zero (the null hypothesis) using likelihood-ratio tests (LRT). Commonly, the tests are presented as “logarithm of the odds” (LOD scores), where \(LOD = \frac{LRT}{2 \times \log_e(10)}\).

In order to evaluate significance (declare a QTL), empirical LOD thresholds are computed for each trait using permutations as proposed by Churchill and Doerge (1994).

LOD threshold permutations

Using the same object ‘data’ from REMIM analyses, one can first compute the thresholds for all or specific traits. The number of simulations is given by n.sim, and 1,000 permutations are generally used:

Once parallel analyses using 4 clusters (n.clusters = 4) are done, one may print the results, specifying a vector probabilities of interest. By default, probs = c(0.95, 0.90), so that 95% and 90% quantiles are shown:

As probs = c(0.95, 0.90) were also given as default in permutations() function, we already stored the genome-wide significance LOD threshold for 0.05 and 0.10 levels, respectively. In order to use the 95% quantiles in the subsequent FEIM analyses, one can do:

Interval mapping

feim function tests every position from the specified step size from ‘data’ – here, every 1 cM. Besides the sig.lod vector containing the thresholds for each trait, one needs to provide a window size (e.g., w.size = 15), which will be used to select QTL peaks within the same linkage group with a minimum distance of the given window size:

A print function shows detailed information on the detected QTL:

Remember, in this case, one should not sum adjusted \(R^2\) from the same trait, as each was obtained from a single-QTL model.

Finally, one may want to plot the profiles and compare with the plot profiles from REMIM:

Notice that one QTL on LG 1 was not detected for the trait ‘T32’ (false negative), while one QTL on LG 3 for the trait ‘T45’ was wrongly assigned (false positive). This exemplifies the power of multiple-QTL models over the single-QTL ones. Therefore, the FEIM model may be recommended only as a first, quick approach, but not as the ultimate model for detecting QTL in autopolyploid species.

Acknowledgments

This package has been developed as part of the Genomic Tools for Sweetpotato Improvement (GT4SP) project, funded by Bill & Melinda Gates Foundation.

References

Churchill, G A, and R W Doerge. 1994. “Empirical threshold values for quantitative trait mapping.” Genetics 138 (3): 963–71.

Covarrubias-Pazaran, Giovanny. 2016. “Genome-assisted prediction of quantitative traits using the R package sommer.” PLoS ONE 11 (6): 1–15. https://doi.org/10.1371/journal.pone.0156744.

Hackett, Christine A., John E. Bradshaw, and Glenn J. Bryan. 2014. “QTL mapping in autotetraploids using SNP dosage information.” Theoretical and Applied Genetics 127 (9): 1885–1904. https://doi.org/10.1007/s00122-014-2347-2.

Kao, C-H, Z-B Zeng, and R D Teasdale. 1999. “Multiple interval mapping for quantitative trait loci.” Genetics 152 (3): 1203–16.

Lander, E S, and P Green. 1987. “Construction of multilocus genetic linkage maps in humans.” Proceedings of the National Academy of Sciences of the United States of America 84: 2363–7.

Mollinari, Marcelo, and Antonio Augusto Franco Garcia. 2019. “Linkage Analysis and Haplotype Phasing in Experimental Autopolyploid Populations with High Ploidy Level Using Hidden Markov Models.” G3: Genes|Genomes|Genetics 9 (10): 3297–3314. https://doi.org/10.1534/g3.119.400378.

Pereira, Guilherme da Silva, Dorcus C. Gemenet, Marcelo Mollinari, Bode A. Olukolu, Joshua C. Wood, Federico Diaz, Veronica Mosquera, et al. 2020. “Multiple QTL Mapping in Autopolyploids: A Random-Effect Model Approach with Application in a Hexaploid Sweetpotato Full-Sib Population.” Genetics 215 (3): 579–95. https://doi.org/10.1534/genetics.120.303080.

Qu, Long, Tobias Guennel, and Scott L. Marshall. 2013. “Linear score tests for variance components in linear mixed models and applications to genetic association studies.” Biometrics 69 (4): 883–92. https://doi.org/10.1111/biom.12095.

Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer.