########################################
#### Genetic indicators in Bulgaria ####
########################################

## This is the script supporting the papar
# "Are we ready for applying genetic diversity indicators? A cross-taxa pilot assessment of wild animal and plant populations in Bulgaria"

## Authors: Violeta Zhelyazkova & Alicia Mastretta-Yanes 
# To start with, we used the code provided here: https://github.com/AliciaMstt/GeneticIndicators

# Check the R version
R.version

# remove all objects previously existing
rm(list = ls())

# set the working directory (has to be changed by each user)
setwd("D:/Science/Genetic indicators")

# load necessary libraries
library(tidyr)
library(dplyr)
library(utile.tools)
library(stringr)
library(janitor)
library(ggplot2)
library(viridis)
library(cowplot)
library(ggsankey)
library(forcats)


# load the custom functions (also attached as supplementary files; they need to be in the same working directory)
source("get_indicatorNe_data_BG.R")
source("transform_to_Ne_BG.R")
source("estimate_indicatorNe_BG.R")

# not in
"%!in%" <- function(x, y)!('%in%'(x,y))


########### PART 1: DATA FORMATTING AND INDICATOR CALCULATION 


# select the file to use
kobo_file="Supplement1_Raw_data_BG.csv"


# open the file and define taxon
kobo_output<-read.csv(file=kobo_file, sep=",", header=TRUE)

# View(kobo_output)


## format the file

kobo_BG<-kobo_output %>% 
  
  # create a taxon column
  mutate(taxon=(utile.tools::paste(genus, species, subspecies_variety, na.rm=TRUE))) %>%
  
  # remove white space at the end of the name
  mutate(taxon=str_trim(taxon, "right")) %>% 
  
  # extract the year of assessment
  mutate(year_assesment=substr(end,1,4)) %>%
  
  # replace empty cells or answers -999 with NA
  mutate(across(everything(),
                ~ {
                  if (is.character(.) || is.factor(.)) {
                    x <- as.character(.)
                    x <- na_if(x, "")
                    na_if(x, "-999")
                  } else if (is.numeric(.)) {
                    na_if(., -999)
                  } else {
                    .
                  }
                }
  )) %>%
  
  # remove columns entirely composed of NA
  remove_empty(which = "cols", cutoff = 1, quiet = F) %>% 
  
  # remove columns that you don't need
  dplyr::select(-start, -end, -username) %>% 
  
  # create a new column that includes the custom Nc/Ne ratios when provided
  # at the end, this was not used in the paper
  mutate(ratio_formatted = case_when(ratio_species_related == "0.06 - 0.14" ~ (0.06+0.14)/2, # these are manually corrected
                                     ratio_species_related == "0.24-0.78" ~ (0.24+0.78)/2,   # to allow calculations
                                     ratio_species_related == "<0.3" ~ 0.29,
                                     ratio_species_related == "0,16" ~ 0.16,
                                     # if it is an integer → copy it
                                     suppressWarnings(!is.na(as.integer(ratio_species_related))) ~ 
                                       as.numeric(ratio_species_related),
                                     TRUE ~ 0.1))

# check if there are indeed NAs in ratio_formatted
unique(is.na(kobo_BG$ratio_formatted))

# check column names so that "_" doesn't interfere with the function get_indicatorNe_data_BG()
names(kobo_BG) <- gsub("transboundary_percentage","transboundaryPercent",names(kobo_BG))
names(kobo_BG) <- gsub("transboundary_countries","transboundaryCountries",names(kobo_BG))
  

# create a shorter dataframe that contains just the species info and not population info
kobo_BG_short <- kobo_BG %>% 
  select(taxon, genus:source_national_taxonID, year_assesment, 
         defined_populations:source_popsize_ratios, ratio_formatted,
         species_comments, gen_studies:source_status_distribution)

# colnames(kobo_BG_short)

# View(kobo_BG_short)

# get Ne data
indNe_BG <- get_indicatorNe_data_BG(kobo_BG)

