Genetic algorithm using the evola package

2a) Best parents for the next generation

One situation that occurs in plant and animal breeding is the so called ‘optimal contribution’ problem where we want to pick a set of parents that can maximize the gain while managing genetic variance as much as possible. In the following example we take a population of 363 possible parents (which will become the genes) and pick the best 20 while conserving genetic variance (group relationship).

data(DT_cpdata)
DT <- DT_cpdata
head(DT)

##        id Row Col Year      color  Yield FruitAver Firmness Rowf Colf occ
## P003 P003   3   1 2014 0.10075269 154.67     41.93  588.917    3    1   0
## P004 P004   4   1 2014 0.13891940 186.77     58.79  640.031    4    1   1
## P005 P005   5   1 2014 0.08681502  80.21     48.16  671.523    5    1   1
## P006 P006   6   1 2014 0.13408561 202.96     48.24  687.172    6    1   1
## P007 P007   7   1 2014 0.13519278 174.74     45.83  601.322    7    1   1
## P008 P008   8   1 2014 0.17406685 194.16     44.63  656.379    8    1   1

Our surrogate of fitness will be the Yield trait and we will have a second trait to control the number of individuals we can select. We will set a constraint for the occurrence (occ) trait to 20 but the only trait contributing to fitness will be Yield (using the traitWeight argument).

# get best 20 individuals weighting variance by 0.5
res<-evolafit(cbind(Yield, occ)~id, dt= DT, 
              # constraints: if sum is greater than this ignore 
              constraintsUB = c(Inf,20), 
              # constraints: if sum is smaller than this ignore
              constraintsLB= c(-Inf,-Inf), 
              # weight the traits for the selection
              traitWeight = c(1,0), 
              # population parameters
              nCrosses = 100, nProgeny = 10, 
              # coancestry parameters
              A=A, lambda= (30*pi)/180 , nQTLperInd = 2, 
              # selection parameters
              propSelBetween = 0.5, propSelWithin =0.5, 
              nGenerations = 20, verbose=FALSE) 

We then use the bestSol() function to extract the solution that maximized our fitness function and constraints.

best = bestSol(res)["pop","Yield"];
sum(res$M[best,]) # total # of inds selected

## [1] 20

We can use the pmonitor() function to see if convergence was achieved between the best and the average solutions.

pmonitor(res)
plot(DT$Yield, col=as.factor(res$M[best,]), 
     pch=(res$M[best,]*19)+1)

2b) Obtaining optimal N crosses from a population for a given trait/feature

A variation of the same problem is when we want to pick the best crosses instead of the best parents to directly find the optimal solution for a crossing block. In the following example we use a dataset of crosses with marker and phenotype information to show how to optimize this problem.

data(DT_technow)
DT <- DT_technow
DT$occ <- 1; DT$occ[1]=0
M <- M_technow
A <- A.mat(M)
head(DT)

##    hybrid dent flint     GY    GM      hy occ
## 1 518.298  518   298  -8.04 -0.85 518:298   0
## 2 518.305  518   305 -11.10  1.70 518:305   1
## 3 518.306  518   306 -16.85  2.24 518:306   1
## 4 518.316  518   316   2.08 -1.33 518:316   1
## 5 518.323  518   323   5.65 -2.71 518:323   1
## 6 518.327  518   327 -16.95 -0.52 518:327   1

The way to specify this problem is exactly the same than with the optimization of parents but the input information is at the level of predicted crosses instead of individuals (genes).

# run the genetic algorithm
  res<-evolafit(formula = c(GY, occ)~hy, dt= DT, 
                # constraints: if sum is greater than this ignore
                constraintsUB = c(Inf,100), 
                # constraints: if sum is smaller than this ignore
                constraintsLB= c(-Inf,-Inf),
                # weight the traits for the selection
                traitWeight = c(1,0), 
                # population parameters
                nCrosses = 100, nProgeny = 10, 
                # coancestry parameters
                A=A, lambda= (20*pi)/180 , nQTLperInd = 100, 
                # selection parameters
                propSelBetween = 0.5, propSelWithin =0.5, 
                nGenerations = 10, verbose=FALSE) 
best = bestSol(res)["pop","GY"]
sum(res$M[best,]) # total # of inds selected

## [1] 98

You can use the pmonitor() or pareto() functions to see the evolution of the solution and see the performance of the solution selected.

pmonitor(res)

plot(DT$GY, col=as.factor(res$M[best,]), 
       pch=(res$M[best,]*19)+1)

Notice that we have maximized te GY variable which is the pure phenotype, but alternatively you could use the GEBVs to get crosses with maximum GEBV or predict the TGV (total genetic value) to maximize the F1 performance (including the dominance).

3a) Optimizing a subsample of size N to be representative of its own

Since the objective is to select a set of 100 lines that represent best the training set (tp) of ~400 lines we will subset a relationship matrix for that training set (As).

As <- A[tp,tp]
DT2 <- DT[rownames(As),]

For this particular case there is no trait to optimize (x’a) but we just want to make sure that we maintain as much variation in the sample as possible (x’Ax). We then just create a dummy trait in the dataset (dummy) to put all the weight into the group relationship (x’Ax) using the lambda argument. The trait for occurrence we will use it as before to control the number of individuals to be in the sample.

