Updates on genomic selection

Day 2: Genomic selection under the genotype-by-environment interaction

Data Science for Plant Breeding Lab

Piracicaba, September 2^nd

Objectives

To understand the factors that influence the performance of genomic selection models

To show how the genotype-by-environment interaction can be incorporated into genomic selection models

To run multi-environment models and validate them in different scenarios

On the accuracy of the genomic prediction

Genomic selection

Wallace et al. (2018)

Genomic selection

Hickey et al. (2017)

Genomic selection

Gaynor et al. (2017)

Genomic selection

Werner et al. (2023)

Genomic selection

Lenz et al. (2019)

What impacts model performance?

Following Grattapaglia (2022):

Theory from Population and Quantitative Genetics

Effective population size \(\left(N_e\right)\)
Genotyping density (marker-QTL LD)
Size and composition of training populations
Trait heritability
Genetic architecture of the target trait

Practical factors

Phenotyping quality
Genotyping platform (quality and cost)
Genotype-by-environment interaction
Model recycling
Adaptation of the breeding pipeline
Software

What impacts model performance?

Effective population size

Guilherme Pereira

Size and composition of the training population

Guilherme Pereira

What impacts model performance?

Size and composition of the training population

Couto et al. (2024)

What impacts model performance?

Size and composition of the training population

Fernández-Gonzaléz et al. (2023)

What impacts model performance?

Trait heritability

Gemenet et al. (2020)

What impacts model performance?

Alemu et al. (2024)

Genotype-by-environment interaction

Genotype-by-environment interaction

Aubé & Landry (2024)

GEI is a multi-dimensional phenomenon
Association with fitness
- An individual’s ability to survive and leave descendants
Dependent on trait and germplasm

Multi-environment trials

Bančič et al. (2024)

Trials laid-out in multiple locations/years/seasons
Different managements
Sample of GEI effects
Integral part of breeding programs
Value for Cultivation and Use (VCU) trials

Multi-environment trials

Historical data

Krause et al. (2023)

Stability

GEI and GS

A model trained in one environment may not work in another. An “environment” can be a location, a season, a year or even a different management.

Certainly, different environmental conditions
Probably, different genotypes

\[ \Delta G = \frac{i \times \color{red}{r_A} \times \sigma_A}{L \times C} \]

GEI and GS

Real life: unbalanced data

Show code

library(ggpubr)
library(tidyverse)

dat = readRDS("../Data/G2F.RDS")

geno_per_trial = as.data.frame(table(dat$gen, dat$loc)) |> 
  mutate(Freq = ifelse(Freq != 0, 1, 0))
geno_per_year = as.data.frame(table(dat$gen, dat$year)) |> 
  mutate(Freq = ifelse(Freq != 0, 1, 0))


ggarrange(
  ggplot(data = geno_per_year, aes(
    x = Var2, y = Var1, fill = as.factor(Freq)
  )) +
    geom_tile() +
    theme_minimal() +
    theme(
      axis.text.y = element_blank(),
      legend.title = element_blank(),
      legend.position = 'top'
    ) +
    scale_fill_manual(
      values = c('#b2bc63', '#10342d'),
      labels = c("Abscent", "Present")
    ) +
    scale_alpha_manual(
      values = c(.4, 1),
      labels = c("Abscent", "Present")
    ) +
    
    labs(x = "Year", y = "Hybrid"),
  
  ggplot(data = geno_per_trial, aes(
    x = Var2, y = Var1, fill = as.factor(Freq)
  )) +
    geom_tile() +
    theme_minimal() +
    theme(
      axis.text.y = element_blank(),
      legend.title = element_blank(),
      legend.position = 'top',
      axis.text.x = element_text(
        angle = 90,
        vjust = .5,
        hjust = 1
      )
    ) +
    scale_fill_manual(
      values = c('#b2bc63', '#10342d'),
      labels = c("Abscent", "Present")
    ) +
    scale_alpha_manual(
      values = c(.5, 1),
      labels = c("Abscent", "Present")
    ) +
    
    labs(x = "Location", y = "Hybrid"),
  common.legend = TRUE,
  labels = "AUTO",
  nrow = 1
)

GEI and GS

Sparse-testing

Jarquin et al. (2020)

GEI and GS

Sparse-testing

Jarquin et al. (2020)

GEI and GS

Sparse-testing

Atanda et al. (2022)

GEI and GS

Leveraging historical data

Skovbjerg et al. (2025)

GEI and GS

Leveraging historical data

Skovbjerg et al. (2025)

GEI and GS

Leveraging historical data

Atanda et al. (2021)

GEI and GS

Defining breeding zones

Chen et al. (2026)

Multi-environment genomic prediction

Dimensionality

Marker-based models are harder to fit
Dimensionality increase
More appropriate (with equivalent results): Genotype-based models
GBLUP (or RKHS, in BGLR)

Multi-environment GBLUP model:

\[ \boldsymbol y = \boldsymbol X \boldsymbol e + \boldsymbol Z_1 \boldsymbol g + \boldsymbol Z_2 \boldsymbol {ge} + \boldsymbol \varepsilon \]

GBLUP

\[ \begin{bmatrix} \boldsymbol{X}^\prime \boldsymbol X & \boldsymbol X^\prime \boldsymbol Z_1 & \boldsymbol X^\prime \boldsymbol Z_2 \\ \boldsymbol Z_1^\prime \boldsymbol X & \left(\boldsymbol Z_1^\prime \boldsymbol Z_1 + \lambda_1 \boldsymbol{G}^{-1} \right) & \boldsymbol Z_1^\prime\boldsymbol Z_2 \\ \boldsymbol Z_2^\prime \boldsymbol X & \boldsymbol Z^\prime_2\boldsymbol Z_1 & \left(\boldsymbol Z_2^\prime \boldsymbol Z_2 + \lambda_2 \boldsymbol{G}^{-1} \right) \end{bmatrix} \begin{bmatrix} \boldsymbol e \\ \boldsymbol g \\ \boldsymbol {ge} \end{bmatrix} = \begin{bmatrix} \boldsymbol X^\prime \boldsymbol y \\ \boldsymbol Z_1^\prime \boldsymbol y \\ \boldsymbol Z_2^\prime \boldsymbol y \end{bmatrix} \]