# Check if some taxons were discarded (this is if they didn't have any population information)
n_distinct(kobo_BG$taxon)
n_distinct(indNe_BG$taxon)

# check which are these taxons
View(kobo_BG[which(kobo_BG$taxon %!in% indNe_BG$taxon),])


# estimate total country population size for each taxon

national_NcPoint <- indNe_BG %>% 
  group_by(taxon) %>% 
  filter(all(!is.na(NcPoint))) %>% 
  summarise(national_pop_size <- sum(NcPoint))

# check for how many species national population size is available
nrow(national_NcPoint)

# View(national_NcPoint)


# calculate Ne from Nc (this creates a separate row for each population)

indNe_BG_red <-transform_to_Ne_BG(indNe_BG, ratio = 0.1) 
# View(indNe_BG_red)
nrow(indNe_BG_red) # this is the total number of populations

# select columns to use for per population calculations and to export for a supplement
Pop_numbers <- indNe_BG_red %>% 
  select(taxonomic_group, taxon, population, Name, Origin, 
         transboundary, transboundaryPercent, transboundaryCountries,
         NcRange, NcRangeDetails, NcPoint, NcLower, NcUpper, NcType, Ne, Ne_combined_standart) %>% 
  arrange(taxonomic_group, taxon)
  

# View(Pop_numbers)

nrow(Pop_numbers)

# more formatting
Pop_numbers$transboundary[is.na(Pop_numbers$transboundary)] <- "unknown"

# export the file
# write.csv(Pop_numbers, "Supplement7_Population-data.csv")

# estimate the Ne indicator (the proportion of populations with Ne > 500)
Ne_BG <- estimate_indicatorNe_BG(indNe_data_red = indNe_BG_red)
# View(Ne_BG)
nrow(Ne_BG)

# View(Ne_BG)

# Add Ne and the national population size to the metadata
ready_data_BG <- kobo_BG_short %>%
  left_join(Ne_BG, by = "taxon") %>%
  left_join(national_NcPoint, by = "taxon")

# View(ready_data_BG)

# estimate the PM indicator (proportion of populations maintained) and add it to the main data
ready_data_BG$PM_BG <- (ready_data_BG$n_extant_populations / (ready_data_BG$n_extant_populations + ready_data_BG$n_extinct_populations))
# View(ready_data_BG)

# Add the data on the Bulgarian Red Book and the Anexes of the Bulgarian Biodiversity Act
# that was already manually updated

RB_BDA <- read.csv("Supplement2_Red-Book-BDA.csv")

View(RB_BDA)
# View(ready_data_BG)

ready_data_BG_merged <- merge(ready_data_BG, RB_BDA, by.x = "taxon", by.y = "Taxon")

# View(ready_data_BG)


