The purpose of this script is to predict genomic breeding values (GEBVs) using a random regression model. The RR model is \[ PSA_{tjk} = \sum_{k=0}^{2}\phi(t)_{jk}\beta_k + \sum_{k=0}^{2}\phi(t)_{jk} u_{jk} + \sum_{k=0}^{1}\phi(t)_{jk} s_{jk} + e_{tjk}\]. Where \(\beta\) is the fixed second-order Legendre polynomial to model the overall trend in the trait overtime, \(u_{jk}\) and \(s_{jk}\) are the \(k^{th}\) random regression coefficients for additive genetic effect and random experiment of line \(j\), and \(e_{tjk}\) is the random residual. This model was seelcted through comparisons between several other RR models. The complete process is decribed in our previous paper in Plant Direct and on bioRxiv.

Since asreml-R does not allow the use of Legendre polynomials all analysis were done with the standalone version of ASREML. Below the “.as” file is provided as well as the call from the commandline.

Fitting the RR model.

Here is the .as file.

RR model selection
 Exp !A
 DOI 20 !I
 PSA !/100000
G2.grm !SKIP 1 !MAXITER 1000 !EXTRA 100 !WORKSPACE 6144 !ASUV !SLOW 10 !DOPART 1

1.37E-03 1.11E-03 1.87E-03 3.00E-03 4.81E-03 7.08E-03 1.04E-02 1.33E-02 1.78E-02 2.40E-02 3.52E-02 5.10E-02 7.23E-02 0.101765 0.1486 0.203647 0.308825 0.484353 0.743772 1.11587
0.456320 0.360961
0.117960 0.945203E-01 0.251287E-01
0.400317E-01 0.267477E-01
PSA ~ leg(DOI,2) !r us(leg(DOI,1) $USp).Exp us(leg(DOI,2) $USg).grm(NID) !f mv
residual id(1071).idh(DOI $USe)

Here is how it was run.

## bash: line 0: cd: RR_GP/: No such file or directory
##  ASReml 4.1 [28 Dec 2014]  mr [31 May 2017]  22 Jan 2019 14:17:12
## Registered to: University of Nebraska Lincoln
## Serial Number: 40216022,   expiry: 28-feb-2019         
##   Unable to open file                                                     

Solving for GEBVs at each time point

GEBVs at each time point were be obtained following to Mrode (2014). For line \(j\) at time \(t\), the GEBVs can be obtained by \(\text{gBLUP}_{jt} = \phi_t\hat{u}_j\); where \(\phi_t\) is the row vector of the matrix of Legendre polynomials of order 2.

The functions below were adapted from Mrode (2005) by Gota Morota

Here, we’ll take the sln files, extract the breeding values for the RR coefficients, and calcualte the gBLUPs at each time point.