inst/doc/F_Quantitative_Genetics.R

In SIMplyBee: 'AlphaSimR' Extension for Simulating Honeybee Populations and Breeding Programmes

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  include = TRUE
)

## ----founderGenomes-----------------------------------------------------------
library(package = "SIMplyBee")
library(package = "ggplot2")
founderGenomes <- quickHaplo(nInd = 20, nChr = 16, segSites = 1000)

## ----SimParamBee_mean_and_varA------------------------------------------------
# Global simulation parameters
SP <- SimParamBee$new(founderGenomes)

nQtlPerChr <- 100

# Genetic parameters for queen and workers effects - each represented by a trait
mean <- c(10, 10 / SP$nWorkers)
varA <- c(1, 1 / SP$nWorkers)

## ----SimParamBee_corA_and_addTrait--------------------------------------------
corA <- matrix(data = c( 1.0, -0.5, 
                        -0.5,  1.0), nrow = 2, byrow = TRUE)
SP$addTraitA(nQtlPerChr = nQtlPerChr, mean = mean, var = varA, corA = corA,
             name = c("queenTrait", "workersTrait"))

## ----SimParamBee_varE_and_corR------------------------------------------------
varE <- c(3, 3 / SP$nWorkers)
corE <- matrix(data = c(1.0, 0.3, 
                        0.3, 1.0), nrow = 2, byrow = TRUE)
SP$setVarE(varE = varE, corE = corE)

## ----basePop_virgin_queens, echo = FALSE, fig.height = 5, fig.width = 6-------
# Base population virgin queens
basePop <- createVirginQueens(founderGenomes, n = 20)
head(basePop@gv)
head(basePop@pheno)
oldpar <- par(mfrow=c(2,2))
limQ <- range(c(basePop@gv[, "queenTrait"], basePop@pheno[, "queenTrait"]))
brkQ <- seq(from = limQ[1], to = limQ[2], length.out = 10)
limW <- range(c(basePop@gv[, "workersTrait"], basePop@pheno[, "workersTrait"]))
brkW <- seq(from = limW[1], to = limW[2], length.out = 10)
hist(basePop@gv[, "queenTrait"], xlab = "Genetic value", main = "Queen effect", xlim = limQ, breaks = brkQ)
hist(basePop@gv[, "workersTrait"], xlab = "Genetic value", main = "Workers effect", xlim = limW, breaks = brkW)
hist(basePop@pheno[, "queenTrait"], xlab = "Phenotypic value", main = "Queen effect", xlim = limQ, breaks = brkQ)
hist(basePop@pheno[, "workersTrait"], xlab = "Phenotypic value", main = "Workers effect", xlim = limW, breaks = brkW)
par(oldpar)

## ----basePop_drones, echo = FALSE, fig.height = 2.7, fig.width = 5------------
# Base population drones
drones <- createDrones(x = basePop[1:5], nInd = 3)
head(drones@gv)
oldpar <- par(mfrow=c(1,2))
hist(drones@gv[, "queenTrait"], xlab = "Genetic value", main = "Queen effect")
hist(drones@gv[, "workersTrait"], xlab = "Genetic value", main = "Workers effect")
par(oldpar)

## ----create_colony------------------------------------------------------------
colony <- createColony(x = basePop[6])
colony <- cross(x = colony, drones = drones, checkCross = "warning")
colony <- addWorkers(x = colony, nInd = 50)
colony

## ----getGv_and_getPheno-------------------------------------------------------
getGv(colony, caste = "queen")
getGv(colony, caste = "workers") |> head(n = 4)

getPheno(colony, caste = "queen")
getPheno(colony, caste = "workers") |> head(n = 4)

## ----getGv_and_getPheno_caste-------------------------------------------------
getQueenGv(colony)
getWorkersGv(colony) |> head(n = 4)

getQueenPheno(colony)
getWorkersPheno(colony) |> head(n = 4)

## ----build_up_colony----------------------------------------------------------
# Check if colony is productive
isProductive(colony)

# Build-up the colony and check the production status again
colony <- buildUp(colony)
colony
isProductive(colony)

## ----data.frame---------------------------------------------------------------
# Collate genetic and phenotypic values of workers
df <- data.frame(id = colony@workers@id,
                 mother = colony@workers@mother,
                 father = colony@workers@father,
                 gvQueenTrait = colony@workers@gv[, "queenTrait"],
                 gvWorkersTrait = colony@workers@gv[, "workersTrait"],
                 pvQueenTrait =  colony@workers@pheno[, "queenTrait"],
                 pvWorkersTrait = colony@workers@pheno[, "workersTrait"])