# format the data and add new columns
ready_data_BG_merged <- ready_data_BG_merged %>% 
  mutate(taxonomic_group = factor(taxonomic_group,
                                  levels = c("mammal", "bird", "reptile", "amphibian", "fish", "invertebrate",
                                             "angiosperm", "gymnosperm", "pteridophytes"),
                                  labels = c("Mammal", "Bird", "Reptile", "Amphibian", "Fish", "Invertebrate",
                                             "Angiosperm", "Gymnosperm", "Pteridophyte"))) %>% 
  mutate(gen_studies = factor(gen_studies,
                              levels = c("phylo", "pop", "phylo_pop", "no"),
                              labels = c("Phylgen", "Popgen", "Both", "No data"))) %>% 
  mutate(ne_pops_exists = factor(ne_pops_exists,
                                 levels = c("ne_available", "other_genetic_info", "no_genetic_data"),
                                 labels = c("Ne available", "Other genetic data", "No data"))) %>% 
  mutate(temp_gen_monitoring = factor(temp_gen_monitoring,
                                      levels = c("yes", "no", "unknown"),
                                      labels = c("Yes", "No", "No data"))) %>% 
  mutate(Nc_type_species = factor(Nc_type_species,
                                  levels = c("Point_data", "Range_data", "Mixed data", "No data"),
                                  labels = c("Point", "Range", "Mixed", "No data"))) %>% 
  mutate(transboundary_populations = factor(transboundary_populations,
                                            levels = c("yes", "no", "not_considered"),
                                            labels = c("Yes", "No", "Not considered"))) %>% 
 
  mutate(n_pops_factor = case_when(n_pops == 1 ~ "1",
                                   n_pops == 2 ~ "2",
                                   n_pops %in% 3:9 ~ "3-9",
                                   n_pops %in% 11:25 ~ "11-25")) %>% 
  mutate(n_pops_factor = factor(n_pops_factor,
                                levels = c("1", "2", "3-9", "11-25", "No data"))) %>% 

  mutate(indicatorNe_standart = round(indicatorNe_standart, 2),
         PM_BG = round(PM_BG, 2)) %>% 
  mutate(In_which_Annex_of_th_s_the_species_listed = factor(In_which_Annex_of_th_s_the_species_listed,
                                                            levels = c("annex_i_birds_directive",
                                                                       "annex_i_birds_directive annex_ii_birds_directive",
                                                                       "annex_i_birds_directive annex_ii_birds_directive annex_iii_birds_directive",
                                                                       "annex_ii_birds_directive",
                                                                       "annex_ii_habitat_directive",
                                                                       "annex_ii_habitat_directive annex_iv_habitat_directive",
                                                                       "annex_ii_habitat_directive annex_v_habitat_directive",
                                                                       "annex_iv_habitat_directive",
                                                                       "annex_v_habitat_directive" ),
                                                            labels = c("BD, Anex 1",
                                                                       "BD, Anex 1 & 2",
                                                                       "BD, Nex 1, 2 & 3",
                                                                       "BD, Anex 2",
                                                                       "HD, Anex 2",
                                                                       "HD, Anex 2 & 4",
                                                                       "HD, Anex 2 & 5",
                                                                       "HD, Anex 4",
                                                                       "HD, Anex 5")))


ready_data_BG_merged$ne_pops_exists[is.na(ready_data_BG_merged$ne_pops_exists)] <- "No data"
ready_data_BG_merged$Nc_type_species[is.na(ready_data_BG_merged$Nc_type_species)] <- "No data"
ready_data_BG_merged$n_pops_factor[is.na(ready_data_BG_merged$n_pops_factor)] <- "No data"
ready_data_BG_merged$PM_BG <- fct_explicit_na(as.factor(ready_data_BG_merged$PM_BG), na_level = "No data")
ready_data_BG_merged$indicatorNe_standart <- fct_explicit_na(as.factor(ready_data_BG_merged$indicatorNe_standart), na_level = "No data")
ready_data_BG_merged$RedBook <- fct_explicit_na(as.factor(ready_data_BG_merged$RedBook), na_level = "Not in RB")
ready_data_BG_merged$In_which_Annex_of_th_s_the_species_listed <- fct_explicit_na(as.factor(ready_data_BG_merged$In_which_Annex_of_th_s_the_species_listed), na_level = "Not in Dir")

# save the file
# write.csv(ready_data_BG_merged, "Ready_data_BG_merged_20260927.csv")


#### Data summaries ######


# check the possible population definitions (to list in the Methods)
levels(as.factor(ready_data_BG_merged$defined_populations))

# number of assessed species
n_distinct(ready_data_BG_merged$taxon)

# number of assessed taxonomic groups
n_distinct(ready_data_BG_merged$taxonomic_group)

# number of species per taxonomic groups
summary(ready_data_BG_merged$taxonomic_group)

# number of species with just one or two populations
ready_data_BG_merged %>% filter(n_extant_populations == "1" | n_extant_populations == "2") %>% nrow()

# species with > 10 populations
ready_data_BG_merged %>% mutate(n_extant_populations = as.numeric(n_extant_populations)) %>% 
  filter(n_extant_populations > 10) %>% summarise(unique(taxon))


