Quantitative genetics
2023-12-22
Source:vignettes/F_Quantitative_Genetics.Rmd
F_Quantitative_Genetics.RmdIntroduction
This vignette describes and demonstrates how SIMplyBee implements quantitative genetics principles for honeybees. Specifically, it describes three different examples where we simulate:
Honey yield - a single colony trait,
Honey yield and Calmness - two colony traits, and
Colony strength and Honey yield - two colony traits where one trait impacts the other one via the number of workers.
We start by loading SIMplyBee and quickly simulating genomes for some founder honeybees. Specifically, we will simulate genomes for 20 individuals with 16 chromosomes and 1000 segregating sites per chromosome.
library(package = "SIMplyBee")
#> Loading required package: AlphaSimR
#> Loading required package: R6
#>
#> Attaching package: 'SIMplyBee'
#> The following object is masked from 'package:base':
#>
#> split
library(package = "ggplot2")
founderGenomes <- quickHaplo(nInd = 20, nChr = 16, segSites = 1000)Honey yield
This section shows how to simulate one colony trait, honey yield, that is influenced by the queen and workers as well as the environment. We will achieve this by:
- setting base population quantitative genetic parameters,
- inspecting individual values in the base population,
- inspecting individual values in a colony,
- calculating colony value,
- calculating multi-colony values, and
- selecting on colony values.
Base population quantitative genetic parameters
AlphaSimR, and hence SIMplyBee, simulates each individual with its
corresponding genome, and quantitative genetic and phenotypic values. To
enable this simulation, we must set base population quantitative genetic
parameters for the traits of interest in the global simulation
parameters via SimParamBee. We must set:
- the number of traits,
- the number of quantitative trait loci (QTL) that affect the traits,
- the distribution of QTL effects,
- trait means, and
- trait genetic and environmental variances - if we simulate multiple traits, we must also specify genetic and environmental covariances between the traits.
In honeybees, the majority of traits are influenced by the queen and workers. There are many biological mechanisms for these queen and workers effects. Depending on which caste is the main driver of the trait (the queen or workers), we also talk about direct and indirect effects. For example, for honey yield, workers directly affect honey yield by foraging, while the queen indirectly affects honey yield by stimulating workers via pheromone production. The queen and workers effects for a trait can be genetically and environmentally independent or correlated (usually negatively).
Here, we will simulate two traits to represent the queen and workers
effects on honey yield. From this point onward we will use the terms the
queen effect and queen trait interchangeably. The same applies to
workers effect and workers trait. These two effects (=traits) will give
rise to honey yield trait. We will assume that colony honey yield is
approximately normally distributed with the mean of 20 kg and variance
of 4 \(kg^2\), which implies that most
colonies will have honey yield between 14 kg and 26 kg (see
hist(rnorm(n = 1000, mean = 20, sd = sqrt(4)))). Traits
like honey yield have a complex polygenic genetic architecture, so we
will assume that this trait is influenced by 100 QTL per chromosome
(with 16 chromosomes, this gives us 1600 QTL in total).
We will first initiate global simulation parameters and set the mean of queen effects to 10 kg with genetic variance of 1 \(kg^2\), while we will set the mean of workers effects to 10 kg with genetic variance of 1 \(kg^2\). The mean and the variance for the worker effect are proportionally scaled by the expected number of workers in a colony. The mean and variance for the queen effect is assumed larger than for the workers effect, because there is one queen and many workers in colony and we assume that workers effects “accumulate”. Deciding how to split the colony mean between queen and workers effects will depend on the individual to colony mapping function, which we will describe in the Colony value sub-section.
# 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)We next set genetic correlation between the queen and workers effects
to -0.5 to reflect the commonly observed antagonistic relationship
between these effects. With all the quantitative genetic parameters
defined, we now add two additive traits to global simulation parameters
and name them queenTrait and workerTrait.
These parameters drive the simulation of QTL effects. Read about all the
other trait simulation options in AlphaSimR via:
vignette(topic = "traits", package="AlphaSimR").
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"))Finally, we set the environmental variance of the queen and workers effects to 3 \(kg^2\) and we again scale the worker variance by the expected number of workers. Contrary to the negative genetic correlation, we here assume that environmental correlation between the queen and workers effects is slightly positive, 0.3. This is just an example! These parameters should be based on literature or simulation scenarios of interest.
Individual values in the base population
Now we create a base population of virgin queens. Since we defined
two traits, all honeybees in the simulation will have genetic and
phenotypic values for both traits. The genetic values are stored in the
gv slot of each Pop object, while phenotypic
values are stored in the pheno slot.
#> queenTrait workersTrait
#> [1,] 10.454034 0.07548599
#> [2,] 9.528805 0.20444512
#> [3,] 9.839773 0.08518108
#> [4,] 10.389981 0.22269740
#> [5,] 10.505698 0.03370364
#> [6,] 10.642709 -0.02208124
#> queenTrait workersTrait
#> [1,] 10.103301 0.1121241
#> [2,] 5.344579 0.3129457
#> [3,] 10.355808 0.3726329
#> [4,] 10.550549 0.3413543
#> [5,] 11.656893 0.4479464
#> [6,] 12.847070 0.0422259
Note that these are virgin queens, yet we obtained queen and workers effect values for them! Is this wrong? No! Virgin queens carry DNA with genes that are differentially expressed in different castes, which would be only showed in their phenotype. Hence, virgin queens have genetic values for the queen and worker effects, but they might never actually express these effects. In this simulation virgin queens also obtained phenotypic values for both of the effects. This is technically incorrect because virgin queens don’t express genes for the worker effect at all, and they also do not express the queen effect, not until they become the queen of a colony. We can treat these phenotypic values for virgin queens as values that we could see if these virgin queens would express these traits. We will show later in the Colony value sub-section how we use these traits from different castes. If existence of these phenotypic values for certain castes is a hindrance, we can always remove them for population or colony objects by modifying the corresponding slots as required.
As with the virgin queens, drones also carry DNA with genes that are expressed in different castes. Therefore, drones will also have the queen and workers effect genetic (and phenotypic values) for honey yield even though they do not contribute to this trait in a colony.
#> queenTrait workersTrait
#> [1,] 10.467976 0.10551154
#> [2,] 11.002152 0.01920906
#> [3,] 10.171420 0.02618648
#> [4,] 8.024650 0.32761118
#> [5,] 7.552598 0.20489027
#> [6,] 9.973988 0.14110013
Individual values in a colony
We continue by creating a colony from one base population virgin queen, crossing it, and adding some workers.
colony <- createColony(x = basePop[6])
colony <- cross(x = colony, drones = drones, checkCross = "warning")
colony <- addWorkers(x = colony, nInd = 50)
colony
#> An object of class "Colony"
#> Id: 1
#> Location: 0 0
#> Queen: 6
#> Number of fathers: 15
#> Number of workers: 50
#> Number of drones: 0
#> Number of virgin queens: 0
#> Has split: FALSE
#> Has swarmed: FALSE
#> Has superseded: FALSE
#> Has collapsed: FALSE
#> Is productive: FALSEWe can access the genetic and phenotypic values of colony members
with functions getGv() and getPheno(), both of
which have the caste argument (see more via
help(getGv)).
getGv(colony, caste = "queen")
#> queenTrait workersTrait
#> 6 10.64271 -0.02208124
getGv(colony, caste = "workers") |> head(n = 4)
#> queenTrait workersTrait
#> 36 9.234616 0.18025244
#> 37 9.003230 0.09253210
#> 38 10.235667 -0.03916070
#> 39 10.992221 -0.02080473
getPheno(colony, caste = "queen")
#> queenTrait workersTrait
#> 6 12.84707 0.0422259
getPheno(colony, caste = "workers") |> head(n = 4)
#> queenTrait workersTrait
#> 36 6.525535 0.22300672
#> 37 7.231879 -0.03854496
#> 38 10.965199 0.04827966
#> 39 12.027848 0.19901663For convenience, there are also alias functions for accessing the genetic and phenotypic values of each caste directly.
getQueenGv(colony)
#> queenTrait workersTrait
#> 6 10.64271 -0.02208124
getWorkersGv(colony) |> head(n = 4)
#> queenTrait workersTrait
#> 36 9.234616 0.18025244
#> 37 9.003230 0.09253210
#> 38 10.235667 -0.03916070
#> 39 10.992221 -0.02080473
getQueenPheno(colony)
#> queenTrait workersTrait
#> 6 12.84707 0.0422259
getWorkersPheno(colony) |> head(n = 4)
#> queenTrait workersTrait
#> 36 6.525535 0.22300672
#> 37 7.231879 -0.03854496
#> 38 10.965199 0.04827966
#> 39 12.027848 0.19901663Some phenotypes, such as honey yield, are only expressed if colony is
at full size. This is achieved by the buildUp() colony
event function that adds worker and drones and hence turns on the
production status of the colony (to TRUE).
SIMplyBee includes a function ìsProductive() to check the
production status of a colony.
# Check if colony is productive
isProductive(colony)
#> [1] FALSE
# Build-up the colony and check the production status again
colony <- buildUp(colony)
colony
#> An object of class "Colony"
#> Id: 1
#> Location: 0 0
#> Queen: 6
#> Number of fathers: 15
#> Number of workers: 100
#> Number of drones: 100
#> Number of virgin queens: 0
#> Has split: FALSE
#> Has swarmed: FALSE
#> Has superseded: FALSE
#> Has collapsed: FALSE
#> Is productive: TRUE
isProductive(colony)
#> [1] TRUEFor the ease of further demonstration, we now combine workers’ values into a single 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)
#> id mother father gvQueenTrait gvWorkersTrait pvQueenTrait pvWorkersTrait
#> 1 86 6 28 10.346503 0.043455224 8.442814 0.12631675
#> 2 87 6 34 10.095304 0.001673821 12.649532 -0.04915628
#> 3 88 6 27 9.249778 0.112582559 11.563785 0.18416347
#> 4 89 6 21 10.746130 0.061059150 11.896217 0.24176906
#> 5 90 6 26 10.372937 -0.011849426 7.787357 -0.27137811
#> 6 91 6 24 9.254217 0.127753126 7.588014 -0.06526353To visualise correlation between queen and workers effects in workers, we plot these effect values against each other.
# 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)
In SIMplyBee, we know genetic values of all individuals, including drones that the queen mated with (=fathers in a colony)!
# Variation in patriline genetic values
getFathersGv(colony)
#> queenTrait workersTrait
#> 21 10.467976 0.105511542
#> 22 11.002152 0.019209059
#> 23 10.171420 0.026186477
#> 24 8.024650 0.327611177
#> 25 7.552598 0.204890271
#> 26 9.973988 0.141100130
#> 27 9.012092 0.154409907
#> 28 10.506579 0.012161614
#> 29 9.330237 0.114271654
#> 30 11.268878 0.155746156
#> 31 8.800179 0.284058746
#> 32 10.975785 0.132860958
#> 33 12.386634 -0.180618064
#> 34 9.882475 -0.013788111
#> 35 9.986798 -0.007105668Knowing the father of each worker, we inspect variation in the distribution of genetic values of worker by the patriline (workers from a single father drone) for the workers effect.
Colony value
However, in honeybees we usually don’t observe values on individuals,
but on a colony. SIMplyBee provides functions for mapping individual
values to a colony value. The general function for this is
calcColonyValue(), which can combine any value and trait
from any caste. There are also aliases calcColonyGv() and
calcColonyPheno(). These functions require users to specify
the so-called mapping function (via the FUN argument). The
mapping function specifies queen and workers traits (potentially also
drone traits) and what function we want to apply to each of them before
mapping them to the colony value(s). We can also specify whether the
colony value(s) depend on the production status. For example, if a
colony is not productive, its honey yield would be 0 or unobserved.
SIMplyBee provides a general mapping function
mapCasteToColonyValue() and aliases
mapCasteToColonyGv() and
mapCasteToColonyPheno(). These functions have arguments to
cater for various situations. By default, they first calculate caste
values: leave the queen’s value as it is, sum workers’ values,
potentially sum drones’ values, and lastly sum all these caste values
together into a colony value. Users can provide their own mapping
function(s) too!
We now calculate honey yield for our colony - a single value for the colony.
# Colony phenotype value
calcColonyPheno(colony, queenTrait = "queenTrait", workersTrait = "workersTrait")
#> [,1]
#> [1,] 17.57193
help(calcColonyPheno)
help(mapCasteToColonyPheno)These colony values are not stored in a colony, because they change as colony changes due to various events. For example, reducing the number of workers will reduce the colony honey yield.
# Colony phenotype value from a reduced colony
removeWorkers(colony, p = 0.5) |>
calcColonyPheno(queenTrait = "queenTrait", workersTrait = "workersTrait")
#> [,1]
#> [1,] 15.09647Please note that we assumed that the queen contributes half to colony honey yield and workers contribute the other half. This means that removing workers will still give a non-zero honey yield! This shows that we have to design the mapping between individual, caste, and colony values with care!
# Colony phenotype value from a reduced colony
removeWorkers(colony, p = 0.99) |>
calcColonyPheno(queenTrait = "queenTrait", workersTrait = "workersTrait")
#> [,1]
#> [1,] 12.81607Finally, note that SIMplyBee currently does not provide functionality
for breeding values, dominance deviations, and epistatic deviations at
caste and colony levels, despite the availabiliy of AlphaSimR
bv(), dd(), and aa() functions.
This is because we have to check or develop theory on how to calculate
these values across active colonies and hence we currently advise
against the use of AlphaSimR bv(), dd(), and
aa() functions with SIMplyBee as the output of these
functions could be easily misinterpreted.
MultiColony values
The same functions can be used on a MultiColony class
object. Let’s create an apiary.
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)We can extract the genetic and phenotypic values from multiple
colonies in the same manner as from a single colony, by using
get*Gv() and get*Pheno() functions. The output
of these function is a named list with values for each colony or a
single matrix if we set the collapse argument to
TRUE.
getQueenGv(apiary) |> head(n = 4)
#> $`2`
#> queenTrait workersTrait
#> 7 9.448344 0.2707414
#>
#> $`3`
#> queenTrait workersTrait
#> 8 10.70298 0.05659256
#>
#> $`4`
#> queenTrait workersTrait
#> 9 7.586921 0.1220963
#>
#> $`5`
#> queenTrait workersTrait
#> 10 8.422281 0.1070055
getQueenGv(apiary, collapse = TRUE) |> head(n = 4)
#> queenTrait workersTrait
#> 7 9.448344 0.27074140
#> 8 10.702985 0.05659256
#> 9 7.586921 0.12209631
#> 10 8.422281 0.10700554In a similar manner, we can calculate colony value for all the colonies in our apiary, where the row names of the output represent colony IDs.
colonyGv <- calcColonyGv(apiary)
colonyPheno <- calcColonyPheno(apiary)
data.frame(colonyGv, colonyPheno)
#> colonyGv colonyPheno
#> 2 28.77589 27.360001
#> 3 20.20733 23.000192
#> 4 16.76474 13.166057
#> 5 20.95291 20.315861
#> 6 11.82838 6.928691
#> 7 26.17827 29.480943
#> 8 21.87154 19.032215
#> 9 19.99128 18.581427
#> 10 21.08818 24.630987
#> 11 14.28482 14.892091
#> 12 22.24619 21.222769
#> 13 29.26271 30.036844
#> 14 20.51424 22.559457
#> 15 21.88229 22.216354Selection on colony values
Since the aim of selection is to select the best individuals or
colonies for the reproduction, we could select the best colony in our
apiary based on either genetic or phenotypic value for grafting the new
generation of virgin queens. We can use the function
selectColonies() that takes a matrix of colony values (the
output of calcColonyValue() function). The default behavior
is to select the colonies with the highest value (argument
selectTop set to TRUE), but you can also
select the colonies with the lowest values (argument
selectTop set to FALSE).
# Select the best colony based on gv
selectColonies(apiary, n = 1, by = colonyGv)
#> An object of class "MultiColony"
#> Number of colonies: 1
#> Are empty: 0
#> Are NULL: 0
#> Have split: 0
#> Have swarmed: 0
#> Have superseded: 0
#> Have collapsed: 0
#> Are productive: 0
# Select the best colony based on phenotype
selectColonies(apiary, n = 1, by = colonyPheno)
#> An object of class "MultiColony"
#> Number of colonies: 1
#> Are empty: 0
#> Are NULL: 0
#> Have split: 0
#> Have swarmed: 0
#> Have superseded: 0
#> Have collapsed: 0
#> Are productive: 0The same functionality is implemented in pullColonies()
and removeColonies().
Honey yield and Calmness
In this section we expand simulation to two uncorrelated colony traits with queen and workers effects, honey yield and calmness. We follow the same recipe as in the previous section where we simulated only one colony trait.
We first reinitialize the global simulation parameters because we will define new traits. For honey yield we will use the same parameters as before, while for calmness trait we will assume that the trait is scored continuously in such a way that negative values are undesirable and positive values are desirable with zero being population mean. We will further assume the same variances for calmness as for honey yield, and a genetic (and environmental) correlation between the queen and workers effects of -0.4 (and 0.2) for calmness. We assume no genetic or environmental correlation between honey yield and calmness. Beware, this is just an example to show you how to simulate multiple colony traits - we have made up these parameters - please use literature estimates in your simulations!
# 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)We continue by creating a base population of virgin queens and from them an apiary with 10 full-sized colonies.
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
#> An object of class "MultiColony"
#> Number of colonies: 15
#> Are empty: 0
#> Are NULL: 0
#> Have split: 0
#> Have swarmed: 0
#> Have superseded: 0
#> Have collapsed: 0
#> Are productive: 15We can again inspect the genetic (and phenotypic) values of all
individuals in each colony and whole apiary with get*Gv()
and get*Pheno() functions. Now, the output contains four
traits representing the queen and workers effect for honey yield and
calmness. These functions also take an nInd argument to
sample a number of individuals along with their values.
getQueenGv(apiary) |> head(n = 4)
#> $`1`
#> yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 6 9.933572 0.05710443 0.9857495 -0.03510868
#>
#> $`2`
#> yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 7 7.66001 0.3691692 0.8552613 -0.1114546
#>
#> $`3`
#> yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 8 8.725186 0.1610305 -0.4964208 0.1783041
#>
#> $`4`
#> yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 9 9.031904 0.08626311 1.276438 -0.04014436
getWorkersPheno(apiary, nInd = 3) |> head(n = 4)
#> $`1`
#> yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 521 11.641016 -0.11109681 -0.6877905 0.34316991
#> 522 8.026872 0.15567208 0.5348766 0.01965276
#> 523 7.845169 -0.01430521 0.6852825 0.28820129
#>
#> $`2`
#> yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 721 6.693340 0.01507982 0.2445978 -0.001231914
#> 722 8.397202 0.51950866 2.4859809 -0.100405795
#> 723 10.291019 0.48789481 1.0408610 -0.008253788
#>
#> $`3`
#> yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 921 10.90046 0.08663074 -1.418556 0.15404511
#> 922 11.00305 0.54165852 3.530549 0.20706255
#> 923 10.23667 0.34108318 -2.043590 0.08833385
#>
#> $`4`
#> yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 1121 7.630843 -0.1305089 -2.0353731 -0.4029159
#> 1122 13.192528 0.3641863 0.4024319 0.1122139
#> 1123 11.200786 -0.1157791 -0.9056648 -0.1681684Now, we calculate colony genetic and phenotypic values for all
colonies in the apiary. Since we are simulating two traits, honey yield
and calmness, we have two ways to calculate corresponding colony values.
The first way is to use the default mapCasteToColony*()
function in calcColony*() and only define additional
arguments as shown here:
colonyValues <- calcColonyPheno(apiary,
queenTrait = c("yieldQueenTrait", "calmQueenTrait"),
workersTrait = c("yieldWorkersTrait", "calmWorkersTrait"),
traitName = c("yield", "calmness"),
checkProduction = c(TRUE, FALSE)) |> as.data.frame()
colonyValues
#> yield calmness
#> 1 19.50487 0.1444758
#> 2 32.55951 -5.6222774
#> 3 21.32290 6.6425657
#> 4 17.78464 -8.3617643
#> 5 22.36657 -10.4064331
#> 6 15.75597 3.3906119
#> 7 15.51778 -4.4352700
#> 8 18.48535 5.2294245
#> 9 14.49118 -2.1968434
#> 10 14.29748 -3.5123717
#> 11 14.14055 1.4324777
#> 12 26.90330 1.9021851
#> 13 26.05325 -1.8664478
#> 14 25.64785 4.9913782
#> 15 15.96541 6.4746797The second way is to create our own mapping function. An equivalent
outcome to the above is shown below just to demonstrate use of your own
function, but we are simply just reusing
mapCasteToColonyPheno() twice;)
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
#> yield calmness
#> 1 19.50487 0.1444758
#> 2 32.55951 -5.6222774
#> 3 21.32290 6.6425657
#> 4 17.78464 -8.3617643
#> 5 22.36657 -10.4064331
#> 6 15.75597 3.3906119
#> 7 15.51778 -4.4352700
#> 8 18.48535 5.2294245
#> 9 14.49118 -2.1968434
#> 10 14.29748 -3.5123717
#> 11 14.14055 1.4324777
#> 12 26.90330 1.9021851
#> 13 26.05325 -1.8664478
#> 14 25.64785 4.9913782
#> 15 15.96541 6.4746797Again, we can now select the best colony based on the best phenotypic
value for either yield, calmness, or an index of both. Let’s say that
both traits are equally important so we select on a weighted sum of both
of them - we will use the AlphaSimR selIndex() function
that enables this calculation along with scaling. We will represent the
index such that it has a mean of 100 and standard deviation of 10
units.
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)
#> [1] 14We see that we selected colony with ID “4”, but we would be selecting a different colony based on different selection criteria (yield, calmness, or index).
Strength and honey yield
In this section we change simulation to two traits where the phenotype realisation of the first trait affects the phenotype realisation of the second trait. Specifically, we will assume that queen’s fecundity, and hence the number of workers, is under the genetic affect of the queen and her environment. Furthermore, we will assume as before that colony honey yield is due to the queen effect and workers effect. Since the value of the workers effect depends on then number of workers, we obtain correlation between fecundity and honey yield, even if these traits would be uncorrelated on the queen level. We emphasise that this is just an example and the biology of these traits might be different.
We follow the same logic as before and simulate three traits that will contribute to two colony traits, queen’s fecundity, that is colony strength, and honey yield. We assume that fecundity is only due to the queen (and not the workers), hence we simulate only the queen effect for this trait. For honey yield we again assume that both the queen and workers contribute to the colony value. For speed of simulation we only simulate 100 workers per colony on average and split honey yield mean between the queen and workers. We measure fecundity with the number of workers, which is a count variable and for such variables Poisson distribution is a good model. This distribution has just one parameter (lambda) that represents both the mean and variance of the variable. To this end we set phenotypic variance to 100 and split it into 25 for genetic and 65 for environmental variance. As before we warn that these are just exemplary values to demonstrate the code functionality and do not necessarily reflect published values!
# 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)We continue by creating an apiary with 10 colonies.
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")Let’s explore queen’s genetic and phenotypic values for fecundity and
honey yield. The below printouts show quite some variation in fecundity
between queens at the genetic, but particularly phenotypic level. This
is a small example, so we should not put too much into correlations
between these three variables. However, if you restart this simulation
many times, you will notice zero correlation on average between
fecundityQueenTrait and the other two traits and negative
correlation on average between yieldQueenTrait and
yieldWorkersTrait. Just like we defined in the global
simulation parameters.
#> fecundityQueenTrait yieldQueenTrait yieldWorkersTrait
#> 6 97.67154 12.376273 -0.048847064
#> 7 100.88150 8.275271 0.248410499
#> 8 100.63731 8.450588 0.154642277
#> 9 92.66094 10.235552 0.121927900
#> 10 91.08507 10.779584 0.122417518
#> 11 103.51349 9.477871 0.247591019
#> 12 101.47314 8.755354 0.242193354
#> 13 107.55632 9.320233 -0.039866221
#> 14 108.24194 10.200024 -0.001463130
#> 15 96.34383 10.507022 -0.006119282
#> 16 96.31990 9.957808 -0.015126904
#> 17 94.33700 9.696282 0.239329927
#> 18 98.66068 10.569513 0.085806929
#> 19 104.21648 10.502979 0.031686311
#> 20 104.53985 10.323235 0.102521413
#> fecundityQueenTrait yieldQueenTrait yieldWorkersTrait
#> fecundityQueenTrait 1.00000000 0.00158083 0.2795059
#> yieldQueenTrait 0.00158083 1.00000000 0.2709895
#> yieldWorkersTrait 0.27950586 0.27098953 1.0000000
We next build-up colonies in the apiary. But instead of building them
all up to the same fixed number of workers, we build them up according
to queen’s fecundity. For that we use the sampling function
nWorkersColonyPhenotype(), that samples the number of
workers based on phenotypes of colony members, in our case
fecundityQueenTrait in queens. Correspondingly, each colony
will have a different number of workers. Read more about this function
in it’s help page.
apiary <- buildUp(apiary, nWorkers = nWorkersColonyPhenotype,
queenTrait = "fecundityQueenTrait")
cbind(nWorkers = nWorkers(apiary), queenPheno)
#> nWorkers fecundityQueenTrait yieldQueenTrait yieldWorkersTrait
#> 1 88 88.35647 8.915364 -0.140314346
#> 2 94 97.86836 9.698944 0.438368728
#> 3 93 93.40605 7.896838 0.077074693
#> 4 87 87.37309 8.436295 0.180431645
#> 5 84 83.71399 11.845748 0.230197150
#> 6 108 107.69697 9.278812 0.298073295
#> 7 107 106.58665 10.491046 0.490729138
#> 8 108 108.38055 9.757189 -0.001839875
#> 9 93 93.35867 10.408413 -0.059982182
#> 10 95 94.80427 8.602538 -0.304760548
#> 11 90 89.73196 10.682694 -0.024289681
#> 12 94 93.81233 9.329215 0.372988877
#> 13 91 99.69227 9.793778 -0.017172019
#> 14 93 97.71532 14.287510 0.330001152
#> 15 91 94.63905 11.499326 -0.035917514
help(nWorkersColonyPhenotype)To compute the colony value for honey yield, we again employ the
calcColonyPheno() function. Correlating the queen and
colony values we will now see a positive correlation because our
individual to colony mapping function sums workers effect across all
workers and the more workers there are the larger the sum.
#> nWorkers fecundityQueenTrait yieldQueenTrait
#> nWorkers 1.0000000 0.93712454 -0.11375783
#> fecundityQueenTrait 0.9371245 1.00000000 0.00158083
#> yieldQueenTrait -0.1137578 0.00158083 1.00000000
#> yieldWorkersTrait 0.2746720 0.27950586 0.27098953
#> yield 0.2911165 0.37189536 0.13431821
#> yieldWorkersTrait yield
#> nWorkers 0.2746720 0.2911165
#> fecundityQueenTrait 0.2795059 0.3718954
#> yieldQueenTrait 0.2709895 0.1343182
#> yieldWorkersTrait 1.0000000 0.8180122
#> yield 0.8180122 1.0000000