Dealing with sources of variability in the data-analysis of phenotyping experiments with transgenic rice

Abstract

The analysis of phenotyping experiments for transgenics deserves special attention. Experiments set up for the detection of interesting phenotypes among transgenic plants have to screen several primary events obtained by transforming with a particular transgene, since expression levels of the transgene differ considerably. Agronomically most interesting lines might have an intermediate level of transgene expression. Therefore, attention should be paid to all transformants and how their expression levels differ. Experimental design and the analysis of the data have to focus on the variability among lines and have to be able to detect small differences in quantitative traits. The mixed model is the most adequate approach to analyse data of phenotyping experiments because it reflects the structure and provides the researcher with important measures to allow broader inferences. The paper explains the model and illustrates it using a screening experiment carried out by the high-throughput phenotyping method of TraitMill^TM. Besides inference for a particular experiment and a particular set of lines, the output allows more general predictions for a wider set of non-tested lines. It quantifies the various sources of variability encountered and helps to understand the underlying process. It also helps to optimise the experimental set-up of future experiments. The model presented here has been implemented in the R-language and SAS. The scripts are attached.

Appendix

The variance-covariance matrix of the mixed model

The important features of the mixed model are given in the main text. Some more explanation should be provided on the variance structure. The matrix V is a n ×  n matrix (n is the total number of plants) with a block-withina-block structure. Each larger block represents the variance-covariance between all plants of a line and is composed of 4 subblocks, indicating either the variance-covariance between nulls of a line, the variance-covariance between transgenic plants of a line, or the covariance between transgenics and nulls of a line. The rest of the matrix V is filled by zeros. For instance, the block for the first line (V₁) is given by:

$$ {\mathbf{V}}_1 = \left[ \begin{array}{ll} V_a& V_{ag}\\ V_{ag} &V_g\\ \end{array}\right] $$

where V_a and V_g are square matrices with dimensions equal to the number of nulls (in V_a) or transgenics (in V_g) of line 1, with the structure:

$$ {\mathbf{V}}_{\rm a}=\left[ \begin{array}{lllll} \sigma^2_a+\sigma^2_e & \sigma^2_a & \sigma^2_a& \cdots & \sigma^2_a\\ \sigma^2_a &\sigma^2_a+\sigma^2_e&\sigma^2_a&\cdots&\sigma^2_a\\ \sigma^2_a&\sigma^2_a&\sigma^2_a+\sigma^2_e&\cdots&\sigma^2_a\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \sigma^2_a&\sigma^2_a&\sigma^2_a&\cdots&\sigma^2_a+\sigma^2_e\\ \end{array} \right] $$

$$ {\mathbf{V}}_{\rm g}=\left[ \begin{array}{lllll} \sigma^2_g+\sigma^2_e & \sigma^2_g & \sigma^2_g& \cdots & \sigma^2_g\\ \sigma^2_g &\sigma^2_g+\sigma^2_e&\sigma^2_g&\cdots&\sigma^2_g\\ \sigma^2_g&\sigma^2_g&\sigma^2_g+\sigma^2_e&\cdots&\sigma^2_g\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \sigma^2_g&\sigma^2_g&\sigma^2_g&\cdots&\sigma^2_g+\sigma^2_e\\ \end{array} \right] $$

where as V_ag is rectangular and all elements equal to σ _ag :

$$ {\mathbf{V}}_{\rm ag}=\left[ \begin{array}{ccc} \sigma_{ag}&\cdots&\sigma_{ag}\\ \vdots&&\vdots\\ \sigma_{ag}&\cdots&\sigma_{ag}\\ \end{array}\right] $$

Estimation of the parameters

Estimation of variance components

The estimation of σ ² _e , σ ² _a , σ ² _g and σ _ag cannot be achieved with the usual least squares procedure. Several methods have been proposed for this purpose, but Maximum Likelihood (ML) approaches are commonly used with the Restricted Maximum Likelihood (REML (Patterson and Thompson, 1971)) approach specifically recommended and found as the default method in most statistical software. ML-methods can be used when the form of the probability of the distribution is known (or assumed). Through iterative algorithms, ML chooses as estimates for the parameters those values that are most consistent with the sample data. For mixed models ML yields biased estimates because it does not take in to account fact that fixed effects are in fact estimates when estimating the variance components. REML is similar to ML but takes into account the fixed effects estimates through the use of appropriate linear combinations of observations to obtain unbiased estimates.