\(\lambda_1\) and \(\lambda_2 \to\) Shrinkage factors
\(\boldsymbol G \to\) Relationship matrix

Relationship matrices

Numerator relationship matrix

Derived from the pedigree

Diagonal: \(1 + F_i\) (inbreeding coefficient)

Off-diagonal: Numerator of the coefficient of relationship:

\[ r_{ii^\prime} = \frac{a_{ii^\prime}}{\sqrt{\left(1 + F_i\right)\left(1 + F_{i^\prime}\right)}} \]

Twice the coefficient of coancestry (The probability that an allele in X is IBD to an allele in the same locus in Y)

\[ a_{ii^\prime} \begin{cases} \mbox{Half-siblings} \to \frac{1}{4} \\ \mbox{Full siblings} \to \frac{1}{2} \\ \mbox{Parent-offspring} \to \frac{1}{2} \end{cases} \]

Relationship matrices

Numerator relationship matrix: Average (expected) relationship

Relationship matrices

Genomic relationship matrix

Derived from molecular data

Diagonal: Homozygosity in relation to the reference population

Off-diagonal: Realized relationship (shared alleles) between individuals \(i\) and \(i^\prime \to\) Mendelian sampling

VanRaden (2008):

\[ \boldsymbol{G} = \frac{\boldsymbol{\dot{W}}\boldsymbol{\dot{W}}^\prime}{2\sum_m p_m \left(1-p_m\right)} \]

Relationship matrices

Genomic relationship matrix: Realized relationship

Multi-environment GS model

Reaction norm model via RKHS (Reproducing Kernel Hilbert Space)

\[ \bar{\boldsymbol y} = \boldsymbol 1 \mu + \boldsymbol X \boldsymbol e + \boldsymbol Z_1 \boldsymbol g + \boldsymbol Z_2 \boldsymbol {ge} + \boldsymbol \varepsilon \]

\(\boldsymbol e \to\) Fixed effect of environment
\(\boldsymbol g \to\) Random main genetic effect
- \(\boldsymbol g \sim \mathcal N \left(\boldsymbol 0, \boldsymbol K_g \sigma^2_g\right)\)
\(\boldsymbol {ge} \to\) Random effect of the genotype-by-environment interaction
- \(\boldsymbol {ge} \sim \mathcal N \left(\boldsymbol 0, \boldsymbol K_{ge} \sigma^2_{ge}\right)\)
\(\boldsymbol \varepsilon \to\) IID residual

Cross-validations

CV1: Ability to predict untested individuals. The genetic covariance between relatives in the sole source of information (more challenging)

Cross-validations

CV2: Sparse-testing scenario. There are two sources of information in this case: Covariance between relatives and the individual’s performance in tested environments

Concluding remarks

Genomic selection is strongly influenced by practical and theoretical factors
Placing the genomic selection into a established breeding pipeline involves considering all factors above
The GEI is a rule in plant breeding, and must not be ignored
Breaking the TPE into zones where GEI is small may be the way to go
GBLUP models are flexible enough to incorporate multi-level information and provide multi-environment predictions

Practical session

The dataset

G2F

Phenotypic data

library(tidyverse)

dat = readRDS("../Data/G2F.RDS")
str(dat)

'data.frame':   85949 obs. of  7 variables:
 $ gen      : chr  "B37/C103" "B37/H95" "B73/MO17" "B73/PHG35" ...
 $ year     : int  2014 2014 2014 2014 2014 2014 2014 2014 2014 2014 ...
 $ loc      : chr  "GAH1" "GAH1" "GAH1" "GAH1" ...
 $ env      : chr  "GAH1_2014" "GAH1_2014" "GAH1_2014" "GAH1_2014" ...
 $ blue     : num  4.69 5.98 6.29 9.73 8.05 ...
 $ std.error: num  0.828 0.828 0.828 0.828 0.828 ...
 $ wt       : num  1.91 1.91 1.91 1.91 1.91 ...

Phenotypic data

Show code

ggplot(data = as.data.frame(table(dat$gen, dat$env)),
       aes(x = Var2, y = Var1, fill = as.factor(Freq))) + 
  geom_tile(aes(alpha = as.factor(Freq))) + 
  theme_minimal() + 
  theme(axis.text = element_blank(), legend.title = element_blank(), 
        legend.position = 'top') +
  scale_fill_manual(values = c('#b2bc63','#10342d'), labels = c("Abscent", "Present")) +
  scale_alpha_manual(values = c(.6, 1), labels = c("Abscent", "Present")) + 
  labs(x = "Environment (Year-Location combination)", y = "Hybrid")

Genomic data

geno.code = readRDS(file = "../Data/mrk_filtered.RDS")
str(geno.code)

 int [1:5899, 1:2105] 1 0 0 1 1 1 1 0 0 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:5899] "01CSI6/LH287" "01DIB2/LH287" "01DIB2/PHP02" "2369/DK3IIH6" ...
  ..$ : chr [1:2105] "S1_1007742" "S1_1020677" "S1_2018002" "S1_2275970" ...

any(is.na(geno.code))

[1] FALSE

Data management

Filtering for didactic purposes!