head(df)

## ----plot_queen_vs_worker_values----------------------------------------------
# Covariation between queen and workers effect genetic values in workers
p <- ggplot(data = df, aes(x = gvQueenTrait, y = gvWorkersTrait)) +
  xlab("Genetic value for the queen effect") +
  ylab("Genetic value for the workers effect") +
  geom_point() +
  theme_classic()
print(p)

## ----fathers_values-----------------------------------------------------------
# Variation in patriline genetic values
getFathersGv(colony)

## ----distribution_by_patriline, echo = FALSE, fig.height = 4.5, fig.width = 6----
# Variation in workers effect genetic values by patriline in workers
p <- ggplot(data = df, aes(x = gvWorkersTrait, colour = father)) +
  xlab("Genetic value for the workers effect") +
  geom_density() +
  theme_classic()
print(p)

## ----colony_pheno-------------------------------------------------------------
# Colony phenotype value
calcColonyPheno(colony, queenTrait = "queenTrait", workersTrait = "workersTrait")
help(calcColonyPheno)
help(mapCasteToColonyPheno)

## ----colony_pheno_change------------------------------------------------------
# Colony phenotype value from a reduced colony
removeWorkers(colony, p = 0.5) |>
  calcColonyPheno(queenTrait = "queenTrait", workersTrait = "workersTrait")

## ----colony_pheno_change1-----------------------------------------------------
# Colony phenotype value from a reduced colony
removeWorkers(colony, p = 0.99) |>
  calcColonyPheno(queenTrait = "queenTrait", workersTrait = "workersTrait")

## ----multicolony--------------------------------------------------------------
apiary <- createMultiColony(basePop[7:20])
drones <- createDrones(basePop[1:5], nInd = 100)
droneGroups <- pullDroneGroupsFromDCA(drones, n = nColonies(apiary), nDrones = 15)
apiary <- cross(x = apiary, drones = droneGroups, checkCross = "warning")
apiary <- buildUp(apiary)

## ----multicolony_gv-----------------------------------------------------------
getQueenGv(apiary) |> head(n = 4)
getQueenGv(apiary, collapse = TRUE) |> head(n = 4)

## ----multicolony_pheno--------------------------------------------------------
colonyGv <- calcColonyGv(apiary)
colonyPheno <- calcColonyPheno(apiary)
data.frame(colonyGv, colonyPheno)

## ----multicolony_selection----------------------------------------------------
# Select the best colony based on gv
selectColonies(apiary, n = 1, by = colonyGv)
# Select the best colony based on phenotype
selectColonies(apiary, n = 1, by = colonyPheno)

## ----SimParamBee_2------------------------------------------------------------
# Global simulation parameters
SP <- SimParamBee$new(founderGenomes)

nQtlPerChr <- 100

# Quantitative genetic parameters - for two traits, each with the queen and workers effects
meanP <- c(10, 10 / SP$nWorkers, 0, 0)
varA <- c(1, 1 / SP$nWorkers, 1, 1 / SP$nWorkers)
corA <- matrix(data = c( 1.0, -0.5,  0.0,  0.0, 
                        -0.5,  1.0,  0.0,  0.0,
                         0.0,  0.0,  1.0, -0.4, 
                         0.0,  0.0, -0.4,  1.0), nrow = 4, byrow = TRUE)
SP$addTraitA(nQtlPerChr = 100, mean = meanP, var = varA, corA = corA,
             name = c("yieldQueenTrait", "yieldWorkersTrait",
                      "calmQueenTrait", "calmWorkersTrait"))

varE <- c(3, 3 / SP$nWorkers, 3, 3 / SP$nWorkers)
corE <- matrix(data = c(1.0, 0.3, 0.0, 0.0,
                        0.3, 1.0, 0.0, 0.0,
                        0.0, 0.0, 1.0, 0.2,
                        0.0, 0.0, 0.2, 1.0), nrow = 4, byrow = TRUE)
SP$setVarE(varE = varE, corE = corE)

