05-explore.Rmd

```{r, echo=FALSE}
library(mia)
library(miaViz)
library(dplyr)
tse <- readRDS("data/Tengeler2020/tse.rds")
```



# Microbiome data exploration

Now we have loaded the data set into R. Next, let us walk through some
basic operations for data exploration to confirm that the data has all
the necessary components.

Chapters 4-6 in
[OMA](https://microbiome.github.io/OMA/datamanipulation.html) provide
many more examples on exploring and manipulating microbiome data.


## Data structure

Let us now investigate how taxonomic profiling data is organized in R.

Dimensionality tells us how many taxa and samples the data
contains. As we can see, there are `r dim(tse)[1]` taxa and `r dim(tse)[2]`
samples.

```{r}
dim(tse)
```

The `rowData` slot contains a taxonomic table. This includes taxonomic
information for each of the `r nrow(tse)` entries. With the `head()`
command, we can print just the beginning of the table.

The `rowData` seems to contain information from `r ncol(rowData(tse))`
different taxonomy classes.

```{r}
knitr::kable(head(rowData(tse))) %>% 
  kableExtra::kable_styling("striped", 
                            latex_options="scale_down") %>% 
  kableExtra::scroll_box(width = "100%")
```

The colData slot contains sample metadata. It contains information for all `r ncol(tse)` samples.
However, here only the 6 first samples are shown as we use the `head()` command. There
are `r ncol(colData(tse))` columns, that contain information, e.g., about patients' status, and cohort.

```{r}
knitr::kable(head(colData(tse))) %>% 
  kableExtra::kable_styling("striped", 
                            latex_options="scale_down") %>% 
  kableExtra::scroll_box(width = "100%")
```

From here, we can draw summaries of the sample (column) data, for
instance to see what is the patient status distribution.

The command `colData(tse)$patient_status` fetches the data from the
column, and `table()` creates a table that shows how many times each
class is present, and `sort()` sorts the table to ascending order.

There are `r length(colData(tse)$patient_status[colData(tse)$patient_status == "ADHD"])` 
samples from patients having ADHD, 
and `r length(colData(tse)$patient_status[colData(tse)$patient_status == "Control"])` control samples.

```{r}
sort(table(colData(tse)$patient_status))
```




### Transformations

Microbial abundances are typically 'compositional' (relative) in the
current microbiome profiling data sets. This is due to technical
aspects of the data generation process (see e.g. [Gloor et al.,
2017](https://www.frontiersin.org/articles/10.3389/fmicb.2017.02224/full)).

The next example calculates relative abundances as these are usually easier to
interpret than plain counts. For some statistical models we need to
transform the data into other formats as explained in above link (and
as we will see later).

```{r}
# Calculates relative abundances, and stores the table to assays
tse <- transformCounts(tse, method = "relabundance")
```

A variety of standard transformations for microbiome data are available for `TSE` data objects through [mia R package](https://microbiome.github.io/mia/reference/transformCounts.html).


### Aggregation

Microbial species can be called at multiple taxonomic resolutions. We
can easily agglomerate the data based on taxonomic ranks. Here, we
agglomerate the data at Phylum level.

```{r 05explore_agglomerateByRank}
tse_phylum <- agglomerateByRank(tse, rank = "Phylum")

# Show dimensionality
dim(tse_phylum)
```



Now there are `r dim(tse_phylum)[1]` taxa and `r dim(tse_phylum)[2]`
samples, meaning that there are `r dim(tse_phylum)[1]` different
Phylum level taxonomic groups. Looking at the `rowData` after
agglomeration shows all Firmicutes are combined together, and all
lower rank information is lost.

From the assay we can see that all abundances of taxa that belong to
Firmicutes are summed up.

```{r}
knitr::kable(head(rowData(tse_phylum))) %>% 
  kableExtra::kable_styling("striped", 
                            latex_options="scale_down") %>% 
  kableExtra::scroll_box(width = "100%")
```

If you are sharp, you have by now noticed that all the aggregated
values in the above example are NA's (missing data). This is because
the agglomeration is missing abundances for certain taxa, and in that
case the sum is not defined by default (`na.rm = FALSE`). We can
ignore the missing values in summing up the data by setting `na.rm =
TRUE`; then the taxa that do not have information in specified level
will be removed. Those taxa that do not have information in specified
level are agglomerated at lowest possible level that is left after
agglomeration.

```{r}
temp <- rowData(agglomerateByRank(tse, rank = "Genus"))

# Prints those taxa that do not have information at the Genus level
knitr::kable(head(temp[temp$Genus == "",])) %>% 
  kableExtra::kable_styling("striped", 
                            latex_options="scale_down") %>% 
  kableExtra::scroll_box(width = "100%")
```

Here agglomeration is done similarly, but na.rm = TRUE

```{r}
temp2 <- rowData(agglomerateByRank(tse, rank = "Genus", na.rm = TRUE))

print(paste0("Agglomeration with na.rm = FALSE: ", dim(temp)[1], " taxa."))
print(paste0("Agglomeration with na.rm = TRUE: ", dim(temp2)[1], " taxa."))
```

The [mia
package](https://microbiome.github.io/mia/reference/index.html)
contains further examples on various data agglomeration and splitting
options.

## Visualization

The [miaViz package](https://microbiome.github.io/miaViz/) facilitates
data visualization. Let us plot the Phylum level abundances.

```{r}

# Here we specify "relabundance" to be abundance table that we use for plotting.
# Note that we can use agglomerated or non-agglomerated tse as an input, because
# the function agglomeration is built-in option. 

# Legend does not fit into picture, so its height is reduced.
plot_abundance <- plotAbundance(tse, abund_values="relabundance", rank = "Phylum") +
  theme(legend.key.height = unit(0.5, "cm")) +
  scale_y_continuous(label = scales::percent)

plot_abundance 
```

**Density plot** shows the overall abundance distribution for a given
  taxonomic group. Let us check the relative abundance of Firmicutes
  across the sample collection. The density plot is a smoothened
  version of a standard histogram.

The plot shows peak abundances around 30 %.

```{r}
# Subset data by taking only Firmicutes
tse_firmicutes <- tse_phylum["Firmicutes"]

# Gets the abundance table
abundance_firmicutes <- assay(tse_firmicutes, "relabundance")

# Creates a data frame object, where first column includes abundances
firmicutes_abund_df <- as.data.frame(t(abundance_firmicutes))
# Rename the first and only column
colnames(firmicutes_abund_df) <- "abund"

# Creates a plot. Parameters inside feom_density are optional. With 
# geom_density(bw=1000), it is possible to adjust bandwidth.
firmicutes_abund_plot <- ggplot(firmicutes_abund_df, aes(x = abund)) + 
  geom_density(color="darkred", fill="lightblue") + 
  labs(x = "Relative abundance", title = "Firmicutes") +
  theme_classic() + # Changes the background
  scale_x_continuous(label = scales::percent)

firmicutes_abund_plot

```

```{r, echo=FALSE}
# # Does the same thing but differently
# # Calculates the density. Bandwidth can be adjusted; here, it is 0.065.
# # density() is from stats package
# density_firmicutes <- density(abundance_firmicutes, bw = 0.065)
# 
# # Plots the density
# plot(density_firmicutes,
#      xlab="Relative abundance",
#      ylab="Density",
#      main=paste0("Firmicutes (",density_firmicutes$n, " obs, ", 
#      density_firmicutes$bw, " bw)"))
```


For more visualization options and examples, see the [miaViz vignette](https://microbiome.github.io/miaViz/articles/miaViz.html).



## Exercises (optional)

Explore some of the following questions on your own by following
[online examples](https://microbiome.github.io/OMA/). Prepare a
reproducible report (Rmarkdown), and include the code that you use to
import the data and generate the analyses.

 * **Abundance table** Retrieve the taxonomic abundance table from the
   example data set (TSE object). Tip: check "assays" in [data import
   section](https://microbiome.github.io/OMA/data-introduction.html#loading-experimental-microbiome-data)
   
 * How many different samples and genus-level groups this phyloseq
   object has? Tips: see dim(), rowData()

 * What is the maximum abundance of Akkermansia in this data set? Tip:
   aggregate the data to Genus level with agglomerateByRank, pick
   abundance assay, and check a given genus (row) in the assay

 * Draw a histogram of library sizes (total number of reads per
   sample). Tip: Library size section in
   [OMA](https://microbiome.github.io/OMA/quality-control.html). You
   can use the available function, or count the sum of reads per
   sample by using the colSums command applied on the abundance
   table. Check [Vandeputte et
   al. 2017](https://www.nature.com/articles/nature24460) for further
   discussion on the differences between absolute and relative
   quantification of microbial abundances.

 * **Taxonomy table** Retrieve the taxonomy table and print out the
   first few lines of it with the R command head(). Investigate how
   many different phylum-level groups this phyloseq object has? Tips:
   rowData, taxonomicRanks in
   [OMA](https://microbiome.github.io/OMA/taxonomic-information.html#functions-to-access-taxonomic-information).

 * **Sample metadata** Retrieve sample metadata. How many patient
     groups this data set has? Draw a histogram of sample
     diversities. Tips: colData

 * **Subsetting** Pick a subset of the data object including only
     ADHD individuals from Cohort 1. How many there are? Tips: subsetting in [OMA](https://microbiome.github.io/OMA/datamanipulation.html#subsetting)

 * **Transformations** The data contains read counts. We can convert
  these into relative abundances and other formats. Compare abundance
  of a given taxonomic group using the example data before and after
  the compositionality transformation (with a cross-plot, for
  instance). You can also compare the results to CLR-transformed data
  (see e.g. [Gloor et
  al. 2017](https://www.frontiersin.org/articles/10.3389/fmicb.2017.02224/full))

 * **Visual exploration** Visualize the population distribution of
   abundances for certain taxonomic groups. Do the same for
   CLR-transformed abundances. Tip: assays, transformCounts
   
 * Experiment with other data manipulation tools from
   [OMA](https://microbiome.github.io/OMA/taxonomic-information.html#functions-to-access-taxonomic-information).





## Exploring data: example solutions

 * **Abundance table** Retrieve the taxonomic abundance table from the
   example data set (TSE object).
   
```{r}
# We show only part of it
assays(tse)$counts[1:6,1:6]
```
 * How many different samples and genus-level groups this phyloseq
   object has?

```{r}
dim(rowData(tse))
```

 * What is the maximum abundance of Akkermansia in this data set?

```{r}
# Agglomerating to Genus level
tse_genus <- agglomerateByRank(tse,rank="Genus")
# Retrieving the count for Akkermansia
Akkermansia_abund <- assays(tse_genus)$count["Genus:Akkermansia",]
max(Akkermansia_abund)
```
 * Draw a histogram of library sizes (total number of reads per
   sample).

```{r}
library(scater)
tse <- addPerCellQC(tse) # adding a new column "sum" to colData, of total counts/sample.
hist(colData(tse)$sum, xlab = "Total number of reads per sample", main = "Histogram of library sizes")

```
 * **Taxonomy table** Retrieve the taxonomy table and print out the
   first few lines of it with the R command head(). Investigate how
   many different phylum-level groups this phyloseq object has?

```{r}
head(rowData(tse))
taxonomyRanks(tse) # The taxonomic ranks available
unique(rowData(tse)["Phylum"]) # phylum-level groups
length(unique(rowData(tse)["Phylum"])[,1]) # number of phylum-level groups
```

 * **Sample metadata** Retrieve sample metadata. How many patient
     groups this data set has? Draw a histogram of sample
     diversities.
     
```{r}
colData(tse) # samples metadata
unique(colData(tse)$patient_status) # patient groups
# Example of sample diversity using Shannon index
tse <- mia::estimateDiversity(tse, 
                             abund_values = "counts",
                             index = "shannon", 
                             name = "shannon")
hist(colData(tse)$shannon, xlab = "Shannon index", main = "Histogram of sample diversity")
```

 * **Subsetting** Pick a subset of the data object including only
     ADHD individuals from Cohort 1. How many there are?

```{r}
sub_cohort_1 <- tse[, colData(tse)$cohort=="Cohort_1"]
sub_cohort_1_ADHD <- sub_cohort_1[, colData(sub_cohort_1)$patient_status=="ADHD"]
colData(sub_cohort_1_ADHD)
```

 * **Transformations** The data contains read counts. We can convert
  these into relative abundances and other formats. Compare abundance
  of a given taxonomic group using the example data before and after
  the compositionality transformation (with a cross-plot, for
  instance). You can also compare the results to CLR-transformed data
  (see e.g. [Gloor et
  al. 2017](https://www.frontiersin.org/articles/10.3389/fmicb.2017.02224/full))
  
```{r}
tse <- transformCounts(tse, method = "relabundance")
tse <- transformCounts(tse, method = "clr", abund_values = "counts",pseudocount = 1)
# Lets compare with taxa: A29
taxa <- "A29"
df <- as.data.frame(list(
                          counts=assays(tse)$counts[,taxa],
                          relabundance=assays(tse)$relabundance[,taxa],
                          clr=assays(tse)$clr[,taxa])
                        )
ggplot(df, aes(x=counts,y=relabundance))+
  geom_point()+
  geom_smooth()

ggplot(df, aes(x=counts,y=clr))+
  geom_point()+
  geom_smooth()

ggplot(df, aes(x=relabundance,y=clr))+
  geom_point()+
  geom_smooth()
```
 * **Visual exploration** Visualize the population distribution of
   abundances for certain taxonomic groups. Do the same for
   CLR-transformed abundances.
   
``` {r}
# Same taxa used as earlier
ggplot(df, aes(x=counts,colour="blue")) +
    geom_density(alpha=.2)+
  theme(legend.position = "none")+
  labs(title =paste("Distribution of counts for",taxa, collapse = ": "))
ggplot(df, aes(x=clr,colour="red")) + 
    geom_density(alpha=.2)+
  theme(legend.position = "none")+
  labs(title =paste("Distribution of clr for",taxa, collapse = ": "))
```