dat = droplevels(dat[which(dat$year > 2021 & dat$loc %in% c("GAH1","IAH1","IAH2","ILH1", "INH1", "MNH1", "MOH1", "MOH2")),])
dat = dat[which(dat$gen %in% rownames(geno.code)),]
set.seed(7)
dat = dat[which(dat$gen %in% sample(unique(dat$gen),150)),]
geno.code = geno.code[which(rownames(geno.code) %in% dat$gen),]
all(unique(dat$gen) %in% rownames(geno.code))

[1] TRUE

Data management

Working dataset

Show code

ggplot(data = dat, aes(
  x = reorder(env, blue),
  y = reorder(gen, blue),
  fill = blue
)) +
  geom_tile() +
  theme_minimal() +
  theme(axis.text.y = element_blank(), axis.text.x = element_text(angle = 90),
        legend.position = "top", legend.title = element_blank()) +
  scale_fill_viridis_c(option = "inferno") +
  labs(x = paste("Environments: ", length(unique(dat$env))),
       y = paste("Hybrid: ", length(unique(dat$gen))))

GRM

\[ \boldsymbol{G} = \frac{\boldsymbol{\dot{W}}\boldsymbol{\dot{W}}^\prime}{2\sum_m p_m \left(1-p_m\right)} \]

library(AGHmatrix)
Gmat = Gmatrix(SNPmatrix = geno.code, method = "VanRaden")

Initial data: 
    Number of Individuals: 150 
    Number of Markers: 2105 

Missing data check: 
    Total SNPs: 2105 
     0 SNPs dropped due to missing data threshold of 0.5 
    Total of: 2105  SNPs 

MAF check: 
    No SNPs with MAF below 0 

Heterozigosity data check: 
    No SNPs with heterozygosity, missing threshold of =  0 

Summary check: 
    Initial:  2105 SNPs 
    Final:  2105  SNPs ( 0  SNPs removed) 
 
Completed! Time = 0.34  seconds

GRM

Show code

library(ComplexHeatmap)
library(circlize)
col_fun = colorRamp2(c(-0.2, 0, 1.5), c("blue4", "white", "green4"))

Heatmap(
  Gmat,
  col = col_fun,
  show_row_names = F,
  show_column_names = F,
  heatmap_legend_param = list(title = "IBS")
)

GRM

library(plotly)

pc = prcomp(Gmat) 
pc.df = rbind(
  as.data.frame(pc$rotation) |> 
    rownames_to_column("ID") |> 
    mutate(comp = "load")
)

plot_ly(
  pc.df,
  x = ~ PC1,
  y = ~ PC2,
  z = ~ PC3
)  |> add_markers(
  text = ~ paste0("Hybrid: ", ID),
  hoverinfo = "text",
  marker = list(color = '#b2bc63', opacity = .62)
) |>  layout(scene = list(
  xaxis = list(
    title = paste0("PC1 (", round(pc$sdev[1]^2/sum(pc$sdev^2)*100,2), "%)"),
    showgrid = TRUE,
    range = c(min(c(
      range(pc.df$pc1), range(pc.df$pc2), range(pc.df$pc3)
    )), max(c(
      range(pc.df$pc1), range(pc.df$pc2), range(pc.df$pc3)
    )))
  ),
  yaxis = list(
    title = paste0("PC2 (", round(pc$sdev[2]^2/sum(pc$sdev^2)*100,2), "%)"),
    showgrid = TRUE,
    range = c(min(c(
      range(pc.df$pc1), range(pc.df$pc2), range(pc.df$pc3)
    )), max(c(
      range(pc.df$pc1), range(pc.df$pc2), range(pc.df$pc3)
    )))
  ),
  zaxis = list(
    title = paste0("PC3 (", round(pc$sdev[3]^2/sum(pc$sdev^2)*100,2), "%)"),
    showgrid = TRUE,
    range = c(min(c(
      range(pc.df$pc1), range(pc.df$pc2), range(pc.df$pc3)
    )), max(c(
      range(pc.df$pc1), range(pc.df$pc2), range(pc.df$pc3)
    )))
  )
))

Multi-environment GS model

Reaction norm model via RKHS (Reproducing Kernel Hilbert Space)

\[ \bar{\boldsymbol y} = \boldsymbol 1 \mu + \boldsymbol Z_e \boldsymbol e + \boldsymbol Z_g \boldsymbol g + \boldsymbol Z_{ge} \boldsymbol {ge} + \boldsymbol \varepsilon \]

\(\boldsymbol e \to\) Fixed effect of environment
\(\boldsymbol g \to\) Random main genetic effect
- \(\boldsymbol g \sim \mathcal N \left(\boldsymbol 0, \boldsymbol K_g \sigma^2_g\right)\)
\(\boldsymbol {ge} \to\) Random effect of the genotype-by-environment interaction
- \(\boldsymbol {ge} \sim \mathcal N \left(\boldsymbol 0, \boldsymbol K_{ge} \sigma^2_{ge}\right)\)
\(\boldsymbol \varepsilon \to\) IID residual

Multi-environment GS model

Expanding the data frame

nrow(dat)

[1] 2027

grid = expand.grid(
  gen = unique(dat$gen),
  env = unique(dat$env),
  stringsAsFactors = FALSE
)
dat = merge(grid, dat, by = c("gen", "env"), all.x = TRUE)
nrow(dat)

[1] 2250

dat$gen = factor(dat$gen, levels = rownames(Gmat), ordered = TRUE)
dat = dat[order(dat$env, dat$gen),]

Multi-environment GS model

Expanding the data frame

Show code

