# =========================================================
# R version 4.3.2 (2023-10-31)  #
# Operating System: arm64    #
# Code for statistical analysis   #
# Supplement to: Gabriel Fénisse*, Manuel Chevalier, Odile Peyron, David Vincent Bekaert, Pierre-Henri Blard: Advancing Last Glacial Maximum paleoclimate reconstructions in Europe using pollen data: a multi-method (mega)biomization approach, EGUsphere [preprint], https://doi.org/10.5194/egusphere-2026-1590, 2026. #
# Contact: *Gabriel Fénisse (gabriel.fenisse@univ-lorraine.fr)[Centre de Recherches Pétrographiques et Géochimiques, Vandoeuvre-les-Nancy, 2025] #
# =========================================================

###Biome reconstruction in R
# Biome affinity score of pollen-based reconstructions 

############################################################
# 1. DATA IMPORT
############################################################

# Load pollen percentage dataset - Modern
pollen <- read.csv(
  "EMPD2-Set.txt",
  sep = ",",
  dec = ".",
  header = TRUE,
  check.names = FALSE
)
NameAssCal <- read.delim(
  "NameAssCal_Megabiome.txt",
  sep = "\t",
  header = TRUE,
  check.names = FALSE
)

taxa <- pollen[, -(1:2)]
colnames(taxa) <- NameAssCal$Names[-c(1,2)]

keep <- !grepl("Delete", colnames(taxa), ignore.case = TRUE)
taxa <- taxa[, keep, drop = FALSE]

rs <- rowSums(taxa, na.rm = TRUE)
taxa <- sweep(taxa, 1, ifelse(rs == 0, NA, rs), FUN = "/") * 100

merge_duplicated_cols <- function(df) {
  groups <- split(seq_along(names(df)), names(df))
  
  merged <- lapply(groups, function(idx) {
    if (length(idx) == 1) {
      df[[idx]]
    } else {
      rowSums(df[, idx, drop = FALSE], na.rm = TRUE)
    }
  })
  out <- as.data.frame(merged)
  names(out) <- names(groups)
  out
}
taxa <- merge_duplicated_cols(taxa)
names(taxa) <- gsub("\\.[0-9]+$", "", names(taxa))
pollen_final <- cbind(pollen[, 1:2], taxa)
pollen <- pollen_final

# Load pollen percentage dataset - Fossil
pollen <- read.csv2(
    "CoolingPollen.csv", # "CoolingPollen.csv", "AssFLDB-2.csv"
    sep = ";",
    dec = ","
  )
NameAssCal <- read.csv2(
  "./Name_CoolingPollen_Harmonized_Megabiome.txt", # "./Name_CoolingPollen_Harmonized_Megabiome.txt", "./NameLDB_Harmonized.csv"
  sep = ";",
  dec = "."
)
pollen[, -(1:2)] <- pollen[, -(1:2)] * 100
pollen_taxa <- pollen[, -(1:2)]
names(pollen_taxa) <- NameAssCal$Names
# merge duplicated taxa names
pollen_taxa <- as.data.frame(
  sapply(
    split.default(pollen_taxa, names(pollen_taxa)),
    rowSums,
    na.rm = TRUE
  )
)
# rebuild full dataset
pollen <- cbind(
  pollen[, 1:2],
  pollen_taxa
)

# Load Taxa → Plant Functional Types (PFT) assignment table
pft_Int <- read.csv2("Taxon_PFT_Europe.csv", sep=";", stringsAsFactors=FALSE, check.names=FALSE)  # In Europe 
pft_Int <- type.convert(pft_Int, as.is = TRUE) 
colnames(pft_Int)[1] <- "taxa"
#setdiff(colnames(pollen),pft_Int$taxa)

# Load PFTs-Biomes data
bio <- read.csv2("PFT_Megabiomes_Europe.csv", header=TRUE, stringsAsFactors=FALSE, check.names=FALSE, sep=";")     # In Europe
bio <- type.convert(bio, as.is = TRUE) 
colnames(bio)[1] <- "biome" 

# Load taxa-specific parameters (thresholds and weights)
variable <- read.csv2("pollen_variable.csv", dec=",", sep=";", stringsAsFactors=FALSE)   # In Europe
variable <- type.convert(variable, as.is = TRUE) 
colnames(variable)[1] <- "taxa" # Rename the first column to "taxa"

############################################################
# 2. HELPER FUNCTION
############################################################

maxN <- function(x, N=1){
  len <- length(x)
  if(N > len){
    warning('N greater than length(x).  Setting N=length(x)')
    N <- length(x)
  }
  sort(x,partial=len-N+1)[len-N+1]
} 

############################################################
# 3. BIOME SCORE CALCULATION (SITE LOOP)
############################################################

# Unique site identifiers
list_ID <- unique(pollen$ID)
# Initialize output table
resulttable_complete <- NULL

##########################################################
# 3.1 Subset site data
##########################################################

