Introduction

In this exercise, you will simulate one generation of single-trait selection and evaluate the response to selection under different:

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

Base population

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

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

# Load AlphaSimR, simulate founder genomes, define a trait, and simulate a base population
library(AlphaSimR)
## Loading required package: R6
founderGenomes = runMacs(nInd = 100,
                         nChr = 10,
                         segSites = 100,
                         species = "MAIZE")
SP = SimParam$new(founderGenomes)
SP$addTraitA(nQtlPerChr = 100, mean = 10, var = 1)
basePop = newPop(founderGenomes)

Accuracy of selection (one generation)

First, vary the accuracy of selection by varying the heritability of phenotypes used to identify the superior individuals. Assume heritability of 0.5 and select 50 individuals with the highest phenotype values.

# Define parameters for the scenario with trait heritability of 0.5 and 50 selected
heritability = 0.5
nSelected = 50

# Phenotype base population
basePop = setPheno(basePop, h2 = heritability)

# Select the best performing individuals according to their phenotype
basePopSelected = selectInd(pop = basePop,
                            nInd = nSelected,
                            use = "pheno")

# Create a new population and phenotype it
newPop = randCross(basePopSelected, nCrosses = nInd(basePop))
newPop = setPheno(newPop, h2 = heritability)

# Evaluate observed response to selection between generations
# (as the difference between the mean of genetic values)
(deltaGObserved_h205 = meanG(newPop) - meanG(basePop))
##    Trait1 
## 0.6647666
# Quantify the intensity of selection
# (using the function selInt())
(selInt_h205 = selInt(p = nSelected / nInd(basePop)))
## [1] 0.7978846
# Quantify the accuracy of selection
# (as the correlation between genetic value and selection criterion)
(selAcc_h205 = cor(gv(basePop), pheno(basePop)))
##           Trait1
## Trait1 0.6453998
# Quantify the standard deviation of genetic values
(sdGen_h205 = sqrt(varG(basePop)))
##        Trait1
## Trait1      1
# Predict expected response to selection between generations
(deltaGExpected_h205 = selInt_h205 * selAcc_h205 * sdGen_h205)
##           Trait1
## Trait1 0.5149545

Now repeat the above code, but change trait heritability to 0.1 and save the results into objects deltaGObserved_h201, selInt_h201, selAcc_h201, sdGen_h201, and deltaGExpected_h201. Then repeat the above code by changing trait heritability to 0.9 and save the results into objects deltaGObserved_h209, selInt_h209, selAcc_h209, sdGen_h209, and deltaGExpected_h209.

Start with heritability of 0.1.

# Define parameters for the scenario with trait heritability of 0.1 and 50 selected
heritability = 0.1
nSelected = 50

# Phenotype base population
basePop = setPheno(basePop, h2 = heritability)

# Select the best performing individuals according to their phenotype
basePopSelected = selectInd(pop = basePop,
                            nInd = nSelected,
                            use = "pheno")

# Create a new population and phenotype it
newPop = randCross(basePopSelected, nCrosses = nInd(basePop))
newPop = setPheno(newPop, h2 = heritability)

# Evaluate observed response to selection between generations
# (as the difference between the mean of genetic values)
(deltaGObserved_h201 = meanG(newPop) - meanG(basePop))
##    Trait1 
## 0.2586856
# Quantify the intensity of selection
# (using the function selInt())
(selInt_h201 = selInt(p = nSelected / nInd(basePop)))
## [1] 0.7978846
# Quantify the accuracy of selection
# (as the correlation between genetic value and selection criterion)
(selAcc_h201 =  cor(gv(basePop), pheno(basePop)))
##           Trait1
## Trait1 0.2834074
# Quantify the standard deviation of genetic values
(sdGen_h201 = sqrt(varG(basePop)))
##        Trait1
## Trait1      1
# Predict expected response to selection between generations
(deltaGExpected_h201 = selInt_h201 * selAcc_h201 * sdGen_h201)
##           Trait1
## Trait1 0.2261264

