Assistance with Genomic Selection scenarios with Multi-Trait index

[cecilia249]

Hello,
I’m reaching out here regarding the use of AlphaSimR in breeding program simulations in rice. We’re setting up a new scenario that includes genomic selection in the first evaluation stage (E1, i.e. PYT) and incorporates a selection index (using selIndex) based on yield, grain quality and disease resistance traits.

In the phenotypic selection scenario, we successfully ran this code:
E1 = selectInd(pop = F6, nInd = 100, simParam = SP, trait = selIndex, b = wE1)

where wE1 is our vector of trait weights. However, when attempting genomic selection with the following code:

E1 = selectInd(pop = E1, nInd = 100, use = "ebv", trait = selIndex, b = wE1, simParam = SP)

we got an error: "ncol(response) == 1 is not TRUE.", which seems to indicate that the response (trait) must be a single trait. All references we've found in GitHub repositories and papers seem to focus on single-trait rather than multi-trait GS simulations.

Could you advise on whether genomic selection can be implemented with a multi-trait index, and if so, any insights into resolving this error?

Thank you very much for the time and assistance!

[gaynorr]

Hi @cecilia249

Here is a script showing the various steps with comments. Please ask questions if anything is unclear and I'd be glad to elaborate.

I also showed how to run a multi-trait model using the bWGR package. AlphaSimR's implementation of a multi-trait model really slows down with large datasets. Whereas the multi-trait model in bWGR scales nicely to large data. See this link for a good example.

library(AlphaSimR)

## Create founder haplotypes ----
founderPop = runMacs(nInd = 50, 
                     nChr = 12, 
                     segSites = 1100, 
                     inbred = TRUE)

## Set simulation parameters ----

# Initialize simulation parameters
SP = SimParam$new(founderPop)

# Choose 100 loci per chromosome to function as markers being genotyped
# These markers won't have an effect on the traits and they will be prevented
# from being selected as QTL
SP$addSnpChip(nSnpPerChr = 100)

# Adding two traits that are highly polygenic (1000 QTL per chromosome), 
# pleiotropic, and have moderate genetic correlation (0.5)
SP$addTraitA(nQtlPerChr = 1000, 
             mean = c(0,0), 
             var = c(1,1), 
             corA = 0.5*diag(2)+0.5)

# Adding a third trait controlled by fewer QTL and no inherent correlation with 
# the other two
SP$addTraitA(nQtlPerChr = 10,
             mean = 0, 
             var = 1)

# Setting heritabilities for each trait
# Note that I'm not seeting corE, so residuals are uncorrelated
SP$setVarE(h2 = c(0.2, 0.6, 0.4))

# Assign selection index weights for multitrait selection
b = c(1, -0.5, 0.3)

## Simulate a population ----

# Create parents
Parents = newPop(founderPop)

# Make 100 crosses
F1 = randCross(Parents, nCrosses = 100)

# Self to create populations of 20 F2s per cross
F2 = self(F1, nProgeny = 20)

# Self one more time to have 20 F3s per cross
F3 = self(F2)

# Assume the F3s are used to derive lines and increased for yield trials 
# I'll use the automatically generated phenotypes for the F3s to represent 
# the values coming from these yield trials



## Fit genomic prediction models ----


## Three univariate models using the built-in solver
uniTrait1 = RRBLUP(F3, traits = 1)
uniTrait2 = RRBLUP(F3, traits = 2)
uniTrait3 = RRBLUP(F3, traits = 3)

# Assign BLUPs to population, giving a new name to keep it separate from other models
# This requires using setEBV once for each model. Note that the order you apply the 
# models needs to match the order of the traits and you have to use append=TRUE.
# When append=FALSE (default), any existing EBVs are deleted.
E1_uni = setEBV(pop = F3, 
                solution = uniTrait1)
E1_uni = setEBV(pop = E1_uni, 
                solution = uniTrait2, 
                append = TRUE)
E1_uni = setEBV(pop = E1_uni, 
                solution = uniTrait3, 
                append = TRUE)

# Measure accuracy
diag(cor(gv(E1_uni), ebv(E1_uni)))

# Make selections
E1_uni_sel = selectInd(E1_uni, 
                       nInd = 100, 
                       trait = selIndex, 
                       use = "ebv", 
                       b = b)

# Measure selection response
meanG(E1_uni_sel) - meanG(E1_uni)

# Economic response
(meanG(E1_uni_sel) - meanG(E1_uni)) %*% b



## One multivariate model using the built-in solver
# This model can be very slow with large data
mtASR = RRBLUP(F3, traits = 1:3)

# Assign BLUPs
E1_mtASR = setEBV(F3, solution = mtASR)

# Measure accuracy
diag(cor(gv(E1_mtASR), ebv(E1_mtASR)))

# Make selections
E1_mtASR_sel = selectInd(E1_mtASR, 
                         nInd = 100, 
                         trait = selIndex, 
                         use = "ebv", 
                         b = b)

# Measure selection response
meanG(E1_mtASR_sel) - meanG(E1_mtASR)

# Economic response
(meanG(E1_mtASR_sel) - meanG(E1_mtASR)) %*% b


## Using the external, PEGS based multivariate solver from bGWR package
# This model scales much better to large data than AlphaSimR's internal model
# at the cost of only slightly less accurate estimation of variance components
library(bWGR)

# Extract phenotype and genotype data
Y = pheno(F3)
M = pullSnpGeno(F3)

# Fit the model
mtBWGR = mrr(Y=Y, X=M)

# Assign BLUPs to population
# Note that to do this correctly I should be centering the markers, but neglecting 
# this step only shifts all BLUPs by a constant without any adverse effects on selection
E1_mtBWGR = F3
E1_mtBWGR@ebv = M %*% mtBWGR$b

# Measure accuracy
diag(cor(gv(E1_mtBWGR), ebv(E1_mtBWGR)))

# Make selections
E1_mtBWGR_sel = selectInd(E1_mtBWGR, 
                          nInd = 100, 
                          trait = selIndex, 
                          use = "ebv", 
                          b = b)

# Measure selection response
meanG(E1_mtBWGR_sel) - meanG(E1_mtBWGR)

# Economic response
(meanG(E1_mtBWGR_sel) - meanG(E1_mtBWGR)) %*% b

[cecilia249]

Thank you very much for all the information @gaynorr

I will use this with my data then and see if this code can solve the problem I'm been having.
Kind regards