Single Locus Breeding Values

Disclaimer

The dataset with one locus is used to compute breeding values based on one locus. This requires to first estimate the genotypic values.

Read the dataset

The dataset is read, with an explict specification of column datatypes.

s_gen_phen_data <- ifelse(b_online, 
                          "https://charlotte-ngs.github.io/lbgfs2022/data/p1_mrk_one_locus.csv",
                          file.path(here::here(), "docs", "data", "p1_mrk_one_locus.csv"))
library(readr)
tbl_gen_phen <- read_delim(file = s_gen_phen_data,
                           delim = ";",
                           col_types = cols(
                             ID = col_integer(),
                             Genotype = col_integer(),
                             Phenotype = col_double()
                           ))
head(tbl_gen_phen)

Descriptive Statistics for Genotypes

What are the available genotypes

library(dplyr)
tbl_gen_phen %>% 
  group_by(Genotype) %>%
  summarize(geno_count = n())

Genotypic Values

The genotypic values are estimated with a regression through the homozygous genotypes. For that, we create a new tibble with only the homozygous genotypes

tbl_homozygous_gen <- tbl_gen_phen %>%
  filter(Genotype < 3)
head(tbl_homozygous_gen)

Checking with the counts

tbl_homozygous_gen %>%
  group_by(Genotype) %>%
  summarise(geno_count = n())

Whether, we need a regression with or without intercept can be determined by the average value of all animals with Genotype \(0\).

Because, the mean of all animals with genotype \(0\) is different from \(0\), we use a regression with intercept.

The regression is modelled via lm()

lm_geno_val <- lm(Phenotype ~ Genotype, data = tbl_homozygous_gen)
(sry_geno_val <- summary(lm_geno_val))

The genotypic value \(a\) for the homozygous animals is given by the estimate of the regression coefficient and is

(geno_val_a <- sry_geno_val$coefficients["Genotype", "Estimate"])

[1] 0.04541911

Assigning genotypic values to the homozygous genotypes requires to subtract the intercept and the genotypic value \(a\) from the current mean phenotypes for the two genoypes

geno_val_hom0 <- tbl_hom_gen_mean$mean_phen[tbl_hom_gen_mean$Genotype == 0] - 
  sry_geno_val$coefficients["(Intercept)", "Estimate"] -
  geno_val_a
geno_val_hom0

[1] -0.04541911

For genotype \(2\)

geno_val_hom2 <- tbl_hom_gen_mean$mean_phen[tbl_hom_gen_mean$Genotype == 2] - 
  sry_geno_val$coefficients["(Intercept)", "Estimate"] -
  geno_val_a
geno_val_hom2

[1] 0.04541911

The genotypic value \(d\) is determined with the heterozygous animals. First we select them into a separate tibble

tbl_heterozygous_gen <- tbl_gen_phen %>%
  filter(Genotype > 2)
head(tbl_heterozygous_gen)

The genotype of the heterozygous animals is re-assiged to be \(1\) and no longer \(3\) or \(4\)

tbl_heterozygous_gen$Genotype <- 1L
head(tbl_heterozygous_gen)

The mean of the heterozygous animals is

n_mean_phen_het <- mean(tbl_heterozygous_gen$Phenotype)
n_mean_phen_het

The meaan of the heterozygous animals must also corrected for the intercept and for the genotypic value \(a\).

Hence the value for \(d\) is

(gen_val_d <- n_mean_phen_het - - 
  sry_geno_val$coefficients["(Intercept)", "Estimate"] -
  geno_val_a)

[1] 0.2190996

Plots

Try some plots which might be instructive

library(ggplot2)
p_hom <- ggplot(data = tbl_homozygous_gen, aes(x = Genotype, y = Phenotype)) + 
  geom_point()
p_hom

All genotypes

tbl_geno_recode <- bind_rows(tbl_homozygous_gen, tbl_heterozygous_gen)
p_all <- ggplot(data = tbl_geno_recode, mapping = aes(x = Genotype, y = Phenotype)) +
  geom_point()
p_all

Adding the regression line to the plot

p_all_reg <- p_all + 
  geom_abline(intercept = sry_geno_val$coefficients["(Intercept)", "Estimate"], 
              slope = sry_geno_val$coefficients["Genotype", "Estimate"])
p_all_reg