res<-evolafit(cbind(dummy, occ)~id, dt= DT2, 
                # constraints: if sum is greater than this ignore 
                constraintsUB = c(Inf, 100), 
                # constraints: if sum is smaller than this ignore
                constraintsLB= c(-Inf, -Inf), 
                # weight the traits for the selection
                traitWeight = c(1,0), 
                # population parameters
                nCrosses = 100, nProgeny = 10, 
                # coancestry parameters
                A=As,
                lambda=(60*pi)/180, nQTLperInd = 80, 
                # selection parameters
                propSelBetween = 0.5, propSelWithin =0.5, 
                nGenerations = 15, verbose = FALSE)

best = bestSol(res)["pop","dummy"]
sum(res$M[best,]) # total # of inds selected

## [1] 100

You can see which individuals were selected.

cex <- rep(0.5,nrow(PCWheat))
names(cex) <- rownames(PCWheat)
cex[names(which(res$M[best,]==1))]=2
plot(PCWheat[,1], PCWheat[,2], col = TreeWheat, cex=cex,
     pch = TreeWheat, xlab = "pc1", ylab = "pc2")

Gender in breeding

In this case is better if you only create the cross combinations that are possible (e.g., where male and female can couple) and you handed them to the evolutionary algorithm. That means, the rows of the crosses to be in the searching space only include the realistic ones.

Number of times a parent should be used

In this case you can modify the fitness function to set to a low value the fitness of solutions that have used too many times the same parent. Using the DT_technow dataset this would be done the following way:

First you create an incidence matrix for parents in columns and hybrids in crosses:

data(DT_technow)
DT <- DT_technow
DT$occ <- 1; DT$occ[1]=0
M <- M_technow
A <- A.mat(M)

Z=with(DT,overlay(dent,flint) )#  Matrix::sparse.model.matrix(~dent-1, data=DT)
rownames(Z) <- DT$hy # needed to link to the QTL matrix

the secons step is to create a new fitness function for the genetic algorithm. Our objective function to be maximized is normally of the form Yb - d, where Y is the trait values, b is the trait weights, and d is the group relationship x’Ax. We then are going to put some additional constrait that parents of the crosses can’t show up more than twice. This can be done in the following way:

# regular fitness function
fitnessf <-function (Y, b, d, Q, Z) {
  fit <- Y %*% b - d
  return(fit)
}
# new fitness function with constraint
fitnessf <-function (Y, b, d, Q, Z) {
  X=Q%*%Z[colnames(Q),]
  bad <- as.vector( apply(X,1, function(x){length(which(x > 5))}) ) 
  bad <- which(bad > 0)
  fit <- Y %*% b - d
  if(length(bad) > 0){fit[bad,1]=min(fit[,1])}
  return(fit)
}

Notice that we have added a matrix product Q%*%Z to see how may times each parent is used in the proposed solution of crosses. The next step would be to provide the new fitness function to the evolafit() function and the additional argument Z which is the overlay matrix formed in the first chunck of code:

res<-evolafit(formula = c(GY, occ)~hy,
              dt= DT, 
              # constraints: if sum is greater than this ignore
              constraintsUB = c(Inf,50), 
              # constraints: if sum is smaller than this ignore
              constraintsLB= c(-Inf,-Inf),
              # weight the traits for the selection
              traitWeight = c(1,0), 
              # population parameters
              nCrosses = 100, nProgeny = 10, 
              # coancestry parameters
              A=A, lambda= (10*pi)/180 , nQTLperInd = 40, 
              # new fitness function and additional args
              fitnessf = fitnessf, Z=Z,
              # selection parameters
              propSelBetween = 0.5, propSelWithin =0.5, 
              nGenerations = 15, verbose=FALSE) 

best = bestSol(res)["pop","GY"]
xa = (res$M %*% DT$GY)[best,]; xa 

##  30051 
## 327.28

xAx = res$M[best,] %*% A %*% res$M[best,]; xAx 

##          [,1]
## [1,] 33.48015

sum(res$M[best,]) # total # of inds selected

## [1] 34

Now, last but not least we check how many times each parent was used:

# check how many times an individual was used in the final crosses
crosses <- data.frame(cross=names(which( res$M[best,] == 1)))
table(unlist(strsplit(crosses$cross,":")))

## 
## 1863 3205 3207 3208 3212 3214 3216 3217 3218  336  338  340  345  348  353  354 
##    1    1    1    2    1    1    1    2    1    1    1    2    2    1    1    1 
##  359 3606 3607 3608 3609 3613 3614 3618 3619  367  377  392  394  396  398  417 
##    1    1    2    2    1    1    2    1    1    2    1    1    1    2    1    1 
##  427  430  431  437  444  491  536  560  578  607  608  614  617  618  632  641 
##    1    1    1    1    1    1    1    1    2    1    1    1    1    1    1    1 
##  647  652  666  671  677  716  755  797 
##    1    1    1    2    2    1    1    1

# check performance of crosses selected
plot(DT$GY, col=as.factor(res$M[best,]), 
       pch=(res$M[best,]*19)+1)