#load additional packages
library(BiocManager)
library(phyloseq)
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionImporting the data
# convert qiime artifact directly to phyloseq
library(phyloseq)
library(remotes)
library(qiime2R)
S4phy <- qza_to_phyloseq ("ch-filtered.qza",
"rooted-tree.qza",
"classification.qza",
"microbiome-meta.txt")
#Manual import
##Importing ASVs abundance file
ASVs <- read_qza ("ch-filtered.qza")
##Importing metadata
metadata <- read.table ("microbiome-meta.txt", , sep='\t', header = T, row.names= 1, comment= "")
metadata <- metadata[-1,] # remove the second line that specifies the data type
##Importing tree
tree <- read_qza("rooted-tree.qza")
##importing taxonomy
taxonomy <- read_qza("classification.qza")
tax_table <- do.call(rbind, strsplit(as.character(taxonomy$data$Taxon),";"))## Warning in (function (..., deparse.level = 1) : number of columns of result is
## not a multiple of vector length (arg 2)colnames(tax_table) <- c("Kingdom", "Phylum", "Class", "Order", "Family","Genus", "Species")
rownames(tax_table) <- taxonomy$data$Feature.ID
#creating phyloseq object
physeq <- phyloseq(otu_table(ASVs$data, taxa_are_rows = T),
phy_tree(tree$data),
tax_table(tax_table),
sample_data(metadata))
#merge the data
ps <- physeq
print(ps)## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3617 taxa and 69 samples ]
## sample_data() Sample Data: [ 69 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 3617 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3617 tips and 3614 internal nodes ]physeq.m = sample_data(physeq)Alpha Diversity
plot_richness(physeq, x="treatment", measures=c("Observed", "Shannon", "Simpson"), color="stage") + theme_bw()## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.
# First we will create a dataframe containing all diversity measures from above using the `estimate_richness()` function.
adiv <- estimate_richness(physeq, measures = c("Observed", "Shannon", "Simpson", "Chao1", "InvSimpson"))## Warning in estimate_richness(physeq, measures = c("Observed", "Shannon", : The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.# Create boxplots of the alpha diversity measures
simpsonbox <- ggplot(as.data.frame(physeq.m),
aes(x = treatment,
y = adiv$Simpson,
fill = stage)) +
labs(fill = "Stage", x = "Treatment", y = "Simpson's Diversity Index", title ="Effect of Treatment Type & Stage on Alpha Diversity") +
geom_boxplot() +
theme(title = element_text(size = 12)) +
theme_update(plot.title = element_text(hjust = 0.5))
# Creates Faith boxplot
shannonbox <- ggplot(as.data.frame(physeq.m), aes(x = treatment,
y = adiv$Observed,
fill = stage)) +
labs(fill = "Stage", x = "Treatment", y = "Shannon Diversity Index") +
geom_boxplot() +
theme(title = element_text(size = 12)) # makes titles smaller
# Puts them into same picture
gridExtra::grid.arrange(shannonbox, simpsonbox, nrow = 2)
# Then we can run an analysis of variance test using the `aov()` function.
summary(aov(adiv$Observed ~ physeq.m$treatment))## Df Sum Sq Mean Sq F value Pr(>F)
## physeq.m$treatment 3 1092 364 0.266 0.849
## Residuals 65 88797 1366Beta Diversity
Distance between samples
# Convert from raw to relative abundance
pn = transform_sample_counts(physeq, function(x) 100 * x/sum(x))
# Create Bray Curtis Distance Matrix
iDist <- distance(pn, method = "bray")
#Create NMDS ordination
pn.nmds = ordinate(pn,
method = "NMDS",
distance = iDist)## Run 0 stress 0.2780104
## Run 1 stress 0.2853044
## Run 2 stress 0.2815991
## Run 3 stress 0.2817538
## Run 4 stress 0.2873967
## Run 5 stress 0.2813362
## Run 6 stress 0.282923
## Run 7 stress 0.2978807
## Run 8 stress 0.2876882
## Run 9 stress 0.2862005
## Run 10 stress 0.2912546
## Run 11 stress 0.2832514
## Run 12 stress 0.2886704
## Run 13 stress 0.2813683
## Run 14 stress 0.2891912
## Run 15 stress 0.2962385
## Run 16 stress 0.2787542
## Run 17 stress 0.2815047
## Run 18 stress 0.2879669
## Run 19 stress 0.2815357
## Run 20 stress 0.2843201
## *** No convergence -- monoMDS stopping criteria:
## 20: stress ratio > sratmaxplot.pn.nmds = plot_ordination(pn, pn.nmds, justDF = TRUE)
#Display plot
print(ggplot(plot.pn.nmds, aes(x = NMDS1, y = NMDS2, color = outcome, shape = treatment)) +
geom_point( size = 4, alpha = 0.75)) 
#Display plot
print(ggplot(plot.pn.nmds, aes(x = NMDS1, y = NMDS2, color = treatment, shape = stage)) +
geom_point( size = 4, alpha = 0.75)) 
PERMANOVA on beta diversity scores
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-2p.df = as(sample_data(pn), "data.frame")
p.d = distance(pn, method = "bray")
p.adonis = adonis2(p.d ~ treatment + stage + outcome + treatment*stage + origin * treatment*outcome, p.df)
p.adonis## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = p.d ~ treatment + stage + outcome + treatment * stage + origin * treatment * outcome, data = p.df)
## Df SumOfSqs R2 F Pr(>F)
## treatment 3 1.5414 0.06451 2.3240 0.001 ***
## stage 1 0.3051 0.01277 1.3801 0.125
## outcome 1 1.4002 0.05860 6.3335 0.001 ***
## origin 13 8.0160 0.33549 2.7890 0.001 ***
## treatment:stage 3 0.5743 0.02404 0.8659 0.764
## treatment:origin 3 1.4544 0.06087 2.1928 0.001 ***
## outcome:origin 1 1.0950 0.04583 4.9528 0.001 ***
## Residual 43 9.5066 0.39788
## Total 68 23.8930 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(pairwiseAdonis) pairwise.adonis(p.d, factors = pn$treatment, sim.function = “vegdist”, sim.method = “bray”, p.adjust.m = “bonferroni”, reduce = NULL, perm = 999)
library(RVAideMemoire)## *** Package RVAideMemoire v 0.9-81-2 ***pairwise.perm.manova(p.d, p.df$treatment, nperm = 999, progress = TRUE, p.method = "fdr", F = FALSE, R2 = FALSE)## 'adonis' will be deprecated: use 'adonis2' instead## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead##
## Pairwise comparisons using permutation MANOVAs on a distance matrix
##
## data: p.d by p.df$treatment
## 999 permutations
##
## Chloramphenicol Combination Doxycycline
## Combination 0.087 - -
## Doxycycline 0.452 0.070 -
## None 0.313 0.057 0.057
##
## P value adjustment method: fdrlibrary(RVAideMemoire)
pairwise.perm.manova(p.d, p.df$outcome, nperm = 999, progress = TRUE, p.method = "fdr", F = FALSE, R2 = FALSE)## 'adonis' will be deprecated: use 'adonis2' instead##
## Pairwise comparisons using permutation MANOVAs on a distance matrix
##
## data: p.d by p.df$outcome
## 999 permutations
##
## EUTHANISED
## SURVIVED 0.001
##
## P value adjustment method: fdrhow do i do post-hoc tests on the permanova testing the interactions?
###post hoc ideas?? PAIRWISE.PERM.MANOVA BETADISPER
- SIMPER?
- ANOSIM -MANTEL
Differential Abundance Testing
#Differential Abundance Testing used to see which sequences differ between treatments and stages of treatment.