ggplot(data = dat, aes(
  x = reorder(env, blue),
  y = reorder(gen, blue),
  fill = blue
)) +
  geom_tile() +
  theme_minimal() +
  theme(axis.text.y = element_blank(), axis.text.x = element_text(angle = 90),
        legend.position = "top", legend.title = element_blank(), 
        text = element_text(size = 18)) +
  scale_fill_viridis_c(option = "inferno", na.value = "lightgrey") +
  labs(x = paste("Environments: ", length(unique(dat$env))),
       y = paste("Hybrid: ", length(unique(dat$gen))))

Multi-environment GS model

Incidence matrices

ZG = model.matrix(~-1 + gen, data = dat)
rownames(ZG) = dat$gen
colnames(ZG) = gsub("gen", "", colnames(ZG))
ZG[1:5, 1:5]

                  B37/H95 B73/TX779 BGEM-0039-N/LH244 BGEM-0063-N/LH244
B37/H95                 1         0                 0                 0
B73/TX779               0         1                 0                 0
BGEM-0039-N/LH244       0         0                 1                 0
BGEM-0063-N/LH244       0         0                 0                 1
GEMN-0096/LH244         0         0                 0                 0
                  GEMN-0096/LH244
B37/H95                         0
B73/TX779                       0
BGEM-0039-N/LH244               0
BGEM-0063-N/LH244               0
GEMN-0096/LH244                 1

Multi-environment GS model

Incidence matrices

ZE = model.matrix(~-1 + env, data = dat)
rownames(ZE) = dat$env
colnames(ZE) = gsub("env", "", colnames(ZE))
ZE[1:5,1:5]

          GAH1_2022 GAH1_2023 IAH1_2022 IAH1_2023 IAH2_2022
GAH1_2022         1         0         0         0         0
GAH1_2022         1         0         0         0         0
GAH1_2022         1         0         0         0         0
GAH1_2022         1         0         0         0         0
GAH1_2022         1         0         0         0         0

Multi-environment GS model

\[ \boldsymbol K_g = \boldsymbol Z_g \boldsymbol G \boldsymbol Z_g^\prime \]

ZGZ = ZG %*% Gmat %*% t(ZG)
rownames(ZGZ) -> aux1
dim(ZGZ)

[1] 2250 2250

Multi-environment GS model

Expanded \(\boldsymbol G\) (for three environments)

Show code

aux = do.call(c, lapply(split(dat, dat$env), nrow))
submat = ZGZ[1:sum(aux[1:3]),1:sum(aux[1:3])]
col_fun = colorRamp2(c(min(submat), 0, max(submat)), c("blue4", "white", "green4"))

Heatmap(
  submat,
  show_row_dend = FALSE,
  show_row_names = FALSE,
  show_column_dend = FALSE,
  show_column_names = FALSE,
  column_split = rep(names(aux)[1:3], aux[1:3]),
  row_split = rep(names(aux)[1:3], aux[1:3]),
  col = col_fun,
  show_heatmap_legend = FALSE,
  border = TRUE,
  border_gp = gpar(col = "black", lty = 2)
)

Multi-environment GS model

\[ \boldsymbol K_e = \boldsymbol Z_e \boldsymbol Z_e^\prime \]

ZEZ = tcrossprod(ZE)
dim(ZEZ)

[1] 2250 2250

Multi-environment GS model

Only 1s and 0s!

Show code

submat = ZEZ[1:sum(aux[1:4]), 1:sum(aux[1:4])]
col_fun = colorRamp2(c(min(submat), 0, max(submat)), c("blue4", "white", "green4"))

Heatmap(
  submat,
  show_row_dend = FALSE,
  show_row_names = FALSE,
  show_column_dend = FALSE,
  show_column_names = FALSE,
  column_split = rep(names(aux)[1:4], aux[1:4]),
  row_split = rep(names(aux)[1:4], aux[1:4]),
  col = col_fun,
  show_heatmap_legend = FALSE,
  border = TRUE,
  border_gp = gpar(col = "black", lty = 2)
)

Multi-environment GS model

\[ \boldsymbol K_{ge} = \boldsymbol K_g \odot \boldsymbol K_e \]

GEI = ZGZ * ZEZ
aux3 = rownames(GEI)
dim(GEI)

[1] 2250 2250

Multi-environment GS model

Show code

submat = GEI[1:sum(aux[1:5]),1:sum(aux[1:5])]
col_fun = colorRamp2(c(min(submat), 0, max(submat)), c("blue4", "white", "green4"))

Heatmap(
  submat,
  show_row_dend = FALSE,
  show_row_names = FALSE,
  show_column_dend = FALSE,
  show_column_names = FALSE,
  column_split = rep(names(aux)[1:5], aux[1:5]),
  row_split = rep(names(aux)[1:5], aux[1:5]),
  col = col_fun,
  show_heatmap_legend = FALSE,
  border = TRUE,
  border_gp = gpar(col = "black", lty = 2),
  row_title_gp = gpar(fontsize = 8)
)

Multi-environment GS model

Eigendecomposition (efficient computation)

\[ \begin{gather} \boldsymbol K_{ge} = \boldsymbol V_{ge} \boldsymbol D_{ge} \boldsymbol V_{ge} \\ \boldsymbol K_{g} = \boldsymbol V_{g} \boldsymbol D_{g} \boldsymbol V_{g} \end{gather} \]

GEI = eigen(GEI)
rownames(GEI$vectors) = aux3
ZGZ = eigen(ZGZ)
rownames(ZGZ$vectors) = aux1

Multi-environment GS model

Lesser number of iterations required (than BayesB)

library(BGLR)
ETA = list(
  env = list(X = ZE, model = 'FIXED'),
  hyb = list(V = ZGZ$vectors, d = ZGZ$values, model = 'RKHS'),
  gei = list(V = GEI$vectors, d = GEI$values, model = 'RKHS')
)

model1 = BGLR(y = dat$blue, ETA = ETA, nIter = 12000, burnIn = 2000, 
              thin = 10, verbose = TRUE, saveAt = "complete_met1_")

