AlphaSimR: Implementing a breeding program II
This vignette we will perform a simulation of a maize breeding program to measure the effect of switching from initial evaluation of DH lines on testcross performance. It was originally developed by Dr. Chris Gaynor, and it was made available at link. We used the main script and we made some punctual changes.
The simulation models 15 years of burn-in for a breeding using the testcross scheme followed by 20 years of breeding in either the testcross or per se scenarios. The design of the breeding program is based on the breeding program on page 211 of “Essentials of Plant Breeding” (Bernardo 2023). The rationale for this comparison is based on a recommendation of Troyer and Wellin (2009). The number of QTL and the values for dominance degree have been chosen to approximately match the long-term rates of genetic gain shown in figure 2 of that paper (Troyer and Wellin 2009).
To run this simulation, run this main vignette. Assuming all other scripts are in your working directly, this script will automatically call all necessary scripts. It will, at the end, save the results to files called “Scenario1.rds” and “Scenario2.rds”.
This is the first script of the simulation. All other scripts are called by this script.
## Loading required package: R6
# Load global parameters
source("maizeBP/GlobalParameters.R")
nReps = 2 #Number of replications for experiment
burninYears = 15
futureYears = 20
# Initialize variables for results
hybridCorr = inbredMean = hybridMean = inbredVar = hybridVar =
rep(NA_real_,burninYears+futureYears)
output = list(inbredMean=NULL,
inbredVar=NULL,
hybridMean=NULL)
# File to store the results
saveRDS(output,"results/Maize_Scenario1.rds")
saveRDS(output,"results/Maize_Scenario2.rds")
Two distinct scenarios will be implemented, both based on phenotypic truncated selection on a reciprocal recurrently selection program.
Scenario 1: maize breeding program with two heterotic groups
(Male and females). The parental selection is made in the YT3/YT4
populations.
Scenario 2: maize breeding program with two heterotic groups
(Male and females). The early parental selection is made, and the
candidate to parents come from the YT2 population.
# Loop through replications
for(rep in 1:nReps){
cat("Working on repetition:",rep,"\n")
# Create initial parents and set testers and hybrid parents
source("maizeBP/CreateParents.R")
# Fill breeding pipeline with unique individuals from initial parents
source("maizeBP/FillPipeline.R")
# p-values for GxY effects
P = runif(burninYears+futureYears)
# Cycle years
for(year in 1:burninYears){ #Change to any number of desired years
cat("Working on year:",year,"\n")
p = P[year]
source("maizeBP/UpdateParents.R") #Pick new parents based on last year's data
source("maizeBP/UpdateTesters.R") #Pick new testers and hybrid parents
source("maizeBP/AdvanceYear.R") #Advances yield trials by a year
source("maizeBP/UpdateResults.R") #Track summary data
}
#Save burn-in to load later for scenario 2
save.image("results/BurnIn_rep.rda")
####>>>-------------------------
#######-- Scenario 1
####>>>-------------------------
cat("Working on Scenario 1\n")
for(year in (burninYears+1):(burninYears+futureYears)){
cat("Working on year:",year,"\n")
p = P[year]
source("maizeBP/UpdateParents_Scen1.R") #Pick new parents based on last year's data
source("maizeBP/UpdateTesters.R") #Pick new testers and hybrid parents
source("maizeBP/AdvanceYear.R") #Advances yield trials by a year
source("maizeBP/UpdateResults.R") #Track summary data
}
# Report results for scenario 1
output = readRDS("results/Maize_Scenario1.rds")
output = list(inbredMean=rbind(output$inbredMean,inbredMean),
inbredVar=rbind(output$inbredVar,inbredVar),
hybridMean=rbind(output$hybridMean,hybridMean))
saveRDS(output,"results/Maize_Scenario1.rds")
####>>>-------------------------
#######-- Scenario 2
####>>>-------------------------
load("results/BurnIn_rep.rda")
cat("Working on Scenario 2\n")
for(year in (burninYears+1):(burninYears+futureYears)){
cat("Working on year:",year,"\n")
p = P[year]
source("maizeBP/UpdateParents_Scen2.R") #Pick new parents based on last year's data
source("maizeBP/UpdateTesters.R") #Pick new testers and hybrid parents
source("maizeBP/AdvanceYear.R") #Advances yield trials by a year
source("maizeBP/UpdateResults.R") #Track summary data
}
# Report results for scenario 2
output = readRDS("results/Maize_Scenario2.rds")
output = list(inbredMean=rbind(output$inbredMean,inbredMean),
inbredVar=rbind(output$inbredVar,inbredVar),
hybridMean=rbind(output$hybridMean,hybridMean))
saveRDS(output,"results/Maize_Scenario2.rds")
#Delete tmp file
file.remove("results/BurnIn_rep.rda")
}
## Working on repetition: 1
## FillPipeline year: 1 of 6
## FillPipeline year: 2 of 6
## FillPipeline year: 3 of 6
## FillPipeline year: 4 of 6
## FillPipeline year: 5 of 6
## FillPipeline year: 6 of 6
## Working on year: 1
## Working on year: 2
## Working on year: 3
## Working on year: 4
## Working on year: 5
## Working on year: 6
## Working on year: 7
## Working on year: 8
## Working on year: 9
## Working on year: 10
## Working on year: 11
## Working on year: 12
## Working on year: 13
## Working on year: 14
## Working on year: 15
## Working on Scenario 1
## Working on year: 16
## Working on year: 17
## Working on year: 18
## Working on year: 19
## Working on year: 20
## Working on year: 21
## Working on year: 22
## Working on year: 23
## Working on year: 24
## Working on year: 25
## Working on year: 26
## Working on year: 27
## Working on year: 28
## Working on year: 29
## Working on year: 30
## Working on year: 31
## Working on year: 32
## Working on year: 33
## Working on year: 34
## Working on year: 35
## Working on Scenario 2
## Working on year: 16
## Working on year: 17
## Working on year: 18
## Working on year: 19
## Working on year: 20
## Working on year: 21
## Working on year: 22
## Working on year: 23
## Working on year: 24
## Working on year: 25
## Working on year: 26
## Working on year: 27
## Working on year: 28
## Working on year: 29
## Working on year: 30
## Working on year: 31
## Working on year: 32
## Working on year: 33
## Working on year: 34
## Working on year: 35
## Working on repetition: 2
## FillPipeline year: 1 of 6
## FillPipeline year: 2 of 6
## FillPipeline year: 3 of 6
## FillPipeline year: 4 of 6
## FillPipeline year: 5 of 6
## FillPipeline year: 6 of 6
## Working on year: 1
## Working on year: 2
## Working on year: 3
## Working on year: 4
## Working on year: 5
## Working on year: 6
## Working on year: 7
## Working on year: 8
## Working on year: 9
## Working on year: 10
## Working on year: 11
## Working on year: 12
## Working on year: 13
## Working on year: 14
## Working on year: 15
## Working on Scenario 1
## Working on year: 16
## Working on year: 17
## Working on year: 18
## Working on year: 19
## Working on year: 20
## Working on year: 21
## Working on year: 22
## Working on year: 23
## Working on year: 24
## Working on year: 25
## Working on year: 26
## Working on year: 27
## Working on year: 28
## Working on year: 29
## Working on year: 30
## Working on year: 31
## Working on year: 32
## Working on year: 33
## Working on year: 34
## Working on year: 35
## Working on Scenario 2
## Working on year: 16
## Working on year: 17
## Working on year: 18
## Working on year: 19
## Working on year: 20
## Working on year: 21
## Working on year: 22
## Working on year: 23
## Working on year: 24
## Working on year: 25
## Working on year: 26
## Working on year: 27
## Working on year: 28
## Working on year: 29
## Working on year: 30
## Working on year: 31
## Working on year: 32
## Working on year: 33
## Working on year: 34
## Working on year: 35
# Loading the results
Scenario1 = readRDS("results/Maize_Scenario1.rds")
Scenario2 = readRDS("results/Maize_Scenario2.rds")
# Getting the mean across repetitions
MeanScen1 = lapply(Scenario1, function(x) apply(x,2,mean))
MeanScen2 = lapply(Scenario2, function(x) apply(x,2,mean))
# Plot hybrid mean of genetic values over time
meanRanges = range(c(MeanScen1$inbredMean[15:35], MeanScen2$inbredMean[15:35]))
plot(x = 15:35, y = MeanScen1$inbredMean[15:35], type = "l", col = "black", lwd = 3,
xlab = "Year", ylab = "Mean of genetic values", ylim = meanRanges)
lines(x = 15:35, y = MeanScen2$inbredMean[15:35], type = "l", col = "blue", lwd = 3)
legend(x = "topleft", legend = c('Conv', 'EarlySel'), title = "Scenarios",
lwd = 3, lty = c(1, 2), col = c("black", "blue"), bty = "n")
University of Florida, mresende@ufl.edu↩︎
University of Florida, deamorimpeixotom@ufl.edu↩︎