Understanding photothermal interactions can help expand production range and increase genetic diversity of lentil (Lens culinaris Medik.)
Plants, People, Planet. 00:1-11.
Introduction
This vignette contains the code and analysis done for the paper:
https://github.com/derekmichaelwright/AGILE_LDP_Phenology
This work done as part of the AGILE project at the University of Saskatchewan.
Data Preparation
Load the nessesary R packages, Prepare the data for analysis.
#install.packages(c("tidyverse","scales","rworldmap","ggrepel","magick",
# "GGally","ggpubr,"ggbeeswarm","FactoMineR","plot3D","stringr",
# "plotly","leaflet","leaflet.minicharts","htmlwidgets"))
# Load libraries
library(tidyverse) # data wrangling
library(scales) # rescale()
library(rworldmap) # mapBubbles()
library(ggrepel) # geom_text_repel() + geom_label_repel()
library(magick) # image editing
library(GGally) # ggpairs() + ggmatrix()
library(ggpubr) # ggarrange()
library(ggbeeswarm) # geom_quasirandom()
library(FactoMineR) # PCA() & HCPC()
library(plot3D) # 3D plots
library(stringr) # str_pad()
library(plotly) # plot_ly()
library(leaflet) # leaflet()
library(leaflet.minicharts) # addMinicharts()
library(htmlwidgets) # saveWidget()
# General color palettes
<- c("darkred", "darkorange3", "darkgoldenrod2", "deeppink3",
colors "steelblue", "darkorchid4", "cornsilk4", "darkgreen")
# Expts color palette
<- c("lightgreen", "palegreen4", "darkgreen", "darkolivegreen3",
colors_Expt "darkolivegreen4", "springgreen4", "orangered2", "orangered4",
"palevioletred", "mediumvioletred", "orange2", "orange4",
"slateblue1", "slateblue4", "aquamarine3", "aquamarine4",
"deepskyblue3", "deepskyblue4" )
# Locations
<- c("Rosthern, Canada", "Sutherland, Canada", "Central Ferry, USA",
names_Location "Bhopal, India", "Jessore, Bangladesh", "Bardiya, Nepal",
"Cordoba, Spain", "Marchouch, Morocco", "Metaponto, Italy" )
# Experiments
<- c("Rosthern, Canada 2016", "Rosthern, Canada 2017",
names_Expt "Sutherland, Canada 2016", "Sutherland, Canada 2017",
"Sutherland, Canada 2018", "Central Ferry, USA 2018",
"Bhopal, India 2016", "Bhopal, India 2017",
"Jessore, Bangladesh 2016", "Jessore, Bangladesh 2017",
"Bardiya, Nepal 2016", "Bardiya, Nepal 2017",
"Cordoba, Spain 2016", "Cordoba, Spain 2017",
"Marchouch, Morocco 2016", "Marchouch, Morocco 2017",
"Metaponto, Italy 2016", "Metaponto, Italy 2017" )
# Experiment short names
<- c("Ro16", "Ro17", "Su16", "Su17", "Su18", "Us18",
names_ExptShort "In16", "In17", "Ba16", "Ba17", "Ne16", "Ne17",
"Sp16", "Sp17", "Mo16", "Mo17", "It16", "It17" )
# Macro-Environments
<- c("Temperate", "South Asia", "Mediterranean")
names_MacroEnvs # Lentil Diversity Panel metadata
<- read.csv("data/data_ldp.csv")
ldp # Country info
<- read.csv("data/data_countries.csv") %>% filter(Country %in% ldp$Origin)
ct # ggplot theme
<- theme_bw() +
theme_AGL theme(strip.background = element_rect(colour = "black", fill = NA, size = 0.5),
panel.background = element_rect(colour = "black", fill = NA, size = 0.5),
panel.border = element_rect(colour = "black", size = 0.5),
panel.grid = element_line(color = alpha("black", 0.1), size = 0.5),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())
# Create scaling function
<- function(x, trait) {
traitScale <- rep(NA, nrow(x))
xout for(i in unique(x$Expt)) {
<- x %>% filter(Expt == i) %>% pull(trait) %>% min(na.rm = T)
mn <- x %>% filter(Expt == i) %>% pull(trait) %>% max(na.rm = T)
mx <- ifelse(x$Expt == i, rescale(x %>% pull(trait), c(1,5), c(mn,mx)), xout)
xout
}
xout
}# Prep data
# Notes:
# - DTF2 = non-flowering genotypes <- group_by(Expt) %>% max(DTF)
<- read.csv("data/data_raw.csv") %>%
rr mutate(Rep = factor(Rep),
Year = factor(Year),
PlantingDate = as.Date(PlantingDate),
Location = factor(Location, levels = names_Location),
Expt = factor(Expt, levels = names_Expt),
ExptShort = plyr::mapvalues(Expt, names_Expt, names_ExptShort),
ExptShort = factor(ExptShort, levels = names_ExptShort),
DTF2_scaled = traitScale(., "DTF2"),
RDTF = round(1 / DTF2, 6),
VEG = DTF - DTE,
REP = DTM - DTF)
# Average raw data
<- rr %>%
dd group_by(Entry, Name, Expt, ExptShort, Location, Year) %>%
summarise_at(vars(DTE, DTF, DTS, DTM, VEG, REP, RDTF, DTF2),
funs(mean), na.rm = T) %>%
ungroup() %>%
mutate(DTF2_scaled = traitScale(., "DTF2"))
# Prep environmental data
<- read.csv("data/data_env.csv") %>%
ee mutate(Date = as.Date(Date),
ExptShort = plyr::mapvalues(Expt, names_Expt, names_ExptShort),
ExptShort = factor(ExptShort, levels = names_ExptShort),
Expt = factor(Expt, levels = names_Expt),
Location = factor(Location, levels = names_Location),
MacroEnv = factor(MacroEnv, levels = names_MacroEnvs),
DayLength_rescaled = rescale(DayLength, to = c(0, 40)) )
# Prep field trial info
<- dd %>% group_by(Expt) %>%
xx summarise_at(vars(DTE, DTF, DTS, DTM), funs(min, mean, max), na.rm = T) %>%
ungroup()
<- read.csv("data/data_info.csv") %>%
ff mutate(Start = as.Date(Start) ) %>%
left_join(xx, by = "Expt")
for(i in unique(ee$Expt)) {
<- ee %>%
ee filter(Expt != i | (Expt == i & DaysAfterPlanting <= ff$DTM_max[ff$Expt == i]))
}<- ee
xx for(i in unique(ee$Expt)) {
<- xx %>%
xx filter(Expt != i | (Expt == i & DaysAfterPlanting <= ff$DTF_max[ff$Expt == i]))
} <- xx %>%
xx group_by(Location, Year) %>%
summarise(T_mean = mean(Temp_mean, na.rm = T), T_sd = sd(Temp_mean, na.rm = T),
P_mean = mean(DayLength, na.rm = T), P_sd = sd(DayLength, na.rm = T) ) %>%
ungroup() %>%
mutate(Expt = paste(Location, Year)) %>%
select(-Location, -Year)
<- ff %>%
ff left_join(xx, by = "Expt") %>%
mutate(ExptShort = plyr::mapvalues(Expt, names_Expt, names_ExptShort),
ExptShort = factor(ExptShort, levels = names_ExptShort),
Expt = factor(Expt, levels = names_Expt),
Location = factor(Location, levels = names_Location),
MacroEnv = factor(MacroEnv, levels = names_MacroEnvs),
T_mean = round(T_mean, 1),
P_mean = round(P_mean, 1))
ldp
= Lentil Diversity Panel Metadatarr
= Raw Phenotype Datadd
= Averaged Phenotype Dataee
= Environmental Dataff
= Field Trial Infoct
= Country Info
Materials & Methods
Supplemental Table 1: LDP
Table S1: Genotype entry number, name, common synonyms, origin and source of lentil genotypes used in this study. These genotypes are gathered from the University of Saskatchewan (USASK), Plant Gene Resources of Canada (PGRC), United States Department of Agriculture (USDA), International Center for Agricultural Research in the Dry Areas (ICARDA).
<- select(ldp, Entry, Name, Origin, Source, Synonyms)
s1 write.csv(s1, "Supplemental_Table_01.csv", row.names = F)
Supplemental Table 2: Field Trial Info
Table S2: Details of the field trials used in this study, including location information, planting dates, mean temperature and photoperiods and details on plot type and number of seeds sown.
<- ff %>%
s2 select(Location, Year, `Short Name`=ExptShort, Latitude=Lat, Longitude=Lon,
`Planting Date`=Start, `Temperature (mean)`=T_mean, `Photoperiod (mean)`=P_mean,
`Number of Seeds Sown`=Number_of_Seeds_Sown, `Plot Type`=Plot_Type)
write.csv(s2, "Supplemental_Table_02.csv", row.names = F)
<- ee %>%
xx select(ExptShort, Temp_min, Temp_max, DayLength) %>%
group_by(ExptShort) %>%
summarise(Temp_min = min(Temp_min),
Temp_max = max(Temp_max),
Photoperiod_min = min(DayLength),
Photoperiod_max = max(DayLength))
::kable(xx) knitr
ExptShort | Temp_min | Temp_max | Photoperiod_min | Photoperiod_max |
---|---|---|---|---|
Ro16 | -2.251 | 41.123 | 14.51 | 16.62 |
Ro17 | 0.107 | 42.541 | 14.89 | 16.62 |
Su16 | -3.720 | 46.528 | 14.51 | 16.52 |
Su17 | -2.200 | 37.419 | 14.91 | 16.52 |
Su18 | -0.686 | 41.940 | 13.94 | 16.52 |
Us18 | 1.670 | 36.670 | 12.49 | 15.64 |
In16 | 0.900 | 37.500 | 10.58 | 12.14 |
In17 | 5.380 | 35.690 | 10.58 | 11.59 |
Ba16 | 7.900 | 33.100 | 10.57 | 12.02 |
Ba17 | 5.390 | 39.460 | 10.57 | 12.15 |
Ne16 | 2.000 | 38.600 | 10.20 | 12.64 |
Ne17 | NA | NA | 10.19 | 12.82 |
Sp16 | -2.900 | 34.000 | 9.37 | 14.18 |
Sp17 | -3.100 | 31.100 | 9.37 | 14.17 |
Mo16 | -1.700 | 36.900 | 9.77 | 13.59 |
Mo17 | -1.000 | 28.100 | 9.77 | 14.18 |
It16 | -5.832 | 31.250 | 9.11 | 14.73 |
It17 | -2.366 | 32.750 | 9.11 | 14.59 |
Figure 1: Field Trial Info
Figure 1: Growing Environments. (a) Locations of field trials conducted in the summer and winter of 2016, 2017 and 2018, along with (b) mean temperature and photoperiod of each field trial: Rosthern, Canada 2016 and 2017 (Ro16, Ro17), Sutherland, Canada 2016, 2017 and 2018 (Su16, Su17, Su18), Central Ferry, USA 2018 (Us18), Metaponto, Italy 2016 and 2017 (It16, It17), Marchouch, Morocco 2016 and 2017 (Mo16, Mo17), Cordoba, Spain 2016 and 2017 (Sp16, Sp17), Bhopal, India 2016 and 2017 (In16, In17), Jessore, Bangladesh 2016 and 2017 (Ba16, Ba17), Bardiya, Nepal 2016 and 2017 (Ne16, Ne17).
# Plot (a) Map
invisible(png("Additional/Temp/Temp_F01_1.png",
width = 4200, height = 1575, res = 600))
par(mai = c(0,0,0,0), xaxs = "i", yaxs = "i")
mapBubbles(dF = ff, nameX = "Lon", nameY = "Lat",
nameZColour = "MacroEnv", nameZSize = "Year",
symbolSize = 0.5, pch = 20, fill = F, addLegend = F,
colourPalette = c("darkgreen", "darkorange3", "darkblue"),
addColourLegend = F,
xlim = c(-140,110), ylim = c(10,35),
oceanCol = "grey90", landCol = "white", borderCol = "black")
invisible(dev.off())
# Plot (b) mean T and P
<- ggplot(ff, aes(x = T_mean, y = P_mean)) +
mp geom_point(aes(color = MacroEnv), size = 3, alpha = 0.8) +
geom_text_repel(aes(label = ExptShort)) +
scale_x_continuous(breaks = 11:22) + scale_y_continuous(breaks = 11:16) +
scale_color_manual(name = "Macro-Environment",
values = c("darkgreen","darkorange3","darkblue")) +
+
theme_AGL theme(legend.position = "top", legend.margin = unit(c(0,0,0,0), "cm")) +
labs(x = expression(paste("Mean temperature (", degree, "C)", sep = "")),
y = "Mean photoperiod (h)")
ggsave("Additional/Temp/Temp_F01_2.png", mp,
width = 7, height = 3.25, dpi = 600)
# Labels were added to "Additional/Temp/Temp_F1_1.png" in image editing software
# Append (a) and (b)
<- image_read("Additional/Temp/Temp_F01_1_1.png") %>%
im1 image_annotate("(a)", size = 35)
<- image_read("Additional/Temp/Temp_F01_2.png") %>%
im2 image_scale("1200x") %>%
image_annotate("(b)", size = 35)
<- image_append(c(im1, im2), stack = T)
im image_write(im, "Figure_01.png")
Additional Figure 1: LDP Origin Map
# Prep data
<- ldp %>% filter(Origin != "Unknown") %>%
x1 mutate(Origin = recode(Origin, "ICARDA"="Syria", "USDA"="USA")) %>%
group_by(Origin) %>%
summarise(Count = n()) %>%
left_join(select(ct, Origin = Country, Lat, Lon), by = "Origin") %>%
ungroup() %>%
as.data.frame()
is.na(x1)] <- 0
x1[# Plot
invisible(png("Additional/Additional_Figure_01.png",
width = 3600, height = 2055, res = 600))
par(mai = c(0,0,0,0), xaxs = "i",yaxs = "i")
mapBubbles(dF = x1, nameX = "Lon", nameY = "Lat",
nameZSize = "Count", nameZColour = "darkred",
xlim = c(-140,110), ylim = c(5,20),
oceanCol = "grey90", landCol = "white", borderCol = "black")
invisible(dev.off())
Supplemental Figure 1: DTF Scaling
Figure S1: Distribution of days from sowing to flowering for raw data (top) and scaled data (1-5) (bottom) for all 18 field trials: Rosthern, Canada 2016 and 2017 (Ro16, Ro17), Sutherland, Canada 2016, 2017 and 2018 (Su16, Su17, Su18), Central Ferry, USA 2018 (Us18), Metaponto, Italy 2016 and 2017 (It16, It17), Marchouch, Morocco 2016 and 2017 (Mo16, Mo17), Cordoba, Spain 2016 and 2017 (Sp16, Sp17), Bhopal, India 2016 and 2017 (In16, In17), Jessore, Bangladesh 2016 and 2017 (Ba16, Ba17), Bardiya, Nepal 2016 and 2017 (Ne16, Ne17). Genotypes which did not flower were given a scaled value of 5.
# Prep data
<- c("Days from sowing to flower (days)", "Scaled (1-5)")
levs <- dd %>% select(Entry, Expt, ExptShort, DTF, DTF2_scaled) %>%
xx left_join(select(ff, Expt, MacroEnv), by = "Expt") %>%
gather(Trait, Value, DTF, DTF2_scaled) %>%
mutate(Trait = plyr::mapvalues(Trait, c("DTF", "DTF2_scaled"), levs),
Trait = factor(Trait, levels = levs) )
# Plot
<- ggplot(xx, aes(x = ExptShort, y = Value)) +
mp geom_violin(fill = "grey", alpha = 0.3, color = NA) +
geom_quasirandom(aes(color = MacroEnv), size = 0.1, alpha = 0.5) +
scale_color_manual(values = c("darkgreen","darkorange3","darkblue")) +
facet_grid(Trait ~ MacroEnv, scales = "free") +
+
theme_AGL theme(legend.position = "none",
panel.grid.major.x = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(x = NULL, y = NULL)
ggsave("Supplemental_Figure_01.png", mp, width = 8, height = 5, dpi = 600)
Phenology
Figure 2: Data Overview
Figure 2: Variations in temperature, day length and phenological traits across contrasting environment for a lentil (Lens culinaris Medik.) diversity panel. (a) Daily mean temperature (red line) and day length (blue line) from seeding to full maturity of all genotypes. The shaded ribbon represents the daily minimum and maximum temperature. The shaded area between the vertical bars corresponds to the windows of flowering. (b) Distribution of mean days from sowing to: flowering (DTF), swollen pods (DTS) and maturity (DTM), and (c) vegetative (VEG) and reproductive periods (REP) of 324 genotypes across 18 site-years. Rosthern, Canada 2016 and 2017 (Ro16, Ro17), Sutherland, Canada 2016, 2017 and 2018 (Su16, Su17, Su18), Central Ferry, USA 2018 (Us18), Metaponto, Italy 2016 and 2017 (It16, It17), Marchouch, Morocco 2016 and 2017 (Mo16, Mo17), Cordoba, Spain 2016 and 2017 (Sp16, Sp17), Bhopal, India 2016 and 2017 (In16, In17), Jessore, Bangladesh 2016 and 2017 (Ba16, Ba17), Bardiya, Nepal 2016 and 2017 (Ne16, Ne17).
# Create plot function
<- function(x, mybreaks, nr = 2, nc = 3) {
ggEnvPlot <- ff %>%
yy filter(Expt %in% unique(x$Expt)) %>%
select(ExptShort, Location, Year, min=DTF_min, max=DTF_max) %>%
mutate(Trait = "DTF Window")
ggplot(x) +
geom_rect(data = yy, aes(xmin = min, xmax = max, fill = Trait),
ymin = -Inf, ymax = Inf, alpha = 0.1, color = "black") +
geom_line(aes(x = DaysAfterPlanting, y = DayLength_rescaled, color = "Day Length")) +
geom_line(aes(x = DaysAfterPlanting, y = Temp_mean, color = "Temperature") ) +
geom_ribbon(aes(x = DaysAfterPlanting, ymin = Temp_min, ymax = Temp_max),
fill = "darkred", alpha = 0.3) +
facet_wrap(ExptShort ~ ., scales = "free_x", dir = "v", nrow = 2, ncol = 3) +
scale_x_continuous(breaks = mybreaks) +
scale_color_manual(name = NULL, values = c("darkblue", "darkred")) +
scale_fill_manual(name = NULL, values = "darkgreen") +
guides(colour = guide_legend(order = 1, override.aes = list(size = 1.25)),
fill = guide_legend(order = 2)) +
coord_cartesian(ylim=c(0, 40)) +
+
theme_AGL theme(plot.margin = unit(c(0,0,0,0), "cm"),
legend.text = element_text(size = 12)) +
labs(y = NULL, x = NULL)
}# Plot (a) T and P
.1 <- ggEnvPlot(ee %>% filter(MacroEnv == "Temperate"), c(25, 50, 75)) +
mp1labs(y = expression(paste("Degrees Celcius (", degree, "C)"))) +
theme(strip.background = element_rect(alpha("darkgreen", 0.3)),
plot.margin = unit(c(0,0,0,0.155), "cm"))
.2 <- ggEnvPlot(ee %>% filter(MacroEnv == "South Asia"), c(25, 75, 125)) +
mp1labs(x = "Days After Planting") +
theme(strip.background = element_rect(fill = alpha("darkgoldenrod2", 0.3)),
axis.text.y = element_blank(), axis.ticks.y = element_blank())
.3 <- ggEnvPlot(ee %>% filter(MacroEnv == "Mediterranean"), c(50, 100, 150)) +
mp1scale_y_continuous(sec.axis = sec_axis(~ (16.62 - 9.11) * . / (40 - 0) + 9.11,
name = "Hours (h)", breaks = c(10, 12, 14, 16))) +
theme(strip.background = element_rect(fill = alpha("darkblue", 0.3)),
plot.margin = unit(c(0,0.17,0,0), "cm"),
axis.text.y.left = element_blank(), axis.ticks.y.left = element_blank())
<- ggarrange(mp1.1, mp1.2, mp1.3, nrow = 1, ncol = 3, align = "h",
mp1 legend = "top", common.legend = T)
# Prep data
<- dd %>% select(Entry, Year, Expt, ExptShort, Location, DTF, DTS, DTM) %>%
xx left_join(select(ff, Expt, MacroEnv), by = "Expt") %>%
gather(Trait, Value, DTF, DTS, DTM) %>%
mutate(Trait = factor(Trait, levels = c("DTF", "DTS", "DTM")))
# Create plot function
<- function(x) {
ggDistroDTF ggplot(x, aes(x = Trait, y = Value) ) +
geom_violin(color = NA, fill = "grey", alpha = 0.3) +
geom_quasirandom(size = 0.3, aes(color = Trait)) +
facet_wrap(ExptShort ~ ., scales = "free_x", dir = "v", ncol = 3, nrow = 2) +
scale_color_manual(values = c("darkgreen", "darkred", "darkgoldenrod2")) +
scale_y_continuous(limits = c(30,190), breaks = seq(25,175, 25)) +
+ labs(y = NULL, x = NULL) +
theme_AGL theme(plot.margin = unit(c(0.1,0,0.3,0), "cm"))
}# Plot (b) DTF, DTS and DTM
.1 <- ggDistroDTF(xx %>% filter(MacroEnv == "Temperate")) +
mp2labs(y = "Days After Planting") +
theme(strip.background = element_rect(fill = alpha("darkgreen", 0.3)),
panel.grid.major.x = element_blank())
.2 <- ggDistroDTF(xx %>% filter(MacroEnv == "South Asia")) +
mp2theme(strip.background = element_rect(fill = alpha("darkgoldenrod2", 0.3)),
panel.grid.major.x = element_blank(),
axis.text.y = element_blank(), axis.ticks.y = element_blank())
.3 <- ggDistroDTF(xx %>% filter(MacroEnv == "Mediterranean")) +
mp2scale_y_continuous(limits = c(30,190), breaks = seq(25,175, 25),
sec.axis = sec_axis(~ ., name = "Days After Planting",
breaks = seq(25,175, 25))) +
theme(strip.background = element_rect(fill = alpha("darkblue", 0.3)),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_line(),
axis.text.y.left = element_blank(),
axis.ticks.y.left = element_blank() )
# Append
<- ggarrange(mp2.1, mp2.2, mp2.3, nrow = 1, ncol = 3, align = "h", legend = "none")
mp2 # Prep data
<- dd %>% select(Entry, Name, Expt, ExptShort, Location, Year, VEG, REP) %>%
xx left_join(select(ff, Expt, MacroEnv), by = "Expt") %>%
gather(Trait, Value, VEG, REP) %>%
mutate(Trait = factor(Trait, levels = c("VEG", "REP")))
# Create plot function
<- function(x) {
ggDistroREP ggplot(x, aes(x = Trait, y = Value)) +
geom_violin(color = NA, fill = "grey", alpha = 0.3) +
geom_quasirandom(size = 0.3, aes(color = Trait)) +
facet_wrap(ExptShort ~ ., scales = "free_x", dir = "v", ncol = 3, nrow = 2) +
scale_color_manual(values = c("steelblue", "purple")) +
scale_y_continuous(limits = c(25,135), breaks = seq(25,125, 25)) +
+ labs(x = NULL, y = NULL) +
theme_AGL theme(plot.margin = unit(c(0,0,0.3,0), "cm"))
}# Plot (c) REP and VEG
.1 <- ggDistroREP(xx %>% filter(MacroEnv == "Temperate")) + labs(y = "Days") +
mp3theme(strip.background = element_rect(fill = alpha("darkgreen", 0.3)),
panel.grid.major.x = element_blank())
.2 <- ggDistroREP(xx %>% filter(MacroEnv == "South Asia")) +
mp3theme(strip.background = element_rect(fill = alpha("darkgoldenrod2", 0.3)),
axis.text.y = element_blank(), axis.ticks.y = element_blank(),
panel.grid.major.x = element_blank())
.3 <- ggDistroREP(xx %>% filter(MacroEnv == "Mediterranean")) +
mp3scale_y_continuous(limits = c(25,135), breaks = seq(25,125, 25),
sec.axis = sec_axis(~ ., name = "Days", breaks = seq(25,125, 25))) +
theme(strip.background = element_rect(fill = alpha("darkblue", 0.3)),
axis.text.y.left = element_blank(), axis.ticks.y.left = element_blank(),
panel.grid.major.x = element_blank())
# Append
<- ggarrange(mp3.1, mp3.2, mp3.3, nrow = 1, ncol = 3, align = "h", legend = "none")
mp3 # Save
ggsave("Additional/Temp/Temp_F02_1.png", mp1, width = 12, height = 4, dpi = 500, bg = "white")
ggsave("Additional/Temp/Temp_F02_2.png", mp2, width = 12, height = 4, dpi = 500)
ggsave("Additional/Temp/Temp_F02_3.png", mp3, width = 12, height = 4, dpi = 500)
# Append (a), (b) and (c)
<- image_read("Additional/Temp/Temp_F02_1.png") %>% image_annotate("(a)", size = 150)
mp1 <- image_read("Additional/Temp/Temp_F02_2.png") %>% image_annotate("(b)", size = 150)
mp2 <- image_read("Additional/Temp/Temp_F02_3.png") %>% image_annotate("(c)", size = 150)
mp3 <- image_append(c(mp1, mp2, mp3), stack = T)
mp image_write(mp, "Figure_02.png")
Additional Figures: Entry Phenology
# Create plotting function
<- function(x, xE, colnums) {
ggPhenology <- c("darkgreen", "darkorange3", "darkblue")
mycols ggplot(xE, aes(x = Trait, y = Value, group = Entry, color = MacroEnv)) +
geom_line(data = x, color = "grey", alpha = 0.5) +
geom_line() +
geom_point() +
facet_grid(MacroEnv ~ ExptShort) +
scale_color_manual(values = mycols[colnums]) +
+
theme_AGL theme(legend.position = "none", panel.grid.major.x = element_blank()) +
ylim(c(min(x$Value, na.rm = T), max(x$Value, na.rm = T))) +
labs(x = NULL, y = "Days")
}# Prep data
<- dd %>% select(Entry, Name, ExptShort, DTF, DTS, DTM) %>%
xx left_join(select(ff, ExptShort, MacroEnv), by = "ExptShort") %>%
gather(Trait, Value, DTF, DTS, DTM) %>%
mutate(Trait = factor(Trait, levels = c("DTF","DTS","DTM")))
<- xx %>% filter(MacroEnv == "Temperate")
x1 <- xx %>% filter(MacroEnv == "South Asia")
x2 <- xx %>% filter(MacroEnv == "Mediterranean")
x3 # Create PDF
pdf("Additional/pdf_Phenology.pdf", width = 8, height = 6)
for(i in 1:324) {
<- xx %>% filter(Entry == i, !is.na(Value), MacroEnv == "Temperate")
xE1 <- xx %>% filter(Entry == i, !is.na(Value), MacroEnv == "South Asia")
xE2 <- xx %>% filter(Entry == i, !is.na(Value), MacroEnv == "Mediterranean")
xE3 <- ggPhenology(x1, xE1, 1)
mp1 <- ggPhenology(x2, xE2, 2)
mp2 <- ggPhenology(x3, xE3, 3)
mp3 <- paste("Entry", str_pad(i, 3, "left", "0"), "|", unique(xE1$Name))
figlab <- ggarrange(mp1, mp2, mp3, nrow = 3, ncol = 1) %>%
mp annotate_figure(top = figlab)
print(mp)
ggsave(paste0("Additional/Entry_Phenology/Phenology_Entry_",
str_pad(i, 3, "left", "0"), ".png"),
width = 8, height = 6, dpi = 600)
mp,
}dev.off()
Additional Figure 2: DTF DTS DTM
# Prep data
<- dd %>% select(Entry, Expt, ExptShort, DTF, DTS, DTM) %>%
xx left_join(select(ff, Expt, MacroEnv), by = "Expt") %>%
gather(Trait, Value, DTF, DTS, DTM) %>%
mutate(Trait = factor(Trait, levels = c("DTF", "DTS", "DTM")) )
# Plot
<- ggplot(xx, aes(x = ExptShort, y = Value)) +
mp geom_violin(fill = "grey", alpha = 0.25, color = NA) +
geom_quasirandom(size = 0.1, alpha = 0.5, aes(color = MacroEnv)) +
facet_grid(Trait ~ MacroEnv, scales = "free") +
scale_color_manual(values = c("darkgreen", "darkorange3", "darkblue")) +
+
theme_AGL theme(legend.position = "none",
panel.grid.major.x = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(x = NULL, y = "Days After Planting")
ggsave("Additional/Additional_Figure_02.png", mp, width = 8, height = 6, dpi = 600)
Additional Figure 3: ggridges
# Prep data
<- dd %>% select(Expt, DTF, DTS, DTM) %>%
xx gather(Trait, Value, DTF, DTS, DTM) %>%
mutate(Trait = factor(Trait, levels = c("DTF", "DTS", "DTM")))
# Plot
<- ggplot(xx, aes(x = Value, y = Expt, fill = Trait)) +
mp ::geom_density_ridges(alpha = 0.7) +
ggridgesscale_fill_manual(name = NULL, values = c("darkgreen", "darkred", "darkgoldenrod2")) +
+
theme_AGL theme(legend.position = "top", legend.margin = unit(c(0,0,0,0), "cm")) +
labs(y = NULL, x = "Days After Sowing")
ggsave("Additional/Additional_Figure_03.png", mp, width = 6, height = 4, dpi = 600)
Additional Figure 4: MacroEnv Phenology
# Prep data
<- ee %>% filter(ExptShort %in% c("Su17", "Ba17", "It17"))
xx <- ff %>% filter(Expt %in% unique(xx$Expt)) %>%
yy mutate(DTF_min = Start + DTF_min, DTF_max = Start + DTF_max,
DTM_min = Start + DTM_min, DTM_max = Start + DTM_max)
<- select(yy, Expt, Location, Year, MacroEnv, min = DTF_min, max = DTF_max) %>%
y1 mutate(Trait = "DTF")
<- select(yy, Expt, Location, Year, MacroEnv, min = DTM_min, max = DTM_max) %>%
y2 mutate(Trait = "DTM")
<- bind_rows(y1, y2)
yy # Plot
<- ggplot(xx) +
mp geom_rect(data = yy, aes(xmin = min, xmax = max, fill = Trait),
ymin = 0, ymax = 40, alpha = 0.4) +
geom_line(aes(x = Date, y = DayLength_rescaled, color = "Blue")) +
geom_line(aes(x = Date, y = Temp_mean, color = "darkred") ) +
geom_ribbon(aes(x = Date, ymin = Temp_min, ymax = Temp_max),
fill = alpha("darkred", 0.25), color = alpha("darkred", 0.25)) +
facet_grid(Location + MacroEnv ~ ., scales = "free_x", space = "free_x") +
scale_color_manual(name = NULL, values = c("Blue", "darkred"),
labels = c("Day length", "Temperature") ) +
scale_fill_manual(name = NULL, values = c("darkgreen", "darkgoldenrod2")) +
coord_cartesian(ylim = c(0,40)) +
+
theme_AGL theme(legend.position = "bottom",
legend.text = element_text(size = 12),
axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_x_date(breaks = "1 month", labels = date_format("%b")) +
scale_y_continuous(sec.axis = sec_axis(~ (16.62 - 9.11) * . / (40 - 0) + 9.11,
breaks = c(10, 12, 14, 16), name = "Hours")) +
guides(colour = guide_legend(order = 1, override.aes = list(size = 1.25)),
fill = guide_legend(order = 2)) +
labs(title = "2017 - 2018", x = NULL,
y = expression(paste(degree, "Celcius")))
ggsave("Additional/Additional_Figure_04.png", mp, width = 8, height = 6, dpi = 600)
Additional Figures: Phenology + EnvData
# Plotting function
<- function(exptshort = "Ro17") {
ggPhenologyEnvD # Prep data
<- dd %>% filter(ExptShort == exptshort) %>%
x1 select(Entry, Name, ExptShort, Expt, DTF, DTS, DTM) %>%
gather(Trait, Value, DTF, DTS, DTM) %>%
mutate(Trait = factor(Trait, levels = c("DTF", "DTS", "DTM")))
<- ee %>% filter(ExptShort == exptshort)
x2 <- max(x2$DaysAfterPlanting)
myMax <- c("steelblue", "darkgreen", "darkgoldenrod2")
myColors # Plot
<- ggplot(x1, aes(x = Value, fill = Trait)) +
mp1 geom_density(alpha = 0.7) +
facet_grid(. ~ Expt) +
scale_fill_manual(name = NULL, values = myColors) +
coord_cartesian(xlim = c(0,myMax)) +
+
theme_AGL theme(legend.position = "top",
axis.text.x = element_blank(),
axis.ticks.x = element_blank()) +
labs(x = NULL, y = NULL)
<- ggplot(x2) +
mp2 geom_line(aes(x = DaysAfterPlanting, y = DayLength_rescaled, color = "Blue")) +
geom_line(aes(x = DaysAfterPlanting, y = Temp_mean, color = "darkred") ) +
geom_ribbon(aes(x = DaysAfterPlanting, ymin = Temp_min, ymax = Temp_max),
fill = alpha("darkred", 0.25), color = alpha("darkred", 0.25)) +
scale_color_manual(name = NULL, values = c("Blue", "darkred"),
labels = c("Day length", "Temperature") ) +
coord_cartesian(ylim = c(0,40)) +
+
theme_AGL theme(legend.position = "bottom",
legend.text = element_text(size = 12)) +
scale_x_continuous(limits = c(0, myMax)) +
scale_y_continuous(sec.axis = sec_axis(~ (16.62 - 9.11) * . / (40 - 0) + 9.11,
breaks = c(10, 12, 14, 16), name = "Hours")) +
guides(colour = guide_legend(order = 1, override.aes = list(size = 1.25)),
fill = guide_legend(order = 2)) +
labs(y = expression(paste(degree, "Celcius"), x = NULL))
ggarrange(mp1, mp2, nrow = 2, align = "v", heights = c(1, 1.2))
}for(i in names_ExptShort) {
<- ggPhenologyEnvD(i)
mp ggsave(paste0("Additional/Expt/1_", i, ".png"), width = 6, height = 6)
}<- image_read(paste0("Additional/Expt/1_", names_ExptShort[1:6], ".png"))
im <- image_append(c(image_append(im[1:3]), image_append(im[4:6])), stack = T)
im image_write(im, "Additional/Expt/2_Temperate.png")
<- image_read(paste0("Additional/Expt/1_", names_ExptShort[7:12], ".png"))
im <- image_append(c(image_append(im[1:3]), image_append(im[4:6])), stack = T)
im image_write(im, "Additional/Expt/2_SouthAsia.png")
<- image_read(paste0("Additional/Expt/1_", names_ExptShort[13:18], ".png"))
im <- image_append(c(image_append(im[1:3]), image_append(im[4:6])), stack = T)
im image_write(im, "Additional/Expt/2_Mediterranean.png")
Supplemental Figure 2: Missing Data
Figure S2: Percentage of lentil genotypes reaching key phenological time points in South Asian locations. Days from sowing to: flowering (DTF), swollen pods (DTS) and maturity (DTM).
# Prep data
<- dd %>%
xx filter(Location %in% c("Bhopal, India", "Jessore, Bangladesh", "Bardiya, Nepal")) %>%
mutate(DTF = ifelse(is.na(DTF), 0, 1),
DTS = ifelse(is.na(DTS), 0, 1),
DTM = ifelse(is.na(DTM), 0, 1) ) %>%
group_by(Expt, Location, Year) %>%
summarise_at(vars(DTF, DTS, DTM), funs(sum), na.rm = T) %>%
ungroup() %>%
gather(Trait, Flowered, DTF, DTS, DTM) %>%
mutate(Total = ifelse(Expt == "Bardiya, Nepal 2016", 323, 324),
# One accession was not planted in Bardiya, Nepal 2016
DidNotFlower = Total - Flowered,
Percent = round(100 * Flowered / Total),
Label = paste0(Percent, "%"),
Trait = factor(Trait, levels = c("DTF", "DTS", "DTM")))
# Plot
<- ggplot(xx, aes(x = Trait, y = Percent, fill = Trait)) +
mp geom_bar(stat = "identity", color = "black", alpha = 0.7) +
geom_text(aes(label = Label), nudge_y = -3, size = 3.5) +
facet_grid(. ~ Location + Year) +
scale_fill_manual(values = c("darkgreen", "darkred", "darkgoldenrod2")) +
scale_y_continuous(limits = c(0,100), expand = c(0,0)) +
+
theme_AGL theme(legend.position = "none",
panel.grid.major.x = element_blank() ) +
labs(x = NULL, y = "Percent of accessions reaching key phenological time points")
ggsave("Supplemental_Figure_02.png", width = 10, height = 5, dpi = 600)
Supplemental Figure 3: Correlation Plots
Figure S3: Correlations along with the corresponding correlation coefficients (R2) between days from sowing to: flowering (DTF), swollen pod (DTS) and maturity (DTM), in temperate (top), South Asian (middle) and Mediterranean (bottom) locations.
# Prep data
<- dd %>%
xx left_join(select(ff, Expt, MacroEnv), by = "Expt") %>%
select(Entry, Expt, MacroEnv, DTF, DTS, DTM)
# Create plotting function
<- function(x, legend.title, colNums) {
ggCorPlot # Plot (a)
<- round(cor(x$DTF, x$DTS, use ="complete", method = "pearson")^2, 2)
r2 <- ggplot(x) +
tp1 geom_point(aes(x = DTF, y = DTS, color = Expt, shape = Expt), alpha = 0.8) +
geom_label(x = -Inf, y = Inf, hjust = 0, vjust = 1, parse = T,
label = paste("italic(R)^2 == ", r2) ) +
scale_color_manual(name = legend.title, values = colors_Expt[colNums]) +
scale_shape_manual(name = legend.title, values = c(15,16,17,15,16,17)) +
theme_AGL# Plot (b)
<- round(cor(x$DTF, x$DTM, use ="complete.obs", method = "pearson")^2, 2)
r2 <- ggplot(x) +
tp2 geom_point(aes(x = DTF, y = DTM, color = Expt, shape = Expt), alpha = 0.8) +
geom_label(x = -Inf, y = Inf, hjust = 0, vjust = 1, parse = T,
label = paste("italic(R)^2 == ", r2) ) +
scale_color_manual(name = legend.title, values = colors_Expt[colNums]) +
scale_shape_manual(name = legend.title, values = c(15,16,17,15,16,17)) +
theme_AGL# Plot (c)
<- round(cor(x$DTS, x$DTM, use = "complete", method = "pearson")^2, 2)
r2 <- ggplot(x) +
tp3 geom_point(aes(x = DTS, y = DTM, color = Expt, shape = Expt), alpha = 0.8) +
geom_label(x = -Inf, y = Inf, hjust = 0, vjust = 1, parse = T,
label = paste("italic(R)^2 == ", r2) ) +
scale_color_manual(name = legend.title, values = colors_Expt[colNums]) +
scale_shape_manual(name = legend.title, values = c(15,16,17,15,16,17)) +
theme_AGL# Append (a), (b) and (c)
<- ggarrange(tp1, tp2, tp3, nrow = 1, ncol = 3, common.legend = T, legend = "right")
mp
mp
}# Plot
<- ggCorPlot(xx %>% filter(MacroEnv == "Temperate"), "Temperate", 1:6 )
mp1 <- ggCorPlot(xx %>% filter(MacroEnv == "South Asia"), "South Asia", 7:12)
mp2 <- ggCorPlot(xx %>% filter(MacroEnv == "Mediterranean"), "Mediterranean", 13:18)
mp3 <- ggarrange(mp1, mp2, mp3, nrow = 3, ncol = 1, common.legend = T, legend = "right")
mp ggsave("Supplemental_Figure_03.png", mp, width = 10, height = 8, dpi = 600)
Additional Figures: Correlations
# Prep data
<- dd %>%
xx left_join(select(ff, Expt, MacroEnv), by = "Expt") %>%
mutate(DTE = ifelse(Location == "Cordoba, Spain", NA, DTE))
<- xx %>% filter(MacroEnv == "Temperate")
x1 <- xx %>% filter(MacroEnv == "South Asia")
x2 <- xx %>% filter(MacroEnv == "Mediterranean")
x3 # Create plotting functions
<- function(data, mapping, cols = colors_Expt, ...) {
my_lower ggplot(data = data, mapping = mapping) +
geom_point(alpha = 0.5, size = 0.3, aes(color = Expt)) +
scale_color_manual(values = cols) +
theme_bw() +
theme(axis.text = element_text(size = 7.5))
}<- function(data, mapping, cols = colors_Expt, ...) {
my_middle ggplot(data = data, mapping = mapping) +
geom_density(alpha = 0.5) +
scale_color_manual(name = NULL, values = cols) +
scale_fill_manual(name = NULL, values = cols) +
guides(color = F, fill = guide_legend(nrow = 3, byrow = T)) +
theme_bw() +
theme(axis.text = element_text(size = 7.5))
}# See: https://github.com/ggobi/ggally/issues/139
<- function(data, mapping, color = I("black"), sizeRange = c(1,5), ...) {
my_upper # Prep data
<- eval_data_col(data, mapping$x)
x <- eval_data_col(data, mapping$y)
y #
<- cor(x, y, method = "pearson", use = "complete.obs")^2
r2 <- format(r2, digits = 2)[1]
rt <- max(sizeRange)
cex <- as.character(rt)
tt # plot the cor value
<- ggally_text(label = tt, mapping = aes(), color = color,
p xP = 0.5, yP = 0.5, size = 6, ... ) + theme_bw()
# Create color palette
<- RColorBrewer::brewer.pal(n = 10, name = "RdBu")[2:9]
corColors if (r2 <= -0.9) { corCol <- alpha(corColors[1], 0.5)
else if (r2 >= -0.9 & r2 <= -0.6) { corCol <- alpha(corColors[2], 0.5)
} else if (r2 >= -0.6 & r2 <= -0.3) { corCol <- alpha(corColors[3], 0.5)
} else if (r2 >= -0.3 & r2 <= 0) { corCol <- alpha(corColors[4], 0.5)
} else if (r2 >= 0 & r2 <= 0.3) { corCol <- alpha(corColors[5], 0.5)
} else if (r2 >= 0.3 & r2 <= 0.6) { corCol <- alpha(corColors[6], 0.5)
} else if (r2 >= 0.6 & r2 <= 0.9) { corCol <- alpha(corColors[7], 0.5)
} else { corCol <- alpha(corColors[8], 0.5) }
} # Plot
<- p +
p theme(panel.background = element_rect(fill = corCol),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p
}# Plot Correlations for each Expt
for(i in 1:length(names_Expt)) {
<- ggpairs(xx %>% filter(Expt == names_Expt[i]),
mp columns = c("DTE", "DTF", "DTS", "DTM", "REP"),
upper = list(continuous = my_upper),
diag = list(continuous = my_middle),
lower = list(continuous = wrap(my_lower, cols = "black")),
title = names_Expt[i]) +
theme(strip.background = element_rect(fill = "White"))
ggsave(paste0("Additional/Corr/Corr_", str_pad(i, 2, "left", "0"),
"_", names_ExptShort[i], ".png"),
width = 6, height = 6, dpi = 600)
mp,
}# Plot (a) Temperate
<- ggpairs(x1, columns = c("DTE", "DTF", "DTS", "DTM", "REP"),
mp1 aes(color = Expt, fill = Expt),
upper=list(continuous = my_upper),
diag =list(continuous = wrap(my_middle, cols = colors_Expt[1:6])),
lower=list(continuous = wrap(my_lower, cols = colors_Expt[1:6])),
title = "(a) Temperate",
legend = c(2,2)) +
theme(strip.background = element_rect(fill = "White"),
legend.position = "bottom")
ggsave("Additional/Corr/Corr_Temperate.png", mp1, width = 6, height = 6, dpi = 600)
# Plot (b) South Asia
<- ggpairs(x2, columns = c("DTE", "DTF", "DTS", "DTM", "REP"),
mp2 aes(color = Expt, fill = Expt),
upper = list(continuous = my_upper),
diag = list(continuous = wrap(my_middle, cols = colors_Expt[7:12])),
lower = list(continuous = wrap(my_lower, cols = colors_Expt[7:12])),
title = "(b) South Asia",
legend = c(2,2)) +
theme(strip.background = element_rect(fill = "White"),
legend.position = "bottom")
ggsave("Additional/Corr/Corr_SouthAsia.png", mp2, width = 6, height = 6, dpi = 600)
# Plot (c) Mediterranean
<- ggpairs(x3, columns = c("DTE", "DTF", "DTS", "DTM", "REP"),
mp3 aes(color = Expt, fill = Expt),
upper = list(continuous = my_upper),
diag = list(continuous = wrap(my_middle, cols = colors_Expt[13:18])),
lower = list(continuous = wrap(my_lower, cols = colors_Expt[13:18])),
title = "(c) Mediterranean",
legend = c(2,2)) +
theme(strip.background = element_rect(fill = "White"),
legend.position = "bottom")
ggsave("Additional/Corr/Corr_Mediterranean.png", mp3, width = 6, height = 6, dpi = 600)
# Plot All
<- ggpairs(xx, columns = c("DTE", "DTF", "DTS", "DTM", "REP"),
mp4 aes(color = ExptShort, fill = ExptShort),
upper = list(continuous = my_upper),
diag = list(continuous = my_middle),
lower = list(continuous = my_lower),
title = "ALL",
legend = c(2,2)) +
theme(strip.background = element_rect(fill = "White"),
legend.position = "bottom")
ggsave("Additional/Corr/Corr_All.png", mp4, width = 6, height = 6, dpi = 600)
PCA
Figure 3: PCA Clusters
Figure 3: Clustering of a lentil (Lens culinaris Medik.) diversity panel based days from sowing to flower (DTF). (a) Principal Component Analysis on DTF, scaled from 1-5, and hierarchical k-means clustering into eight groups. (b) Mean scaled DTF (1-5) for each cluster group across all field trials: Rosthern, Canada 2016 and 2017 (Ro16, Ro17), Sutherland, Canada 2016, 2017 and 2018 (Su16, Su17, Su18), Central Ferry, USA 2018 (Us18), Metaponto, Italy 2016 and 2017 (It16, It17), Marchouch, Morocco 2016 and 2017 (Mo16, Mo17), Cordoba, Spain 2016 and 2017 (Sp16, Sp17), Bhopal, India 2016 and 2017 (In16, In17), Jessore, Bangladesh 2016 and 2017 (Ba16, Ba17), Bardiya, Nepal 2016 and 2017 (Ne16, Ne17). Shaded areas represent one standard deviation from the mean. Dashed, vertical bars separate temperate, South Asian and Mediterranean macro-environments. (c) Composition of cluster groups in genotypes by country of origin. Pie size is relative to the number of genotypes originating from that country.
# Prep data
<- dd %>%
xx select(Entry, Expt, DTF2_scaled) %>%
spread(Expt, DTF2_scaled)
<- xx %>% column_to_rownames("Entry") %>% as.matrix()
xx # PCA
<- PCA(xx, ncp = 10, graph = F)
mypca # Heirarcical clustering
<- HCPC(mypca, nb.clust = 8, graph = F)
mypcaH <- round(mypca[[1]][,2], 1)
perc <- mypcaH[[1]] %>% rownames_to_column("Entry")
x1 <- mypca[[3]]$coord %>% as.data.frame() %>% rownames_to_column("Entry")
x2 <- left_join(x1, x2, by = "Entry") %>%
pca mutate(Entry = as.numeric(Entry)) %>%
left_join(select(ldp, Entry, Name, Origin), by = "Entry") %>%
left_join(select(ct, Origin=Country, Region), by = "Origin") %>%
select(Entry, Name, Origin, Region, Cluster=clust,
PC1=Dim.1, PC2=Dim.2, PC3=Dim.3, PC4=Dim.4, PC5=Dim.5,
PC6=Dim.6, PC7=Dim.7, PC8=Dim.8, PC9=Dim.9, PC10=Dim.10)
# Save data
write.csv(pca, "data/data_pca_results.csv", row.names = F)
# Prep data
<- dd %>%
x2 left_join(select(pca, Entry, Cluster), by = "Entry") %>%
group_by(Expt, ExptShort, Cluster) %>%
summarise(mean = mean(DTF2_scaled, na.rm = T),
sd = sd(DTF2_scaled, na.rm = T)) %>%
ungroup() %>%
mutate(ClusterNum = plyr::mapvalues(Cluster, as.character(1:8), summary(pca$Cluster)))
<- pca %>%
x3 count(Cluster) %>%
mutate(Cluster = factor(Cluster, levels = rev(levels(Cluster))), y = n/2)
for(i in 2:nrow(x3)) { x3$y[i] <- sum(x3$n[1:(i-1)]) + (x3$n[i]/2) }
# Plot (a) PCA 1v2
<- function(df) df[chull(df[,"PC1"], df[,"PC2"]), ]
find_hull <- plyr::ddply(pca, "Cluster", find_hull) %>% mutate(Cluster = factor(Cluster))
polys .1 <- ggplot(pca) +
mp1geom_polygon(data = polys, alpha = 0.15, aes(x = PC1, y = PC2, fill = Cluster)) +
geom_point(aes(x = PC1, y = PC2, colour = Cluster)) +
scale_fill_manual(values = colors) +
scale_color_manual(values = colors) +
+
theme_AGL theme(legend.position = "none", panel.grid = element_blank()) +
labs(x = paste0("PC1 (", perc[1], "%)"),
y = paste0("PC2 (", perc[2], "%)"))
# Plot (a) PCA 1v3
<- function(df) df[chull(df[,"PC1"], df[,"PC3"]), ]
find_hull <- plyr::ddply(pca, "Cluster", find_hull) %>% mutate(Cluster = factor(Cluster))
polys .2 <- ggplot(pca) +
mp1geom_polygon(data = polys, alpha = 0.15, aes(x = PC1, y = PC3, fill = Cluster)) +
geom_point(aes(x = PC1, y = PC3, colour = Cluster)) +
scale_fill_manual(values = colors) +
scale_color_manual(values = colors) +
+
theme_AGL theme(legend.position = "none", panel.grid = element_blank()) +
labs(x = paste0("PC1 (", perc[1], "%)"),
y = paste0("PC3 (", perc[3], "%)"))
# Plot (a) PCA 2v3
<- function(df) df[chull(df[,"PC2"], df[,"PC3"]), ]
find_hull <- plyr::ddply(pca, "Cluster", find_hull) %>% mutate(Cluster = factor(Cluster))
polys .3 <- ggplot(pca) +
mp1geom_polygon(data = polys, alpha = 0.15, aes(x = PC2, y = PC3, fill = Cluster)) +
geom_point(aes(x = PC2, y = PC3, colour = Cluster)) +
scale_fill_manual(values = colors) +
scale_color_manual(values = colors) +
+
theme_AGL theme(legend.position = "none", panel.grid = element_blank()) +
labs(x = paste0("PC2 (", perc[2], "%)"),
y = paste0("PC3 (", perc[3], "%)"))
# Append
<- ggarrange(mp1.1, mp1.2, mp1.3, nrow = 1, ncol = 3, hjust = 0)
mp1 # Plot (b) DTF
<- ggplot(x2, aes(x = ExptShort, y = mean, group = Cluster)) +
mp2 geom_point(aes(color = Cluster)) +
geom_vline(xintercept = 6.5, lty = 2) +
geom_vline(xintercept = 12.5, lty = 2) +
geom_ribbon(aes(ymin = mean - sd, ymax = mean + sd, fill = Cluster),
alpha = 0.2, color = NA) +
geom_line(aes(color = Cluster), size = 1) +
scale_color_manual(values = colors) +
scale_fill_manual(values = colors) +
coord_cartesian(ylim = c(0.95,5.05), expand = F) +
+
theme_AGL theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none", panel.grid = element_blank()) +
labs(y = "DTF (scaled 1-5)", x = NULL)
#
ggsave("Additional/Temp/Temp_F03_1.png", mp1, width = 8, height = 1*6/2.5, dpi = 600)
ggsave("Additional/Temp/Temp_F03_2.png", mp2, width = 8, height = 1.5*6/2.5, dpi = 600)
# Plot (c)
<- ldp %>% left_join(select(pca, Entry, Cluster), by = "Entry")
xx <- xx %>%
xx filter(!Origin %in% c("ICARDA","USDA","Unknown")) %>%
group_by(Origin, Cluster) %>%
summarise(Count = n()) %>%
spread(Cluster, Count) %>%
left_join(select(ct, Origin=Country, Lat, Lon), by = "Origin") %>%
ungroup() %>%
as.data.frame()
is.na(xx)] <- 0
xx[invisible(png("Additional/Temp/Temp_F03_3.png", width = 4800, height = 2200, res = 600))
par(mai = c(0,0,0,0), xaxs = "i", yaxs = "i")
mapPies(dF = xx, nameX = "Lon", nameY = "Lat", zColours = colors,
nameZs = c("1","2","3","4","5","6","7","8"), lwd = 1,
xlim = c(-140,110), ylim = c(0,20), addCatLegend = F,
oceanCol = "grey90", landCol = "white", borderCol = "black")
legend(-139.5, 15.5, title = "Cluster", legend = 1:8, col = colors,
pch = 16, cex = 1, pt.cex = 2, box.lwd = 2)
invisible(dev.off())
# Append (a), (b) and (c)
<- image_read("Additional/Temp/Temp_F03_1.png") %>%
im1 image_annotate("(a)", size = 120, location = "+0+0")
<- image_read("Additional/Temp/Temp_F03_2.png") %>%
im2 image_annotate("(b)", size = 120, location = "+0+0")
<- image_read("Additional/Temp/Temp_F03_3.png") %>%
im3 image_annotate("(c)", size = 120, location = "+0+0")
<- image_append(c(im1, im2, im3), stack = T)
im image_write(im, "Figure_03.png")
Additional Figure 5: PCA Plot
# Prep data
<- read.csv("data/data_pca_results.csv") %>%
xx mutate(Cluster = factor(Cluster),
PC2 = PC2*2.5, PC3 = PC3*2.5,
myColors = plyr::mapvalues(Cluster, 1:8, colors))
# Plot png
<- ggplot(xx, aes(x = PC1, y = PC2, color = Cluster)) +
mp scale_color_manual(values = colors) +
geom_point() +
+
theme_AGL theme(panel.grid.major.x = element_blank(),
panel.grid.major.y = element_blank(),
axis.text = element_blank())
ggsave("Additional/Additional_Figure_05.png", mp, width = 6, height = 4, dpi = 600)
# Plot html
<- plot_ly(xx, x = ~PC1, y = ~PC2, z = ~PC3,
mp color = ~Cluster, colors = colors,
text = ~paste(Entry, "|", Name, "\nOrigin:", Origin)) %>%
add_markers()
saveWidget(mp, file="Additional/Additional_Figure_05.html")
# Prep data
<- read.csv("data/data_pca_results.csv") %>%
xx mutate(Cluster = factor(Cluster),
PC2 = PC2*2.5, PC3 = PC3*2.5,
myColors = plyr::mapvalues(Cluster, 1:8, colors))
# Create plotting function
<- function(x = xx, i = 1) {
ggPCA # Plot
par(mar=c(1.5, 2.5, 1.5, 0.5))
scatter3D(x = x$PC1, y = x$PC3, z = x$PC2, pch = 16, cex = 0.75,
col = alpha(colors,0.8), colvar = as.numeric(x$Cluster), colkey = F,
col.grid = "gray90", bty = "u",
theta = i, phi = 15,
#ticktype = "detailed", cex.lab = 1, cex.axis = 0.5,
xlab = "PC1", ylab = "PC3", zlab = "PC2", main = "LDP - DTF")
}# Plot each angle
for (i in seq(0, 360, by = 2)) {
png(paste0("Additional/PCA/PCA_", str_pad(i, 3, pad = "0"), ".png"),
width = 1000, height = 1000, res = 300)
ggPCA(xx, i)
dev.off()
}# Create animation
<- list.files("Additional/PCA")
lf <- image_read(paste0("Additional/PCA/", lf))
mp <- image_animate(mp, fps = 10)
animation image_write(animation, "Additional/Animation_PCA.gif")
Additional Figure 6: PCA by Region
# Prep data
<- c("Africa", "Asia", "Europe", "Americas", "ICARDA", "USDA")
levs <- read.csv("data/data_pca_results.csv") %>%
xx filter(Origin != "Unknown") %>%
mutate(Cluster = factor(Cluster),
PC2 = PC2*2.5, PC3 = PC3*2.5,
Origin = as.character(Origin),
Region = as.character(Region),
Region = ifelse(Origin %in% levs[5:6], Origin, Region),
Region = factor(Region, levels = levs))
# Plot png
<- ggplot(xx, aes(x = PC1, y = PC2, color = Region, shape = Region)) +
mp geom_point(alpha = 0.7, size = 2) +
scale_color_manual(values = colors[c(1,3,8,5,7,7)]) +
scale_shape_manual(values = c(16,18,15,17,13,10)) +
+
theme_AGL theme(panel.grid.major.x = element_blank(),
panel.grid.major.y = element_blank(),
axis.text = element_blank())
ggsave("Additional/Additional_Figure_06.png", mp, width = 6, height = 4, dpi = 600)
# Plot html
<- plot_ly(xx, x = ~PC1, y = ~PC2, z = ~PC3, color = ~Region,
mp colors = colors[c(1,3,8,5,7,7)],
text = ~paste(Entry, "|", Name, "\nOrigin:", Origin,
"\nCluster:", Cluster)) %>%
add_markers()
saveWidget(mp, file="Additional/Additional_Figure_06.html")
Additional Figure 7: DTF By Cluster
# Prep data
<- read.csv("data/data_pca_results.csv") %>%
x1 mutate(Cluster = factor(Cluster))
<- c("India", "Bangladesh", "Ethiopia", "Pakistan",
x2 "Jordan", "Iran", "Turkey", "Syria", "Chile",
"Spain", "Czech Republic", "Canada" )
<- dd %>%
xx left_join(ldp, by = "Entry") %>%
filter(ExptShort %in% c("Ro17", "Ba17", "Ne17", "It17"), Origin != "Unknown") %>%
left_join(select(x1, Entry, Cluster), by = "Entry") %>%
mutate(Origin = factor(Origin, levels = unique(Origin)[rev(order(unique(Origin)))])) %>%
filter(Origin %in% x2) %>%
mutate(Origin = factor(Origin, levels = x2))
# Plot
<- ggplot(xx, aes(y = DTF2, x = Origin, group = Expt)) +
mp geom_quasirandom(aes(color = Cluster)) +
facet_grid(Location + Year ~ ., scales = "free_y") +
scale_color_manual(values = colors) +
+
theme_AGL theme(legend.position = "top",
panel.grid.major.x = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)) +
guides(colour = guide_legend(nrow = 1, override.aes = list(size = 3))) +
labs(y = "Days from sowing to flower")
ggsave("Additional/Additional_Figure_07.png", mp, width = 5, height = 7, dpi = 600)
Additional Figure 8: Cluster Origins
# Prep data
<- read.csv("data/data_pca_results.csv") %>%
xx mutate(Cluster = factor(Cluster))
<- ldp %>%
xx left_join(select(xx, Entry, Cluster), by = "Entry")
<- xx %>% filter(Origin != "ICARDA") %>%
x1 group_by(Origin, Cluster) %>%
summarise(Count = n()) %>%
spread(Cluster, Count) %>%
left_join(select(ct, Origin=Country, Lat, Lon), by = "Origin") %>%
ungroup() %>%
as.data.frame()
is.na(x1)] <- 0
x1[# Plot
<- ggplot(xx, aes(x = Origin, fill = Cluster)) +
mp geom_bar(stat = "count") +
facet_grid(. ~ Cluster, scales = "free", space = "free") +
scale_fill_manual(values = colors) +
+
theme_AGL theme(legend.position = "none",
panel.grid.major.x = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(y = "Count", x = NULL)
ggsave("Additional/Additional_Figure_08.png", width = 16, height = 4, dpi = 600)
Additional Figure 9: LDP Origins By Cluster
# Prep data
<- read.csv("data/data_pca_results.csv") %>%
xx mutate(Cluster = factor(Cluster))
<- ldp %>%
xx select(Entry, Name, Lat, Lon, Origin) %>%
left_join(select(xx, Entry, Cluster), by = "Entry") %>%
left_join(select(ct, Origin=Country, cLat=Lat, cLon=Lon), by = "Origin") %>%
mutate(Lat = ifelse(is.na(Lat), cLat, Lat),
Lon = ifelse(is.na(Lon), cLon, Lon),
Lat = ifelse(duplicated(Lat), jitter(Lat, 1, 1), Lat),
Lon = ifelse(duplicated(Lon), jitter(Lon, 1, 1), Lon),
Size = 1)
# Plot png
invisible(png("Additional/Additional_Figure_09.png",
width = 7200, height = 4200, res = 600))
par(mai = c(0,0,0,0), xaxs = "i", yaxs = "i")
mapBubbles(dF = xx, nameX = "Lon", nameY = "Lat", nameZColour = "Cluster",
pch = 20, nameZSize = "Size", symbolSize = 0.5,
xlim = c(-140,110), ylim = c(0,20),
colourPalette = colors[1:8], addColourLegend = F, addLegend = F,
oceanCol = "grey90", landCol = "white", borderCol = "black")
legend(-139.5, 15.5, title = "Cluster", legend = 1:8, col = colors,
pch = 16, cex = 1.75, pt.cex = 3, box.lwd = 2)
invisible(dev.off())
# Plot html
<- colorFactor(colors, domain = 1:8)
pal <- leaflet() %>%
mp addProviderTiles(providers$CartoDB.Positron) %>%
addCircles(lng = xx$Lon, lat = xx$Lat, weight = 10,
color = pal(xx$Cluster), opacity = 1, fillOpacity = 1,
popup = paste(xx$Entry,"|",xx$Name,"\nCluster:",xx$Cluster)) %>%
addLegend("bottomleft", pal = pal, values = xx$Cluster,
title = "Cluster", opacity = 1)
saveWidget(mp, file="Additional/Additional_Figure_09.html")
Modeling DTF
Additional Figures: Entry Regressions
<- alpha(c("darkgreen", "darkorange3", "darkblue"), 0.5)
myfills <- min(rr$RDTF, na.rm = T)
mymin <- max(rr$RDTF, na.rm = T)
mymax <- list()
mp for(i in 1:324) {
<- rr %>%
xx filter(Entry == i) %>%
left_join(select(ff, Expt, MacroEnv, T_mean, P_mean), by = "Expt") %>%
mutate(myfill = MacroEnv)
<- xx %>% filter(MacroEnv != "South Asia")
x1 <- xx %>% filter(MacroEnv != "Temperate")
x2 <- xx %>% filter(MacroEnv != "Mediterranean")
x3 <- paste("Entry", str_pad(i, 3, "left", "0"), "|", unique(xx$Name))
figlab # Plot (a) 1/f = a + bT
<- ggplot(xx, aes(x = T_mean, y = RDTF)) +
mp1 geom_point(aes(shape = MacroEnv, color = MacroEnv)) +
geom_smooth(data = x1, method = "lm", se = F, color = "black", lty = 3) +
geom_smooth(data = x2, method = "lm", se = F, color = "black") +
scale_y_continuous(sec.axis = dup_axis(~ 1/., name = NULL, breaks = c(35,50,100,150)),
trans = "reverse", breaks = c(0.01,0.02,0.03), limits = c(mymax, mymin)) +
scale_x_continuous(breaks = c(11,13,15,17,19,21)) +
scale_shape_manual(name = "Macro-environment", values = c(16,15,17)) +
scale_color_manual(name = "Macro-environment", values = myfills) +
+
theme_AGL labs(title = figlab, y = "1 / DTF",
x = expression(paste("Temperature (", degree, "C)", sep="")))
# Plot (b) 1/f = a + cP
<- ggplot(xx, aes(x = P_mean, y = RDTF)) +
mp2 geom_point(aes(shape = MacroEnv, color = MacroEnv)) +
geom_smooth(data = x1, method = "lm", se = F, color = "black", lty = 3) +
geom_smooth(data = x3, method = "lm", se = F, color = "black") +
scale_y_continuous(sec.axis = dup_axis(~ 1/., name="DTF", breaks = c(35,50,100,150)),
trans = "reverse", breaks = c(0.01,0.02,0.03), limits = c(mymax, mymin)) +
scale_x_continuous(breaks = c(11,12,13,14,15,16)) +
scale_shape_manual(name = "Macro-environment", values = c(16,15,17)) +
scale_color_manual(name = "Macro-environment", values = myfills) +
+
theme_AGL labs(title = "", y = NULL,x = "Photoperiod (hours)")
#
<- ggarrange(mp1, mp2, ncol = 2, common.legend = T, legend = "bottom")
mp[[i]] ggsave(paste0("Additional/Entry_TP/TP_Entry_", str_pad(i, 3, pad = "0"), ".png"),
width = 8, height = 4, dpi = 600)
mp[[i]],
}pdf("Additional/pdf_TP.pdf", width = 8, height = 4)
for(i in 1:324) { print(mp[[i]]) }
dev.off()
Additional Figures: PhotoThermal Plane
# Prep data
<- rr %>%
xx filter(!is.na(RDTF)) %>%
left_join(select(ff, Expt, T_mean, P_mean, MacroEnv), by = "Expt")
# Create plotting function
<- function(x, i) {
ggPTplane <- x %>% filter(Entry == i)
x1 <- x1$T_mean
x <- x1$P_mean
y <- x1$RDTF
z <- lm(z ~ x + y)
fit # Create PhotoThermal plane
<- predict(fit)
fitp <- 12
grid.lines <- seq(min(x), max(x), length.out = grid.lines)
x.p <- seq(min(y), max(y), length.out = grid.lines)
y.p <- expand.grid(x = x.p, y = y.p)
xy <- matrix(predict(fit, newdata = xy), nrow = grid.lines, ncol = grid.lines)
z.p <- plyr::mapvalues(x1$Expt, names_Expt, c(rep(16,6), rep(15,6), rep(17,6))) %>%
pchs as.character() %>% as.numeric()
# Plot with regression plane
par(mar=c(1.5, 2.5, 1.5, 0.5))
scatter3D(x, y, z, pch = pchs, cex = 2, main = unique(x1$Name),
col = alpha(c("darkgreen", "darkorange3", "darkblue"),0.5),
colvar = as.numeric(x1$MacroEnv), colkey = F, col.grid = "gray90",
bty = "u", theta = 40, phi = 25, ticktype = "detailed",
cex.lab = 1, cex.axis = 0.5, zlim = c(0.005,0.03),
xlab = "Temperature", ylab = "Photoperiod", zlab = "1 / DTF",
surf = list(x = x.p, y = y.p, z = z.p, col = "black", facets = NA, fit = fitp) )
}# Plot each Entry
for (i in 1:324) {
png(paste0("Additional/Entry_3D/3D_Entry_", str_pad(i, 3, pad = "0"), ".png"),
width = 1000, height = 1000, res = 300)
ggPTplane(xx, i)
dev.off()
}# Create PDF
pdf("Additional/pdf_3D.pdf")
par(mar=c(1.5, 2.5, 1.5, 0.5))
for (i in 1:324) {
ggPTplane(xx, i)
}dev.off()
#
# Plot ILL 5888 & ILL 4400 & Laird
<- xx %>% mutate(Name = gsub(" AGL", "", Name))
xx for (i in c(235, 94, 128)) {
png(paste0("Additional/Temp/3D_Entry_", str_pad(i, 3, pad = "0"), ".png"),
width = 2000, height = 2000, res = 600)
ggPTplane(xx, i)
dev.off()
}# Create animation
<- read.csv("data/model_t+p_coefs.csv") %>% arrange(b, c)
xx <- list.files("Additional/Entry_3D")[xx$Entry]
lf <- image_read(paste0("Additional/Entry_3D/", lf))
mp <- image_animate(mp, fps = 10)
animation image_write(animation, "Additional/Animation_3D.gif")
Supplemental Figure 4: Regressions
Figure S4: Effects of mean temperature and photoperiod on the rate of progress towards flowering (1 / DTF) in three contrasting selected genotypes. (a) Effect of temperature on 1 / DTF, (b) effect of photoperiod on 1 / DTF, and (c) effect of temperature and photoperiod on 1 / DTF modelled using equation 1. For (a) and (b), solid lines represent regressions among locations of relatively constant photoperiod or temperature, respectively, while dotted lines indicate a break in the assumption of constant photoperiod or temperature, respectively, across environments (see Figure 1). (d) Scaled DTF (1-5) of each genotype (grey lines) across all site-years with ILL5888, PI 420925 LSP and Laird highlighted according to their corresponding cluster group, 1, 5 and 8 respectively. ILL 5888 is an early maturing, genotype from Bangladesh. PI 420925 LSP is a landrace from Jordan with medium maturity. Laird is a late maturing, Canadian cultivar.
# Prep data
<- alpha(c("darkgreen", "darkorange3", "darkblue"), 0.5)
myfills <- c("ILL 5888 AGL", "PI 420925 LSP AGL", "Laird AGL")
yy <- rr %>%
xx filter(Name %in% yy, !is.na(DTF)) %>%
left_join(select(ff, Expt, MacroEnv, T_mean, P_mean), by = "Expt") %>%
mutate(Name = gsub(" AGL", "", Name), myfill = MacroEnv,
Name = factor(Name, levels = gsub(" AGL", "", yy)))
<- xx %>% filter(MacroEnv != "South Asia")
x1 <- xx %>% filter(MacroEnv != "Temperate")
x2 <- xx %>% filter(MacroEnv != "Mediterranean")
x3 # Plot (a) 1/f = a + bT
<- ggplot(xx, aes(x = T_mean, y = RDTF)) +
mp1 geom_point(aes(shape = MacroEnv, color = MacroEnv)) +
geom_smooth(data = x1, method = "lm", se = F, color = "black", lty = 3) +
geom_smooth(data = x2, method = "lm", se = F, color = "black") +
scale_y_continuous(trans = "reverse", breaks = c(0.01,0.02,0.03),
sec.axis = dup_axis(~ 1/., name = "DTF", breaks = c(35,50,100,150))) +
scale_x_continuous(breaks = c(11,13,15,17,19,21)) +
scale_shape_manual(name = "Macro-environment", values = c(16,15,17)) +
scale_color_manual(name = "Macro-environment", values = myfills) +
+ facet_grid(. ~ Name) +
theme_AGL theme(axis.title.y = element_text(size = 9),
plot.margin = unit(c(0,0,0,0), "cm")) +
guides(colour = guide_legend(override.aes = list(size = 3))) +
labs(y = "1 / DTF", x = expression(paste("Temperature (", degree, "C)", sep = "")))
# Plot (b) 1/f = a + cP
<- ggplot(xx, aes(x = P_mean, y = RDTF)) +
mp2 geom_point(aes(shape = MacroEnv, color = MacroEnv)) +
geom_smooth(data = x1, method = "lm", se = F, color = "black", lty = 3) +
geom_smooth(data = x3, method = "lm", se = F, color = "black") +
scale_y_continuous(trans = "reverse", breaks = c(0.01,0.02,0.03),
sec.axis = dup_axis(~ 1/., name="DTF", breaks = c(35,50,100,150))) +
scale_x_continuous(breaks = c(11,12,13,14,15,16)) +
scale_shape_manual(name = "Macro-environment", values = c(16,15,17)) +
scale_color_manual(name = "Macro-environment", values = myfills) +
+ facet_grid(. ~ Name) +
theme_AGL theme(axis.title.y = element_text(size = 9),
plot.margin = unit(c(0,0,0,0), "cm")) +
guides(colour = guide_legend(override.aes = list(size = 3))) +
labs(y = "1 / DTF",x = "Photoperiod (h)")
# Append (a) and (b)
<- ggarrange(mp1, mp2, ncol = 1, common.legend = T, legend = "bottom")
mp ggsave("Additional/Temp/Temp_SF04_1.png", mp, width = 6, height = 3.75, dpi = 600, bg = "white")
# Append (c)s
<- image_read("Additional/Temp/3D_Entry_094.png")
im1 <- image_read("Additional/Temp/3D_Entry_235.png")
im2 <- image_read("Additional/Temp/3D_Entry_128.png")
im3 <- image_append(c(im1, im2, im3)) %>% image_scale("3600x")
im image_write(im, "Additional/Temp/Temp_SF04_2.png")
# Prep data
<- dd %>%
xx filter(Name %in% yy) %>%
mutate(Name = gsub(" AGL", "", Name),
Name = factor(Name, levels = gsub(" AGL", "", yy)))
# Plot D)
<- ggplot(dd, aes(x = ExptShort, y = DTF2_scaled, group = Name)) +
mp3 geom_line(color = "black", alpha = 0.1) +
geom_vline(xintercept = 6.5, lty = 2) +
geom_vline(xintercept = 12.5, lty = 2) +
geom_line(data = xx, aes(color = Name), size = 2) +
scale_color_manual(name = "Genotype", values = colors[c(1,5,8)]) +
+ labs(y = "DTF (scaled 1-5)", x = NULL) +
theme_AGL theme(legend.position = "bottom",
legend.margin = unit(c(0,0,0,0), "cm"),
plot.margin = unit(c(0,0.2,0,0.5), "cm"),
panel.grid = element_blank(),
axis.title.y = element_text(size = 9),
axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("Additional/Temp/Temp_SF04_3.png", mp3, width = 6, height = 2.5, dpi = 600)
# Append (a), (b), C) and D)
<- image_read("Additional/Temp/Temp_SF04_1.png") %>%
im1 image_annotate("(a)", size = 100, location = "+0+0") %>%
image_annotate("(b)", size = 100, location = "+0+1000")
<- image_read("Additional/Temp/Temp_SF04_2.png") %>%
im2 image_annotate("(c)", size = 100)
<- image_read("Additional/Temp/Temp_SF04_3.png") %>%
im3 image_annotate("(d)", size = 100)
<- image_append(c(im1, im2, im3), stack = T)
im image_write(im, "Supplemental_Figure_04.png")
Modeling DTF - Functions
# Create functions
# Plot Observed vs Predicted
<- function(x, title = NULL, type = 1,
ggModel1 mymin = min(c(x$DTF,x$Predicted_DTF)) - 2,
mymax = max(c(x$DTF,x$Predicted_DTF)) + 2 ) {
<- x %>% mutate(Flowered = ifelse(is.na(DTF), "Did not Flower", "Flowered"))
x # Prep data
if(type == 1) {
<- "DTF"; myy <- "Predicted_DTF"
myx <- x %>% filter(!is.na(DTF))
x
}if(type == 2) {
<- "RDTF"; myy <- "Predicted_RDTF"
myx <- x %>% filter(!is.na(RDTF))
x
}<- colors_Expt[names_Expt %in% unique(x$Expt)]
myPal <- round(modelR2(x = x[,myx], y = x[,myy]), 3)
r2 # Plot
<- ggplot(x) +
mp geom_point(aes(x = get(myx), y = get(myy), color = Expt)) +
geom_abline() +
geom_label(x = mymin, y = mymax, hjust = 0, vjust = 1, parse = T,
label = paste("italic(R)^2 == ", r2)) +
scale_x_continuous(limits = c(mymin, mymax), expand = c(0, 0)) +
scale_y_continuous(limits = c(mymin, mymax), expand = c(0, 0)) +
scale_color_manual(name = NULL, values = myPal) +
guides(colour = guide_legend(override.aes = list(size = 2))) +
theme_AGL if(type == 1) {
<- mp + labs(y = "Predicted DTF (days)", x = "Observed DTF (days)")
mp
}if(type == 2) {
<- mp +
mp scale_x_reverse(limits = c(mymax, mymin), expand = c(0, 0)) +
scale_y_reverse(limits = c(mymax, mymin), expand = c(0, 0)) +
labs(y = "Predicted RDTF (1/days)", x = "Observed RDTF (1/days)")
}if(!is.null(title)) { mp <- mp + labs(title = title) }
mp
}# Facets by Expt
<- function(x, myX = "DTF", myY = "Predicted_DTF", title = NULL,
ggModel2 x1 = 30, x2 = 30, y1 = 145, y2 = 120, legend.pos = "bottom") {
# Prep data
<- read.csv("data/data_pca_results.csv") %>%
pca select(Entry, Cluster) %>%
mutate(Cluster = factor(Cluster))
<- x %>%
x filter(!is.na(get(myX))) %>%
left_join(pca, by = "Entry")
<- x %>%
xf group_by(Expt) %>%
summarise(Mean = mean(DTF)) %>%
ungroup() %>%
mutate(r2 = NULL, RMSE = NULL)
for(i in 1:nrow(xf)) {
<- x %>% filter(Expt == xf$Expt[i])
xi "r2"] <- round(modelR2(x = xi[,myX], y = xi[,myY]), 2)
xf[i,"RMSE"] <- round(modelRMSE(x = xi[,myX], y = xi[,myY]), 1)
xf[i,
}# Plot
<- ggplot(x, aes(x = get(myX), y = get(myY))) +
mp geom_point(aes(color = Cluster), size = 0.75, alpha = 0.7) + geom_abline() +
geom_text(x = x1, y = y1, color = "black", hjust = 0, vjust = 0, size = 3,
aes(label = paste("RMSE = ", RMSE, sep = "")), data = xf) +
geom_text(x = x2, y = y2, color = "black", hjust = 0, vjust = 0, size = 3,
aes(label = paste("italic(R)^2 == ", r2)), parse = T, data = xf) +
facet_wrap(Expt ~ ., ncol = 6, labeller = label_wrap_gen(width = 17)) +
scale_x_continuous(limits = c(min(x[,myX]), max(x[,myX]))) +
scale_y_continuous(limits = c(min(x[,myX]), max(x[,myX]))) +
scale_color_manual(values = colors) +
+
theme_AGL theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = legend.pos,
legend.margin = unit(c(0,0,0,0), "cm"),
panel.grid = element_blank()) +
guides(colour = guide_legend(nrow = 1, override.aes = list(size = 3))) +
labs(y = "Predicted DTF (days)", x = "Observed DTF (days)")
if(!is.null(title)) { mp <- mp + labs(title = title) }
mp
}# R^2 function
<- function(x, y) {
modelR2 1 - ( sum((x - y)^2, na.rm = T) / sum((x - mean(x, na.rm = T))^2, na.rm = T))
}# RMSE function
<- function(x, y) {
modelRMSE sqrt(sum((x-y)^2, na.rm = T) / length(x))
}
\(R^2=1-\frac{SS_{residuals}}{SS_{total}}=1-\frac{\sum (x-y)^2}{\sum (x-\bar{x})}\)
\(RMSE=\frac{\sum (y-x)^2}{n}\)
Modeling DTF (T + P) - All Site-years
###########################
# 1/f = a + bT + cP (All) #
###########################
# Prep data
<- rr %>%
xx filter(!is.na(RDTF)) %>%
left_join(select(ff, Expt, T_mean, P_mean), by = "Expt") %>%
select(Plot, Entry, Name, Rep, Expt, ExptShort, T_mean, P_mean, RDTF, DTF)
<- NULL
mr <- NULL
md <- ldp %>% select(Entry, Name) %>%
mc mutate(a = NA, b = NA, c = NA, RR = NA, Environments = NA, aP = NA, bP = NA, cP = NA)
# Model
for(i in 1:324) {
# Prep data
<- xx %>% filter(Entry == i)
xri <- xri %>%
xdi group_by(Entry, Name, Expt, ExptShort) %>%
summarise_at(vars(DTF, RDTF, T_mean, P_mean), funs(mean), na.rm = T) %>%
ungroup()
# Train Model
<- lm(RDTF ~ T_mean + P_mean, data = xri)
mi # Predict DTF
<- xri %>% mutate(Predicted_RDTF = predict(mi),
xri Predicted_DTF = 1 / predict(mi))
<- xdi %>% mutate(Predicted_RDTF = predict(mi, newdata = xdi),
xdi Predicted_DTF = 1 / predict(mi, newdata = xdi))
# Save to table
<- bind_rows(mr, xri)
mr <- bind_rows(md, xdi)
md # Save coefficients
c("a","b","c")] <- mi$coefficients
mc[i,# Calculate rr and # of environments used
"RR"] <- 1 - sum((xri$DTF - xri$Predicted_DTF)^2, na.rm = T) /
mc[i,sum((xri$Predicted_DTF - mean(xri$DTF, na.rm = T))^2, na.rm = T)
"Environments"] <- length(unique(xri$Expt[!is.na(xri$DTF)]))
mc[i,"aP"] <- summary(mi)[[4]][1,4]
mc[i,"bP"] <- summary(mi)[[4]][2,4]
mc[i,"cP"] <- summary(mi)[[4]][3,4]
mc[i,
}<- mr %>% mutate(Expt = factor(Expt, levels = names_Expt))
mr <- md %>% mutate(Expt = factor(Expt, levels = names_Expt))
md # Save Results
write.csv(mr, "data/model_t+p.csv", row.names = F)
write.csv(md, "data/model_t+p_d.csv", row.names = F)
write.csv(mc, "data/model_t+p_coefs.csv", row.names = F)
#
# Plot Each Entry
<- list()
mp for(i in 1:324) {
<- ggModel1(mr %>% filter(Entry == i), paste("Entry", i, "| DTF"),
mp1 mymin = min(c(mr$Predicted_DTF, mr$DTF), na.rm = T),
mymax = max(c(mr$Predicted_DTF, mr$DTF), na.rm = T))
<- ggModel1(mr %>% filter(Entry == i), paste("Entry", i, "| RDTF"), type = 2,
mp2 mymin = min(c(mr$Predicted_RDTF, mr$RDTF)) - 0.001,
mymax = max(c(mr$Predicted_RDTF, mr$RDTF)) + 0.001)
<- ggarrange(mp1, mp2, ncol = 2, common.legend = T, legend = "right")
mp[[i]] <- paste0("Additional/entry_model/model_entry_", str_pad(i, 3, pad = "0"), ".png")
fname ggsave(fname, mp[[i]], width = 10, height = 4.5, dpi = 600)
}pdf("Additional/pdf_model.pdf", width = 10, height = 4.5)
for (i in 1:324) { print(mp[[i]]) }
dev.off()
# Prep data
<- read.csv("data/model_t+p_d.csv") %>%
xx mutate(Expt = factor(Expt, levels = names_Expt))
# Plot Observed vs Predicted
<- ggModel1(xx, title = "Model = T + P")
mp ggsave("Additional/Model/Model_1_1.png", mp, width = 7, height = 5, dpi = 600)
# Plot Observed vs Predicted
<- ggModel2(xx, title = "Model = T + P")
mp ggsave("Additional/Model/Model_2_1.png", mp, width = 8, height = 5.5, dpi = 600)
# Plot Observed vs Predicted for Temperate Locations
<- c("Ro16","Ro17","Su16","Su17","Su18","Us18")
myexpts <- ggModel2(xx %>% filter(ExptShort %in% myexpts))
mp ggsave("Additional/Model/Model_3_1.png", mp, width = 8, height = 3, dpi = 600)
Modeling DTF (T x P) - All Site-years
###########################
# 1/f = a + bT + cP (All) #
###########################
# Prep data
<- rr %>%
xx filter(!is.na(RDTF)) %>%
left_join(select(ff, Expt, T_mean, P_mean), by = "Expt") %>%
select(Plot, Entry, Name, Rep, Expt, ExptShort, T_mean, P_mean, RDTF, DTF)
<- NULL
mr <- NULL
md <- ldp %>% select(Entry, Name) %>%
mc mutate(a = NA, b = NA, c = NA, d = NA, RR = NA, Environments = NA,
aP = NA, bP = NA, cP= NA, dP = NA)
# Model
for(i in 1:324) {
# Prep data
<- xx %>% filter(Entry == i)
xri <- xri %>%
xdi group_by(Entry, Name, Expt, ExptShort) %>%
summarise_at(vars(DTF, RDTF, T_mean, P_mean), funs(mean), na.rm = T) %>%
ungroup()
# Train Model
<- lm(RDTF ~ T_mean * P_mean, data = xri)
mi # Predict DTF
<- xri %>% mutate(Predicted_RDTF = predict(mi),
xri Predicted_DTF = 1 / predict(mi))
<- xdi %>% mutate(Predicted_RDTF = predict(mi, newdata = xdi),
xdi Predicted_DTF = 1 / predict(mi, newdata = xdi))
# Save to table
<- bind_rows(mr, xri)
mr <- bind_rows(md, xdi)
md # Save coefficients
c("a","b","c","d")] <- mi$coefficients
mc[i,# Calculate rr and # of environments used
"RR"] <- 1 - sum((xri$DTF - xri$Predicted_DTF)^2, na.rm = T) /
mc[i,sum((xri$Predicted_DTF - mean(xri$DTF, na.rm = T))^2, na.rm = T)
"Environments"] <- length(unique(xri$Expt[!is.na(xri$DTF)]))
mc[i,"aP"] <- summary(mi)[[4]][1,4]
mc[i,"bP"] <- summary(mi)[[4]][2,4]
mc[i,"cP"] <- summary(mi)[[4]][3,4]
mc[i,"dP"] <- summary(mi)[[4]][4,4]
mc[i,
}<- mr %>% mutate(Expt = factor(Expt, levels = names_Expt))
mr <- md %>% mutate(Expt = factor(Expt, levels = names_Expt))
md # Save Results
write.csv(mr, "data/model_txp.csv", row.names = F)
write.csv(md, "data/model_txp_d.csv", row.names = F)
write.csv(mc, "data/model_txp_coefs.csv", row.names = F)
# Prep data
<- read.csv("data/model_txp_d.csv") %>%
xx mutate(Expt = factor(Expt, levels = names_Expt))
# Plot Observed vs Predicted
<- ggModel1(xx, title = "Model = T x P")
mp ggsave("Additional/Model/Model_1_2.png", mp, width = 7, height = 5, dpi = 600)
# Plot Observed vs Predicted
<- ggModel2(xx, title = "Model = T x P")
mp ggsave("Additional/Model/Model_2_2.png", mp, width = 8, height = 5.5, dpi = 600)
Supplemental Table 3: Model Constants
Table S3: Values of the constants derived from equations 1 and 2 using data from all site-years, for each of the genotypes used in this study.
# Prep data
<- read.csv("data/model_t+p_coefs.csv")
x1 <- read.csv("data/model_txp_coefs.csv")
x2 <- bind_rows(x1, x2) %>%
xx arrange(Entry) %>%
select(Entry, Name, a, b, c, d, RR, Environments,
a_p.value=aP, b_p.value=bP, c_p.value=cP, d_p.value=dP)
# Save
write.csv(xx, "Supplemental_Table_03.csv", na = "", row.names = F)
Additional Figure 10: significant T x P interactions
# Prep data
<- read.csv("data/model_t+p_coefs.csv")
x1 <- read.csv("data/model_txp_coefs.csv")
x2 <- bind_rows(x1, x2) %>%
xx arrange(Entry) %>%
select(Entry, Name, a, b, c, d, RR)
<- x2 %>% filter(dP < 0.05) %>% pull(Entry)
entries <- xx %>% filter(Entry %in% entries)
xx <- xx %>%
xx arrange(RR) %>%
mutate(Entry = factor(Entry, levels = unique(Entry)),
Model = ifelse(is.na(d), "T + P", "T x P"))
# Plot
<- ggplot(xx, aes(x = Entry, y = RR, fill = Model)) +
mp geom_bar(stat = "identity", position = "dodge") +
scale_y_continuous(breaks = seq(0,1,0.1), minor_breaks = seq(0,1,0.01)) +
+
theme_AGL theme(legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("Additional/Additional_Figure_10.png", mp,
width = 7, height = 5, dpi = 600)
#
length(entries)
[1] 30
Supplemental Figure 5: Model T + P vs. T x P
Figure S5: Comparison of observed and predicted values for days from sowing to flowering using (a) equation 1 and (b) equation 2.
# Prep data
<- read.csv("data/model_t+p_d.csv") %>%
xx mutate(Expt = factor(Expt, levels = names_Expt))
# Plot Observed vs Predicted
<- ggModel1(xx, title = expression(paste("(a) ", italic(" 1 / f = a + bT + cP"))))
mp1 # Prep data
<- read.csv("data/model_txp_d.csv") %>%
xx mutate(Expt = factor(Expt, levels = names_Expt))
# Plot Observed vs Predicted
<- ggModel1(xx, title = expression(paste("(b) ", italic(" 1 / f = a + bT + cP + dTP"))))
mp2 # Append
<- ggarrange(mp1, mp2, ncol = 2, common.legend = T, legend = "right")
mp ggsave("Supplemental_Figure_05.png", mp, width = 12, height = 5, dpi = 600, bg = "white")
Additional Figure 11: Constants
# Prep data
<- read.csv("data/data_pca_results.csv")
pca <- read.csv("data/model_t+p_coefs.csv") %>% mutate(Model = "T+P")
x1 <- read.csv("data/model_txp_coefs.csv") %>% mutate(Model = "TxP")
x2 <- bind_rows(x1, x2) %>%
xx gather(Coef, Value, a, b, c, d) %>%
left_join(pca, by = c("Entry","Name")) %>%
mutate(Cluster = factor(Cluster), Coef = factor(Coef))
# Plot
<- ggplot(xx %>% filter(Model=="T+P", !is.na(Value)),
mp1 aes(x = Cluster, y = Value * 10000, color = Cluster)) +
geom_quasirandom() +
+
theme_AGL theme(strip.text.y = element_text(face = "italic")) +
facet_grid(Coef ~ Model, scales = "free_y", drop = F) +
scale_color_manual(values = colors) + labs(y = "x 10,000")
<- ggplot(xx %>% filter(Model=="TxP"),
mp2 aes(x = Cluster, y = Value * 10000, color = Cluster)) +
geom_quasirandom() +
facet_grid(Coef ~ Model, scales = "free_y") +
+
theme_AGL theme(strip.text.y = element_text(face = "italic"),
panel.grid.major.x = element_blank()) +
scale_color_manual(values = colors) + labs(y = "x 10,000")
<- ggarrange(mp1, mp2, ncol = 2, legend = "none")
mp ggsave("Additional/Additional_Figure_11.png", mp, width = 8, height = 5, dpi = 600)
Additional Figure 12: Coefficient p-values
# Prep data
<- read.csv("data/data_pca_results.csv") %>% mutate(Cluster = factor(Cluster))
pca <- read.csv("data/model_t+p_coefs.csv") %>% mutate(Model = "T+P")
x1 <- read.csv("data/model_txp_coefs.csv") %>% mutate(Model = "TxP")
x2 <- bind_rows(x1, x2) %>%
xx gather(Coef, Value, aP,bP,cP,dP) %>%
left_join(pca, by = "Entry") %>%
mutate(Cluster = factor(Cluster), Coef = factor(Coef))
# Plot
<- ggplot(xx %>% filter(Model=="T+P", !is.na(Value)),
mp1 aes(x = Cluster, y = Value, color = Cluster)) +
geom_quasirandom() +
facet_grid(Coef ~ Model, scales = "free_y", drop = F) +
scale_color_manual(values = colors) +
+
theme_AGL theme(strip.text.y = element_text(face = "italic"),
panel.grid.major.x = element_blank()) +
labs(y = expression(paste(italic("p"), " value")))
<- ggplot(xx%>%filter(Model=="TxP"),
mp2 aes(x = Cluster, y = Value, color = Cluster)) +
geom_quasirandom() +
facet_grid(Coef ~ Model, scales = "free_y") +
+
theme_AGL theme(strip.text.y = element_text(face = "italic"),
panel.grid.major.x = element_blank()) +
scale_color_manual(values = colors) +
labs(y = expression(paste(italic("p"), " value")))
<- ggarrange(mp1,mp2,ncol=2, legend = "none")
mp ggsave("Additional/Additional_Figure_12.png", mp, width = 8, height = 5, dpi = 600)
Additional Figure 13: p-values b c
# Prep data
<- read.csv("data/data_pca_results.csv") %>%
pca mutate(Cluster = factor(Cluster))
<- read.csv("data/model_t+p_coefs.csv") %>%
xx mutate(Sig = factor(ifelse(bP > 0.0001, "Sig","Less Sig"))) %>%
select(Entry, Sig, p.value=bP, b) %>%
gather(Trait, Value, p.value, b) %>%
left_join(pca, by = "Entry")
<- xx %>% filter(Sig == "Sig", Trait == "p.value")
x1 # Plot (a)
<- ggplot(xx, aes(x = "", y = Value, group = "")) +
mp1 geom_quasirandom(aes(color = Cluster)) +
geom_point(data = x1, size = 2.5, pch = 25,
color = "black", fill = "darkblue") +
geom_text_repel(data = x1, aes(label = Entry)) +
facet_wrap(.~Trait, scales = "free_y") +
scale_color_manual(values = colors) +
guides(colour = guide_legend(nrow = 1, override.aes = list(size = 3))) +
+
theme_AGL theme(panel.grid.major.x = element_blank()) +
labs(title = "(a)", x = NULL, y = NULL)
# Prep data
<- read.csv("data/model_t+p_coefs.csv") %>%
xx mutate(Sig = factor(ifelse(cP > 0.000001, "Sig","Less Sig"))) %>%
select(Entry, Sig, p.value=cP, c) %>%
gather(Trait, Value, p.value, c) %>%
left_join(pca, by = "Entry")
<- xx %>% filter(Sig == "Sig")
x1 # Plot B)
<- ggplot(xx, aes(x = "", y = Value, group = "")) +
mp2 geom_quasirandom(aes(color = Cluster)) +
geom_point(data = x1, size = 2.5, pch = 25,
color = "black", fill = "darkred") +
geom_text_repel(data = x1 %>% filter(Trait =="p.value"), aes(label = Entry)) +
facet_wrap(.~Trait, scales = "free_y") +
scale_color_manual(values = colors) +
guides(colour = guide_legend(nrow = 1, override.aes = list(size = 3))) +
+
theme_AGL theme(panel.grid.major.x = element_blank()) +
labs(title = "(b)", x = NULL, y = NULL)
# Append (a) and (b)
<- ggarrange(mp1, mp2, ncol = 2, common.legend = T, legend = "bottom")
mp ggsave("Additional/Additional_Figure_13.png", mp, width = 8, height = 4, dpi = 600, bg = "white")