Abundance_plots_metabarcoding.Rmd

---
title: "Abundance_metabarcoding_090519"
author: "Georgina Brennan" 
date: "5/10/2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
Load packages
```{r, results='hide',message=FALSE}
library(ggplot2)
library(reshape2)
library(tidyr)
#library(tidyverse)
library(RColorBrewer)
#library(cowplot)
#library(grid)
library(plyr)


```
load tree height metabarcoding count data and the mapping file or metadata

```{r}

OTUdata<-read.csv("bact_alldata_taxatable_wTax.csv", row.names = 1, check.names =FALSE)
str(OTUdata)

meta<-read.table("bact_alldata_mapfile.txt")
names(meta)
```
Remove taxonomy before transposing the data
```{r}
OTUdata$taxonomy<-NULL

```
Transpose the data to have sample names on rows

```{r}
OTUt<-as.data.frame(t(OTUdata))

```
Apply proportion normalisation

```{r}
OTUdata<-OTUt/rowSums(OTUt)
OTUdata<-OTUdata[,order(colSums(OTUdata),decreasing=TRUE)]
```
#Extract list of top N Taxa

```{r}
N<-20
taxa_list<-colnames(OTUdata)[1:N]
N<-length(taxa_list)
```

Generate a new table with everything in the top N list

```{r}

abundOTU<-data.frame(OTUdata[,colnames(OTUdata) %in% taxa_list])
names(abundOTU)

```
Decide how you want to filter your data for visualization 
Here I am looking at the most the most abundant OTUS - top 20 abundant OTUs across all sampling ID's

```{r}
abundOTU$Index<-rownames(abundOTU)
OTUt.prop <- merge(meta, abundOTU, by.y = 'Index', by.x = "SampleNo", all=FALSE)

# check for errors in merging data by looking for missing data
sum(is.na(OTUt.prop$Project)) # 0 missing data 

```
make data frame into a long format for plotting and statistical analysis 
Here we have the OTUs in columns 15 to the Nth column in the data.frame. This should be edited if you change the input file

```{r}
names(OTUt.prop)
OTU.long <- gather(OTUt.prop, Taxa, proportion, 15:34, factor_key=TRUE)

names(OTU.long)

```

Read in the tax ID - we removed this at the begining before we tranposed the data

```{r}
taxalist <- read.csv("taxalist.csv", header=T)

```
Merge the tax ID with the data based on the OTU ID

```{r}
OTU.long<-merge(OTU.long, taxalist, by.x = "Taxa", by.y = "OTU")
```
TIP: if you have a time series with dates that you want to plot use the following code using the as.Date function 
# > df$Date <- as.Date(df$Date, "%d/%m/%Y") 

```{r}
```

Load some pretty colours

```{r}
colours29 <- c("#404142", "#2f4b7c", "#a05195", "#ff7c43", "#665191", "#e0e084",
               "#c68d23", "#8DD3C7", "#9894ae", "#FB8072", "#80B1D3", "#f95d6a", "#B3DE69", "#FCCDE5", 
               "#D9D9D9", "#BC80BD", "#CCEBC5", "#d45087", "#FFED6F", "#1F78B4", "#33A02C", "#FB9A99", 
               "#FDBF6F", "#FF7F00", "#a51213", "#CAB2D6", "#77a070", "#B15928", "#003f5c")
```


Start plotting data 

1) abundance graph to visualise the abundant taxa across different tree heights.