## PM indicator summary

# per taxonomic group (not used in the manuscript)
PM_taxonomic_group_summary <- ready_data_BG_merged %>%
  group_by(taxonomic_group) %>%
  summarise(
    Total = n(),
    `PM = 1 (n)` = sum(PM_BG == "1", na.rm = TRUE),
    `PM < 1 (n)` = sum(PM_BG != "1" & PM_BG != "No data" & PM_BG != "0", na.rm = TRUE),
    `PM = 0 (n)` = sum(PM_BG == "0", na.rm = TRUE),
    `No data (n)` = sum(PM_BG == "No data", na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    `PM = 1 (%)` = round(100 * `PM = 1 (n)` / Total),
    `PM < 1 (%)` = round(100 * `PM < 1 (n)` / Total),
    `PM = 0 (%)` = round(100 * `PM = 0 (n)` / Total),
    `No data (%)` = round(100 * `No data (n)` / Total)
  ) %>%
  select(
    `Taxonomic group` = taxonomic_group,
    `PM = 1 (n)`, `PM = 1 (%)`,
    `PM < 1 (n)`, `PM < 1 (%)`,
    `PM = 0 (n)`, `PM = 0 (%)`,
    `No data (n)`, `No data (%)`,
    Total
  )

# View(PM_taxonomic_group_summary)


# for all taxonomic groups togehter (used in the abstract)

PM_summary <- ready_data_BG_merged %>%
  summarise(
    Total = n(),
    `PM = 1 (n)` = sum(PM_BG == "1", na.rm = TRUE),
    `PM < 1 (n)` = sum(PM_BG != "1" & PM_BG != "No data" & PM_BG != "0", na.rm = TRUE),
    `PM = 0 (n)` = sum(PM_BG == "0", na.rm = TRUE),
    `No data (n)` = sum(PM_BG == "No data", na.rm = TRUE)
  ) %>%
  mutate(
    `PM = 1 (%)` = round(100 * `PM = 1 (n)` / Total),
    `PM < 1 (%)` = round(100 * `PM < 1 (n)` / Total),
    `PM = 0 (%)` = round(100 * `PM = 0 (n)` / Total),
    `No data (%)` = round(100 * `No data (n)` / Total)
  ) %>%
  select(
    `PM = 1 (n)`, `PM = 1 (%)`,
    `PM < 1 (n)`, `PM < 1 (%)`,
    `PM = 0 (n)`, `PM = 0 (%)`,
    `No data (n)`, `No data (%)`,
    Total
  )

# View(PM_summary)

## Ne indicator summary

# per taxonomic group (not used in the manuscript)
Ne_taxonomic_group_summary <- ready_data_BG_merged %>%
  group_by(taxonomic_group) %>%
  summarise(
    Total = n(),
    `Ne = 1 (n)` = sum(indicatorNe_standart == "1", na.rm = TRUE),
    `Ne < 1 (n)` = sum(indicatorNe_standart != "1" & indicatorNe_standart != "No data", na.rm = TRUE),
    `No data (n)` = sum(indicatorNe_standart == "No data", na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    `Ne = 1 (%)` = round(100 * `Ne = 1 (n)` / Total),
    `Ne < 1 (%)` = round(100 * `Ne < 1 (n)` / Total),
    `No data (%)` = round(100 * `No data (n)` / Total)
  ) %>%
  select(
    `Taxonomic group` = taxonomic_group,
    `Ne = 1 (n)`, `Ne = 1 (%)`,
    `Ne < 1 (n)`, `Ne < 1 (%)`,
    `No data (n)`, `No data (%)`,
    Total
  )

# View(Ne_taxonomic_group_summary)


# for all taxonomic groups togehter (used in the abstract)

Ne_summary <- ready_data_BG_merged %>%
  summarise(
    Total = n(),
    `Ne = 1 (n)` = sum(indicatorNe_standart == "1", na.rm = TRUE),
    `Ne < 1 (n)` = sum(indicatorNe_standart != "1" & indicatorNe_standart != "No data" & indicatorNe_standart != "0", na.rm = TRUE),
    `Ne = 0 (n)` = sum(indicatorNe_standart == "0", na.rm = TRUE),
    `No data (n)` = sum(indicatorNe_standart == "No data", na.rm = TRUE)
  ) %>%
  mutate(
    `Ne = 1 (%)` = round(100 * `Ne = 1 (n)` / Total),
    `Ne < 1 (%)` = round(100 * `Ne < 1 (n)` / Total),
    `Ne = 0 (%)` = round(100 * `Ne = 0 (n)` / Total),
    `No data (%)` = round(100 * `No data (n)` / Total)
  ) %>%
  select(
    `Ne = 1 (n)`, `Ne = 1 (%)`,
    `Ne < 1 (n)`, `Ne < 1 (%)`,
    `Ne = 0 (n)`, `Ne = 0 (%)`,
    `No data (n)`, `No data (%)`,
    Total
  )

# View(Ne_summary)

# count the number of species with genetic studies
ready_data_BG_merged %>% filter(gen_studies != "No data") %>% summarise(n_distinct(taxon))

# look at how many of the taxons where genetic studies were conducted are mammals
ready_data_BG_merged %>% filter(gen_studies != "No data") %>% 
  filter(taxonomic_group == "Mammal") %>%  summarise(n_distinct(taxon))

## calculate whether population numbers were given as range or point data

# number of populations with range data
Pop_numbers %>% filter(NcType == "Nc_range") %>% nrow()
177/253 # as a percentage of the total number of populations

# number of species with range data
Pop_numbers %>% filter(NcType == "Nc_range") %>% summarise(n_distinct(taxon))
71/96 # as a percentage of the total number of species

# number of populations with point data
Pop_numbers %>% filter(NcType == "Nc_point") %>% nrow()
60/253 # as a percentage of the total number of populations

# number of species with pint data
Pop_numbers %>% filter(NcType == "Nc_point") %>% summarise(n_distinct(taxon))
26/96 # as a percentage of the total number of species

# number of populations with no data
Pop_numbers %>% filter(NcType == "unknown") %>% nrow()
16/253 # as a percentage of the total number of populations


## transboundary populations

# total number of transboundary populations
Pop_numbers %>% filter(transboundary == "yes_transboundary") %>% nrow()
79/253 # as a percentage of the total number of populations

# number of species with transboundary populations
Pop_numbers %>% filter(transboundary == "yes_transboundary") %>% summarise(n_distinct(taxon))
61/96 # as a percentage of the total number of species

# see the species for which Ne data is available
ready_data_BG_merged %>% 
  filter(ne_pops_exists == "Ne available") %>% 
  select(taxon)

# see species where genetic monitoring has been conducted
ready_data_BG_merged %>% 
  filter(temp_gen_monitoring == "Yes") %>% 
  select(taxon)


# summarise conservation status
summary(as.factor(ready_data_BG_merged$RedBook))


# Create a table with the values of both indicators per species

ready_data_BG_merged <- ready_data_BG_merged %>%
  rename(national_pop_size = `national_pop_size <- sum(NcPoint)`)

Indicators_species <- ready_data_BG_merged %>% 
  select(taxonomic_group, taxon, n_extant_populations, n_extinct_populations, defined_populations, 
         national_pop_size, transboundary_populations, n_pops_more_500_standart,
         indicatorNe_standart, PM_BG, IUCN, In_which_Annex_of_th_s_the_species_listed, RedBook, BDA) %>% 
  arrange(taxonomic_group, taxon)


colnames(Indicators_species)<- c("Group", "Taxon", "ExtantPop", "ExtinctPop", "PopDefiniton",
                    "NatPopSize", "Transboundary", "BigPops", "Ne", "PM", "IUCN", "EUDir", "RedBook", "BDA")

View(Indicators_species)

# view species where information on extinct populations is missing
Indicators_species %>% filter(PM == "No data") %>% select(Taxon, ExtinctPop, Ne) %>% View()

# export the table
# write.csv(Indicators_species, "Table1_Indicators_species.csv")


##### Figure 1 #####

dat <- ready_data_BG_merged %>% 
  mutate(Ne_cathegory = case_when(indicatorNe_standart == "1" ~ "1",
                                  indicatorNe_standart != "1" & indicatorNe_standart != "0" & indicatorNe_standart != "No data" ~ "0 < Ne < 1",
                                  indicatorNe_standart == "0" ~ "0",
                                  indicatorNe_standart == "No data" ~ "No data")) %>% 
  mutate(Ne_cathegory = factor(Ne_cathegory,
                               levels = c("1", "0", "0 < Ne < 1", "No data"))) %>% 
  mutate(PM_cathegory = case_when(PM_BG == "1" ~ "1",
                                  PM_BG != "1" & indicatorNe_standart != "0" & indicatorNe_standart != "No data" ~ "0 < PM < 1",
                                  PM_BG == "0" ~ "0",
                                  PM_BG == "No data" ~ "No data")) %>% 
  mutate(PM_cathegory = factor(PM_cathegory,
                               levels = c("1", "0", "0 < PM < 1", "No data")))

gen_studies <- ggplot(dat, aes(x=taxonomic_group, fill = as.factor(gen_studies))) + 
  geom_bar(stat = "count", position = position_stack(reverse = TRUE)) +
  scale_fill_viridis(discrete=TRUE) +
  guides(fill = guide_legend(nrow = 1)) +
  ylim(0,40) +
  theme_bw() +
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6),
        legend.position = "top",
        legend.box.margin = margin(b = -10),
        legend.title = element_blank(),
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.3, "cm"),
        plot.title = element_text(hjust = 0.5, size = 8)) +
        ggtitle("Genetic studies conducted")


