diversity_agroecoseqC

rm(list = ls()) 
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

spXplot <- readRDS("species_x_plots.RDS")
plotXsp <- as.data.frame(t(spXplot))
plotXenv <- readRDS("plot_X_env.RDS")
dataset_names <- c("CS1_CREA", "CS2_INRAE_CF", "CS3_INRAE", "CS5_CSIC_INIA", "CS7_LAMMC", "CS8_Tagem", "CS9_CRAW")

PLOT NON AGGREGATI, ossia subplot

# Esegui la fusione basata sui rownames
agroecoseqC_2.1 <- merge(plotXenv, plotXsp, by = "row.names", all = TRUE)

# Se desideri ripristinare i rownames originali
rownames(agroecoseqC_2.1) <- agroecoseqC_2.1$Row.names
agroecoseqC_2.1$Row.names <- NULL  # Rimuovi la colonna aggiuntiva dei rownames

# Funzione per trasformare le colonne (tranne le prime due) in numerico
trasforma_colonne_numeriche <- function(df) {
  df <- df %>%
    mutate_at(vars(-(1:3)), as.numeric)
  return(df)
}

# lista_dataframes <- lapply(lista_dataframes, trasforma_colonne_numeriche)

lista_dataframes <- split(agroecoseqC_2.1, agroecoseqC_2.1$exp)
list2env(lista_dataframes, envir = .GlobalEnv)

## <environment: R_GlobalEnv>

library(vegan)

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.6-8

# Elenco dei tuoi dataframe
dataframes <- list(CS1_CREA = CS1_CREA, CS2_INRAE_CF = CS2_INRAE_CF, CS3_INRAE = CS3_INRAE, CS5_CSIC_INIA = CS5_CSIC_INIA, CS7_LAMMC = CS7_LAMMC, CS8_Tagem = CS8_Tagem, CS9_CRAW = CS9_CRAW)

# Calcola gli indici per ciascun dataframe e crea un nuovo dataframe con i risultati
for(name in names(dataframes)){
  # Escludi le prime tre colonne
  df <- dataframes[[name]][, -c(1:5)] #elimino le colonne non numeriche
  
  # Calcola gli indici
  shannon_index <- diversity(df, index = "shannon")
  simpson_index <- diversity(df, index = "simpson")
  evenness <- shannon_index / log(specnumber(df))
  richness <- species_richness <- rowSums(df > 0)
  
  # Crea un nuovo dataframe con le prime tre colonne del dataframe originale e i nuovi indici
  new_df_name <- paste0(name, "_div")
  assign(new_df_name, data.frame(dataframes[[name]][, 1:3], Shannon_Index = shannon_index, Simpson_Index = simpson_index, Evenness = evenness, Species_richness = richness), envir = .GlobalEnv)
}


merged_df_div <- bind_rows(CS1_CREA_div, CS2_INRAE_CF_div, CS3_INRAE_div,
                       CS5_CSIC_INIA_div, CS7_LAMMC_div, CS8_Tagem_div, CS9_CRAW_div)

plotXdiv <- merged_df_div

library(ggrepel)

## Loading required package: ggplot2

library(factoextra)

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library(dplyr)

# Filtra il dataset in base ai valori di "replica_temp"
plotXsp_rep1 <- plotXsp[plotXenv$replica_temp == "1", ]
plotXenv_rep1 <- plotXenv[plotXenv$replica_temp == "1", ]
plotXsp_rep2 <- plotXsp[plotXenv$replica_temp == "2", ]
plotXenv_rep2 <- plotXenv[plotXenv$replica_temp == "2", ]

# Funzione per rimuovere colonne costanti o con variabilità nulla
remove_zero_variance <- function(df) {
  df <- df[, apply(df, 2, sd) != 0]
  return(df)
}

# Applicare la funzione ai tuoi dataset
plotXsp_rep1_clean <- remove_zero_variance(plotXsp_rep1)
plotXsp_rep2_clean <- remove_zero_variance(plotXsp_rep2)