Now with heritability of 0.9.

# Define parameters for the scenario with trait heritability of 0.9 and 50 selected
heritability = 0.9
nSelected = 50

# Phenotype base population
basePop = setPheno(basePop, h2 = heritability)

# Select the best performing individuals according to their phenotype
basePopSelected = selectInd(pop = basePop,
                            nInd = nSelected,
                            use = "pheno")

# Create a new population and phenotype it
newPop = randCross(basePopSelected, nCrosses = nInd(basePop))
newPop = setPheno(newPop, h2 = heritability)

# Evaluate observed response to selection between generations
# (as the difference between the mean of genetic values)
(deltaGObserved_h209 = meanG(newPop) - meanG(basePop))
##    Trait1 
## 0.7964692
# Quantify the intensity of selection
# (using the function selInt())
(selInt_h209 = selInt(p = nSelected / nInd(basePop)))
## [1] 0.7978846
# Quantify the accuracy of selection
# (as the correlation between genetic value and selection criterion)
(selAcc_h209 =  cor(gv(basePop), pheno(basePop)))
##           Trait1
## Trait1 0.9591685
# Quantify the standard deviation of genetic values
(sdGen_h209 = sqrt(varG(basePop)))
##        Trait1
## Trait1      1
# Predict expected response to selection between generations
(deltaGExpected_h209 = selInt_h209 * selAcc_h209 * sdGen_h209)
##           Trait1
## Trait1 0.7653058

Intensity of selection (one generation)

Now repeat the above code, by changing trait heritability back to 0.5 and changing the number of selected individuals to 20 and save the results into objects deltaGObserved_n20, selInt_n20, selAcc_n20, sdGen_n20, and deltaGExpected_n20. Then change the number of selected individuals to 10 and save the results into objects deltaGObserved_n10, selInt_10, selAcc_n10, sdGen_n10, and deltaGExpected_n10.

Start with 20 selected individuals.

# Define parameters for the scenario with trait heritability of 0.5 and 20 selected
heritability = 0.5
nSelected = 20

# Phenotype base population
basePop = setPheno(basePop, h2 = heritability)

# Select the best performing individuals according to their phenotype
basePopSelected = selectInd(pop = basePop,
                            nInd = nSelected,
                            use = "pheno")

# Create a new population and phenotype it
newPop = randCross(basePopSelected, nCrosses = nInd(basePop))
newPop = setPheno(newPop, h2 = heritability)

# Evaluate observed response to selection between generations
# (as the difference between the mean of genetic values)
(deltaGObserved_n20 = meanG(newPop) - meanG(basePop))
##   Trait1 
## 1.157193
# Quantify the intensity of selection
# (using the function selInt())
(selInt_n20 = selInt(p = nSelected / nInd(basePop)))
## [1] 1.39981
# Quantify the accuracy of selection
# (as the correlation between genetic value and selection criterion)
(selAcc_n20 =  cor(gv(basePop), pheno(basePop)))
##           Trait1
## Trait1 0.7535116
# Quantify the standard deviation of genetic values
(sdGen_n20 = sqrt(varG(basePop)))
##        Trait1
## Trait1      1
# Predict expected response to selection between generations
(deltaGExpected_n20 = selInt_n20 * selAcc_n20 * sdGen_n20)
##          Trait1
## Trait1 1.054773

Now with 10 selected individuals.

# Define parameters for the scenario with trait heritability of 0.5 and 10 selected
heritability = 0.5
nSelected = 10

# Phenotype base population
basePop = setPheno(basePop, h2 = heritability)

# Select the best performing individuals according to their phenotype
basePopSelected = selectInd(pop = basePop,
                            nInd = nSelected,
                            use = "pheno")

# Create a new population and phenotype it
newPop = randCross(basePopSelected, nCrosses = nInd(basePop))
newPop = setPheno(newPop, h2 = heritability)