Ne_av <- ggplot(dat, aes(x=taxonomic_group, fill = as.factor(ne_pops_exists))) + 
  geom_bar(stat = "count", position = position_stack(reverse = TRUE)) +
  scale_fill_viridis(discrete=TRUE, na.translate = F) +
  guides(fill = guide_legend(nrow = 1)) +
  ylim(0,40) +
  theme_bw() +
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6),
        legend.position = "top",
        legend.box.margin = margin(b = -10),
        legend.title = element_blank(),
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.3, "cm"),
        plot.title = element_text(hjust = 0.5, size = 8)) +
  ggtitle("Effective population size data")


gen_mon <- ggplot(dat, aes(x=taxonomic_group, fill = as.factor(temp_gen_monitoring))) + 
  geom_bar(stat = "count", position = position_stack(reverse = TRUE)) +
  scale_fill_viridis(discrete=TRUE, na.translate = F) +
  guides(fill = guide_legend(nrow = 1)) +
  ylim(0,40) +
  theme_bw() +
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6),
        legend.position = "top",
        legend.box.margin = margin(b = -10),
        legend.title = element_blank(),
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.3, "cm"),
        plot.title = element_text(hjust = 0.5, size = 8)) +
  ggtitle("Genetic monitoring conducted")


Nc_data <- ggplot(dat, aes(x=taxonomic_group, fill = as.factor(Nc_type_species))) + 
  geom_bar(stat = "count", position = position_stack(reverse = TRUE)) +
  scale_fill_viridis(discrete=TRUE, na.translate = F) +
  guides(fill = guide_legend(nrow = 1)) +
  ylim(0,40) +
  theme_bw() +
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6),
        legend.position = "top",
        legend.box.margin = margin(b = -10),
        legend.title = element_blank(),
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.3, "cm"),
        plot.title = element_text(hjust = 0.5, size = 8)) +
  ggtitle("Census population size data")

