Until now we have been using the breeders equation to try to compare different breeding strategies. Because we are using an equation to compute rates of genetic gain, this is called deterministic simulation. Stochastic simulation takes a different approach in that it models each process including recombiantion, phenotyping, genotyping etc. For most of the processes, there is random chance involved, thus each simulation run will give slightly different results. To be able to compute the rate of genetic gain using stochastic simulation, simulations will need to be repeated hundreds or even thousands of times. In this example we will only be doing a single simulation run.

Load the AlphaSimR package

Currently AlphaSimR and BreedingSchemeLanguage are the two best R packages for breeding program simulation. Each one has its advantages and disadvantages. AlphaSimR may be faster and its development is more active, but BreedingSchemeLanguage has very useful functions for plant breeding which makes simulating realistic breeding programs a bit easier.

library(AlphaSimR)
## Loading required package: R6

Create a founder population

This creates a set of founder haplotypes according to a certain demographic history using the MaCS software (Chen, Marjoram, and Wall 2009). The demographic history is basically how the effective population size has changed over time. The function allows users to input their own population demographic history or to select from some pre-programmed population histories. We will use the population history for Wheat, and we set inbred to TRUE so that the founder individuals are inbred. The base population would be the founders prior to inbreeding. We can set the founder population size and the number of chromosomes. For testing, its convenient to use fewer chromosomes than the actual number.

founderPop<-  runMacs(nInd=100, nChr=5, inbred=TRUE, species='WHEAT')

Create the parameters object

A simulation parameters object needs to be created and then populated. In this example the simulation paramater object is SP. To view all the parameters that can be set, use str(SP)

SP<- SimParam$new(founderPop)#create the simulation parameters object

Populate the parameters object

SP$addTraitA(nQtlPerChr=100) #add trait, by default the additive genetic variance = 1
SP$setSexes("no") #all individuals are hermaphrodites 
SP$setTrackPed(TRUE) #keeps pedigree information
SP$setTrackRec(TRUE) #keeps records of all individuals created

SP$addSnpChip(nSnpPerChr=50) #creates a marker set with 50 SNP per chromosome
SP$setVarE(h2=0.4) #adds error to achieve a heritability of 0.4

Make the generation 0 population equal to the founder population

This step seems redundant but it is needed in order to create a population object that can be populated with phenotypic information and used later for selection and recombination

gen0<- newPop(founderPop)
str(gen0, max.level = 2)
## Formal class 'Pop' [package "AlphaSimR"] with 17 slots
##   ..@ id     : chr [1:100] "1" "2" "3" "4" ...
##   ..@ mother : chr [1:100] "0" "0" "0" "0" ...
##   ..@ father : chr [1:100] "0" "0" "0" "0" ...
##   ..@ sex    : chr [1:100] "H" "H" "H" "H" ...
##   ..@ nTraits: int 1
##   ..@ gv     : num [1:100, 1] -0.923 0.557 1.195 0.272 0.341 ...
##   ..@ pheno  : num [1:100, 1] -0.0233 1.1176 0.2494 0.52 0.2364 ...
##   ..@ ebv    : num[1:100, 0 ] 
##   ..@ gxe    :List of 1
##   ..@ fixEff : int [1:100] 1 1 1 1 1 1 1 1 1 1 ...
##   ..@ reps   : num [1:100] 1 1 1 1 1 1 1 1 1 1 ...
##   ..@ misc   :List of 100
##   ..@ nInd   : int 100
##   ..@ nChr   : int 5
##   ..@ ploidy : int 2
##   ..@ nLoci  : int [1:5] 180434 181579 185492 182241 178928
##   ..@ geno   :List of 5
##   .. ..- attr(*, "dim")= int [1:2] 5 1

Phenotype generation 0

This creates phenotypes by adding error to the breeding values, just like what we have been doing in this class

gen0<- setPheno(pop=gen0)

Select individals

This selects individuals in the upper tail. The ‘use’ parameter indicates what is used for selection. The other options are estimated breeding values, true breeding values, and random

gen0sel<- selectInd(pop=gen0, nInd=20, use="pheno")

Randomly cross the selected individuals

This function randomly intermates individuals in a population. In this example we are simulating 100 crosses with 1 progeny per cross.

gen1_F1<- randCross(pop=gen0sel, nCrosses=100, nProgeny = 1)
str(gen1_F1, max.level = 2)
## Formal class 'Pop' [package "AlphaSimR"] with 17 slots
##   ..@ id     : chr [1:100] "101" "102" "103" "104" ...
##   ..@ mother : chr [1:100] "85" "85" "85" "85" ...
##   ..@ father : chr [1:100] "59" "35" "3" "91" ...
##   ..@ sex    : chr [1:100] "H" "H" "H" "H" ...
##   ..@ nTraits: int 1
##   ..@ gv     : num [1:100, 1] 1.675 0.477 1.127 1.52 1.816 ...
##   ..@ pheno  : num [1:100, 1] 2.221 0.668 2.182 1.843 0.741 ...
##   ..@ ebv    : num[1:100, 0 ] 
##   ..@ gxe    :List of 1
##   ..@ fixEff : int [1:100] 1 1 1 1 1 1 1 1 1 1 ...
##   ..@ reps   : num [1:100] 1 1 1 1 1 1 1 1 1 1 ...
##   ..@ misc   :List of 100
##   ..@ nInd   : int 100
##   ..@ nChr   : int 5
##   ..@ ploidy : int 2
##   ..@ nLoci  : int [1:5] 180434 181579 185492 182241 178928
##   ..@ geno   :List of 5
##   .. ..- attr(*, "dim")= int [1:2] 5 1

Self pollinate

Simulates a single generation of self-pollination. In this example we simulate just one progeny from self-pollination.

gen1_F2<- self(pop=gen1_F1, nProgeny = 1)

Create double-haploids

Simulates the generation of double haploid individuals. In this example we simulate just double haploid generated per individual in the F2 generation.

gen1_DH<- makeDH(pop=gen1_F2, nDH = 1)

Compare breeding values across populations

I use a loop to get the breeding values from each generation and combine it into a table. Then I create a boxplot using that table of results

pops<- c('gen0', 'gen0sel', 'gen1_F1', 'gen1_F2', 'gen1_DH')
for(i in 1:length(pops)){
  df<- data.frame(pop=pops[i], bv=get(pops[i])@gv)
  if(i==1){
    dfall<- df
  }else{
    dfall<-rbind(dfall, df)
  }
}
boxplot(dfall$bv~dfall$pop, las=3, col='grey', xlab="", ylab='breeding value')

Exercises

  1. Think of a scenario/other set of parameters that you would like to compare with what we just simulated and alter the demonstration code to simulate that scenario.

  2. Sketch a basic breeding scheme, write down each step in the process, and then simulate the breeding scheme that you described.