Load packages

library("ampvis")

Load data

data(DNAext_1.0)

Check that we were able to remove spike in DNA from e.coli

ecoli <- subset_samples(V13, Exp.PMA == "YES") %>%
  rarefy_even_depth(sample.size = 25000, rngseed = 712) %>%
  subset_taxa(Family == "f__Enterobacteriaceae")
## `set.seed(712)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(712); .Random.seed` for the full vector
## ...
## 7943OTUs were removed because they are no longer 
## present in any sample after random subsampling
## 
## ...

We did sequence some of the spiked in DNA.

amp_heatmap(ecoli, group = c("PMA","E.coli"), tax.aggregate = "Family", scale.seq = 25000) +
  scale_x_discrete(labels = c("0%", "20%","40%","50%","0%\nPMA","20%\nPMA","40%\nPMA","50%\nPMA")) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 0, hjust = 0.3))

Overall differences between samples

Remove spiked in e.coli

Remove e.coli and rarefy to even depth in order to use ordination afterwards. In addition, low abundant OTUs are removed.

edna_clean_rare <- subset_samples(V13, Exp.PMA == "YES") %>%
  subset_taxa(Family != "f__Enterobacteriaceae") %>%
  rarefy_even_depth(sample.size=25000, rngseed=712) %>%
  filter_taxa(function(x) max(x) >= 10, TRUE)
## `set.seed(712)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(712); .Random.seed` for the full vector
## ...
## 7958OTUs were removed because they are no longer 
## present in any sample after random subsampling
## 
## ...

Figure S3C: PCA colored by PMA treatment

pca <- amp_ordinate(data = edna_clean_rare, 
             plot.color = "PMA",             
             plot.point.size = 3,
             envfit.factor = "PMA",
             envfit.show = F,
             output = "complete",
             plot.theme = "clean"
             )

It looks like there is significant groupings.

pca$plot +
    theme(legend.key.height = unit(3, "mm"))

ggsave("plots/S4C.eps", width = 80, height = 55, units = "mm")

The model reports a p-value of 0.001, hence there is a significant effect of PMA treatment.

pca$eff.model
## 
## ***FACTORS:
## 
## Centroids:
##            PC1     PC2
## PMANo   1.9826 -0.4351
## PMAYes -1.9826  0.4351
## 
## Goodness of fit:
##         r2 Pr(>r)    
## PMA 0.4693  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

Figure S4D: Cluster analysis of beta diversity using Bray-Curtis

The Bray-Curtis dissimilarity index is used as an alternative method to test for significant groupings in the dataset.

beta <- amp_test_cluster(data = edna_clean_rare, 
                         group = "PMA", 
                         method = "bray", 
                         plot.color = "PMA", 
                         plot.label = "PMA",
                         plot.theme = "clean")

Using adonis we also see a significant effect of eDNA removal as the p-value is 0.001.

beta$adonis
## 
## Call:
## adonis(formula = test_formula, data = sample) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)    
## PMA        1  0.082235 0.082235  11.429 0.34189  0.001 ***
## Residuals 22  0.158293 0.007195         0.65811           
## Total     23  0.240528                  1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Clustering also shows that there is a significant effect.

beta$plot_cluster +
  theme(legend.position = "none")

ggsave("plots/S3D.eps", width = 60, height = 55, units = "mm")

Figure S3E: Variance compared to time-series data

edna_time <- subset_samples(V13, (Exp.PMA == "YES" & E.coli == "0%") | Exp.time == "YES") %>%
  rarefy_even_depth(sample.size=25000, rngseed=712) %>%
  filter_taxa(function(x) max(x) > 10, TRUE)
## `set.seed(712)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(712); .Random.seed` for the full vector
## ...
## 7775OTUs were removed because they are no longer 
## present in any sample after random subsampling
## 
## ...

There seem to a quite large effect also compared to the varience over time.

amp_ordinate(data = edna_time, 
             plot.color = "Date",             
             plot.point.size = 3,
             plot.theme = "clean",
             plot.shape = "PMA"
             ) +
  scale_color_discrete(name = "Sampling date") +
  theme(legend.key.height = unit(3, "mm"))

ggsave("plots/S3E.eps", width = 90, height = 55, units = "mm")

Figure S3F: Using clustering to estimate classification resolution

beta_time <- amp_test_cluster(data = edna_time, 
                              group = "Date", 
                              method = "bray", 
                              plot.color = "Date", 
                              plot.label = "PMA",
                              plot.theme = "clean")

The PMA treated samples are still most closely related to the samples from the same timepoint. However, they are clearly disinct and the difference seem to be much larger than the weekly variance.

beta_time$plot_cluster +
  theme(legend.position = "none")

ggsave("plots/S3F.eps", width = 60, height = 55, units = "mm")

Figure S3B: Impact of eDNA removal on Alpha-diversity

Subset to the relevant dataset and remove e.coli spike in.

edna_alpha <- subset_samples(V13, Exp.PMA == "YES") %>%
  subset_taxa(Family != "f__Enterobacteriaceae") %>%
  rarefy_even_depth(sample.size = 25000, rngseed = 712)
## `set.seed(712)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(712); .Random.seed` for the full vector
## ...
## 7958OTUs were removed because they are no longer 
## present in any sample after random subsampling
## 
## ...

Calculate alpha diverisity.

alpha <- cbind.data.frame(estimate_richness(edna_alpha), sample_data(edna_alpha))

Make a simple t-test.

t.test(alpha$Observed~alpha$PMA)
## 
##  Welch Two Sample t-test
## 
## data:  alpha$Observed by alpha$PMA
## t = 5.0512, df = 18.926, p-value = 7.174e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   77.87494 188.12506
## sample estimates:
##  mean in group No mean in group Yes 
##          2235.667          2102.667

Plot the data.

ggplot(data = alpha, aes(y = Observed, x = PMA)) +
  geom_boxplot() + 
  geom_jitter(color = "darkred", size = 1) +
  ylab("Observed number of OTUs") +
  xlab("PMA treatment") +
  theme(legend.position = "none",
      text = element_text(size = 8, color = "black"),
      axis.text = element_text(size = 8, color = "black"),
      axis.text.y = element_text(hjust = 1, size = 6),
      plot.margin = unit(c(0,0,0,0), "mm"),
      axis.line = element_line(color = "black"),
      panel.grid = element_blank(),
      axis.ticks = element_line(color = "black"),
      panel.background = element_blank())

ggsave("plots/S3B.eps", width = 60, height = 55, units = "mm")