rm(list=ls())
library(readxl)
library(bipartite)
library(dplyr)
library(ggplot2)
library(ggrepel)

# Read the matrix from Excel (we use Sheet2 as you mentioned)
matriz <- read_excel("data-METANETWORKS-QS.xlsx", sheet = "aggregated") #change sheet for others

# Step 2: Convert to a numeric matrix
numerica <- as.matrix(matriz)

# Initialize the variable to store the highest modularity value
max_modularity <- -Inf
best_result <- NULL

# Step 3: Run DIRT_LPA_wb_plus 100 times
for (i in 1:100) {
  
  # Run DIRT_LPA_wb_plus to obtain the modules
  res <- DIRT_LPA_wb_plus(numerica)
  
  # Extract modularity
  modularity_value <- res$modularity
  
  # If the current modularity is the highest, store the result
  if (modularity_value > max_modularity) {
    max_modularity <- modularity_value
    best_result <- res
  }
}

# Step 4: Convert the best result to a 'moduleWeb' object for plotting
modulos <- convert2moduleWeb(numerica, best_result)

# Step 5: Plot the modular network using plotModuleWeb
plotModuleWeb(modulos)

# Step 4: Extract c and z values from the best result
cz_values <- czvalues(modulos, weighted = TRUE)
c_values <- cz_values$c
z_values <- cz_values$z

# Step 5: Calculate the thresholds (95th percentile) for c and z
c_threshold <- quantile(c_values, 0.95, na.rm = TRUE)
z_threshold <- quantile(z_values, 0.95, na.rm = TRUE)

# Show the calculated thresholds
cat("Threshold c (95th percentile):", c_threshold, "\n")
cat("Threshold z (95th percentile):", z_threshold, "\n")

# Step 6: Identify species that exceed the c and z thresholds
higher_c <- names(c_values)[c_values > c_threshold]
higher_z <- names(z_values)[z_values > z_threshold]

# Show connector and hub species
cat("\nConnector species (c > threshold):\n", higher_c, "\n")
cat("\nHub species (z > threshold):\n", higher_z, "\n")

# Step 7: Plot c and z values with thresholds
# Specify output file and high-quality parameters
tiff("cz_plot_pbirds_high_res.tiff", width = 2000, height = 2000, res = 300, compression = "lzw")
plot(c_values, z_values, pch = 16, xlab = "Among-module connectivity, c", ylab = "Within-module degree, z", xlim = c(0, 1), las = 1, cex = 0.8, 
     main = "Parana-Birds")
abline(v = c_threshold, col = "blue", lty = 2)
abline(h = z_threshold, col = "red", lty = 2)
text(c_values, z_values, labels = names(c_values), pos = 4, cex = 0.7)
legend("right", legend = c("Threshold c (95%)", "Threshold z (95%)"), cex = 0.7,
       col = c("blue", "red"), lty = 2, bty = "n")
# Quadrant labels
text(0.1, 2.05, "Module hubs", cex = 0.8, col = "blue")
text(0.9, 2.05, "Network hubs", cex = 0.8, col = "red")
text(0.1, -0.8, "Peripherals", cex = 0.8, col = "violet")
text(0.9, -0.8, "Connectors", cex = 0.8, col = "green")

# Close the graphic device
dev.off()

# Create the data frame
data <- na.omit(data.frame(
  c = c_values,
  z = z_values,
  name = names(c_values)
))

# Separate the data into peripherals and non-peripherals
data_peripherals <- data %>%
  filter(c < c_threshold & z < z_threshold)

data_main <- data %>%
  filter(!(c < c_threshold & z < z_threshold))

# Create the plot
cz_plot <- ggplot() +
  # Peripheral points without labels
  geom_point(data = data_peripherals, aes(x = c, y = z), size = 2, color = "gray") +
  # Other category points with labels
  geom_point(data = data_main, aes(x = c, y = z), size = 3) +
  geom_text_repel(data = data_main, aes(x = c, y = z, label = name), size = 3, max.overlaps = 20) +
  # Threshold lines
  geom_vline(xintercept = c_threshold, color = "blue", linetype = "dashed") +
  geom_hline(yintercept = z_threshold, color = "red", linetype = "dashed") +
  # Threshold annotations
  annotate("text", x = 0.9, y = -1, label = "Cut-off c (95%)", color = "blue", size = 4, hjust = 1) +
  annotate("text", x = 1.05, y = z_threshold + 0.2, label = "Cut-off z (95%)", color = "red", size = 4, hjust = 0) +
  # Quadrant names
  annotate("text", x = 0.2, y = 3.5, label = "Module hubs", color = "blue", size = 4) +
  annotate("text", x = 0.9, y = 3.5, label = "Network hubs", color = "red", size = 4) +
  annotate("text", x = 0.2, y = -1.2, label = "Peripherals", color = "purple", size = 4) +
  annotate("text", x = 0.9, y = -1.2, label = "Connectors", color = "green", size = 4) +
  # Labels and theme
  labs(
    x = "Among-module connectivity, c",
    y = "Within-module degree, z",
    title = "Aggregated metanetwork interaction roles"
  ) +
  theme_minimal() +
  scale_x_continuous(expand = c(0.1, 0.1)) +
  scale_y_continuous(expand = c(0.1, 0.1)) +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA)
  )

# Save as TIFF file
ggsave(
  filename = "aggregatedczvalues.tiff", 
  plot = cz_plot, 
  device = "tiff", 
  width = 10, height = 10, units = "in", 
  dpi = 300, compression = "lzw",
  bg = "white"
)