#Loading packages
library(picante)
library(ggplot2)
library(phytools)
library(dplyr)
library(grid)
library(PMCMRplus)
library(gridExtra)
library(multcomp)
library(multcompView)
## Folders, Themes, colors
source("prelude.R")
# Setting up directories
data_dir <-paste0(data_dir_path,"xptransmission")
output_dir <- paste0(output_dir_path,"Figure6_transmission")
### Rarefaction analysis : importing data
setwd(data_dir)
comm_abundance_1<-read.table("nmt3_abondance_tpn_notax_xp_transa.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 a pdf file
#pdf("plot_rarefaction1_trans.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 a pdf file
#pdf("plot_rarefaction2_trans.pdf",width=5,height=4)
#print(rarecurve(comm_abundance_1, step = 20, sample = raremax, col = "blue", cex = 0.6))
#dev.off()
###Diversity analysis
# Importing data
# Defining the directory containing the data to import
setwd(data_dir)
# Working on abundance data from 100 rarefactions
allcom<-read.table("xptrans_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
#Importating metadata
metadata<-read.csv("nmt3_metadata_XP_transa.csv",sep=";",header=TRUE,row.names=1)
metadata$XP_Tissue<-as.factor(metadata$XP_Tissue)
#Importating 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)
## CRT1E CRT1FFE CRT1FP CRT1MFE CRT1MP CRT2E CRT2FFE CRT2FP CRT2MFE CRT2MP
## 1 1 1 1 1 1 1 1 1 1
## CRT3E CRT3FFE CRT3FP CRT3MFE CRT3MP CRT4E CRT4FFE CRT4FP CRT4MFE CRT4MP
## 1 1 1 1 1 1 1 1 1 1
## TW2FFE TW2FP TW2MFE TW2MP XTTW2E TW3FFE TW3FP TW3MP XTTW3E TW4FP
## 1 1 1 1 1 1 1 1 1 1
## TW4MFE TW4MP XTTW4E TW5FP TW5MFE TW5MP XTTW5E
## 1 1 1 1 1 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$XP_Tissue,Cross=metadata$Cross)
# 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("Eggs","Skin","Feces"), labels=c("Eggs","Skin","Feces"))+theme_npgray()
#ggsave("Plot_transm_1.pdf",width=10,height=8,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$XP_Tissue)
# 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("Eggs","Skin","Feces"), labels=c("Eggs","Skin","Feces"))+theme_npgray()
#ggsave("Plot_transm_2.pdf",width=10,height=8,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("Eggs","Skin","Feces"))+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("Eggs","Skin","Feces"), labels=c("Eggs","Skin","Feces"))+theme_nb()+labs(y="Phylogenetic diversity",x="")
# Saving plot to a pdf file
#pdf("plot12_transm.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$XP_Tissue)
plotH<-ggplot(dfH,aes(x=Sample,y=H))+scale_x_discrete(limits=c("Eggs","Skin","Feces"), labels=c("Eggs","Skin","Feces"))+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
comm.D<-diversity(comm, index="simpson")
dfD<-data.frame(D=comm.D,Sample=metadata$XP_Tissue)
plotD<-ggplot(dfD,aes(x=Sample,y=D))+scale_x_discrete(limits=c("Eggs","Skin","Feces"), labels=c("Eggs","Skin","Feces"))+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
# Diversity using Pielou equitability index
comm.J<-comm.H/log(specnumber(comm), base=exp(1))
dfJ<-data.frame(J= comm.J,Sample=metadata$XP_Tissue)
plotJ<-ggplot(dfJ,aes(x=Sample,y=J))+scale_x_discrete(limits=c("Eggs","Skin","Feces"), labels=c("Eggs","Skin","Feces"))+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
# Boxplots of the alpha diversity metrics
# grid.arrange(plotH, plotD, plotJ,nrow=3,ncol=1)
# Saving plot to a pdf file
#pdf("plotHDJ_Stage.pdf",width=6,height=10)
#print(grid.arrange(plotH, plotD, plotJ,nrow=3,ncol=1))
#dev.off()
#Significance testing of the diversity measures using Kruskal-Wallis test
# We select Kruskal-Wallis test because we have less than 30 samples and thus a potential problem of normality and/or homogeneity
kw_specnumber<-kruskal.test(specnumber(comm)~metadata$XP_Tissue, data= comm)
kw_FaithPD<-kruskal.test(comm.pd$PD~metadata$XP_Tissue, data= comm)
kw_Shannon<-kruskal.test(comm.H~metadata$XP_Tissue, data= comm)
kw_Evenness<-kruskal.test(comm.J~metadata$XP_Tissue, data= comm)
#Pairwise comparisons using Dunn-Bonferroni test
pp_spec<-kwAllPairsDunnTest(x=specnumber(comm), metadata$XP_Tissue,p.adjust.method="bonferroni")
pp_pd<-kwAllPairsDunnTest(x=comm.pd$PD, metadata$XP_Tissue,p.adjust.method="bonferroni")
pp_shannon<-kwAllPairsDunnTest(x=comm.H, metadata$XP_Tissue,p.adjust.method="bonferroni")
pp_Evenness<-kwAllPairsDunnTest(x=comm.J, metadata$XP_Tissue,p.adjust.method="bonferroni")
# Making letters corresponding to the significance level
mymat_spec <-tri.to.squ(pp_spec$p.value)
mymat_pd <-tri.to.squ(pp_pd$p.value)
myletters_spec <- multcompLetters(mymat_spec,compare="<=",threshold=0.05,Letters=letters)
myletters_pd <- multcompLetters(mymat_pd,compare="<=",threshold=0.05,Letters=letters)
# Store the letters in a dataframe
myletters_spec_df <- data.frame(Sample=names(myletters_spec$Letters),letter = myletters_spec$Letters )
myletters_pd_df <- data.frame(Sample=names(myletters_pd$Letters),letter = myletters_pd$Letters )
#Plot for Figure
plot_div_1<-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("Eggs","Skin","Feces"),labels=c("Eggs","Skin","Feces"))+theme_nb()+xlab("Samples")+geom_text(data = myletters_spec_df, aes(label = letter, y = 390 ), colour="black", size=5,fontface="bold")
plot_div_2<-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("Eggs","Skin","Feces"), labels=c("Eggs","Skin","Feces"))+theme_nb()+xlab("Samples")+geom_text(data = myletters_spec_df, aes(label = letter, y = 20.5 ), colour="black", size=5,fontface="bold")
grid.arrange(plot_div_1, plot_div_2, ncol=2)
# Saving plot as a pdf file
#pdf("plot_div_12_transm.pdf",width=10,height=8)
#print(grid.arrange(plot_div_1, plot_div_2, ncol=2))
#dev.off()
#Comparison of dissimilarity measures
rankindex(comm,comm,c("euc","man","bray","jac","kul"))
## euc man bray jac kul
## -0.3040496 0.4716522 0.4716522 0.4716522 0.4716522
# 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 eggs, skin and feces")
# Saving plot as a pdf file
#pdf("plot_spec_accum_transm.pdf",width=7,height=5)
#print(plot(specaccum(comm), xlab = "# of samples", ylab = "# of species",main="Species accumulation in eggs, skin and feces"))
#dev.off()
#Computing dissimilarity using Bray-Curtis (weighted)
comm.bc.dist <- vegdist(comm, method = "bray")
#Clustering to analyse the samples - Bray
comm.bc.clust <- hclust(comm.bc.dist, method = "average")
plot(comm.bc.clust, ylab = "Bray-Curtis dissimilarity")
# Saving plot as a pdf file
#pdf("plot_bc_clust_transm.pdf",width=8,height=6)
#print(plot(comm.bc.clust, ylab = "Bray-Curtis dissimilarity"))
#dev.off()
#Computing dissimilarity using Sorenson (unweighted)
comm.so.dist<-vegdist(comm,method="bray",binary="TRUE")
#Clustering to analyse the samples - Sorenson
comm.so.clust <- hclust(comm.so.dist, method = "average")
plot(comm.so.clust, ylab = "Sorenson dissimilarity")
# Saving plot as a pdf file
#pdf("plot_so_clust_transm.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.1315412
## Run 1 stress 0.138584
## Run 2 stress 0.1385899
## Run 3 stress 0.1408388
## Run 4 stress 0.1315412
## ... Procrustes: rmse 1.250927e-05 max resid 4.078501e-05
## ... Similar to previous best
## Run 5 stress 0.1470581
## Run 6 stress 0.1315413
## ... Procrustes: rmse 3.994118e-05 max resid 0.000123488
## ... Similar to previous best
## Run 7 stress 0.1727038
## Run 8 stress 0.1315412
## ... Procrustes: rmse 2.913916e-05 max resid 0.000136542
## ... Similar to previous best
## Run 9 stress 0.1315412
## ... Procrustes: rmse 8.551574e-06 max resid 2.96676e-05
## ... Similar to previous best
## Run 10 stress 0.131639
## ... Procrustes: rmse 0.00353371 max resid 0.01744158
## Run 11 stress 0.1316389
## ... Procrustes: rmse 0.003527948 max resid 0.01746212
## Run 12 stress 0.1315412
## ... Procrustes: rmse 6.965893e-06 max resid 1.987238e-05
## ... Similar to previous best
## Run 13 stress 0.1408385
## Run 14 stress 0.1415515
## Run 15 stress 0.1315412
## ... Procrustes: rmse 9.90658e-06 max resid 2.538888e-05
## ... Similar to previous best
## Run 16 stress 0.1315414
## ... Procrustes: rmse 0.0001173342 max resid 0.0005387483
## ... Similar to previous best
## Run 17 stress 0.1315412
## ... Procrustes: rmse 4.418682e-05 max resid 0.0001644307
## ... Similar to previous best
## Run 18 stress 0.1315414
## ... Procrustes: rmse 4.643509e-05 max resid 0.0001627682
## ... Similar to previous best
## Run 19 stress 0.3958534
## Run 20 stress 0.1316389
## ... Procrustes: rmse 0.0035248 max resid 0.01745808
## *** Solution reached
#Checking the NMDS
# Saving plot as a pdf file
stressplot(comm.bc.mds)
#pdf("plot_nmds_stressplot.pdf",width=8,height=6)
#print(stressplot(comm.bc.mds))
#dev.off()
# We plot the NMDS ordination using Bray-Curtis
#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$XP_Tissue,show.groups="Eggs",col="green")
points(mds.fig2, "sites", pch = 19, col = "green", select = metadata$XP_Tissue == "Eggs")
ordispider(comm.bc.mds,groups=metadata$XP_Tissue,show.groups="Skin",col="blue")
points(mds.fig2, "sites", pch = 19, col = "blue", select = metadata$XP_Tissue == "Skin")
ordispider(comm.bc.mds,groups=metadata$XP_Tissue,show.groups="Feces",col="black")
points(mds.fig2, "sites", pch = 19, col = "black", select = metadata$XP_Tissue == "Feces")
ordiellipse(comm.bc.mds, metadata$XP_Tissue, conf = 0.95, label = FALSE,lty = 'dotted')
# Saving plot as a pdf file
#pdf("plot3_nmds.pdf",width=8,height=6)
#print(mds.fig2 <- ordiplot(comm.bc.mds, type = "none",xlim=c(-2,2),ylim=c(-2,2)))
#print(points(mds.fig2, "species", pch = 1,cex=0.5, col = "grey"))
#print(ordispider(comm.bc.mds,groups=metadata$XP_Tissue,show.groups="Eggs",col="green"))
#print(points(mds.fig2, "sites", pch = 19, col = "green", select = metadata$XP_Tissue == "Eggs"))
#print(ordispider(comm.bc.mds,groups=metadata$XP_Tissue,show.groups="Skin",col="blue"))
#print(points(mds.fig2, "sites", pch = 19, col = "blue", select = metadata$XP_Tissue == "Skin"))
#print(ordispider(comm.bc.mds,groups=metadata$XP_Tissue,show.groups="Feces",col="black"))
#print(points(mds.fig2, "sites", pch = 19, col = "black", select = metadata$XP_Tissue == "Feces"))
#print(ordiellipse(comm.bc.mds, metadata$XP_Tissue, conf = 0.95, label = FALSE,lty = 'dotted'))
#dev.off()
################################
# Ordination beta-dispersion #
################################
# PCoA of beta-dispersion (PERMDISP)
comm.betadisp=betadisper(comm.bc.dist, metadata$XP_Tissue, type="median")
comm.betadisp
##
## Homogeneity of multivariate dispersions
##
## Call: betadisper(d = comm.bc.dist, group = metadata$XP_Tissue, type =
## "median")
##
## No. of Positive Eigenvalues: 31
## No. of Negative Eigenvalues: 5
##
## Average distance to median:
## Eggs Feces Skin
## 0.4401 0.4609 0.4781
##
## Eigenvalues for PCoA axes:
## (Showing 8 of 36 eigenvalues)
## PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
## 3.7034 1.3696 1.1709 0.7754 0.7448 0.6284 0.5064 0.4412
mds.fig3a<-plot(comm.betadisp)
#pdf("plot1_betadispa.pdf",width=8,height=6)
#print(mds.fig3a<-plot(comm.betadisp))
#dev.off()
# Another Figure of PCoA for beta dispersion
#pdf("plot1_betadisp.pdf",width=8,height=6)
#print(mds.fig3 <- ordiplot(comm.betadisp, type = "none"))
#print(ordispider(comm.betadisp,groups=metadata$XP_Tissue,show.groups="Eggs",col="green"))
#print(points(mds.fig3, "sites", pch = 19, col = "green", select = metadata$XP_Tissue == "Eggs"))
#print(ordispider(comm.betadisp,groups=metadata$XP_Tissue,show.groups="Skin",col="blue"))
#print(points(mds.fig3, "sites", pch = 19, col = "blue", select = metadata$XP_Tissue == "Skin"))
#print(ordispider(comm.betadisp,groups=metadata$XP_Tissue,show.groups="Feces",col="black"))
#print(points(mds.fig3, "sites", pch = 19, col = "black", select = metadata$XP_Tissue == "Feces"))
#print(ordiellipse(comm.betadisp, metadata$XP_Tissue, conf = 0.95, label = FALSE,lty = 'dotted'))
#dev.off()
############################################
# Community structure using MPD 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 = 9)
comm.sesmntd <- ses.mntd(comm, phy.dist, null.model = "richness", abundance.weighted = FALSE, runs = 9)
#Looking at the values
#comm.sesmpd
#comm.sesmntd
# Significance
#Using Kruskal-Wallis test
kw_sesmpd<-kruskal.test(comm.sesmpd$mpd.obs.z ~ metadata$XP_Tissue, data= comm)
kw_sesmntd<-kruskal.test(comm.sesmntd$mntd.obs.z ~ metadata$XP_Tissue, data= comm)
#Pairwise comparisons using Dunn-Bonferroni test implemented in PMCMRplus
pp_sesmpd<-kwAllPairsDunnTest(x=comm.sesmpd$mpd.obs.z, metadata$XP_Tissue,p.adjust.method="bonferroni")
pp_sesmntd<-kwAllPairsDunnTest(x=comm.sesmntd$mntd.obs.z, metadata$XP_Tissue,p.adjust.method="bonferroni")
mymat_sesmpd <-tri.to.squ(pp_sesmpd$p.value)
mymat_sesmntd <-tri.to.squ(pp_sesmntd$p.value)
# Making letters corresponding to the significance level
myletters_sesmpd <- multcompLetters(mymat_sesmpd,compare="<=",threshold=0.05,Letters=letters)
myletters_sesmntd <- multcompLetters(mymat_sesmntd,compare="<=",threshold=0.05,Letters=letters)
# Store the letters in a dataframe
myletters_sesmpd_df <- data.frame(Sample=names(myletters_sesmpd$Letters),letter = myletters_sesmpd$Letters )
myletters_sesmntd_df <- data.frame(Sample=names(myletters_sesmntd$Letters),letter = myletters_sesmntd$Letters )
#Preparing data for the graphical representation of SES(MPD) and SES(MNTD)
df3<-data.frame(SES_MPD = comm.sesmpd$mpd.obs.z, Sample = metadata$XP_Tissue)
df4<-data.frame(SES_MNTD = comm.sesmntd$mntd.obs.z, Sample = metadata$XP_Tissue)
# Plotting the graphical comparisons of SES(MPD) and et SES(MNTD) values
plot3<-ggplot(df3,aes(x=Sample,y=SES_MPD))+geom_boxplot()+scale_x_discrete(limits=c("Eggs","Skin","Feces"), labels=c("Eggs","Skin","Feces"))+theme_nb()+labs(x="",y="SES MPD")+geom_text(data = myletters_sesmpd_df, aes(label = letter, y = 5.5 ),colour="black", size=5,fontface="bold")
plot4<-ggplot(df4,aes(x=Sample,y=SES_MNTD))+geom_boxplot()+scale_x_discrete(limits=c("Eggs","Skin","Feces"), labels=c("Eggs","Skin","Feces"))+theme_nb()+labs(x="",y="SES MNTD")+geom_text(data = myletters_sesmntd_df, aes(label = letter, y = 5.5 ), colour="black", size=5,fontface="bold")
grid.arrange(plot3, plot4, ncol=2)
ggsave("plot34_sesmpdmntd.pdf",width=5,height=4)
#####################################
# Phylogenetic Beta-diversity #
#####################################
# beta diversity: Taxonomic (Bray-Curtis) dissimilarity explained (Permanova)
adonis(comm.bc.dist ~ XP_Tissue, data = metadata)
##
## Call:
## adonis(formula = comm.bc.dist ~ XP_Tissue, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## XP_Tissue 2 4.2090 2.10448 8.7918 0.34088 0.001 ***
## Residuals 34 8.1385 0.23937 0.65912
## Total 36 12.3474 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Significance of beta-dispersion
comm.betadisp.perm=permutest(comm.betadisp, group= metadata$XP_Tissue, 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 2 0.007877 0.0039386 0.6945 999 0.524
## Residuals 34 0.192824 0.0056713
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Eggs Feces Skin
## Eggs 0.58300 0.264
## Feces 0.54427 0.566
## Skin 0.24960 0.55114
comm.betadisp.HSD<-TukeyHSD(comm.betadisp)
plot(comm.betadisp.HSD)
# Testing variance homogeneity using anova
anova(comm.betadisp)
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.007877 0.0039386 0.6945 0.5063
## Residuals 34 0.192824 0.0056713
# beta diversity: compute phylogenetic dissimilarity using MPD
comm.mpd.dist<-comdist(comm,phy.dist,abundance.weighted=TRUE)
adonis(comm.mpd.dist ~ XP_Tissue, data = metadata)
##
## Call:
## adonis(formula = comm.mpd.dist ~ XP_Tissue, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## XP_Tissue 2 2.2369 1.1184 9.0711 0.34794 0.001 ***
## Residuals 34 4.1921 0.1233 0.65206
## Total 36 6.4289 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 ~ XP_Tissue, data = metadata)
##
## Call:
## adonis(formula = comm.mntd.dist ~ XP_Tissue, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## XP_Tissue 2 0.011919 0.0059593 88.134 0.8383 0.001 ***
## Residuals 34 0.002299 0.0000676 0.1617
## Total 36 0.014217 1.0000
## ---
## 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.6583
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0693 0.0969 0.1167 0.1428
## 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.6117
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0538 0.0765 0.0957 0.1236
## Permutation: free
## Number of permutations: 999