---
title: "Figure S10 Primer at Order"
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 order level abundances
```{r order}
pV13 <- amp_heatmap(data=V13n, tax.show = "all", output = "complete", scale.seq = 100, tax.empty="remove", tax.aggregate = "Order",tax.class = "p__Proteobacteria", tax.add = "Phylum")
pV4 <- amp_heatmap(data=V4n, tax.show = "all", output = "complete", scale.seq = 100, tax.empty="remove", tax.aggregate = "Order",tax.class = "p__Proteobacteria", tax.add = "Phylum")
pV34 <- amp_heatmap(data=V34n, tax.show = "all", output = "complete", scale.seq = 100, tax.empty="remove", tax.aggregate = "Order",tax.class = "p__Proteobacteria", tax.add = "Phylum")
pMT <- amp_heatmap(data=MTn, tax.show = "all", output = "complete", scale.seq = 100, tax.empty="remove", tax.aggregate = "Order",tax.class = "p__Proteobacteria", tax.add = "Phylum")
pMG <- amp_heatmap(data=MGn, tax.show = "all", output = "complete", scale.seq = 100, tax.empty="remove", tax.aggregate = "Order",tax.class = "p__Proteobacteria", tax.add = "Phylum")
```

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)))
```

# Order level comparison
## Aggregate the data
Only the 50 most abundant orders are shown.
```{r order_format}
c_order <- rbind.data.frame(pV13d, pV34d, pV4d, pMTd, pMGd) %>%
  group_by(Display, Experiment, Type) %>%
  summarise(Abundance = round(mean(Abundance),1))

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

Total <- group_by(c_order, Display) %>%
  summarise(Mean = mean(Abundance)) %>%
  arrange(desc(Mean)) %>%
  filter(row_number() < 51)

order <- subset(c_order, Display %in% Total$Display)
```

## Figure S10: Primer coverage at order level
```{r FigS10, fig.align='center', fig.width=6, fig.height=7}
ggplot(order, 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", na.value = "#EF8A62") +
  facet_grid(~Type, scales = "free_x", space = "free_x") +  
  theme(legend.position = "none",
        axis.text.x = element_text(size =6, color = "black", hjust = 0.4, angle = 0),
        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")
        )
```

```{r save_S10, eval=FALSE}
ggsave("plots/S10_Fig.eps", width = 100, height = 230, units = "mm")
```