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

1 Introduction

BayPass is a genome-environment association (GEA) method developed by Mathieu Gautier (2015). These models explicitly account for (and can estimate) the covariance structure between population allele frequencies arising from the shared demographic history of populations. BayPass provides three models:
- Core model: an FST-based scan that does not incorporate climatic or environmental variables.
- Auxiliary covariate model: a GEA approach that does not account for population genetic structure.
- Standard covariate model: a GEA approach that accounts for population genetic structure.

In this study, we used the standard covariate model to perform outlier detection. This model estimates a scaled covariance matrix between populations, which is used as a covariate to account for population genetic structure. It evaluates the extent to which each marker is linearly associated with the covariates (i.e., the covariance matrix and environmental variables).

The regression coefficient for each SNP with respect to the covariates is estimated using either MCMC (not detailed here) or an importance sampling (IS) approximation. The IS approach also allows estimation of the Bayes factor, which quantifies support for the association between SNP i and covariate k, by comparing the model with association (βᵢₖ ≠ 0) to the null model (βᵢₖ = 0) (BayPass manual).

Steps to perform outlier detection with BayPass:
- Prepare genomic data (with and without MAC filtering).
- Prepare climatic/environmental data in the required format.
- Run the standard covariate model with IS on Linux.
- Analyse the results following the BayPass manual in R.

#meta data pop
meta_data_pop <- read.csv("Data/Species/Taxus_baccata/Populations/taxus_sample_29pop.csv",h=T,sep=";",dec=",")
#meta data indiv
meta_data_vcf=read.csv("Data/Species/Taxus_baccata/Samples/samples_taxus_baccata_adapcon_gentree.csv",h=T,sep=";",dec=",")

2 Data

The genomic data must be formatted with SNPs in rows, with one row per SNP, and populations in columns (the number of columns should be twice the number of populations), with two columns per population (one for each allele of a SNP). The genomic information is encoded in allele counts at the population level, as shown below:
POP1 POP2 POP3 POP4 POP5 POP6
SNP 1: 71 8 115 0 61 36 51 39 10 91 69 58
SNP 2: 82 0 91 0 84 14 24 57 28 80 18 80

We will need to create two genomic files in this format: one for the GEA and another to calculate the covariance matrix to account for the population genetic structure. The covariance matrix file will be filtered for MAC =< 3 (except for T. baccata), as MAC is important for assessing population genetic structure.

2.1 Format of the genomic data for the GEA analysis (corrected for MAC)

The genomic matrix consists of the imputed and MAC-corrected dataset (for GEA) or the dataset without MAC correction (for the omega matrix), comprising 475 individuals and 8,616 SNPs.

#genomic data for outlier detection
load("Data/Species/Taxus_baccata/genetic_new/Gen_matrix_imp_T_Adapcon_Gentree_475_8616.Rdata")
genomic_data_maf_c <-Gen_matrix_imp_T_Adapcon_Gentree_475_8616

#genomic data with MAC (non corrected) for covariance matrix
load("Data/Species/Taxus_baccata/genetic_new/GEA/Data_geno_no_maf_c_8252SNP_452IND.Rdata")

genomic_data_no_maf_c <- Data_geno_no_maf_c_8252SNP_452IND

genomic_data_maf_c_ID <- rownames_to_column(genomic_data_maf_c)
#name of VCF_ID to merge 
names(genomic_data_maf_c_ID)[names(genomic_data_maf_c_ID) == 'rowname'] <- 'VCF_ID'

#add the population information
genomic_data_MAF_pop <- data.frame(merge(meta_data_vcf[,c(1,3)],genomic_data_maf_c_ID, "VCF_ID"))

#This function enables to calculate the allelic count of the 1st allele of SNP from the genotypic data in format: 0, 1, 2. 
reformat_genotype <- function(allele) {
  allele1 <- case_when(
    allele == " 0" ~ 2,
    allele == " 1" ~ 1,
    allele == " 2" ~ 0,
    allele == "0" ~ 2,
    allele == "1" ~ 1,
    allele == "2" ~ 0,
    TRUE ~ NA_real_
  )
  return(allele1)
}
#Apply the reformat_genotype function to SNP columns and calculate Allele2. Allele1 and 2 in two different dataframe to merge them more easily after. 
#allele1 
df_allele1 <- genomic_data_MAF_pop %>%
  mutate(across(starts_with("Tbac"), ~ reformat_genotype(.))) 