# Evaluate observed response to selection between generations
# (as the difference between the mean of genetic values)
(deltaGObserved_n10 = meanG(newPop) - meanG(basePop))
##   Trait1 
## 1.269964
# Quantify the intensity of selection
# (using the function selInt())
(selInt_n10 = selInt(p = nSelected / nInd(basePop)))
## [1] 1.754983
# Quantify the accuracy of selection
# (as the correlation between genetic value and selection criterion)
(selAcc_n10 = cor(gv(basePop), pheno(basePop)))
##           Trait1
## Trait1 0.6462686
# Quantify the standard deviation of genetic values
(sdGen_n10 = sqrt(varG(basePop)))
##        Trait1
## Trait1      1
# Predict expected response to selection between generations
(deltaGExpected_n10 = selInt_n10 * selAcc_n10 * sdGen_n10)
##          Trait1
## Trait1 1.134191

Finally, combine all the results into a data.frame and discuss the observed and expected response to selection across the scenarios.

results = data.frame(h2 = c(0.5, 0.1, 0.9, 0.5, 0.5),
                     nSelected = c(50, 50, 50, 20, 10),
                     deltaGObserved = c(deltaGObserved_h205, deltaGObserved_h201, deltaGObserved_h209, deltaGObserved_n20, deltaGObserved_n10),
                     deltaGExpected = c(deltaGExpected_h205, deltaGExpected_h201, deltaGExpected_h209, deltaGExpected_n20, deltaGExpected_n10),
                     selInt = c(selInt_h205, selInt_h201, selInt_h209, selInt_n20, selInt_n10),
                     selAcc = c(selAcc_h205, selAcc_h201, selAcc_h209, selAcc_n20, selAcc_n10),
                     sdGen = c(sdGen_h205, sdGen_h201, sdGen_h209, sdGen_n20, sdGen_n10))
print(results)
##    h2 nSelected deltaGObserved deltaGExpected    selInt    selAcc sdGen
## 1 0.5        50      0.6647666      0.5149545 0.7978846 0.6453998     1
## 2 0.1        50      0.2586856      0.2261264 0.7978846 0.2834074     1
## 3 0.9        50      0.7964692      0.7653058 0.7978846 0.9591685     1
## 4 0.5        20      1.1571934      1.0547727 1.3998096 0.7535116     1
## 5 0.5        10      1.2699644      1.1341906 1.7549833 0.6462686     1

What do the results show you about changing these parameters?

Response to selection increased as we increased trait heritability and therefore accuracy of selection on phenotypes (with 10 selected individuals response to selection was 0.2586856 with heritability of 0.1, 0.6647666 with heritability of 0.5, and 0.7964692 with heritability of 0.9).

Response to selection increased as we decreased the number of selected individuals (and therefore increased intensity of selection) (with heritability 0.50 response to selection was 0.6647666 with 50% selected individuals, 1.1571934 with 20% selected individuals, and 1.2699644 with 10% selected individuals)

The observed and expected responses to selection matched well, but there were some deviations between theory and simulation.

EXTRA: Intensity of selection (many generations)

Lastly, use the code below and evaluate the effect of the intensity of selection on response to selection over 50 generations in the following three steps:

  1. Run the simulation with 10 selected individuals and save results into objects meanGAll_n10 and varGAll_n10.

  2. Repeat the simulation with 50 selected individuals and save results into objects meanGAll_n50 and varGAll_n50.

  3. Plot trend in the mean and variance of genetic values over generations for both scenarios and discuss results.

Start with 10 selected individuals.

# Set simulation parameters for the scenario with trait heritability of 0.5 and 10 selected
heritability = 0.5
nSelected = 10

# Phenotype base population
basePop = setPheno(basePop, h2 = heritability)

# Select the best performing individuals according to their phenotype
basePopSelected = selectInd(pop = basePop,
                            nInd = nSelected,
                            use = "pheno")

# Allocate vectors
nGenerations = 50 + 1 # +1 to store starting generation
meanGAll = numeric(nGenerations)
varGAll = numeric(nGenerations)

