Load packages

library("ampvis")

Load data

data(DNAext_1.0)

Species affected by eDNA removal

The DESeq2 package is used to test if any OTUs are in differential abundance after PMA treatment. There is no need to rarefy as DESeq2 handles different sample sizes nice. In addition, all low abundant OTUs are removed.

edna_deseq2 <- subset_samples(V13, Exp.PMA == "YES") %>%
  subset_taxa(Family != "f__Enterobacteriaceae") %>%
  filter_taxa(function(x) max(x) >= 10, TRUE)

To call significant we use an adjusted p-value of 0.001 and an abosolute cutoff of 0.5 log2fold change.

edna_test <- amp_test_species(data = edna_deseq2, 
                              group= "PMA", 
                              tax.add=c("Phylum","Genus"),
                              sig = 0.001, 
                              fold = 0.5, 
                              plot.type = "boxplot", 
                              plot.show = 10,
                              plot.theme = "clean",
                              plot.point.size = 1)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

Figure 2A: Effect of PMA treatment on all OTUs

Looking at the MA plot it can be seen that there are some OTUs that are in significant lower abundance after PMA treatment.

edna_test$plot_MA + 
  ylim(-4,4) +
  scale_x_log10(breaks = c(10, 100, 1000)) +
  scale_color_manual(labels = c("Not significant", "Significant"), values = c("black", "red"), name = "") +
  theme(legend.key.height = unit(3, "mm"),
        legend.position = c(0.3, 0.9),
        text = element_text(size = 10),
        axis.text = element_text(size = 8), 
        panel.grid.major = element_blank()) +
  annotate("text", x = 2500, y = 2.2, label = "+PMA", size = 3, color = "#00BFC4") +
  annotate("text", x = 2500, y = -2.2, label = "-PMA", size = 3, color = "#F8766D")

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

Figure 2B: Effect of PMA treatment on specific OTUs

197 out of 1643 OTUs had significantly changed abundance. Looking at the 10 most significantly changed species.

edna_test$plot_sig +
  scale_y_log10(breaks = c(0.01,0.1, 1)) +
  theme(legend.position = "none",
        text = element_text(size = 10),
        axis.text = element_text(size = 8)) 

ggsave("plots/Fig2B.eps", width = 100, height = 60, units = "mm")

Printing the 40 OTUs with most significant changed abundance.

grid.table(edna_test$clean_res[1:40,],
           show.rownames = F, 
           core.just = "right", 
           gpar.coretext = gpar(fontsize = 10), 
           gpar.coltext = gpar(fontsize = 12),
           padding.v = unit(2, "mm"),
           padding.h = unit(4, "mm"))