Prediction of the random effects

The BLUP of the nulls of line i is given by the expected value of the random effect, conditional on the mean of the observations of the nulls of the i ^th line ($\overline{y}_{i0}.$), but also on the mean of the transgenics of this line ($\overline{y}_{i1}.$):

$$ \hat{a}_i = E(a_i|\overline{y}_{i0.}, \overline{y}_{i1.}) $$

(1)

$$ \hat{g}_i = E(g_i|\overline{y}_{i0.}, \overline{y}_{i1.}) $$

(2)

To solve 1 we need the joint distribution of a _i , $\overline{y}_{i0}.$ and $\overline{y}_{il}.$, which would be fully known if σ ² _a , σ ² _g , σ _ag and σ ² _e would be known. In the absence of the true values for these parameters, the REML estimates are used. Furthermore, â _i and ĝ _i depend also the number of observation, respectively n _io and n _il . The BLUP of real interest in these experiments is the effect of the individual insertions or $\widehat{g_i - a_i}$, which is obtained by subtracting (1) from (2). The prediction of these effects will differ from their estimates in the fixed model, and can also be shrunk or expanded. Examples of these BLUPs and their shrinkage or expansion are given in the main text.

Estimation of fixed effects

The fixed effects in a mixed model are estimated through a generalised least squares estimation much like in fixed ANOVA. This method needs to know the variance components. Since these components are not known precisely, but estimated by REML, the GLS estimator of the fixed effects is biased. A small bias is acceptable since the estimator is more precise than the other candidate estimators (Greenland, 2000).

Simultaneous estimation of fixed and random effects

A method for getting BLUPs and fixed effects simultaneously is through the solution of Henderson’s mixed model equations, given by

$$ \left[\begin{array}{ll} X^{T}R^{-1}X& X^{T}R^{-1}X\\ Z^{T}R^{-1}X& Z^{T}R^{-1}X + G^{-1} \end{array} \right]\left[\begin{array}{l} \hat{\beta}\\ \hat{u}\end{array}\right]=\left[\begin{array}{l}X^{T}R^{-1}Y\\ Z^{T}R^{-1}Y\end{array}\right] $$

(3)

where Y,X and Z are the matrices from equation (5) from the main text. R is the dispersion matrix of the residuals, i.e. σ ² _e I _k . G is the dispersion matrix of the random effects, defined by matrix (9) from the main text. The advantage of this formulation is that the calculations involve the inverting of matrices is far simpler than the inversion of V, which is necessary in the generalised least squares solution.

From equations (3) it can be derived that for the present model the BLUP for the nulls of the first line is:

$$ \begin{array}{ll} \hat{a}_i=&(\sum r_{null}(\frac{\sigma_a^2}{D}+\frac{n}{\sigma_e^2})+ \sum r_{tr} \frac{\sigma_{ag}^2}{D})\\ &*\frac{\sigma_p^2 D^2}{(\sigma_a^2\sigma_e^2+nD)(\sigma_g^2\sigma_e^2+nD)-\sigma_{ag}^2} \end{array} $$

(4)

where ∑ r _null is the sum of all deviations of the null plants in line 1 from the fixed effect, ∑ r _tr the sum of all deviations of the transgenic plants in line 1 from the fixed effect, n the number of null plants in line 1, and

$$ D = \sigma_a^2\sigma_g^2-\sigma_{ag}^2 $$

The BLUP for the transgenics of the first line is similar. It is clear from equation 4 that the BLUP for the nulls does not only involves the observations for the nulls, but also the observations of the transgenics, the variance components and the number of observations. The smaller the covariance σ _ag , the less important are the observations on the transgenics. The BLUP also tends to be more close to the fixed estimate if the number of observations n is large and the residual variance σ ² _e is small. D will tend to zero if the correlation ρ _ag tends to 1. High correlations will cause a strong shrinkage.

Estimation of linear combinations

The mixed model equations (3) can also be used to estimate linear combinations of fixed effects and/or BLUPs. For instance, the BLUP of a particular linear combination is defined as:

$$ L^{T}\widehat{\beta+M^T} u=L^{T}\hat{\beta}+M^{T}\hat{u} $$

(5)

The variance of this linear combination is given by:

