ORPHEE : activités enzymatiques et diversité
Nicolas Fanin
2026-08-30
1 Présentation
Ce document présente l’analyse des activités enzymatiques des sols du dispositif ORPHEE. Les analyses portent sur les effets de la composition et de la richesse spécifique des communautés, ainsi que sur les relations entre les activités enzymatiques et les variables du sol.
Les figures sont générées directement par R lors du Knit du document.
Fichier de données requis :
2026_ORPHEE_Enzyme_CNP.csv, placé dans le même dossier que ce fichier.Rmd.
############################################################
# ANALYSE HFA / ACTIVITES ENZYMATIQUES
# ORPHEE
############################################################
############################################################
# 1. PACKAGES
############################################################
library(nlme)
library(dplyr)
library(tidyr)
library(ggplot2)
library(corrplot)
library(vegan)
############################################################
# 2. CHARGEMENT ET PREPARATION DES DONNEES
############################################################
# Le fichier CSV doit être placé dans le même dossier que ce Rmd.
ORPHEE_data <- read.csv(
"2026_ORPHEE_Enzyme_CNP.csv",
sep = ";",
header = TRUE
)
dim(ORPHEE_data)## [1] 256 31str(ORPHEE_data)## 'data.frame': 256 obs. of 31 variables:
## $ Block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Plot : int 1 2 3 4 1 2 3 4 5 6 ...
## $ SpNo : int 1 1 1 1 1 1 1 1 1 2 ...
## $ Bp : int 1 0 0 0 1 0 0 0 0 1 ...
## $ Pp : int 0 0 0 0 0 0 0 0 1 0 ...
## $ Qp : int 0 0 1 0 0 0 1 0 0 0 ...
## $ Qi : int 0 0 0 1 0 0 0 1 0 0 ...
## $ Qr : int 0 1 0 0 0 1 0 0 0 1 ...
## $ Depth : chr "A" "A" "A" "A" ...
## $ ID_Lab : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Ctot : num 31.4 26.3 27.8 28.9 32.1 ...
## $ Ntot : num 0.799 0.894 0.952 0.954 0.863 ...
## $ Ptot : num 0.003 0.004 0.004 0.004 0.003 0.003 0.003 0.003 0.004 0.004 ...
## $ pH : num 4.2 4.31 4.23 4.22 4.15 4.44 4.28 4.15 4.34 4.4 ...
## $ DM : int 91 89 88 90 92 92 90 87 90 91 ...
## $ BG : int 145 169 187 142 90 91 142 118 180 119 ...
## $ CBH : int 52 61 70 35 18 39 40 31 63 34 ...
## $ XYL : int 46 70 75 32 27 33 57 31 44 39 ...
## $ NAG : int 68 52 97 87 32 28 39 60 127 33 ...
## $ AP : int 208 260 254 266 106 136 230 186 191 176 ...
## $ AG : num 9 7.9 4.8 9.4 4.2 6.8 0.4 7.8 3 3.5 ...
## $ AS : num 5.7 5.3 5 6.1 2.8 6.1 0 6.2 0.1 5.9 ...
## $ LAP : num 2.3 2.5 2.4 2.1 1.8 2 1.8 1.7 1.9 2.7 ...
## $ inc..time: num 3 3 3 3 3 3 3 3 3 3 ...
## $ Cenz : num 0.75 0.92 1.01 0.66 0.42 0.51 0.72 0.56 0.87 0.59 ...
## $ Nenz : num 0.21 0.16 0.3 0.27 0.1 0.09 0.12 0.19 0.39 0.11 ...
## $ Penz : num 0.62 0.78 0.76 0.8 0.32 0.41 0.69 0.56 0.57 0.53 ...
## $ Senz : num 0.03 0.02 0.01 0.03 0.01 0.02 0 0.02 0.01 0.01 ...
## $ Cenz.Nenz: num 3.6 5.6 3.4 2.5 4.1 5.6 5.8 3 2.2 5.4 ...
## $ Nenz.Penz: num 0.3 0.2 0.4 0.3 0.3 0.2 0.2 0.3 0.7 0.2 ...
## $ Cenz.Penz: num 1.2 1.2 1.3 0.8 1.3 1.2 1 1 1.5 1.1 ...# Conversion des variables
ORPHEE_data <- ORPHEE_data %>%
mutate(
across(
c(Block, Plot, SpNo, Bp, Pp, Qp, Qi, Qr, Depth, ID_Lab),
factor
),
across(
c(DM, BG, CBH, XYL, NAG, AP),
as.numeric
)
)
str(ORPHEE_data)## 'data.frame': 256 obs. of 31 variables:
## $ Block : Factor w/ 4 levels "1","3","6","8": 1 1 1 1 1 1 1 1 1 1 ...
## $ Plot : Factor w/ 32 levels "1","2","3","4",..: 1 2 3 4 1 2 3 4 5 6 ...
## $ SpNo : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 2 ...
## $ Bp : Factor w/ 2 levels "0","1": 2 1 1 1 2 1 1 1 1 2 ...
## $ Pp : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 2 1 ...
## $ Qp : Factor w/ 2 levels "0","1": 1 1 2 1 1 1 2 1 1 1 ...
## $ Qi : Factor w/ 2 levels "0","1": 1 1 1 2 1 1 1 2 1 1 ...
## $ Qr : Factor w/ 2 levels "0","1": 1 2 1 1 1 2 1 1 1 2 ...
## $ Depth : Factor w/ 2 levels "A","B": 1 1 1 1 2 2 2 2 1 1 ...
## $ ID_Lab : Factor w/ 256 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Ctot : num 31.4 26.3 27.8 28.9 32.1 ...
## $ Ntot : num 0.799 0.894 0.952 0.954 0.863 ...
## $ Ptot : num 0.003 0.004 0.004 0.004 0.003 0.003 0.003 0.003 0.004 0.004 ...
## $ pH : num 4.2 4.31 4.23 4.22 4.15 4.44 4.28 4.15 4.34 4.4 ...
## $ DM : num 91 89 88 90 92 92 90 87 90 91 ...
## $ BG : num 145 169 187 142 90 91 142 118 180 119 ...
## $ CBH : num 52 61 70 35 18 39 40 31 63 34 ...
## $ XYL : num 46 70 75 32 27 33 57 31 44 39 ...
## $ NAG : num 68 52 97 87 32 28 39 60 127 33 ...
## $ AP : num 208 260 254 266 106 136 230 186 191 176 ...
## $ AG : num 9 7.9 4.8 9.4 4.2 6.8 0.4 7.8 3 3.5 ...
## $ AS : num 5.7 5.3 5 6.1 2.8 6.1 0 6.2 0.1 5.9 ...
## $ LAP : num 2.3 2.5 2.4 2.1 1.8 2 1.8 1.7 1.9 2.7 ...
## $ inc..time: num 3 3 3 3 3 3 3 3 3 3 ...
## $ Cenz : num 0.75 0.92 1.01 0.66 0.42 0.51 0.72 0.56 0.87 0.59 ...
## $ Nenz : num 0.21 0.16 0.3 0.27 0.1 0.09 0.12 0.19 0.39 0.11 ...
## $ Penz : num 0.62 0.78 0.76 0.8 0.32 0.41 0.69 0.56 0.57 0.53 ...
## $ Senz : num 0.03 0.02 0.01 0.03 0.01 0.02 0 0.02 0.01 0.01 ...
## $ Cenz.Nenz: num 3.6 5.6 3.4 2.5 4.1 5.6 5.8 3 2.2 5.4 ...
## $ Nenz.Penz: num 0.3 0.2 0.4 0.3 0.3 0.2 0.2 0.3 0.7 0.2 ...
## $ Cenz.Penz: num 1.2 1.2 1.3 0.8 1.3 1.2 1 1 1.5 1.1 ...############################################################
# 3. EXPLORATION DU DISPOSITIF
############################################################
# Nombre d'observations par niveau de diversité
table(ORPHEE_data$SpNo)##
## 1 2 3 4 5
## 40 80 80 40 16# Nombre d'observations par profondeur
table(ORPHEE_data$Depth)##
## A B
## 128 128# Diversité × profondeur
table(
ORPHEE_data$SpNo,
ORPHEE_data$Depth
)##
## A B
## 1 20 20
## 2 40 40
## 3 40 40
## 4 20 20
## 5 8 8# Nombre d'observations par bloc
table(ORPHEE_data$Block)##
## 1 3 6 8
## 64 64 64 64# Structure Block × SpNo × Depth
with(
ORPHEE_data,
table(Block, SpNo, Depth)
)## , , Depth = A
##
## SpNo
## Block 1 2 3 4 5
## 1 5 10 10 5 2
## 3 5 10 10 5 2
## 6 5 10 10 5 2
## 8 5 10 10 5 2
##
## , , Depth = B
##
## SpNo
## Block 1 2 3 4 5
## 1 5 10 10 5 2
## 3 5 10 10 5 2
## 6 5 10 10 5 2
## 8 5 10 10 5 2# Nombre de plots par traitement
ORPHEE_data %>%
count(Block, SpNo, Depth)## Block SpNo Depth n
## 1 1 1 A 5
## 2 1 1 B 5
## 3 1 2 A 10
## 4 1 2 B 10
## 5 1 3 A 10
## 6 1 3 B 10
## 7 1 4 A 5
## 8 1 4 B 5
## 9 1 5 A 2
## 10 1 5 B 2
## 11 3 1 A 5
## 12 3 1 B 5
## 13 3 2 A 10
## 14 3 2 B 10
## 15 3 3 A 10
## 16 3 3 B 10
## 17 3 4 A 5
## 18 3 4 B 5
## 19 3 5 A 2
## 20 3 5 B 2
## 21 6 1 A 5
## 22 6 1 B 5
## 23 6 2 A 10
## 24 6 2 B 10
## 25 6 3 A 10
## 26 6 3 B 10
## 27 6 4 A 5
## 28 6 4 B 5
## 29 6 5 A 2
## 30 6 5 B 2
## 31 8 1 A 5
## 32 8 1 B 5
## 33 8 2 A 10
## 34 8 2 B 10
## 35 8 3 A 10
## 36 8 3 B 10
## 37 8 4 A 5
## 38 8 4 B 5
## 39 8 5 A 2
## 40 8 5 B 2# Composition des traitements
ORPHEE_data %>%
distinct(
SpNo,
Bp,
Pp,
Qp,
Qi,
Qr
) %>%
arrange(SpNo)## SpNo Bp Pp Qp Qi Qr
## 1 1 1 0 0 0 0
## 2 1 0 0 0 0 1
## 3 1 0 0 1 0 0
## 4 1 0 0 0 1 0
## 5 1 0 1 0 0 0
## 6 2 1 0 0 0 1
## 7 2 1 0 1 0 0
## 8 2 1 0 0 1 0
## 9 2 1 1 0 0 0
## 10 2 0 0 1 0 1
## 11 2 0 0 0 1 1
## 12 2 0 1 0 0 1
## 13 2 0 0 1 1 0
## 14 2 0 1 1 0 0
## 15 2 0 1 0 1 0
## 16 3 1 0 1 0 1
## 17 3 1 0 0 1 1
## 18 3 1 1 0 0 1
## 19 3 1 0 1 1 0
## 20 3 1 1 1 0 0
## 21 3 1 1 0 1 0
## 22 3 0 0 1 1 1
## 23 3 0 1 1 0 1
## 24 3 0 1 0 1 1
## 25 3 0 1 1 1 0
## 26 4 1 0 1 1 1
## 27 4 1 1 1 0 1
## 28 4 1 1 0 1 1
## 29 4 1 1 1 1 0
## 30 4 0 1 1 1 1
## 31 5 1 1 1 1 1############################################################
# 4. VARIABLES ET STATISTIQUES DESCRIPTIVES
############################################################
# Enzymes
enzymes <- c(
"BG",
"CBH",
"XYL",
"AG",
"NAG",
"LAP",
"AP"
)
# Variables CNP du sol
env_vars <- c(
"Ctot",
"Ntot",
"Ptot",
"pH"
)
# Variables CNP enzymatiques
cnp_enz <- c(
"Cenz",
"Nenz",
"Penz"
)
# Statistiques CNP du sol par diversité et profondeur
ORPHEE_data %>%
group_by(
SpNo,
Depth
) %>%
summarise(
n = n(),
Ctot_mean = mean(Ctot, na.rm = TRUE),
Ctot_sd = sd(Ctot, na.rm = TRUE),
Ntot_mean = mean(Ntot, na.rm = TRUE),
Ntot_sd = sd(Ntot, na.rm = TRUE),
Ptot_mean = mean(Ptot, na.rm = TRUE),
Ptot_sd = sd(Ptot, na.rm = TRUE),
pH_mean = mean(pH, na.rm = TRUE),
pH_sd = sd(pH, na.rm = TRUE),
.groups = "drop"
)## # A tibble: 10 × 11
## SpNo Depth n Ctot_mean Ctot_sd Ntot_mean Ntot_sd Ptot_mean Ptot_sd
## <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 A 20 27.0 6.71 0.864 0.216 0.0035 0.000761
## 2 1 B 20 27.8 10.7 0.818 0.312 0.0028 0.00106
## 3 2 A 40 27.9 10.9 0.829 0.340 0.00322 0.00117
## 4 2 B 40 29.8 10.7 0.827 0.315 0.00325 0.00151
## 5 3 A 40 28.0 11.3 0.868 0.336 0.00338 0.00117
## 6 3 B 40 33.1 14.9 0.949 0.404 0.00358 0.00117
## 7 4 A 20 26.4 11.5 0.837 0.428 0.0029 0.000968
## 8 4 B 20 28.6 11.0 0.839 0.377 0.00345 0.00128
## 9 5 A 8 28.3 14.0 0.823 0.389 0.00325 0.00149
## 10 5 B 8 30.9 9.12 0.874 0.317 0.00338 0.000744
## # ℹ 2 more variables: pH_mean <dbl>, pH_sd <dbl># Statistiques des enzymes
ORPHEE_data %>%
group_by(
SpNo,
Depth
) %>%
summarise(
across(
all_of(enzymes),
list(
mean = ~mean(.x, na.rm = TRUE),
sd = ~sd(.x, na.rm = TRUE)
),
.names = "{.col}_{.fn}"
),
.groups = "drop"
)## # A tibble: 10 × 16
## SpNo Depth BG_mean BG_sd CBH_mean CBH_sd XYL_mean XYL_sd AG_mean AG_sd
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 A 110. 39.2 35.8 15.4 34.6 16.4 5.62 3.65
## 2 1 B 97.0 35.6 30.7 12.5 31.6 14.6 5.34 4.28
## 3 2 A 94.9 23.6 30.0 11.8 32.4 12.8 6.65 4.27
## 4 2 B 102. 29.5 32.9 13.8 38.3 15.9 6.79 3.98
## 5 3 A 107. 38.2 34.4 14.2 33.2 12.6 5.82 2.95
## 6 3 B 107. 50.8 35.1 18.0 39.0 20.3 5.75 4.48
## 7 4 A 73 15.4 22 6.40 21.8 9.74 5.96 4.94
## 8 4 B 75.2 27.5 21.8 11.8 25.6 14.2 5.78 6.70
## 9 5 A 77.6 19.7 24.1 9.19 25 11.8 4.04 5.24
## 10 5 B 78.8 24.9 22.2 13.3 28.5 19.3 4.22 7.78
## # ℹ 6 more variables: NAG_mean <dbl>, NAG_sd <dbl>, LAP_mean <dbl>,
## # LAP_sd <dbl>, AP_mean <dbl>, AP_sd <dbl>############################################################
# 5. CORRELATIONS ENTRE VARIABLES
############################################################
vars_corr <- c(
"Ctot",
"Ntot",
"Ptot",
"pH",
"Cenz",
"Nenz",
"Penz",
enzymes
)
cor_mat <- cor(
ORPHEE_data[, vars_corr],
use = "pairwise.complete.obs"
)
corrplot(
cor_mat,
method = "color",
type = "upper",
tl.col = "black",
tl.srt = 45
)
############################################################
# 6. CREATION DES TRAITEMENTS
############################################################
# Code numérique du traitement
ORPHEE_data <- ORPHEE_data %>%
mutate(
Treatment = paste(
Bp,
Pp,
Qp,
Qi,
Qr,
sep = "_"
),
Treatment_name = paste0(
ifelse(Bp == "1", "Bp", ""),
ifelse(Pp == "1", "Pp", ""),
ifelse(Qp == "1", "Qp", ""),
ifelse(Qi == "1", "Qi", ""),
ifelse(Qr == "1", "Qr", "")
)
)
# Vérification des traitements
ORPHEE_data %>%
distinct(
SpNo,
Treatment,
Treatment_name
) %>%
arrange(
SpNo,
Treatment_name
)## SpNo Treatment Treatment_name
## 1 1 1_0_0_0_0 Bp
## 2 1 0_1_0_0_0 Pp
## 3 1 0_0_0_1_0 Qi
## 4 1 0_0_1_0_0 Qp
## 5 1 0_0_0_0_1 Qr
## 6 2 1_1_0_0_0 BpPp
## 7 2 1_0_0_1_0 BpQi
## 8 2 1_0_1_0_0 BpQp
## 9 2 1_0_0_0_1 BpQr
## 10 2 0_1_0_1_0 PpQi
## 11 2 0_1_1_0_0 PpQp
## 12 2 0_1_0_0_1 PpQr
## 13 2 0_0_0_1_1 QiQr
## 14 2 0_0_1_1_0 QpQi
## 15 2 0_0_1_0_1 QpQr
## 16 3 1_1_0_1_0 BpPpQi
## 17 3 1_1_1_0_0 BpPpQp
## 18 3 1_1_0_0_1 BpPpQr
## 19 3 1_0_0_1_1 BpQiQr
## 20 3 1_0_1_1_0 BpQpQi
## 21 3 1_0_1_0_1 BpQpQr
## 22 3 0_1_0_1_1 PpQiQr
## 23 3 0_1_1_1_0 PpQpQi
## 24 3 0_1_1_0_1 PpQpQr
## 25 3 0_0_1_1_1 QpQiQr
## 26 4 1_1_0_1_1 BpPpQiQr
## 27 4 1_1_1_1_0 BpPpQpQi
## 28 4 1_1_1_0_1 BpPpQpQr
## 29 4 1_0_1_1_1 BpQpQiQr
## 30 4 0_1_1_1_1 PpQpQiQr
## 31 5 1_1_1_1_1 BpPpQpQiQr############################################################
# 7. EFFET DE LA COMPOSITION DES COMMUNAUTES
############################################################
# Cette analyse teste si les différents assemblages
# de communautés diffèrent pour Cenz, Nenz et Penz.
#
# Effets fixes :
# Treatment_name
# Depth
# interaction Treatment × Depth
#
# Effet aléatoire :
# Block
mod_Cenz_treatment <- lme(
Cenz ~ Treatment_name * Depth,
random = ~1 | Block,
data = ORPHEE_data,
method = "REML"
)
mod_Nenz_treatment <- lme(
Nenz ~ Treatment_name * Depth,
random = ~1 | Block,
data = ORPHEE_data,
method = "REML"
)
mod_Penz_treatment <- lme(
Penz ~ Treatment_name * Depth,
random = ~1 | Block,
data = ORPHEE_data,
method = "REML"
)
# Résultats
anova(mod_Cenz_treatment)## numDF denDF F-value p-value
## (Intercept) 1 191 176.31037 <.0001
## Treatment_name 30 191 2.27177 0.0005
## Depth 1 191 0.43624 0.5097
## Treatment_name:Depth 30 191 0.51934 0.9823anova(mod_Nenz_treatment)## numDF denDF F-value p-value
## (Intercept) 1 191 686.4263 <.0001
## Treatment_name 30 191 2.1423 0.0011
## Depth 1 191 22.9368 <.0001
## Treatment_name:Depth 30 191 0.3755 0.9989anova(mod_Penz_treatment)## numDF denDF F-value p-value
## (Intercept) 1 191 158.22972 <.0001
## Treatment_name 30 191 2.84152 <.0001
## Depth 1 191 10.95087 0.0011
## Treatment_name:Depth 30 191 0.14009 1.0000############################################################
# 8. VISUALISATION Cenz, Nenz ET Penz
############################################################
ORPHEE_long <- ORPHEE_data %>%
select(
Treatment_name,
Depth,
Cenz,
Nenz,
Penz
) %>%
pivot_longer(
cols = c(
Cenz,
Nenz,
Penz
),
names_to = "Element",
values_to = "Activity"
)
ggplot(
ORPHEE_long,
aes(
x = Treatment_name,
y = Activity,
fill = Depth
)
) +
geom_boxplot(
position = position_dodge(width = 0.8),
width = 0.7
) +
scale_fill_manual(
values = c(
"A" = "#6BAED6",
"B" = "#F4A261"
)
) +
facet_grid(
Element ~ .,
scales = "free_y",
switch = "y"
) +
labs(
x = "Species composition",
y = expression(
paste(
"Enzyme activity (mmol ",
kg^{-1},
")"
)
),
fill = "Depth"
) +
theme_classic() +
theme(
axis.text.x = element_text(
angle = 90,
hjust = 1,
vjust = 0.5,
colour = "black"
),
axis.text.y = element_text(
colour = "black"
),
axis.title = element_text(
colour = "black"
),
axis.line = element_line(
colour = "black",
linewidth = 0.8
),
panel.border = element_rect(
colour = "black",
fill = NA,
linewidth = 0.8
),
strip.placement = "outside",
strip.background = element_rect(
colour = "black",
fill = "white",
linewidth = 0.8
),
strip.text.y.left = element_text(
angle = 90,
colour = "black"
)
)
############################################################
# 9. EFFET DE LA DIVERSITE
############################################################
# On considère ici SpNo comme un facteur.
#
# Question :
# l'activité change-t-elle avec le nombre d'espèces ?
#
# L'interaction avec Depth permet de vérifier si cette
# réponse diffère entre les deux profondeurs.
mod_Cenz_diversity <- lme(
Cenz ~ SpNo * Depth,
random = ~1 | Block,
data = ORPHEE_data,
method = "REML"
)
mod_Nenz_diversity <- lme(
Nenz ~ SpNo * Depth,
random = ~1 | Block,
data = ORPHEE_data,
method = "REML"
)
mod_Penz_diversity <- lme(
Penz ~ SpNo * Depth,
random = ~1 | Block,
data = ORPHEE_data,
method = "REML"
)
# ANOVA
anova(mod_Cenz_diversity)## numDF denDF F-value p-value
## (Intercept) 1 243 176.31037 <.0001
## SpNo 4 243 8.67985 <.0001
## Depth 1 243 0.44555 0.5051
## SpNo:Depth 4 243 0.71914 0.5796anova(mod_Nenz_diversity)## numDF denDF F-value p-value
## (Intercept) 1 243 686.4263 <.0001
## SpNo 4 243 4.7049 0.0011
## Depth 1 243 22.6774 <.0001
## SpNo:Depth 4 243 0.4246 0.7908anova(mod_Penz_diversity)## numDF denDF F-value p-value
## (Intercept) 1 243 158.22972 <.0001
## SpNo 4 243 4.89779 0.0008
## Depth 1 243 10.31091 0.0015
## SpNo:Depth 4 243 0.36705 0.8320############################################################
# 10. VISUALISATION DE L'EFFET DE LA DIVERSITE
############################################################
ORPHEE_diversity <- ORPHEE_data %>%
select(
SpNo,
Depth,
Cenz,
Nenz,
Penz
) %>%
pivot_longer(
cols = c(
Cenz,
Nenz,
Penz
),
names_to = "Element",
values_to = "Activity"
)
ggplot(
ORPHEE_diversity,
aes(
x = factor(SpNo),
y = Activity,
fill = Depth
)
) +
geom_boxplot(
position = position_dodge(width = 0.8),
width = 0.7
) +
scale_fill_manual(
values = c(
"A" = "#6BAED6",
"B" = "#F4A261"
)
) +
facet_grid(
Element ~ .,
scales = "free_y"
) +
labs(
x = "Number of species",
y = expression(
paste(
"Enzyme activity (mmol ",
kg^{-1},
")"
)
),
fill = "Depth"
) +
theme_classic() +
theme(
axis.text = element_text(
colour = "black"
),
axis.title = element_text(
colour = "black"
),
axis.line = element_line(
colour = "black",
linewidth = 0.8
),
panel.border = element_rect(
colour = "black",
fill = NA,
linewidth = 0.8
),
strip.background = element_rect(
colour = "black",
fill = "white",
linewidth = 0.8
),
strip.text = element_text(
colour = "black"
)
)
############################################################
# 11. PIN COMME REFERENCE
############################################################
# Ici on change de question.
#
# Le Pin (Pp) est pris comme référence.
#
# La monoculture de Pin correspond à :
#
# Pp = 1
# Bp = 0
# Qp = 0
# Qi = 0
# Qr = 0
#
# Pour chaque assemblage contenant Pp, on calcule :
#
# activité de l'assemblage
# -
# activité de la monoculture de Pin
#
# Une valeur négative signifie donc une activité plus faible
# que celle de la monoculture de Pin.
############################################################
# 11.1 Monoculture de Pin
############################################################
pin_mono <- ORPHEE_data %>%
filter(
Pp == "1",
Bp == "0",
Qp == "0",
Qi == "0",
Qr == "0"
)
############################################################
# 11.2 Valeurs de référence du Pin
############################################################
pin_reference <- pin_mono %>%
group_by(
Depth
) %>%
summarise(
Cenz_ref = mean(
Cenz,
na.rm = TRUE
),
Nenz_ref = mean(
Nenz,
na.rm = TRUE
),
Penz_ref = mean(
Penz,
na.rm = TRUE
),
.groups = "drop"
)
############################################################
# 11.3 Assemblages contenant le Pin
############################################################
pin_assemblages <- ORPHEE_data %>%
filter(
Pp == "1"
) %>%
group_by(
Treatment_name,
SpNo,
Depth
) %>%
summarise(
Cenz_mean = mean(
Cenz,
na.rm = TRUE
),
Nenz_mean = mean(
Nenz,
na.rm = TRUE
),
Penz_mean = mean(
Penz,
na.rm = TRUE
),
Cenz_sd = sd(
Cenz,
na.rm = TRUE
),
Nenz_sd = sd(
Nenz,
na.rm = TRUE
),
Penz_sd = sd(
Penz,
na.rm = TRUE
),
n_Cenz = sum(
!is.na(Cenz)
),
n_Nenz = sum(
!is.na(Nenz)
),
n_Penz = sum(
!is.na(Penz)
),
.groups = "drop"
)
############################################################
# 11.4 Difference par rapport au Pin monoculture
############################################################
pin_effects <- pin_assemblages %>%
left_join(
pin_reference,
by = "Depth"
) %>%
mutate(
Cenz_effect = Cenz_mean - Cenz_ref,
Nenz_effect = Nenz_mean - Nenz_ref,
Penz_effect = Penz_mean - Penz_ref
)
############################################################
# 11.5 Format long
############################################################
pin_effects_long <- pin_effects %>%
select(
Treatment_name,
SpNo,
Depth,
Cenz_effect,
Nenz_effect,
Penz_effect
) %>%
pivot_longer(
cols = c(
Cenz_effect,
Nenz_effect,
Penz_effect
),
names_to = "Element",
values_to = "Effect"
) %>%
mutate(
Element = case_when(
Element == "Cenz_effect" ~ "Cenz",
Element == "Nenz_effect" ~ "Nenz",
Element == "Penz_effect" ~ "Penz"
)
)
############################################################
# 12. EFFET DE LA DIVERSIFICATION RELATIF AU PIN MONOCULTURE
############################################################
library(dplyr)
library(tidyr)
library(ggplot2)
############################################################
# 12.1 IDENTIFICATION DU PIN MONOCULTURE
############################################################
# Le Pin monoculture correspond à :
# Pp = 1
# Bp = 0
# Qp = 0
# Qi = 0
# Qr = 0
pin_data <- ORPHEE_data %>%
filter(
as.character(Pp) == "1",
as.character(Bp) == "0",
as.character(Qp) == "0",
as.character(Qi) == "0",
as.character(Qr) == "0"
)
############################################################
# 12.2 VÉRIFICATION DU PIN MONOCULTURE
############################################################
cat(
"Nombre d'observations du Pin monoculture :",
nrow(pin_data),
"\n"
)## Nombre d'observations du Pin monoculture : 8cat("\nRépartition du Pin monoculture par profondeur :\n")##
## Répartition du Pin monoculture par profondeur :print(
table(
pin_data$Depth,
useNA = "ifany"
)
)##
## A B
## 4 4############################################################
# 12.3 MOYENNE DU PIN MONOCULTURE PAR PROFONDEUR
############################################################
# La référence est calculée séparément pour chaque profondeur.
#
# Cela permet de comparer chaque communauté à la valeur
# moyenne du Pin monoculture dans la même profondeur.
pin_reference <- pin_data %>%
group_by(Depth) %>%
summarise(
Cenz_Pin = mean(
Cenz,
na.rm = TRUE
),
Nenz_Pin = mean(
Nenz,
na.rm = TRUE
),
Penz_Pin = mean(
Penz,
na.rm = TRUE
),
.groups = "drop"
)
cat("\nValeurs de référence du Pin monoculture :\n")##
## Valeurs de référence du Pin monoculture :print(pin_reference)## # A tibble: 2 × 4
## Depth Cenz_Pin Nenz_Pin Penz_Pin
## <fct> <dbl> <dbl> <dbl>
## 1 A 0.548 0.182 0.512
## 2 B 0.478 0.0975 0.412############################################################
# 12.4 CALCUL DE L'EFFET RELATIF AU PIN
############################################################
# Pour chaque observation :
#
# Effet = valeur de la communauté - valeur moyenne du Pin
#
# Une valeur :
# > 0 = activité supérieure au Pin
# < 0 = activité inférieure au Pin
# = 0 = activité équivalente au Pin
pin_effects <- ORPHEE_data %>%
# On retire la monoculture de Pin elle-même
filter(
!(
as.character(Pp) == "1" &
as.character(Bp) == "0" &
as.character(Qp) == "0" &
as.character(Qi) == "0" &
as.character(Qr) == "0"
)
) %>%
# Ajout de la référence correspondant à la profondeur
left_join(
pin_reference,
by = "Depth"
) %>%
mutate(
Cenz_effect = Cenz - Cenz_Pin,
Nenz_effect = Nenz - Nenz_Pin,
Penz_effect = Penz - Penz_Pin
)
############################################################
# 12.5 VÉRIFICATION DE LA RICHESSE SPÉCIFIQUE
############################################################
cat("\nRépartition de SpNo avant nettoyage :\n")##
## Répartition de SpNo avant nettoyage :print(
table(
pin_effects$SpNo,
useNA = "ifany"
)
)##
## 1 2 3 4 5
## 32 80 80 40 16############################################################
# 12.6 PASSAGE AU FORMAT LONG
############################################################
# On conserve uniquement les communautés de richesse 2 à 5.
#
# La richesse 1 n'est pas utilisée car elle correspond
# au Pin monoculture de référence.
#
# Les éventuelles valeurs NA de SpNo sont supprimées.
pin_effects_long <- pin_effects %>%
filter(
!is.na(SpNo),
as.character(SpNo) %in% c("2", "3", "4", "5")
) %>%
select(
SpNo,
Depth,
Cenz_effect,
Nenz_effect,
Penz_effect
) %>%
pivot_longer(
cols = c(
Cenz_effect,
Nenz_effect,
Penz_effect
),
names_to = "Element",
values_to = "Effect"
) %>%
mutate(
Element = recode(
Element,
Cenz_effect = "Cenz",
Nenz_effect = "Nenz",
Penz_effect = "Penz"
),
SpNo = factor(
as.character(SpNo),
levels = c("2", "3", "4", "5")
)
)
############################################################
# 12.7 VÉRIFICATION DES DONNÉES FINALES
############################################################
cat(
"\nNombre de lignes finales :",
nrow(pin_effects_long),
"\n"
)##
## Nombre de lignes finales : 648cat("\nNombre d'observations par richesse :\n")##
## Nombre d'observations par richesse :print(
table(
pin_effects_long$SpNo,
useNA = "ifany"
)
)##
## 2 3 4 5
## 240 240 120 48cat("\nNombre d'observations par profondeur :\n")##
## Nombre d'observations par profondeur :print(
table(
pin_effects_long$Depth,
useNA = "ifany"
)
)##
## A B
## 324 324cat("\nNombre d'observations par élément :\n")##
## Nombre d'observations par élément :print(
table(
pin_effects_long$Element,
useNA = "ifany"
)
)##
## Cenz Nenz Penz
## 216 216 216############################################################
# 12.8 CALCUL DES MOYENNES ET DES ERREURS STANDARD
############################################################
# Pour chaque combinaison :
# richesse × profondeur × élément
#
# Mean = moyenne de l'effet relatif au Pin
#
# SD = écart-type
#
# SE = SD / sqrt(n)
pin_effects_summary <- pin_effects_long %>%
group_by(
SpNo,
Depth,
Element
) %>%
summarise(
n = sum(
!is.na(Effect)
),
Mean = mean(
Effect,
na.rm = TRUE
),
SD = sd(
Effect,
na.rm = TRUE
),
SE = SD / sqrt(n),
.groups = "drop"
)
############################################################
# 12.9 AFFICHER LE TABLEAU RÉCAPITULATIF
############################################################
print(
pin_effects_summary
)## # A tibble: 24 × 7
## SpNo Depth Element n Mean SD SE
## <fct> <fct> <chr> <int> <dbl> <dbl> <dbl>
## 1 2 A Cenz 40 -0.0555 0.142 0.0224
## 2 2 A Nenz 40 -0.031 0.0784 0.0124
## 3 2 A Penz 40 -0.0485 0.143 0.0226
## 4 2 B Cenz 40 0.0637 0.172 0.0271
## 5 2 B Nenz 40 0.0165 0.0385 0.00608
## 6 2 B Penz 40 0.0160 0.163 0.0257
## 7 3 A Cenz 40 -0.00525 0.188 0.0297
## 8 3 A Nenz 40 -0.0155 0.0678 0.0107
## 9 3 A Penz 40 0.0228 0.236 0.0373
## 10 3 B Cenz 40 0.0820 0.264 0.0418
## # ℹ 14 more rows############################################################
# 12.10 FIGURE : EFFET RELATIF AU PIN MONOCULTURE
############################################################
ggplot(
pin_effects_summary,
aes(
x = SpNo,
y = Mean,
fill = Depth
)
) +
# Référence = Pin monoculture
geom_hline(
yintercept = 0,
linetype = "dashed",
colour = "black",
linewidth = 0.6
) +
# Moyenne
geom_col(
position = position_dodge(
width = 0.8
),
width = 0.65
) +
# Erreur standard
geom_errorbar(
aes(
ymin = Mean - SE,
ymax = Mean + SE
),
position = position_dodge(
width = 0.8
),
width = 0.15,
linewidth = 0.6
) +
# Couleurs des profondeurs
scale_fill_manual(
values = c(
"A" = "#6BAED6",
"B" = "#F4A261"
)
) +
# Un panneau par élément
facet_grid(
Element ~ .,
scales = "free_y"
) +
labs(
x = "Species richness",
y = "Effect relative to Pin monoculture\n(mmol kg⁻¹)",
fill = "Depth"
) +
theme_classic() +
theme(
axis.text = element_text(
colour = "black"
),
axis.title = element_text(
colour = "black"
),
axis.line = element_line(
colour = "black",
linewidth = 0.8
),
panel.border = element_rect(
colour = "black",
fill = NA,
linewidth = 0.8
),
strip.background = element_rect(
colour = "black",
fill = "white",
linewidth = 0.8
),
strip.text = element_text(
colour = "black"
)
)
######################################################
### 13. RDA : enzymes ~ CNP + pH du sol
######################################################
library(vegan)
library(dplyr)
######################################################
### 13.1 Variables
######################################################
# Variables enzymatiques = variables réponses
enzymes <- c(
"BG",
"CBH",
"XYL",
"AG",
"NAG",
"LAP",
"AP"
)
# Variables environnementales = construisent la RDA
env_vars <- c(
"Ctot",
"Ntot",
"Ptot",
"pH"
)
# Variables supplémentaires
# Elles ne participent pas à la construction des axes
supp_vars <- c(
"Cenz",
"Nenz",
"Penz"
)
######################################################
### 13.2 Préparation des données
######################################################
RDA_data <- ORPHEE_data %>%
select(
Block,
Depth,
SpNo,
Treatment_name,
all_of(
c(
enzymes,
env_vars,
supp_vars
)
)
) %>%
na.omit()
######################################################
### 13.3 Matrice des enzymes
######################################################
Y <- RDA_data %>%
select(
all_of(enzymes)
)
######################################################
### 13.4 RDA
######################################################
# Les enzymes sont expliquées par :
# Ctot, Ntot, Ptot et pH.
#
# scale = TRUE :
# les enzymes sont standardisées avant l'analyse.
mod_RDA <- rda(
Y ~ Ctot + Ntot + Ptot + pH,
data = RDA_data,
scale = TRUE
)
######################################################
### 13.5 Résumé de la RDA
######################################################
summary(mod_RDA)##
## Call:
## rda(formula = Y ~ Ctot + Ntot + Ptot + pH, data = RDA_data, scale = TRUE)
##
## Partitioning of correlations:
## Inertia Proportion
## Total 7.000 1.0000
## Constrained 1.495 0.2136
## Unconstrained 5.505 0.7864
##
## Eigenvalues, and their contribution to the correlations
##
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 PC1 PC2 PC3
## Eigenvalue 1.0032 0.44422 0.03913 0.008485 2.8939 0.9242 0.7591
## Proportion Explained 0.1433 0.06346 0.00559 0.001212 0.4134 0.1320 0.1084
## Cumulative Proportion 0.1433 0.20677 0.21236 0.213574 0.6270 0.7590 0.8675
## PC4 PC5 PC6 PC7
## Eigenvalue 0.50393 0.27985 0.09071 0.053347
## Proportion Explained 0.07199 0.03998 0.01296 0.007621
## Cumulative Proportion 0.93944 0.97942 0.99238 1.000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## RDA1 RDA2 RDA3 RDA4
## Eigenvalue 1.003 0.4442 0.03913 0.008485
## Proportion Explained 0.671 0.2971 0.02618 0.005675
## Cumulative Proportion 0.671 0.9681 0.99432 1.000000######################################################
### 13.6 Significativité globale
######################################################
anova(mod_RDA)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = Y ~ Ctot + Ntot + Ptot + pH, data = RDA_data, scale = TRUE)
## Df Variance F Pr(>F)
## Model 4 1.495 17.041 0.001 ***
## Residual 251 5.505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1######################################################
### 13.7 Significativité des axes
######################################################
anova(
mod_RDA,
by = "axis"
)## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = Y ~ Ctot + Ntot + Ptot + pH, data = RDA_data, scale = TRUE)
## Df Variance F Pr(>F)
## RDA1 1 1.0032 45.7401 0.001 ***
## RDA2 1 0.4442 20.3349 0.001 ***
## RDA3 1 0.0391 1.7985 0.319
## RDA4 1 0.0085 0.3915 0.802
## Residual 251 5.5050
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1######################################################
### 13.8 Significativité des variables environnementales
######################################################
anova(
mod_RDA,
by = "terms"
)## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = Y ~ Ctot + Ntot + Ptot + pH, data = RDA_data, scale = TRUE)
## Df Variance F Pr(>F)
## Ctot 1 0.6741 30.7360 0.001 ***
## Ntot 1 0.5542 25.2678 0.001 ***
## Ptot 1 0.1330 6.0632 0.005 **
## pH 1 0.1338 6.0984 0.002 **
## Residual 251 5.5050
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1######################################################
### 13.9 R2 et R2 ajusté
######################################################
RsquareAdj(mod_RDA)## $r.squared
## [1] 0.2135741
##
## $adj.r.squared
## [1] 0.2010414######################################################
### 13.10 Variables supplémentaires
######################################################
# Cenz, Nenz et Penz sont projetées dans l'espace
# de la RDA mais ne participent pas à la construction
# des axes.
supp <- RDA_data %>%
select(
all_of(supp_vars)
)
fit_supp <- envfit(
mod_RDA,
supp,
permutations = 999
)
# Résultats des variables supplémentaires
fit_supp##
## ***VECTORS
##
## RDA1 RDA2 r2 Pr(>r)
## Cenz 0.84437 -0.53576 0.9191 0.001 ***
## Nenz 0.99950 -0.03154 0.4382 0.001 ***
## Penz 0.89493 0.44620 0.8450 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999######################################################
### 13.11 Coordonnées des observations
######################################################
site_scores <- scores(
mod_RDA,
display = "sites",
scaling = 2
)
######################################################
### 13.12 Coordonnées des enzymes
######################################################
enzyme_scores <- scores(
mod_RDA,
display = "species",
scaling = 2
)
######################################################
### 13.13 Coordonnées des variables environnementales
######################################################
env_scores <- scores(
mod_RDA,
display = "bp",
scaling = 2
)
######################################################
### 13.14 Coordonnées des variables supplémentaires
######################################################
supp_scores <- scores(
fit_supp,
display = "vectors"
)
######################################################
### 13.15 Pourcentage de variance expliquée
######################################################
eig <- summary(mod_RDA)$cont$importance
RDA1_percent <- round(
eig[2, 1] * 100,
1
)
RDA2_percent <- round(
eig[2, 2] * 100,
1
)
######################################################
### 13.16 Graphique RDA
######################################################
plot(
mod_RDA,
scaling = 2,
type = "n",
xlab = paste0(
"RDA1 (",
RDA1_percent,
"%)"
),
ylab = paste0(
"RDA2 (",
RDA2_percent,
"%)"
),
main = "RDA : activités enzymatiques et environnement du sol"
)
######################################################
### 13.17 Points : profondeur + richesse spécifique
######################################################
# Formes selon le nombre d'espèces
#
# 1 espèce = cercle
# 2 espèces = carré
# 3 espèces = losange
# 4 espèces = triangle
# 5 espèces = triangle inversé
pch_species <- c(
"1" = 21,
"2" = 22,
"3" = 23,
"4" = 24,
"5" = 25
)
# Couleur de remplissage selon la profondeur
point_col <- ifelse(
RDA_data$Depth == "A",
"#6BAED6",
"#F4A261"
)
# Ajout des points
points(
site_scores[, 1],
site_scores[, 2],
pch = pch_species[
as.character(RDA_data$SpNo)
],
bg = point_col,
col = "black",
cex = 0.9
)
######################################################
### 13.18 Variables environnementales
######################################################
arrows(
0,
0,
env_scores[, 1],
env_scores[, 2],
length = 0.08,
col = "blue",
lwd = 1.8
)
text(
env_scores[, 1],
env_scores[, 2],
labels = rownames(env_scores),
col = "blue",
pos = 3,
cex = 1.1
)
######################################################
### 13.19 Enzymes
######################################################
arrows(
0,
0,
enzyme_scores[, 1],
enzyme_scores[, 2],
length = 0.08,
col = "black",
lwd = 1.5
)
text(
enzyme_scores[, 1],
enzyme_scores[, 2],
labels = rownames(enzyme_scores),
col = "black",
pos = 3,
cex = 1
)
######################################################
### 13.20 Variables supplémentaires
######################################################
arrows(
0,
0,
supp_scores[, 1],
supp_scores[, 2],
length = 0.08,
col = "red",
lwd = 1.5,
lty = 2
)
text(
supp_scores[, 1],
supp_scores[, 2],
labels = rownames(supp_scores),
col = "red",
pos = 3,
cex = 1
)
######################################################
### 13.21 Légende profondeur
######################################################
legend(
"topright",
legend = c(
"Depth A",
"Depth B",
"CNP + pH du sol",
"Enzymes",
"Cenz / Nenz / Penz"
),
pch = c(
21,
21,
NA,
NA,
NA
),
pt.bg = c(
"#6BAED6",
"#F4A261",
NA,
NA,
NA
),
col = c(
"black",
"black",
"blue",
"black",
"red"
),
lty = c(
NA,
NA,
1,
1,
2
),
lwd = c(
NA,
NA,
1.8,
1.5,
1.5
),
bty = "n"
)
######################################################
### 13.22 Légende richesse spécifique
######################################################
legend(
"bottomright",
legend = c(
"1 species",
"2 species",
"3 species",
"4 species",
"5 species"
),
pch = c(
21,
22,
23,
24,
25
),
pt.bg = "white",
col = "black",
bty = "n",
title = "Species richness"
)