# Set seed for reproducibility
set.seed(1234)
n_species <- 30
# Create groups
control <- matrix(
rpois(n_species * 20, lambda = runif(n_species, 10, 50)),
nrow = 20, ncol = n_species
)
case <- control
up_species <- sample(1:n_species, size = n_species/2)
case[, up_species] <- case[, up_species] + rpois(length(up_species) * 20, lambda = 30)
case[, -up_species] <- pmax(0, case[, -up_species] - rpois(length(-up_species) * 20, lambda = 10))
# Combine groups
otu_table <- rbind(control, case)
rownames(otu_table) <- paste0("Sample", 1:40)
colnames(otu_table) <- paste0("Species", 1:n_species)
# Metadata
metadata <- data.frame(
SampleID = rownames(otu_table),
Treatment = rep(c("control", "case"), each = 20)
)
## Warning: package 'vegan' was built under R version 4.5.2
## Warning: package 'permute' was built under R version 4.5.2
p <- apply(otu_table, 1, function(x){
x=x[x!=0]
x/sum(x)
})
shannon <- sapply(p, function(x){
(-sum(x*log(x)))
})
diversity_tbl <- cbind(metadata, shannon)
plot<- ggplot(diversity_tbl, aes(x=Treatment, y=shannon))+
geom_boxplot(outlier.shape = 2)+
theme(axis.text.x = element_text(angle = 45, hjust = 1, face = c("italic")))+
labs(y= "Relative abundance (%)")+
geom_jitter(aes(color = shannon))+
theme(legend.position = "none")
plot
compare <- list(c("case", "control"))
plot + stat_compare_means(method = "t.test", comparisons = compare)
ordination <- metaMDS(otu_table, distance = "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1492806
## Run 1 stress 0.1546364
## Run 2 stress 0.1543374
## Run 3 stress 0.1543374
## Run 4 stress 0.1492804
## ... New best solution
## ... Procrustes: rmse 8.715526e-05 max resid 0.0004763319
## ... Similar to previous best
## Run 5 stress 0.1500269
## Run 6 stress 0.1490274
## ... New best solution
## ... Procrustes: rmse 0.01415732 max resid 0.08202732
## Run 7 stress 0.1500269
## Run 8 stress 0.19911
## Run 9 stress 0.1543374
## Run 10 stress 0.1492802
## ... Procrustes: rmse 0.01453476 max resid 0.08434747
## Run 11 stress 0.1490274
## ... Procrustes: rmse 2.661767e-05 max resid 6.914301e-05
## ... Similar to previous best
## Run 12 stress 0.1492802
## ... Procrustes: rmse 0.01455681 max resid 0.08447234
## Run 13 stress 0.1490274
## ... New best solution
## ... Procrustes: rmse 7.906731e-06 max resid 4.03314e-05
## ... Similar to previous best
## Run 14 stress 0.1492802
## ... Procrustes: rmse 0.01443135 max resid 0.08375218
## Run 15 stress 0.1492803
## ... Procrustes: rmse 0.01465414 max resid 0.08501035
## Run 16 stress 0.1490274
## ... Procrustes: rmse 2.105977e-05 max resid 0.000107871
## ... Similar to previous best
## Run 17 stress 0.1490274
## ... Procrustes: rmse 1.144579e-05 max resid 5.73852e-05
## ... Similar to previous best
## Run 18 stress 0.1500269
## Run 19 stress 0.1500269
## Run 20 stress 0.1500269
## *** Best solution repeated 3 times
site <- as.data.frame(scores(ordination, display = "sites"))
site <- cbind(site, Treatment = metadata$Treatment)
theme_set(theme_bw())
nmds <- ggplot(site, aes(x=NMDS1, y=NMDS2)) +
geom_point(aes(NMDS1, NMDS2, colour = factor(site$Treatment),shape = factor(site$Treatment)), size = 6) +
labs(colour = "Treatment", shape = "Treatment")
nmds
## Warning: Use of `site$Treatment` is discouraged.
## ℹ Use `Treatment` instead.
## Use of `site$Treatment` is discouraged.
## ℹ Use `Treatment` instead.
dist <- vegdist(otu_table, method = "bray")
permanova <- adonis2(dist ~ Treatment, data = site, permutations = 999)
permanova
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist ~ Treatment, data = site, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.7537 0.39923 25.252 0.001 ***
## Residual 38 1.1342 0.60077
## Total 39 1.8879 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1