extended_example.Rmd

---
title: "Plotting the distribution of taxa"
author: "Sur Herrera Paredes"
date: "`r Sys.Date()`"
output:
  github_document:
    toc: TRUE
---

# Introduction

In this extended example I go through every step to produce a relative
abundance barplot that represents bacterial communities living in
individual hosts. 

Bacterial communities are everywhere, and when we characterize them
it is important to describe they overall taxonomic structure as that
gives us clues as to what types of functions might be performed by the
community. At the same time, it is important to show the variability
of these communities and thus it is useful to plot them at the lowest
aggregation level possible.

In this example, I utilize data from a big experiment that
was published [here](https://www.nature.com/articles/nature11237). In
that experiment, we planted individual *Arabidopsis thaliana* plants in
individual pots. The pots each had one of two types of natural soil,
each from two different seasons. We planted eight different accessions in
each of those soils, and plants were harvested at two developmental stages.
Additionally we had unplanted soil only pots, these are the "soil" samples
in the example. For each individual plant, we harvested two fractions
(i.e. E & R), one which we called the Endophytic Compartment (E) and
corresponds to the interior of the root after removing the outer cell wall,
and another which we called Rhizosphere (R) which is the soil within 1mm of
the plant root. So **E** samples contains bacteria inside the root, and **R**
samples contain bacteria immediately surrounding the root.

Plant-bacteria interactions in the root are incredibly important because
the root is both the gut and the brain of the plant. Microbes there
can benefit of the products of the plant photosynthesis as sources
of nutrition, and can also provide chemistry that the plant couldn't perform
by itself. However, the story is more complicated because microbial competition
and the plant immune system also provide a fertile evolutionary environment 
for antagonistic interactions.

Ultimately understanding and being able to manipulate plant-bacteria
interactions has a lot of implications as hunger is one of the most pressing
problems of humanity with 800 million people living in hunger. Agriculture
is our only sustainable tool against hunger, and even though it currently 
employs a quarter of the World population it has not been enough to tackle
this challenge. 

# Getting ready

First you need to get the data. If you haven't check the 
[README](https://github.com/surh/scip_barplot/blob/master/README.md) file
of the GitHub repository of the workshop. It is also recommended that
you watch the YouTube [video](https://www.youtube.com/watch?v=siIoupAnILk)
that runs through the example. Finally, you will need to install the
`tidyverse` package.

Once you have everything you need, start an R session and load the tidyverse
package:

```{r}
library(tidyverse)
```
 
# Read data

First read the OTU table. You may need to change the file path
to wherever you downloaded the files in your machine.

```{r}
Tab <- read_tsv("data/rhizo/otu_table.tsv")
Tab
```

The code above reads the file into a `tibble`, which is a type
of `data.frame` that has some neat additional properties. You
don't need to concern yourself too much with the differences. 

The code above also produces a warning, indicating that `read_tsv`
tried to guess the types of data in each column of the table. It
guessed correctly but you should always specify the expected columns
with the option `col_types` (use `?read_tsv` for additional details).


```{r}
Tab <- read_tsv("data/rhizo/otu_table.tsv",
                col_types = cols(otu_id = col_character(),
                                 .default = col_number()))
Tab
```

# Basic barplot 

We need to think back to the original figure and reformat our data to have one
column for the x-axis and another for the y-axis. This is a requirement for
`ggplot2`. We can to that with `pivot_longer`, a function of the `tidyverse`.

```{r}
Tab %>%
  pivot_longer(-otu_id, names_to = "sample_id", values_to = "count")
```

In the code above the options `samples_to` and `names_to` indicate the
names of the new columns in the new tibble.

Lets create a smaller subset of the data to make some basic plots

```{r}
dat <- Tab %>%
  pivot_longer(-otu_id, names_to = "sample_id", values_to = "count") %>%
  filter(otu_id %in% c("OTU_14834", "OTU_18567", "OTU_14402", "OTU_14822"))
print(dat)
```

Now we create a basic barplot that shows relative abundances of the 4
selected OTUs
```{r}
dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  geom_bar(aes(fill = otu_id), stat = "identity", position = "fill")
```

Now we re-make the plot with some beautification

```{r}
dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  geom_bar(aes(fill = otu_id), stat = "identity", position = "fill") +
  scale_y_continuous(name = "Relative abundance", labels = scales::percent) +
  theme(axis.text.x = element_text(angle = 90))
```

# Adding sample metadata

Read sample metadata

```{r}
Meta <- read_tsv("data/rhizo/sample_metadata.tsv",
                 col_types = cols(.default = col_character()))
Meta
```

We **join** our long format dataset with the sample metadata table

```{r}
print(dat)
dat <- dat %>%
  left_join(Meta, by = "sample_id")
dat
```
We add facts with `facet_grid`

```{r}
p1 <- dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  geom_bar(aes(fill = otu_id), stat = "identity", position = "fill") +
  facet_grid(~ fraction) +
  scale_y_continuous(name = "Relative abundance", labels = scales::percent) +
  theme(axis.text.x = element_text(angle = 90),
        axis.text.y = element_text(color = "black"))
p1
```
We set `scales="free_x"` to plot a different set of samples on each
facet.

```{r}
p1 <- dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  geom_bar(aes(fill = otu_id), stat = "identity", position = "fill") +
  facet_grid(~ fraction, scales = "free_x") +
  scale_y_continuous(name = "Relative abundance", labels = scales::percent) +
  theme(axis.text.x = element_text(angle = 90),
        axis.text.y = element_text(color = "black"))
p1
```

We repeat the plot with some beautification

```{r}
p1 <- dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  geom_bar(aes(fill = otu_id), stat = "identity", position = "fill", width = 1) +
  facet_grid(~ fraction, scales = "free_x") +
  scale_y_continuous(name = "Relative abundance", labels = scales::percent) +
  theme(axis.text.x = element_text(angle = 90),
        axis.text.y = element_text(color = "black"),
        strip.text = element_text(face = "bold"),
        strip.background = element_blank())
p1
```
We can add an extra variable to the facet

```{r}
p1 <- dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  geom_bar(aes(fill = otu_id), stat = "identity", position = "fill", width = 1) +
  facet_grid(~ fraction + soil, scales = "free_x") +
  scale_y_continuous(name = "Relative abundance", labels = scales::percent) +
  theme(axis.text.x = element_text(angle = 90),
        axis.text.y = element_text(color = "black"),
        strip.text = element_text(face = "bold"),
        strip.background = element_blank())
p1
```

```{r}
p1 <- dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  geom_bar(aes(fill = otu_id), stat = "identity", position = "fill", width = 1) +
  facet_grid(~ fraction + soil, scales = "free_x", space = "free_x") +
  scale_y_continuous(name = "Relative abundance", labels = scales::percent) +
  theme(axis.text.x = element_text(angle = 90),
        axis.text.y = element_text(color = "black"),
        strip.text = element_text(face = "bold"),
        strip.background = element_blank())
p1
```

## Excercise 1

Re-make the plot above but using accession instead of soil type in the facets

## Excercise 2

Re-make the plot above but show all three variables (fraction, accession and
soil) in the facets. What happens if you change the order of the terms
in the facet formulas? is any of the orders better?

# Adding OTU taxonomy

```{r}
Tax <- read_tsv("data/rhizo/otu_taxonomy.tsv",
                col_types = cols(otu_id = col_character(),
                                 Phylum = col_character()))
Tax
```
Let's include all OTUs in our data table now

```{r}
dat <-Tab %>%
  pivot_longer(-otu_id, names_to = "sample_id",
               values_to = "count")
print(dat)
```

As we did before with sample metadata, we use `left_join` to add OTU
taxonomic information to our data table

```{r}
dat <- dat %>%
  left_join(Tax, by = "otu_id")
print(dat)
```

Let's make our basic plot, but using `Phylum` instead of `otu_id`

```{r}
p1 <- dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  geom_bar(aes(fill = Phylum), stat = "identity", position = "fill", width = 1) +
  scale_y_continuous(name = "Relative abundance", labels = scales::percent) +
  theme(axis.text.x = element_text(angle = 90),
        axis.text.y = element_text(face = "bold")) 
p1
```

Can we find better colors

```{r}
p1 <- dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  geom_bar(aes(fill = Phylum), stat = "identity", position = "fill", width = 1) +
  scale_fill_brewer(palette = "Paired") +
  scale_y_continuous(name = "Relative abundance", labels = scales::percent) +
  theme(axis.text.x = element_text(angle = 90),
        axis.text.y = element_text(face = "bold")) 
p1
```
I want phyla to be ordered by abundance and "unclassified" to be
at the end

Code to get order not shown but results are in `mean_freqs`

```{r, echo=FALSE}
mean_freqs <- dat %>%
  group_by(sample_id, Phylum) %>%
  summarise(phyl_count = sum(count),
            .groups = 'drop') %>%
  group_by(sample_id) %>%
  summarise(phyl_freq = phyl_count / sum(phyl_count),
            Phylum = Phylum,
            .groups = 'drop') %>%
  group_by(Phylum) %>%
  summarise(mean = mean(phyl_freq),
            .groups = 'drop') %>%
  arrange(desc(mean))
print(mean_freqs)
```
Create a vector with the desired order of phyla and use it to convert
phylum column to a factor.

```{r}
phyla_order <- c("Proteobacteria",                                                                                                                      
                 "Actinobacteria",                                                                                                        
                 "Bacteroidetes",                                                                                        
                 "Acidobacteria",                                                                       
                 "Firmicutes",                                                         
                 "Cyanobacteria",                                     
                 "Verrucomicrobia",                 
                 "Gemmatimonadetes",                                             
                 "Armatimonadetes",
                 "Chloroflexi",
                 "unclassified")

dat <- dat %>%
  mutate(Phylum = factor(Phylum, levels = phyla_order))
dat

```


```{r}
p1 <- dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  geom_bar(aes(fill = Phylum), stat = "identity", position = "fill", width = 1) +
  scale_fill_brewer(palette = "Paired") +
  scale_y_continuous(name = "Relative abundance", labels = scales::percent) +
  theme(axis.text.x = element_text(angle = 90),
        axis.text.y = element_text(face = "bold")) 
p1
```
## Excercise

Use `scale_color_manual` to manually select a good set of colors for this plot

# Combining sample metadata and OTU taxonomy

We can use `lef_join` twice in a row to get a combinde data table with all
the information

```{r}
dat <- Tab %>%
  pivot_longer(-otu_id, names_to = "sample_id", values_to = "count") %>%
  left_join(Meta, by = "sample_id") %>%
  left_join(Tax, by = "otu_id")
print(dat)
```

To get the same order as above we convert the Phylum column to a factor

```{r}
dat <- dat %>%
  mutate(Phylum = factor(Phylum, levels = phyla_order))
dat
```


Now we can combine facets and phylum color and ordering in just one plot

```{r}
p1 <- dat %>%
  ggplot(aes(x = sample_id, y = count)) +
  facet_grid(~ fraction + soil, scales = "free_x", space = "free_x") +
  geom_bar(aes(fill = Phylum), stat = "identity", position = "fill", width = 1) +
  scale_y_continuous(name = "Relative abundance", labels = scales::percent) +
  scale_fill_brewer(palette = "Paired") +
  theme(axis.text.x = element_text(angle = 90, size = 6),
        axis.text.y = element_text(color = "black"),
        strip.text = element_text(face = "bold"),
        strip.background = element_blank())
p1
```


We save the plot to a file

```{r}
ggsave("rhizo_phylo_distribution.png", p1, width = 8, height = 4)
```



# Extra excercises

Look at the files at [data/hmp_v13](data/hmp_v13) which contain much bigger
data tables generated from the Human Microbiome Project (HMP).

Can you make similar plots illustrating the bacterial taxonomic distributions
of the Stool and Saliva? Are there any differences by sex?

Can you order the samples by Proteobacteria relative abundance in the
rhizosphere dataset? what about the Bacteroidetes abundance in the
HMP dataset? (HINT: Use the same approach we used to sort taxonomic groups)


# Session info

```{r}
sessionInfo()
```