#allele 2
df_allele2 <- df_allele1 %>%
  mutate_at(vars(starts_with("Tbac")), ~ 2 - .)

# Group by Population and summarise allele counts
#allele1
allele_counts_allele1 <- df_allele1 %>%
  group_by(Population) %>%
  summarize(
    across(starts_with("Tbac"), ~ sum(.)),
  )

#allele2
allele_counts_allele2 <- df_allele2 %>%
  group_by(Population) %>%
  summarize(
    across(starts_with("Tbac"), ~ sum(.)),
  )

#rename allele_counts_2 to merge them with allele1
allele_counts_allele2$Population <- paste0(allele_counts_allele2$Population,"_Allele_2")

#final dataset with both allele 1 and 2 
final_dtf <- rbind(allele_counts_allele1,allele_counts_allele2);row.names(final_dtf) <- final_dtf$Population

#order the population to have the format where allele1 and allele 2 for each pop are beside
final_r_order <- final_dtf[order(row.names(final_dtf)), ]

#allele in row and population*2 in columns side by side
data_allele_count_BAYPASS_MAF_c <- data.frame(t(final_r_order))

#export the data in txt
write.table(x=data_allele_count_BAYPASS_MAF_c,
  file = "Data/Species/Taxus_baccata/GEA/BAYPASS/data_allele_count_BAYPASS_MAF_c.txt", 
            sep = " ",
            row.names = F, 
            col.names = F)

The data_allele_count_BAYPASS_MAF_c file contains the genotypic data to run the BAYPASS standard covariate model with importance sampling (IS).

2.2 Format of the genomic data for the covariance matrix (non corrected for MAC)

genomic_data_maf_no_c_ID <- rownames_to_column(genomic_data_no_maf_c)
#name of VCF_ID to merge 
names(genomic_data_maf_no_c_ID)[names(genomic_data_maf_no_c_ID) == 'rowname'] <- 'VCF_ID'

#add the population information
genomic_data_no_MAF_pop <- data.frame(merge(meta_data_vcf[,c(1,3)],genomic_data_maf_no_c_ID, "VCF_ID"))

#allele1 
df_allele1_no_mac <- genomic_data_no_MAF_pop %>%
  mutate(across(starts_with("Tbac"), ~ reformat_genotype(.))) 

#allele 2
df_allele2_no_mac <- df_allele1_no_mac %>%
  mutate_at(vars(starts_with("Tbac")), ~ 2 - .)

# Group by Population and summarise allele counts
#allele1
allele_counts_allele1_no_mac <- df_allele1_no_mac %>%
  group_by(Population) %>%
  summarize(
    across(starts_with("Tbac"), ~ sum(.,na.rm = TRUE)),
  )

#allele2
allele_counts_allele2_no_mac <- df_allele2_no_mac %>%
  group_by(Population) %>%
  summarize(
    across(starts_with("Tbac"), ~ sum(.,na.rm = TRUE)),
  )

#rename allele_counts_2 to merge them with allele1
allele_counts_allele2_no_mac$Population <- paste0(allele_counts_allele2_no_mac$Population,"_Allele_2")

#final dataset with both allele 1 and 2 
final_dtf_no_mac <- rbind(allele_counts_allele1_no_mac,allele_counts_allele2_no_mac);row.names(final_dtf_no_mac) <- final_dtf_no_mac$Population

#order the population to have the format where allele1 and allele 2 for each pop are beside
final_r_order_no_mac <- final_dtf_no_mac[order(row.names(final_dtf_no_mac)), ]

#allele in row and population*2 in columns side by side
data_allele_count_BAYPASS_MAF_no_c <- data.frame(t(final_r_order_no_mac))

#export the data in txt
write.table(x=data_allele_count_BAYPASS_MAF_no_c,
  file = "Data/Species/Taxus_baccata/GEA/BAYPASS/data_allele_count_BAYPASS_MAF_no_c.txt", 
            sep = " ",
            row.names = F, 
            col.names = F)

data_allele_count_BAYPASS_MAF_no_c contains the genomic data used to run the core model of BAYPASS to calculate the omega matrix.

