Simulation using AlphaSimR

Julius Mugambe

2022-02-09

Basic simple simulation using AlphaSimR

Lets install and load the package

To install use:

# install.packages('AlphaSimR')

The development version of AlphaSimR (potentially unstable) can be accessed from the devel branch on GitHub.

To install use:

# devtools::install_github(repo="gaynorr/AlphaSimR@devel")

To install with vignettes use:

# devtools::install_github(repo="gaynorr/AlphaSimR@devel", build_vignettes=TRUE)

To load the library use:

require("AlphaSimR")

In this example, I will consider an animal breeding program.

  • Discrete generations
    • about 1000 animals per generation
  • Selection will be based on true genetic value
    • select all females
    • select the best 25 males

First we clean the working environment

rm(list=ls())

Then set up the working environment

Set the location of the working directory

setwd("D:/Qgene/Simulation/AlphaSimR/Test00/")
getwd()
## [1] "D:/Qgene/Simulation/AlphaSimR/Test00"

Simulation Structure

  • Creating Founder population
  • Setting Simulation Parameters
  • Modeling the Breeding Program
  • Examining the Results

Creating the founder population

While creating the founder population, you can decide to use external SNP data that is either phased (for haplotypes) or not. You will also need a genetic map for external data.

This software uses the Markovian Coalescent Simulator (MaCS) software i.e it used the Coalescent model. This is the preferred method for the simulation of founder populations

# Creating Founder population
nInd      = 1000
nChr      = 10
segSites  = 1000

# founder population haplotypes
founderPop = runMacs(nInd,
                     nChr,
                     segSites)

Setting Simulation parameters

First, you set the founder population Add traits ( scaled by the founder population) Add SNP Chips Several other options

# Setting Simulation Parameters 

SP   = SimParam$ 
       new(founderPop)$                           # Initialize SimParam
       addTraitA(nQtlPerChr=1000, mean=0, var=1)$  # Add trait
       setSexes("yes_sys")                        # Set gender usage

Modelling the Breeding program

This is highly variable depending on the breeding program.
There is no correct way to write the code, It is by design.

# Modeling the Breeding Program

nFemale   = 500
nMale     = 25

# Initialize the population

pop       = newPop(founderPop)        #  Create initial population
for(cycle in 1:20){                   #  20 generations of breeding
  
# Select and mate
pop = selectCross(pop = pop, 
                  nFemale,
                  nMale,
                  use="gv", 
                  nCrosses=1000)
}

Examining the Results

Here you can use other packages in R

  • Plots
  • Statistical analysis
  • This step requires planning and may be an iterative process
meanG(pop)             # View new population mean
## [1] 11.88981
varG(pop)              # View new population variance
##           [,1]
## [1,] 0.4314468