RDA_Genomic_offset
Thomas Francisco
2024-07-18
meta_data_pop <- read.csv("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Populations/taxus_sample_29pop.csv",h=T,sep=";",dec=",")
#alphabetic order
meta_data_pop_order <- meta_data_pop[order(meta_data_pop$Population),]1 Introduction
This script calculates the Euclidean distances between observed and predicted population scores along the retained RDA axes, and normalises the values using min-max scaling to a zero-to-one scale. (later called genomic discrepancy index) and the genomic offset index using the RDA method. The genomic discrepancy index is a metric we introduce to visualise how far the observed genomic composition of populations deviates from that predicted by the GEA model under the same climatic conditions. Populations are unlikely to have exactly the genomic composition predicted by the model’s linear local adaptation relationships. Deviations may suggest that the degree of local adaptation is lower than expected, or that adaptation relies on a different genetic basis than captured by the model. Possible explanations include neutral processes such as drift (fixing potentially adaptive alleles), restricted gene flow (leading to absence of adaptive variants in some populations), or populations experiencing adaptation lag or pre-adaptation to future climates. However, alternative explanations such as conditional neutrality of alleles or model misspecification should also be considered.
The genomic offset index, first proposed by Fitzpatrick & Keller (2015), quantifies the amount of genomic change required for populations (or individuals) to maintain fitness under future climatic conditions. Specifically, genomic offset measures the distance between the estimated genomic composition under current/reference climates and that predicted under future climates, focusing on SNPs potentially associated with adaptation. To compute genomic offset, potentially climate-adaptive SNPs could first be identified (alternatively, all SNPs can be used under the assumption that local adaptation is widespread). GEA methods are then applied to establish associations between allelic variation and climatic variables. Using these relationships, we can estimate past/present and future genomic compositions from the corresponding climate data, effectively interpolating/extrapolating the GEA models in space and time. The genomic offset is then quantified as the Euclidean distance between genomic compositions under the reference and future climates.
2 Data
2.1 Past/present and future climatic data
#climatic/structure/IBD data
load("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/GEA_new_var/variance_partitioning/Climatic_data_RDA_pRDA.Rdata") Here we will also load the raster data not scaled but with the mean and sd values of the past climatic data to scaled the data together to follow Capblancq and Forester (2021).
list_models_clim<- c("GFDL_ESM4","IPSL_CM6A_LR","MPI_ESM1_2_HR","MRI_ESM2_0","UKESM1_0_LL")
#load past and future raster data
climatic_variables_to_keep <-c("bio1", "bio12", "bio15","bio2","bio4","bio9") #subset of retained clim var, see script processing clim data
list_raster <- c("past_climatic_data_raster","present_climatic_data_raster","future_climatic_data_raster")
for(i in 1:length(list_models_clim)){
model <- list_models_clim[i]
for(x in 1:length(list_raster)){
name <- c("past","present","future")
name_clim <- name[x]#name final raster
var <- list_raster[x]#name of the raster
if(var == "future_climatic_data_raster"){
var= paste0("future_climatic_data_raster_",model)
load(paste0("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Genomic_offset/RDA/",var,".Rdata")) #load data
raster_clim <- get(var)[[climatic_variables_to_keep]]#keep only var in clim_variables to keep
names(raster_clim)=c("Annual_Tc","Annual_P","P_Seasonality","Diurnal_range_Tc","Tc_Seasonality","Tc_driest_quarter") #change the name to match the names of the scaled and center values
assign(paste0("raster_",name_clim,"_clim_",model),raster_clim)
} else {
load(paste0("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Genomic_offset/RDA/",var,".Rdata")) #load data
raster_clim <- get(var)[[climatic_variables_to_keep]]#keep only var in clim_variables to keep
names(raster_clim)=c("Annual_Tc","Annual_P","P_Seasonality","Diurnal_range_Tc","Tc_Seasonality","Tc_driest_quarter") #change the name to match the names of the scaled and center values
assign(paste0("raster_",name_clim,"_clim"),raster_clim)
}
}
}We also need the scale and center values:
load("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Climatic_data/new_selection/scale_env_value_new_cli.Rdata")
load("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/Climatic_data/new_selection/center_env_value_new_cli.Rdata")list_name <- c("past","present","future")
coords <- data.frame(apply(meta_data_pop_order[,c(5:4)], 2, as.numeric))#we need to have longitude then latitude
for(i in 1:length(list_models_clim)){
model <- list_models_clim[i]
for(x in 1: length(list_name)){
name <- list_name[x]
if(name == "future"){
raster_data <- get(paste0("raster_",name,"_clim_",model))
clim_data <- raster::extract(raster_data, coords)
clim_data_order <- data.frame(clim_data[,c(1,4,5,6,2,3)])#put it in the same order as center_env and scale object
# Standardization of the environmental variables
clim_data_scale <- as.data.frame(scale(clim_data_order, center=center_env_value_new_cli, scale=scale_env_value_new_cli))
assign(paste0("clim_df_",name,"_scale_",model),clim_data_scale)
} else {
raster_data <- get(paste0("raster_",name,"_clim"))
clim_data <- raster::extract(raster_data, coords)
clim_data_order <- data.frame(clim_data[,c(1,4,5,6,2,3)])#put it in the same order as center_env and scale object
# Standardization of the environmental variables
clim_data_scale <- as.data.frame(scale(clim_data_order, center=center_env_value_new_cli, scale=scale_env_value_new_cli))
assign(paste0("clim_df_",name,"_scale"),clim_data_scale)
}
}
}2.2 Genomic data
The genomic data set includes allele frequency corrected for MAC in the 29 populations.
#genomic data
load("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/genetic_new/data_allelic_frequencies_29pop_adapcon_gentree_475_8616.Rdata")
genomic_matrix <- data_allelic_frequencies_29pop_adapcon_gentree_475_8616 2.3 Outliers sets
We created 5 set of outliers:
2.3.1 Random SNPs
There are 2 sets of randoms SNPs:
- 1 neutral set with loci randomly selected from a set of SNPs not including the candidate ones:
load("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GEA_new_var/outliers/random_neutral_set_SNPs_T_adapcon_gentree_bis.Rdata")
list_random_SNPs <- colnames(random_neutral_set_SNPs_T_adapcon_gentree_bis)
#write.csv(list_random_SNPs,"C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GEA_new_var/outliers/list_random_SNPs.csv")- 1 neutral set with loci selected at random, while maintaining the same number of SNPs for each class of allelic frequencies as in the outlier set (LC).
load("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GEA_new_var/outliers/random_set_taxus.Rdata")
list_random_SNPs_same_AF <- random_set_taxus$name_snps
#write.csv(list_random_SNPs_same_AF,"C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GEA_new_var/outliers/list_random_SNPs_same_AF.csv")2.3.2 Less conservatives thresholds
#set of less conservative thresholds
load("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GEA_new_var/outliers/outliers_set_final_overlapping_no_LD_LC_new_var.Rdata")2.3.3 More conservatives thresholds
load("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GEA_new_var/outliers/outliers_set_final_overlapping_no_LD_new_var.Rdata")2.3.4 Overlapping snps LC with CG populations
load("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/genetic_new/gen_matrix_imp_CG.Rdata")
CG_common_set <- intersect(outliers_set_final_overlapping_no_LD_LC_new_var, colnames(gen_matrix_imp_CG))2.3.5 another set of random SNPs following the allelic frequency of the outliers
load("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GEA_new_var/outliers/random_set_taxus_not_overlap_both_dataset_V2.Rdata")
list_random_SNPs_same_AF_V2 <- random_set_taxus_not_overlap_both_dataset_V2$name_snps2.3.6 all the outliers
load("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GEA_new_var/outliers/unique_outliers.Rdata")3 Calculate GEA relationship
Our genomic offset analysis is based on the relationship between SNPs
and the climatic variables identified by GEA methods. In this script, we
will use redundancy analysis (RDA). First, we need to perform a GEA
analysis on the SNP set to identify the relationships that can then be
interpolated or extrapolated.
To achieve this, we ran an RDA on the SNP and the retained climatic
variables.
3.1 Models
- all SNPs
RDA_all <- rda(formula = genomic_matrix ~ Annual_Tc+Diurnal_range_Tc+Tc_Seasonality+Tc_driest_quarter+Annual_P+P_Seasonality, data = Climatic_data_RDA_pRDA, scale=F)
RsquareAdj(RDA_all)## $r.squared
## [1] 0.3262662
##
## $adj.r.squared
## [1] 0.1425206- random set of SNPs
RDA_random <- rda(formula = genomic_matrix[list_random_SNPs] ~ Annual_Tc+Diurnal_range_Tc+Tc_Seasonality+Tc_driest_quarter+Annual_P+P_Seasonality, data = Climatic_data_RDA_pRDA, scale=F)
RsquareAdj(RDA_random)## $r.squared
## [1] 0.2976353
##
## $adj.r.squared
## [1] 0.1060813- random set keeping the same AF
RDA_random_same_AF <- rda(formula = genomic_matrix[list_random_SNPs_same_AF] ~ Annual_Tc+Diurnal_range_Tc+Tc_Seasonality+Tc_driest_quarter+Annual_P+P_Seasonality, data = Climatic_data_RDA_pRDA, scale=F)
RsquareAdj(RDA_random_same_AF)## $r.squared
## [1] 0.2871292
##
## $adj.r.squared
## [1] 0.09270992- Set of overlapping SNPs with common garden genomic data
RDA_CG <- rda(formula = genomic_matrix[CG_common_set] ~ Annual_Tc+Diurnal_range_Tc+Tc_Seasonality+Tc_driest_quarter+Annual_P+P_Seasonality, data = Climatic_data_RDA_pRDA, scale=F)
RsquareAdj(RDA_CG)## $r.squared
## [1] 0.5201462
##
## $adj.r.squared
## [1] 0.389277- Less conservative outlier set (LC)
RDA_outliers_LC <- rda(formula = genomic_matrix[outliers_set_final_overlapping_no_LD_LC_new_var] ~ Annual_Tc+Diurnal_range_Tc+Tc_Seasonality+Tc_driest_quarter+Annual_P+P_Seasonality, data = Climatic_data_RDA_pRDA, scale=F)
RsquareAdj(RDA_outliers_LC)## $r.squared
## [1] 0.5173666
##
## $adj.r.squared
## [1] 0.3857393- More conservative outlier set (MC)
RDA_outliers_MC <- rda(formula = genomic_matrix[outliers_set_final_overlapping_no_LD_new_var] ~ Annual_Tc+Diurnal_range_Tc+Tc_Seasonality+Tc_driest_quarter+Annual_P+P_Seasonality, data = Climatic_data_RDA_pRDA, scale=F)
RsquareAdj(RDA_outliers_MC)## $r.squared
## [1] 0.5947706
##
## $adj.r.squared
## [1] 0.4842534- V2 set random same AF (random_AF_V2)
RDA_random_AF_V2 <- rda(formula = genomic_matrix[list_random_SNPs_same_AF_V2] ~ Annual_Tc+Diurnal_range_Tc+Tc_Seasonality+Tc_driest_quarter+Annual_P+P_Seasonality, data = Climatic_data_RDA_pRDA, scale=F)
RsquareAdj(RDA_random_AF_V2)## $r.squared
## [1] 0.3046185
##
## $adj.r.squared
## [1] 0.114969- All outliers (all_outliers)
RDA_all_outliers <- rda(formula = genomic_matrix[unique_outliers] ~ Annual_Tc+Diurnal_range_Tc+Tc_Seasonality+Tc_driest_quarter+Annual_P+P_Seasonality, data = Climatic_data_RDA_pRDA, scale=F)
RsquareAdj(RDA_all_outliers)## $r.squared
## [1] 0.4189863
##
## $adj.r.squared
## [1] 0.260528We can save the models to use them later to evaluate the models:
list_models <- c("RDA_all","RDA_random","RDA_random_same_AF","RDA_CG","RDA_outliers_LC","RDA_outliers_MC","RDA_random_AF_V2","RDA_all_outliers")
for(x in 1:length(list_models)){
name_model<- list_models[x]
file_save<- get(name_model)
name <- paste0("model_GEA_",name_model)
assign(name,file_save)
save(list = paste0("model_GEA_",name_model),file=paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/data/models/",name,".Rdata"))
}To conduct further analysis, we will examine the different sets to see if the same patterns emerge.
list_outliers_set <- c("list_all_SNPs","list_random_SNPs","list_random_SNPs_same_AF","outliers_set_final_overlapping_no_LD_LC_new_var","outliers_set_final_overlapping_no_LD_new_var","CG_common_set")
list_RDA <- c("RDA_all","RDA_random","RDA_random_same_AF","RDA_outliers_LC","RDA_outliers_MC","RDA_CG")
list_outliers_set_V2 <- c("list_all_SNPs","list_random_SNPs","list_random_SNPs_same_AF","outliers_set_final_overlapping_no_LD_LC_new_var","outliers_set_final_overlapping_no_LD_new_var","CG_common_set","list_random_SNPs_same_AF_V2","unique_outliers")
list_RDA_V2 <- c("RDA_all","RDA_random","RDA_random_same_AF","RDA_outliers_LC","RDA_outliers_MC","RDA_CG","RDA_random_AF_V2","RDA_all_outliers")3.2 Graphical visualization
We can plot how the outliers are distributed along the RDA axes and the climatic variables.
for(x in 1: length(list_RDA_V2)){
RDA <- get(list_RDA_V2[x])
name <- c("all","random_SNPs","RDA_random_same_AF","outliers_Lc","Outliers_Mc","CG_set","random_AF_V2","all_outliers")
#screeplot
plot(RDA$CCA$eig, option="screeplot")
#in hist
screeplot_RDA <- screeplot(RDA,main = paste0("Hist screeplot ",name[x]))
#explained variance along each RDA axis
explain_variance <-RDA$CCA$eig*100/sum(RDA$CCA$eig)
assign(paste0("explain_variance",name[x]),explain_variance)
print(get(paste0("explain_variance",name[x])))
#save
#pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/screeplot/Screeplot_RDA_",name[x],".pdf"));screeplot(RDA,main = paste0("Hist screeplot ",name[x]));dev.off()
}
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 59.299725 24.467288 5.571946 4.712949 3.804574 2.143519As can be seen, the first two RDA axes appear to explain most of the genetic variation of outliers across all sets. We will use these two axes in the further analysis of each set.
We can also plot the loci and climatic variables in the RDA space to visualise their relationships:
for(x in 1:length(list_RDA_V2)){
RDA <- get(list_RDA_V2[x])
name <- c("all","random_SNPs","RDA_random_same_AF","outliers_Lc","Outliers_Mc","CG_set","random_AF_V2","all_outliers")
score_loci <- as.data.frame(scores(RDA, choices=c(1:2), display="species", scaling="none"))
score_climatic_var <- as.data.frame(scores(RDA, choices=c(1:2), display="bp"))
explained_variance_round <- round(RDA$CCA$eig*100/sum(RDA$CCA$eig),digits=1)
#Biplot with SNPs and climatic variables along the two first RDA axis.
biplot_outliers_RDA<- ggplot() +
geom_hline(yintercept=0, linetype="dashed", color = gray(.80), size=0.6) +
geom_vline(xintercept=0, linetype="dashed", color = gray(.80), size=0.6) +
geom_point(data = score_loci, aes(x=RDA1*15, y=RDA2*15,color="locus"), size = 1.4) +
geom_segment(data = score_climatic_var, aes(xend=RDA1, yend=RDA2, x=0, y=0), colour="black", size=0.15, linetype=1, arrow=arrow(length = unit(0.02, "npc"))) +
geom_text(data = score_climatic_var, aes(x=1.1*RDA1, y=1.1*RDA2, label = row.names(score_climatic_var)), size = 2.5) +
xlab(paste0("RDA 1 (",explained_variance_round[1],"%)")) +
ylab(paste0("RDA 2 (",explained_variance_round[2],"%)")) +
ggtitle(paste0("Biplot RDA ",name[x])) +
guides(color=guide_legend(title="Locus type")) +
scale_color_manual(values ="#F9A242FF") +
theme_bw(base_size = 11) +
theme(legend.position="right", panel.grid = element_blank(), strip.text = element_text(size=11),plot.title = element_text(hjust = 0.5,color = "Black",face="italic"))
print(biplot_outliers_RDA)
#save plots
#pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/biplot/biplot_outliers_RDA_",name[x],".pdf"));print(biplot_outliers_RDA);dev.off()
}
Along the first axis, we observe a pattern of continentality/elevation, with colder temperatures, greater temperature seasonality and a higher diurnal temperature range on the left, which is associated with continental or high-altitude regions. On the right, there are outliers associated with higher temperatures, particularly during the driest quarter, as well as increased precipitation seasonality. These are linked to less continental or lower-altitude areas. For the MC set, the pattern is similar but reversed along the first RDA axis.
We can plot the populations in this space to see if their distribution follows the expected climatic pattern. Here, the population positions represent observed values rather than predicted ones, as in other studies.
for(x in 1:length(list_RDA_V2)){
RDA <- get(list_RDA_V2[x])
name <- c("all","random_SNPs","RDA_random_same_AF","outliers_Lc","Outliers_Mc","CG_set","random_AF_V2","all_outliers")
#score along the 2 first RDA axis
score_climatic_var <- as.data.frame(scores(RDA, choices=c(1:2), display="bp"))
#Score_population <- data.frame(RDA$CCA$u[,c(1,2)]) #predicted population scores
observed_scores_populations <- as.data.frame(scores(RDA,choices=c(1,2), display = "sites",scaling="none"))
#explained variance
explained_variance_round <- round(RDA$CCA$eig*100/sum(RDA$CCA$eig),digits=1)
#merge for country info
Score_population_bis <- rownames_to_column(observed_scores_populations,"Population")
score_with_country_info <- merge(Score_population_bis,meta_data_pop_order[,c(1,2)],"Population")
score_with_country_info$Country <- as.factor(score_with_country_info$Country)
group_palette <- c("Bosnia"="orangered3", "France"="gold2","Germany"= "darkorchid3", "Greece"="navyblue", "Italy"="turquoise2", "Norway"="green3", "Slovakia"="blue", "Slovenia"="red", "Spain"="black", "Sweden"="gray", "Switzerland"="orange", "UK"="darkgreen")
##Biplot with populations and climatic variables along the 2 first RDA axis
biplot_populations <- ggplot() +
geom_hline(yintercept = 0, linetype = "dashed", color = gray(0.80), size = 0.6) +
geom_vline(xintercept = 0, linetype = "dashed", color = gray(0.80), size = 0.6) +
geom_point(data = score_with_country_info, aes(x = RDA1 * 3, y = RDA2 * 3, colour = Country), size = 2, alpha = 0.8) +
geom_segment(data = score_climatic_var, aes(xend = RDA1, yend = RDA2, x = 0, y = 0), colour = "black", size = 0.15, linetype = 1, arrow = arrow(length = unit(0.02, "npc"))) +
geom_text(data = score_climatic_var, aes(x=1.1*RDA1, y=1.1*RDA2, label = row.names(score_climatic_var)), size = 2.5)+
xlab(paste0("RDA 1 (",explained_variance_round[1],"%)")) +
ylab(paste0("RDA 2 (",explained_variance_round[2],"%)")) +
ggtitle(paste0("Biplot RDA Populations ",name[x])) +
scale_color_manual(name = "Countries", values = group_palette, labels = levels(score_with_country_info$Country)) +
theme_bw(base_size = 11) +
theme(legend.position="right", panel.grid = element_blank(), strip.text = element_text(size=11),plot.title = element_text(hjust = 0.5,color = "Black",face="italic"))+
labs(color = "Country")
print(biplot_populations)
#save plots
#pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/biplot/biplot_population_RDA_",name[x],".pdf"));print(biplot_populations);dev.off()
}
As before, we observe that the SNP pattern aligns here: more continental or high-altitude populations are located on the left, with the opposite pattern evident for the LC set on the right. Similarly, the pattern along the first RDA axis is reversed for the MC outlier set, though the information remains consistent.
4 Interpolate/extrapolate the GEA relationship to space and time to calculate the genomic composition
We can now calculate the predicted genomic composition for past, present and future climates. Please note that, unlike in the adaptive_index function, climatic values must be scaled prior to using this function. These results will then be used to calculate the GDI index and genomic offset.
list_rda <- c("RDA_all","RDA_random","RDA_random_same_AF","RDA_outliers_LC","RDA_outliers_MC","RDA_CG","RDA_random_AF_V2","RDA_all_outliers")
#list_clim <- c("clim_df_past_scale","clim_df_present_scale","clim_df_future_scale")
list_period <- c("past","present","future")
list_set <- c("all","random","random_same_AF","LC","MC","CG","random_AF_V2","all_outliers")
K=2 #number of rda axes retained
for(i in 1: length(list_rda)){
rda_model <- get(list_rda[i])
set <- list_set[i]
for(x in 1: length(list_period)){
name <- list_period[x]
if(name == "future"){
for(a in 1:length(list_models_clim)){
models_clim <- list_models_clim[a]
clim <- get(paste0("clim_df_",name,"_scale_",models_clim))
predictions <- data.frame(Population=meta_data_pop_order$Population,predict(rda_model,clim, type = "lc", rank = K, scaling = "none"))
assign(paste0("adaptive_values_",set,"_",name,"_",models_clim),predictions)
}
}else{
clim <- get(paste0("clim_df_",name,"_scale"))
#model
predictions <- data.frame(Population=meta_data_pop_order$Population,predict(rda_model,clim, type = "lc", rank = K, scaling = "none"))
assign(paste0("adaptive_values_",set,"_",name),predictions)
}
}
}We compared the predicted values obtained using the Predict and Loading methods.
Although the two methods do not produce identical values, we observe a high degree of correlation between them. This suggests that we should select one method and use it consistently when calculating the different adaptive compositions and subsequently the genomic offset.
5 Genomic discrepancy index (GDI)
The GDI is an index that calculates the difference between the predicted and observed genomic compositions. This enables us to identify populations that differ more from what would be expected based on local adaptation relationships. This suggests that these populations may have a lower degree of local adaptation, or that there may be a non-linear relationship (e.g. conditional neutrality), or that the SNPs used may not fully inform us about local adaptation that could arise due to the absence of local adaptation.
We applied the same method used for genomic offset to calculate this: measuring the distance between two genomic compositions, here, the observed and predicted values at the time of establishment (1901–1950, referred to as ‘past’).
5.1 Calculation
genomic_discrepancy_index <- function(RDA, K, observed_score, predicted_score,meta_data){
# Weights based on axis eigen values
weights <- RDA$CCA$eig/sum(RDA$CCA$eig)
# Weighing the current and future adaptive indices based on the eigen values of the associated axes
Proj_offset_observed <- as.data.frame(do.call(cbind, lapply(1:K, function(x) observed_score[,x]*weights[x])))
predicted_score_df <- predicted_score[,-1]
Proj_offset_predicted <- as.data.frame(do.call(cbind, lapply(1:K, function(x) predicted_score_df[,x]*weights[x])))
#Now we want to calculate the distance between present and future for each RDA axis before doing it for both axis simultaneously
Proj_offset <- list()
for(i in 1:K){
Proj_offset[[i]] <- abs(Proj_offset_observed[[i]] - Proj_offset_predicted[[i]])
names(Proj_offset)[i] <- paste0("RDA", as.character(i))
}
# Predict a global genetic offset, incorporating the K first axes weighted by their eigen values
ras <- Proj_offset[[1]] #we reused the format of the previous distance per RDA axis
ras[!is.na(ras)] <- unlist(lapply(1:nrow(Proj_offset_observed), function(x) dist(rbind(Proj_offset_observed[x,], Proj_offset_predicted[x,]), method = "euclidean"))) #calculation of the euclidean distance on the non Na values of the previous distance -> that why we used the format of the previous distance, to be sure to only select the rows without Nas because they are not deal by euclidean distance,
#the euclidean distance is still calculated on the weighted data (not the previous distance but on the genomic composition weighted)
names(ras) <- "GDI"
Proj_offset_global <- ras
# Return prediction of genetic offset for each RDA axis and a global genetic offset for each population
return(list(Population=meta_data$Population,Proj_offset = Proj_offset, Proj_offset_global = Proj_offset_global, weights = weights[1:K]))
}list_rda <- c("RDA_all","RDA_random","RDA_random_same_AF","RDA_outliers_LC","RDA_outliers_MC","RDA_CG","RDA_random_AF_V2","RDA_all_outliers")
list_predicted_adaptive_values <- c("adaptive_values_all_past","adaptive_values_random_past","adaptive_values_random_same_AF_past","adaptive_values_LC_past","adaptive_values_MC_past","adaptive_values_CG_past","adaptive_values_random_AF_V2_past","adaptive_values_all_outliers_past")
for(x in 1: length(list_rda)){
rda <- get(list_rda[x])
observed_values <- scores(rda, display = "sites",scaling="none")
predicted_values <- data.frame(get(list_predicted_adaptive_values[x]))
set <- list_set[x]
calculated_GDI <- genomic_discrepancy_index(rda,K,observed_values,predicted_values,meta_data_pop_order)
genomic_offset_df<- data.frame(Population=unlist(calculated_GDI$Population),Genomic_offset_random=unlist(calculated_GDI$Proj_offset_global))
colnames(genomic_offset_df) <- c("Population",paste0("values_",set))
assign(paste0("GDI_index_",set),genomic_offset_df)
}We can save the values of GDI.
Compare the GDI with the residuals output by the residual function. The residual function provides a residual for each SNP in a population, so we need to summarise this information.
res <- residuals(RDA_random_same_AF)
res_df <- as.data.frame(res)
#change residual value per SNP to numeric
res_df_numeric <- res_df %>%
dplyr::mutate(across(everything(), as.numeric))
#add population info
res_df_numeric$Population <- row.names(res_df_numeric)
#calculate the mean residual across SNPs for each population
res_df_with_mean <- res_df_numeric %>%
mutate(row_mean = rowMeans(dplyr::select(., where(is.numeric)), na.rm = TRUE))
#save the mean value of residual for each pop
res_df_with_mean_final <- data.frame(Population=row.names(res_df_numeric),residual=res_df_with_mean[,c(102)])
#load GDI value
load("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GDI/RDA/data/GDI_index_random_same_AF.Rdata")
#calculate correlation
cor(res_df_with_mean_final$residual,GDI_index_random_same_AF$values_random_same_AF)## [1] 0.31553385.1.1 Relative GDI
Not used in this study.
genomic_discrepancy_index <- function(RDA, K, observed_score, predicted_score, meta_data) {
# Weights based on axis eigenvalues
weights <- RDA$CCA$eig / sum(RDA$CCA$eig)
# Weighted observed and predicted scores
Proj_offset_observed <- as.data.frame(do.call(cbind, lapply(1:K, function(x) observed_score[,x] * weights[x])))
predicted_score_df <- predicted_score[,-1]
Proj_offset_predicted <- as.data.frame(do.call(cbind, lapply(1:K, function(x) predicted_score_df[,x] * weights[x])))
# Calculate the absolute differences for each axis
Proj_offset <- list()
for (i in 1:K) {
Proj_offset[[i]] <- abs(Proj_offset_observed[[i]] - Proj_offset_predicted[[i]])
names(Proj_offset)[i] <- paste0("RDA", as.character(i))
}
# Calculate the global GDI using Euclidean distance across all axes
Proj_offset_global <- unlist(lapply(1:nrow(Proj_offset_observed), function(x) {
dist(rbind(Proj_offset_observed[x,], Proj_offset_predicted[x,]), method = "euclidean")
}))
names(Proj_offset_global) <- "GDI"
# Calculate the relative index as a Euclidean distance of relative differences
relative_change <- unlist(lapply(1:nrow(Proj_offset_observed), function(x) {
# Extract the observed and predicted scores for this row as numeric vectors
observed_row <- as.numeric(Proj_offset_observed[x, ])
predicted_row <- as.numeric(Proj_offset_predicted[x, ])
# Calculate relative differences for all axes
relative_diff <- (predicted_row - observed_row) / observed_row
relative_diff[is.nan(relative_diff) | is.infinite(relative_diff)] <- 0 # Handle division by zero or NA
# Compute Euclidean distance of relative differences
dist(rbind(rep(0, K), relative_diff), method = "euclidean")
}))
names(relative_change) <- "Relative_Index"
# Return results
return(list(
Population = meta_data$Population,
Proj_offset = Proj_offset,
Proj_offset_global = Proj_offset_global,
Relative_change = relative_change,
weights = weights[1:K]
))
}list_rda <- c("RDA_all", "RDA_random", "RDA_random_same_AF", "RDA_outliers_LC", "RDA_outliers_MC", "RDA_CG")
list_predicted_adaptive_values <- c("adaptive_values_all_past", "adaptive_values_random_past", "adaptive_values_random_same_AF_past", "adaptive_values_LC_past", "adaptive_values_MC_past", "adaptive_values_CG_past")
for (x in 1:length(list_rda)) {
rda <- get(list_rda[x])
observed_values <- scores(rda, display = "sites", scaling = "none")
predicted_values <- data.frame(get(list_predicted_adaptive_values[x]))
set <- list_set[x]
# Calculate GDI and Relative Change using the updated function
calculated_GDI <- genomic_discrepancy_index(rda, K, observed_values, predicted_values, meta_data_pop_order)
# Create a dataframe for GDI and combined relative index values
genomic_offset_df <- data.frame(
Population = unlist(calculated_GDI$Population),
Genomic_offset = unlist(calculated_GDI$Proj_offset_global),
Relative_Index = unlist(calculated_GDI$Relative_change) # Add the combined relative index
)
# Rename columns to be more informative
colnames(genomic_offset_df) <- c("Population", paste0("values_", set), paste0("Relative_Index_", set))
# Assign the dataframe to a new variable with a dynamic name
assign(paste0("GDI_index_bis_", set), genomic_offset_df)
}5.2 Plot the GDI values
We can plot the GDI values for each population. As before, we can use the ‘cut_interval’ function (never use ‘cut_number’!) to rescale the values so that they match the colour scale used for other genomic offsets, such as GF.
shapefile <- st_read("C:/Users/tfrancisco/Documents/Thesis/Data/Species/Taxus_baccata/shapefiles/Taxus_baccata_plg_clip.shp")
list_period <- c("present","future")
list_set <- c("all","random","random_same_AF","LC","MC","CG","random_AF_V2","all_outliers")
# Loop to find the global min and max values across all data frames to cut the values and be able to compare the values across model.
#for(set in list_set){
# GDI_df <- get(paste0("GDI_index_",set))
# global_min <- min(global_min, min(GDI_df[,2], na.rm = TRUE))
# global_max <- max(global_max, max(GDI_df[,2], na.rm = TRUE))
#}
for(x in 1: length(list_set)){
set <- list_set[x]
GDI_df <- get(paste0("GDI_index_",set))
#first, we need to add the coordinates
GDI_coord <- merge(GDI_df,meta_data_pop_order[,c(2,4,5)],"Population")
#transform longitude and latitude to numeric variables
GDI_coord <- GDI_coord %>% mutate(Longitude=as.numeric(Longitude),Latitude=as.numeric(Latitude))
colors <- c( "darkgreen", "#FDF7B3","#FC4E2A","#BD0026","darkorchid4")
colors <- c("#005000","#85B22C","#FFFFBF","#F76D5E","#E51932")
#background map
admin <- ne_countries(scale = "medium", returnclass = "sf")
plot <- ggplot(data = GDI_coord) +
geom_sf(data = admin, fill = gray(0.92), size = 0) +
geom_sf(data = shapefile, fill = "lightgrey", color = "black", size = 0.2) + # Adding shapefile
geom_point(aes(x = Longitude, y = Latitude, fill = cut_interval(GDI_coord[,2], n = 5)), shape = 21,size=3, color = "black") +
#geom_point(aes(x = Longitude, y = Latitude, fill = cut(GDI_coord[,2], breaks=seq(global_min, global_max, length.out = 6), include.lowest = T)), shape = 21,size=3, color = "black") + #if we want to class with the same values the 4 models to compare color
scale_fill_manual(
values = colors,
labels = c("0","","","","1"),
drop = FALSE,na.translate = FALSE)+ # Ensure all levels are shown in the legend
geom_sf(data = admin, fill = NA, size = 0.1) +
coord_sf(xlim = c(-10, 30), ylim = c(35, 62), expand = FALSE) +
xlab("Longitude") + ylab("Latitude") +
guides(fill = guide_legend(title = NULL)) +
ggtitle(paste0("Adaptive genomic disruption across populations ",set)) +
annotation_scale(location = "br", width_hint = 0.3) + # Scale bar at bottom left
annotation_north_arrow(location = "tl", which_north = "true",
pad_x = unit(-0.2, "cm"), pad_y = unit(0.08, "cm"),
style = north_arrow_fancy_orienteering) + # Compass rose
theme_bw(base_size = 11) +
theme(legend.position="right", panel.grid = element_blank(), strip.text = element_text(size=11),plot.title = element_text(hjust = 0.5,color = "Black",face="italic"))
print(plot)
#save
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GDI/RDA/figures/GDI_across_populations ",set,"_newcolors.pdf"));print(plot);dev.off()
}
Save the GDI using the set of SNPs containing all the outliers for Chapter 2.
write.csv(GDI_index_random_same_AF,file="C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GDI/RDA/predictions/GDI_index_random_same_AF.csv")list_set <- c("all","random","random_same_AF","LC","MC","CG")
# Combine all GDI values into one dataframe
combined_data <- data.frame()
for (set in list_set) {
GDI_df <- get(paste0("GDI_index_", set))
colnames(GDI_df)[2] <- "GDI" # Ensure second column is named GDI
GDI_df$Set <- set
combined_data <- rbind(combined_data, GDI_df[, c("GDI", "Set")])
}
# Density plot
density_plot <- ggplot(combined_data, aes(x = GDI, fill = Set)) +
geom_density(alpha = 0.7, color = "black") +
scale_fill_manual(values = c("darkred", "#E69F00", "#56B4E9", "darkgreen", "#009E73", "#CC79A7")) +
#geom_vline(xintercept = 1, color = "red", linetype = "dashed", size = 1) + # Optional threshold
facet_wrap(~ Set, scales = "free_y", ncol = 2) +
labs(
title = "Density of Genomic Discrepancy Index (GDI) across datasets",
x = "Genomic Discrepancy Index (GDI)",
y = "Density"
) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5, face = "italic"),
strip.text = element_text(size = 10),
legend.position = "none"
)
print(density_plot)
# Save to PDF
pdf("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GDI/RDA/figures/GDI_density_faceted.pdf")
print(density_plot)
dev.off()## png
## 2list_set <- c("all","random","random_same_AF","LC","MC","CG")
# Combine all GDI values into one dataframe
combined_data <- data.frame()
for (set in list_set) {
GDI_df <- get(paste0("GDI_index_", set))
colnames(GDI_df)[2] <- "GDI" # Ensure second column is named GDI
GDI_df$Set <- set
combined_data <- rbind(combined_data, GDI_df[, c("GDI", "Set")])
}
# Plot: histogram of GDI values, faceted by set
plot <- ggplot(combined_data, aes(x = GDI)) +
geom_histogram(aes(fill = Set), binwidth = 0.005, color = "black", alpha = 0.8) +
scale_fill_manual(values = c("darkred", "#E69F00", "#56B4E9", "darkgreen", "#009E73", "#CC79A7")) +
#geom_vline(xintercept = 1, color = "red", linetype = "dashed", size = 1) + # Adjust or remove if not relevant
facet_wrap(~ Set, scales = "free_y", ncol = 2) +
labs(
title = "Distribution of Genomic Discrepancy Index (GDI) across datasets",
x = "Genomic Discrepancy Index (GDI)",
y = "Number of Populations"
) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5, face = "italic"),
strip.text = element_text(size = 10),
legend.position = "none"
)
print(plot)
# #Save
pdf("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GDI/RDA/figures/GDI_distribution_histogram_faceted.pdf")
print(plot)
dev.off()## png
## 2# Prepare an empty data frame to collect means
mean_GDI_df <- data.frame(Set = character(), Mean_GDI = numeric(), stringsAsFactors = FALSE)
for (set in list_set) {
GDI_df <- get(paste0("GDI_index_", set))
colnames(GDI_df)[2] <- "GDI"
mean_val <- mean(GDI_df$GDI, na.rm = TRUE)
mean_GDI_df <- rbind(mean_GDI_df, data.frame(Set = set, Mean_GDI = mean_val))
}
# Plot
barplot <- ggplot(mean_GDI_df, aes(x = Set, y = Mean_GDI, fill = Set)) +
geom_bar(stat = "identity", color = "black") +
ylab("Mean Genomic Discrepancy Index") +
xlab("Set") +
ggtitle("Mean GDI per Dataset") +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(hjust = 0.5, face = "italic"),
legend.position = "none"
)
print(barplot)
# Save plot
ggsave(
filename = "C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GDI/RDA/figures/Mean_GDI_barplot_across_sets.pdf",
plot = barplot, width = 8, height = 6
)Example LC: Overall, an interesting pattern emerges, with higher GDI values observed in populations that appear to be isolated from others, such as those in Greece. However, this pattern is inconsistent, as high GDI values are not observed in isolated northern populations (e.g. Sweden and Norway). We can also compare GDI values across models and with respect to elevation.
#df GDI and elevation
df_GDI_cor <- data.frame(GDI_all=GDI_index_all$values_all,GDI_random=GDI_index_random$values_random,GDI_random_same_AF=GDI_index_random_same_AF$values_random_same_AF,GDI_LC=GDI_index_LC$values_LC,GDI_MC=GDI_index_MC$values_MC,GDI_CG=GDI_index_CG$values_CG,elevation=as.numeric(meta_data_pop_order$Elevation.DEM_90m.),Latitude=as.numeric(meta_data_pop_order$Latitude))
#correlation
correlation_GDI <- cor(df_GDI_cor)
corrplot(correlation_GDI, method = "number", addrect = 2, col = c("darkorange","darkred"), type = "lower", tl.col = "black", tl.cex = 0.6, number.cex = 0.6)
pdf("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/GDI/RDA/figures/comparison/correlation_GDI_values_RDA_taxus.pdf");corrplot(correlation_GDI, method = "number", addrect = 2, col = c("darkorange","darkred"), type = "lower", tl.col = "black", tl.cex = 0.6, number.cex = 0.6);dev.off()## png
## 26 Genomic offset
In our study, genomic offset is defined as the Euclidean distance between future and present genomic compositions. This could also represent the distance between the same population in different common gardens, or between populations in different locations and time periods. To calculate this distance, we consider genomic composition to be estimated using the first two RDA axes, and weight the values for each population based on the genetic variation explained by these axes.
We adapted the function from Capblancq and Forester (2021) so that our analysis focuses solely on calculating the genomic offset for populations without interpolating or extrapolating across space.
6.1 Genomic offset with predicted genomic composition
6.1.1 Model
Due to the significant computational time required by the function used in Capblancq and Forester (2021) to calculate the genomic offset (GO) for all raster points, we modified the function to calculate the GO exclusively for our sampled points.
genomic_offset_pop <- function(RDA, K, Past_score, Future_score,meta_data){
# Weights based on axis eigen values
weights <- RDA$CCA$eig/sum(RDA$CCA$eig)
# Weighing the current and future adaptive indices based on the eigen values of the associated axes
Past_score_df <- Past_score[,-c(1)]
Proj_offset_past <- as.data.frame(do.call(cbind, lapply(1:K, function(x) Past_score_df[,x]*weights[x])))
Future_score_df <- Future_score[,-c(1)]
Proj_offset_fut <- as.data.frame(do.call(cbind, lapply(1:K, function(x) Future_score_df[,x]*weights[x])))
#Now we want to calculate the distance between present and future for each RDA axis before doing it for both axis simultaneously
Proj_offset <- list()
for(i in 1:K){
Proj_offset[[i]] <- abs(Proj_offset_past[[i]] - Proj_offset_fut[[i]])
names(Proj_offset)[i] <- paste0("RDA", as.character(i))
}
# Predict a global genetic offset, incorporating the K first axes weighted by their eigen values
ras <- Proj_offset[[1]] #we reused the format of the previous distance per RDA axis
ras[!is.na(ras)] <- unlist(lapply(1:nrow(Proj_offset_past), function(x) dist(rbind(Proj_offset_past[x,], Proj_offset_fut[x,]), method = "euclidean"))) #calculation of the euclidean distance on the non Na values of the previous distance -> that why we used the format of the previous distance, to be sure to only select the rows without Nas because they are not deal by euclidean distance,
#the euclidean distance is still calculated on the weighted data (not the previous distance but on the genomic composition weighted)
names(ras) <- "Global_offset"
Proj_offset_global <- ras
# Return prediction of genetic offset for each RDA axis and a global genetic offset for each population
return(list(Population=meta_data$Population,Proj_offset = Proj_offset, Proj_offset_global = Proj_offset_global, weights = weights[1:K]))
}We can run the GO function for all the set of genetic markers for present and future climate:
list_RDA <- c("RDA_all","RDA_random", "RDA_random_same_AF", "RDA_outliers_LC", "RDA_outliers_MC", "RDA_CG","RDA_random_AF_V2","RDA_all_outliers")
list_past_score <- c("adaptive_values_all_past","adaptive_values_random_past", "adaptive_values_random_same_AF_past", "adaptive_values_LC_past", "adaptive_values_MC_past", "adaptive_values_CG_past","adaptive_values_random_AF_V2_past","adaptive_values_all_outliers_past")
list_present_score <- c("adaptive_values_all_present","adaptive_values_random_present", "adaptive_values_random_same_AF_present", "adaptive_values_LC_present", "adaptive_values_MC_present", "adaptive_values_CG_present","adaptive_values_random_AF_V2_present","adaptive_values_all_outliers_present")
#list_future_score <- c("adaptive_values_all_future","adaptive_values_random_future", "adaptive_values_random_same_AF_future", "adaptive_values_LC_future", "adaptive_values_MC_future", "adaptive_values_CG_future")
list_name_present <- c("all_present","random_present", "random_same_AF_present", "LC_present", "MC_present", "CG_present","random_AF_V2_present","all_outliers_present")
list_name_future <- c("all_future","random_future", "random_same_AF_future", "LC_future", "MC_future", "CG_future","random_AF_V2_future","all_outliers_future")
# Loop through each RDA model
for (i in 1:length(list_RDA)) {
# Get the RDA model, past score, present score, and future score
RDA <- get(list_RDA[i])
adaptive_past_score <- get(list_past_score[i])
present_score <- get(list_present_score[i])
# Calculate genomic offset for the present scenario
name_present <- list_name_present[i]
Run_genomic_offset_pop_present <- genomic_offset_pop(RDA, 2, adaptive_past_score, present_score, meta_data_pop_order)
# Create a data frame for the present scenario genomic offset
genomic_offset_df_present <- data.frame(Population = unlist(Run_genomic_offset_pop_present$Population),
Genomic_offset_random = unlist(Run_genomic_offset_pop_present$Proj_offset_global))
# Order the data frame by population and rename columns
Genomic_offset_Taxus_Adapcon_Gentree_RDA_present <- genomic_offset_df_present[order(genomic_offset_df_present$Population), ]
colnames(Genomic_offset_Taxus_Adapcon_Gentree_RDA_present) <- c("Population", paste0("GO_", name_present))
# Assign the present result to a new variable
assign(paste0("Genomic_offset_Taxus_Adapcon_Gentree_RDA_", name_present), Genomic_offset_Taxus_Adapcon_Gentree_RDA_present)
# Calculate genomic offset for the future scenario
for (j in 1:length(list_models_clim)) { # Loop through 5 future models
name_clim <- list_models_clim[j]
name_future <- paste0(list_name_future[i],"_", name_clim)
future_score <- get(paste0("adaptive_values_",list_name_future[i],"_", name_clim))
Run_genomic_offset_pop_future <- genomic_offset_pop(RDA, 2, adaptive_past_score, future_score, meta_data_pop_order)
# Create a data frame for the future scenario genomic offset
genomic_offset_df_future <- data.frame(Population = unlist(Run_genomic_offset_pop_future$Population),
Genomic_offset_random = unlist(Run_genomic_offset_pop_future$Proj_offset_global))
# Order the data frame by population and rename columns
Genomic_offset_Taxus_Adapcon_Gentree_RDA_future <- genomic_offset_df_future[order(genomic_offset_df_future$Population), ]
colnames(Genomic_offset_Taxus_Adapcon_Gentree_RDA_future) <- c("Population", paste0("GO_", name_future))
# Assign the future result to a new variable
assign(paste0("Genomic_offset_Taxus_Adapcon_Gentree_RDA_",name_future), Genomic_offset_Taxus_Adapcon_Gentree_RDA_future)
}
}We can save the genomic offset values for downstream analysis.
save present predictions for chapter 2 for the random set of SNPs
save(GO_T_Adapcon_Gentree_RDA_random_same_AF_present,file="C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/predictions/GO_T_Adapcon_Gentree_RDA_random_same_AF_present.Rdata")
write.csv(GO_T_Adapcon_Gentree_RDA_random_same_AF_present,file="C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/predictions/GO_T_Adapcon_Gentree_RDA_random_same_AF_present.csv")6.1.2 Comparison GO predictions across climatic models
list_set <- c("all","random","random_same_AF","LC","MC","CG","random_AF_V2","all_outliers")
for(x in 1:length(list_set)){
set <- list_set[x]
df_tot_GO_predictions_clim_models_future <- data.frame(Population=GO_T_Adapcon_Gentree_RDA_random_future_GFDL_ESM4[,1],GFDL_ESM4=get(paste0("GO_T_Adapcon_Gentree_RDA_",set,"_future_GFDL_ESM4"))[,2],IPSL_CM6A_LR=get(paste0("GO_T_Adapcon_Gentree_RDA_",set,"_future_IPSL_CM6A_LR"))[,2],MPI_ESM1_2_HR=get(paste0("GO_T_Adapcon_Gentree_RDA_",set,"_future_MPI_ESM1_2_HR"))[,2],MRI_ESM2_0=get(paste0("GO_T_Adapcon_Gentree_RDA_",set,"_future_MRI_ESM2_0"))[,2],UKESM1_0_LL=get(paste0("GO_T_Adapcon_Gentree_RDA_",set,"_future_UKESM1_0_LL"))[,2])#elevation=as.numeric(meta_data_pop_order$Elevation.DEM_90m.
correlation_future <- cor(df_tot_GO_predictions_clim_models_future[,-1])
corrplot(correlation_future, method = "number", addrect = 2, col = c("darkorange","darkred"), type = "lower", tl.col = "black", tl.cex = 0.6, number.cex = 0.6)
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/comparison/GCMs/correlation_future_predictions_GCMs_",set,".pdf"));corrplot(correlation_future, method = "number", addrect = 2, col = c("darkorange","darkred"), type = "lower", tl.col = "black", tl.cex = 0.6, number.cex = 0.6);dev.off()
}
list_models_clim<- c("GFDL_ESM4","IPSL_CM6A_LR","MPI_ESM1_2_HR","MRI_ESM2_0","UKESM1_0_LL")
for(x in 1:length(list_set)){
set <- list_set[x]
all_GO <- data.frame()
for(i in 1:length(list_models_clim)){
model <- list_models_clim[i]
GO_data <- get(paste0("Genomic_offset_Taxus_Adapcon_Gentree_RDA_", set, "_future_", model))
colnames(GO_data) <- c("Population","GO")
all_GO <- rbind(all_GO,GO_data)
}
mean_GO <- aggregate(GO ~ Population, data = all_GO, FUN = mean, na.rm = TRUE)
colnames(mean_GO) <- c("Population", "Mean_GO") # Rename the columns for clarity
assign(paste0("GO_T_Adapcon_Gentree_RDA_", set, "_future_mean"),mean_GO)
#save
file_save<- get(paste0("GO_T_Adapcon_Gentree_RDA_", set, "_future_mean"))
name <- paste0("GO_T_Adapcon_Gentree_RDA_", set, "_future_mean")
save(list = name,file=paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/data/",name,".Rdata"))
}list_models_clim<- c("mean","GFDL_ESM4","IPSL_CM6A_LR","MPI_ESM1_2_HR","MRI_ESM2_0","UKESM1_0_LL")
# List of SNP sets
list_sets <- c("all","random","random_same_AF","LC","MC","CG","random_AF_V2","all_outliers")
# Initialize an empty list to store data for all SNP sets
plot_data <- list()
# Load the RColorBrewer package for dynamic color palettes
library(RColorBrewer)
color_palette <- c("#FF69B4","red","blue","orange","black","green")
# Initialize a vector to keep track of top populations across all GCMs
meta_data_pop$pop_abbrev <- c("VIS","BAU","SDF","PAT","CHO","VOU","OLY","ANC","PDT","UMB","VAL","BRA","HAR","GOR","BOM","UNS","BUJA","BUJG","CAR","RAS","SUE","SAV","BLI","HHR","JUR","WAL","CED","RWM","YEW")
meta_data_f <- meta_data_pop[order(meta_data_pop$Population),]
# Loop over each SNP set
for (x in 1:length(list_sets)) {
set <- list_sets[x]
# Initialize an empty data frame for all models in the current SNP set
ranks_df <- data.frame()
all_top_populations <- character(0)
# Loop over each climatic model
for (i in 1:length(list_models_clim)) {
model <- list_models_clim[i]
# Get the genomic offset data for the current SNP set and model
GO_data <- get(paste0("GO_T_Adapcon_Gentree_RDA_", set, "_future_", model))
colnames(GO_data) <- c("Population", "GO")
GO_data_V2 <- merge(GO_data,meta_data_f[,c(1,2,7)],"Population")
# Rank the populations by genomic offset (highest GO gets rank 1)
GO_data_V2$Rank <- rank(-GO_data_V2$GO, ties.method = "first")
# Add columns for the climatic model and SNP set
GO_data_V2$Model <- model
GO_data_V2$Set <- set
# Append to ranks dataframe
ranks_df <- rbind(ranks_df, GO_data_V2)
# Find the top 1 population for the current model (highest GO = rank 1)
top_1_population <- GO_data_V2 %>% top_n(1, GO) %>% pull(Population)
# Add this top 1 population to the list of all top populations (avoid duplicates)
all_top_populations <- unique(c(all_top_populations, top_1_population))
}
# Convert 'Model' and 'Set' to factors for plotting
ranks_df$Model <- factor(ranks_df$Model, levels = list_models_clim)
ranks_df$Set <- factor(ranks_df$Set, levels = list_sets)
# Adjust the number of colors if there are more top populations than the current color palette
num_top_populations <- length(all_top_populations)
if (num_top_populations > length(color_palette)) {
color_palette <- colorRampPalette(brewer.pal(9, "Set1"))(num_top_populations) # Use a different palette if needed
}
# Assign a unique color to each of the top populations (those ranked 1 in any model)
population_colors <- setNames(color_palette[1:num_top_populations], all_top_populations)
# Create a new column to store color assignments (highlight only rank 1 populations for each GCM)
ranks_df$Highlight <- ifelse(ranks_df$Population %in% all_top_populations, ranks_df$Population, "Other")
# Map colors for top populations and grey for the rest
ranks_df$Color <- ifelse(ranks_df$Highlight == "Other", "grey", population_colors[ranks_df$Highlight])
# Store the data frame for this SNP set in the plot_data list
plot_data[[set]] <- ranks_df
# Modify the plot generation code
p <- ggplot(ranks_df, aes(x = Model, y = Rank, group = Population)) +
geom_line(aes(color = Highlight), alpha = 0.7, size = 0.8) + # Color only top populations
geom_point(aes(color = Highlight), size = 2) + # Color only top populations
geom_text(
data = ranks_df %>% filter(Model == list_models_clim[1]), # Label only in first column
aes(label = meta_data_f$pop_abbrev[match(Population, meta_data_f$Population)]),
hjust = 1, # Align to the left
vjust = 0,
size = 3.8, # Increase text size
nudge_x = -0.25,
fontface= "italic"
# Move text further left
) +
scale_y_reverse() + # Reverse rank (1 = top)
scale_color_manual(values = c(population_colors, "Other" = "grey")) + # Keep top populations colored, others grey
labs(
title = paste("RDA Rank Variability for SNP Set:", set),
x = "Climatic Models", y = "Population Rank"
) +
theme_minimal() +
theme(
legend.position = "none",
panel.grid.minor = element_blank(),
strip.text = element_text(size = 10, face = "bold"),
legend.title = element_text(size = 10),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 12),
legend.text = element_text(size = 12)
)
print(p) # Display plot for this SNP set
#save
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/comparison/GCMs/GO_rank_comparison_",set,".pdf"));print(p);dev.off()
}
We are going to use the mean of the GO predictions from the five GCMs as the future GO.
6.1.3 Comparison GO predictions across SNPs sets
6.1.4 Graphical visualization
We can plot the genomic offset values for each population. As before, we can use the ‘cut_interval’ function to rescale the values so that they align with the colour scale used for other genomic offsets, such as GF.
list_period <- c("present","future")
list_set <- c("all","random","random_same_AF","LC","MC","CG","random_AF_V2","all_outliers")
for(i in 1:length(list_period)){
period <- list_period[i]
for(x in 1: length(list_set)){
set <- list_set[x]
if(period == "future"){
for (j in 1:length(list_models_clim)) {
model_clim <- list_models_clim[j]
name_file <- paste0(set,"_",period,"_",model_clim)
GO_df <- get(paste0("GO_T_Adapcon_Gentree_RDA_",name_file))
#first, we need to add the coordinates
Genomic_offset_coord <- merge(GO_df,meta_data_pop_order[,c(2,4,5)],"Population")
#transform longitude and latitude to numeric variables
Genomic_offset_coord <- Genomic_offset_coord %>% mutate(Longitude=as.numeric(Longitude),Latitude=as.numeric(Latitude))
colors <- c("#005000","#85B22C","#FFFFBF","#F76D5E","#E51932")
#background map
admin <- ne_countries(scale = "medium", returnclass = "sf")
plot <- ggplot(data = Genomic_offset_coord) +
geom_sf(data = admin, fill = gray(0.92), size = 0) +
geom_point(aes(x = Longitude, y = Latitude, fill = cut_interval(Genomic_offset_coord[,2], n = 5)), shape = 21,size=3, color = "black") +
scale_fill_manual(
values = colors,
labels = c("low values","","","","high values"),
drop = FALSE,na.translate = FALSE)+ # Ensure all levels are shown in the legend
geom_sf(data = admin, fill = NA, size = 0.1) +
coord_sf(xlim = c(-10, 30), ylim = c(35, 62), expand = FALSE) +
xlab("Longitude") + ylab("Latitude") +
guides(fill = guide_legend(title = "Genomic offset")) +
ggtitle(paste0("Genomic offset across populations ",set," ",period,"_",model_clim)) +
theme_bw(base_size = 11) +
theme(legend.position="right", panel.grid = element_blank(), strip.text = element_text(size=11),plot.title = element_text(hjust = 0.5,color = "Black",face="italic"))
print(plot)
#save
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/GO/GCMs_separate_and_mean/future/GO_RDA_standard_across_populations_",set,"_",period,"_",model_clim,".pdf"));print(plot);dev.off()
}
} else {
name_file <- paste0(set,"_",period)
GO_df <- get(paste0("GO_T_Adapcon_Gentree_RDA_",name_file))
#first, we need to add the coordinates
Genomic_offset_coord <- merge(GO_df,meta_data_pop_order[,c(2,4,5)],"Population")
#transform longitude and latitude to numeric variables
Genomic_offset_coord <- Genomic_offset_coord %>% mutate(Longitude=as.numeric(Longitude),Latitude=as.numeric(Latitude))
colors <- c("#005000","#85B22C","#FFFFBF","#F76D5E","#E51932")
#background map
admin <- ne_countries(scale = "medium", returnclass = "sf")
plot <- ggplot(data = Genomic_offset_coord) +
geom_sf(data = admin, fill = gray(0.92), size = 0) +
geom_point(aes(x = Longitude, y = Latitude, fill = cut_interval(Genomic_offset_coord[,2], n = 5)), shape = 21,size=3, color = "black") +
scale_fill_manual(
values = colors,
labels = c("low values","","","","high values"),
drop = FALSE,na.translate = FALSE)+ # Ensure all levels are shown in the legend
geom_sf(data = admin, fill = NA, size = 0.1) +
coord_sf(xlim = c(-10, 30), ylim = c(35, 62), expand = FALSE) +
xlab("Longitude") + ylab("Latitude") +
guides(fill = guide_legend(title = "Genomic offset")) +
ggtitle(paste0("Genomic offset across populations ",set," ",period)) +
theme_bw(base_size = 11) +
theme(legend.position="right", panel.grid = element_blank(), strip.text = element_text(size=11),plot.title = element_text(hjust = 0.5,color = "Black",face="italic"))
print(plot)
#save
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/GO/GCMs_separate_and_mean/present/GO_RDA_standard_across_populations_",set,"_",period,".pdf"));print(plot);dev.off()
}
}
}
6.1.5 Correlation GO values across models
#df
df_tot_GO_predictions_present <- data.frame(Population=GO_T_Adapcon_Gentree_RDA_random_present[,1],GO_all_pres=GO_T_Adapcon_Gentree_RDA_all_present[,2],GO_rand_pres=GO_T_Adapcon_Gentree_RDA_random_present[,2],GO_rand_same_AF_pres=GO_T_Adapcon_Gentree_RDA_random_same_AF_present[,2],GO_LC_pres=GO_T_Adapcon_Gentree_RDA_LC_present[,2],GO_MC_pres=GO_T_Adapcon_Gentree_RDA_MC_present[,2],GO_CG_pres=GO_T_Adapcon_Gentree_RDA_CG_present[,2],elevation=as.numeric(meta_data_pop_order$Elevation.DEM_90m.))
correlation_present <- cor(df_tot_GO_predictions_present[,-1])
corrplot(correlation_present, method = "number", addrect = 2, col = c("darkorange","darkred"), type = "lower", tl.col = "black", tl.cex = 0.6, number.cex = 0.6)
df_tot_GO_predictions_future <- data.frame(Population=GO_T_Adapcon_Gentree_RDA_random_present[,1],GO_all_Fut=GO_T_Adapcon_Gentree_RDA_all_future_mean[,2],GO_rand_Fut=GO_T_Adapcon_Gentree_RDA_random_future_mean[,2],GO_rand_same_AF_Fut=GO_T_Adapcon_Gentree_RDA_random_same_AF_future_mean[,2],GO_LC_Fut=GO_T_Adapcon_Gentree_RDA_LC_future_mean[,2],GO_MC_Fut=GO_T_Adapcon_Gentree_RDA_MC_future_mean[,2],GO_CG_Fut=GO_T_Adapcon_Gentree_RDA_CG_future_mean[,2],GO_random_AF_V2_Fut=GO_T_Adapcon_Gentree_RDA_random_AF_V2_future_mean[,2],GO_all_outliers_Fut=GO_T_Adapcon_Gentree_RDA_all_outliers_future_mean[,2],elevation=as.numeric(meta_data_pop_order$Elevation.DEM_90m.))
correlation_future <- cor(df_tot_GO_predictions_future[,-1])
corrplot(correlation_future, method = "number", addrect = 2, col = c("darkorange","darkred"), type = "lower", tl.col = "black", tl.cex = 0.6, number.cex = 0.6)
pdf("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/comparison/correlation_GO_present_values_RDA_T_set_period.pdf");corrplot(correlation_present, method = "number", addrect = 2, col = c("darkorange","darkred"), type = "lower", tl.col = "black", tl.cex = 0.6, number.cex = 0.6);dev.off()## png
## 2pdf("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/comparison/correlation_GO_future_values_RDA_T_set_period.pdf");corrplot(correlation_future, method = "number", addrect = 2, col = c("darkorange","darkred"), type = "lower", tl.col = "black", tl.cex = 0.6, number.cex = 0.6);dev.off()## png
## 26.1.6 Temporal trend GO
We can also examine the GO trend over time to see if it is consistent.
list_set <- c("all","random","random_same_AF","LC","MC","CG")
for(i in 1:length(list_set)){
set <- list_set[i]
period_1 <- c("present")
period_2 <- c("future_mean")
df_pres <- get(paste0("GO_T_Adapcon_Gentree_RDA_",set,"_",period_1))
df_fut <- get(paste0("GO_T_Adapcon_Gentree_RDA_",set,"_",period_2))
#df
df_f_pres <- data.frame(Population = df_pres[,1], Values = df_pres[,2])
df_f_fut <- data.frame(Population = df_fut[,1], Values = df_fut[,2])
#add period
df_f_pres$Period <- "Present"
df_f_fut$Period <- "Future"
#combine them
df_GO <- rbind(df_f_fut,df_f_pres)
df_GO$Period <- factor(df_GO$Period, levels = c("Present", "Future"))
#plot
trend_plot <- ggplot(df_GO, aes(x = Period, y = Values, group = Population)) +
geom_line(aes(color = Population)) + # Line plot, connecting the points for each population
geom_point(aes(color = Population)) + # Adding points for each value
theme_minimal() +
labs(title = paste0("Genomic offset trend across present and future ",set),
x = "Period",
y = "Genomic offset prediction") +
theme(
legend.position = "none",
strip.text = element_text(size = 11),
plot.title = element_text(hjust = 0.5, color = "Black", face = "italic")
)
plot(trend_plot)
#save
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/comparison/GO_pres_future/Scatterplot_comparison_trend_GO_present_future_",set,"_Taxus.pdf"));print(trend_plot);dev.off()
}
While the values may be similar, we should consider that the ranking of populations could differ. We will investigate this further using graphical visualisation of the ranks.
list_period <- c("present","future_mean")
list_set <- c("all","random","random_same_AF","LC","MC","CG")
for(i in 1:length(list_period)){
period <- list_period[i]
# Loop through each set and compare it with all subsequent sets
for(x in 1:(length(list_set)-1)){
for(y in (x+1):length(list_set)){
set1 <- list_set[x]
set2 <- list_set[y]
# Merge the two datasets
GO_RDA_merge <- merge(get(paste0("GO_T_Adapcon_Gentree_RDA_",set1,"_",period)),
get(paste0("GO_T_Adapcon_Gentree_RDA_",set2,"_",period)),
by = "Population")
# Create a dataframe for comparison
GO_RDA_set_df <- data.frame(Population = GO_RDA_merge$Population,
GO_RDA_1 = GO_RDA_merge[,2],
GO_RDA_2 = GO_RDA_merge[,3])
# Rank the values
GO_RDA_set_df$rank_1 <- rank(GO_RDA_set_df$GO_RDA_1)
GO_RDA_set_df$rank_2 <- rank(GO_RDA_set_df$GO_RDA_2)
# Merge with metadata to include country information
GO_RDA_set_df_meta <- merge(GO_RDA_set_df, meta_data_pop_order[,c(1,2)], by = "Population")
GO_RDA_set_df_meta$Country <- as.factor(GO_RDA_set_df_meta$Country)
# Create the scatter plot
Scatterplot <- ggplot(GO_RDA_set_df_meta, aes(x = rank_1, y = rank_2)) +
geom_point(aes(color = Country), size = 3) +
scale_colour_manual(name = "Countries",
values = c("orangered3","gold2","darkorchid3","navyblue","turquoise2","green3","blue","red","black","gray","orange","darkgreen")) +
geom_abline(intercept = 0, slope = 1, color = "gray60") +
ggtitle(paste0("Comparison GO rank of populations RDA ", set1, "/", set2, " ", period)) +
theme_bw()
# Print the plot
print(Scatterplot)
# Save the plot as a PDF
pdf(paste0("C:/Users/tfrancisco/Documents/Thesis/Results/species/taxus/Genomic_offset/RDA/figures/comparison/Scatterplot_comparison_rank_pop_",set1,"_",set2,"_",period,"_Taxus.pdf"));print(Scatterplot);dev.off()
}
}
}
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.3.2 (2023-10-31 ucrt)
## os Windows 11 x64 (build 26100)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate French_France.utf8
## ctype French_France.utf8
## tz Europe/Paris
## date 2025-09-22
## pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-8 2024-09-12 [1] CRAN (R 4.3.3)
## backports 1.4.1 2021-12-13 [1] CRAN (R 4.3.1)
## base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.3.1)
## broom 1.0.8 2025-03-28 [1] CRAN (R 4.3.3)
## bslib 0.6.1 2023-11-28 [1] CRAN (R 4.3.2)
## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.1)
## car 3.1-2 2023-03-30 [1] CRAN (R 4.3.2)
## carData 3.0-5 2022-01-06 [1] CRAN (R 4.3.2)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.3.2)
## class 7.3-22 2023-05-03 [2] CRAN (R 4.3.2)
## classInt 0.4-10 2023-09-05 [1] CRAN (R 4.3.2)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.1)
## cluster 2.1.6 2023-12-01 [2] CRAN (R 4.3.2)
## codetools 0.2-19 2023-02-01 [2] CRAN (R 4.3.2)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.1)
## corrplot * 0.92 2021-11-18 [1] CRAN (R 4.3.2)
## curl 5.2.0 2023-12-08 [1] CRAN (R 4.3.2)
## data.table 1.15.0 2024-01-30 [1] CRAN (R 4.3.2)
## DBI 1.2.2 2024-02-16 [1] CRAN (R 4.3.2)
## DEoptimR 1.1-3 2023-10-07 [1] CRAN (R 4.3.1)
## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.3.3)
## digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.1)
## dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.2)
## e1071 1.7-14 2023-12-06 [1] CRAN (R 4.3.2)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.1)
## evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.2)
## farver 2.1.2 2024-05-13 [1] CRAN (R 4.3.3)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.1)
## fit.models * 0.64 2020-08-02 [1] CRAN (R 4.3.2)
## fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.1)
## generics 0.1.4 2025-05-09 [1] CRAN (R 4.3.2)
## ggplot2 * 3.5.2 2025-04-09 [1] CRAN (R 4.3.2)
## ggspatial * 1.1.9 2023-08-17 [1] CRAN (R 4.3.3)
## glue 1.7.0 2024-01-09 [1] CRAN (R 4.3.2)
## gtable 0.3.6 2024-10-25 [1] CRAN (R 4.3.3)
## highr 0.10 2022-12-22 [1] CRAN (R 4.3.1)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.2)
## htmltools 0.5.7 2023-11-03 [1] CRAN (R 4.3.2)
## htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.3.2)
## httpuv 1.6.13 2023-12-06 [1] CRAN (R 4.3.2)
## httr 1.4.7 2023-08-15 [1] CRAN (R 4.3.2)
## import 1.3.2 2024-01-21 [1] CRAN (R 4.3.3)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.1)
## jsonlite 1.8.8 2023-12-04 [1] CRAN (R 4.3.2)
## KernSmooth 2.23-22 2023-07-10 [2] CRAN (R 4.3.2)
## knitr 1.45 2023-10-30 [1] CRAN (R 4.3.2)
## labeling 0.4.3 2023-08-29 [1] CRAN (R 4.3.1)
## later 1.3.2 2023-12-06 [1] CRAN (R 4.3.2)
## lattice * 0.22-5 2023-10-24 [2] CRAN (R 4.3.2)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.3.2)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.2)
## lubridate * 1.9.3 2023-09-27 [1] CRAN (R 4.3.2)
## magrittr * 2.0.3 2022-03-30 [1] CRAN (R 4.3.1)
## markdown 1.12 2023-12-06 [1] CRAN (R 4.3.2)
## MASS 7.3-60.0.1 2024-01-13 [2] CRAN (R 4.3.2)
## Matrix 1.6-5 2024-01-11 [1] CRAN (R 4.3.2)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.1)
## mgcv 1.9-1 2023-12-21 [2] CRAN (R 4.3.2)
## mime 0.12 2021-09-28 [1] CRAN (R 4.3.1)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.3.2)
## mnormt 2.1.1 2022-09-26 [1] CRAN (R 4.3.1)
## munsell 0.5.1 2024-04-01 [1] CRAN (R 4.3.3)
## mvtnorm 1.3-1 2024-09-03 [1] CRAN (R 4.3.3)
## nlme 3.1-164 2023-11-27 [2] CRAN (R 4.3.2)
## patchwork 1.2.0 2024-01-08 [1] CRAN (R 4.3.2)
## pcaPP 2.0-4 2023-12-07 [1] CRAN (R 4.3.2)
## permute * 0.9-7 2022-01-27 [1] CRAN (R 4.3.2)
## pillar 1.10.2 2025-04-05 [1] CRAN (R 4.3.3)
## pkgbuild 1.4.3 2023-12-10 [1] CRAN (R 4.3.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.2)
## pkgload 1.3.4 2024-01-16 [1] CRAN (R 4.3.2)
## plotly 4.10.4 2024-01-13 [1] CRAN (R 4.3.2)
## plyr 1.8.9 2023-10-02 [1] CRAN (R 4.3.2)
## png 0.1-8 2022-11-29 [1] CRAN (R 4.3.1)
## profvis 0.3.8 2023-05-02 [1] CRAN (R 4.3.2)
## promises 1.2.1 2023-08-10 [1] CRAN (R 4.3.2)
## proxy 0.4-27 2022-06-09 [1] CRAN (R 4.3.2)
## psych 2.4.1 2024-01-18 [1] CRAN (R 4.3.2)
## purrr 1.0.2 2023-08-10 [1] CRAN (R 4.3.2)
## qvalue * 2.34.0 2023-10-24 [1] Bioconductor
## R6 2.6.1 2025-02-15 [1] CRAN (R 4.3.3)
## radiant.data * 1.6.3 2023-12-17 [1] CRAN (R 4.3.3)
## ragg 1.2.7 2023-12-11 [1] CRAN (R 4.3.2)
## randomizr 1.0.0 2023-08-10 [1] CRAN (R 4.3.3)
## raster * 3.6-26 2023-10-14 [1] CRAN (R 4.3.2)
## RColorBrewer * 1.1-3 2022-04-03 [1] CRAN (R 4.3.1)
## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.1)
## readr 2.1.5 2024-01-10 [1] CRAN (R 4.3.2)
## readxl 1.4.3 2023-07-06 [1] CRAN (R 4.3.2)
## remotes 2.5.0 2024-03-17 [1] CRAN (R 4.3.3)
## reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.3.2)
## rgdal * 1.6-7 2023-05-31 [1] CRAN (R 4.3.1)
## rlang 1.1.3 2024-01-10 [1] CRAN (R 4.3.2)
## rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.3.1)
## rnaturalearth * 1.0.1 2023-12-15 [1] CRAN (R 4.3.2)
## rnaturalearthdata 1.0.0 2024-02-09 [1] CRAN (R 4.3.2)
## robust * 0.7-4 2024-02-04 [1] CRAN (R 4.3.2)
## robustbase 0.99-2 2024-01-27 [1] CRAN (R 4.3.2)
## rrcov 1.7-5 2024-01-30 [1] CRAN (R 4.3.2)
## rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.2)
## sass 0.4.8 2023-12-06 [1] CRAN (R 4.3.2)
## scales * 1.3.0 2023-11-28 [1] CRAN (R 4.3.2)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.2)
## sf * 1.0-15 2023-12-18 [1] CRAN (R 4.3.2)
## shiny 1.8.0 2023-11-17 [1] CRAN (R 4.3.2)
## shinyAce 0.4.2 2022-05-06 [1] CRAN (R 4.3.3)
## shinyFiles 0.9.3 2022-08-19 [1] CRAN (R 4.3.3)
## sp * 2.1-3 2024-01-30 [1] CRAN (R 4.3.2)
## stringi 1.8.3 2023-12-11 [1] CRAN (R 4.3.2)
## stringr 1.5.1 2023-11-14 [1] CRAN (R 4.3.2)
## systemfonts 1.0.5 2023-10-09 [1] CRAN (R 4.3.2)
## terra * 1.7-71 2024-01-31 [1] CRAN (R 4.3.2)
## textshape * 1.7.5 2024-04-01 [1] CRAN (R 4.3.3)
## textshaping 0.3.7 2023-10-09 [1] CRAN (R 4.3.2)
## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.2)
## tidyr * 1.3.1 2024-01-24 [1] CRAN (R 4.3.2)
## tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.3.3)
## timechange 0.3.0 2024-01-18 [1] CRAN (R 4.3.2)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.3.2)
## units 0.8-5 2023-11-28 [1] CRAN (R 4.3.2)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.3.2)
## usethis 2.2.3 2024-02-19 [1] CRAN (R 4.3.2)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.2)
## vegan * 2.6-4 2022-10-11 [1] CRAN (R 4.3.2)
## viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.2)
## withr 3.0.2 2024-10-28 [1] CRAN (R 4.3.3)
## writexl * 1.5.0 2024-02-09 [1] CRAN (R 4.3.2)
## xfun 0.41 2023-11-01 [1] CRAN (R 4.3.2)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.2)
## yaml 2.3.8 2023-12-11 [1] CRAN (R 4.3.2)
##
## [1] C:/Users/tfrancisco/AppData/Local/R/win-library/4.3
## [2] C:/Users/tfrancisco/AppData/Local/Programs/R/R-4.3.2/library
##
## ──────────────────────────────────────────────────────────────────────────────