2.3 Climatic data

We load the non scale climatic data:

#climatic data 
climatic_data <- read.csv("Data/Species/Taxus_baccata/Climatic_data/new_selection/Past_new_6_Climatic_data_scale_df.csv",sep=";",dec=",")

climatic_data_BAYPASS <- data.frame(t(climatic_data))

We need to format the data for BayPass, which can be done in two ways: either by creating a separate text file for each climatic variable with populations in columns, or by creating a single text file containing all climatic variables in rows and populations in columns.

write.table(x=climatic_data_BAYPASS[-c(1,2),],
  file = "Data/Species/Taxus_baccata/GEA_new_var/BAYPASS/climatic_data_BAYPASS.txt", 
            sep = " ",
            row.names = F, 
            col.names = F)

climatic_data_BAYPASS contains the climatic data used to run the standard covariate model with IS in BAYPASS.

3 Covariate model

We ran the core model to estimate the covariance matrix used to correct for population structure. We used genomic data without MAC correction: 8,252 SNPs, 29 populations, with allele counts at the population level (data_allele_count_BAYPASS_MAF_no_c). The core model was run with the following parameters:
i) 20 pilot runs of 500 iterations (to adjust proposal distributions)
ii) a burn-in period of 5000 iterations
iii) final MCMC sampling of 1000 parameter values sampled every 20 iterations (i.e., a total of 20 000 iterations)

The covariance matrix used to correct for population genetic structure is scaled directly within the model, so no additional scaling is required.

The output results are:

omega <- as.matrix(read.table("Data/Species/Taxus_baccata/GEA/BAYPASS/covari_matrix_core_model/out_mat_omega.out"))
country_names <- climatic_data$Country
dimnames(omega)=list(country_names,country_names)

Now, we can visualise the covariance matrix using plots, corrplot, and heatmaps, as proposed by Mathieu Gautier (BAYPASS Manual) and following Archambeau et al. (2025).

plot.omega <- function(omega, PC = c(1,2), country, main = expression(Omega*" matrix"), 
                       group_palette, pch = 16, expand_factor = 1.5) {
  om.svd <- svd(omega)
  eig <- om.svd$d
  pcent.var <- 100 * eig / sum(eig)
  # Assign colors based on country
  col <- group_palette[country]
  # Get min and max values for PC axes
  x_range <- range(om.svd$u[, PC[1]])
  y_range <- range(om.svd$u[, PC[2]])
  # Expand range symmetrically
  x_expand <- diff(x_range) * (expand_factor - 1) / 2
  y_expand <- diff(y_range) * (expand_factor - 1) / 2
  xlim <- x_range + c(-x_expand, x_expand)
  ylim <- y_range + c(-y_expand, y_expand)

  # Adjust margins for legend space
  par(mar = c(5, 5, 4, 6))
  plot(om.svd$u[, PC], main = main, pch = pch, col = col, 
       xlim = xlim, ylim = ylim,  
       xlab = paste0("PC", PC[1], " (", round(pcent.var[PC[1]], 2), "%)"),
       ylab = paste0("PC", PC[2], " (", round(pcent.var[PC[2]], 2), "%)"),
       cex.lab = 1.5, cex.axis = 1.3, cex.main = 1.8  
  )
  # Fix legend position 
  legend("topright", inset = c(-0.15, 0), legend = names(group_palette), 
         col = group_palette, pch = 16, title = "Country", 
         cex = 1, bty = "n", xpd = TRUE, ncol = 2)  
}
dev.new(width = 10, height = 8)  

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")

#order country in the matrix
country_vector <- c("Italy", "Switzerland","Slovenia","Norway","Spain","Spain","Spain","UK","Italy","Italy","Slovenia","Slovakia","Switzerland","Switzerland","Greece","Greece","Greece","Germany","Spain","UK","France","Sweden","France","Spain","Slovenia","Italy","Bosnia","Switzerland","UK") # 

# Run the function
plot.omega(omega, PC = c(1,2), country = country_vector, group_palette = group_palette)

# as a correlation plot
cor_mat <- cov2cor(omega)
corrplot::corrplot(cor_mat)

# corrplot(cor_mat,method="color",mar=c(2,1,2,2)+0.1,
# main=expression("Correlation map based on"~hat(Omega)))

