rm(list = ls())
knitr::opts_chunk$set(
    message = FALSE,
    warning = FALSE,
    cache = FALSE
)

1 Introduction

In this script, we performed outlier detection using redundancy analysis (RDA), which does not account for population genetic structure, and partial redundancy analysis (pRDA), which does.

RDA is a multivariate canonical analysis. Unlike multiple regression, which models a single response variable, RDA handles a matrix of response variables. It is considered a canonical method because it produces canonical axes analogous to those obtained in PCA.

The procedure begins with linear regressions between the explanatory and response variables. PCA is then performed on the fitted values from these regressions to derive the canonical axes (RDA axes), which reduce dimensionality. These axes represent linear combinations of the explanatory variables. In summary, RDA combines features of PCA and multiple regression, producing PCA-like outputs based on fitted values, which allow interpretation under a linear relationship framework. Thus, RDA is a multivariate method based on a linear hypothesis.

pRDA is an extension of RDA that accounts for covariates in the analysis. The steps are similar to those of RDA, with the addition of an initial step that controls for variables that may influence either the response or explanatory variables. pRDA first regresses the response variables against the conditioning variables (covariates). The explanatory variables are also regressed against the same covariates. RDA is then performed on the residuals of both the response and explanatory variables. This procedure removes the effects of the conditioned variables, allowing a clearer assessment of the relationship between explanatory and response variables.