# transboundary populations
trans_data <- ggplot(dat, aes(x=taxonomic_group, fill = transboundary_populations)) + 
  geom_bar(stat = "count", position = position_stack(reverse = TRUE)) +
  scale_fill_viridis(discrete=TRUE, na.translate = F) +
  guides(fill = guide_legend(nrow = 1)) +
  ylim(0,40) +
  theme_bw() +
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6),
        legend.position = "top",
        legend.box.margin = margin(b = -10),
        legend.title = element_blank(),
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.3, "cm"),
        plot.title = element_text(hjust = 0.5, size = 8)) +
  ggtitle("Transboundary populations")

BDA_data <- ggplot(dat, aes(x=taxonomic_group, fill = BDA_Anex)) + 
  geom_bar(stat = "count", position = position_stack(reverse = TRUE)) +
  scale_fill_viridis(discrete=TRUE, na.translate = F) +
  guides(fill = guide_legend(nrow = 1)) +
  ylim(0,40) +
  theme_bw() +
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6),
        legend.position = "top",
        legend.box.margin = margin(b = -10),
        legend.title = element_blank(),
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.3, "cm"),
        plot.title = element_text(hjust = 0.5, size = 8)) +
  ggtitle("Bulgarian Biodiversity Act")

n_pop <- ggplot(dat, aes(x=taxonomic_group, fill = n_pops_factor)) + 
  geom_bar(stat = "count", position = position_stack(reverse = TRUE)) +
  scale_fill_viridis(discrete=TRUE, na.translate = F) +
  guides(fill = guide_legend(nrow = 1)) +
  ylim(0,40) +
  theme_bw() +
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6),
        legend.position = "top",
        legend.box.margin = margin(b = -10),
        legend.title = element_blank(),
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.3, "cm"),
        plot.title = element_text(hjust = 0.5, size = 8)) +
  ggtitle("Number of populations")