# as a heatmap and hierarchical clustering tree (using the average agglomeration method)
##we use the population names
pop_names <- climatic_data$Population
dimnames(omega)=list(pop_names,pop_names)

#cor matrix 
cor_mat <- cov2cor(omega)

hclust_ave <- function(x) hclust(x, method="average")
Heatmap_omega_matrix<-heatmap(1-cor_mat,hclustfun = hclust_ave,
main=expression("Heatmap of "~hat(Omega)~"("*d[ij]*"=1-"*rho[ij]*")"))

4 standard covariate model

4.1 Run under Linux

We can run the IS models multiple times to ensure that the outputs are consistent across runs. The models are executed using the IS method via the Linux command-line interface of BAYPASS. We will interpret the association between SNPs and climatic predictors using the Bayes factor (BF)

4.2 Results

First, we extracted the BF from the BAYPASS output

#names snp
names_snps <- colnames(genomic_data_maf_c)

#name climatic data
  climatic_variables <- colnames(climatic_data)
  COVARIABLE <- c(1,2,3,4,5,6)
climatic_variables_merge <- data.frame(climatic_variables[-c(1,2)], COVARIABLE)

#only the name of the variables 
names_climatic_variables <- climatic_variables[-c(1,2)]

for(x in 1:5){
  seed <- x+11
  output_BAYPASS_run <- read.table(paste0("Results/species/taxus/GEA_new_var/BAYPASS/runs/out_BAYPASS_outliers_ide_",seed,"_summary_betai_reg.out"),h=T)

  BAYPASS_clim <- merge(climatic_variables_merge,output_BAYPASS_run,"COVARIABLE")

BAYPASS_clim$climatic_variables..c.1..2.. <- as.factor(BAYPASS_clim$climatic_variables..c.1..2..)

subset_name <- paste0("RUN_",x, "_BAYPASS_results")

  Data_results <- BAYPASS_clim %>% 
    group_by(climatic_variables..c.1..2..) %>% 
    mutate(names_snps = names_snps) %>%
    ungroup()
  
  assign(subset_name, Data_results)
}

names_subset_data_frame <- c("RUN_1_BAYPASS_results","RUN_2_BAYPASS_results","RUN_3_BAYPASS_results","RUN_4_BAYPASS_results","RUN_5_BAYPASS_results")

We also added the run info for BF.

## we also add the RUN info
    for (x in 1:5) {
      var <- paste0("RUN_",x,"_BAYPASS_results")
    # Create the new dataframe name
    names_final <- paste0(var, "_final")
    
    # Extract the dataframe using get
    data <- get(var)
    
    name_BF.dB. <- paste0("BF_RUN",x)
    # Mutate the data
    names(data)[names(data)== "BF.dB."] <- name_BF.dB.
    
    # Assign the mutated data to a new dataframe
    assign(names_final, data)
    }

Finally, we merged all the runs in on dataframe

#merge all runs
data_tot_results_allRUNs <- cbind(RUN_1_BAYPASS_results_final[,c(12,2,8)],RUN_2_BAYPASS_results_final[,c(8)],RUN_3_BAYPASS_results_final[,c(8)],RUN_4_BAYPASS_results_final[,c(8)],RUN_5_BAYPASS_results_final[,c(8)])

4.3 Results across runs

Next, we assessed whether the runs produced consistent results. To do this, we calculated the correlation of the Bayes factors (BF) across runs for each climatic variable.

First, we subset the dataset of all runs by climatic variable to create one dataset per variable.

#subset the results at the climatic variable scale to compare values between runs
for(var in names_climatic_variables){
  subset_name <- paste0(var, "_BAYPASS_all_runs") #name the output dataframe of the loop
  assign(subset_name, subset(data_tot_results_allRUNs, climatic_variables..c.1..2.. == var))#assign the name to the dataframe in a loop
}

#list of dataset with data of all run for each climatic data
final_names_dataframe <- c("Annual_Tc_BAYPASS_all_runs","Diurnal_range_Tc_BAYPASS_all_runs","Tc_Seasonality_BAYPASS_all_runs","Tc_driest_quarter_BAYPASS_all_runs","Annual_P_BAYPASS_all_runs","P_Seasonality_BAYPASS_all_runs")
#name climatic data
names_climatic_variables <- colnames(climatic_data[,-c(1,2)])