# Save the starting values
meanGAll[1] = meanG(basePop)
varGAll[1] = varG(basePop)

# To make the for loop below simpler we will make a copy of the object basePopSelected
newPopSelected = basePopSelected

# Selection over many generations
for (generation in 1:(nGenerations - 1)) {
  # Cross parents, phenotype progeny, and select new parents
  newPop = randCross(newPopSelected, nCrosses = nInd(basePop))
  newPop = setPheno(newPop, h2 = heritability)
  newPopSelected = selectInd(pop = newPop,
                             nInd = nSelected,
                             use = "pheno")
  # Save summaries
  meanGAll[1 + generation] = meanG(newPop)
  varGAll[1 + generation] = varG(newPop)
}

# Now save these outputs by copying the objects
meanGAll_n10 = meanGAll
varGAll_n10 = varGAll

Now with 50 selected individuals.

# Set simulation parameters for the scenario with trait heritability of 0.5 and 50 selected
heritability = 0.5
nSelected = 50

# Phenotype base population
basePop = setPheno(basePop, h2 = heritability)

# Select the best performing individuals according to their phenotype
basePopSelected = selectInd(pop = basePop,
                            nInd = nSelected,
                            use = "pheno")

# Allocate vectors
nGenerations = 50 + 1 # +1 to store starting generation
meanGAll = numeric(nGenerations)
varGAll = numeric(nGenerations)

# Save the starting values
meanGAll[1] = meanG(basePop)
varGAll[1] = varG(basePop)

# To make the for loop below simpler we will make a copy of the object basePopSelected
newPopSelected = basePopSelected

# Selection over many generations
for (generation in 1:(nGenerations - 1)) {
  # Cross parents, phenotype progeny, and select new parents
  newPop = randCross(newPopSelected, nCrosses = nInd(basePop))
  newPop = setPheno(newPop, h2 = heritability)
  newPopSelected = selectInd(pop = newPop,
                             nInd = nSelected,
                             use = "pheno")
  # Save summaries
  meanGAll[1 + generation] = meanG(newPop)
  varGAll[1 + generation] = varG(newPop)
}

# Now save these outputs by copying the objects
meanGAll_n50 = meanGAll
varGAll_n50 = varGAll
par(mfrow = c(2, 1),
    mar = c(4, 4, 1, 1))

# Plot mean of genetic values over time
meanRanges = range(c(meanGAll_n10, meanGAll_n50))
plot(x = 1:nGenerations, y = meanGAll_n10, type = "l", col = "black", lwd = 3,
     xlab = "Generation", ylab = "Mean of genetic values", ylim = meanRanges)
lines(x = 1:nGenerations, y = meanGAll_n50, type = "l", col = "purple", lty = 2, lwd = 3)
legend(x = "topleft", legend = c(10, 50), title = "nSelected",
       lwd = 3, lty = c(1, 2), col = c("black", "purple"), bty = "n")

# Plot variance of genetic values over time
varRanges = range(c(varGAll_n10, varGAll_n50))
plot(x = 1:nGenerations, y = varGAll_n10, type = "l", col = "black", lwd = 3,
     xlab = "Generation", ylab = "Variance of genetic values", ylim = varRanges)
lines(x = 1:nGenerations, y = varGAll_n50, type = "l", col = "purple", lty = 2, lwd = 3)
legend(x = "topright", legend = c(10, 50), title = "nSelected",
       lwd = 3, lty = c(1, 2), col = c("black", "purple"), bty = "n")

What do the results show you about changing these parameters?

The results show that mean genetic value increased in both scenarios, when we selected 10 or 50 individuals. While both scenarios achieved comparable mean genetic value after 50 generations, the response was initially much faster in the scenario with the higher intensity of selection and then the response to selection started to level of and both scenarios ended at comparable means.

Trends in the variance of genetic values shows that the scenario with the higher intensity of selection depleted genetic variance faster, which slowed down response to selection in the long-term compared to the scenario with the lower intensity of selection.