Sources: Legendre and Legendre 2012; Capblancq and Forester 2021; workshop on redundancy analysis by the Quebec Center for Biodiversity Science (https://r.qcbs.ca/workshop10/book-en/redundancy-analysis.html) and Brenna Forester tutorial (https://bookdown.org/hhwagner1/LandGenCourse_book/WE_11.html).

2 Data

The genomic data consist of allele frequencies at the population level (29 populations) from the imputed dataset (475 individuals, 8,616 SNPs) with minor allele count (MAC) correction.

#Genomic data
load("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  

#climatic/structure/IBD
load("Data/Species/Taxus_baccata/GEA_new_var/variance_partitioning/Climatic_data_RDA_pRDA.Rdata")

3 Models

We perform RDA using the genomic matrix corrected for MAC and the scaled climatic variables of interest. For pRDA, we also include the first two PC axes of the genetic PCA.

#model
RDA_env <- 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_env)
## $r.squared
## [1] 0.3262662
## 
## $adj.r.squared
## [1] 0.1425206
#other analysis to test the model
#summary(RDA_env)
#significativity of the model
#anova.cca(RDA_env)
#test the significativity of the axis and the climatic variables using permutation
#anova.cca(RDA_env, step = 1000, by = "axis")
#anova.cca(RDA_env, step = 1000, by = "term")
#model
pRDA_env <- rda(formula = genomic_matrix ~  Annual_Tc+Diurnal_range_Tc+Tc_Seasonality+Tc_driest_quarter+Annual_P+P_Seasonality + Condition(PC1+PC2), data = Climatic_data_RDA_pRDA, scale=F)
RsquareAdj(pRDA_env)
## $r.squared
## [1] 0.1801329
## 
## $adj.r.squared
## [1] 0.03423138
#other analysis to test the model
#summary(pRDA_env)
#significativity of the model
#anova.cca(pRDA_env)
#test the significativity of the axis and the climatic variables using permutation
#anova.cca(pRDA_env, step = 1000, by = "axis")
#anova.cca(pRDA_env, step = 1000, by = "term")

Next, we conducted candidate selection based on these models. To achieve this, we followed the procedure developed by Capblancq and Forester (2021) and identified outliers based on their positions along a distribution of Mahalanobis distances. These distances are estimated between each locus and the center of the RDA space using a specified number of canonical axes (K).

First, we need to determine the number of canonical axes (RDA axes) to retain.

list_models <- c("RDA","pRDA")

for(i in 1:length(list_models)){
  
  model <- get(paste0(list_models[i],"_env"))
#screeplot  
explained_variance <- model$CCA$eig * 100 / sum(model$CCA$eig)
  
scree_data <- data.frame(Axis = factor(1:length(explained_variance), 
                                       labels = paste0("RDA", 1:length(explained_variance))),
                         Variance = explained_variance)
# Plot with ggplot2
plot <- ggplot(scree_data, aes(x = Axis, y = Variance)) +
  geom_bar(stat = "identity", fill = "darkgrey") +
  labs(title = paste("Screeplot -", list_models[i]), 
       x = "RDA Axis", y = "Explained Variance (%)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),panel.background = element_blank())

print(plot)
#explained variance along each RDA axis
model$CCA$eig*100/sum(model$CCA$eig)

pdf(paste0("Results/species/taxus/GEA_new_var/Redundancy_analyses/screeplot_",list_models[i],".pdf"))
print(plot);dev.off()
}

An interesting aspect of redundancy analyses is that we can visualise populations (or individuals or genotypes, depending on the scale of the study), loci, and climatic variables in the same space. First, we plotted the populations and climatic variables in a biplot within the RDA space.

list_models <- c("RDA","pRDA")

for(i in 1:length(list_models)){
   model <- get(paste0(list_models[i],"_env"))
#score along the 2 first RDA axis
score_climatic_var <- as.data.frame(scores(model, choices=c(1:2), display="bp"))
score_pop_var <- data.frame(model$CCA$u[,c(1,2)])

#meta_data
meta_data_pop <- read.csv("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),]

#merge for country info
Score_population_bis <- rownames_to_column(score_pop_var,"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)

#explained variance along each RDA axis
explained_variance <- data.frame(model$CCA$eig)*100/sum(model$CCA$eig) # variance of each axis*100 divided by the sum of all the explained variance for all axis
explained_variance_round <- round(explained_variance$model.CCA.eig,digits=1)

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 = 3, alpha = 0.8) +
  geom_segment(data = score_climatic_var, aes(xend = RDA1*2, yend = RDA2*2, 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=2.1*RDA1, y=2.5*RDA2, label = row.names(score_climatic_var)), size = 3.5)+
  xlab(paste0("RDA 1 (",explained_variance_round[1],"%)")) + 
  ylab(paste0("RDA 2 (",explained_variance_round[2],"%)")) +
  ggtitle(paste0("Biplot ", list_models[i]," Populations")) +
  scale_color_manual(name = "Countries", values = group_palette, labels = levels(score_with_country_info$Country)) +
   scale_x_continuous(expand = expansion(mult = 0.14)) +
  scale_y_continuous(expand = expansion(mult = 0.1)) + 
  guides(color = guide_legend(override.aes = list(size = 3)))+ 
  theme_bw(base_size = 14) +
  theme(legend.position="right", panel.grid = element_blank(), strip.text = element_text(size=14),plot.title = element_text(hjust = 0.5,color = "Black",face="italic"),legend.text = element_text(size = 13),)+
  labs(color = "Country")

 print(biplot_populations)
 pdf(paste0("Results/species/taxus/GEA_new_var/Redundancy_analyses/biplot_",list_models[i],"_populations.pdf"));print(biplot_populations);dev.off()
}

4 Outlier identification

The next step is to run the genome scan to calculate the Mahalanobis distances. This will allow us to compute the associated p-values and q-values, which will be used as thresholds for candidate selection.

4.1 Calculating p-values and q-values

We used the rdadapt function from Capblancq and Forester (2021):

rdadapt<-function(rda,K)
{
  zscores<-rda$CCA$v[,1:as.numeric(K)]
  resscale <- apply(zscores, 2, scale)
  resmaha <- covRob(resscale, distance = TRUE, na.action= na.omit, estim="pairwiseGK")$dist
  lambda <- median(resmaha)/qchisq(0.5,df=K)
  reschi2test <- pchisq(resmaha/lambda,K,lower.tail=FALSE)
  qval <- qvalue(reschi2test)
  q.values_rdadapt<-qval$qvalues
  return(data.frame(p.values=reschi2test, q.values=q.values_rdadapt))
}
for(i in 1:length(list_models)){
  
   model <- get(paste0(list_models[i],"_env"))

#Perform the function to calculate the mahalanobis distance and then pvalues/qvalues
genome_scan <- rdadapt(model,K=2) 

assign(paste0("genome_scan_",list_models[i]), genome_scan)

#save the pvalues for each snp for retaining SNPs in LD with the greater signal in the outliers identification script
pvalues_model_snp <- data.frame(snp_names=colnames(genomic_matrix),pvalues= genome_scan$p.values)

assign(paste0("pvalues_", list_models[i],"_snp"), pvalues_model_snp)

   # Now save the object with a descriptive filename
   save(list = paste0("pvalues_", list_models[i],"_snp"), file = paste0("Data/Species/Taxus_baccata/genetic_new/GEA_new_var/pvalues_", list_models[i], "_snp.Rdata"))

load(paste0("Data/Species/Taxus_baccata/genetic_new/GEA_new_var/pvalues_",list_models[i],"_snp.Rdata"))

write_xlsx(get(paste0("pvalues_", list_models[i],"_snp")),paste0("Data/Species/Taxus_baccata/genetic_new/GEA_new_var/pvalues_",list_models[i],"_snp.xlsx"))

#plot the distribution of pvalues 
Histogram_of_Pvalues_RDA<- hist(genome_scan$p.values,
     main= paste0("Histogram of ",list_models[i]," P-values"),
     xlab= "P-values")

#save the histogram
 png(filename=paste0("Results/species/taxus/GEA_new_var/Redundancy_analyses/Histogram_of_Pvalues_",list_models[i],".png"));hist(genome_scan$p.values,
     main= paste0("Histogram of ",list_models[i]," P-values"),
     xlab= "P-values");dev.off()
}

4.2 Threshold selection

Using the calculated p-values and q-values, we can now proceed with candidate selection. This involves setting thresholds to identify significant loci based on their statistical significance. We used FDR 10% or 5% depending on the species (see Supplementary Methods).

for(i in 1:length(list_models)){
   genome_scan_model <- get(paste0("genome_scan_",list_models[i]))
# qvalues < 0.05
thres_qvalues <- 0.05

outliers_qvalues_5perc <- data.frame(Loci = colnames(genomic_matrix)[which(genome_scan_model$q.values<thres_qvalues)], p.value = genome_scan_model$q.values[which(genome_scan_model$q.values<thres_qvalues)], contig = unlist(lapply(strsplit(colnames(genomic_matrix)[which(genome_scan_model$q.values<thres_qvalues)], split = "_"), function(x) x[1])))

nrow(outliers_qvalues_5perc)

assign(paste0("outliers_qvalues_5perc_", list_models[i]), outliers_qvalues_5perc)

# qvalues < 0.10
thres_qvalues <- 0.10

outliers_qvalues_10perc <- data.frame(Loci = colnames(genomic_matrix)[which(genome_scan_model$q.values<thres_qvalues)], p.value = genome_scan_model$q.values[which(genome_scan_model$q.values<thres_qvalues)], contig = unlist(lapply(strsplit(colnames(genomic_matrix)[which(genome_scan_model$q.values<thres_qvalues)], split = "_"), function(x) x[1])))

nrow(outliers_qvalues_10perc)
assign(paste0("outliers_qvalues_10perc_", list_models[i]), outliers_qvalues_10perc)

#threshold pvalues < 0.05 with Bonferonni correction
thres_pvalues <- 0.05/length(genome_scan_model$p.values)

outliers_pvalues <- data.frame(Loci = colnames(genomic_matrix)[which(genome_scan_model$p.values<thres_pvalues)], p.value = genome_scan_model$p.values[which(genome_scan_model$p.values<thres_pvalues)], contig = unlist(lapply(strsplit(colnames(genomic_matrix)[which(genome_scan_model$p.values<thres_pvalues)], split = "_"), function(x) x[1])))

nrow(outliers_pvalues)
assign(paste0("outliers_pvalues_", list_models[i]), outliers_pvalues)
}

4.3 Graphic representations of the outliers

4.3.1 FDR 5%

We can visualise both candidate and presumed neutral loci within the RDA space, alongside the climatic variables.

for(i in 1:length(list_models)){

   model <- get(paste0(list_models[i],"_env"))
   outliers_qvalues_5perc <- get(paste0("outliers_qvalues_5perc_", list_models[i]))
score_loci <- as.data.frame(scores(model, choices=c(1:2), display="species", scaling="none"))
score_loci_outliers <- data.frame(names = row.names(score_loci), score_loci)
score_loci_outliers$FDR5 <- "other SNPs"
score_loci_outliers$FDR5[score_loci_outliers$names%in%outliers_qvalues_5perc$Loci] <- "candidate SNPs FDR 5%"
score_loci_outliers$FDR5 <- factor(score_loci_outliers$FDR5, levels = c("candidate SNPs FDR 5%","other SNPs"))
score_loci_outliers <- score_loci_outliers[order(score_loci_outliers$FDR5),]

score_climatic_var <- as.data.frame(scores(model, choices=c(1:2), display="bp"))

explained_variance <- data.frame(model$CCA$eig)*100/sum(model$CCA$eig) # variance of each axis*100 divided by the sum of all the explained variance for all axis
explained_variance_round <- round(explained_variance$model.CCA.eig,digits=1)

#Biplot with SNPs and climatic variables along the two first RDA axis. 
 biplot_outliers<- 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_outliers, aes(x=RDA1*15, y=RDA2*15,colour=FDR5), 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 = 3.7) +
  xlab(paste0("RDA 1 (",explained_variance_round[1],"%)")) + 
  ylab(paste0("RDA 2 (",explained_variance_round[2],"%)")) +
   ggtitle(paste0(list_models[i]," space: candidates FDR 5%")) +
  scale_color_manual(values = c("#F9A242FF","lightblue")) +
    scale_x_continuous(expand = expansion(mult = 0.18)) +  
  scale_y_continuous(expand = expansion(mult = 0.16)) +  
   guides(color = guide_legend(title="SNP type",override.aes = list(size = 4)))+ 
  theme_bw(base_size = 14) +
  theme(legend.position="right", panel.grid = element_blank(), strip.text = element_text(size=14),plot.title = element_text(hjust = 0.5,color = "Black",face="italic"))
 
 print(biplot_outliers)
 
 pdf(paste0("Results/species/taxus/GEA_new_var/Redundancy_analyses/biplot_outliers_",list_models[i],"_FDR_5perc.pdf"));print(biplot_outliers);dev.off()
}

4.3.2 FDR 10%

We can visualise the SNPs identified under a FDR threshold of 10%

for(i in 1:length(list_models)){

   model <- get(paste0(list_models[i],"_env"))
   outliers_qvalues_10perc <- get(paste0("outliers_qvalues_10perc_", list_models[i]))
score_loci <- as.data.frame(scores(model, choices=c(1:2), display="species", scaling="none"))
score_loci_outliers <- data.frame(names = row.names(score_loci), score_loci)
score_loci_outliers$FDR5_10 <- "other SNPs"
score_loci_outliers$FDR5_10[score_loci_outliers$names%in%outliers_qvalues_10perc$Loci] <- "candidate SNPs FDR 10%"
score_loci_outliers$FDR5_10 <- factor(score_loci_outliers$FDR5_10, levels = c("candidate SNPs FDR 10%","other SNPs"))
score_loci_outliers <- score_loci_outliers[order(score_loci_outliers$FDR5_10),]

score_climatic_var <- as.data.frame(scores(model, choices=c(1:2), display="bp"))

explained_variance <- data.frame(model$CCA$eig)*100/sum(model$CCA$eig)
explained_variance_round <- round(explained_variance$model.CCA.eig,digits=1)

#Biplot with SNPs and climatic variables along the two first RDA axis including the less conservative threshold
 biplot_outliers_RDA_LC<- 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_outliers, aes(x=RDA1*15, y=RDA2*15,colour=FDR5_10), size = 1.4) +
  geom_segment(data = score_climatic_var, aes(xend=RDA1, yend=RDA2, x=0, y=0), colour="black", size=0.3, 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 = 3.9) +
  xlab(paste0("RDA 1 (",explained_variance_round[1],"%)")) + 
  ylab(paste0("RDA 2 (",explained_variance_round[2],"%)")) +
   ggtitle(paste0(list_models[i]," space: candidates FDR 10%")) +
  scale_color_manual(values = c("orange","lightblue")) +
  scale_x_continuous(expand = expansion(mult = 0.18)) +  
  scale_y_continuous(expand = expansion(mult = 0.16)) +  
   guides(color = guide_legend(override.aes = list(size = 4)))+ 
  guides(color = guide_legend(title="SNP type",override.aes = list(size = 4)))+ 
  theme_bw(base_size = 14) +
  theme(legend.position="right", panel.grid = element_blank(), strip.text = element_text(size=14),plot.title = element_text(hjust = 0.5,color = "Black",face="italic"))
 
 print(biplot_outliers_RDA_LC)
 
  pdf(paste0("Results/species/taxus/GEA_new_var/Redundancy_analyses/biplot_outliers_",list_models[i],"_FDR_10perc.pdf"));print(biplot_outliers_RDA_LC);dev.off()
}

We can also visualise the top candidate loci, selected based on a false discovery rate (FDR) of 5%, alongside the presumably neutral loci within the RDA space.

for(i in 1:length(list_models)){

   model <- get(paste0(list_models[i],"_env"))
   outliers_qvalues_10perc <- get(paste0("outliers_qvalues_10perc_", list_models[i]))
   outliers_qvalues_5perc <- get(paste0("outliers_qvalues_5perc_", list_models[i]))

score_loci <- as.data.frame(scores(model, choices=c(1:2), display="species", scaling="none"))
score_loci_outliers <- data.frame(names = row.names(score_loci), score_loci)
score_loci_outliers$FDR5_10 <- "other SNPs"
score_loci_outliers$FDR5_10[score_loci_outliers$names%in%outliers_qvalues_10perc$Loci] <- "candidate SNPs FDR 10%"
score_loci_outliers$FDR5_10[score_loci_outliers$names%in%outliers_qvalues_5perc$Loci] <- "candidate SNPs FDR 5%"
score_loci_outliers$FDR5_10 <- factor(score_loci_outliers$FDR5_10, levels = c("candidate SNPs FDR 10%","candidate SNPs FDR 5%","other SNPs"))
score_loci_outliers <- score_loci_outliers[order(score_loci_outliers$FDR5_10),]

score_climatic_var <- as.data.frame(scores(model, choices=c(1:2), display="bp"))

explained_variance <- data.frame(model$CCA$eig)*100/sum(model$CCA$eig)
explained_variance_round <- round(explained_variance$model.CCA.eig,digits=1)

#Biplot with SNPs and climatic variables along the two first RDA axis including the less conservative threshold
 biplot_outliers_RDA_LC<- 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_outliers, aes(x=RDA1*15, y=RDA2*15,colour=FDR5_10), 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 = 3.7) +
  xlab(paste0("RDA 1 (",explained_variance_round[1],"%)")) + 
  ylab(paste0("RDA 2 (",explained_variance_round[2],"%)")) +
   ggtitle(paste0(list_models[i]," space: candidates FDR 10%")) +
  scale_color_manual(values = c("darkgreen","orange","lightblue")) +
    scale_x_continuous(expand = expansion(mult = 0.18)) +  
  scale_y_continuous(expand = expansion(mult = 0.16)) +  
    guides(color = guide_legend(title="SNP type",override.aes = list(size = 4)))+ 
  theme_bw(base_size = 14) +
  theme(legend.position="right", panel.grid = element_blank(), strip.text = element_text(size=14),plot.title = element_text(hjust = 0.5,color = "Black",face="italic"))
 
print(biplot_outliers_RDA_LC)
 
  pdf(paste0("Results/species/taxus/GEA_new_var/Redundancy_analyses/biplot_outliers_",list_models[i],"_FDR_10perc+5perc.pdf"));print(biplot_outliers_RDA_LC);dev.off()
}