Then, we wanted to see if these values were consistent across runs by testing the correlation of the values across runs:

#correlation between runs for BF 
list_title <- c("BIO 1","BIO 2","BIO 4","BIO 9","BIO 12", "BIO 15")

for(x in 1:6){
  var <- final_names_dataframe[x]
  title_name <- list_title[x]
  #names of the corr matrix for each biovariable
  names_corr <- paste0("correlation_",var,"_BF")
  #title of corrplot
  title <- (paste0(title_name," all runs"))
  
  #group for each bio var with the values of Bayes factor only 
corr_bio <- cor(get(var)[, grepl("BF", names(get(var)))]) 

#name the corr_bio with names_corr
assign(names_corr,corr_bio)

#plot corrplot
corr_plot <- corrplot(get(names_corr), method = "number", addrect = 2, col = c("red", "white", "red"), type = "lower", tl.col = "black", tl.cex = 1, number.cex = 1, title = title,mar=c(0,0,1,0) )

#save
pdf(paste0("Results/species/taxus/GEA_new_var/BAYPASS/figures/correlation_BF_values_across_runs_",names_climatic_variables[x],".pdf"));corrplot(get(names_corr), method = "number", addrect = 2, col = c("red", "white", "red"), type = "lower", tl.col = "black", tl.cex = 1, number.cex = 1, title = title,mar=c(0,0,1,0) );dev.off()
}

5 Candidate detection

5.1 Thresholds

5.1.1 Bayes factor values

To select candidate SNPs, we will set a threshold based on the Bayes factor (BF). Following Mathieu Gautier (2015), we will apply Jeffreys’ rule (Jeffreys, 1961) to quantify the strength of evidence in favour of the association between the SNPs and the covariates, using the following unit scale:
- 10 < BF < 15 = strong evidence
- 15 < BF < 20 = very strong evidence
- BF > 20 = conclusive evidence

We used threshold of BF > 8 or 10 depending on the species.

To begin this selection process, we calculate the mean BF values across runs for each SNP.

final_names_dataframe <- c("Annual_Tc_BAYPASS_all_runs","Diurnal_range_Tc_BAYPASS_all_runs","Tc_Seasonality_BAYPASS_all_runs","Tc_driest_quarter_BAYPASS_all_runs","Annual_P_BAYPASS_all_runs","P_Seasonality_BAYPASS_all_runs")

#mean BF values 
for(data in final_names_dataframe){
  #names of the corr matrix for each biovariable
  names_corr <- paste0("mean_",data,"_BF_values")
  
  #select only the BF values
  dataset <- get(data)[, grepl("BF", names(get(data)))]
  
  #mean
  mean_BF_values <- data.frame(names_snps,rowMeans(dataset)); colnames(mean_BF_values)=c("names_snps","BF_values")
    
#name the corr_bio with names_corr
assign(names_corr,mean_BF_values)
}
list_mean_BF_clim <- c("mean_Annual_Tc_BAYPASS_all_runs_BF_values","mean_Diurnal_range_Tc_BAYPASS_all_runs_BF_values","mean_Tc_Seasonality_BAYPASS_all_runs_BF_values","mean_Tc_driest_quarter_BAYPASS_all_runs_BF_values","mean_Annual_P_BAYPASS_all_runs_BF_values","mean_P_Seasonality_BAYPASS_all_runs_BF_values")

Next, we will identify candidates by applying a Bayes factor threshold of 10 for each climatic variable.

thres_BF <- 10

for(x in 1:6){
  #names of the corr matrix for each biovariable
  names_climatic_var <- names_climatic_variables[x]
  data <- get(list_mean_BF_clim[x])
  
 selection_outliers <- data.frame(Loci=data$names_snps[which(data$BF_values>thres_BF)],BF = data$BF_values[which(data$BF_values>thres_BF)], climatic_variables=names_climatic_var)
 
 assign(paste0("names_outliers_",names_climatic_var),selection_outliers)
 # Count the number of candidates for each climatic variable
  count <- data.frame(Climatic_variable=names_climatic_var,Number_outliers=nrow(selection_outliers))
assign(paste0("outliers_",names_climatic_var),count)
}
  # Combine the results for each climatic variable to see the number of candidates
  outliers_summary_BF10 <- rbind(outliers_Annual_Tc,outliers_Diurnal_range_Tc,outliers_Tc_Seasonality,outliers_Tc_driest_quarter,outliers_Annual_P,outliers_P_Seasonality)

