Estimate Single Locus Breeding Values from Data

Disclaimer

Breeding values are estimated from data

Read Data

This is the notebook developed in class. It shows how genotype frequencies can be determined using R.

We want to read a dataset with genotype and phenotype information.

s_gen_phen_data <- "https://charlotte-ngs.github.io/lbgfs2022/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

library(dplyr)
tbl_gen_phen %>%
  group_by(Genotype) %>%
  summarise(geno_count = n()) %>%
  summarise(geno_freq = geno_count / sum(geno_count))

NA

Genotypic Values

Regression of phenotypes on Genotypes using homozygous genotypes

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

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

Regression model, with or without intercept? Compute the mean phenotypic value for all animals with Genotype \(0\)

tbl_hom_gen0 <- tbl_hom_gen %>%
  filter(Genotype < 1)
mean(tbl_hom_gen0$Phenotype)

[1] 0.1384745

Plot the data

library(ggplot2)
p_g0 <- ggplot(data = tbl_hom_gen, mapping = aes(x = Genotype, y = Phenotype)) +
  geom_point()
p_g0

From the plot and from the mean of the animals with genotype \(0\), it becomes clear that we need to fit a regression with intercept

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

Call:
lm(formula = Phenotype ~ Genotype, data = tbl_hom_gen)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1200 -0.6895  0.0200  0.7068  3.8270 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.13847    0.17801   0.778    0.437
Genotype     0.04542    0.08966   0.507    0.613

Residual standard error: 1.007 on 2182 degrees of freedom
Multiple R-squared:  0.0001176, Adjusted R-squared:  -0.0003407 
F-statistic: 0.2566 on 1 and 2182 DF,  p-value: 0.6125

Add the regression line to the plot

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

The genotypic value \(a\) is

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

[1] 0.04541911

The genotypic values for heterozygous animals

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

n_geno_d <- mean(tbl_het_gen$Phenotype) - sry_geno_val$coefficients["(Intercept)", "Estimate"] - n_geno_a
n_geno_d

[1] -0.05784937

tbl_het_gen$Genotype <- 1L
head(tbl_het_gen)

tbl_gen_all <- bind_rows(tbl_hom_gen, tbl_het_gen)
head(tbl_gen_all)

p_g_all <- ggplot(data = tbl_gen_all, mapping = aes(x=Genotype,  y=Phenotype))+
  geom_point()+ 
  geom_abline(slope = sry_geno_val$coefficients["Genotype", "Estimate"],
              intercept = sry_geno_val$coefficients["(Intercept)", "Estimate"])
p_g_all

Now that the genotypic values \(a\) and \(d\) are computed, we can do the transformations of the phenotypic values such that the three genotypes have the correct genotypic value

tbl_trans_gen <- tbl_gen_all
tbl_trans_gen$Phenotype <- tbl_trans_gen$Phenotype - sry_geno_val$coefficients["(Intercept)", "Estimate"] - n_geno_a

To check the result of the transformation, we can compute the means of the transformed phenotypic values

tbl_trans_gen %>%
  group_by(Genotype) %>%
  summarise(mean_phen = mean(Phenotype))