Xenopus metatetard XPDEV
# Analysis using phyloseq
library(phyloseq)
library(ggplot2)
library(dplyr)
library(phytools)
library(philr)
library(picante)
## Folders, Themes, colors
source("prelude.R")
# Setting up directories
data_dir <-paste0(data_dir_path,"xptransmission")
output_dir <- paste0(output_dir_path,"Figure_late")
# Importing data
# Defining the directory containing the data to import
setwd(data_dir)
# Importing abundance data
comm<-read.table("nmt3_abondance_tpn_notax_xp_transa.csv",header=TRUE,row.names=1,sep=";")
# 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_metadata_XP_transa.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")
#Rooting the tree using midpoint rooting
phy<-midpoint.root(phy)
physeq.ini = merge_phyloseq(OTU,sample_data,TAX,phy)
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
# Filtering OTUs
MTDEV <- filter_taxa(physeq.ini, function(x) sum(x > 3) > (0.05*length(x)), TRUE)
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
# We use pseudocounts of 1 to compute logs
MTDEV <- transform_sample_counts(MTDEV, function(x) x+1)
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
#Check that phylogeny is rooted
is.rooted(phy_tree(MTDEV))
## [1] TRUE
is.binary.tree(phy_tree(MTDEV))
## [1] TRUE
# Naming tree nodes
phy_tree(MTDEV) <- makeNodeLabel(phy_tree(MTDEV), method="number", prefix='n')
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
# Reformatting OTU table for philr
otu.table <- t(otu_table(MTDEV))
tree <- phy_tree(MTDEV)
metadata <- sample_data(MTDEV)
tax <- tax_table(MTDEV)
# Using philR metric
gp.philr <- philr(otu.table, tree, part.weights='enorm.x.gm.counts',ilr.weights='blw.sqrt')
## Building Sequential Binary Partition from Tree...
## Building Contrast Matrix...
## Transforming the Data...
## Calculating ILR Weights...
#Setting the working directory for output files
setwd(output_dir)
#PCoA Ordination
gp.dist <- dist(gp.philr, method="euclidean")
gp.pcoa <- ordinate(MTDEV, 'PCoA', distance=gp.dist)
ploto<-plot_ordination(MTDEV, gp.pcoa, shape='XP_Tissue') +geom_point(size=3,alpha=0.5,fill="grey")+scale_shape_manual(values=c(15,0,16,1,6))+labs(shape="Sample")
ploto<-ploto+ggtitle("Community composition")+theme_npgray()+ stat_ellipse(type="norm",linetype=2)+theme(legend.position="top")
ploto
ggsave("PCOA_transmission_philr.pdf",width=8,height=6)
#perMANOVA (adonis) analysis
adonis(gp.dist ~ XP_Tissue, as(sample_data(MTDEV), "data.frame"))
##
## Call:
## adonis(formula = gp.dist ~ XP_Tissue, data = as(sample_data(MTDEV), "data.frame"))
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## XP_Tissue 2 5993.0 2996.52 23.145 0.57653 0.001 ***
## Residuals 34 4401.9 129.47 0.42347
## Total 36 10395.0 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1