Updates on genomic selection
Day 3: Genomic selection models enriched with environmental data
Data Science for Plant Breeding Lab
Piracicaba, September 3rd
Genotype-by-environment interaction
Genotype-by-environment interaction
GEI from the environmental perspective
Environmics
An omics approach that uses the full spectrum of environmental variations to understand and predict G × E interaction — Fritsche-Neto et al. (2025)
Environmental covariates
Data extrinsic to the tested genotypes
Environmic marker
An index vector derived from environmental covariates (ECs) which reflects the potential expression of a phenotype in a given space and time — Marcatti & Resende (2026)
Environmic marker
Marker engineering
Environmics
A new paradigm
Environmics
Plant breeding triangle
Reaction norm models
Cross-validations
CV0: predict the performance of tested genotypes in untested environments
Cross-validations
CV00: predict the performance of untested genotypes in untested environments (most challenging scenario of them all!)
Cross-validations
Overview
The dataset
Environmental data
Environmental profile of each environment
Show code
ggplot(data = soil |> pivot_longer(colnames(soil[-c(1, 2)])),
aes(x = Env, y = value)) +
facet_wrap(. ~ name, scales = "free_y") +
geom_col(color = "black", fill = '#b2bc63') +
theme_bw() +
theme(axis.text.x = element_text(
angle = 90,
vjust = .5,
hjust = 1
), axis.title.y = element_blank()) +
labs(x = "Environments")
Environmental data
Environmental profile of environment “GAH1_2022”
Show code
weather$Date = lubridate::ymd(weather$Date)
ggplot(data = weather |> pivot_longer(colnames(weather[-c(1, 2)])) |>
filter(Env == "GAH1_2022"),
aes(x = Date, y = value)) +
facet_wrap(name ~ ., scales = "free_y") +
geom_point(color = '#10342d', alpha = .75) +
geom_line(aes(group = 1), color = '#10342d', alpha = .8) +
theme_bw() +
theme(text = element_text(size = 15), axis.title.y = element_blank())
Environmental markers
Overview of the matrix of environmental markers
# any(is.na(envdata))
envdata = scale(envdata) ## Very important!!!Show code
library(ComplexHeatmap)
library(circlize)
col_fun = colorRamp2(c(min(envdata), 0, max(envdata)), c("blue4", "white", "red4"))
Heatmap(matrix = envdata,
show_row_dend = FALSE, show_column_dend = FALSE,
column_order = colnames(envdata),
column_split = ifelse(grepl("\\.", colnames(envdata)), "SOIL", gsub("_.*", "", colnames(envdata))),
show_column_names = FALSE,
column_title_rot = 45,
border = TRUE,
border_gp = gpar(col = "black", lty = 2),
col = col_fun,
show_heatmap_legend = FALSE)
Environmental markers
Environmental relationship matrix
Show code
col_fun = colorRamp2(c(min(Emat), 0, max(Emat)), c("blue4", "white", "green4"))
Heatmap(matrix = Emat,
col = col_fun,
show_heatmap_legend = TRUE,
heatmap_legend_param = list(title = "Relationship"))
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")
)
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
Show code
aux = do.call(c, lapply(split(dat, dat$env), nrow))
submat = ZOZ[1:sum(aux[1:6]), 1:sum(aux[1:6])]
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:6], aux[1:6]),
row_split = rep(names(aux)[1:6], aux[1:6]),
col = col_fun,
show_heatmap_legend = FALSE,
border = TRUE,
border_gp = gpar(col = "black", lty = 2),
row_title_gp = gpar(fontsize = 10)
)
Multi-environment GS model
Show code
submat = GEI_EC[1:sum(aux[1:5]),1:sum(aux[1:5])]
col_fun = colorRamp2(c(min(submat), 0, 3.8), 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)
)
Multi-environment GS model
Factorized matrix of environmental relationship
Show code
Heatmap(
L_E,
show_row_dend = FALSE,
show_row_names = FALSE,
show_column_dend = FALSE,
show_column_names = FALSE,
row_order = rownames(L_E),
column_order = colnames(L_E),
col = col_fun,
show_heatmap_legend = FALSE,
row_split = rownames(L_E),
row_title_gp = gpar(fontsize = 5)
)
Multi-environment GS model
Factorized covariance matrix of environmental covariates
Show code
Heatmap(
X_ec,
show_row_dend = FALSE,
show_row_names = FALSE,
show_column_dend = FALSE,
show_column_names = FALSE,
row_order = rownames(X_ec),
column_order = colnames(X_ec),
row_split = rep(names(aux), aux),
col = col_fun,
show_heatmap_legend = FALSE,
row_title_gp = gpar(fontsize = 5)
)
Multi-environment GS model
Factorized genotype-by-environmental covariate kernel
Show code
submat = X_ge_ec
rownames(submat) = colnames(submat) = rep(names(aux), aux)
col_fun = colorRamp2(c(min(submat), 0, max(submat)), c("blue4", "white", "green4"))
Heatmap(
X_ge_ec,
show_row_dend = FALSE,
show_row_names = FALSE,
show_column_dend = FALSE,
show_column_names = FALSE,
row_order = 1:nrow(X_ge_ec),
column_order = 1:ncol(X_ge_ec),
col = col_fun,
show_heatmap_legend = FALSE,
# row_split = rep(names(aux), aux),
# row_title_gp = gpar(fontsize = 5)
)
Multi-environment GS model
Ellapsed time
Show code
ellapsed = data.frame(
model = rep(c("GxEC", "GxEC + GxE"), each = 2),
mechanism = rep(c("RKHS", "CHOL"), times = 2),
ellapsed = c(model1$ellapsed, model2$ellapsed, model3$ellapsed, model4$ellapsed)
)
ggplot(data = ellapsed, aes(x = model, y = ellapsed, fill = mechanism)) +
geom_col(position = position_dodge(), color = "black") +
theme_bw() +
theme(text = element_text(size = 18), legend.position = 'top',
legend.title = element_blank()) +
scale_fill_manual(values = c('#b2bc63','#10342d')) +
labs(x = "Model", y = "Ellapsed time (min)")
Multi-environment GS model
Model parameters
Show code
modparam |> mutate(control = paste(model, mechanism, sep = "_")) %>%
group_by(control) |>
mutate(varT = ifelse(grepl("GxE_", control),residual + ec + g + gxec + gxe,
residual + ec + g + gxec),
PvR = residual/varT,
PvEC = ec/varT,
PvG = g/varT,
PvGxEC = gxec/varT,
PvGxE = gxe/varT) |>
pivot_longer(PvR:PvGxE) |>
ggplot(aes(x = model, y = value)) +
facet_wrap(.~mechanism) +
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('PvR' = "varE", 'PvG' = "varG", 'PvGxE' = "varGEI",
"PvGxEC" = "varE", "PvEC" = "varEC")) +
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,
'GxEC_RKHS' = model1$yHat, 'GxEC_CHOL' = model2$yHat,
'GxE_GxEC_RKHS' = model3$yHat, 'GxE_GxEC_CHOL' = model4$yHat)Show code
GGally::ggpairs(res[,-c(1:3)]) + theme_bw() +
theme(text = element_text(size = 15),
strip.text.y.right = element_text(size = 10))
Cross-validations
CV0: prediction of untested environments
Cross-validations
CV00: prediction of untested genotypes in untested environments
Cross-validations
Show code
rbind(res_cv1_mean, res_cv2_mean, res_cv00_mean) |>
ggplot(aes(x = CV, y = value, fill = model)) +
facet_wrap(.~metric, scales = "free_y", labeller = labeller(.cols = c("pa" = "Predictive ability", "rmspe" = "RMSPE"))) +
geom_boxplot() +
theme_bw() +
theme(axis.title.y = element_blank(), text = element_text(size = 18),
legend.position = 'top') +
labs(x = "Cross-validation scenario", fill = "Model") +
scale_fill_manual("Model", values = c("#b2bc63", "#10342d"),
labels = c("GxE + GxEC", "GxE")) +
scale_color_manual(values = c("#b2bc63", "#10342d"),
labels = c("GxE + GxEC", "GxE"), guide = 'none') +
geom_linerange(
data = res_cv0_mean, aes(x = CV, ymin = 0, ymax = value, color = model),
position = position_dodge(width = 0.6),
show.legend = FALSE,
linewidth = 1
) +
geom_point(
data = res_cv0_mean, aes(x = CV, y = value, color = model),
position = position_dodge(width = 0.6),
size = 2.5,
show.legend = FALSE,
shape = 17
)
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"
}
for (i in 1:length(res_cv00)){
res_cv00[[i]]$rept = i
res_cv00[[i]]$CV = "CV00"
}
rbind(
do.call(rbind, res_cv1),
do.call(rbind, res_cv2),
do.call(rbind, res_cv00)
) |>
ggplot(aes(x = env, y = value, fill = model)) +
facet_grid(
metric ~ 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),
legend.position = "top"
) +
geom_boxplot() +
labs(x = "Environments", fill = "Model") +
scale_fill_manual(
values = c("#b2bc63", "#10342d"),
labels = c("GxE + GxEC", "GxE")
) +
scale_color_manual(
values = c("#b2bc63", "#10342d"),
labels = c("GxE + GxEC", "GxE"),
guide = 'none'
) +
geom_linerange(
data = res_cv0 |> mutate(CV = "CV0"),
aes(x = env, ymin = 0, ymax = value, color = model, group = model),
position = position_dodge(width = 0.6),
show.legend = FALSE,
linewidth = 1, alpha = .7
) +
geom_point(
data = res_cv0 |> mutate(CV = "CV0"),
aes(x = env, y = value, color = model, group = model),
position = position_dodge(width = 0.6),
size = 2.5,
show.legend = FALSE,
shape = 17,
alpha = .7
)