# Xenopus metatetard XPDEV
# Analysis using phyloseq
library(phyloseq)
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(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(phytools)
## Loading required package: maps
## 
## Attaching package: 'phytools'
## The following object is masked from 'package:vegan':
## 
##     scores
## Folders, Themes, Colors
source("prelude.R")

# Setting up directories
data_dir <-paste0(data_dir_path,"xpdev")
output_dir <- paste0(output_dir_path,"Figure2_ordination")


##################################
#   Ordination on rarefied data  #
##################################
# Defining the directory containing the data to import
setwd(data_dir)
# Working on abundance data from 100 rarefactions
comm<-read.table("xpdev_fichier_allcom.txt",header=TRUE,row.names=1)
# Getting rid of OTUs with zero abundance
comm<-comm %>% select_if((function(col) is.numeric(col) && sum(col) > 0))
comm<-t(comm)
#Importing metadata
metadata<-read.csv("nmt3_2_metadata_xp_dev.csv",sep=";",header=TRUE,row.names=1)
taxmat<-read.csv("nmt3_Gal18_abondance_tax.csv",sep=";",header=TRUE,row.names=1)
taxmat<-as.matrix(taxmat)
OTU = otu_table(comm, taxa_are_rows = TRUE)
TAX = tax_table(taxmat)
sample_data = sample_data(metadata)
#Importing phylogeny
phy<-read.tree("nmt3_sequences_mpr.tree")
#Tree midpoint rooting
phy<-midpoint.root(phy)
physeq.ini = merge_phyloseq(OTU,sample_data,TAX,phy)

# Testing the significance of unifrac distances across Life Stages

# beta diversity: compute unweighted unifrac distance
uudist<-distance(physeq.ini,method="uunifrac")
adonis(uudist~LifeStage,data=metadata)
## 
## Call:
## adonis(formula = uudist ~ LifeStage, data = metadata) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## LifeStage  4    2.4826 0.62065  9.5308 0.57655  0.001 ***
## Residuals 28    1.8234 0.06512         0.42345           
## Total     32    4.3059                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# beta diversity: compute weighted unifrac distance
uwdist<-distance(physeq.ini,method="wunifrac")
adonis(uwdist~LifeStage,data=metadata)
## 
## Call:
## adonis(formula = uwdist ~ LifeStage, data = metadata) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## LifeStage  4    2.2626 0.56566  11.822 0.62811  0.001 ***
## Residuals 28    1.3397 0.04785         0.37189           
## Total     32    3.6023                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LifeStage significantly impacts the communities both in structure and in composition.

#############################
#   Ordination on raw data  #
#############################
# Defining the directory containing the data to import
setwd(data_dir)
# Working on abundance data from 100 rarefactions
comm<-read.table("nmt3_abondance_tpn_notax_XP_Dev.asc",header=TRUE,row.names=1)
# Getting rid of OTUs with zero abundance
comm<-comm %>% select_if((function(col) is.numeric(col) && sum(col) > 0))
comm<-t(comm)
#Importing metadata
metadata<-read.csv("nmt3_2_metadata_xp_dev.csv",sep=";",header=TRUE,row.names=1)
taxmat<-read.csv("nmt3_Gal18_abondance_tax.csv",sep=";",header=TRUE,row.names=1)
taxmat<-as.matrix(taxmat)
OTU = otu_table(comm, taxa_are_rows = TRUE)
TAX = tax_table(taxmat)
sample_data = sample_data(metadata)
#Importing phylogeny
phy<-read.tree("nmt3_sequences_mpr.tree")
#Tree midpoint rooting
phy<-midpoint.root(phy)
physeq.ini.raw = merge_phyloseq(OTU,sample_data,TAX,phy)

# Testing the significance of unifrac distances across Life Stages

# beta diversity: compute unweighted unifrac distance
uudist<-distance(physeq.ini.raw,method="uunifrac")
adonis(uudist~LifeStage,data=metadata)
## 
## Call:
## adonis(formula = uudist ~ LifeStage, data = metadata) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## LifeStage  4    2.4826 0.62065  9.5308 0.57655  0.001 ***
## Residuals 28    1.8234 0.06512         0.42345           
## Total     32    4.3059                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# beta diversity: compute weighted unifrac distance
uwdist<-distance(physeq.ini.raw,method="wunifrac")
adonis(uwdist~LifeStage,data=metadata)
## 
## Call:
## adonis(formula = uwdist ~ LifeStage, data = metadata) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## LifeStage  4    2.2627 0.56567  11.823 0.62812  0.001 ***
## Residuals 28    1.3396 0.04784         0.37188           
## Total     32    3.6023                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1