Funmap2: an R package for QTL mapping using longitudinal phenotypes

Statistical methods

Functional mapping is a statistical framework aimed at identifying QTLs that are associated significantly with a longitudinal phenotype of interest in an experimental population, such as recombinant inbred lines (RIL) or a doubled-haploid (DH) population. Functional mapping computes maximum likelihood estimation (MLE) of mixture models that integrate the likelihood over QTL genotypes. The model also accounts for the (1) longitudinal trend of the trait using a continuous curve (i.e., growth trajectory for time-dependent traits and reaction norm for environment-dependent traits), and (2) internal correlation of traits across longitudinal measurements using a covariance matrix.

The model assumes that the longitudinal trend of the trait follows a particular mathematical curve. In a statistical setting, the phenotypes of the trait at all time points follow a multivariate normal distribution, which can be described by the following Eq. (1): (1)${\displaystyle {\mathrm{f}}_{\mathrm{j}}\left({\mathrm{y}}_{\mathrm{i}}\right)=\frac{1}{{\left(2\pi \right)}^{\frac{\mathrm{m}}{2}}{\left|\sum \right|}^{\frac{1}{2}}}exp\left[-{\left({\mathrm{y}}_{\mathrm{i}}-\alpha {\mathrm{X}}_{\mathrm{i}}-{\mathrm{g}}_{\mathrm{j}}\right)}^{\mathrm{T}}{\sum }^{-1}\left({\mathrm{y}}_{\mathrm{i}}-\alpha {\mathrm{X}}_{\mathrm{i}}-{\mathrm{g}}_{\mathrm{j}}\right)∕2\right].}$

Assuming that y_i is a vector of measured values for individual i at m time points, which describes the phenotypic values, g_j denotes the overall mean vector for genotype j that is described as a mathematical curve, X_i denotes the covariate for individual i, α is a vector of coefficient value for each covariate, and∑ is the covariance matrix. Therefore, ${\mathrm{f}}_{\mathrm{j}}\left({\mathrm{y}}_{\mathrm{i}}\right)$ is the probability density function that relates the measured traits of individual i to the combination of the genetic effects contributed by genotype j and covariate effects of individual i. In Functional mapping, we assume the genetic effects of each genotype vary from time to time and follow a trajectory which can be defined by a mathematical curve (g_j), such as the logistic, Legendre or other type curves (Ma, Casella & Wu, 2002). For the logistic curve, (g_j) can be described by (2)${\displaystyle {\mathrm{g}}_{\mathrm{j}}\left(t_i\right)=\frac{a_j}{1+b_j{\mathrm{e}}^{-r_j{t}_i}}}$ where a_j, b_j and r_j are the parameters for genotype j.

The estimated likelihood of genotype j for individual i are summed, weighted by the corresponding conditional probability of the QTL genotype, given the adjacent marker and inbred type. Therefore, the functional mapping method formulates the likelihood calculation of the mixture model as: (3)${\displaystyle \text{logL}\left(\hat{\Omega }\right)=\sum _{\mathrm{i}=1}^{\mathrm{N}}log\left[\sum _{\mathrm{j}=1}^{\mathrm{J}}{\mathrm{p}}_{\mathrm{ij}}{\mathrm{f}}_{\mathrm{j}}\left({\mathrm{y}}_{\mathrm{i}}\right)\right]}$ where each element p_ij indicates the genotypic possibility of subject i for gene j (QQ, Qq, or qq), N denotes the individual number in the experimental population, and $\hat{\Omega }$ denotes the variables in this log-likelihood function that include the covariate coefficient, curve parameters for multiple genotypes, and covariance parameters. In Eq. (3), we prefer to use the log-likelihood function.

The likelihood function for the null hypothesis model, in which QTL does not affect the trait, is built as follows, (4)${\displaystyle \text{logL}\left(\tilde{\Omega }\right)=\sum _{\mathrm{i}=1}^{\mathrm{N}}log\left({\mathrm{f}}_0\left({\mathrm{y}}_{\mathrm{i}}\right)\right)}$ where ${\mathrm{f}}_0\left({\mathrm{y}}_{\mathrm{i}}\right)$ differs from (1) in g _j by assuming the same longitudinal curve for all genotypes. $\tilde{\Omega }$ denotes the estimated parameters in this log-likelihood function.

The goal of functional mapping is to compute the log-likelihood ratio (LR) for each QTL, defined by the following equation, and then to choose the significant QTL with a high value of LR. (5)${\displaystyle \mathrm{LR}2=-2log\left[\frac{L\left(\tilde{\Omega }\right)}{L\left(\hat{\Omega }\right)}\right].}$ Intuitively, the null hypothesis H₀ is that there is no gene that controls the growth process, and the alternative hypothesis H₁ is that growth processes are different across the QTL genotypes. $\tilde{\Omega }$ and $\hat{\Omega }$ are the maximum likelihood estimates of parameters under the hypotheses H₀ and H₁, respectively.

To conduct a log-likelihood ratio test on complicated statistical models like the one used in functional mapping, a permutation test is usually used to derive and to compare against the null distribution (Churchill & Doerge, 1994). Because the permutation test is generally applicable to various models, it is intensive to implement computationally. To address this, we proposed a new approach called the filtering method to improve the computational efficiency of a permutation test. We first quantified the correlation between QTL and the longitudinal data using a genotype-oriented curve clustering method. Then, the QTLs that are highly correlated with the outcome were computed in the improved permutation tests (Wang et al., 2017). As a result, this reduced the amount of computation in permutation tests significantly and sped up the computation for data analysis in functional mapping.

Package workflow

The Funmap2 package is an open source package for R with automated data analysis for identifying the significant QTLs for the longitudinal traits that were measured. The Funmap2 pipeline includes modules for data import, curve fitting, QTL scanning, MLE computation, hypothesis testing, and data visualization (Fig. 1).

The user may either run the entire Funmap2 workflow in one function call or step-by-step with customization. Often it is more convenient to call the main function FM2.pipe, which automatically implements all tasks and outputs the summary information and figures in PDF format. In this function, sub-modules are called successively. Functions include data loading (FM2.load.data), data estimation (FM2.estimate.data) for curve fitting and covariance selection, QTL scanning (FM2.qtlscan), and permutation (FM.permutation), and report generation (FM2.report). The typical calling is shown below.

Funmap2 requires users to provide experimental data and to specify several parameters. The phenotype file that contains longitudinal traits, one genotype marker file for the experimental population, and one genetic marker information file are required. A covariate file is optional and is not provided in the above example. In addition to these experimental data, users need to specify the cross type of QTL mapping. Four available types are provided in Funmap2: Backcrossing, F2, RIL, and DH. Note that in this automated run, the user may either let Funmap2 choose the optimal curve type and covariance structure or the user may specify their values as functional arguments.

Alternatively, the data analysis can be conducted by customizing the workflow by running sub-modules of Funmap2 successively, as illustrated below.