AlphaSimR: Implementing a breeding program II
## Recriprocal recurrently selection program: Maize breeding 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](https://github.com/gaynorr/AlphaSimR_Examples). 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" [@bernardo2023essentials]. 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 [@troyer2009heterosis]. 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". ## Breeding program This is the first script of the simulation. All other scripts are called by this script. ```{r} # Packages library(AlphaSimR) # 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") ``` ## Scenarios implementations 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. ```{r} # 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") } ``` ## Plotting the outputs ```{r} # 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") ``` ## References ::: {#refs} :::