Multi-environment GS model

Reparameterized model: Cholesky decomposition of \(\boldsymbol G\)

\[ \boldsymbol G = \boldsymbol L \boldsymbol L^\prime \]

L_G = t(chol(Gmat + diag(x=0.0001, nrow = nrow(Gmat))))

Show code

Heatmap(
  L_G,
  show_row_dend = FALSE,
  show_row_names = FALSE,
  show_column_dend = FALSE,
  show_column_names = FALSE,
  row_order = rownames(L_G),
  column_order = colnames(L_G),
  col = col_fun,
  show_heatmap_legend = FALSE
)

Multi-environment GS model

\[ \begin{gather} \boldsymbol X_g = \boldsymbol Z_g \times \boldsymbol L \\ \boldsymbol X_{ge} = \boldsymbol X_g:\boldsymbol Z_e \end{gather} \]

X_g = ZG %*% L_G
dim(X_g)

[1] 2250  150

X_ge = model.matrix(~ -1 + X_g:dat$env)
dim(X_ge)

[1] 2250 2250

Multi-environment GS model

\(\boldsymbol X_g\):

Show code

submat = X_g[1:sum(aux[1:2]),]
col_fun = colorRamp2(c(min(submat), 0, max(submat)), c("blue4", "white", "green4"))

Heatmap(
  submat,
  show_row_dend = FALSE,
  show_row_names = FALSE,
  show_column_dend = FALSE,
  show_column_names = FALSE,
  row_order = 1:nrow(submat),
  column_order = colnames(submat),
  col = col_fun,
  # row_split = rep(names(aux)[1:2], aux[1:2]),
  show_heatmap_legend = FALSE
)

Multi-environment GS model

\(\boldsymbol X_ge\):

Multi-environment GS model

ETA2 = list(
  env = list(X = ZE, model = 'FIXED'),
  hyb = list(X = X_g, model = 'BRR'),
  gei = list(X = X_ge, model = 'BRR')
)

model2 = BGLR(y = dat$blue, ETA = ETA2, nIter = 12000, burnIn = 2000, 
             thin = 10, verbose = TRUE, saveAt = "complete_met2_")

Multi-environment GS model

Model parameters

modparam = data.frame(varE = c(model1$varE, model2$varE),
                      varG = c(model1$ETA$hyb$varU, model2$ETA$hyb$varB),
                      varGE = c(model1$ETA$gei$varU, model2$ETA$gei$varB),
                      model = c("RKHS", "CHOL (BRR)"));modparam

      varE     varG    varGE      model
1 1.882919 1.553343 1.697239       RKHS
2 1.876174 1.724512 1.703228 CHOL (BRR)

Multi-environment GS model

Model parameters

Show code

modparam |> group_by(model) |> 
  mutate(varT = varE + varG + varGE,
         PvE = varE/varT,
         PvG = varG/varT,
         PvGE = varGE/varT) |> relocate(model, .before = varE) |> 
  pivot_longer(varE:PvGE) |> filter(!grepl("var", name)) |> 
  ggplot(aes(x = model, y = value)) +
  geom_col(aes(fill = name), color = "black") + 
  theme_bw() + 
  labs(x = "Model", y = "Percentage from total variance") +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_viridis_d(option="mako", labels = c("varE", "varG", "varGEI")) +
  theme(text = element_text(size = 18), legend.position = "top", 
        legend.title = element_blank())

Multi-environment GS model

Predictions

res = data.frame(gen = dat$gen, env = dat$env, blue = dat$blue,
                 pred1 = model1$yHat, pred2 = model2$yHat)

Show code

ggplot(data = res, aes(x = pred1, y = pred2)) + 
  geom_point(alpha = .6, color = '#10342d') +
  labs(x = "RKHS", y = "CHOL (BRR)",
       caption = paste("Correlation = ", round(cor(res$pred2, res$pred1), 4))) + 
  theme_bw() +
  theme(text = element_text(size = 18))

Multi-environment GS model

Convergence

library(coda)
trace_varE = read.table("../Saves/dat/complete_met2_varE.dat")[,1]
trace_varG = read.table("../Saves/dat/complete_met2_ETA_hyb_varU.dat")[,1]
trace_varGEI = read.table("../Saves/dat/complete_met2_ETA_gei_varU.dat")[,1]
varE_chain = mcmc(trace_varE[201:length(trace_varE)])
varG_chain = mcmc(trace_varG[201:length(trace_varG)])
varGEI_chain = mcmc(trace_varGEI[201:length(trace_varGEI)])
plot(varE_chain)

Multi-environment GS model

Convergence

effectiveSize(varE_chain)

    var1 
534.0472

effectiveSize(varG_chain)

var1 
1000

effectiveSize(varGEI_chain)

   var1 
229.046

geweke.diag(varE_chain)


Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5 

   var1 
-0.6276

geweke.diag(varG_chain)


Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5 

  var1 
-1.196

geweke.diag(varGEI_chain)


Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5 

  var1 
-0.852

Cross-validations

nfolds = 5
nrept = 5

Cross-validations

CV2: Sparse-testing

sets = list()
i = 1
repeat{
  set.seed(987 * i)
  sets[[i]] = sample(rep(1:nfolds, length.out = nrow(dat)))
  i = i + 1
  if(i > nrept) break
}
cvdata = lapply(sets, function(x){
  cvdata = dat
  cvdata$set = x
  return(cvdata)
})
for (i in 1:length(cvdata)) cvdata[[i]]$rept = i

Cross-validations

CV2: Sparse-testing

