Load packages

library("ampvis")

Load data

data(DNAext_1.0)

Effect of annealing temperature

Subset to the relevant dataset.

pcr_temperature <- subset_samples(V13, Exp.PCR.temp == "YES") %>%
  rarefy_even_depth(sample.size = 25000, rngseed = 712) %>%
  filter_taxa(function(x) max(x) >= 10, TRUE)

Calculate PCA.

pca_temperature <- amp_ordinate(data = pcr_temperature, 
             plot.color = "PCR.temp", 
             plot.point.size = 5,
             plot.group = "chull",
             plot.group.label = "PCR.temp",
             envfit.factor = "PCR.temp",
             output = "complete"
             )

Plot the PCA. It looks like there might be some significant grouping.

pca_temperature$plot

The model reports a p-value of 0.006, hence there is a small effect of the amount of temperature.

pca_temperature$eff.model
## 
## ***FACTORS:
## 
## Centroids:
##                  PC1     PC2
## PCR.temp52 C  3.1952  0.7337
## PCR.temp56 C -0.7427 -2.5089
## PCR.temp58 C -2.4524  1.7752
## 
## Goodness of fit:
##              r2 Pr(>r)   
## PCR.temp 0.7933  0.006 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

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

beta_temperature <- amp_test_cluster(data = pcr_temperature, group = "PCR.temp", method = "bray", plot.color = "PCR.temp", plot.label = "PCR.temp")

Using adonis we also a small significant effect of temperature as the p-value is 0.003.

beta_temperature$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)   
## PCR.temp   2  0.016720 0.0083598  1.8118 0.37653  0.003 **
## Residuals  6  0.027685 0.0046142         0.62347          
## Total      8  0.044405                   1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Clustering the data also shows that there might be a small effect. At least the 52 C and 58 C samples are different.

beta_temperature$plot_cluster

Species abundance 52 C vs. 58C

pcr_temperature_test <- subset_samples(V13, Exp.PCR.temp == "YES" & PCR.temp != "56 C") %>%
  filter_taxa(function(x) max(x) >= 10, TRUE)

The DESeq2 package is used to test if any OTUs are in differential abundance.

temperature_species_test <- amp_test_species(data = pcr_temperature_test, 
                                             group = "PCR.temp", 
                                             tax.add = c("Phylum","Family"), 
                                             sig = 0.001, 
                                             fold = 0.5, 
                                             plot.type = "point", 
                                             plot.show = 10, 
                                             tax.class = "p__Proteobacteria", 
                                             plot.point.size = 1, 
                                             plot.theme = "clean")
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

Looking at the MA plot it can be seen that there are some OTUs that are in differential abundance.

temperature_species_test$plot_MA + ylim(-5,5)

Looking at the 10 most significantly changed species. Most of these are undetected at the high annealing temperature!

temperature_species_test$plot_sig

Figure 6B: MA plot of differential abundant species as function of annealing temperature

temperature_species_test$plot_MA +
  xlab("Base mean (Counts)") +
  ylim(-5,5) +
  scale_color_manual(labels = c("Not significant", "Significant"), values=c("black", "red")) + 
  annotate("text", x = 5000, y = -3, label = "52*' '*degree*C", size = 2, size = 3, color = "#00B6EB", parse = T) +
  annotate("text", x = 5000, y = 3, label = "58*' '*degree*C", size = 2, size = 3, color = "#A58AFF", parse = T) +
  theme(legend.position = c(0.3,0.9),
        legend.title = element_blank(),
        text = element_text(size = 8, color = "black"),
        axis.text = element_text(size = 6, color = "black"),
        legend.key.height = unit(3, "mm"),
        legend.key.width = unit(2, "mm"))
## Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Save the plot

ggsave("plots/Fig6B.eps", width = 50, height = 50, units = "mm")

Figure 6C: point plot of differential abundant species as function of annealing temperature

id <- levels(temperature_species_test$plot_sig$data$Tax)
id2 <- data.frame(do.call('rbind', str_split(id,';',3))) 
id3 <- paste(id2[,1], id2[,2], sep = ";")
species <- gsub(" ", "", id2[,3])

pcr_n <- subset_samples(V13, Exp.PCR.temp == "YES") %>%
  transform_sample_counts(function(x) x / sum(x) * 100) %>%
  prune_taxa(taxa = species)

out <- amp_rabund(pcr_n, group = "PCR.temp", plot.type = "point", tax.aggregate = "OTU", scale.seq = 100, output = "complete")

out$data$Display <- factor(out$data$Display, levels = species)

ggplot(out$data, aes(x = Abundance, y = Display, color = Group)) +
  geom_point(size = 1) +
  scale_x_continuous(breaks = c(0,0.1,0.2,0.3)) +
  scale_y_discrete(labels = id3) +
  scale_color_manual(values=c("#00B6EB","#FB61D7","#A58AFF")) +
  ylab("") +
  xlab("Read abundance (%)") +
  theme(legend.position = "none",
        axis.ticks.length = unit(1, "mm"),
        axis.ticks = element_line(color = "black"),
        text = element_text(size = 8, color = "black"),
        axis.text = element_text(size = 6, color = "black"),
        plot.margin = unit(c(0,0,0,0), "mm"),
        legend.key.width = unit(3, "mm"),
        legend.key.height = unit(3, "mm"),
        legend.key = element_blank(),
        legend.title = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(color = "black")
        )

Save the plot

ggsave("plots/Fig6C.eps", width = 75, height = 50, units = "mm")