```{r}
P1<-ggplot(OTU.long,aes(Vertposition,proportion,fill=Order, col=Order)) +
  geom_density( stat = "identity", alpha = 0.2)+
  #geom_line(data = df.data, aes(x = Date2, y= counts, col = Site), linetype = "dashed")+
  xlab("Vertical Hieght (m)")+
  ylab("Proportion of Sequence Reads")+
  #ylim(0,4.5e+06)+
  facet_wrap( . ~ transecttree, ncol=2, scales = "free_y") +
  theme_bw() +   # remove grey background
  #theme(legend.position="none")+
  scale_fill_manual(values= colours29)+
  scale_colour_manual(values= colours29)+
  theme(panel.grid.minor = element_blank())+   # remove minor lines on plot
  theme(panel.grid.major = element_blank())+   # remove major lines on plot
  theme(axis.text.x = element_text(colour="grey20",size=12,angle=0,hjust=.5,vjust=.5,face="plain"),  #horizontal axis text, grey, size 16, no angle etc.
        axis.text.y = element_text(colour="grey20",size=12,angle=0,hjust=1,vjust=0,face="plain"),     #vertical axis text, grey, size 16, no angle etc.
        axis.title.x = element_text(colour="grey20",size=14,angle=0,hjust=.5,vjust=0,face="plain"),   #horizontal axis label, grey, size 20, no angle etc.
        axis.title.y = element_text(colour="grey20",size=14,angle=90,hjust=.5,vjust=.5,face="plain")) #vertical axis label, grey, size 20, no angle etc.

P1

```

2) Have a look at the biological replicates - one option is to use the facet function:

```{r}

P2<-ggplot(OTU.long,aes(Vertposition,proportion,fill=Order, col=Order)) +
  geom_density( stat = "identity", alpha = 0.2)+
  #geom_line(data = df.data, aes(x = Date2, y= counts, col = Site), linetype = "dashed")+
  xlab("Vertical Hiegh (m)")+
  ylab("Proportion of Sequence Reads")+
 # ylim(0,1)+
  facet_grid(transecttree ~ Order, scales = "free") +
  theme_bw() +   # remove grey background
  theme(legend.position="none")+
  theme(strip.text.x = element_text(size=8, angle=75),
          strip.text.y = element_text(size=12, face="bold"),
          strip.background = element_rect(colour="black", fill="#CCCCFF"))+
  scale_fill_manual(values= colours29)+
  scale_colour_manual(values= colours29)+
  theme(panel.grid.minor = element_blank())+   # remove minor lines on plot
  theme(panel.grid.major = element_blank())+   # remove major lines on plot
  theme(axis.text.x = element_text(colour="grey20",size=8,angle=0,hjust=.5,vjust=.5,face="plain"),  #horizontal axis text, grey, size 16, no angle etc.
        axis.text.y = element_text(colour="grey20",size=8,angle=0,hjust=1,vjust=0,face="plain"),     #vertical axis text, grey, size 16, no angle etc.
        axis.title.x = element_text(colour="grey20",size=12,angle=0,hjust=.5,vjust=0,face="plain"),   #horizontal axis label, grey, size 20, no angle etc.
        axis.title.y = element_text(colour="grey20",size=12,angle=90,hjust=.5,vjust=.5,face="plain")) #vertical axis label, grey, size 20, no angle etc.

P2
```


This data contains 10 biological replicates (samples collected from 10 different trees).

Using the summarize function here we calculated the average proportion of sequences, the variance and standard deviation at each height across all biolgical replicates. 


```{r}
## summary statistcs
names(OTU.long)
str(OTU.long)

summary<-ddply(OTU.long,  .(Vertposition, Order), summarize,
               average.proportion = mean(proportion, na.rm=TRUE), sd.proportion = (sd(proportion, na.rm=TRUE)),var.proportion = (var(proportion)))

str(summary)
```


3) Plot summary statistics: average ± the biological variation
Here we can see that most orders found a lower vertical distances have similar proportion across all replicate trees.