## ----base_pop_and_colony------------------------------------------------------
basePop <- createVirginQueens(founderGenomes)
drones <- createDrones(x = basePop[1:5], nInd = 100)
apiary <- createMultiColony(basePop[6:20])
droneGroups <- pullDroneGroupsFromDCA(drones, nColonies(apiary), nDrones = 15)
apiary <- cross(x = apiary, drones = droneGroups, checkCross = "warning")
apiary <- buildUp(apiary)
apiary

## -----------------------------------------------------------------------------
getQueenGv(apiary) |> head(n = 4)
getWorkersPheno(apiary, nInd = 3) |> head(n = 4)

## ----colony_pheno_2a----------------------------------------------------------
colonyValues <- calcColonyPheno(apiary,
                                queenTrait = c("yieldQueenTrait", "calmQueenTrait"),
                                workersTrait = c("yieldWorkersTrait", "calmWorkersTrait"),
                                traitName = c("yield", "calmness"),
                                checkProduction = c(TRUE, FALSE)) |> as.data.frame()
colonyValues

## ----colony_pheno_2b----------------------------------------------------------
myMapCasteToColonyPheno <- function(colony) {
  yield <- mapCasteToColonyPheno(colony,
                                 queenTrait = "yieldQueenTrait",
                                 workersTrait = "yieldWorkersTrait",
                                 traitName = "yield",
                                 checkProduction = TRUE)
  calmness <- mapCasteToColonyPheno(colony,
                                    queenTrait = "calmQueenTrait",
                                    workersTrait = "calmWorkersTrait",
                                    traitName = "calmness",
                                    checkProduction = FALSE)
  return(cbind(yield, calmness))
}
colonyValues <- calcColonyPheno(apiary, FUN = myMapCasteToColonyPheno) |> as.data.frame()
colonyValues

## -----------------------------------------------------------------------------
colonyValues$Index <- selIndex(Y = colonyValues, b = c(0.5, 0.5), scale = TRUE) * 10 + 100
bestColony <- selectColonies(apiary, n = 1, by = colonyValues$Index)
getId(bestColony)

## ----SimParamBee_3------------------------------------------------------------
# Global simulation parameters
SP <- SimParamBee$new(founderGenomes)

# Quantitative genetic parameters
# - the first trait has only the queen effect
# - the second trait has both the queen and workers effects
nWorkers <- 100
mean <- c(nWorkers, 10, 10 / nWorkers)
varA <- c(25, 1, 1 / nWorkers)
corA <- matrix(data = c(1.0,  0.0,  0.0,
                        0.0,  1.0, -0.5, 
                        0.0, -0.5,  1.0), nrow = 3, byrow = TRUE)
SP$addTraitA(nQtlPerChr = 100, mean = mean, var = varA, corA = corA,
             name = c("fecundityQueenTrait", "yieldQueenTrait", "yieldWorkersTrait"))

varE <- c(75, 3, 3 / nWorkers)
corE <- matrix(data = c(1.0, 0.0, 0.0,
                        0.0, 1.0, 0.3,
                        0.0, 0.3, 1.0), nrow = 3, byrow = TRUE)
SP$setVarE(varE = varE, corE = corE)

## ----base_pop_and_colony_2----------------------------------------------------
basePop <- createVirginQueens(founderGenomes)
drones <- createDrones(x = basePop[1:5], nInd = 100)
apiary <- createMultiColony(basePop[6:20])
droneGroups <- pullDroneGroupsFromDCA(drones, nColonies(apiary), nDrones = 15)
apiary <- cross(x = apiary, drones = droneGroups, checkCross = "warning")

## ----queen_values, echo = FALSE, fig.height = 5, fig.width = 6----------------
getQueenGv(apiary, collapse = TRUE)
queenPheno <- getQueenPheno(apiary, collapse = TRUE) |> as.data.frame()
cor(queenPheno)
plot(queenPheno)

## ----colony_strength----------------------------------------------------------
apiary <- buildUp(apiary, nWorkers = nWorkersColonyPhenotype,
                  queenTrait = "fecundityQueenTrait")
cbind(nWorkers = nWorkers(apiary), queenPheno)
help(nWorkersColonyPhenotype)

## ----colony_pheno_3, echo = FALSE, fig.height = 5.5, fig.width = 6.5----------
colonyValuesPheno <- calcColonyPheno(apiary,
                                     queenTrait = "yieldQueenTrait",
                                     workersTrait = "yieldWorkersTrait")
pheno <- cbind(nWorkers = nWorkers(apiary), queenPheno, yield = colonyValuesPheno)
cor(pheno)
plot(pheno)