---
title: "Figure 5A Primer Comparison"
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)
```

## Subset and normalise the data
```{r subset}
V13n <- subset_samples(V13,Exp.beadbeating == "YES" & Beadbeating == "160s6ms") %>%
  transform_sample_counts(function(x) x / sum(x) * 100)

V34n <- transform_sample_counts(V34, function(x) x / sum(x) * 100)

V4n <- subset_samples(V4,Exp.beadbeating == "YES" & Beadbeating == "160s6ms") %>%
  transform_sample_counts(function(x) x / sum(x) * 100)

MTn <- subset_taxa(MT, Kingdom == "k__Bacteria") %>%
  transform_sample_counts(function(x) x / sum(x) * 100)

MGn <- subset_taxa(MG, Kingdom == "k__Bacteria") %>%
  transform_sample_counts(function(x) x / sum(x) * 100)
```

## Extract phylum level abundances
```{r phylum}
pV13 <- amp_heatmap(data=V13n, tax.show = "all", output = "complete", scale.seq = 100, tax.empty="remove", tax.aggregate = "Phylum",tax.class = "p__Proteobacteria")
pV4 <- amp_heatmap(data=V4n, tax.show = "all", output = "complete", scale.seq = 100, tax.empty="remove", tax.aggregate = "Phylum",tax.class = "p__Proteobacteria")
pV34 <- amp_heatmap(data=V34n, tax.show = "all", output = "complete", scale.seq = 100, tax.empty="remove", tax.aggregate = "Phylum",tax.class = "p__Proteobacteria")
pMT <- amp_heatmap(data=MTn, tax.show = "all", output = "complete", scale.seq = 100, tax.empty="remove", tax.aggregate = "Phylum",tax.class = "p__Proteobacteria")
pMG <- amp_heatmap(data=MGn, tax.show = "all", output = "complete", scale.seq = 100, tax.empty="remove", tax.aggregate = "Phylum",tax.class = "p__Proteobacteria")
```

Convert to a data frames for later aggregation.
```{r to_dataframe}
pV13d <- cbind.data.frame(pV13$data[,c(1,3,5)], 
                          Experiment = rep("V1-3", nrow(pV13$data)), 
                          Type = rep("Amplicon", nrow(pV13$data)))
pV34d <- cbind.data.frame(pV34$data[,c(1,3,5)], 
                          Experiment = rep("V3-4", nrow(pV34$data)), 
                          Type = rep("Amplicon", nrow(pV34$data)))
pV4d <- cbind.data.frame(pV4$data[,c(1,3,5)], 
                          Experiment = rep("V4", nrow(pV4$data)), 
                          Type = rep("Amplicon", nrow(pV4$data)))
pMTd <- cbind.data.frame(pMT$data[,c(1,3,5)], 
                          Experiment = rep("MT", nrow(pMT$data)), 
                          Type = rep("Meta-", nrow(pMT$data)))
pMGd <- cbind.data.frame(pMG$data[,c(1,3,5)], 
                          Experiment = rep("MG", nrow(pMG$data)), 
                          Type = rep("Meta-", nrow(pMG$data)))
```

## Theoretical coverage using Silva
### Load data from the three different sets of primers
The data was generated using Silva's `test prime` function.

### Extract the relevant phyla
```{r mm_extract}
pshow <- c("Bacteria;Proteobacteria;Betaproteobacteria;",
  "Bacteria;Actinobacteria;",
  "Bacteria;Bacteroidetes;",
  "Bacteria;Proteobacteria;Alphaproteobacteria;",
  "Bacteria;Firmicutes;",
  "Bacteria;Chloroflexi;",
  "Bacteria;Proteobacteria;Gammaproteobacteria;",
  "Bacteria;Proteobacteria;Deltaproteobacteria;",
  "Bacteria;Nitrospirae;",
  "Bacteria;Chlorobi;",
  "Bacteria;Acidobacteria;")

pTH <- filter(TH, Mismatch == 1) %>% 
  filter(taxonomy %in% pshow)

pTH$taxonomy <- gsub("Proteobacteria;", "", pTH$taxonomy)
pTH$taxonomy <- gsub("Bacteria;", "", pTH$taxonomy)
pTH$taxonomy <- gsub(";", "", pTH$taxonomy)
```

Finally, the extracted data is also converted to a dataframe.
```{r to_dataframe2}
pTHd <- data.frame(Display = pTH$taxonomy, 
                   Group = rep("TH", nrow(pTH)), 
                   Abundance = pTH$coverage, 
                   Experiment = pTH$Primer, 
                   Type = "Theoretical")
```

## Aggregate the data
Only the 11 most abundant phyla are shown.
```{r Phylum_format}
c_phylum <- rbind.data.frame(pV13d, pV34d, pV4d, pMTd, pMGd, pTHd) %>%
  group_by(Display, Experiment, Type) %>%
  summarise(Abundance = round(mean(Abundance),1))

c_phylum$Experiment <- factor(c_phylum$Experiment, levels = c("V1-3","V3-4","V4", "MG", "MT", "TH"))

Total <- filter(c_phylum, Type != "Theoretical") %>%
  group_by(Display) %>%
  summarise(Mean = mean(Abundance)) %>%
  arrange(desc(Mean)) %>%
  filter(row_number() < 12)

phylum <- subset(c_phylum, Display %in% Total$Display)
```

## Figure 5A_1: Primer coverage at Phylum level
The figure is build in 2 seperate layers to facilitate different scales on the heatmap.
```{r Fig5A_1, fig.align='center', fig.width=3.5, fig.height=2.8}
ggplot(phylum, aes(x = Experiment, y = Display, label = Abundance)) + 
  geom_tile(aes(fill = Abundance), colour = "white", size = 0.5) + 
  labs(x = "", y = "", fill = "Abundance") + 
  geom_text(colour = "black", size = 2) +
  scale_fill_gradientn(colours = brewer.pal(3, "RdBu"), trans = "log10", limits = c(1,30)) +
  facet_grid(~Type, scales = "free_x", space = "free_x") +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 8, color = "black", vjust = 0.4, hjust = 1, angle = 90),
        axis.text.y = element_text(size = 8, color = "black"),
        axis.title = element_blank(),
        axis.ticks.length = unit(1, "mm"),
        strip.text = element_text(size = 8, color = "black"),
        strip.background = element_blank(),
        plot.margin = unit(c(0,0,0,0), "mm")
        )
```

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

## Figure 5A_2: Primer coverage at Phylum level
```{r Fig5A_2, fig.align='center', fig.width=3.5, fig.height=2.8}
ggplot(phylum, aes(x = Experiment, y = Display, label = Abundance)) + 
  geom_tile(aes(fill = Abundance), colour = "white", size = 0.5) + 
  labs(x = "", y = "", fill = "Abundance") + 
  geom_text(colour = "black", size = 2) +
  scale_fill_gradientn(colours = brewer.pal(3, "RdBu"), trans = "log10", limits = c(80,100)) +
  facet_grid(~Type, scales = "free_x", space = "free_x") +  
  theme(legend.position = "none",
        axis.text.x = element_text(size = 8, color = "black", vjust = 0.4, hjust = 1, angle = 90),
        axis.text.y = element_text(size = 8, color = "black"),
        axis.title = element_blank(),
        axis.ticks.length = unit(1, "mm"),
        strip.text = element_text(size = 8, color = "black"),
        strip.background = element_blank(),
        plot.margin = unit(c(0,0,0,0), "mm")
        )
```

## Save the plot
```{r save2, eval=FALSE}
ggsave("plots/Fig5A_2.eps", width = 90, height = 70, units = "mm")
```
