# Loading packages
library(picante)
## Loading required package: ape
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
## Loading required package: nlme
library(ggplot2)
library(phytools)
## Loading required package: maps
##
## Attaching package: 'phytools'
## The following object is masked from 'package:vegan':
##
## scores
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:nlme':
##
## collapse
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(grid)
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(ggstatsplot)
## Registered S3 methods overwritten by 'broom.mixed':
## method from
## augment.lme broom
## augment.merMod broom
## glance.lme broom
## glance.merMod broom
## glance.stanreg broom
## tidy.brmsfit broom
## tidy.gamlss broom
## tidy.lme broom
## tidy.merMod broom
## tidy.rjags broom
## tidy.stanfit broom
## tidy.stanreg broom
## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
library(PMCMRplus)
library(multcompView)
## Folders, Themes, colors
source("prelude.R")
# Setting up directories
data_dir <-paste0(data_dir_path,"xpgut")
output_dir <- paste0(output_dir_path,"Figure3_gut")
# Defining the directory containing the data to import
setwd(data_dir)
###Rarefaction analysis
comm_abundance_1<-read.table("nmt3_abondance_tpn_notax_XP_AdultTissue.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)
setwd(output_dir)
plot(S, Srare, xlab = "Observed No. of Species", ylab = "Rarefied No. of Species")
abline(0, 1)
rarecurve(comm_abundance_1, step = 20, sample = raremax, col = "blue", cex = 0.6)
#pdf("plot_rarefaction1_gut.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()
#pdf("plot_rarefaction2_gut.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("xpgut_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_AdultTissue.csv",sep=";",header=TRUE,row.names=1)
metadata$Tissue<-as.factor(metadata$Tissue)
#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)
## Lab-col-01 Lab-col-02 Lab-col-03 TW2FFE TW2MFE TW3FFE
## 1 1 1 1 1 1
## TW4MFE TW5MFE CRT1FFE CRT1MFE CRT2FFE CRT2MFE
## 1 1 1 1 1 1
## CRT3FFE CRT3MFE CRT4FFE CRT4MFE Lab-feces-01 Lab-feces-02
## 1 1 1 1 1 1
## Lab-feces-03 XTLSM1 XTLSM2 XTLSM3 XTLTW2M XTLTW3F
## 1 1 1 1 1 1
## XTLTW3M XTLTW4F XTLTW5F Lab-int-01 Lab-int-02 Lab-int-03
## 1 1 1 1 1 1
## Lab-est-01 Lab-est-02 Lab-est-03 TW2FP TW2MP TW3FP
## 1 1 1 1 1 1
## TW3MP TW4FP TW4MP TW5FP TW5MP CRT1FP
## 1 1 1 1 1 1
## CRT1MP CRT2FP CRT2MP CRT3FP CRT3MP CRT4FP
## 1 1 1 1 1 1
## CRT4MP Lab-swab-01 Lab-swab-02
## 1 1 1
#######################################
# Analysis of taxonomic diversity #
#######################################
# Defining the output directory
setwd(output_dir)
#For the boxplot with ggplot2
# OTU richness
df<-data.frame(Nb.OTU= specnumber(comm),Sample=metadata$Tissue)
#Analysis of phylogenetic diversity
comm.pd <- pd(comm, phy)
df2<-data.frame(Faith.PD= comm.pd$PD,Sample=metadata$Tissue)
# 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("Skin","Feces","Colon","Intestine","Gut","Stomach"),labels=c("Ski","Fec","Col","Int","Gut","Sto"))+theme_nb()
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
ggsave("Plot_AdultTissue_1.pdf",width=10,height=13,units="cm")
# 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("Skin","Feces","Colon","Intestine","Gut","Stomach"),labels=c("Ski","Fec","Col","Int","Gut","Sto"))+theme_nb()
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
ggsave("Plot_AdultTissue_2.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("Skin","Feces","Colon","Intestine","Gut","Stomach"),labels=c("Ski","Fec","Col","Int","Gut","Sto"))+theme_nb()+labs(y="OTU richness",x="")
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
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("Skin","Feces","Colon","Intestine","Gut","Stomach"),labels=c("Ski","Fec","Col","Int","Gut","Sto"))+theme_nb()+labs(y="Phylogenetic diversity",x="")
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
grid.arrange(plot1, plot2, ncol=2)
#pdf("plot12_AdultTissue.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$Tissue)
plotH<-ggplot(dfH,aes(x=Sample,y=H))+scale_x_discrete(limits=c("Skin","Feces","Colon","Intestine","Gut","Stomach"),labels=c("Ski","Fec","Col","Int","Gut","Sto"))+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="")
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
plotH
# Diversity using Simpson index
comm.D<-diversity(comm, index="simpson")
dfD<-data.frame(D= comm.D,Sample=metadata$Tissue)
plotD<-ggplot(dfD,aes(x=Sample,y=D))+scale_x_discrete(limits=c("Skin","Feces","Colon","Intestine","Gut","Stomach"),labels=c("Ski","Fec","Col","Int","Gut","Sto"))+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="")
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
plotD
# Diversity using Pielou equitability index
comm.J<-comm.H/log(specnumber(comm), base=exp(1))
dfJ<-data.frame(J= comm.J,Sample=metadata$Tissue)
plotJ<-ggplot(dfJ,aes(x=Sample,y=J))+scale_x_discrete(limits=c("Skin","Feces","Colon","Intestine","Gut","Stomach"),labels=c("Ski","Fec","Col","Int","Gut","Sto"))+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="")
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
plotJ
# Boxplots of the alpha diversity metrics
#pdf("plotHDJ_AdultTissue.pdf",width=6,height=10)
#library("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.2256845 0.1501432 0.1501393 0.1501404 0.1501436
# Checking that we have a correct sampling of the biodiversity in our samples
plot(specaccum(comm), xlab = "# of samples", ylab = "# of species",main="Species accumulation in adult tissues microbiomes " )
#pdf("plot_spec_accum_AdultTissue.pdf",width=7,height=5)
#print(plot(specaccum(comm), xlab = "# of samples", ylab = "# of species",main="Species accumulation in adult tissues microbiomes " ))
#dev.off()
# Testing significance 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$Tissue, data= comm)
kw_FaithPD<-kruskal.test(comm.pd$PD~metadata$Tissue, data= comm)
kw_Shannon<-kruskal.test(comm.H~metadata$Tissue, data= comm)
kw_Simpson<-kruskal.test(comm.D~metadata$Tissue, data= comm)
kw_Evenness<-kruskal.test(comm.J~metadata$Tissue, data= comm)
#Pairwise comparisons using Dunn-Bonferroni test
pp_spec<-kwAllPairsDunnTest(x=specnumber(comm), metadata$Tissue,p.adjust.method="bonferroni")
## Warning in kwAllPairsDunnTest.default(x = specnumber(comm), metadata$Tissue, :
## Ties are present. z-quantiles were corrected for ties.
pp_pd<-kwAllPairsDunnTest(x=comm.pd$PD, metadata$Tissue,p.adjust.method="bonferroni")
pp_shannon<-kwAllPairsDunnTest(x=comm.H, metadata$Tissue,p.adjust.method="bonferroni")
pp_Evenness<-kwAllPairsDunnTest(x=comm.J, metadata$Tissue,p.adjust.method="bonferroni")
#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")
#pdf("plot_bc_clust_AdultTissue.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")
#pdf("plot_so_clust_AdultTissue.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.1596667
## Run 1 stress 0.1777503
## Run 2 stress 0.1571168
## ... New best solution
## ... Procrustes: rmse 0.01637732 max resid 0.1048842
## Run 3 stress 0.2147483
## Run 4 stress 0.1882596
## Run 5 stress 0.1720643
## Run 6 stress 0.1782901
## Run 7 stress 0.1571169
## ... Procrustes: rmse 1.816624e-05 max resid 7.098485e-05
## ... Similar to previous best
## Run 8 stress 0.1816818
## Run 9 stress 0.1747871
## Run 10 stress 0.1697258
## Run 11 stress 0.1683307
## Run 12 stress 0.1641456
## Run 13 stress 0.1747239
## Run 14 stress 0.1953012
## Run 15 stress 0.2080515
## Run 16 stress 0.1571169
## ... Procrustes: rmse 1.935651e-05 max resid 7.546089e-05
## ... Similar to previous best
## Run 17 stress 0.157117
## ... Procrustes: rmse 4.528465e-05 max resid 0.0001190438
## ... Similar to previous best
## Run 18 stress 0.2123465
## Run 19 stress 0.1872079
## Run 20 stress 0.1695201
## *** Solution reached
#Checking the NMDS
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")
#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()
# We plot the NMDS ordination using Bray-Curtis
plot.new()
mds.fig2 <- ordiplot(comm.bc.mds, type = "none",xlim=c(-2,2),ylim=c(-2,2))
points(mds.fig2, "species", pch = 1,cex=0.5, col = "darkgrey")
ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Feces",col="#984ea3")
points(mds.fig2, "sites", pch = 19, col = "#984ea3", select = metadata$Tissue == "Feces")
ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Skin",col="black")
points(mds.fig2, "sites", pch = 19, col = "black", select = metadata$Tissue == "Skin")
ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Colon",col="#e41a1c")
points(mds.fig2, "sites", pch = 19, col = "#e41a1c", select = metadata$Tissue == "Colon")
ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Gut",col="#ff7f00")
points(mds.fig2, "sites", pch = 19, col = "#ff7f00", select = metadata$Tissue == "Gut")
ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Stomach",col="green")
points(mds.fig2, "sites", pch = 19, col = "green", select = metadata$Tissue == "Stomach")
ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Intestine",col="#377eb8")
points(mds.fig2, "sites", pch = 19, col = "#377eb8", select = metadata$Tissue == "Intestine")
ordiellipse(comm.bc.mds, metadata$Tissue, conf = 0.95, label = FALSE,lty = 'dotted')
#plot.new()
#par(new=TRUE)
#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 = "darkgrey"))
#print(ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Feces",col="#984ea3"))
#print(points(mds.fig2, "sites", pch = 19, col = "#984ea3", select = metadata$Tissue == "Feces"))
#print(ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Skin",col="black"))
#print(points(mds.fig2, "sites", pch = 19, col = "black", select = metadata$Tissue == "Skin"))
#print(ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Colon",col="#e41a1c"))
#print(points(mds.fig2, "sites", pch = 19, col = "#e41a1c", select = metadata$Tissue == "Colon"))
#print(ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Gut",col="#ff7f00"))
#print(points(mds.fig2, "sites", pch = 19, col = "#ff7f00", select = metadata$Tissue == "Gut"))
#print(ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Stomach",col="green"))
#print(points(mds.fig2, "sites", pch = 19, col = "green", select = metadata$Tissue == "Stomach"))
#print(ordispider(comm.bc.mds,groups=metadata$Tissue,show.groups="Intestine",col="#377eb8"))
#print(points(mds.fig2, "sites", pch = 19, col = "#377eb8", select = metadata$Tissue == "Intestine"))
#print(ordiellipse(comm.bc.mds, metadata$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$Tissue, type="median")
comm.betadisp
##
## Homogeneity of multivariate dispersions
##
## Call: betadisper(d = comm.bc.dist, group = metadata$Tissue, type =
## "median")
##
## No. of Positive Eigenvalues: 45
## No. of Negative Eigenvalues: 5
##
## Average distance to median:
## Colon Feces Gut Intestine Skin Stomach
## 0.1845 0.5245 0.3612 0.1792 0.5039 0.2860
##
## Eigenvalues for PCoA axes:
## (Showing 8 of 50 eigenvalues)
## PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
## 4.1560 1.9949 1.5860 1.1028 0.9317 0.8787 0.8269 0.7887
mds.fig3a<-plot(comm.betadisp)
#pdf("plot1_betadispa.pdf",width=8,height=6)
#print(mds.fig3a<-plot(comm.betadisp))
#dev.off()
# Significance testing using a permutest
comm.betadisp.perm=permutest(comm.betadisp, group= metadata$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 5 0.70807 0.14161 10.592 999 0.001 ***
## Residuals 45 0.60165 0.01337
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## Colon Feces Gut Intestine Skin Stomach
## Colon 2.0000e-03 1.0000e-03 9.7500e-01 1.0000e-03 0.441
## Feces 5.4853e-04 4.0000e-03 3.0000e-03 5.8300e-01 0.023
## Gut 3.0198e-04 3.0634e-03 3.4000e-02 3.0000e-03 0.306
## Intestine 9.6827e-01 1.4523e-03 3.0538e-02 1.0000e-03 0.559
## Skin 2.1829e-05 6.0663e-01 4.8687e-04 1.8056e-04 0.011
## Stomach 4.2684e-01 1.6190e-02 2.8527e-01 5.4793e-01 4.7323e-03
############################################
# 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 = 999)
comm.sesmntd <- ses.mntd(comm, phy.dist, null.model = "richness", abundance.weighted = FALSE, runs = 999)
# Looking at the values
#comm.sesmpd
#comm.sesmntd
# Significance
# Using Kruskal-Wallis test
kw_sesmpd<-kruskal.test(comm.sesmpd$mpd.obs.z ~ metadata$Tissue, data= comm)
kw_sesmntd<-kruskal.test(comm.sesmntd$mntd.obs.z ~ metadata$Tissue, data= comm)
# Pairwise comparisons using Dunn-Bonferroni test implemented in PMCMRplus
pp_sesmpd<-kwAllPairsDunnTest(x=comm.sesmpd$mpd.obs.z, metadata$Tissue,p.adjust.method="bonferroni")
pp_sesmntd<-kwAllPairsDunnTest(x=comm.sesmntd$mntd.obs.z, metadata$Tissue,p.adjust.method="bonferroni")
# Making letters corresponding to the significance level
mymat_sesmpd <-tri.to.squ(pp_sesmpd$p.value)
mymat_sesmntd <-tri.to.squ(pp_sesmntd$p.value)
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$Tissue)
df4<-data.frame(SES_MNTD = comm.sesmntd$mntd.obs.z, Sample = metadata$Tissue)
# Plotting the graphical comparisons of SES(MPD) and et SES(MNTD) values
plot3<-ggplot(df3,aes(x=Sample,y=SES_MPD))+geom_jitter(width=0.1,height=0.01,size=3,colour="darkgrey",alpha=0.9)+geom_boxplot(outlier.colour=NA,alpha=0.1)+scale_x_discrete(limits=c("Skin","Stomach","Intestine","Colon","Gut","Feces"),labels=c("Ski","Fec","Col","Int","Gut","Sto"))+geom_text(data = myletters_sesmpd_df, aes(label = letter, y = 5 ),fontface="bold",size=7)+theme_nb()
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
plot4<-ggplot(df4,aes(x=Sample,y=SES_MNTD))+geom_jitter(width=0.1,height=0.01,size=3,colour="darkgrey",alpha=0.9)+geom_boxplot(outlier.colour=NA,alpha=0.1)+scale_x_discrete(limits=c("Skin","Stomach","Intestine","Colon","Gut","Feces"),labels=c("Ski","Fec","Col","Int","Gut","Sto"))+geom_text(data = myletters_sesmntd_df, aes(label = letter, y = 1.3 ),fontface="bold",size=7)+theme_nb()
## Warning: `axis.ticks.margin` is deprecated. Please set `margin` property of
## `axis.text` instead
grid.arrange(plot3, plot4, ncol=2)
ggsave("plot_ses_adulttissues.pdf",width=6,height=4)
#####################################
# Phylogenetic Beta-diversity #
#####################################
# beta diversity: Taxonomic (Bray-Curtis) dissimilarity explained (Permanova)
adonis(comm.bc.dist ~ Tissue, data = metadata)
##
## Call:
## adonis(formula = comm.bc.dist ~ Tissue, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tissue 5 8.412 1.68240 6.9278 0.43495 0.001 ***
## Residuals 45 10.928 0.24285 0.56505
## Total 50 19.340 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 ~ Tissue, data = metadata)
##
## Call:
## adonis(formula = comm.mpd.dist ~ Tissue, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tissue 5 3.5845 0.71690 4.9223 0.35355 0.001 ***
## Residuals 45 6.5540 0.14565 0.64645
## Total 50 10.1386 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 ~ Tissue, data = metadata)
##
## Call:
## adonis(formula = comm.mntd.dist ~ Tissue, data = metadata)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Tissue 5 0.057063 0.0114126 24.1 0.7281 0.001 ***
## Residuals 45 0.021310 0.0004736 0.2719
## Total 50 0.078373 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.4656
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0481 0.0644 0.0808 0.1013
## 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.4881
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0477 0.0631 0.0758 0.0856
## Permutation: free
## Number of permutations: 999