Short code for demo

Please, copy-and-paste the code below to an R script file in RStudio: File > New File > R Script. See installing the qtlpoly package section if you have not installed it yet.


Quantitative trait loci (QTL) mapping studies main goal is to investigate the genetic architecture of traits of interest. Single- and multiple-QTL models have been available for inbred, diploid mapping populations for quite some time now (see Da Costa E Silva and Zeng (2010) for a comprehensive review). However, only recently these methods became available for outbred, polyploid mapping populations.

The R package qtlpoly (v. 0.2.1) is an under development software to map multiple QTL in full-sib families of outcrossing, autopolyploid species (G. da S. Pereira, Gemenet, et al. 2020). It is based on the partition of the phenotypic variance (\(\sigma^2_p\)) into variance due to QTL (\(\sigma^2_q\)) in addition to the residual variance (\(\sigma^2_e\)) as follows:

\[\sigma^2_p = \sum_{q=1}^Q \boldsymbol{G}^{(i,i')}_q\sigma^2_q + \sigma^2_{e}\]

where \(\boldsymbol{G}^{(i,i')}_q\) is the additive relationship between full-sibs \(i\) and \(i'\), whose computation is based on the genotype conditional probabilities of QTL \(q\). See how to estimate these genotype probabilities using mappoly in Mollinari et al. (2020).

Because we only need to estimate one parameter per QTL (the very variance component associated with it), it is relatively easy to look for additional QTL and add them to the variance component model, without ending up with an overparameterized model. A multiple-QTL model is known to have increased power when compared to a single-QTL model, with ability to detect minor or separate linked QTL (G. da S. Pereira, Gemenet, 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 was first developed for the Polyploid Tools Training Workshop (December 12-15, 2021), and tested in R version 4.0.3 (2020-10-10) running on either Ubuntu 18.04.1 LTS (64-bit) or Windows 10.

Installing the qtlpoly package

In order to run this tutorial, you will need to install the qtlpoly package, which is available at GitHub. You can install all needed packages within R using the commands below:

If you have installed qtlpoly before, please make sure to run the last line of the installation process again, so that you will have all the last updates. Then, use the function library() – or require() – to load the package:

Tetraploid potato data

A cross between two potato (Solanum tuberosum, 2n = 4x = 48) cultivars, ‘Atlantic’ \(\times\) ‘B1829-5’, resulted in 156 full-sibs. The population has been phenotyped (4-year evaluation) and genotyped (8k SNP array), and analyses have been performed to call SNP dosage, build a genetic map and map QTL (G. da S. Pereira, Mollinari, et al. 2020).

For brevity’s sake, we have selected three phenotypes (foliage maturity evaluated in years 2007, 2008 and 2014, i.e. FM07, FM08 and FM14) and three linkage groups (LGs, namely 1, 5 and 7) for this demo. Foliage maturity is measured as the “area under the curve” of foliage color along the plant cycle. All analyses can be found at this GitHub page, though.

Use the function data() to preload the genotype probabilities computed in the LGs 1, 5 and 7 (now called 1, 2 and 3, respectively) as well as the phenotypic data:

The genoprob4x object contains the genotype probability along the LGs for each individual, reflecting all recombination events. For example, individuals 1 and 2 show the following patterns of genotype probabilities, and so will every other individual:

The pheno4x object contains the phenotypic values of foliage maturity for years 2007, 2008 and 2014:

In a region where a QTL exists, the more alleles a pair of individuals share, the more their phenotypes will look alike. In fact, this is the basis of QTL detection.

Performing QTL detection

A algorithm proposed for QTL model selection (Kao, Zeng, and Teasdale 1999) is available within the software, and it can be described as follows:

  1. Null model: for each trait, a model starts with no QTL \[\sigma^2_p = \sigma^2_{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 genome-wide significance level (e.g., \(\alpha < 0.20\)) \[\sigma^2_p = \sum_{q=1}^Q \boldsymbol{G}^{(i,i')}_q \sigma^2_q + \sigma^2_{e}\]
  3. Model optimization: each QTL \(r\) is tested again conditional to the remaining one(s) in the model under a more stringent genome-wide significance level (e.g., \(\alpha < 0.05\)) \[\sigma^2_p = \boldsymbol{G}^{(i,i')}_r \sigma^2_r + \sum_{q \neq r} \boldsymbol{G}^{(i,i')}_q \sigma^2_q + \sigma^2_{e}\] Steps 1 and 2 iterate 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 round, the following forward searches use the more stringent threshold (e.g., \(\alpha < 0.05\)).
  4. QTL profiling: score statistics for the whole genome are updated conditional to the final set of selected QTL.

Score-based resampling to assess genome-wide significance

Rather than guessing pointwise significance levels for declaring QTL, you can use the score-based resampling method to assess the genome-wide significance level (Zou et al. 2004). This method is relatively intensive, because it involves score statistics computation for every position in the map repeated 1,000 times (resampling):

Therefore, we ran it in advance and learned that genome-wide significance levels of \(\alpha=0.20\) and \(\alpha=0.05\) match \(P < 0.0011493\) and \(P < 2.2844648\times 10^{-4}\), respectively, which will be used next. Remember you are supposed to run the resampling method for all the 12 linkage groups together (for our potato example).

Searching for QTL

Now, we are to use the remim function for building our multiple QTL model. One should include the data object, the window size (w.size = 15), the value \(d\) that is decreased from the log of \(P\)-value to compute the support interval (d.sint = 1.5 for ~95%), and the number of cores to be used (n.clusters = 6).

In case you have computed the ‘score.null’ object, both forward search (sig.fwd = 0.20) and backward elimination (sig.bwd = 0.05) will reflect the desired genome-wide significance levels (\(\alpha\)):

Otherwise, you can just include the pointwise significance level computed before based on our resampling method, which is the option we are going to use here:

Use print() and a summary table for each trait will be shown:

Since support intervals have been calculated, you can print them as well by specifying the sint argument properly:

Visualizing results graphically

Given the final profiled models, you can plot either individual or joint \(LOP\) profiles using the function plot_profile(). Triangles show where the mapped QTL peaks are located:

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

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

Fitting multiple QTL models

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

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

Another way of visualization you can use is the one provided by the function plot_qtl():

Dots are located on the respective LG 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.

Estimating 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: