---
title: "Figure 6B and 6C PCR settings"
author: "Mads Albertsen"
date: "`r format(Sys.time(), '%d-%m-%Y')`"
output: html_document
---

## Load packages
```{r Load_packages, message=FALSE, warning=FALSE, results='hide'}
library("ampvis")
```

## Load data
```{r load_data}
data(DNAext_1.0)
```

## Effect of annealing temperature
Subset to the relevant dataset.
```{r subset_temperature, message=FALSE}
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.
```{r PCA_temperature}
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. 
```{r pca_temperature_plot, fig.align='center', fig.height=5, fig.width=6}
pca_temperature$plot
```

The model reports a p-value of `r pca_temperature$eff.model$factors$pvals`, hence there is a small effect of the amount of temperature.
```{r pca_temperature_model}
pca_temperature$eff.model
```

The Bray-Curtis dissimilarity index is used as an alternative method to test for significant groupings in the dataset. 
```{r beta_temperature}
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 `r beta_temperature$adonis$aov.tab$"Pr(>F)"[1]`.
```{r beta_temperature_adonis}
beta_temperature$adonis
```

Clustering the data also shows that there might be a small effect. At least the 52 C and 58 C samples are different.
```{r beta_temperature_cluster, fig.align='center', fig.height=3, fig.width=6}
beta_temperature$plot_cluster
```

## Species abundance 52 C vs. 58C
```{r deseq2_subset}
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.
```{r temperature_sig}
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")
```

Looking at the MA plot it can be seen that there are some OTUs that are in differential abundance.
```{r temperature_sig_MA, fig.align='center', fig.height=5, fig.width=6}
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!
```{r temperature_sig_box, fig.align='center', fig.height=5, fig.width=12}
temperature_species_test$plot_sig
```

## Figure 6B: MA plot of differential abundant species as function of annealing temperature
```{r Fig6B, fig.align='center', fig.width=2, fig.height=2, warning=FALSE}
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"))
```

## Save the plot
```{r save1, eval=FALSE}
ggsave("plots/Fig6B.eps", width = 50, height = 50, units = "mm")
```

## Figure 6C: point plot of differential abundant species as function of annealing temperature
```{r Fig6C, fig.align='center', fig.width=3.5, fig.height=2}
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
```{r save2, eval=FALSE}
ggsave("plots/Fig6C.eps", width = 75, height = 50, units = "mm")
```