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