Ne_value <- ggplot(dat, aes(x=taxonomic_group, fill = Ne_cathegory)) + 
  geom_bar(stat = "count", position = position_stack(reverse = TRUE)) +
  scale_fill_viridis(discrete=TRUE, na.translate = F) +
  guides(fill = guide_legend(nrow = 1)) +
  ylim(0,40) +
  theme_bw() +
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6),
        legend.position = "top",
        legend.box.margin = margin(b = -10),
        legend.title = element_blank(),
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.3, "cm"),
        plot.title = element_text(hjust = 0.5, size = 8)) +
  ggtitle("Ne 500 value")

PM_value <- ggplot(dat, aes(x=taxonomic_group, fill = PM_cathegory)) + 
  geom_bar(stat = "count", position = position_stack(reverse = TRUE)) +
  scale_fill_viridis(discrete=TRUE, na.translate = F) +
  guides(fill = guide_legend(nrow = 1)) +
  ylim(0,40) +
  theme_bw() +
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6),
        legend.position = "top",
        legend.box.margin = margin(b = -10),
        legend.title = element_blank(),
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.3, "cm"),
        plot.title = element_text(hjust = 0.5, size = 8)) +
  ggtitle("PM value")

windows(height = 1200, width = 800)
plot_grid(Nc_data, Ne_av, n_pop, gen_studies, trans_data, gen_mon, Ne_value, PM_value,
          nrow = 4, ncol = 2)
ggsave("Figure1-final-20260715.png", dpi = 600, bg = "white")



##### Figure 2 #####

# View(ready_data_BG_merged)

# select columns you need for the graph
sankey_dat <- ready_data_BG_merged %>% 
  select(taxonomic_group, indicatorNe_standart, RedBook, PM_BG, taxonomic_group)
# View(sankey_dat)
sankey_dat$taxonomic_group2 <- sankey_dat$taxonomic_group

