PCoA and PERMANOVA

Pre Requirements: The data should be structured as follows in a CSV file.

# Load required libraries

library(FD)            # For functional diversity and Gower distance calculations

library(ade4)          # For performing PCoA (Principal Coordinate Analysis)

library(vegan)         # For ecological analysis (PERMANOVA, Betadisper)

library(tidyverse)     # For data manipulation and plotting

library(ggrepel)       # For improved text label placement in ggplot2

library(cluster)       # For calculating Gower distance

library(RColorBrewer)  # For using color palettes in plots

# Load dataset

raptor <- read.csv("C:/Users/Rounak Choudhary/Desktop/raptor.csv",  # Read CSV file containing your data data

                   header=TRUE, sep=",", stringsAsFactors=TRUE, row.names=1)  # First column as row names, factors as strings

# Convert categorical variables to factors

raptor$Habitat <- as.factor(raptor$Habitat)     # Convert Habitat column to factor

raptor$Niche <- as.factor(raptor$Niche)         # Convert Niche column to factor

raptor$Lifestyle <- as.factor(raptor$Lifestyle) # Convert Lifestyle column to factor

raptor$Active <- as.factor(raptor$Active)       # Convert Active column to factor

# Compute Gower’s distance (handles mixed data types)

distances.raptor <- gowdis(raptor)              # Calculate Gower distance matrix suitable for mixed variables

# Perform Principal Coordinate Analysis (PCoA)

pcoa.all <- dudi.pco(sqrt(distances.raptor), scannf = FALSE, nf = 4)  # Perform PCoA on the square-root of the distance matrix, extract 4 axes

# Calculate percentage variance explained by each axis

eig_vals <- pcoa.all$eig                                  # Extract eigenvalues from PCoA result

variance_explained <- eig_vals / sum(eig_vals) * 100      # Calculate % variance explained

var1 <- round(variance_explained[1], 2)                   # Round variance of axis 1 to 2 decimal places

var2 <- round(variance_explained[2], 2)                   # Round variance of axis 2 to 2 decimal places

# Create a dataframe for plotting

pcoa_df <- data.frame(PCoA1 = pcoa.all$li[,1],            # Scores of axis 1

                      PCoA2 = pcoa.all$li[,2],            # Scores of axis 2

                      Species = rownames(raptor),         # Species names from row names

                      Lifestyle = raptor$Lifestyle)       # Lifestyle grouping to color points by Lifestyle

# Ensure Lifestyle is a factor

pcoa_df$Lifestyle <- as.factor(pcoa_df$Lifestyle)         # Make sure Lifestyle is treated as factor (categorical)

# Define colors using RColorBrewer

num_lifestyles <- length(unique(pcoa_df$Lifestyle))       # Count how many unique Lifestyle groups

palette <- brewer.pal(min(num_lifestyles, 8), "Set1")     # Pick appropriate color palette (max 8 colors from 'Set1')

# Plot the PCoA Biplot with species grouped by Lifestyle and variance explained on axes

pcoa_plot <- ggplot(pcoa_df, aes(x = PCoA1, y = PCoA2, color = Lifestyle, label = Species)) +  # Set plot aesthetics

  geom_point(size = 3) +                                # Add points for each species

  geom_text_repel(size = 3, max.overlaps = Inf, box.padding = 0.1, point.padding = 0.1) +  # Add labels with repel to avoid overlap

  scale_color_manual(values = palette) +                # Apply manual color palette

  geom_hline(yintercept = 0, linetype = "dashed", color = "black", linewidth = 0.8) +  # Draw horizontal dashed line at y=0

  geom_vline(xintercept = 0, linetype = "dashed", color = "black", linewidth = 0.8) +  # Draw vertical dashed line at x=0

  labs(title = "PCoA Biplot with Species Grouped by Lifestyle",  # Plot title

       x = paste0("PCoA Axis 1 (", var1, "%)"),                 # X-axis label with % variance

       y = paste0("PCoA Axis 2 (", var2, "%)")) +              # Y-axis label with % variance

  theme_minimal()                                              # Use minimal ggplot theme

# Save as high-resolution PNG

ggsave("PCoA_Biplot_HighRes.png", plot = pcoa_plot, width = 10, height = 8, dpi = 300)  # Save plot as high-res PNG file (10x8 inches, 300 dpi)

#Statistics

# PERMANOVA (adonis2 from vegan package)

# Perform PERMANOVA to test if "Lifestyle" explains variation in the trait space

permanova_result <- adonis2(distances.raptor ~ Lifestyle, data = raptor, permutations = 999, method = "euclidean")  

# View PERMANOVA results

print(permanova_result)  # Print the PERMANOVA result (R², F value, p-value)

# Betadisper (Homogeneity of Dispersion Test)

# Check if group dispersions are homogeneous (important assumption check for PERMANOVA)

dispersion <- betadisper(distances.raptor, raptor$Lifestyle)         # Calculate multivariate dispersion for Lifestyle groups

permutest(dispersion, permutations = 999)                            # Permutation test to check if dispersion is significantly different

# A non-significant result means group differences are not due to spread differences (good for PERMANOVA)

# Display PERMANOVA Result on PCoA Plot (Optional)

# Extract R² and p-value from PERMANOVA result to annotate on the plot

R2 <- round(permanova_result$R2[1], 3)                               # Extract and round R² (proportion of variance explained)

p_val <- permanova_result$`Pr(>F)`[1]                                # Extract p-value from PERMANOVA result

# Add PERMANOVA statistics to the PCoA plot

pcoa_plot1 <- ggplot(pcoa_df, aes(x = PCoA1, y = PCoA2, color = Lifestyle, label = Species)) +  # Set aesthetics

  geom_point(size = 3) +                                # Plot points

  geom_text_repel(size = 3, max.overlaps = Inf, box.padding = 0.1, point.padding = 0.1) +  # Add text labels with repel

  scale_color_manual(values = palette) +                # Apply color palette

  geom_hline(yintercept = 0, linetype = "dashed", color = "black", linewidth = 0.8) +  # Horizontal line at y=0

  geom_vline(xintercept = 0, linetype = "dashed", color = "black", linewidth = 0.8) +  # Vertical line at x=0

  labs(title = "PCoA Biplot with Species Grouped by Lifestyle",   # Title

       x = paste0("PCoA Axis 1 (", var1, "%)"),                  # X-axis label with variance explained

       y = paste0("PCoA Axis 2 (", var2, "%)")) +                # Y-axis label with variance explained

  annotate("text",                                               # Add annotation text (PERMANOVA result)

           x = max(pcoa_df$PCoA1),                               # Position text at maximum X value (right side)

           y = max(pcoa_df$PCoA2),                               # Position text at maximum Y value (top side)

           label = paste0("PERMANOVA: R² = ", R2, ", p = ", p_val),  # Text to display

           hjust = 1,                                            # Align text to the right

           size = 4) +                                           # Text size

  theme_minimal()                                                # Minimal ggplot theme

# Save as high-resolution PNG

ggsave("PCoA_Biplot_HighRes.png", plot = pcoa_plot1, width = 10, height = 8, dpi = 300)  # Save plot as high-res PNG file (10x8 inches, 300 dpi)