```{r}
P3<-ggplot(summary,aes(x = Vertposition, y = average.proportion,fill=Order, col=Order, group = Order)) +
  geom_point(shape = 21, size = 2, alpha = 0.3, position = position_dodge(width = 0.3))+
  #facet_wrap( ~ site, ncol=3, scales = "free_y")+
  geom_errorbar(data = summary, aes(x = Vertposition, y = average.proportion, ymin = average.proportion - var.proportion, ymax = average.proportion + var.proportion, width = 0.7),col = "black", position = position_dodge(width = 0.3)) +
  #geom_errorbar(aes(x = Vertposition, y = average.proportion, ymin = average.proportion - var.proportion, ymax = average.proportion + var.proportion, group = Taxa), col = "black", position = position_dodge(width = 0.3)) +
  xlab("Vertical Height")+
  ylab("Proportion of Sequence reads")+
  theme_bw() +   # remove grey background
  #theme(legend.position="none")+
  theme(panel.grid.minor = element_blank())+   # remove minor lines on plot
  theme(panel.grid.major = element_blank())+   # remove major lines on plot
  theme(axis.text.x = element_text(colour="grey20",size=12,angle=0,hjust=.5,vjust=.5,face="plain"),  #horizontal axis text, grey, size 16, no angle etc.
  axis.text.y = element_text(colour="grey20",size=12,angle=0,hjust=1,vjust=0,face="plain"),     #vertical axis text, grey, size 16, no angle etc.
  axis.title.x = element_text(colour="grey20",size=14,angle=0,hjust=.5,vjust=0,face="plain"),   #horizontal axis label, grey, size 20, no angle etc.
  axis.title.y = element_text(colour="grey20",size=14,angle=90,hjust=.5,vjust=.5,face="plain"))  #vertical axis label, grey, size 20, no angle etc.  

P3

```



Here, even with a subset of the data, we can plot the mean variance relationship which shows a postive relationship between the mean proportion of sequences and the variance of the proportion of sequences across all samples i.e. there are many zeros and many low abundant taxa therefore producing lower variance and few taxa present at high proportions. This also suggests that the whole dataset best fits a negative binomial distribution due to many zeros in the dataset, which is common for abundance data. 


Simple linear model to support the mean and variance relationship:
```{r}
m1<-lm(log(average.proportion) ~ log(var.proportion), data = summary)
anova(m1)
summary(m1)
```
Make a funtion to add the model to the plot
```{r}
lm_eqn = function(m1) {
  
  l <- list(a = format(coef(m1)[1], digits = 2),
            b = format(abs(coef(m1)[2]), digits = 2),
            r2 = format(summary(m1)$r.squared, digits = 2));
  
  if (coef(m1)[2] >= 0)  {
    eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,l)
  } else {
    eq <- substitute(italic(y) == a - b %.% italic(x)*","~~italic(r)^2~"="~r2,l)    
  }
  
  as.character(as.expression(eq));                 
}
```


4) Creat a mean variance plot


```{r}

P4<-ggplot(summary,aes(x = log(average.proportion), y = log(var.proportion),fill=Order, col=Order, group = Order)) +
  geom_point(shape = 21, size = 2, alpha = 0.3)+
  #geom_smooth(method = "lm", se = FALSE)+
  xlab("Mean proportion")+
  ylab("Variation of the proportion")+
   geom_text(aes(x = -6, y = 0, label = lm_eqn(lm(log(average.proportion) ~ log(var.proportion)))), parse = TRUE, col = "black")+
    annotate("text",  x = -7, y = -2, label="p < 0.0001" )+
  theme_bw() +   # remove grey background
  #theme(legend.position="none")+
  theme(panel.grid.minor = element_blank())+   # remove minor lines on plot
  theme(panel.grid.major = element_blank())+   # remove major lines on plot
  theme(axis.text.x = element_text(colour="grey20",size=12,angle=0,hjust=.5,vjust=.5,face="plain"),  #horizontal axis text, grey, size 16, no angle etc.
  axis.text.y = element_text(colour="grey20",size=12,angle=0,hjust=1,vjust=0,face="plain"),     #vertical axis text, grey, size 16, no angle etc.
  axis.title.x = element_text(colour="grey20",size=14,angle=0,hjust=.5,vjust=0,face="plain"),   #horizontal axis label, grey, size 20, no angle etc.
  axis.title.y = element_text(colour="grey20",size=14,angle=90,hjust=.5,vjust=.5,face="plain"))  #vertical axis label, grey, size 20, no angle etc.  

P4
```