Load library
library(vegan)
Load test OTU table and evironmental data with replicates
nifhotu_rep <- read.delim("~/Desktop/nifhotu_with_reps.txt", row.names=1)
env_rep <- read.delim("~/Desktop/env_with_reps.txt", row.names=1)
Load test OTU table and evironmental data without replicates
nifhotu_norep <- read.delim("~/Desktop/nifhotu_no_reps.txt", row.names=1)
env_norep <- read.delim("~/Desktop/env_no_reps.txt", row.names=1)
RDA with data keeping replicates
#distance calculation
bray_dist_rep <- vegdist(nifhotu_rep,method = "bray")
#rda analysis
otu_dbrda_rep <- dbrda(sqrt(bray_dist_rep)~.,data=env_rep,add = TRUE)
anova(otu_dbrda_rep,permutations = how(nperm = 999))
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = sqrt(bray_dist_rep) ~ pH + tn + tp + tk + ep + ek + emo + cec + den + nl + fl + sl + wc, data = env_rep, add = TRUE)
## Df SumOfSqs F Pr(>F)
## Model 12 8.9011 1.9271 0.001 ***
## Residual 63 24.2492
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#forward selection
step_forward_rep <- ordistep(dbrda(sqrt(bray_dist_rep)~1,data=env_rep,add = TRUE),scope = formula(otu_dbrda_rep),direction="forward")
##
## Start: sqrt(bray_dist_rep) ~ 1
##
## Df AIC F Pr(>F)
## + pH 1 264.42 4.6703 0.005 **
## + tp 1 266.44 2.6099 0.005 **
## + tn 1 267.02 2.0293 0.005 **
## + nl 1 267.09 1.9555 0.005 **
## + emo 1 267.23 1.8201 0.005 **
## + fl 1 267.33 1.7211 0.005 **
## + den 1 267.46 1.5851 0.005 **
## + tk 1 267.47 1.5788 0.005 **
## + cec 1 267.51 1.5349 0.005 **
## + sl 1 267.58 1.4695 0.005 **
## + ep 1 267.65 1.3958 0.005 **
## + wc 1 267.53 1.5162 0.010 **
## + ek 1 267.88 1.1758 0.070 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH
##
## Df AIC F Pr(>F)
## + tn 1 264.30 2.0697 0.005 **
## + emo 1 264.64 1.7356 0.005 **
## + den 1 264.68 1.6936 0.005 **
## + tp 1 264.72 1.6565 0.005 **
## + wc 1 264.76 1.6177 0.005 **
## + fl 1 264.76 1.6168 0.005 **
## + tk 1 264.76 1.6140 0.005 **
## + cec 1 264.80 1.5721 0.005 **
## + sl 1 264.82 1.5542 0.005 **
## + nl 1 264.84 1.5320 0.005 **
## + ep 1 264.99 1.3925 0.010 **
## + ek 1 265.18 1.2016 0.030 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn
##
## Df AIC F Pr(>F)
## + emo 1 264.46 1.7612 0.005 **
## + den 1 264.52 1.7025 0.005 **
## + fl 1 264.63 1.5938 0.005 **
## + tk 1 264.64 1.5894 0.005 **
## + sl 1 264.64 1.5886 0.005 **
## + cec 1 264.65 1.5806 0.005 **
## + wc 1 264.65 1.5798 0.005 **
## + tp 1 264.69 1.5388 0.005 **
## + nl 1 264.79 1.4428 0.005 **
## + ep 1 264.81 1.4273 0.005 **
## + ek 1 265.01 1.2264 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn + emo
##
## Df AIC F Pr(>F)
## + den 1 264.67 1.6899 0.005 **
## + tk 1 264.75 1.6178 0.005 **
## + fl 1 264.80 1.5673 0.005 **
## + cec 1 264.80 1.5666 0.005 **
## + sl 1 264.81 1.5581 0.005 **
## + tp 1 264.81 1.5545 0.005 **
## + wc 1 264.82 1.5467 0.005 **
## + ep 1 264.93 1.4471 0.005 **
## + nl 1 264.92 1.4515 0.010 **
## + ek 1 265.17 1.2181 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn + emo + den
##
## Df AIC F Pr(>F)
## + sl 1 264.92 1.6384 0.005 **
## + fl 1 264.92 1.6339 0.005 **
## + tp 1 264.98 1.5733 0.005 **
## + cec 1 265.01 1.5518 0.005 **
## + tk 1 265.01 1.5512 0.005 **
## + nl 1 265.09 1.4754 0.005 **
## + ep 1 265.10 1.4636 0.005 **
## + wc 1 265.11 1.4558 0.005 **
## + ek 1 265.34 1.2361 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn + emo + den + sl
##
## Df AIC F Pr(>F)
## + cec 1 265.17 1.6037 0.005 **
## + tp 1 265.19 1.5856 0.005 **
## + tk 1 265.23 1.5514 0.005 **
## + nl 1 265.29 1.4937 0.005 **
## + fl 1 265.29 1.4937 0.005 **
## + ep 1 265.31 1.4732 0.005 **
## + wc 1 265.31 1.4687 0.005 **
## + ek 1 265.56 1.2453 0.015 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn + emo + den + sl + cec
##
## Df AIC F Pr(>F)
## + tk 1 265.43 1.5708 0.005 **
## + nl 1 265.50 1.5069 0.005 **
## + fl 1 265.50 1.5069 0.005 **
## + wc 1 265.55 1.4680 0.005 **
## + tp 1 265.55 1.4634 0.005 **
## + ep 1 265.56 1.4582 0.005 **
## + ek 1 265.61 1.4074 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn + emo + den + sl + cec + tk
##
## Df AIC F Pr(>F)
## + nl 1 265.73 1.5160 0.005 **
## + fl 1 265.73 1.5160 0.005 **
## + ep 1 265.78 1.4734 0.005 **
## + wc 1 265.79 1.4672 0.005 **
## + tp 1 265.82 1.4388 0.005 **
## + ek 1 265.84 1.4185 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn + emo + den + sl + cec + tk + nl
##
## Df AIC F Pr(>F)
## + tp 1 266.09 1.4401 0.005 **
## + ek 1 266.11 1.4222 0.005 **
## + ep 1 266.20 1.3411 0.010 **
## + wc 1 266.28 1.2770 0.020 *
## + fl 0 265.73
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn + emo + den + sl + cec + tk + nl + tp
##
## Df AIC F Pr(>F)
## + ek 1 266.46 1.4102 0.005 **
## + ep 1 266.52 1.3535 0.010 **
## + wc 1 266.60 1.2848 0.015 *
## + fl 0 266.09
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn + emo + den + sl + cec + tk + nl + tp + ek
##
## Df AIC F Pr(>F)
## + ep 1 266.84 1.3825 0.005 **
## + wc 1 266.97 1.2697 0.015 *
## + fl 0 266.46
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn + emo + den + sl + cec + tk + nl + tp + ek + ep
##
## Df AIC F Pr(>F)
## + wc 1 267.31 1.2775 0.01 **
## + fl 0 266.84
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_rep) ~ pH + tn + emo + den + sl + cec + tk + nl + tp + ek + ep + wc
##
## Df AIC F Pr(>F)
## + fl 0 267.31
RDA with data deleting replicates
#distance calculation
bray_dist_norep <- vegdist(nifhotu_norep,method = "bray")
#rda analysis
otu_dbrda_norep <- dbrda(sqrt(bray_dist_norep)~.,data=env_norep,add = TRUE)
anova(otu_dbrda_norep,permutations = how(nperm = 999))
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: dbrda(formula = sqrt(bray_dist_norep) ~ pH + tn + tp + tk + ep + ek + emo + cec + den + nl + fl + sl + wc, data = env_norep, add = TRUE)
## Df SumOfSqs F Pr(>F)
## Model 12 5.5331 1.185 0.001 ***
## Residual 15 5.8365
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#forward selection
step_forward_norep <- ordistep(dbrda(sqrt(bray_dist_norep)~1,data=env_norep,add = TRUE),scope = formula(otu_dbrda_norep),direction="forward")
##
## Start: sqrt(bray_dist_norep) ~ 1
##
## Df AIC F Pr(>F)
## + pH 1 68.259 2.7227 0.005 **
## + tp 1 69.368 1.6078 0.005 **
## + nl 1 69.635 1.3461 0.015 *
## + tn 1 69.723 1.2602 0.040 *
## + fl 1 69.810 1.1756 0.085 .
## + den 1 69.927 1.0618 0.210
## + emo 1 69.891 1.0965 0.220
## + wc 1 69.954 1.0363 0.310
## + tk 1 69.954 1.0362 0.345
## + cec 1 69.967 1.0233 0.350
## + sl 1 70.006 0.9860 0.425
## + ep 1 70.137 0.8600 0.950
## + ek 1 70.188 0.8112 0.995
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_norep) ~ pH
##
## Df AIC F Pr(>F)
## + tn 1 68.817 1.3213 0.005 **
## + tp 1 68.984 1.1648 0.045 *
## + den 1 69.007 1.1437 0.045 *
## + emo 1 69.037 1.1160 0.105
## + tk 1 69.057 1.0968 0.135
## + fl 1 69.045 1.1083 0.145
## + wc 1 69.059 1.0947 0.185
## + cec 1 69.073 1.0824 0.200
## + nl 1 69.073 1.0820 0.240
## + sl 1 69.105 1.0522 0.330
## + ep 1 69.271 0.8986 0.850
## + ek 1 69.309 0.8633 0.930
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: sqrt(bray_dist_norep) ~ pH + tn
##
## Df AIC F Pr(>F)
## + den 1 69.530 1.1296 0.065 .
## + tp 1 69.520 1.1382 0.070 .
## + emo 1 69.529 1.1304 0.095 .
## + fl 1 69.582 1.0827 0.180
## + cec 1 69.573 1.0903 0.225
## + sl 1 69.595 1.0708 0.240
## + tk 1 69.602 1.0648 0.255
## + wc 1 69.585 1.0799 0.270
## + nl 1 69.693 0.9834 0.630
## + ep 1 69.768 0.9164 0.800
## + ek 1 69.817 0.8734 0.905
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1