cv2 = lapply(cvdata, function(x){
  res.list = list()
  for (i in unique(x$set)) {
    
    yNA = x$blue
    yNA[x$set == i] = NA
    
    mod_cv = BGLR(y = yNA, ETA = ETA2, nIter = 8000, burnIn = 1000, 
                  thin = 10, verbose = FALSE, saveAt = "modCV")
    unlink(list.files(pattern = "modCV"))
    
    res.list[[i]] = data.frame(
      gen = x$gen,
      env = x$env,
      blue = x$blue,
      set = x$set,
      yNA = yNA,
      yhat = mod_cv$yHat
    ) |> filter(set == i)
    
  }
  res.list
})

cv2 = lapply(cv2, function(x) do.call(rbind, x))

Cross-validations

CV1: Untested genotypes

sets = list()
i = 1
repeat{
  set.seed(7 * i)
  sets[[i]] = sample(rep(1:nfolds, length.out = nlevels(dat$gen)))
  i = i + 1
  if(i > nrept) break
}
cvdata = lapply(sets, function(x){
  cvdata = dat
  aux = data.frame(gen = levels(dat$gen), set = x)
  cvdata$set = aux$set[match(cvdata$gen, aux$gen)]
  return(cvdata)
})
for (i in 1:length(cvdata)) cvdata[[i]]$rept = i

Cross-validations

CV1: Untested genotypes

cv1 = lapply(cvdata, function(x){
  res.list = list()
  for (i in unique(x$set)) {
    
    yNA = x$blue
    yNA[x$set == i] = NA
    
    mod_cv = BGLR(y = yNA, ETA = ETA2, nIter = 8000, burnIn = 1000, 
                  thin = 10, verbose = FALSE, saveAt = "modCV")
    unlink(list.files(pattern = "modCV"))
    
    res.list[[i]] = data.frame(
      gen = x$gen,
      env = x$env,
      blue = x$blue,
      set = x$set,
      yNA = yNA,
      yhat = mod_cv$yHat
    ) |> filter(set == i)
    
  }
  res.list
})
cv1 = lapply(cv1, function(x) do.call(rbind, x))

Cross-validations

res_cv1 = lapply(cv1, function(x){
  x |> reframe(rmspe = sqrt(mean((blue - yhat)^2, na.rm = T)),
               pa = cor(blue, yhat),
               .by = env)
})
res_cv1

[[1]]
         env    rmspe        pa
1  GAH1_2022 1.889566 0.1702812
2  IAH1_2022 2.022620 0.4196013
3  IAH2_2022 3.367848 0.3187834
4  ILH1_2022 1.874160 0.2719688
5  INH1_2022 1.423681 0.4246822
6  MNH1_2022 1.558176 0.4690019
7  MOH2_2022 1.650287 0.2446381
8  GAH1_2023 2.035345 0.2008071
9  IAH1_2023 1.836036 0.6045529
10 IAH2_2023 1.221760 0.6056442
11 ILH1_2023 1.404764 0.4389203
12 INH1_2023 1.892373 0.3624826
13 MNH1_2023 1.519305 0.5680988
14 MOH1_2023 1.388520 0.1090614
15 MOH2_2023 1.346523 0.3344103

[[2]]
         env    rmspe        pa
1  GAH1_2022 1.895644 0.1692424
2  IAH1_2022 2.022137 0.4211831
3  IAH2_2022 3.348890 0.3302070
4  ILH1_2022 1.802512 0.3632502
5  INH1_2022 1.467618 0.3716514
6  MNH1_2022 1.646410 0.3707182
7  MOH2_2022 1.605652 0.3188448
8  GAH1_2023 1.968891 0.2829812
9  IAH1_2023 1.820997 0.6173261
10 IAH2_2023 1.244461 0.5841689
11 ILH1_2023 1.412772 0.4303447
12 INH1_2023 1.892635 0.3607798
13 MNH1_2023 1.566753 0.5314768
14 MOH1_2023 1.360501 0.1643439
15 MOH2_2023 1.353044 0.3290723

[[3]]
         env    rmspe        pa
1  GAH1_2022 1.882118 0.1596959
2  IAH1_2022 1.960602 0.4768718
3  IAH2_2022 3.402141 0.2905243
4  ILH1_2022 1.841893 0.3115065
5  INH1_2022 1.523762 0.2920739
6  MNH1_2022 1.535121 0.4934392
7  MOH2_2022 1.655805 0.2327191
8  GAH1_2023 2.024429 0.2070284
9  IAH1_2023 1.835144 0.6052883
10 IAH2_2023 1.195456 0.6312664
11 ILH1_2023 1.336070 0.5263628
12 INH1_2023 1.837789 0.4195306
13 MNH1_2023 1.497860 0.5918761
14 MOH1_2023 1.348625 0.1581148
15 MOH2_2023 1.292592 0.4124761

[[4]]
         env    rmspe        pa
1  GAH1_2022 1.858716 0.2065210
2  IAH1_2022 2.048789 0.3957298
3  IAH2_2022 3.484799 0.2365388
4  ILH1_2022 1.891632 0.2530135
5  INH1_2022 1.533387 0.2818449
6  MNH1_2022 1.662090 0.3516716
7  MOH2_2022 1.665324 0.2225250
8  GAH1_2023 1.999532 0.2455527
9  IAH1_2023 1.924328 0.5472414
10 IAH2_2023 1.298657 0.5312214
11 ILH1_2023 1.454470 0.3739427
12 INH1_2023 1.965409 0.2801419
13 MNH1_2023 1.618175 0.4766807
14 MOH1_2023 1.400039 0.1025656
15 MOH2_2023 1.406617 0.2493828

[[5]]
         env    rmspe        pa