sankey_dat <- sankey_dat %>% 
  mutate(RedBook = factor(RedBook,
                                    levels = c("LC", "NT", "VU", "EN", "CR", "RE", "Not in RB"))) %>% 
  mutate(indicatorNe_standart = factor(indicatorNe_standart,
                                       levels = c("1", "0.6", "0.5", "0.33", "0.27",
                                                  "0.25", "0.14", "0.12", "0", "No data"),
                                       labels = c("1", "0.1-0.6", "0.1-0.6", "0.1-0.6", "0.1-0.6", "0.1-0.6", "0.1-0.6", "0.1-0.6", "0", "No data"))) %>% 
  mutate(PM_BG = factor(PM_BG,
                        levels = c("1", "0.96", "0.8", "0.67",
                                   "0.5", "0.25", "0.2", "0", "No data"),
                        labels = c("1", "0.5-0.9", "0.5-0.9", "0.5-0.9", "0.5-0.9",
                                   "0.2-0.3", "0.2-0.3", "0", "No data")))


sankey_long <-sankey_dat %>% make_long(taxonomic_group,
                                      indicatorNe_standart,
                                      RedBook,
                                      PM_BG,
                                      taxonomic_group2)


sankey_grouped <- sankey_long %>% dplyr::group_by(node)%>%tally

sankey_red<-merge(sankey_long, sankey_grouped, by.x="node", by.y="node", all.x=TRUE)
# View(sankey_red)

labels_sankey<-c("Taxonomic group",
                 "Ne 500 value",
                 "Bulgarian Red Book",
                 "PM value",
                 "Taxonomic group")


sankey_red$node <- factor(sankey_red$node, 
                             levels = c("0", "0.1-0.6", "0.2-0.3", "0.5-0.9", "1", "No data",
                                        "Pteridophyte", "Gymnosperm","Angiosperm", "Invertebrate",
                                        "Fish", "Amphibian", "Reptile", "Bird", "Mammal",
                                        "RE", "CR", "EN", "VU", "NT", "LC", "Not in RB"))



sankey_plot <-
  ggplot(sankey_red, aes(x=x,
                                   next_x=next_x,
                                   node=node,
                                   next_node=next_node,
                                   fill=factor(node),
                                   label=paste0(node, " ")))+
  geom_sankey(flow.alpha=0.4,
              node.color="black",
              width = 0.05,
              show.legend = FALSE)+
  geom_sankey_label(size=4, color = "black", fill = "white")+
  scale_x_discrete(labels = labels_sankey)+
  theme_bw()+
  scale_fill_viridis(discrete = T)+
  theme(axis.title = element_blank(),
        axis.text.x = element_text(size = 15, colour = "black"),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.grid = element_blank())


windows(height = 500, width = 900)
sankey_plot
ggsave("Figure_2-final-20260715-new.png", dpi = 600, bg = "white")


# see the species where discrepancies exist between national conservation status and indicator values

levels(as.factor(Indicators_species$RedBook))

# species with unfavourable values of the indicators but favourable Red Book conservation status
Indicators_species_discrepancies <- Indicators_species %>% 
  filter(PM != "1" | Ne != "1") %>% 
  filter(RedBook == "LC" | RedBook == "Not in RB") %>% 
  select(Group, Taxon, ExtantPop, NatPopSize, Ne, PM, Transboundary, RedBook)

View(Indicators_species_discrepancies)

# species with favourable values of the indicators but unfavourable Red Book conservation status
Indicators_species_discrepancies_2 <- Indicators_species %>% 
  filter(PM == "1" & Ne == "1") %>% 
  filter(RedBook != "LC" & RedBook != "Not in RB") %>% 
  select(Group, Taxon, ExtantPop, NatPopSize, Ne, PM, Transboundary, RedBook)

View(Indicators_species_discrepancies_2)

# species that are listed in the Bulgarian Red Book as endangered or critically endangered but not in EU directives
Indicators_species_discrepancies_3 <- Indicators_species %>% 
  filter(RedBook == "EN" | RedBook == "CR") %>% 
  filter(EUDir == "Not in Dir")

View(Indicators_species_discrepancies_3)