{
  for (ID_Value in list_ID) { #add loop for lots of sites
    
    print(paste0("+++++ SITE ",ID_Value," | (",which(list_ID == ID_Value),"/",length(list_ID),") +++++"))
    
    # Separate metadata (first 2 columns) and pollen data
    subset_site <- subset(pollen, ID == ID_Value) 
    
    # Ensure taxa are present in both dataset and calibration table
    pollen.no   <- nrow(variable)
    pollennames <- variable$taxa
    
    # Split site data into metadata (first 2 columns) and taxa-specific data
    subset_site_info <- subset_site[, c(1:2)]
    subset_site_taxa <- subset_site[, -c(1:2)]
      
    # clean duplicate suffix
    names(subset_site_taxa) <- gsub("\\.[0-9]+$", "", names(subset_site_taxa))
    
    # merge duplicates properly
    subset_site_taxa <- as.data.frame(
      t(sapply(split.default(subset_site_taxa, names(subset_site_taxa)),
               rowSums, na.rm = TRUE))
    )
    names(subset_site_taxa) <- gsub("\\.[0-9]+$", "", names(subset_site_taxa))
    
    # ensure full pollen structure
    missing <- setdiff(pollennames, names(subset_site_taxa))
    subset_site_taxa[missing] <- 0
    
    subset_site_taxa <- subset_site_taxa[, pollennames]
    
    #subset_site_taxa <- subset_site[ ,names(subset_site) %in% pollennames]
    # Update pollen taxa
    pollennames <- intersect(pollennames, colnames(subset_site_taxa))
    pollen.no <- length(pollennames)
    
    ##########################################################
    # 3.2 Thresholding + weighting of pollen data
    ##########################################################
    taxa_corrected <- data.frame(matrix(NA, nrow=dim(subset_site_taxa)[1], ncol=dim(subset_site_taxa)[2], byrow=TRUE))
    colnames(taxa_corrected) <- names(subset_site_taxa)
    
    ##########################################################
    # 3.3 Biome affinity scoring
    ##########################################################
    
    for(k in 1:pollen.no){
      a <- as.vector(unlist(subset_site_taxa[pollennames[k]]))
      a[(a < variable[k,"threshold"])] <- 0  
      taxa_corrected[ ,pollennames[k]] <- a * variable[k,"weight"] 
    } 
    
    sample.no <- nrow(subset_site) 
    info.no   <- ncol(subset_site_info) 
    
    biome.no   <- nrow(bio) # Number of biomes
    biomenames <- bio$biome # Biome names
    
    res         <- data.frame(matrix(NA, nrow=sample.no, ncol=info.no+biome.no+1, byrow=TRUE))
    res[,1:info.no] <- subset_site_info
    colnames(res)   <- c(names(subset_site_info), biomenames, c('Best'))  ##add "Best_Column", "Second_Best_Column", 
  
    # Loop over biomes
    for(j in 1:biome.no){
        
        # Identify the PFTs from the biome
        biome1 <- biomenames[j]                       
        pft1   <- colnames(bio)[which(bio[j, ] == 1)] 
       #pft=subset(pft_Int, pft_Int$taxa %in% pollennames)
        pft <- pft_Int[pft_Int$taxa %in% pollennames, ]

        # Identify taxa that belong to the selected PFTs
        if(length(pft1) > 1){
          
          # PFTs associated with biomes
          #pollen_pft <- pollennames[rowSums(pft[ ,which(names(pft) %in% pft1)]) > 0]
          #pollen_pft <- pft$taxa[rowSums(pft[, which(names(pft) %in% pft1), drop = FALSE]) > 0]
          pollen_pft <- pft$taxa[rowSums(pft[, pft1, drop = FALSE]) > 0]
          
        }else{
          
          # Identify pollen taxa associated with biome PFTs
          pollen_pft <- pollennames[pft[ ,which(names(pft) %in% pft1)] > 0]
          
        }  
        
        # Biome affinity score
        #cat("\nBiome:", biome1, "\n")
        #cat("Taxa used:\n")
        #active_taxa <- pollen_pft[ taxa_corrected[pollen_pft] > 0 ]
        #print(active_taxa)
        res[,biome1] <- rowSums(
          sqrt(taxa_corrected[, pollen_pft, drop = FALSE]),
          na.rm = TRUE
        )}
        
        #res[,biome1] <- sum(sqrt(taxa_corrected[pollen_pft])) }
     # 
        
    ##########################################################
    # 3.4 Best biome attribution
    ##########################################################
    vals <- as.numeric(res[, biomenames])
    res[, "Best"] <- biomenames[which.max(vals)]
    resulttable_complete <- rbind(resulttable_complete, res)                
    
  }}

#For modern samples  
# Save the final results as a CSV file
############################################################
# 4. ADD COORDINATES + CLIMATE VARIABLES
############################################################
  
# Load auxiliary datasets
set_data <- read.csv2("EMPD2-Set.txt", sep = ",", dec = ".")
count_data <- read.csv2(
    "EMPD2_count_v1.csv",
    sep = ";",
    dec = ".",
    fileEncoding = "latin1"
  )
  
climate_data <- read.csv2(
    "CalibrationDataset_EMPD2_B.csv",
    sep = ";",
    dec = ".",
    fileEncoding = "latin1"
  )
  
# Standardize site identifiers
set_data$Place <- trimws(set_data$Place)
count_data$Place <- trimws(count_data$Place)
climate_data$Place <- trimws(climate_data$Place)
  
# Ensure consistent column naming
names(count_data)[names(count_data) == "Place"] <- "Place"
names(climate_data)[names(climate_data) == "Place"] <- "Place"
  
############################################################
# Merge coordinates
############################################################
  
resulttable_complete <- merge(
    resulttable_complete,
    count_data[, c("Place", "Longitude", "Latitude")],
    by = "Place",
    all.x = TRUE
)
  
############################################################
# Merge climate variables (TANN, MTWA, MTCO)
############################################################
  
resulttable_complete <- merge(
    resulttable_complete,
    climate_data[, c("Place", "TANN", "MTWA", "MTCO")],
    by = "Place",
    all.x = TRUE
)

# Export
write.csv(
  resulttable_complete,
  "resulttable_complete_with_coords.csv",
  row.names = FALSE
)

# =========================================================
# END
# =========================================================