1  GAH1_2022 1.873404 0.1869904
2  IAH1_2022 2.031339 0.4127986
3  IAH2_2022 3.371061 0.3154021
4  ILH1_2022 1.794050 0.3727625
5  INH1_2022 1.453547 0.3866245
6  MNH1_2022 1.517166 0.5110332
7  MOH2_2022 1.632968 0.2727572
8  GAH1_2023 2.038913 0.1886679
9  IAH1_2023 1.831827 0.6083938
10 IAH2_2023 1.201414 0.6253156
11 ILH1_2023 1.418105 0.4232666
12 INH1_2023 1.858152 0.3986870
13 MNH1_2023 1.449322 0.6340561
14 MOH1_2023 1.319553 0.2265163
15 MOH2_2023 1.318570 0.3765377

Cross-validations

res_cv1_mean = do.call(rbind, lapply(res_cv1, function(x){
  reframe(x, rmspe = mean(rmspe), pa = mean(pa))
}))
res_cv1_mean$rept = 1:5
res_cv1_mean$CV = "CV1"
res_cv1_mean

     rmspe        pa rept  CV
1 1.762064 0.3695290    1 CV1
2 1.760595 0.3763727    2 CV1
3 1.744627 0.3872516    3 CV1
4 1.814131 0.3169716    4 CV1
5 1.740626 0.3959873    5 CV1

Cross-validations

res_cv2 = lapply(cv2, function(x){
  x |> reframe(rmspe = sqrt(mean((blue - yhat)^2, na.rm = T)),
               pa = cor(blue, yhat),
               .by = env)
})
res_cv2

[[1]]
         env    rmspe        pa
1  GAH1_2022 1.879606 0.2662937
2  IAH1_2022 1.797687 0.5931229
3  IAH2_2022 3.210433 0.4251466
4  ILH1_2022 1.725585 0.4573075
5  INH1_2022 1.177796 0.6625441
6  MNH1_2022 1.479689 0.5490231
7  MOH2_2022 1.530060 0.4449176
8  GAH1_2023 2.058755 0.2433487
9  IAH1_2023 1.539022 0.7588185
10 IAH2_2023 1.035859 0.7375820
11 ILH1_2023 1.244886 0.6049751
12 INH1_2023 1.613326 0.6118256
13 MNH1_2023 1.279705 0.7223321
14 MOH1_2023 1.256384 0.4034762
15 MOH2_2023 1.131159 0.6052655

[[2]]
         env    rmspe        pa
1  GAH1_2022 1.850584 0.2928406
2  IAH1_2022 1.737497 0.6325008
3  IAH2_2022 3.343739 0.3347174
4  ILH1_2022 1.740162 0.4444828
5  INH1_2022 1.202679 0.6439494
6  MNH1_2022 1.415093 0.5985186
7  MOH2_2022 1.591965 0.3748018
8  GAH1_2023 2.028339 0.2760836
9  IAH1_2023 1.528469 0.7627030
10 IAH2_2023 1.024830 0.7437002
11 ILH1_2023 1.208967 0.6339923
12 INH1_2023 1.607658 0.6158269
13 MNH1_2023 1.317938 0.7009847
14 MOH1_2023 1.250440 0.4028758
15 MOH2_2023 1.133991 0.6028341

[[3]]
         env    rmspe        pa
1  GAH1_2022 1.865127 0.2793878
2  IAH1_2022 1.765478 0.6163797
3  IAH2_2022 3.343894 0.3345007
4  ILH1_2022 1.700967 0.4806699
5  INH1_2022 1.166414 0.6718908
6  MNH1_2022 1.421570 0.5929087
7  MOH2_2022 1.518983 0.4505503
8  GAH1_2023 1.980359 0.3149610
9  IAH1_2023 1.534773 0.7588277
10 IAH2_2023 1.018471 0.7472612
11 ILH1_2023 1.238987 0.6119061
12 INH1_2023 1.645343 0.5850448
13 MNH1_2023 1.294980 0.7139840
14 MOH1_2023 1.244532 0.4129931
15 MOH2_2023 1.121777 0.6125460

[[4]]
         env    rmspe        pa
1  GAH1_2022 1.852358 0.2850538
2  IAH1_2022 1.768035 0.6125090
3  IAH2_2022 3.306263 0.3618793
4  ILH1_2022 1.671418 0.5061276
5  INH1_2022 1.191953 0.6521207
6  MNH1_2022 1.449059 0.5709491
7  MOH2_2022 1.545193 0.4249735
8  GAH1_2023 2.033443 0.2636737
9  IAH1_2023 1.556104 0.7556200
10 IAH2_2023 1.033698 0.7381824
11 ILH1_2023 1.223525 0.6227211
12 INH1_2023 1.524027 0.6711545
13 MNH1_2023 1.320506 0.7034818
14 MOH1_2023 1.275118 0.3930340
15 MOH2_2023 1.110774 0.6239383

[[5]]
         env    rmspe        pa
1  GAH1_2022 1.845270 0.2991219
2  IAH1_2022 1.793866 0.5942029
3  IAH2_2022 3.375256 0.3146511
4  ILH1_2022 1.730234 0.4577725
5  INH1_2022 1.170492 0.6687173
6  MNH1_2022 1.460144 0.5639850
7  MOH2_2022 1.545673 0.4246391
8  GAH1_2023 2.019055 0.2711392
9  IAH1_2023 1.564224 0.7458028
10 IAH2_2023 1.012800 0.7515672
11 ILH1_2023 1.257817 0.5936511
12 INH1_2023 1.566635 0.6440214
13 MNH1_2023 1.303418 0.7095823
14 MOH1_2023 1.251500 0.4103725
15 MOH2_2023 1.113108 0.6220204

Cross-validations

res_cv2_mean = do.call(rbind, lapply(res_cv2, function(x) reframe(x, rmspe = mean(rmspe), pa = mean(pa))))
res_cv2_mean$rept = 1:5
res_cv2_mean$CV = "CV2"
res_cv2_mean

     rmspe        pa rept  CV
