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.
# 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 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)
# 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 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))
# 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))