# Message and warning turned off for rendering
#Loading packages
library(picante)
library(ggplot2)
library(phytools)
library(dplyr)
library(grid)
library(gridExtra)
library(ggstatsplot)
# Warning turned off for rendering
## Folders, Themes, colors
source("prelude.R")
# Setting up directories
data_dir <-paste0(data_dir_path,"xpdiet")
output_dir <- paste0(output_dir_path,"Figure5_diet")
# Defining the directory containing the data to import
setwd(data_dir)
### Rarefaction analysis : importing data
comm_abundance_1<-read.table("nmt3_abondance_xpfood_tpn_notax.csv",sep=";",header=TRUE,row.names=1)
# Getting rid of OTUs with zero abundance
comm_abundance_1<-comm_abundance_1 %>% select_if((function(col) is.numeric(col) && sum(col) > 0))
S<-specnumber(comm_abundance_1)
raremax <- min(rowSums(comm_abundance_1))
Srare <- rarefy(comm_abundance_1, raremax)
# Defining the output dir
setwd(output_dir)
# Rarefaction
plot(S, Srare, xlab = "Observed No. of Species", ylab = "Rarefied No. of Species")
abline(0, 1)
#Saving plot to file
#pdf("plot_rarefaction1_food.pdf",width=5,height=4)
#print(plot(S, Srare, xlab = "Observed No. of Species", ylab = "Rarefied No. of Species"))
#print(abline(0, 1))
#dev.off()
# Rarefaction curve
rarecurve(comm_abundance_1, step = 20, sample = raremax, col = "blue", cex = 0.6)
#Saving plot to file
#pdf("plot_rarefaction2_food.pdf",width=5,height=4)
#print(rarecurve(comm_abundance_1, step = 20, sample = raremax, col = "blue", cex = 0.6))
#dev.off()
###Diversity analysis
#Defining the directory containing the data to import
setwd(data_dir)
# Working on abundance data from 100 rarefactions
allcom<-read.table("xpdiet_fichier_allcom.txt", header=TRUE, row.names=1)
# Getting rid of OTUs with zero abundance
allcom<-allcom %>% select_if((function(col) is.numeric(col) && sum(col) > 0))
commRS<-decostand(allcom, method="total")
comm<-commRS
#Importing metadata
metadata<-read.csv("nmt3_2_metadata_food.csv",sep=";",header=TRUE,row.names=1)
#Importing phylogeny
phy<-read.tree("nmt3_sequences_mpr.tree")
#Rooting the tree using midpoint rooting
phy<-midpoint.root(phy)
#Getting rid of taxa that are missing in the studied samples
phy<-prune.sample(comm, phy)
#Combining the input data
combined <- match.phylo.comm(phy, comm)
phy <- combined$phy
comm <- combined$comm
metadata <- metadata[rownames(comm), ]
#Checking the index
all.equal(rownames(comm), rownames(metadata))
## [1] TRUE
#Checking the normalisation
rowSums(comm)
## MT-D1-1 MT-D1-10 MT-D1-2 MT-D1-3 MT-D1-4 MT-D1-5 MT-D1-6 MT-D1-7
## 1 1 1 1 1 1 1 1
## MT-D1-9 MT-D2-1 MT-D2-2 MT-D2-3 MT-D2-4 MT-D2-5 MT-D2-6 MT-D2-7
## 1 1 1 1 1 1 1 1
## MT-D2-8 MT-D2-9
## 1 1
#######################################
# Analysis of taxonomic diversity #
#######################################
# Defining the output dir
setwd(output_dir)
# For the boxplot with ggplot2
# OTU richness
df<-data.frame(Nb.OTU= specnumber(comm),Sample=metadata$Class)
# OTU richness boxplot
ggplot(df,aes(x=Sample,y=Nb.OTU))+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey",alpha=0.9)+geom_boxplot(outlier.colour=NA,alpha=0.1)+scale_x_discrete(limits=c("Tadpole food 1","Tadpole food 2"), labels=c("Diet 1","Diet 2"))+theme_nb()
ggsave("Plot_stage_1.pdf",width=10,height=13,units="cm")
# OTU richness boxplot including significance
set.seed(623)
ggstatsplot::ggbetweenstats(data=df,x=Sample,y=Nb.OTU,type="np",plot.type="box")
## Note: 95% CI for effect size estimate was computed with 100 bootstrap samples.
##
## Note: Shapiro-Wilk Normality Test for Nb.OTU: p-value = 0.104
##
## Note: Bartlett's test for homogeneity of variances for factor Sample: p-value = 0.396
##
ggsave("Plot_stage_1_stats.pdf",width=10,height=13,units="cm")
#Analysis of phylogenetic diversity
comm.pd <- pd(comm, phy)
# For the boxplot with ggplot2
df2<-data.frame(Faith.PD= comm.pd$PD,Sample=metadata$Class)
# Faith's PD boxplot
ggplot(df2,aes(x=Sample,y=Faith.PD)) +geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey",alpha=0.9)+geom_boxplot(outlier.colour=NA,alpha=0.1)+scale_x_discrete(limits=c("Tadpole food 1","Tadpole food 2"), labels=c("Diet 1","Diet 2"))+theme_nb()
ggsave("Plot_stage_2.pdf",width=10,height=13,units="cm")
# Faith's PD boxplot including stats
ggstatsplot::ggbetweenstats(data=df2,x=Sample,y=Faith.PD,type="np",plot.type="box")
## Note: 95% CI for effect size estimate was computed with 100 bootstrap samples.
##
## Note: Shapiro-Wilk Normality Test for Faith.PD: p-value = 0.307
##
## Note: Bartlett's test for homogeneity of variances for factor Sample: p-value = 0.247
##
ggsave("Plot_stage_2_stats.pdf",width=10,height=13,units="cm")
# OTU richness and Faith's PD box plots side by side
plot1<-ggplot(df,aes(x=Sample,y=Nb.OTU))+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey",alpha=0.9)+geom_boxplot(outlier.colour=NA,alpha=0.1)+scale_x_discrete(limits=c("Tadpole food 1","Tadpole food 2"), labels=c("Diet 1","Diet 2"))+theme_nb()+labs(y="OTU richness",x="")
plot2<-ggplot(df2,aes(x=Sample,y=Faith.PD))+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey",alpha=0.9)+geom_boxplot(outlier.colour=NA,alpha=0.1)+scale_x_discrete(limits=c("Tadpole food 1","Tadpole food 2"), labels=c("Diet 1","Diet 2"))+theme_nb()+labs(y="Phylogenetic diversity",x="")
#Saving plot side to side in a file
#pdf("plot12_stage.pdf",width=10,height=8)
#print(grid.arrange(plot1, plot2, ncol=2))
#dev.off()
# Diversity using Shannon index
comm.H<-diversity(comm, index="shannon", base=exp(1))
dfH<-data.frame(H= comm.H,Sample=metadata$Class)
plotH<-ggplot(dfH,aes(x=Sample,y=H))+scale_x_discrete(limits=c("Tadpole food 1","Tadpole food 2"), labels=c("Diet 1","Diet 2"))+geom_jitter(width=0.1,height=0.01,size=2.5,colour="darkgrey")+geom_boxplot(outlier.colour=NA,alpha=0.1)+theme_nb()+labs(y="Shannon diversity",x="")
plotH
# Shannon index boxplot including stats
set.seed(623)
ggstatsplot::ggbetweenstats(data=dfH,x=Sample,y=H,type="np", mean.plotting = FALSE,ylab="Shannon diversity",xlab="",plot.type="box")
## Note: 95% CI for effect size estimate was computed with 100 bootstrap samples.
##
## Note: Shapiro-Wilk Normality Test for Shannon diversity: p-value = 0.049
##
## Note: Bartlett's test for homogeneity of variances for factor : p-value = 0.219
##
ggsave("plot_class_H_stats.pdf",width=6,height=6)
# Diversity using Simpson index
comm.D<-diversity(comm, index="simpson")
dfD<-data.frame(D= comm.D,Sample=metadata$Class)
plotD<-ggplot(dfD,aes(x=Sample,y=D))+scale_x_discrete(limits=c("Tadpole food 1","Tadpole food 2"), labels=c("Diet 1","Diet 2"))+geom_jitter(width=0.1,height=0.01,size=2.5,colour="darkgrey")+geom_boxplot(outlier.colour=NA,alpha=0.1)+theme_nb()+labs(y="Simpson diversity",x="")
plotD
# Simpson index boxplot including stats
set.seed(623)
ggstatsplot::ggbetweenstats(data=dfD,x=Sample,y=D,type="np",mean.plotting = FALSE,ylab="Simpson diversity",xlab="",plot.type="box")
## Note: 95% CI for effect size estimate was computed with 100 bootstrap samples.
##
## Note: Shapiro-Wilk Normality Test for Simpson diversity: p-value = 0.008
##
## Note: Bartlett's test for homogeneity of variances for factor : p-value = 0.008
##
ggsave("plot_class_I_stats.pdf",width=6,height=6)
# Diversity using Pielou equitability index
comm.J <- comm.H/log(specnumber(comm), base=exp(1))
dfJ<-data.frame(J= comm.J,Sample=metadata$Class)
plotJ<-ggplot(dfJ,aes(x=Sample,y=J))+scale_x_discrete(limits=c("Tadpole food 1","Tadpole food 2"), labels=c("Diet 1","Diet 2"))+geom_jitter(width=0.1,height=0.01,size=2.5,colour="darkgrey")+geom_boxplot(outlier.colour=NA,alpha=0.1)+theme_nb()+labs(y="Pielou evenness",x="")
plotJ
# Pielou equitability index boxplot including stats
set.seed(623)
ggstatsplot::ggbetweenstats(data=dfJ,x=Sample,y=J,type="np",mean.plotting = FALSE,ylab="Pielou equitability",xlab="",plot.type="box")
## Note: 95% CI for effect size estimate was computed with 100 bootstrap samples.
##
## Note: Shapiro-Wilk Normality Test for Pielou equitability: p-value = 0.023
##
## Note: Bartlett's test for homogeneity of variances for factor : p-value = 0.265
##
ggsave("plot_class_J_stats.pdf",width=6,height=6)
# Saving boxplots of the three alpha diversity metrics to a file
#pdf("plotHDJ_Stage.pdf",width=6,height=10)
#require(gridExtra)
#print(grid.arrange(plotH, plotD, plotJ,nrow=3,ncol=1))
#dev.off()
#Comparison of dissimilarity measures
rankindex(comm,comm,c("euc","man","bray","jac","kul"))
## euc man bray jac kul
## 0.3344390 0.5566188 0.5566188 0.5566188 0.5566188
# Checking that we have a correct sampling of the diversity in our samples
plot(specaccum(comm), xlab = "# of samples", ylab = "# of species",main="Species accumulation in tadpole's gut microbiomes fed with two different diets" )
# Saving plot to a file
#pdf("plot_spec_accum_stage.pdf",width=7,height=5)
#print(plot(specaccum(comm), xlab = "# of samples", ylab = "# of species",main="Species accumulation in tadpole's gut microbiomes fed with two different diets" ))
#dev.off()
#Significance testing of the diversity measures using Kruskal-Wallis test
kw_specnumber<-kruskal.test(specnumber(comm)~metadata$Class, data= comm)
kw_FaithPD<-kruskal.test(comm.pd$PD~metadata$Class, data= comm)
kw_Shannon<-kruskal.test(comm.H~metadata$Class, data= comm)
kw_Simpson<-kruskal.test(comm.D~metadata$Class, data= comm)
kw_Evenness<-kruskal.test(comm.J~metadata$Class, data= comm)
# Computing dissimilarity using Bray-Curtis distance (weighted)
comm.bc.dist <- vegdist(comm, method = "bray")
# Clustering to analyse samples - Bray
comm.bc.clust <- hclust(comm.bc.dist, method = "average")
plot(comm.bc.clust, ylab = "Bray-Curtis dissimilarity")
#Saving plot to a file
#pdf("plot_bc_clust_food.pdf",width=8,height=6)
#print(plot(comm.bc.clust, ylab = "Bray-Curtis dissimilarity"))
#dev.off()
#Dissimilarity using Sorenson (unweighted)
comm.so.dist<-vegdist(comm,method="bray",binary="TRUE")
# Clustering to analyse samples - Sorenson
comm.so.clust <- hclust(comm.so.dist, method = "average")
plot(comm.so.clust, ylab = "Sorenson dissimilarity")
#pdf("plot_so_clust_food.pdf",width=8,height=6)
#print(plot(comm.so.clust, ylab = "Sorenson dissimilarity"))
#dev.off()
################################
# Ordination NMDS Bray-Curtis #
################################
# Ordination using NMDS
comm.bc.mds <- metaMDS(comm, dist = "bray")
## Run 0 stress 0.07795151
## Run 1 stress 0.07795152
## ... Procrustes: rmse 2.432252e-05 max resid 5.512452e-05
## ... Similar to previous best
## Run 2 stress 0.07795156
## ... Procrustes: rmse 0.0001076706 max resid 0.0002985161
## ... Similar to previous best
## Run 3 stress 0.1130982
## Run 4 stress 0.1095625
## Run 5 stress 0.07704883
## ... New best solution
## ... Procrustes: rmse 0.07374961 max resid 0.2616592
## Run 6 stress 0.0778795
## Run 7 stress 0.09829115
## Run 8 stress 0.07787946
## Run 9 stress 0.1106998
## Run 10 stress 0.07787941
## Run 11 stress 0.07787967
## Run 12 stress 0.118126
## Run 13 stress 0.07704825
## ... New best solution
## ... Procrustes: rmse 0.0002850224 max resid 0.0005449019
## ... Similar to previous best
## Run 14 stress 0.1186749
## Run 15 stress 0.0770485
## ... Procrustes: rmse 0.0002471256 max resid 0.0005003866
## ... Similar to previous best
## Run 16 stress 0.07704831
## ... Procrustes: rmse 0.0001404113 max resid 0.0002867942
## ... Similar to previous best
## Run 17 stress 0.07704827
## ... Procrustes: rmse 5.610118e-05 max resid 0.000125906
## ... Similar to previous best
## Run 18 stress 0.09510528
## Run 19 stress 0.0982972
## Run 20 stress 0.07704831
## ... Procrustes: rmse 0.0001204728 max resid 0.0002470797
## ... Similar to previous best
## *** Solution reached
# We have a stress value of 0.08, thus the NMDS ordination provides a good representation of the distances between the communities of the tadpole's gut samples.
# Checking the NMDS
# Saving plot to a file
stressplot(comm.bc.mds)
#pdf("plot_nmds_stressplot.pdf",width=8,height=6)
#print(stressplot(comm.bc.mds))
#dev.off()
mds.fig <- ordiplot(comm.bc.mds, type = "none")
points(mds.fig, "sites", pch = 19, col = "grey")
orditorp(comm.bc.mds, display = "sites")
# Saving plot to a file
#pdf("plot1_nmds.pdf",width=8,height=6)
#print(mds.fig <- ordiplot(comm.bc.mds, type = "none"))
#print(points(mds.fig, "sites", pch = 19, col = "grey"))
#print(orditorp(comm.bc.mds, display = "sites"))
#dev.off()
# Just another plot of the same NMDS ordination
mds.fig <- ordiplot(comm.bc.mds, type = "none")
points(mds.fig, "sites", pch = 19, col = "green", select = metadata$Class == "Tadpole food 1")
points(mds.fig, "sites", pch = 19, col = "black", select = metadata$Class == "Tadpole food 2")
ordiellipse(comm.bc.mds, metadata$Class, conf = 0.95, label = FALSE)
ordicluster(comm.bc.mds, comm.so.clust, col = "gray")
# Saving plot to file
#pdf("plot2_nmds.pdf",width=8,height=6)
#print(mds.fig <- ordiplot(comm.bc.mds, type = "none"))
#print(points(mds.fig, "sites", pch = 19, col = "green", select = metadata$Class == "Tadpole food 1"))
#print(points(mds.fig, "sites", pch = 19, col = "black", select = metadata$Class == "Tadpole food 2"))
#print(ordiellipse(comm.bc.mds, metadata$Class, conf = 0.95, label = FALSE))
#print(ordicluster(comm.bc.mds, comm.so.clust, col = "gray"))
#dev.off()
# Plotting NMDS ordination using Bray-Curtis distance
#par(new=TRUE)
mds.fig2 <- ordiplot(comm.bc.mds, type = "none")
points(mds.fig2, "species", pch = 1,cex=0.5, col = "grey")
ordispider(comm.bc.mds,groups=metadata$Class,show.groups="Tadpole food 1",col="green")
points(mds.fig2, "sites", pch = 19, col = "green", select = metadata$Class == "Tadpole food 1")
ordispider(comm.bc.mds,groups=metadata$Class,show.groups="Tadpole food 2",col="black")
points(mds.fig2, "sites", pch = 19, col = "black", select = metadata$Class == "Tadpole food 2")
ordiellipse(comm.bc.mds, metadata$Class, conf = 0.95, label = FALSE,lty = 'dotted')
# Saving plot to a file
#pdf("plot3_nmds.pdf",width=8,height=6)
#print(mds.fig2 <- ordiplot(comm.bc.mds, type = "none",xlim=c(-2,2),ylim=c(-1.5,1.5)))
#print(points(mds.fig2, "species", pch = 1,cex=0.5, col = "grey"))
#print(ordispider(comm.bc.mds,groups=metadata$Class,show.groups="Tadpole food 1",col="green"))
#print(points(mds.fig2, "sites", pch = 19, col = "green", select = metadata$Class == "Tadpole food 1"))
#print(ordispider(comm.bc.mds,groups=metadata$Class,show.groups="Tadpole food 2",col="black"))
#print(points(mds.fig2, "sites", pch = 19, col = "black", select = metadata$Class == "Tadpole food 2"))
#print(ordiellipse(comm.bc.mds, metadata$Class, conf = 0.95, label = FALSE,lty = 'dotted'))
#dev.off()
################################
# Ordination beta-dispersion #
################################
# Plotting the ordination results of beta-dispersion PCoA (PERMDISP)
comm.betadisp=betadisper(comm.bc.dist, metadata$Class, type="median")
comm.betadisp
##
## Homogeneity of multivariate dispersions
##
## Call: betadisper(d = comm.bc.dist, group = metadata$Class, type =
## "median")
##
## No. of Positive Eigenvalues: 14
## No. of Negative Eigenvalues: 3
##
## Average distance to median:
## Tadpole food 1 Tadpole food 2
## 0.3588 0.1979
##
## Eigenvalues for PCoA axes:
## (Showing 8 of 17 eigenvalues)
## PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
## 1.82429 0.53684 0.19737 0.11899 0.07153 0.05827 0.04529 0.03039
mds.fig3a<-plot(comm.betadisp)
#Saving plot to a file
#pdf("plot1_betadispa.pdf",width=8,height=6)
#print(mds.fig3a<-plot(comm.betadisp))
#dev.off()
# Another plot of PCoA for beta dispersion
#plot.new()
#pdf("plot1_betadisp.pdf",width=8,height=6)
#print(mds.fig3 <- ordiplot(comm.betadisp, type = "none",xlim=c(-0.7,0.7),ylim=c(-0.7,0.7)))
#print(ordispider(comm.betadisp,groups=metadata$Class,show.groups="Tadpole food 1",col="green"),xlim=c(-0.7,0.7),ylim=c(-0.7,0.7))
#print(points(mds.fig3, "sites", pch = 19, col = "green", select = metadata$Class == "Tadpole food 1"))
#print(ordispider(comm.betadisp,groups=metadata$Class,show.groups="Tadpole food 2",col="black"),xlim=c(-0.7,0.7),ylim=c(-0.7,0.7))
#print(points(mds.fig3, "sites", pch = 19, col = "black", select = metadata$Class == "Tadpole food 2"))
#print(ordiellipse(comm.betadisp, metadata$Class, conf = 0.95, label = FALSE,lty = 'dotted'),xlim=c(-0.7,0.7),ylim=c(-0.7,0.7))
#dev.off()
# Significance
comm.betadisp.perm=permutest(comm.betadisp, group= metadata$Stage, type="median", permutations=999, pairwise=TRUE)
comm.betadisp.perm
##
## 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 1 0.11661 0.11661 13.419 999 0.001 ***
## Residuals 16 0.13904 0.00869
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Tadpole food 1 Tadpole food 2
## Tadpole food 1 0.001
## Tadpole food 2 0.0020995
#######################################################
# Analysis of community structure using MDP and MNTD #
#######################################################
#Phylogenetic distance to measure the community phylogenetic structure without abundance
phy.dist <- cophenetic(phy)
comm.sesmpd <- ses.mpd(comm, phy.dist, null.model = "richness", abundance.weighted = FALSE, runs = 999)
comm.sesmntd <- ses.mntd(comm, phy.dist, null.model = "richness", abundance.weighted = FALSE, runs = 999)
#Inspecting obtained values
#comm.sesmpd
#comm.sesmntd
#Figure with comparisons SES(MPD) and SES(MNTD)
df3<-data.frame(SES_MPD = comm.sesmpd$mpd.obs.z, Sample = metadata$Class)
plot3<-ggplot(df3,aes(x=Sample,y=SES_MPD))+geom_boxplot()+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey",alpha=0.9)+scale_x_discrete(limits=c("Tadpole food 1","Tadpole food 2"), labels=c("Diet 1","Diet 2"))+theme_nb()+labs(x="",y="SES MPD")
df4<-data.frame(SES_MNTD = comm.sesmntd$mntd.obs.z, Sample = metadata$Class)
plot4<-ggplot(df4,aes(x=Sample,y=SES_MNTD))+geom_boxplot()+geom_jitter(width=0.1,height=0.01,size=5,colour="darkgrey",alpha=0.9)+scale_x_discrete(limits=c("Tadpole food 1","Tadpole food 2"), labels=c("Diet 1","Diet 2"))+theme_nb()+labs(x="",y="SES MNTD")
grid.arrange(plot3, plot4, ncol=2)
# Saving plot to a file
#pdf("plot34_sesmpdmntd.pdf",width=5,height=4)
#print(grid.arrange(plot3, plot4, ncol=2))
#dev.off()
#Significance
#Using Kruskal-Wallis test
kw_sesmpd<-kruskal.test(comm.sesmpd$mpd.obs.z ~ metadata$Class, data= comm)
kw_sesmntd<-kruskal.test(comm.sesmntd$mntd.obs.z ~ metadata$Class, data= comm)
set.seed(623)
ggstatsplot::ggbetweenstats(data=df3,x=Sample,y=SES_MPD,mean.plotting = FALSE,type="np",plot.type = "box")
## Note: 95% CI for effect size estimate was computed with 100 bootstrap samples.
##
## Note: Shapiro-Wilk Normality Test for SES_MPD: p-value = 0.434
##
## Note: Bartlett's test for homogeneity of variances for factor Sample: p-value = 0.105
##
ggsave("plot_SESMPD_stat.pdf",width=5,height=4)
ggstatsplot::ggbetweenstats(data=df4,x=Sample,y=SES_MNTD,mean.plotting = FALSE,type="np",plot.type = "box")
## Note: 95% CI for effect size estimate was computed with 100 bootstrap samples.
##
## Note: Shapiro-Wilk Normality Test for SES_MNTD: p-value = 0.817
##
## Note: Bartlett's test for homogeneity of variances for factor Sample: p-value = 0.335
##
ggsave("plot_SESMNTD_stat.pdf",width=5,height=4)
#####################################
# Phylogenetic Beta-diversity #
#####################################
# beta diversity: Taxonomic (Bray-Curtis) dissimilarity explained (Permanova)
adonis(comm.bc.dist ~ Class, data = metadata)
##
## Call:
## adonis(formula = comm.bc.dist ~ Class, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Class 1 1.2887 1.28873 12.638 0.44131 0.001 ***
## Residuals 16 1.6315 0.10197 0.55869
## Total 17 2.9203 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# beta diversity: compute phylogenetic dissimilarity using MPD
comm.mpd.dist<-comdist(comm,phy.dist,abundance.weighted=TRUE)
adonis(comm.mpd.dist ~ Class, data = metadata)
##
## Call:
## adonis(formula = comm.mpd.dist ~ Class, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Class 1 0.4981 0.49814 2.6887 0.14387 0.002 **
## Residuals 16 2.9644 0.18528 0.85613
## Total 17 3.4626 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# beta diversity: compute phylogenetic dissimilarity using MNTD
comm.mntd.dist <- comdistnt(comm, phy.dist, abundance.weighted = TRUE)
adonis(comm.mntd.dist ~ Class, data = metadata)
##
## Call:
## adonis(formula = comm.mntd.dist ~ Class, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Class 1 1.6146e-04 1.6146e-04 29.135 0.64551 0.001 ***
## Residuals 16 8.8671e-05 5.5420e-06 0.35449
## Total 17 2.5013e-04 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# calculate Mantel correlation for taxonomic Bray-Curtis vs. phylogenetic
# MPD diversity
mantel(comm.bc.dist, comm.mpd.dist)
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = comm.bc.dist, ydis = comm.mpd.dist)
##
## Mantel statistic r: 0.8879
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.109 0.159 0.211 0.278
## Permutation: free
## Number of permutations: 999
#MNTD diversity
mantel(comm.bc.dist, comm.mntd.dist)
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = comm.bc.dist, ydis = comm.mntd.dist)
##
## Mantel statistic r: 0.6867
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.111 0.160 0.197 0.224
## Permutation: free
## Number of permutations: 999