# Esegui la PCA per il gruppo "1"
pca_result_rep1 <- prcomp(plotXsp_rep1_clean, scale. = TRUE)
var_rep1 <- get_pca_var(pca_result_rep1)

# Esegui la PCA per il gruppo "2"
pca_result_rep2 <- prcomp(plotXsp_rep2_clean, scale. = TRUE)
var_rep2 <- get_pca_var(pca_result_rep2)

# Crea il biplot PCA per il gruppo "1"
p_rep1 <- fviz_pca_biplot(pca_result_rep1, 
                          geom.ind = "point", 
                          col.ind = plotXenv_rep1$Treatment, 
                          addEllipses = TRUE, 
                          ellipse.level = 0.95, 
                          palette = "jco", 
                          legend.title = "Treatment",
                          label = "none")

# Aggiungi le etichette per le prime 20 variabili al gruppo "1"
p_rep1 <- p_rep1 + geom_text_repel(data = as.data.frame(var_rep1$coord[1:20, ]),
                                    aes(x = Dim.1, y = Dim.2, label = rownames(var_rep1$coord)[1:20]),
                                    size = 3, 
                                    max.overlaps = Inf,
                                    box.padding = 1.2, 
                                    point.padding = 0.7,
                                    segment.size = 0.2)

# Crea il biplot PCA per il gruppo "2"
p_rep2 <- fviz_pca_biplot(pca_result_rep2, 
                          geom.ind = "point", 
                          col.ind = plotXenv_rep2$Treatment, 
                          addEllipses = TRUE, 
                          ellipse.level = 0.95, 
                          palette = "jco", 
                          legend.title = "Treatment",
                          label = "none")

# Aggiungi le etichette per le prime 20 variabili al gruppo "2"
p_rep2 <- p_rep2 + geom_text_repel(data = as.data.frame(var_rep2$coord[1:20, ]),
                                    aes(x = Dim.1, y = Dim.2, label = rownames(var_rep2$coord)[1:20]),
                                    size = 3, 
                                    max.overlaps = Inf,
                                    box.padding = 1.2, 
                                    point.padding = 0.7,
                                    segment.size = 0.2)

# Visualizza i grafici
print(p_rep1)

print(p_rep2)

PCA per ogni caso separato

library(ggrepel)
library(factoextra)
library(dplyr)

# Funzione per rimuovere colonne costanti o con variabilità nulla
remove_zero_variance <- function(df) {
  df <- df[, apply(df, 2, sd) != 0]
  return(df)
}

# Ottieni i livelli unici della variabile "exp" nel dataset plotXenv
unique_exps <- unique(plotXenv$exp)

# Funzione per eseguire PCA e creare un biplot per un determinato livello di exp e replica_temp
perform_pca <- function(exp_level, replica_temp) {
  # Filtra i dataset in base al valore di "exp" e "replica_temp"
  plotXsp_filtered <- plotXsp[plotXenv$exp == exp_level & plotXenv$replica_temp == replica_temp, ]
  plotXenv_filtered <- plotXenv[plotXenv$exp == exp_level & plotXenv$replica_temp == replica_temp, ]
  
  # Rimuovi colonne con variabilità nulla
  plotXsp_filtered_clean <- remove_zero_variance(plotXsp_filtered)
  
  # Verifica se ci sono abbastanza righe dopo il filtraggio
  if (nrow(plotXsp_filtered_clean) < 2) {
    message(paste("Skipping PCA for", exp_level, "Replica", replica_temp, "- Not enough data"))
    return(NULL)
  }

  # Esegui la PCA
  pca_result <- prcomp(plotXsp_filtered_clean, scale. = TRUE)
  var_result <- get_pca_var(pca_result)
  
  # Determina il numero massimo di variabili da etichettare
  max_vars <- min(20, nrow(var_result$coord))
  
  # Crea il biplot PCA
  p <- fviz_pca_biplot(pca_result, 
                       geom.ind = "point", 
                       col.ind = plotXenv_filtered$Treatment, 
                       addEllipses = TRUE, 
                       ellipse.level = 0.95, 
                       palette = "jco", 
                       legend.title = "Treatment",
                       label = "none")
  
  # Aggiungi le etichette per le prime variabili disponibili
  p <- p + geom_text_repel(data = as.data.frame(var_result$coord[1:max_vars, ]),
                           aes(x = Dim.1, y = Dim.2, label = rownames(var_result$coord)[1:max_vars]),
                           size = 3, 
                           max.overlaps = Inf,
                           box.padding = 1.2, 
                           point.padding = 0.7,
                           segment.size = 0.2)
  
  # Aggiungi un titolo al grafico
  p <- p + ggtitle(paste("PCA for", exp_level, "Replica", replica_temp))
  
  return(p)
}

