Introduction

In this exercise, you will simulate a population and analyse its trait variation in five steps:

For each step, we have given you instructions with an AlphaSimR template code to complete (replace ??? with an appropriate code).

This exercise deliberately goes beyond the material we have covered up to now. To help you on this pathway, we indicate which functions should be used and we point to their documentation. Embracing this growth mindset is important for mastering AlphaSimR and combining it with other R functionality.

Simulate founding genomes

# Clean the working environment
rm(list = ls())

# Set the default plot layout
par(mfrow = c(1, 1))

library(AlphaSimR)
## Loading required package: R6
# Use the following parameters: 
#   * 100 individuals (this simulation will need more time than before)
#   * 10 chromosomes
#   * 100 segregating loci (sites) per chromosomes
#   * MAIZE species
founderGenomes = runMacs(nInd = 100,
                         nChr = 10,
                         segSites = 100,
                         species = "MAIZE")
# Create the simulation parameters object SP from founding genomes
SP = SimParam$new(founderGenomes)

Define traits

# Define four traits with additive genetic effects where:
#   * Trait 1 has mean 10, genetic variance 1, and is affected by   1 QTL/chromosome
#   * Trait 2 has mean 10, genetic variance 1, and is affected by  10 QTL/chromosome
#   * Trait 3 has mean 10, genetic variance 1, and is affected by 100 QTL/chromosome
#   * Trait 4 has mean 10, genetic variance 3, and is affected by 100 QTL/chromosome

# Add trait 1
SP$addTraitA(mean = 10, var = 1, nQtlPerChr =   1)

# Add trait 2
SP$addTraitA(mean = 10, var = 1, nQtlPerChr =  10)

# Add trait 3
SP$addTraitA(mean = 10, var = 1, nQtlPerChr = 100)

# Add trait 4
SP$addTraitA(mean = 10, var = 3, nQtlPerChr = 100)

Inspect traits

# Number of traits
SP$nTraits
## [1] 4
# Distribution of additive QTL effects for each trait
par(mfrow = c(2, 2),
    mar = c(4, 4, 1, 1))
hist(SP$traits[[1]]@addEff, xlab = "Additive effects - trait 1", main = "")
hist(SP$traits[[2]]@addEff, xlab = "Additive effects - trait 2", main = "")
hist(SP$traits[[3]]@addEff, xlab = "Additive effects - trait 3", main = "")
hist(SP$traits[[4]]@addEff, xlab = "Additive effects - trait 4", main = "")

par(mfrow = c(1, 1))

Create a population and phenotype it

# Create a population of individuals from founding genomes
maizePop = newPop(founderGenomes)

# Simulate phenotypes with the following environmental variances:
#   * Trait 1 with varE = 0
#   * Trait 2 with varE = 1
#   * Trait 3 with varE = 1
#   * Trait 4 with varE = 2
maizePop = setPheno(maizePop, varE = c(0, 1, 1, 2))

Summarise genetic and phenotype values

# Calculate mean and variance of genetic values for individuals in the maizePop
# using meanG() and varG() functions. Since we have multiple traits, these
# functions will return vectors and matrices (variances on diagonal and covariances
# on off-diagonal elements).
meanG(maizePop)
## Trait1 Trait2 Trait3 Trait4 
##     10     10     10     10
varG(maizePop)
##              Trait1       Trait2       Trait3     Trait4
## Trait1  1.000000000 -0.034734622 -0.005690818 -0.0170347
## Trait2 -0.034734622  1.000000000 -0.002897671 -0.1049792
## Trait3 -0.005690818 -0.002897671  1.000000000  0.5045344
## Trait4 -0.017034700 -0.104979151  0.504534358  3.0000000
# Calculate mean and variance of phenotype values for individuals in the maizePop
# using meanP() and varP() functions. Since we have multiple traits, these
# functions will return vectors and matrices (variances on diagonal and covariances
# on off-diagonal elements).
meanP(maizePop)
##    Trait1    Trait2    Trait3    Trait4 
## 10.000000  9.957154  9.983308 10.075273
varP(maizePop)
##             Trait1     Trait2    Trait3      Trait4
## Trait1  1.00000000 -0.1514729 0.1600604  0.08981877
## Trait2 -0.15147294  1.9135183 0.0260031 -0.48240633
## Trait3  0.16006040  0.0260031 2.0112513  0.48935422
## Trait4  0.08981877 -0.4824063 0.4893542  4.98945043
# Calculate heritabilities for simulated phenotypes in the maizePop
# Note the use of diag() function to extract only variances from the matrices
VarG = varG(maizePop)
VarP = varP(maizePop)
diag(VarG) / diag(VarP)
##    Trait1    Trait2    Trait3    Trait4 
## 1.0000000 0.5225976 0.4972029 0.6012686
# Plot phenotype values as a function of genetic values in the maizePop
par(mfrow = c(2, 2),
    mar = c(4, 4, 1, 1))
plot(x = gv(maizePop)[, 1], y = pheno(maizePop)[, 1], xlab = "gv1", ylab = "pheno1")
plot(x = gv(maizePop)[, 2], y = pheno(maizePop)[, 2], xlab = "gv2", ylab = "pheno2")
plot(x = gv(maizePop)[, 3], y = pheno(maizePop)[, 3], xlab = "gv3", ylab = "pheno3")
plot(x = gv(maizePop)[, 4], y = pheno(maizePop)[, 4], xlab = "gv4", ylab = "pheno4")

par(mfrow = c(1, 1))