4.4 Save candidates

We can save the candidates:

4.4.1 FDR 5% candidates

for(i in 1:length(list_models)){

  outliers_FDR5perc_T_Adapcon_gentree <- get(paste0("outliers_qvalues_5perc_", list_models[i]))[,-3]
  
  assign(paste0("outliers_",list_models[i],"_FDR5perc_T_Adapcon_gentree"),outliers_FDR5perc_T_Adapcon_gentree)
  
   save(list = paste0("outliers_",list_models[i],"_FDR5perc_T_Adapcon_gentree"), file = paste0("Results/species/taxus/GEA_new_var/outliers/outliers_",list_models[i],"FDR5perc_T_Adapcon_gentree_bis.Rdata"))

write_xlsx(get(paste0("outliers_",list_models[i],"_FDR5perc_T_Adapcon_gentree")),paste0("Results/species/taxus/GEA_new_var/outliers/outliers_",list_models[i],"_FDR5perc_T_Adapcon_gentree.xlsx"))
}

4.4.2 FDR 10% candidates

for(i in 1:length(list_models)){

  outliers_FDR10perc_T_Adapcon_gentree <- get(paste0("outliers_qvalues_10perc_", list_models[i]))[,-3]
  
  assign(paste0("outliers_",list_models[i],"_FDR10perc_T_Adapcon_gentree"),outliers_FDR10perc_T_Adapcon_gentree)
  
   save(list = paste0("outliers_",list_models[i],"_FDR10perc_T_Adapcon_gentree"), file = paste0("Results/species/taxus/GEA_new_var/outliers/outliers_",list_models[i],"FDR10perc_T_Adapcon_gentree_bis.Rdata"))

write_xlsx(get(paste0("outliers_",list_models[i],"_FDR10perc_T_Adapcon_gentree")),paste0("Results/species/taxus/GEA_new_var/outliers/outliers_",list_models[i],"_FDR10perc_T_Adapcon_gentree_bis.xlsx"))
}

Session info

packages <- c(
  "vegan", "dplyr", "robust", "ggplot2",
  "radiant.data", "textshape", "writexl"
)
info <- sessionInfo()
packages_used <- info$otherPkgs[packages]
data.frame(
  Package = names(packages_used),
  Version = sapply(packages_used, function(x) x$Version)
)
##                   Package Version
## vegan               vegan   2.6-4
## dplyr               dplyr   1.1.4
## robust             robust   0.7-4
## ggplot2           ggplot2   3.5.2
## radiant.data radiant.data   1.6.3
## textshape       textshape   1.7.5
## writexl           writexl   1.5.0