#name of the candidates
  outliers_names_summary_BF_10 <- rbind(names_outliers_Annual_Tc,names_outliers_Diurnal_range_Tc,names_outliers_Tc_Seasonality,names_outliers_Tc_driest_quarter,names_outliers_Annual_P,names_outliers_P_Seasonality)

  #finally we can search for same candidates across climatic variables
duplicated_loci_BF10 <- duplicated(outliers_names_summary_BF_10$Loci) | duplicated(outliers_names_summary_BF_10$Loci, fromLast = TRUE)

# Subset the data frame to show only the duplicated loci
duplicated_outliers_BF10 <- outliers_names_summary_BF_10[duplicated_loci_BF10, ]

#Number of candidates: 
nrow(outliers_names_summary_BF_10)

## [1] 67

#at the climatic var scale
outliers_summary_BF10

##   Climatic_variable Number_outliers
## 1         Annual_Tc               9
## 2  Diurnal_range_Tc               7
## 3    Tc_Seasonality              23
## 4 Tc_driest_quarter               8
## 5          Annual_P               3
## 6     P_Seasonality              17

Finally, we saved the set of candidates for downstream analysis

outliers_T_adapcon_gentree_BAYPASS_BF_10 <- outliers_names_summary_BF_10[,-2]

#write_xlsx(outliers_T_adapcon_gentree_BAYPASS_BF_10,"Results/species/taxus/GEA_new_var/outliers/outliers_T_adapcon_gentree_BAYPASS_BF_10.xlsx")
save(outliers_T_adapcon_gentree_BAYPASS_BF_10, file="Results/species/taxus/GEA_new_var/outliers/outliers_T_adapcon_gentree_BAYPASS_BF_10.Rdata")

We can also create a less conservative dataset with a threshold of BF>8:

thres_BF <- 8

for(x in 1:6){
  #names of the corr matrix for each biovariable
  names_climatic_var <- names_climatic_variables[x]
  data <- get(list_mean_BF_clim[x])
  
 selection_outliers <- data.frame(Loci=data$names_snps[which(data$BF_values>thres_BF)],BF = data$BF_values[which(data$BF_values>thres_BF)], climatic_variables=names_climatic_var)
 
 assign(paste0("names_outliers_",names_climatic_var),selection_outliers)
 
 # Count the number of candidates for each climatic variable
  count <- data.frame(Climatic_variable=names_climatic_var,Number_outliers=nrow(selection_outliers))
  
  #name the corr_bio with names_corr
assign(paste0("outliers_",names_climatic_var),count)
}
  # Combine the results for each climatic variable to see the number of candidates
  outliers_summary_BF_8 <- rbind(outliers_Annual_Tc,outliers_Diurnal_range_Tc,outliers_Tc_Seasonality,outliers_Tc_driest_quarter,outliers_Annual_P,outliers_P_Seasonality)

#name of the candidates
  outliers_names_summary_BF8 <- rbind(names_outliers_Annual_Tc,names_outliers_Diurnal_range_Tc,names_outliers_Tc_Seasonality,names_outliers_Tc_driest_quarter,names_outliers_Annual_P,names_outliers_P_Seasonality)
  
#finally we can search for same candidates across climatic variables
duplicated_loci_BF_8 <- duplicated(outliers_names_summary_BF8$Loci) | duplicated(outliers_names_summary_BF8$Loci, fromLast = TRUE)

# Subset the data frame to show only the duplicated loci
duplicated_outliers_BF_8 <- outliers_names_summary_BF8[duplicated_loci_BF_8, ]

#Number of candidates: 
outliers_summary_BF_8

##   Climatic_variable Number_outliers
## 1         Annual_Tc              13
## 2  Diurnal_range_Tc              17
## 3    Tc_Seasonality              35
## 4 Tc_driest_quarter              20
## 5          Annual_P              13
## 6     P_Seasonality              35

#at the climatic var scale
nrow(outliers_names_summary_BF8)

## [1] 133