$$ V ar(L^T\hat{\beta}+M^{T}\hat{u}) =[L^T\quad M^T] [C^{-1}]\left[\begin{array}{l} L\\M\end{array}\right] $$

(6)

where C is short for the leftmost matrix of equations (3). All predictable functions can be constructed this way. A detailed description of prediction with mixed models is given in reference (Welham et al, 2004). An important concept in this context is the prediction space, i.e. the population to which the inference applies. We will elaborate further on this in the next section.

Inference for mixed models

After we have fit a model, we normally want to assess the significance and the precision of the estimates. In these experiments there are only two test that are of interest. First there is the test of whether the fixed effect, i.e. the global effect of the transgene, is significant. Second, we are interested in comparing the BLUPs for nulls and transgenics within a line.

Testing a fixed effect

The global difference between nulls and transgenics over all lines is a fixed effect. Fixed effects or linear combinations of fixed effects in mixed models could theoretically be tested by a Wald test. The Wald statistic is given by

$$ A=\frac{\hat\beta_i^2}{\hat{Va}r(\hat\beta_i)} $$

(7)

Under the null hypothesis H ₀:β _i  = 0 against H ₀: β _i ≠  0, A has asymptotically a χ² -distribution with 1 degree of freedom. However, since the variance is estimated based on small samples, the Wald test is inaccurate, having high Type I error rates. Solutions include (Bartlett adjusted) likelihood ratio tests (Welham and Thompson, 1997), but in most cases the Wald statistic is still used but evaluated against an F _ν1,ν2-distribution. In this empirical method the numerator degrees of freedom ν2 is equal to the number of contrasts tested. For the denominator degrees of freedom ν1 several possible solutions exist. Recommended are the Satterthwaite method and the naive method (Manor and Zucker, 2004). For the experiment described above both methods will reduce to twice the number of transgenic lines tested (2q).

Testing a linear combination of fixed and random effects

Linear combinations can be tested through a t-test or by constructing confidence intervals based on the standard errors obtained from the mixed model equations (6). When calculating standard errors in the context of mixed models, the inference space has to be defined. The prediction space is narrow when the prediction applies to the insertions tested in the experiment. In other words, the prediction for the narrow space is resembling estimates in a fixed effects model. On the other hand, the prediction space is broad when the prediction applies to the entire population of insertions. If some random effects are kept fixed and others are random, the prediction space is intermediate. The prediction space can be set by the matrix $\user2{M}$ in the equation (5) and (6). A prediction in the broad space will have all elements of $\user2{M}$ equal to zero. In this way, the random effect does not contribute to (5) but it does add a variance component to (6). In the narrow space, all or some elements of $\user2{M}$ are non-zero. The point prediction is influenced by non-zero elements but they add no variance. The contrasts between the transgenics and the control plants of a particular event involves both fixed and random effects, but for the random effect a particular level is chosen. In other words, these contrast apply necessarily to the narrow prediction space.

Computer programme

All analyses have been done using R version 2.4.0 and the library nlme. The main analysis statement is

$$ \begin{array}{l} {\tt mod < - lme(data = DS,} \\ {\tt fixed = TKW\ ^{\sim}\ Transgenity,}\\ {\tt random = list(Line = pdLogChol({}^{\sim}Transgenity - 1)))} \end{array} $$

With the more recent R library lme4 the same model is obtained by

$$ \begin{array}{l} {\tt mod\ < -\ lmer(data = DS,}\\ {\tt TKW\ ^{\sim}\ Transgenity + (Transgenity-1|Line)}) \end{array} $$

and most presented data can be extracted from the mod object through summary(mod), fixef(mod), ranef(mod) and n getVarCov(mod). For the calculation of standard errors and p-values on differences between BLUPs, a function had to be written based on the Henderson formulation. This function is available from the authors.

The same model can be obtained with SAS proc mixed with the programme:

$$ \begin{array}{l} {\tt TransgenityR = Transgenity}\\ {\tt proc\ mixed\ data=DS \quad method=reml \quad ;}\\ {\tt class\ Transgenity\ Line;}\\ {\tt model\ TKW=Transgenity\ /\ solution}\\ {\tt ddfm=satterth \quad outp=outpredict \quad ;}\\ {\tt random\ TransgenityR\ /\ subject=Line}\\ {\tt type=UN \ solution\ g \quad ;}\\ {\tt run;} \end{array} $$