In charlotte-ngs/lbgfs2020: Livestock Breeding and Genomics (Fall 2020)

knitr::opts_chunk$set(echo = FALSE, results = 'asis', fig.pos = 'H')
knitr::knit_hooks$set(hook_convert_odg = rmdhelp::hook_convert_odg)

cnt <- rmdhelp::R6ClassCount$new()
cnt$set_prefix(ps_prefix = "## Problem")

cat(cnt$out(ps_suffix = "Marker Effect Model"), "\n")

library(dplyr)
### # fix the number of animals
n_nr_animal <- 10
### # intercept
n_inter_cept <- 150
### # residual standard deviation
n_res_sd <- 9.3
### # vector of genotype value coefficients
vec_geno_value_coeff <- c(-1,0,1)
### # sample genotypes of unlinked SNP randomly
set.seed(5432)
### # fix allele frequency of positive allele
n_prob_snps <- .45
### # genotypic values 
vec_geno_val <- c(17.9, 3.3)
n_nr_snp <- length(vec_geno_val)
### # put together the genotypes into a matrix
mat_geno_snp <- matrix(c(sample(vec_geno_value_coeff, n_nr_animal, prob = c((1-n_prob_snps)^2, 
                                                                   2*(1-n_prob_snps)*n_prob_snps, 
                                                                   n_prob_snps^2), 
                       replace = TRUE),
                       sample(vec_geno_value_coeff, n_nr_animal, prob = c(n_prob_snps^2, 
                                                                   2*(1-n_prob_snps)*n_prob_snps, 
                                                                   (1-n_prob_snps)^2), 
                       replace = TRUE)),
                       nrow = n_nr_snp)
mat_obs_y <- n_inter_cept + crossprod(mat_geno_snp, vec_geno_val) + rnorm(n = n_nr_animal, mean = 0, sd = n_res_sd)
### # combine SNP genotypes into a tibble
geno_code <- tibble::tibble(`SNP A` = mat_geno_snp[1,],
                            `SNP B` = mat_geno_snp[2,])

### # add the data
mat_obs_y_rounded <- round(mat_obs_y, digits = 0)
tbl_obs <- tibble::tibble(Observation = mat_obs_y_rounded[,1])
geno_code %>% bind_cols(tbl_obs) -> tbl_all_data
### # add animal ids
tbl_all_data <- bind_cols(Animal = c(1:n_nr_animal),tbl_all_data)

We are given the dataset that is shown in the table below. This dataset contains gentyping results of r n_nr_animal for r n_nr_snp SNP loci.

knitr::kable(tbl_all_data,
             booktabs = TRUE,
             longtable = TRUE,
#             caption = "Animals With Two SNP Loci A and B Affecting A Quantitative Trait",
             escape = FALSE)

Your Task

• The goal of this problem is to estimate SNP marker effects using a marker effect model. Because we have just r n_nr_snp SNP loci, you can use a fixed effects linear model with the r n_nr_snp loci as fixed effects. Furthermore you can also include a fixed intercept into the model.
• Specify all the model components including the vector of observations, the design matrix $X$, the vector of unknowns and the vector of residuals.
• You can use the R-function lm() to get the solutions for estimates of the unknown SNP effects.

cat(cnt$out(ps_suffix = "Breeding Value Model"), "\n")

Use the same data as in Problem 1 to estimate genomic breeding values using a breeding value model.

Hints

• The only fixed effect in this model is the mean $\mu$ which is the same for all observations.
• You can use the following matrix as the genomic relationship matrix

#' Compute genomic relationship matrix based on data matrix
computeMatGrm <- function(pmatData) {
  matData <- pmatData
  # check the coding, if matData is -1, 0, 1 coded, then add 1 to get to 0, 1, 2 coding
  if (min(matData) < 0) matData <- matData + 1
  # Allele frequencies, column vector of P and sum of frequency products
  freq <- apply(matData, 2, mean) / 2
  P <- 2 * (freq - 0.5)
  sumpq <- sum(freq*(1-freq))
  # Changing the coding from (0,1,2) to (-1,0,1) and subtract matrix P
  Z <- matData - 1 - matrix(P, nrow = nrow(matData), 
                             ncol = ncol(matData), 
                             byrow = TRUE)
  # Z%*%Zt is replaced by tcrossprod(Z)
  return(tcrossprod(Z)/(2*sumpq))
}
matG <-computeMatGrm(pmatData = t(mat_geno_snp))
matG_star <- rvcetools::make_pd_rat_ev(matG, pn_max_ratio = 100)

cat(paste(rmdhelp::bmatrix(pmat = round(matG_star, digits = 3), ps_name = 'G', ps_env = '$$'), collapse = '\n'), '\n')

Your Tasks

• Specify all model components of the linear mixed model, including the expected values and the variance-covariance matrix of the random effects.

charlotte-ngs/lbgfs2020 documentation built on Dec. 20, 2020, 5:39 p.m.