# Esegui PCA e crea grafici per ogni combinazione di "exp" e "replica_temp"
plots <- list()

for (exp_level in unique_exps) {
  for (replica_temp in c("1", "2")) {
    plot_name <- paste("PCA_", exp_level, "_replica_", replica_temp, sep = "")
    plots[[plot_name]] <- perform_pca(exp_level, replica_temp)
  }
}

# Visualizza e salva tutti i grafici non nulli
for (plot_name in names(plots)) {
  if (!is.null(plots[[plot_name]])) {
    print(plots[[plot_name]])
  }
}

## Warning: Computation failed in `stat_ellipse()`.
## Caused by error in `chol.default()`:
## ! the leading minor of order 1 is not positive

## Warning: Computation failed in `stat_ellipse()`.
## Caused by error in `chol.default()`:
## ! the leading minor of order 1 is not positive

## Warning: Computation failed in `stat_ellipse()`.
## Caused by error in `chol.default()`:
## ! the leading minor of order 1 is not positive

MATRICI DI DISSIMILARITA

MIN

library(vegan)

# Funzione per rimuovere righe con solo zeri
remove_zero_rows <- function(df) {
  df <- df[rowSums(df != 0) > 0, ]
  return(df)
}

# Applica la funzione al dataset
plotXsp_rep1_nozero <- remove_zero_rows(plotXsp_rep1_clean)

# Rimuovi righe con valori NA
plotXsp_rep1_nozero <- na.omit(plotXsp_rep1_nozero)

# Calcola la matrice di dissimilarità
dissimilarity_matrix <- vegdist(plotXsp_rep1_nozero, method = "bray")

# Esegui ANOSIM
anosim_result_rep1 <- anosim(dissimilarity_matrix, plotXenv[rownames(plotXsp_rep1_nozero),]$Treatment)
print(anosim_result_rep1)

## 
## Call:
## anosim(x = dissimilarity_matrix, grouping = plotXenv[rownames(plotXsp_rep1_nozero),      ]$Treatment) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.003854 
##       Significance: 0.199 
## 
## Permutation: free
## Number of permutations: 999

MAX

library(vegan)

# Funzione per rimuovere righe con solo zeri
remove_zero_rows <- function(df) {
  df <- df[rowSums(df != 0) > 0, ]
  return(df)
}

# Applica la funzione al dataset
plotXsp_rep2_nozero <- remove_zero_rows(plotXsp_rep2_clean)

# Rimuovi righe con valori NA
plotXsp_rep2_nozero <- na.omit(plotXsp_rep2_nozero)

# Calcola la matrice di dissimilarità
dissimilarity_matrix <- vegdist(plotXsp_rep2_nozero, method = "bray")

# Esegui ANOSIM
anosim_result_rep2 <- anosim(dissimilarity_matrix, plotXenv[rownames(plotXsp_rep2_nozero),]$Treatment)
print(anosim_result_rep2)

## 
## Call:
## anosim(x = dissimilarity_matrix, grouping = plotXenv[rownames(plotXsp_rep2_nozero),      ]$Treatment) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.02162 
##       Significance: 0.002 
## 
## Permutation: free
## Number of permutations: 999