Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
# Aggregate basal GIFT into elements
#Get list of present MAGs
present_MAGs <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
filter(rowSums(.[, -1]) != 0) %>%
rownames()
#Align KEGG annotations with present MAGs and remove all-zero and all-one traits
present_MAGs <- present_MAGs[present_MAGs %in% rownames(genome_gifts)]
genome_gifts_filt <- genome_gifts[present_MAGs,] %>%
select_if(~!all(. == 0)) %>% #remove all-zero modules
select_if(~!all(. == 1)) #remove all-one modules
#Filter count table to only contain present MAGs after KEGG filtering
genome_counts_filt_filt <- genome_counts_filt %>%
column_to_rownames(var = "genome")
genome_counts_filt_filt <- genome_counts_filt_filt[present_MAGs,]
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample))
save(genome_gifts_filt,
richness,
neutral,
phylogenetic,
alpha_div,
file = "data/alpha_div_21032025.Rdata")
alpha_div %>%
pivot_longer(-sample, names_to = "alpha", values_to = "value") %>%
left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
group_by(alpha)%>%
summarise(
Aisa_mean=mean(value[Transect=="Aisa"], na.rm=T),
Aisa_sd=sd(value[Transect=="Aisa"], na.rm=T),
Aran_mean=mean(value[Transect=="Aran"], na.rm=T),
Aran_sd=sd(value[Transect=="Aran"], na.rm=T),
Sentein_mean=mean(value[Transect=="Sentein"], na.rm=T),
Sentein_sd=sd(value[Transect=="Sentein"], na.rm=T),
Tourmalet_mean=mean(value[Transect=="Tourmalet"], na.rm=T),
Tourmalet_sd=sd(value[Transect=="Tourmalet"], na.rm=T))%>%
mutate(
Aisa=str_c(round(Aisa_mean,2),"±",round(Aisa_sd,2)),
Aran=str_c(round(Aran_mean,2),"±",round(Aran_sd,2)),
Sentein=str_c(round(Sentein_mean,2),"±",round(Sentein_sd,2)),
Tourmalet=str_c(round(Tourmalet_mean,2),"±",round(Tourmalet_sd,2))) %>%
arrange(-Aisa_mean) %>%
dplyr::select(alpha,Aisa,Aran,Sentein,Tourmalet) %>%
tt()
| alpha |
Aisa |
Aran |
Sentein |
Tourmalet |
| richness |
198.77±65.73 |
153.32±72.26 |
257.32±71.08 |
245.04±68.73 |
| neutral |
87.81±36.55 |
71.13±36.69 |
110.98±45.84 |
101.5±43.7 |
| phylogenetic |
4.98±1.13 |
4.81±1.12 |
5.41±0.86 |
5.29±0.9 |
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == EHI_number)) %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
ggplot(aes(y = value, x = Transect, group=Transect, color=Transect, fill=Transect)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Transect",
breaks=c("Aisa","Aran","Sentein","Tourmalet"),
labels=c("Aisa","Aran","Sentein","Tourmalet"),
values=c("#e5bd5b", "#6b7398","#e2815a", "#876b96")) +
scale_fill_manual(name="Transect",
breaks=c("Aisa","Aran","Sentein","Tourmalet"),
labels=c("Aisa","Aran","Sentein","Tourmalet"),
values=c("#e5bd5b50", "#6b739850","#e2815a50", "#876b9650")) +
facet_wrap(. ~ metric, scales = "free", ncol=4) +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.text.x = element_text(angle = 45, hjust = 1, size=0)
) +
ylab("Alpha diversity")
Regression plots per transect
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == EHI_number)) %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
ggplot(aes(y = value, x = Elevation, group=Transect, colour=Transect)) +
geom_point(alpha=0.4) +
geom_smooth(method = "loess", se=FALSE) +
scale_colour_manual(name="Transect",
breaks=c("Aisa","Aran","Sentein","Tourmalet"),
labels=c("Aisa","Aran","Sentein","Tourmalet"),
values=c("#e5bd5b", "#6b7398","#e2815a", "#876b96")) +
facet_wrap(. ~ metric, scales = "free", ncol=4) +
xlim(800, 2400)+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.text.x = element_text(angle = 45, hjust = 1)
) +
ylab("Alpha diversity") +
xlab("Elevation (m)")
Mixed models per metric
alpha_div_meta<-alpha_div %>%
left_join(sample_metadata, by = join_by(sample == EHI_number))
Richness diversity
Modelq0GLMMNB<- lmerTest::lmer(richness ~ log(sequencing_depth) + Elevation+ I(Elevation^2)+Transect + (1 | Sampling_point), data = alpha_div_meta, REML = FALSE)
summary(Modelq0GLMMNB)
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: richness ~ log(sequencing_depth) + Elevation + I(Elevation^2) + Transect + (1 | Sampling_point)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
1133.5 1157.4 -557.8 1115.5 96
Scaled residuals:
Min 1Q Median 3Q Max
-3.8856 -0.4150 0.0923 0.6956 2.4063
Random effects:
Groups Name Variance Std.Dev.
Sampling_point (Intercept) 321.5 17.93
Residual 2157.7 46.45
Number of obs: 105, groups: Sampling_point, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.696e+03 2.194e+02 9.792e+01 -7.729 9.59e-12 ***
log(sequencing_depth) 1.098e+02 1.130e+01 1.005e+02 9.724 3.81e-16 ***
Elevation -4.040e-03 1.351e-01 2.377e+01 -0.030 0.9764
I(Elevation^2) 1.536e-05 4.379e-05 2.418e+01 0.351 0.7288
TransectAran -1.279e+00 1.785e+01 2.187e+01 -0.072 0.9435
TransectSentein 4.254e+01 1.990e+01 2.313e+01 2.137 0.0434 *
TransectTourmalet -8.757e+00 2.043e+01 2.698e+01 -0.429 0.6717
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) lg(s_) Elevtn I(E^2) TrnscA TrnscS
lg(sqncng_) -0.879
Elevation -0.472 0.004
I(Elevtn^2) 0.453 0.005 -0.992
TransectArn -0.317 0.219 0.152 -0.134
TransctSntn -0.051 -0.101 0.201 -0.191 0.522
TrnsctTrmlt -0.001 -0.248 0.403 -0.421 0.481 0.537
fit warnings:
Some predictor variables are on very different scales: consider rescaling
emmeans::emmeans(Modelq0GLMMNB, pairwise ~ Transect)
$emmeans
Transect emmean SE df lower.CL upper.CL
Aisa 212 16.2 26.9 179 245
Aran 211 15.1 36.1 180 242
Sentein 255 17.9 28.1 218 291
Tourmalet 203 19.0 32.4 165 242
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Aisa - Aran 1.28 21.1 27.9 0.061 0.9999
Aisa - Sentein -42.54 23.4 27.6 -1.816 0.2873
Aisa - Tourmalet 8.76 23.8 31.3 0.368 0.9827
Aran - Sentein -43.82 21.6 32.4 -2.028 0.1989
Aran - Tourmalet 7.48 22.5 41.2 0.332 0.9872
Sentein - Tourmalet 51.30 22.6 30.1 2.267 0.1286
Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 4 estimates
Neutral diversity
Modelq1n <- nlme::lme(neutral ~ log(sequencing_depth) + Elevation+I(Elevation^2)+Transect,
random = ~ 1 | Sampling_point,
data = alpha_div_meta)
summary(Modelq1n)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
1062.022 1085.287 -522.0111
Random effects:
Formula: ~1 | Sampling_point
(Intercept) Residual
StdDev: 14.25385 34.97356
Fixed effects: neutral ~ log(sequencing_depth) + Elevation + I(Elevation^2) + Transect
Value Std.Error DF t-value p-value
(Intercept) -600.7620 166.44662 82 -3.609337 0.0005
log(sequencing_depth) 37.8840 8.52615 82 4.443264 0.0000
Elevation 0.0506 0.10417 82 0.486083 0.6282
I(Elevation^2) 0.0000 0.00003 82 -0.390369 0.6973
TransectAran -1.1070 13.77988 16 -0.080333 0.9370
TransectSentein 18.5815 15.35036 16 1.210495 0.2437
TransectTourmalet -0.9188 15.73085 16 -0.058407 0.9541
Correlation:
(Intr) lg(s_) Elevtn I(E^2) TrnscA TrnscS
log(sequencing_depth) -0.875
Elevation -0.480 0.004
I(Elevation^2) 0.460 0.005 -0.992
TransectAran -0.313 0.214 0.152 -0.134
TransectSentein -0.055 -0.099 0.201 -0.191 0.524
TransectTourmalet -0.010 -0.243 0.402 -0.421 0.485 0.537
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.29352054 -0.63548808 0.06481815 0.65898243 2.05460514
Number of Observations: 105
Number of Groups: 20
emmeans::emmeans(Modelq1n, pairwise ~ Transect)
$emmeans
Transect emmean SE df lower.CL upper.CL
Aisa 93.8 10.6 19 71.6 116
Aran 92.7 10.0 16 71.5 114
Sentein 112.3 11.8 16 87.4 137
Tourmalet 92.8 12.7 16 66.0 120
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Aisa - Aran 1.107 13.8 16 0.080 0.9998
Aisa - Sentein -18.582 15.4 16 -1.210 0.6293
Aisa - Tourmalet 0.919 15.7 16 0.058 0.9999
Aran - Sentein -19.689 14.3 16 -1.379 0.5294
Aran - Tourmalet -0.188 15.1 16 -0.012 1.0000
Sentein - Tourmalet 19.500 15.0 16 1.304 0.5735
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 4 estimates
Phylogenetic diversity
Modelq1p <- nlme::lme(phylogenetic ~ log(sequencing_depth) + Elevation+I(Elevation^2)+Transect,
random = ~ 1 | Sampling_point,
data = alpha_div_meta)
summary(Modelq1p)
Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
355.8819 379.1466 -168.941
Random effects:
Formula: ~1 | Sampling_point
(Intercept) Residual
StdDev: 0.3549572 0.9586922
Fixed effects: phylogenetic ~ log(sequencing_depth) + Elevation + I(Elevation^2) + Transect
Value Std.Error DF t-value p-value
(Intercept) 9.841044 4.504479 82 2.1847243 0.0318
log(sequencing_depth) -0.058510 0.232730 82 -0.2514084 0.8021
Elevation -0.005423 0.002741 82 -1.9784059 0.0512
I(Elevation^2) 0.000002 0.000001 82 2.0697526 0.0416
TransectAran -0.225597 0.361943 16 -0.6232954 0.5419
TransectSentein 0.309787 0.403678 16 0.7674113 0.4540
TransectTourmalet -0.037213 0.415013 16 -0.0896663 0.9297
Correlation:
(Intr) lg(s_) Elevtn I(E^2) TrnscA TrnscS
log(sequencing_depth) -0.882
Elevation -0.466 0.004
I(Elevation^2) 0.447 0.005 -0.992
TransectAran -0.319 0.223 0.152 -0.134
TransectSentein -0.047 -0.103 0.201 -0.191 0.521
TransectTourmalet 0.005 -0.251 0.403 -0.422 0.478 0.536
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.59328141 -0.57869923 -0.05449676 0.59603078 2.19625157
Number of Observations: 105
Number of Groups: 20
emmeans::emmeans(Modelq1p, pairwise ~ Transect)
$emmeans
Transect emmean SE df lower.CL upper.CL
Aisa 4.84 0.278 19 4.26 5.43
Aran 4.62 0.264 16 4.06 5.18
Sentein 5.15 0.310 16 4.50 5.81
Tourmalet 4.81 0.335 16 4.10 5.52
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Aisa - Aran 0.2256 0.362 16 0.623 0.9231
Aisa - Sentein -0.3098 0.404 16 -0.767 0.8679
Aisa - Tourmalet 0.0372 0.415 16 0.090 0.9997
Aran - Sentein -0.5354 0.377 16 -1.422 0.5047
Aran - Tourmalet -0.1884 0.400 16 -0.471 0.9643
Sentein - Tourmalet 0.3470 0.394 16 0.880 0.8152
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 4 estimates
Beta diversity per individual
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, tree = genome_tree)
save(beta_q0n,
beta_q1n,
beta_q1p,
file = "data/beta_div_21032025.Rdata")
Richness across transects and elevational gradients
betadisper(beta_q0n$S, sample_metadata$Transect) %>% permutest(., pairwise=TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.00102 0.0003414 0.0493 999 0.986
Residuals 101 0.69990 0.0069297
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Aisa Aran Sentein Tourmalet
Aisa 0.81200 0.96800 0.806
Aran 0.81748 0.78500 0.990
Sentein 0.97241 0.79372 0.717
Tourmalet 0.77376 0.98856 0.73109
adonis2(beta_q0n$S ~ Transect+Elevation+I(Elevation^2),
data = sample_metadata %>% arrange(match(EHI_number,labels(beta_q0n$S))),
permutations = 999,
strata = sample_metadata %>% arrange(match(EHI_number,labels(beta_q0n$S))) %>% pull(Sampling_point),
by="terms") %>%
broom::tidy() %>%
tt()
Warning in tidy.anova(.): The following column names in ANOVA output were not recognized or transformed: SumOfSqs, R2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Transect |
3 |
2.3447615 |
0.07197986 |
2.672289 |
1.000 |
| Elevation |
1 |
0.7641052 |
0.02345662 |
2.612518 |
0.815 |
| I(Elevation^2) |
1 |
0.5110076 |
0.01568699 |
1.747163 |
0.901 |
| Residual |
99 |
28.9553706 |
0.88887653 |
NA |
NA |
| Total |
104 |
32.5752449 |
1.00000000 |
NA |
NA |
Richness across transects and sampling points
betadisper(beta_q0n$S, sample_metadata$Sampling_point) %>% permutest(., pairwise=TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 19 0.11500 0.0060525 0.8854 999 0.61
Residuals 85 0.58107 0.0068362
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
B Barranc deth Port1 Bifurcation C Cleta Cram2 D Dolmen E Embalse1 F G
B 0.2980000 0.4000000 0.6390000 0.4390000 0.6030000 1.0000000 0.5420000 0.8900000 0.0780000 0.6230000 0.7030000
Barranc deth Port1 0.3040478 0.9890000 0.7370000 0.5530000 0.1470000 0.3260000 0.1290000 0.2460000 0.0170000 0.5240000 0.2760000
Bifurcation 0.3664890 0.9842991 0.7870000 0.5930000 0.2020000 0.3950000 0.2170000 0.3290000 0.0440000 0.6030000 0.3270000
C 0.6236839 0.7325529 0.7650985 0.9460000 0.3920000 0.6570000 0.3460000 0.5700000 0.1040000 0.8780000 0.5350000
Cleta 0.4711597 0.5553415 0.5943158 0.9429777 0.2440000 0.5070000 0.2020000 0.4140000 0.0120000 0.8840000 0.3220000
Cram2 0.6153734 0.1464259 0.2093474 0.3934670 0.2403920 0.6240000 0.8890000 0.7020000 0.2150000 0.3510000 0.9340000
D 0.9966434 0.3158886 0.3799381 0.6340815 0.4928781 0.6216070 0.5610000 0.8900000 0.0700000 0.6420000 0.7390000
Dolmen 0.5231388 0.1309018 0.1908296 0.3583936 0.1878345 0.8995804 0.5334586 0.6200000 0.2770000 0.3020000 0.8200000
E 0.8808842 0.2539881 0.3164401 0.5565211 0.3740076 0.7063156 0.8830570 0.6064810 0.0790000 0.5340000 0.8160000
Embalse1 0.0617047 0.0176032 0.0373442 0.0878521 0.0103868 0.2149619 0.0717360 0.2579120 0.0757024 0.0340000 0.2290000
F 0.6136488 0.5145329 0.5617449 0.8694554 0.8656326 0.3471196 0.6320820 0.2835542 0.5101655 0.0245587 0.4320000
G 0.7116284 0.2480677 0.3195465 0.5151805 0.3218670 0.9262796 0.7203416 0.8378433 0.8004222 0.2082989 0.4372227
H 0.7626589 0.3099177 0.3915869 0.5712920 0.4292667 0.9546452 0.7662136 0.8837448 0.8325213 0.3402873 0.5363681 0.9860756
I 0.6728377 0.2533762 0.3322664 0.5075293 0.3385669 0.9758627 0.6796767 0.9459786 0.7449481 0.3422403 0.4453611 0.9228037
J 0.6339816 0.5069526 0.5555598 0.8593925 0.8488747 0.3606340 0.6513453 0.2965977 0.5305082 0.0273568 0.9827400 0.4535061
K 0.8591639 0.2758125 0.3454357 0.5692563 0.4166744 0.7693488 0.8610388 0.6822271 0.9585672 0.1375477 0.5428266 0.8625865
Presa1 0.4160770 0.9445048 0.9399240 0.7999216 0.6709318 0.2111940 0.4238643 0.1986303 0.3632638 0.0411080 0.6267952 0.3428332
Rest 0.6461720 0.1780386 0.2459707 0.4328678 0.2743747 0.9888836 0.6523837 0.8949505 0.7332441 0.2362772 0.3833795 0.9404571
Senet 0.7796151 0.4636338 0.5219017 0.7947521 0.7480762 0.4711331 0.7896003 0.4064872 0.6816744 0.0606702 0.8741575 0.5801018
Villaler1 0.3412656 0.5873085 0.6169590 0.9991460 0.8915236 0.1747613 0.3709944 0.1251872 0.2500408 0.0034039 0.7490295 0.2239555
H I J K Presa1 Rest Senet Villaler1
B 0.7780000 0.6720000 0.6360000 0.8710000 0.4120000 0.6480000 0.7730000 0.363
Barranc deth Port1 0.3030000 0.2520000 0.5100000 0.2610000 0.9320000 0.1840000 0.4760000 0.593
Bifurcation 0.4030000 0.3590000 0.5670000 0.3700000 0.9350000 0.2310000 0.5250000 0.632
C 0.5850000 0.5050000 0.8450000 0.5660000 0.8030000 0.4220000 0.7960000 1.000
Cleta 0.4360000 0.3490000 0.8480000 0.4250000 0.6500000 0.2810000 0.7510000 0.901
Cram2 0.9470000 0.9700000 0.3490000 0.7920000 0.1970000 0.9830000 0.4580000 0.168
D 0.7680000 0.6750000 0.6600000 0.8640000 0.4580000 0.6860000 0.7930000 0.373
Dolmen 0.8880000 0.9470000 0.3090000 0.6800000 0.1860000 0.8940000 0.4110000 0.118
E 0.8140000 0.7340000 0.5480000 0.9610000 0.3790000 0.7360000 0.6880000 0.261
Embalse1 0.3470000 0.3400000 0.0280000 0.1450000 0.0460000 0.2450000 0.0580000 0.004
F 0.5720000 0.4300000 0.9810000 0.5570000 0.6550000 0.3930000 0.8910000 0.760
G 0.9850000 0.9330000 0.4570000 0.8640000 0.3340000 0.9320000 0.5890000 0.217
H 0.9570000 0.5450000 0.8690000 0.4100000 0.9670000 0.6640000 0.361
I 0.9489899 0.4820000 0.7980000 0.3310000 0.9710000 0.5670000 0.257
J 0.5485147 0.4583109 0.5590000 0.6180000 0.3940000 0.9020000 0.745
K 0.8747050 0.7999933 0.5594890 0.3400000 0.7750000 0.6840000 0.319
Presa1 0.3808394 0.3296386 0.6182053 0.3691535 0.2590000 0.5600000 0.730
Rest 0.9658193 0.9689599 0.3970776 0.7933412 0.2485741 0.4910000 0.180
Senet 0.6419839 0.5604641 0.8903584 0.6848981 0.5634271 0.5074338 0.641
Villaler1 0.3577428 0.2632633 0.7330438 0.3180445 0.7125919 0.2037697 0.6456652
adonis2(beta_q0n$S ~ Transect + Sampling_point,
data = sample_metadata %>% arrange(match(EHI_number,labels(beta_q0n$S))),
permutations = 999,
by="terms") %>%
broom::tidy() %>%
tt()
Warning in tidy.anova(.): The following column names in ANOVA output were not recognized or transformed: SumOfSqs, R2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Transect |
3 |
2.344762 |
0.07197986 |
2.920065 |
0.001 |
| Sampling_point |
16 |
7.479310 |
0.22960103 |
1.746452 |
0.001 |
| Residual |
85 |
22.751174 |
0.69841911 |
NA |
NA |
| Total |
104 |
32.575245 |
1.00000000 |
NA |
NA |
Neutral across transects and elevational gradients
betadisper(beta_q1n$S, sample_metadata$Transect) %>% permutest(., pairwise=TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.00828 0.0027596 0.3593 999 0.787
Residuals 101 0.77583 0.0076815
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Aisa Aran Sentein Tourmalet
Aisa 0.66200 0.68700 0.659
Aran 0.65706 0.38400 0.963
Sentein 0.66556 0.39226 0.346
Tourmalet 0.64203 0.95713 0.33301
adonis2(beta_q1n$S ~ Transect+Elevation+I(Elevation^2),
data = sample_metadata %>% arrange(match(EHI_number,labels(beta_q1n$S))),
permutations = 999,
strata = sample_metadata %>% arrange(match(EHI_number,labels(beta_q1n$S))) %>% pull(Sampling_point),
by="terms") %>%
broom::tidy() %>%
tt()
Warning in tidy.anova(.): The following column names in ANOVA output were not recognized or transformed: SumOfSqs, R2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Transect |
3 |
2.2626912 |
0.07288578 |
2.744266 |
1.000 |
| Elevation |
1 |
0.8398207 |
0.02705229 |
3.055687 |
0.902 |
| I(Elevation^2) |
1 |
0.7328095 |
0.02360525 |
2.666327 |
0.889 |
| Residual |
99 |
27.2090233 |
0.87645668 |
NA |
NA |
| Total |
104 |
31.0443448 |
1.00000000 |
NA |
NA |
Neutral across transects and sampling points
betadisper(beta_q1n$S, sample_metadata$Sampling_point) %>% permutest(., pairwise=TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 19 0.12321 0.0064848 0.8025 999 0.704
Residuals 85 0.68683 0.0080804
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
B Barranc deth Port1 Bifurcation C Cleta Cram2 D Dolmen E Embalse1 F G
B 0.3270000 0.4510000 0.7380000 0.2300000 0.6800000 0.9130000 0.7960000 0.9210000 0.2210000 0.5050000 0.6780000
Barranc deth Port1 0.3007097 0.8930000 0.6910000 0.7140000 0.3410000 0.2840000 0.1800000 0.3940000 0.0400000 0.6410000 0.4580000
Bifurcation 0.4601109 0.8852910 0.7920000 0.9040000 0.5200000 0.4490000 0.3160000 0.5340000 0.0940000 0.8270000 0.5940000
C 0.7174492 0.6542037 0.7953918 0.7810000 0.8480000 0.6910000 0.5720000 0.7930000 0.2250000 0.9210000 0.8920000
Cleta 0.2237026 0.7000553 0.8943929 0.7559349 0.2850000 0.1760000 0.0740000 0.3280000 0.0020000 0.7450000 0.2180000
Cram2 0.6835331 0.3369360 0.5188639 0.8581008 0.2929148 0.6010000 0.4200000 0.7810000 0.0460000 0.6450000 0.9640000
D 0.9214433 0.2621689 0.4154637 0.6659250 0.1726566 0.5911467 0.8770000 0.8360000 0.2570000 0.4180000 0.6190000
Dolmen 0.7882898 0.1723433 0.3081111 0.5549549 0.0888133 0.4281670 0.8722649 0.7450000 0.2790000 0.2800000 0.4580000
E 0.9197167 0.3493618 0.5139431 0.7747148 0.2981315 0.7879994 0.8431427 0.7079785 0.1800000 0.5700000 0.7900000
Embalse1 0.2217370 0.0360113 0.0920715 0.2174621 0.0037684 0.0501232 0.2586208 0.2812940 0.1944003 0.0440000 0.0540000
F 0.4753077 0.6247201 0.7958773 0.9131011 0.7401189 0.6224064 0.4102863 0.2827218 0.5549272 0.0423364 0.6770000
G 0.6910967 0.4437426 0.6101582 0.9026023 0.2133859 0.9628389 0.6061036 0.4510556 0.7874650 0.0633954 0.6697878
H 0.9334010 0.3694505 0.5292628 0.7811830 0.2547006 0.7724845 0.8579103 0.7245387 0.9881613 0.2009472 0.5439375 0.7598872
I 0.8954398 0.3651875 0.5201032 0.7126711 0.3492999 0.6548591 0.9494667 0.9673597 0.8410326 0.4796369 0.5337036 0.7038724
J 0.7971608 0.1943066 0.3316654 0.5793345 0.0726427 0.4325690 0.8800269 0.9928389 0.7187450 0.2756344 0.2854199 0.4292880
K 0.5003251 0.4939369 0.6746651 0.9951405 0.3615639 0.7155778 0.4195839 0.2762432 0.5998292 0.0238066 0.8318935 0.7085591
Presa1 0.3910121 0.8893452 0.9803571 0.7385192 0.8622092 0.4607750 0.3467997 0.2448057 0.4441974 0.0598451 0.7533186 0.5706481
Rest 0.7693406 0.1638242 0.2965449 0.5420928 0.0775536 0.4062297 0.8537620 0.9811738 0.6892404 0.2840266 0.2661229 0.4264451
Senet 0.6275903 0.4890070 0.6633909 0.9582559 0.4798884 0.8528667 0.5515481 0.4076360 0.7153875 0.0683906 0.7919790 0.8919665
Villaler1 0.1714897 0.8092720 0.9958157 0.6740086 0.6792428 0.2052448 0.1307461 0.0637074 0.2334304 0.0026571 0.5909879 0.1382832
H I J K Presa1 Rest Senet Villaler1
B 0.9250000 0.9140000 0.8110000 0.4980000 0.3960000 0.7900000 0.6460000 0.162
Barranc deth Port1 0.3790000 0.3640000 0.2320000 0.5160000 0.9010000 0.1710000 0.5000000 0.815
Bifurcation 0.5620000 0.5360000 0.3560000 0.6770000 0.9810000 0.2920000 0.6540000 0.997
C 0.8000000 0.7410000 0.5950000 0.9970000 0.7460000 0.5590000 0.9630000 0.686
Cleta 0.2510000 0.3590000 0.0810000 0.3890000 0.8740000 0.0730000 0.4840000 0.696
Cram2 0.7570000 0.6420000 0.4650000 0.7320000 0.4950000 0.4150000 0.8520000 0.202
D 0.8640000 0.9590000 0.8790000 0.4180000 0.3460000 0.8770000 0.5630000 0.128
Dolmen 0.7240000 0.9570000 0.9940000 0.2760000 0.2390000 0.9800000 0.4110000 0.070
E 0.9910000 0.8420000 0.7400000 0.6480000 0.4300000 0.7120000 0.7260000 0.250
Embalse1 0.2010000 0.4890000 0.2830000 0.0180000 0.0580000 0.3110000 0.0640000 0.003
F 0.5630000 0.5640000 0.3100000 0.8470000 0.7620000 0.2690000 0.8130000 0.586
G 0.7570000 0.7110000 0.4370000 0.7140000 0.5670000 0.4440000 0.9120000 0.141
H 0.8710000 0.7060000 0.5750000 0.4710000 0.7050000 0.7060000 0.180
I 0.8592512 0.9770000 0.5760000 0.4410000 0.9670000 0.6210000 0.292
J 0.7296426 0.9733256 0.2640000 0.3220000 0.9800000 0.4120000 0.045
K 0.5624808 0.5652795 0.2536636 0.6550000 0.2710000 0.9010000 0.225
Presa1 0.4681700 0.4276233 0.2745667 0.6375916 0.2370000 0.6110000 0.960
Rest 0.7051011 0.9553734 0.9739851 0.2547324 0.2353413 0.3850000 0.048
Senet 0.7024182 0.6407721 0.4104213 0.9053069 0.6118775 0.3879661 0.361
Villaler1 0.1955515 0.3002374 0.0514435 0.2278863 0.9658622 0.0550681 0.3666947
adonis2(beta_q1n$S ~ Transect + Sampling_point,
data = sample_metadata %>% arrange(match(EHI_number,labels(beta_q1n$S))),
permutations = 999,
by="terms") %>%
broom::tidy() %>%
tt()
Warning in tidy.anova(.): The following column names in ANOVA output were not recognized or transformed: SumOfSqs, R2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Transect |
3 |
2.262691 |
0.07288578 |
3.035979 |
0.001 |
| Sampling_point |
16 |
7.665044 |
0.24690628 |
1.928366 |
0.001 |
| Residual |
85 |
21.116610 |
0.68020794 |
NA |
NA |
| Total |
104 |
31.044345 |
1.00000000 |
NA |
NA |
Phylogenetic cross transects and elevational gradients
betadisper(beta_q1p$S, sample_metadata$Transect) %>% permutest(., pairwise=TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.00719 0.002396 0.2197 999 0.903
Residuals 101 1.10166 0.010908
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Aisa Aran Sentein Tourmalet
Aisa 0.75300 0.66400 0.914
Aran 0.73539 0.46700 0.817
Sentein 0.62637 0.44777 0.576
Tourmalet 0.90540 0.81074 0.54936
adonis2(beta_q1p$S ~ Transect+Elevation+I(Elevation^2),
data = sample_metadata %>% arrange(match(EHI_number,labels(beta_q1p$S))),
permutations = 999,
strata = sample_metadata %>% arrange(match(EHI_number,labels(beta_q1p$S))) %>% pull(Sampling_point),
by="terms") %>%
broom::tidy() %>%
tt()
Warning in tidy.anova(.): The following column names in ANOVA output were not recognized or transformed: SumOfSqs, R2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Transect |
3 |
0.24316585 |
0.051781424 |
1.8852255 |
1 |
| Elevation |
1 |
0.16103933 |
0.034292833 |
3.7455356 |
1 |
| I(Elevation^2) |
1 |
0.03529458 |
0.007515874 |
0.8208996 |
1 |
| Residual |
99 |
4.25650571 |
0.906409869 |
NA |
NA |
| Total |
104 |
4.69600548 |
1.000000000 |
NA |
NA |
Phylogenetic across transects and sampling points
betadisper(beta_q1p$S, sample_metadata$Sampling_point) %>% permutest(., pairwise=TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 19 0.09500 0.0050001 0.444 999 0.984
Residuals 85 0.95726 0.0112619
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
B Barranc deth Port1 Bifurcation C Cleta Cram2 D Dolmen E Embalse1 F G H I
B 0.20000 0.54700 0.82700 0.55700 0.92200 0.38900 0.53200 0.94600 0.92800 0.39200 0.43100 0.84200 0.45100
Barranc deth Port1 0.21404 0.78900 0.48800 0.44600 0.15000 0.55700 0.28800 0.21700 0.13900 0.67700 0.49300 0.21000 0.86100
Bifurcation 0.50900 0.76828 0.78700 0.79500 0.50500 0.89500 0.66800 0.54400 0.47700 0.96000 0.83200 0.51600 0.94400
C 0.79258 0.46553 0.74569 0.86600 0.81300 0.79100 0.99000 0.83800 0.75100 0.71500 0.86700 0.75500 0.70300
Cleta 0.52031 0.43020 0.77252 0.86122 0.56000 0.85300 0.79800 0.57200 0.47200 0.80800 0.99000 0.47800 0.71900
Cram2 0.90170 0.16120 0.46233 0.80739 0.51262 0.33900 0.61400 0.97300 0.82500 0.34800 0.41300 0.74700 0.40100
D 0.35229 0.55427 0.88885 0.74242 0.82092 0.32644 0.60800 0.39200 0.30500 0.91900 0.84400 0.26900 0.83400
Dolmen 0.50949 0.26380 0.61539 0.98645 0.76845 0.56127 0.51724 0.60100 0.46000 0.56900 0.65200 0.20500 0.57200
E 0.94098 0.22291 0.52466 0.81954 0.54792 0.96482 0.36790 0.54539 0.87700 0.42200 0.45500 0.77000 0.48500
Embalse1 0.91480 0.15868 0.43870 0.72829 0.43725 0.80049 0.28031 0.41568 0.85015 0.32200 0.31700 0.94400 0.39600
F 0.36610 0.63976 0.95047 0.70480 0.75512 0.32853 0.91318 0.50288 0.38181 0.29417 0.79600 0.33400 0.88600
G 0.38935 0.46638 0.79679 0.85647 0.98758 0.39646 0.80611 0.61613 0.40050 0.31881 0.75409 0.12000 0.73400
H 0.82062 0.20529 0.48696 0.72490 0.42822 0.69649 0.25307 0.19638 0.73508 0.91460 0.30404 0.13440 0.45000
I 0.44779 0.84684 0.93845 0.68655 0.69002 0.39241 0.80503 0.53290 0.46133 0.37808 0.86975 0.72122 0.43356
J 0.86155 0.19238 0.46959 0.72464 0.45423 0.75390 0.31543 0.43119 0.80540 0.92923 0.32839 0.36104 0.99610 0.41375
K 0.39287 0.56464 0.89048 0.75441 0.83265 0.36478 0.99706 0.56191 0.41063 0.31820 0.91498 0.82835 0.30874 0.80903
Presa1 0.45074 0.97342 0.84987 0.62799 0.62154 0.39184 0.72890 0.51958 0.46411 0.38551 0.77905 0.68147 0.45855 0.90792
Rest 0.89262 0.15122 0.42791 0.71470 0.41833 0.77321 0.26092 0.37775 0.82479 0.97962 0.27872 0.28606 0.93199 0.36804
Senet 0.44607 0.43634 0.77928 0.85334 0.99257 0.44348 0.80897 0.70104 0.46856 0.36581 0.74756 0.99367 0.29676 0.69827
Villaler1 0.61977 0.28018 0.62196 0.99565 0.76779 0.66041 0.55527 0.96173 0.66056 0.52514 0.52806 0.69068 0.42784 0.54278
J K Presa1 Rest Senet Villaler1
B 0.88100 0.42600 0.49200 0.90100 0.48400 0.672
Barranc deth Port1 0.19300 0.60700 0.97400 0.14900 0.46200 0.303
Bifurcation 0.48700 0.90700 0.87900 0.44400 0.81300 0.672
C 0.75800 0.77500 0.66500 0.74200 0.87200 0.998
Cleta 0.48700 0.84400 0.67400 0.45600 0.99400 0.807
Cram2 0.77300 0.38700 0.42300 0.79100 0.49300 0.711
D 0.35800 0.99700 0.77100 0.27300 0.84200 0.607
Dolmen 0.47900 0.61400 0.59200 0.40800 0.75800 0.965
E 0.85000 0.46000 0.50800 0.84300 0.49200 0.695
Embalse1 0.93500 0.34900 0.40900 0.97200 0.39200 0.544
F 0.36100 0.92900 0.82200 0.28600 0.78700 0.553
G 0.39500 0.85700 0.72000 0.29200 0.99900 0.749
H 0.99400 0.34300 0.51400 0.94200 0.31800 0.469
I 0.44600 0.81800 0.91800 0.39100 0.73400 0.584
J 0.41000 0.47400 0.95700 0.41900 0.588
K 0.35067 0.75900 0.32400 0.84600 0.631
Presa1 0.41610 0.73036 0.39600 0.70200 0.553
Rest 0.94455 0.30013 0.37939 0.33900 0.542
Senet 0.39395 0.82596 0.64713 0.34089 0.745
Villaler1 0.52680 0.58814 0.51784 0.49637 0.72374
adonis2(beta_q1p$S ~ Transect + Sampling_point,
data = sample_metadata %>% arrange(match(EHI_number,labels(beta_q1p$S))),
permutations = 999,
by="terms") %>%
broom::tidy() %>%
tt()
Warning in tidy.anova(.): The following column names in ANOVA output were not recognized or transformed: SumOfSqs, R2
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Transect |
3 |
0.2431659 |
0.05178142 |
2.039987 |
0.008 |
| Sampling_point |
16 |
1.0755147 |
0.22902757 |
1.691774 |
0.004 |
| Residual |
85 |
3.3773249 |
0.71919100 |
NA |
NA |
| Total |
104 |
4.6960055 |
1.00000000 |
NA |
NA |
NMDS across transects and elevational gradient
p0<-beta_q0n$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
group_by(Transect, Elevation) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=factor(Elevation), shape=Transect)) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
guides(col = guide_legend(ncol = 2))
p1<-beta_q1n$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
group_by(Transect, Elevation) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=factor(Elevation), shape=Transect)) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
guides(col = guide_legend(ncol = 2))
p2<-beta_q1p$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
group_by(Transect, Elevation) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=factor(Elevation), shape=Transect)) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
guides(col = guide_legend(ncol = 2))
ggarrange(p0, p1, p2, ncol=3, nrow=1, common.legend = TRUE, legend="right")
Spatial distance decay in neutral beta diversity
#Create community composition matrix with elevation difference data and transect assignation
# Convert dissimilarity matrix to long format
dissimilarity_long <- reshape2::melt(as.matrix(beta_q1n$S))
# Rename columns for clarity
colnames(dissimilarity_long) <- c("Sample1", "Sample2", "Dissimilarity")
## Create vector of elevations per transect
elevation_vectors <- sample_metadata %>%
split(.$Transect) %>%
lapply(function(df) setNames(df$Elevation, df$EHI_number))
# Compute pairwise absolute elevation differences per transect
elevation_dist_matrices <- lapply(elevation_vectors, function(ev) as.dist(dist(ev, method = "manhattan")))
elevation_dist_long <- lapply(names(elevation_dist_matrices), function(transect) {
df <- reshape::melt(as.matrix(elevation_dist_matrices[[transect]]))
colnames(df) <- c("Sample1", "Sample2", "ElevationDifference")
df$Transect <- transect # Add transect info
return(df)
})
Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the caller; using TRUE
Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the caller; using TRUE
Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the caller; using TRUE
Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the caller; using TRUE
Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the caller; using TRUE
Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the caller; using TRUE
Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the caller; using TRUE
Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by the caller; using TRUE
#Combine the four transects into a single df
elevation_dist_long <- do.call(rbind, elevation_dist_long)
# Merge and keep only samples that appear in both matrices
distance_decay_data <- merge(dissimilarity_long, elevation_dist_long, by = c("Sample1", "Sample2"))
##Remove self comparisons
distance_decay_data<-distance_decay_data %>%
filter(Dissimilarity!=0)
ggplot(distance_decay_data, aes(x = ElevationDifference, y = Dissimilarity, color = Transect)) +
geom_point(alpha = 0.2) +
geom_smooth(method = "lm", se = TRUE) +
scale_colour_manual(name="Transect",
breaks=c("Aisa","Aran","Sentein","Tourmalet"),
labels=c("Aisa","Aran","Sentein","Tourmalet"),
values=c("#e5bd5b", "#6b7398","#e2815a", "#876b96")) +
theme_classic() +
xlab("Elevation difference (m)") +
ylab("Neutral beta diversity dissimilarity") +
theme(legend.title = element_blank())
model <- lm(Dissimilarity ~ ElevationDifference, data = distance_decay_data)
summary(model)
FALSE
FALSE Call:
FALSE lm(formula = Dissimilarity ~ ElevationDifference, data = distance_decay_data)
FALSE
FALSE Residuals:
FALSE Min 1Q Median 3Q Max
FALSE -0.36620 -0.08205 -0.00152 0.07700 0.26580
FALSE
FALSE Coefficients:
FALSE Estimate Std. Error t value Pr(>|t|)
FALSE (Intercept) 7.342e-01 3.585e-03 204.800 < 2e-16 ***
FALSE ElevationDifference 5.045e-05 7.014e-06 7.193 8.05e-13 ***
FALSE ---
FALSE Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
FALSE
FALSE Residual standard error: 0.118 on 2836 degrees of freedom
FALSE Multiple R-squared: 0.01792, Adjusted R-squared: 0.01757
FALSE F-statistic: 51.74 on 1 and 2836 DF, p-value: 8.055e-13
Beta diversity per population
richness_beta <- genome_counts_filt %>%
pivot_longer(!genome,names_to="sample",values_to = "counts") %>%
left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
mutate(Sampling_point=factor(Sampling_point)) %>%
group_by(Sampling_point) %>%
group_split() %>%
map_dbl(., ~ .x %>%
select(genome, sample, counts) %>%
pivot_wider(names_from = sample, values_from = counts) %>%
column_to_rownames(var = "genome") %>%
tss() %>%
as.data.frame() %>%
hilldiss(data=., metric="S", q = 0)
)
names(richness_beta) <- unique(sample_metadata$Sampling_point) %>% sort()
neutral_beta <- genome_counts_filt %>%
pivot_longer(!genome,names_to="sample",values_to = "counts") %>%
left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
mutate(Sampling_point=factor(Sampling_point)) %>%
group_by(Sampling_point) %>%
group_split() %>%
map_dbl(., ~ .x %>%
select(genome, sample, counts) %>%
pivot_wider(names_from = sample, values_from = counts) %>%
column_to_rownames(var = "genome") %>%
tss() %>%
as.data.frame() %>%
hilldiss(data=., metric="S", q = 1)
)
names(neutral_beta) <- unique(sample_metadata$Sampling_point) %>% sort()
phylogenetic_beta <- genome_counts_filt %>%
pivot_longer(!genome,names_to="sample",values_to = "counts") %>%
left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
mutate(Sampling_point=factor(Sampling_point)) %>%
group_by(Sampling_point) %>%
group_split() %>%
map_dbl(., ~ .x %>%
select(genome, sample, counts) %>%
pivot_wider(names_from = sample, values_from = counts) %>%
column_to_rownames(var = "genome") %>%
tss() %>%
as.data.frame() %>%
hilldiss(data=., metric="S", q = 1, tree = genome_tree)
)
names(phylogenetic_beta) <- unique(sample_metadata$Sampling_point) %>% sort()
# Merge all metrics
beta_div <- bind_rows(richness_beta,neutral_beta,phylogenetic_beta) %>%
t() %>%
as.data.frame() %>%
rownames_to_column(var="Sampling_point") %>%
as_tibble() %>%
rename("richness"=V1,"neutral"=V2,"phylogenetic"=V3)
save(richness_beta,
neutral_beta,
phylogenetic_beta,
beta_div,
file = "data/beta_div_pop_21032025.Rdata")
beta_div %>%
pivot_longer(-Sampling_point, names_to = "metric", values_to = "value") %>%
left_join(sample_metadata %>% select(Transect,Sampling_point, Elevation) %>% unique, by = "Sampling_point") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
ggplot(aes(y = value, x = Transect,group=Transect, fill=Transect)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(aes(color = Elevation), alpha = 0.5, size = 2) +
scale_colour_gradient(low="#a7f0ba", high="#044319") +
scale_fill_manual(name="Transect",
breaks=c("Aisa","Aran","Sentein","Tourmalet"),
labels=c("Aisa","Aran","Sentein","Tourmalet"),
values=c("#e5bd5b50", "#6b739850","#e2815a50", "#876b9650")) +
facet_wrap(. ~ metric, scales = "free", ncol=4) +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=2) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.text.x = element_text(angle = 45, hjust = 1)
) +
ylab("Beta diversity")
beta_div_met<-beta_div %>%
pivot_longer(-Sampling_point, names_to = "metric", values_to = "value") %>%
left_join(sample_metadata %>% select(Transect,Sampling_point, Elevation)) %>%
mutate(Zones = case_when(
Transect == "Aisa" & Elevation == 1450 ~ "Z0",
Transect == "Aisa" & Elevation == 1650 ~ "Z1",
Transect == "Aisa" & Elevation == 1850 ~ "Z2",
Transect == "Aran" & Elevation == 1000 ~ "Z0",
Transect == "Aran" & Elevation == 1500 ~ "Z1",
Transect == "Aran" & Elevation == 1650 ~ "Z2",
Transect == "Aran" & Elevation == 1850 ~ "Z3",
Transect == "Sentein" & Elevation == 941 ~ "Z0",
Transect == "Sentein" & Elevation == 1260 ~ "Z1",
Transect == "Sentein" & Elevation == 1628 ~ "Z2",
Transect == "Sentein" & Elevation == 1873 ~ "Z3",
Transect == "Tourmalet" & Elevation == 953 ~ "Z0",
Transect == "Tourmalet" & Elevation == 1255 ~ "Z1",
Transect == "Tourmalet" & Elevation == 1797 ~ "Z2",
Transect == "Tourmalet" & Elevation == 2306 ~ "Z3",
TRUE ~ NA_character_)) %>%
filter(Zones!="NA") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional")))
beta_div_met$Zones<-as.factor(beta_div_met$Zones)
beta_div_met$Tras<-as.factor(beta_div_met$Zones)
beta_div_met$metric <- factor(beta_div_met$metric, levels = c("richness", "neutral", "phylogenetic"))
bxp1<-beta_div_met%>%
ggplot(aes(y = value, x = Sampling_point, group=Sampling_point, colour=Transect, fill=Transect)) +
geom_boxplot() +
geom_jitter() +
#scale_colour_gradient(low="#a7f0ba", high="#044319") +
scale_colour_manual(name="Transect",
breaks=c("Aisa","Aran","Sentein","Tourmalet"),
labels=c("Aisa","Aran","Sentein","Tourmalet"),
values=c("#e5bd5b", "#6b7398","#e2815a", "#876b96")) +
scale_fill_manual(name="Transect",
breaks=c("Aisa","Aran","Sentein","Tourmalet"),
labels=c("Aisa","Aran","Sentein","Tourmalet"),
values=c("#e5bd5b50", "#6b739850","#e2815a50", "#876b9650")) +
facet_wrap(. ~ metric+Transect, scales = "free", ncol=4) +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=2) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.text.x = element_text(angle = 45, hjust = 1)
) +
ylab("Beta diversity")
# Add significance bars with p-values
bxp2<-bxp1 +
geom_signif(comparisons = list(c("Z0", "Z1"), c("Z1", "Z2"), c("Z2", "Z3"), c("Z0", "Z2"), c("Z1", "Z3")),
map_signif_level = TRUE)
bxp2+geom_signif(comparisons = list(c("Z0", "Z1"), c("Z1", "Z2")),
map_signif_level = TRUE)
Transit zones
beta_div_pop<-beta_div %>%
pivot_longer(-Sampling_point, names_to = "metric", values_to = "value") %>%
left_join(sample_metadata %>% select(Transect,Sampling_point, Elevation) %>% unique, by = "Sampling_point")
beta_div_pop <- beta_div_pop %>%
mutate(Zones = case_when(
Transect == "Aisa" & Elevation == 1450 ~ "Z0",
Transect == "Aisa" & Elevation == 1650 ~ "Z1",
Transect == "Aisa" & Elevation == 1850 ~ "Z2",
Transect == "Aran" & Elevation == 1000 ~ "Z0",
Transect == "Aran" & Elevation == 1500 ~ "Z1",
Transect == "Aran" & Elevation == 1650 ~ "Z2",
Transect == "Aran" & Elevation == 1850 ~ "Z3",
Transect == "Sentein" & Elevation == 941 ~ "Z0",
Transect == "Sentein" & Elevation == 1260 ~ "Z1",
Transect == "Sentein" & Elevation == 1628 ~ "Z2",
Transect == "Sentein" & Elevation == 1873 ~ "Z3",
Transect == "Tourmalet" & Elevation == 953 ~ "Z0",
Transect == "Tourmalet" & Elevation == 1255 ~ "Z1",
Transect == "Tourmalet" & Elevation == 1797 ~ "Z2",
Transect == "Tourmalet" & Elevation == 2306 ~ "Z3",
TRUE ~ NA_character_)) %>%
filter(Zones!="NA")
beta_div_pop %>%
#pivot_longer(-Sampling_point, names_to = "metric", values_to = "value") %>%
#left_join(sample_metadata %>% select(Transect,Sampling_point, Elevation) %>% unique, by = "Sampling_point") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
ggplot(aes(y = value, x = Elevation ,group=Elevation, colour=Transect)) +
#geom_boxplot(outlier.shape = NA) +
geom_point(aes(color = Transect), alpha = 0.5, size = 2) +
geom_line(aes(group = interaction(Transect, metric)), position = position_dodge(width = 0.5))+
#scale_color_viridis_c(name = "Elevation", option = "viridis") +
scale_color_manual(name="Transect",
breaks=c("Aisa","Aran","Sentein","Tourmalet"),
labels=c("Aisa","Aran","Sentein","Tourmalet"),
values=c("#e5bd5b50", "#6b739850","#e2815a50", "#876b9650")) +
facet_wrap(. ~ metric, scales = "free", ncol=4) +
coord_cartesian(xlim = c(800, 2400)) +
#stat_compare_means(size=2) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.text.x = element_text(angle = 45, hjust = 1)
) +
ylab("Beta diversity") +
xlab("Elevation (m)")
Aisa
samples_to_keep_aisa <- sample_metadata %>%
filter(Transect == "Aisa") %>%
filter(Elevation!=1250)%>%
select(EHI_number) %>%
pull()
subset_meta_aisa <- sample_metadata %>%
filter(Transect == "Aisa")%>%
filter(Elevation!=1250)%>%
mutate(Zones = case_when(
Transect == "Aisa" & Elevation == 1450 ~ "Z0",
Transect == "Aisa" & Elevation == 1650 ~ "Z1",
Transect == "Aisa" & Elevation == 1850 ~ "Z2"))
Richness
richness_aisa<-as.matrix(beta_q0n$S)
richness_aisa <- as.dist(richness_aisa[rownames(richness_aisa) %in% samples_to_keep_aisa,
colnames(richness_aisa) %in% samples_to_keep_aisa])
betadisper(richness_aisa, subset_meta_aisa$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.003772 0.0018859 0.2748 999 0.753
Residuals 13 0.089215 0.0068627
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2
Z0 0.63300 0.405
Z1 0.64435 0.933
Z2 0.42066 0.94162
adonis2(richness_aisa ~ Zones,
data = subset_meta_aisa %>% arrange(match(EHI_number,labels(richness_aisa))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
2 |
0.8534401 |
0.1924462 |
1.548999 |
0.006 |
| Residual |
13 |
3.5812551 |
0.8075538 |
NA |
NA |
| Total |
15 |
4.4346952 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(richness_aisa, subset_meta_aisa$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z2 vs Z1 |
1 |
0.3371409 |
1.0712622 |
0.11809406 |
0.379 |
1 |
|
| Z2 vs Z0 |
1 |
0.2646211 |
0.9327129 |
0.08531395 |
0.595 |
1 |
|
| Z1 vs Z0 |
1 |
0.2917045 |
1.0008819 |
0.11119821 |
0.418 |
1 |
|
Neutral
neutral_aisa<-as.matrix(beta_q1n$S)
neutral_aisa <- as.dist(neutral_aisa[rownames(neutral_aisa) %in% samples_to_keep_aisa,
colnames(neutral_aisa) %in% samples_to_keep_aisa])
betadisper(neutral_aisa, subset_meta_aisa$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.001232 0.0006162 0.0722 999 0.932
Residuals 13 0.110978 0.0085368
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2
Z0 0.92900 0.682
Z1 0.92770 0.823
Z2 0.67219 0.81707
adonis2(neutral_aisa ~ Zones,
data = subset_meta_aisa %>% arrange(match(EHI_number,labels(neutral_aisa))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
2 |
0.8375024 |
0.1980923 |
1.605671 |
0.009 |
| Residual |
13 |
3.3903377 |
0.8019077 |
NA |
NA |
| Total |
15 |
4.2278401 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(neutral_aisa, subset_meta_aisa$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z2 vs Z1 |
1 |
0.3707416 |
1.2794356 |
0.13787860 |
0.145 |
0.435 |
|
| Z2 vs Z0 |
1 |
0.2742152 |
1.0032615 |
0.09117856 |
0.446 |
1.000 |
|
| Z1 vs Z0 |
1 |
0.2363446 |
0.8463151 |
0.09566866 |
0.681 |
1.000 |
|
Phylogenetic
phylo_aisa<-as.matrix(beta_q1p$S)
phylo_aisa <- as.dist(phylo_aisa[rownames(phylo_aisa) %in% samples_to_keep_aisa,
colnames(phylo_aisa) %in% samples_to_keep_aisa])
betadisper(phylo_aisa, subset_meta_aisa$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.016704 0.0083522 1.2726 999 0.294
Residuals 13 0.085321 0.0065631
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2
Z0 0.26300 0.181
Z1 0.27329 0.698
Z2 0.17102 0.72653
adonis2(phylo_aisa ~ Zones,
data = subset_meta_aisa %>% arrange(match(EHI_number,labels(phylo_aisa))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
2 |
0.1727056 |
0.26169 |
2.303891 |
0.073 |
| Residual |
13 |
0.4872569 |
0.73831 |
NA |
NA |
| Total |
15 |
0.6599625 |
1.00000 |
NA |
NA |
pairwise <- pairwise.adonis(phylo_aisa, subset_meta_aisa$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z2 vs Z1 |
1 |
0.05218133 |
0.9595539 |
0.1070984 |
0.386 |
1 |
|
| Z2 vs Z0 |
1 |
0.04703713 |
1.1244722 |
0.1010809 |
0.355 |
1 |
|
| Z1 vs Z0 |
1 |
0.03893120 |
1.1039304 |
0.1212587 |
0.344 |
1 |
|
Functional
func_aisa<-as.matrix(beta_q0n$S)
func_aisa <- as.dist(func_aisa[rownames(func_aisa) %in% samples_to_keep_aisa,
colnames(func_aisa) %in% samples_to_keep_aisa])
betadisper(func_aisa, subset_meta_aisa$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.003772 0.0018859 0.2748 999 0.756
Residuals 13 0.089215 0.0068627
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2
Z0 0.63700 0.421
Z1 0.64435 0.947
Z2 0.42066 0.94162
adonis2(func_aisa ~ Zones,
data = subset_meta_aisa %>% arrange(match(EHI_number,labels(func_aisa))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
2 |
0.8534401 |
0.1924462 |
1.548999 |
0.003 |
| Residual |
13 |
3.5812551 |
0.8075538 |
NA |
NA |
| Total |
15 |
4.4346952 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(func_aisa, subset_meta_aisa$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z2 vs Z1 |
1 |
0.3371409 |
1.0712622 |
0.11809406 |
0.358 |
1 |
|
| Z2 vs Z0 |
1 |
0.2646211 |
0.9327129 |
0.08531395 |
0.601 |
1 |
|
| Z1 vs Z0 |
1 |
0.2917045 |
1.0008819 |
0.11119821 |
0.452 |
1 |
|
Aran
samples_to_keep_Aran <- sample_metadata %>%
filter(Transect == "Aran") %>%
filter(Elevation!=1080) %>%
filter(Elevation!=1340) %>%
select(EHI_number) %>%
pull()
subset_meta_Aran <- sample_metadata %>%
filter(Transect == "Aran")%>%
filter(Elevation!=1080) %>%
filter(Elevation!=1340) %>%
mutate(Zones = case_when(
Transect == "Aran" & Elevation == 1000 ~ "Z0",
Transect == "Aran" & Elevation == 1500 ~ "Z1",
Transect == "Aran" & Elevation == 1650 ~ "Z2",
Transect == "Aran" & Elevation == 1850 ~ "Z3"))
Richness
richness_Aran<-as.matrix(beta_q0n$S)
richness_Aran <- as.dist(richness_Aran[rownames(richness_Aran) %in% samples_to_keep_Aran,
colnames(richness_Aran) %in% samples_to_keep_Aran])
betadisper(richness_Aran, subset_meta_Aran$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.01499 0.0049968 0.3525 999 0.796
Residuals 21 0.29765 0.0141739
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.72700 0.93900 0.621
Z1 0.69606 0.53000 0.366
Z2 0.93379 0.55617 0.562
Z3 0.61474 0.35171 0.54452
adonis2(richness_Aran ~ Zones,
data = subset_meta_Aran %>% arrange(match(EHI_number,labels(richness_Aran))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
1.688494 |
0.2111284 |
1.873434 |
0.001 |
| Residual |
21 |
6.308983 |
0.7888716 |
NA |
NA |
| Total |
24 |
7.997477 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(richness_Aran, subset_meta_Aran$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z0 vs Z3 |
1 |
0.2683277 |
0.7494360 |
0.06971863 |
0.744 |
1 |
|
| Z0 vs Z1 |
1 |
0.2291126 |
0.7202304 |
0.06718423 |
0.899 |
1 |
|
| Z0 vs Z2 |
1 |
0.3661773 |
1.1061894 |
0.09137387 |
0.264 |
1 |
|
| Z3 vs Z1 |
1 |
0.3533102 |
1.0446704 |
0.09458593 |
0.384 |
1 |
|
| Z3 vs Z2 |
1 |
0.3555732 |
1.0179833 |
0.08470500 |
0.399 |
1 |
|
| Z1 vs Z2 |
1 |
0.3663377 |
1.1704367 |
0.09617047 |
0.213 |
1 |
|
Neutral
neutral_Aran<-as.matrix(beta_q1n$S)
neutral_Aran <- as.dist(neutral_Aran[rownames(neutral_Aran) %in% samples_to_keep_Aran,
colnames(neutral_Aran) %in% samples_to_keep_Aran])
betadisper(neutral_Aran, subset_meta_Aran$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.00535 0.0017819 0.0935 999 0.959
Residuals 21 0.40001 0.0190481
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.69500 0.77800 0.998
Z1 0.70904 0.86000 0.698
Z2 0.76165 0.86873 0.715
Z3 0.99716 0.69886 0.74882
adonis2(neutral_Aran ~ Zones,
data = subset_meta_Aran %>% arrange(match(EHI_number,labels(neutral_Aran))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
1.651330 |
0.2207076 |
1.982508 |
0.002 |
| Residual |
21 |
5.830648 |
0.7792924 |
NA |
NA |
| Total |
24 |
7.481978 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(neutral_Aran, subset_meta_Aran$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z0 vs Z3 |
1 |
0.2836206 |
0.8488135 |
0.07824022 |
0.650 |
1 |
|
| Z0 vs Z1 |
1 |
0.1899430 |
0.5979893 |
0.05642479 |
0.989 |
1 |
|
| Z0 vs Z2 |
1 |
0.3179402 |
1.0153314 |
0.08450299 |
0.411 |
1 |
|
| Z3 vs Z1 |
1 |
0.3268096 |
1.0309829 |
0.09346247 |
0.430 |
1 |
|
| Z3 vs Z2 |
1 |
0.3075466 |
0.9839904 |
0.08210874 |
0.487 |
1 |
|
| Z1 vs Z2 |
1 |
0.3040944 |
1.0219976 |
0.08501063 |
0.434 |
1 |
|
Phylogenetic
phylo_Aran<-as.matrix(beta_q1p$S)
phylo_Aran <- as.dist(phylo_Aran[rownames(phylo_Aran) %in% samples_to_keep_Aran,
colnames(phylo_Aran) %in% samples_to_keep_Aran])
betadisper(phylo_Aran, subset_meta_Aran$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.01979 0.0065956 0.2464 999 0.885
Residuals 21 0.56207 0.0267651
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.69500 0.93100 0.663
Z1 0.65487 0.56800 0.868
Z2 0.91960 0.52777 0.551
Z3 0.60010 0.85366 0.50253
adonis2(phylo_Aran ~ Zones,
data = subset_meta_Aran %>% arrange(match(EHI_number,labels(phylo_Aran))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
0.2448675 |
0.1537826 |
1.272106 |
0.233 |
| Residual |
21 |
1.3474291 |
0.8462174 |
NA |
NA |
| Total |
24 |
1.5922966 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(phylo_Aran, subset_meta_Aran$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z0 vs Z3 |
1 |
0.07013449 |
0.8683266 |
0.07989515 |
0.532 |
1 |
|
| Z0 vs Z1 |
1 |
0.04948648 |
0.7852274 |
0.07280582 |
0.641 |
1 |
|
| Z0 vs Z2 |
1 |
0.05654666 |
1.2175475 |
0.09965564 |
0.306 |
1 |
|
| Z3 vs Z1 |
1 |
0.03794568 |
0.4142491 |
0.03977714 |
0.881 |
1 |
|
| Z3 vs Z2 |
1 |
0.07841010 |
1.0826500 |
0.08960368 |
0.339 |
1 |
|
| Z1 vs Z2 |
1 |
0.03769619 |
0.6696805 |
0.05738636 |
0.595 |
1 |
|
Functional
func_Aran<-as.matrix(beta_q0n$S)
func_Aran <- as.dist(func_Aran[rownames(func_Aran) %in% samples_to_keep_Aran,
colnames(func_Aran) %in% samples_to_keep_Aran])
betadisper(func_Aran, subset_meta_Aran$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.01499 0.0049968 0.3525 999 0.813
Residuals 21 0.29765 0.0141739
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.69100 0.94000 0.622
Z1 0.69606 0.56500 0.368
Z2 0.93379 0.55617 0.549
Z3 0.61474 0.35171 0.54452
adonis2(func_Aran ~ Zones,
data = subset_meta_Aran %>% arrange(match(EHI_number,labels(func_Aran))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
1.688494 |
0.2111284 |
1.873434 |
0.001 |
| Residual |
21 |
6.308983 |
0.7888716 |
NA |
NA |
| Total |
24 |
7.997477 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(func_Aran, subset_meta_Aran$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z0 vs Z3 |
1 |
0.2683277 |
0.7494360 |
0.06971863 |
0.750 |
1 |
|
| Z0 vs Z1 |
1 |
0.2291126 |
0.7202304 |
0.06718423 |
0.886 |
1 |
|
| Z0 vs Z2 |
1 |
0.3661773 |
1.1061894 |
0.09137387 |
0.266 |
1 |
|
| Z3 vs Z1 |
1 |
0.3533102 |
1.0446704 |
0.09458593 |
0.396 |
1 |
|
| Z3 vs Z2 |
1 |
0.3555732 |
1.0179833 |
0.08470500 |
0.380 |
1 |
|
| Z1 vs Z2 |
1 |
0.3663377 |
1.1704367 |
0.09617047 |
0.222 |
1 |
|
Sentein
samples_to_keep_Sentein <- sample_metadata %>%
filter(Transect == "Sentein") %>%
select(EHI_number) %>%
pull()
subset_meta_Sentein <- sample_metadata %>%
filter(Transect == "Sentein")%>%
mutate(Zones = case_when(
Transect == "Sentein" & Elevation == 941 ~ "Z0",
Transect == "Sentein" & Elevation == 1260 ~ "Z1",
Transect == "Sentein" & Elevation == 1628 ~ "Z2",
Transect == "Sentein" & Elevation == 1873 ~ "Z3"))
Richness
richness_Sentein<-as.matrix(beta_q0n$S)
richness_Sentein <- as.dist(richness_Sentein[rownames(richness_Sentein) %in% samples_to_keep_Sentein,
colnames(richness_Sentein) %in% samples_to_keep_Sentein])
betadisper(richness_Sentein, subset_meta_Sentein$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.000586 0.0001953 0.0274 999 0.996
Residuals 15 0.106954 0.0071302
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.96600 0.87400 0.996
Z1 0.95843 0.85100 0.971
Z2 0.85766 0.83351 0.862
Z3 0.99725 0.96275 0.85693
adonis2(richness_Sentein ~ Zones,
data = subset_meta_Sentein %>% arrange(match(EHI_number,labels(richness_Sentein))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
0.9080638 |
0.2142842 |
1.363624 |
0.004 |
| Residual |
15 |
3.3295968 |
0.7857158 |
NA |
NA |
| Total |
18 |
4.2376606 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(richness_Sentein, subset_meta_Sentein$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z3 vs Z2 |
1 |
0.2375436 |
1.0174160 |
0.11282789 |
0.407 |
1 |
|
| Z3 vs Z1 |
1 |
0.2669639 |
1.1015029 |
0.13596278 |
0.290 |
1 |
|
| Z3 vs Z0 |
1 |
0.2025126 |
0.8581745 |
0.09687938 |
0.759 |
1 |
|
| Z2 vs Z1 |
1 |
0.2427161 |
1.0148385 |
0.12661996 |
0.381 |
1 |
|
| Z2 vs Z0 |
1 |
0.1536841 |
0.6590676 |
0.07611300 |
1.000 |
1 |
|
| Z1 vs Z0 |
1 |
0.2569696 |
1.0617346 |
0.13170051 |
0.309 |
1 |
|
Neutral
neutral_Sentein<-as.matrix(beta_q1n$S)
neutral_Sentein <- as.dist(neutral_Sentein[rownames(neutral_Sentein) %in% samples_to_keep_Sentein,
colnames(neutral_Sentein) %in% samples_to_keep_Sentein])
betadisper(neutral_Sentein, subset_meta_Sentein$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.005606 0.0018688 0.1825 999 0.916
Residuals 15 0.153627 0.0102418
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.89000 0.78100 0.422
Z1 0.87111 0.74300 0.338
Z2 0.78640 0.73283 0.884
Z3 0.41948 0.32948 0.86846
adonis2(neutral_Sentein ~ Zones,
data = subset_meta_Sentein %>% arrange(match(EHI_number,labels(neutral_Sentein))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
1.038080 |
0.2525336 |
1.689264 |
0.002 |
| Residual |
15 |
3.072579 |
0.7474664 |
NA |
NA |
| Total |
18 |
4.110659 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(neutral_Sentein, subset_meta_Sentein$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z3 vs Z2 |
1 |
0.2550295 |
1.1940719 |
0.12987411 |
0.233 |
1 |
|
| Z3 vs Z1 |
1 |
0.2722656 |
1.2040506 |
0.14676295 |
0.241 |
1 |
|
| Z3 vs Z0 |
1 |
0.1802123 |
0.8256212 |
0.09354822 |
0.749 |
1 |
|
| Z2 vs Z1 |
1 |
0.2667350 |
1.0964422 |
0.13542272 |
0.283 |
1 |
|
| Z2 vs Z0 |
1 |
0.1291948 |
0.5538189 |
0.06474523 |
0.992 |
1 |
|
| Z1 vs Z0 |
1 |
0.2304475 |
0.9268351 |
0.11692372 |
0.587 |
1 |
|
Phylogenetic
phylo_Sentein<-as.matrix(beta_q1p$S)
phylo_Sentein <- as.dist(phylo_Sentein[rownames(phylo_Sentein) %in% samples_to_keep_Sentein,
colnames(phylo_Sentein) %in% samples_to_keep_Sentein])
betadisper(phylo_Sentein, subset_meta_Sentein$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.00740 0.0024665 0.2359 999 0.943
Residuals 15 0.15682 0.0104546
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.87300 0.77800 0.951
Z1 0.79532 0.72500 0.898
Z2 0.63893 0.62016 0.724
Z3 0.91344 0.83008 0.61693
adonis2(phylo_Sentein ~ Zones,
data = subset_meta_Sentein %>% arrange(match(EHI_number,labels(phylo_Sentein))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
0.1791132 |
0.2402074 |
1.580743 |
0.033 |
| Residual |
15 |
0.5665474 |
0.7597926 |
NA |
NA |
| Total |
18 |
0.7456606 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(phylo_Sentein, subset_meta_Sentein$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z3 vs Z2 |
1 |
0.09032861 |
1.8196959 |
0.18531082 |
0.092 |
0.552 |
|
| Z3 vs Z1 |
1 |
0.05533407 |
1.9062867 |
0.21403833 |
0.118 |
0.708 |
|
| Z3 vs Z0 |
1 |
0.05336036 |
1.8567636 |
0.18837457 |
0.152 |
0.912 |
|
| Z2 vs Z1 |
1 |
0.03341661 |
0.6235568 |
0.08179343 |
0.894 |
1.000 |
|
| Z2 vs Z0 |
1 |
0.02675707 |
0.5326789 |
0.06242810 |
0.906 |
1.000 |
|
| Z1 vs Z0 |
1 |
0.01873231 |
0.6306449 |
0.08264634 |
0.806 |
1.000 |
|
Functional
func_Sentein<-as.matrix(beta_q0n$S)
func_Sentein <- as.dist(func_Sentein[rownames(func_Sentein) %in% samples_to_keep_Sentein,
colnames(func_Sentein) %in% samples_to_keep_Sentein])
betadisper(func_Sentein, subset_meta_Sentein$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.000586 0.0001953 0.0274 999 0.995
Residuals 15 0.106954 0.0071302
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.96200 0.85600 1.000
Z1 0.95843 0.85100 0.966
Z2 0.85766 0.83351 0.887
Z3 0.99725 0.96275 0.85693
adonis2(func_Sentein ~ Zones,
data = subset_meta_Sentein %>% arrange(match(EHI_number,labels(func_Sentein))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
0.9080638 |
0.2142842 |
1.363624 |
0.005 |
| Residual |
15 |
3.3295968 |
0.7857158 |
NA |
NA |
| Total |
18 |
4.2376606 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(func_Sentein, subset_meta_Sentein$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z3 vs Z2 |
1 |
0.2375436 |
1.0174160 |
0.11282789 |
0.418 |
1 |
|
| Z3 vs Z1 |
1 |
0.2669639 |
1.1015029 |
0.13596278 |
0.267 |
1 |
|
| Z3 vs Z0 |
1 |
0.2025126 |
0.8581745 |
0.09687938 |
0.758 |
1 |
|
| Z2 vs Z1 |
1 |
0.2427161 |
1.0148385 |
0.12661996 |
0.382 |
1 |
|
| Z2 vs Z0 |
1 |
0.1536841 |
0.6590676 |
0.07611300 |
1.000 |
1 |
|
| Z1 vs Z0 |
1 |
0.2569696 |
1.0617346 |
0.13170051 |
0.336 |
1 |
|
Tourmalet
samples_to_keep_Tourmalet <- sample_metadata %>%
filter(Transect == "Tourmalet") %>%
filter(Elevation!=1561)%>%
filter(Elevation!=2065)%>%
filter(Elevation!=2134)%>%
select(EHI_number) %>%
pull()
subset_meta_Tourmalet <- sample_metadata %>%
filter(Transect == "Tourmalet")%>%
filter(Elevation!=1561)%>%
filter(Elevation!=2065)%>%
filter(Elevation!=2134)%>%
mutate(Zones = case_when(
Transect == "Tourmalet" & Elevation == 953 ~ "Z0",
Transect == "Tourmalet" & Elevation == 1255 ~ "Z1",
Transect == "Tourmalet" & Elevation == 1797 ~ "Z2",
Transect == "Tourmalet" & Elevation == 2306 ~ "Z3"))
Richness
richness_Tourmalet<-as.matrix(beta_q0n$S)
richness_Tourmalet <- as.dist(richness_Tourmalet[rownames(richness_Tourmalet) %in% samples_to_keep_Tourmalet,
colnames(richness_Tourmalet) %in% samples_to_keep_Tourmalet])
betadisper(richness_Tourmalet, subset_meta_Tourmalet$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.01082 0.0036066 0.4255 999 0.745
Residuals 14 0.11868 0.0084770
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.55300 0.63900 0.318
Z1 0.54233 0.85200 0.786
Z2 0.63198 0.83988 0.548
Z3 0.32313 0.79049 0.53472
adonis2(richness_Tourmalet ~ Zones,
data = subset_meta_Tourmalet %>% arrange(match(EHI_number,labels(richness_Tourmalet))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
1.644463 |
0.3244252 |
2.241031 |
0.001 |
| Residual |
14 |
3.424388 |
0.6755748 |
NA |
NA |
| Total |
17 |
5.068851 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(richness_Tourmalet, subset_meta_Tourmalet$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z3 vs Z2 |
1 |
0.3159344 |
1.0174080 |
0.12689987 |
0.399 |
1.000 |
|
| Z3 vs Z1 |
1 |
0.2326859 |
0.7306096 |
0.09450866 |
0.952 |
1.000 |
|
| Z3 vs Z0 |
1 |
0.3651518 |
1.2919249 |
0.13903738 |
0.166 |
0.996 |
|
| Z2 vs Z1 |
1 |
0.3395145 |
1.0840568 |
0.15302768 |
0.407 |
1.000 |
|
| Z2 vs Z0 |
1 |
0.2885497 |
1.0570176 |
0.13119217 |
0.291 |
1.000 |
|
| Z1 vs Z0 |
1 |
0.3124845 |
1.1122898 |
0.13711170 |
0.278 |
1.000 |
|
Neutral
neutral_Tourmalet<-as.matrix(beta_q1n$S)
neutral_Tourmalet <- as.dist(neutral_Tourmalet[rownames(neutral_Tourmalet) %in% samples_to_keep_Tourmalet,
colnames(neutral_Tourmalet) %in% samples_to_keep_Tourmalet])
betadisper(neutral_Tourmalet, subset_meta_Tourmalet$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.005366 0.0017887 0.2496 999 0.867
Residuals 14 0.100345 0.0071675
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.73600 0.56500 0.814
Z1 0.72491 0.52600 0.683
Z2 0.55399 0.50036 0.251
Z3 0.78813 0.66695 0.24647
adonis2(neutral_Tourmalet ~ Zones,
data = subset_meta_Tourmalet %>% arrange(match(EHI_number,labels(neutral_Tourmalet))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
1.766135 |
0.3549255 |
2.56764 |
0.001 |
| Residual |
14 |
3.209938 |
0.6450745 |
NA |
NA |
| Total |
17 |
4.976073 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(neutral_Tourmalet, subset_meta_Tourmalet$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z3 vs Z2 |
1 |
0.3265025 |
1.1532417 |
0.1414458 |
0.273 |
1.000 |
|
| Z3 vs Z1 |
1 |
0.2578411 |
0.9476322 |
0.1192345 |
0.477 |
1.000 |
|
| Z3 vs Z0 |
1 |
0.4639250 |
1.8123792 |
0.1847033 |
0.045 |
0.270 |
|
| Z2 vs Z1 |
1 |
0.4306105 |
1.4278859 |
0.1922331 |
0.127 |
0.762 |
|
| Z2 vs Z0 |
1 |
0.3114323 |
1.1164685 |
0.1375559 |
0.291 |
1.000 |
|
| Z1 vs Z0 |
1 |
0.4436572 |
1.6559508 |
0.1913078 |
0.081 |
0.486 |
|
Phylogenetic
phylo_Tourmalet<-as.matrix(beta_q1p$S)
phylo_Tourmalet <- as.dist(phylo_Tourmalet[rownames(phylo_Tourmalet) %in% samples_to_keep_Tourmalet,
colnames(phylo_Tourmalet) %in% samples_to_keep_Tourmalet])
betadisper(phylo_Tourmalet, subset_meta_Tourmalet$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.006258 0.002086 0.7418 999 0.558
Residuals 14 0.039367 0.002812
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.83600 0.24200 0.566
Z1 0.83103 0.37900 0.812
Z2 0.24810 0.36948 0.197
Z3 0.56736 0.83334 0.17818
adonis2(phylo_Tourmalet ~ Zones,
data = subset_meta_Tourmalet %>% arrange(match(EHI_number,labels(phylo_Tourmalet))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
0.1959592 |
0.3774918 |
2.829888 |
0.002 |
| Residual |
14 |
0.3231493 |
0.6225082 |
NA |
NA |
| Total |
17 |
0.5191085 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(phylo_Tourmalet, subset_meta_Tourmalet$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z3 vs Z2 |
1 |
0.02640453 |
0.8924459 |
0.11307597 |
0.505 |
1.000 |
|
| Z3 vs Z1 |
1 |
0.02031821 |
0.8798851 |
0.11166217 |
0.524 |
1.000 |
|
| Z3 vs Z0 |
1 |
0.05542631 |
2.3936078 |
0.23029614 |
0.098 |
0.588 |
|
| Z2 vs Z1 |
1 |
0.05981848 |
1.7757976 |
0.22837498 |
0.137 |
0.822 |
|
| Z2 vs Z0 |
1 |
0.01619581 |
0.5022690 |
0.06694894 |
0.840 |
1.000 |
|
| Z1 vs Z0 |
1 |
0.08582117 |
3.3327983 |
0.32254557 |
0.036 |
0.216 |
|
Functional
func_Tourmalet<-as.matrix(beta_q0n$S)
func_Tourmalet <- as.dist(func_Tourmalet[rownames(func_Tourmalet) %in% samples_to_keep_Tourmalet,
colnames(func_Tourmalet) %in% samples_to_keep_Tourmalet])
betadisper(func_Tourmalet, subset_meta_Tourmalet$Zones) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.01082 0.0036066 0.4255 999 0.738
Residuals 14 0.11868 0.0084770
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Z0 Z1 Z2 Z3
Z0 0.52800 0.62700 0.294
Z1 0.54233 0.83600 0.801
Z2 0.63198 0.83988 0.531
Z3 0.32313 0.79049 0.53472
adonis2(func_Tourmalet ~ Zones,
data = subset_meta_Tourmalet %>% arrange(match(EHI_number,labels(func_Tourmalet))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| Model |
3 |
1.644463 |
0.3244252 |
2.241031 |
0.001 |
| Residual |
14 |
3.424388 |
0.6755748 |
NA |
NA |
| Total |
17 |
5.068851 |
1.0000000 |
NA |
NA |
pairwise <- pairwise.adonis(func_Tourmalet, subset_meta_Tourmalet$Zones, perm=999)
pairwise%>%
tt()
| pairs |
Df |
SumsOfSqs |
F.Model |
R2 |
p.value |
p.adjusted |
sig |
| Z3 vs Z2 |
1 |
0.3159344 |
1.0174080 |
0.12689987 |
0.408 |
1.000 |
|
| Z3 vs Z1 |
1 |
0.2326859 |
0.7306096 |
0.09450866 |
0.945 |
1.000 |
|
| Z3 vs Z0 |
1 |
0.3651518 |
1.2919249 |
0.13903738 |
0.162 |
0.972 |
|
| Z2 vs Z1 |
1 |
0.3395145 |
1.0840568 |
0.15302768 |
0.386 |
1.000 |
|
| Z2 vs Z0 |
1 |
0.2885497 |
1.0570176 |
0.13119217 |
0.316 |
1.000 |
|
| Z1 vs Z0 |
1 |
0.3124845 |
1.1122898 |
0.13711170 |
0.250 |
1.000 |
|