Overview

We will learn how to connect phenotypes and genomics using genomic best linear unbiased prediction (GBLUP).

Load R objects

Use the function read_excel in the readxl package to read the pedigree file (Simdata.xlsx) in a data frame format. Use the function load() to load the SNP matrix (W.Rda), the numerator relationship matrix (A.Rda) and the results from pedigree-based BLUP (day17.Rda) we created in class.

library(readxl)
dat <- read_excel(file.choose())
dim(dat)
head(dat)


load(file.choose())  # A.Rda
dim(A)

load(file.choose())  # W.Rda
dim(W)

load(file.choose())  # day17.Rda

Statistical model

The statistical model we will fit is given by \[ \mathbf{y} = \mathbf{X}\mathbf{b} + \mathbf{Z}\mathbf{u} + \mathbf{e}, \] where \(\mathbf{y}\) is the vector of birth weights, \(\mathbf{X}\) and \(\mathbf{Z}\) are incident matrices for fixed and random effects, respectively, \(\mathbf{b}\) is the vector of regression coefficients for fixed effects, \(\mathbf{u}\) is the vector of regression coefficients for random genetic values, and \(\mathbf{e}\) is the vector of residuals. Recall that BLUP of \(\mathbf{u}\) is given by Henderson (1975) (doi: 10.2307/2529430) \[ \begin{align*} BLUP(\mathbf{u}) &= E(\mathbf{u} | \mathbf{y}) \\ &= Cov(\mathbf{u}, \mathbf{y}') Var(\mathbf{y})^{-1} [\mathbf{y} - E(\mathbf{y})] \\ &= \mathbf{G}\sigma^2_{G} \mathbf{Z}'\mathbf{V}^{-1}(\mathbf{y} - \mathbf{X}\hat{\mathbf{b}}), \end{align*} \] where \[ \begin{align*} \mathbf{V} &= Var(\mathbf{y}) \\ &= \mathbf{Z}\mathbf{G}\sigma^2_G\mathbf{Z}' + \mathbf{I}\sigma^2_e \end{align*} \]

We predict genomic estimated breeding values (GEBV) or genetic values \(\hat{\mathbf{u}}\) in the following two steps.

• fit OLS to estimate fixed effects (\(\hat{\mathbf{b}}\))
• fit BLUP to obtain GEBV (\(\hat{\mathbf{u}}\)) condition on the estimated \(\hat{\mathbf{b}}\)

Later we will discuss how to obtain \(\hat{\mathbf{b}}\) and \(\hat{\mathbf{u}}\) simultaneously.

Computing a genomic relationship matrix (G)

We will compute the \(\mathbf{G}\) matrix defined by VanRaden (2008). The function computeG() accepts two arguments: 1) a genotype matrix and 2) a cutoff point for minor allele frequency (MAF).

computeG <- function(W, maf) {
    p <- colMeans(W)/2
    maf2 <- pmin(p, 1 - p)
    index <- which(maf2 < maf)
    W2 <- W[, -index]
    p2 <- p[-index]
    Wc <- scale(W2, center = TRUE, scale = FALSE)
    G <- tcrossprod(Wc)/(2 * sum(p2 * (1 - p2)))
    return(G)
}

G <- computeG(W, maf = 0.05)
dim(G)

Incidence matrix X

We will now contruct each component one by one. First we will create the incidence matrix \(\mathbf{X}\). We first subset the variable dat by age of dams and sex of calves, and create a new variable dat2. The model.matrix() function returns the incidence matrix \(X\).

dat2 <- dat[, c("AgeDam", "Sex")]
head(dat2)
X <- model.matrix(~dat2$AgeDam + dat2$Sex)
dim(X)
head(X)
tail(X)

Incidence matrix Z

Next we will create the incidence matrix \(\mathbf{Z}\).

Z <- diag(nrow(G))
dim(Z)
diag(Z)

Incidence matrix I

The third incidence matrix we create is \(\mathbf{I}\).

I <- diag(nrow(G))
dim(I)
diag(I)

Variance components estimation using the regress package

The regress() function in the regress package fits a liner mixed model and estimate variance components (\(\sigma^2_G\) and \(\sigma^2_e\)) through a restricted maximum likelihood (REML). The variance-covariance structure of \(\mathbf{u}\) is given by the genomic relationship matrix \(\mathbf{G}\).

install.packages("regress")
library(regress)
`?`(regress)

y <- dat$BW_205d
varcomp <- regress(y ~ -1 + X, ~G)
summary(varcomp)
varcomp$sigma

sigma2G <- varcomp$sigma[1]  # additive genetic variance
sigma2G
sigma2e <- varcomp$sigma[2]  # residual variance 
sigma2e

Inverse of V

Here we compute the inverse of \(V\) matrix. We will use 1) the multiplication operator %*% and 2) the function solve() to obtain the inverse of matrix.

V <- Z %*% G %*% t(Z) * sigma2G + I * sigma2e
dim(V)
Vinv <- solve(V)
dim(Vinv)

Ordinary least squares

We will use the lm() function to fit OLS and estimate fixed effect coefficients. This function offers an intuitive formula syntax. The summary() function summarizes a model fit.

fit <- lm(BW_205d ~ AgeDam + factor(Sex), data = dat)
summary(fit)
residuals(fit)

GBLUP of genetic values

Now we compute GEBV (\(\hat{\mathbf{u}}\)).

uhatG <- sigma2G * G %*% t(Z) %*% Vinv %*% (matrix(y) - matrix(predict(fit)))
uhatG2 <- sigma2G * G %*% t(Z) %*% Vinv %*% (matrix(y) - matrix(X %*% fit$coefficients))  # alternative
table(uhatG == uhatG2)
head(uhatG)
tail(uhatG)

Let’s plot a histogram of GEBV.

hist(uhatG, col = "lightblue", main = "Histogram", xlab = "Genomic estimated breeding values")

Save as R objects

save(uhatG, h2G, sigma2G, file = "day18.Rda")
save(G, file = "G.Rda")