1 1.597330 0.5390653    1 CV2
2 1.598823 0.5373875    2 CV2
3 1.590777 0.5455875    3 CV2
4 1.590765 0.5456946    4 CV2
5 1.600633 0.5380831    5 CV2

Cross-validations

Show code

rbind(res_cv1_mean, res_cv2_mean) |> 
  pivot_longer(rmspe:pa) |> 
  ggplot(aes(x = CV, y = value)) +
  facet_wrap(.~name, scales = "free_y", labeller = labeller(.cols = c("pa" = "Predictive ability", "rmspe" = "RMSPE")))+
  geom_boxplot(fill = "#b2bc63") +
  theme_bw() + 
  theme(axis.title.y = element_blank(), text = element_text(size = 18))

Cross-validations

Show code

for (i in 1:length(res_cv1)){
  res_cv1[[i]]$rept = i
  res_cv1[[i]]$CV = "CV1"
} 
for (i in 1:length(res_cv2)){
  res_cv2[[i]]$rept = i
  res_cv2[[i]]$CV = "CV2"
} 
rbind(
  do.call(rbind, res_cv1),
  do.call(rbind, res_cv2)
) |> pivot_longer(rmspe:pa) |>
  ggplot(aes(x = env, y = value)) +
  # facet_wrap(.~name, scales = "free_y", labeller = labeller(.cols = c("pa" = "Predictive ability", "rmspe" = "RMSPE"))) +
  facet_grid(name~CV, scales = 'free_y',
             labeller = labeller(.rows = c("pa" = "Predictive ability", 
                                           "rmspe" = "RMSPE")))+
  theme_bw() + 
  theme(axis.title.y = element_blank(), text = element_text(size = 14),
        axis.text.x = element_text(angle = 90, vjust = .5, hjust = 1)) +
  geom_boxplot(fill = "#b2bc63")

References

Couto, E. G. O., Chaves, S. F. S., Dias, K. O. G., Morales-Marroqu, J. A., Alves-Pereira, A., Colombo, C. A., & Zucchi, M. I. (2024). Training set optimization is a feasible alternative for perennial orphan crop domestication and germplasm management: an Acrocomia aculeata example. Frontiers in Plant Science, 15, 1441683. https://doi.org/10.3389/fpls.2024.1441683

de los Campos, G., Hickey, J. M., Pong-Wong, R., Daetwyler, H. D., & Calus, M. P. L. (2013). Whole-genome regression and prediction methods applied to plant and animal breeding. Genetics, 193(2), 327–345. https://doi.org/10.1534/genetics.112.143313

Fernández-González, J., Akdemir, D., & Isidro y Sánchez, J. (2023). A comparison of methods for training population optimization in genomic selection. Theoretical and Applied Genetics, 136(3), 30. https://doi.org/10.1007/s00122-023-04265-6

Grattapaglia, D. (2022). Twelve years into genomic selection in forest trees: climbing the slope of enlightenment of marker assisted tree breeding. Forests, 13(10), 1554. https://doi.org/10.3390/f13101554

Hickey, J. M., Chiurugwi, T., Mackay, I., & Powell, W. (2017). Genomic prediction unifies animal and plant breeding programs to form platforms for biological discovery. Nature Genetics, 49(9), 1297–1303. https://doi.org/10.1038/ng.3920

Jarquin, D., Crossa, J., Lacaze, X., Du Cheyron, P., Daucourt, J., Lorgeou, J., Piraux, F., Guerreiro, L., Pérez, P., Calus, M., Burgueño, J., & de los Campos, G. (2014). A reaction norm model for genomic selection using high-dimensional genomic and environmental data. Theoretical and applied genetics, 127(3), 595–607. https://doi.org/10.1007/s00122-013-2243-1

Jarquin, D., Howard, R., Crossa, J., Beyene, Y., Gowda, M., Martini, J. W. R., Covarrubias Pazaran, G., Burgueño, J., Pacheco, A., Grondona, M., Wimmer, V., & Prasanna, B. M. (2020). Genomic prediction enhanced sparse testing for multi-environment trials. G3 GenesGenomesGenetics, 10(8), 2725–2739. https://doi.org/10.1534/g3.120.401349

Misztal, I., Aguilar, I., Lourenco, D., Ma, L., Steibel, J. P., & Toro, M. (2021). Emerging issues in genomic selection. Journal of Animal Science, 99(6), skab092. https://doi.org/10.1093/jas/skab092

Pérez, P., & de los Campos, G. (2014). Genome-wide regression and prediction with the BGLR statistical package. Genetics, 198(2), 483–495. https://doi.org/10.1534/genetics.114.164442

VanRaden, P. M. (2008). Efficient methods to compute genomic predictions. Journal of Dairy Science, 91(11), 4414–4423. https://doi.org/10.3168/jds.2007-0980

Voss-Fels, K. P., Cooper, M., & Hayes, B. J. (2019). Accelerating crop genetic gains with genomic selection. Theoretical and Applied Genetics, 132(3), 669–686. https://doi.org/10.1007/s00122-018-3270-8

Wallace, J. G., Rodgers-Melnick, E., & Buckler, E. S. (2018). On the road to breeding 4.0: Unraveling the good, the bad, and the boring of crop quantitative genomics. Annual Review of Genetics, 52(1), 421–444. https://doi.org/10.1146/annurev-genet-120116-024846

Werner, C. R., Gaynor, R. C., Sargent, D. J., Lillo, A., Gorjanc, G., & Hickey, J. M. (2023). Genomic selection strategies for clonally propagated crops. Theoretical and Applied Genetics, 136(4), 74. https://doi.org/10.1007/s00122-023-04300-6