---
title: "Understanding RNA-seq Data through Proportionality Analysis"
author: "Thomas Quinn"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Understanding RNA-seq Data through Proportionality Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## Introduction

In this vignette, we use a real dataset to show how we can apply proportionality analysis to understand RNA-seq count data. We place a particular emphasis here on documenting the visualization tools included with this package. Although the user may feel eager to start here, we strongly recommend first reading the companion vignette, "Calculating the Proportionality Coefficients of Compositional Data".

## Big Counts

As a use case, we analyze raw RNA-seq counts from a published study on cane toad evolution and adaptation (Rollins 2015). This dataset contains the transcript counts for 20 toads, sampled from two locations in the wild. Sugar cane farmers introduced cane toads to Australia in 1935 as a pest control measure, but these toads quickly became pests themselves. Starting in Queensland (QLD), they have since spread out into parts of Western Australia (WA). The two locations sampled, which we will treat as the experimental groups, include the early settlement region in QLD and the front of expansion into WA. In this analysis, we try to tease apart the genomic differences between the regional and invasive toads. To begin, we load the `propr` library along with the example data.

```{r, eval = FALSE}
library(propr)
data(caneToad.counts)
data(caneToad.groups)
```

```{r, echo = FALSE}
library(propr)
data(caneToad.groups)
data(top)
```

Next, we calculate proportionality between all transcripts in the dataset. When working with a large number of transcripts (57,580 in this case), we may first want to filter non-informative transcripts to minimize the computational burden of analysis. However, while the variance of the log-ratio transformed feature pair (VLR) (i.e., the numerator portion of $\phi$ and $\rho$) remains fixed regardless of pre-filtering, the individual log-ratio transformed variances (VLS) (i.e., the denominator portion of $\phi$ and $\rho$) do change. In other words, although VLR is sub-compositionally coherent, VLS is not. Therefore, pre-filtering *before* calculating $\rho$ will impact the final result.

That said, in some circumstances, it is simply infeasible to calculate a proportionality matrix using all features. For example, calculating pairwise proportionality in this case requires at least 32 GB of RAM, not counting the additional RAM required for indexing and visualization. However, beginning with version 2.0.3, we can now pre-filter *while* calculating $\rho$, which does not impact the final result. For this pre-filter step, we remove all transcripts that do not have at least 10 counts in at least 10 samples. We base this on the intuition that lowly-expressed transcripts likely have high relative error in cross-sample comparisons due to random noise from the assay technology.

```{r, eval = FALSE}
keep <- apply(caneToad.counts, 2, function(x) sum(x >= 10) >= 10)
```

The proportionality metric $\rho$ has two advantages compared to $\phi$. First, $\rho$ is measured on a scale from $[-1, 1]$, reinforcing its analogy to correlation. Second, $\rho$ adjusts the measure of dependence between two features by the individual variance of *both* features. This symmetrizes the proportionality matrix and may help make proportionality more robust (Erb 2016).

```{r, eval = FALSE}
rho <- perb(caneToad.counts, select = keep)
```

Next, we index the most highly proportional pairs based on an arbitrary threshold. In the absence of any statistical testing framework, we set this threshold at $\rho>0.99$ to include only "very proportional" transcript pairs. Alternatively, we could set this threshold at $\rho<-0.99$ to include only "very unproportional" pairs. However, take note that we use a more stringent threshold here so that the vignette renders more quickly.

```{r, eval = FALSE}
best <- rho[">", .995]
```

## Index-aware Plots

After indexing, it is a good idea to check the pair-wise distribution of the log-ratio transformed data using the `plot` (or, equivalently, `smear`) function. As an "index-aware" function, `plot` only plots feature pairs indexed with `[`. When interpreting this figure, a smear of straight diagonal lines confirms that the pairs are proportional. Intuitively, this means that as one transcript increases in log-ratio transformed expression, so does the other. If you see lines deviating considerably from the $y=x$ diagonal, you likely set your index threshold too low.

```{r, eval = FALSE}
plot(best)
```

```{r, echo = FALSE, results = "hide", fig.show = "hold", fig.keep = "last"}
plot(top)
```

You can also inspect how the indexed pairs cluster using the `dendrogram` function. Specifically, this function clusters features based on a hierarchical clustering of the dissimiliarity measure of the proportionality matrix. In this package, we define the dissimiliarity measure of a proportionality matrix as `as.dist(1-abs(rho@matrix))`. Like `plot`, this function is "index-aware" and only plots feature pairs indexed with `[`. Take note that this dendrogram matches the dendrogram used to label co-clusters in the `prism`, `bokeh`, and `bucket` plots, described below. However, it differs from the dendrogram produced by the `snapshot` plot. Heatmap intensity is not scaled.

```{r, eval = FALSE}
dendrogram(best)
```

```{r, echo = FALSE, results = "hide"}
dendrogram(top)
```

## Index-naive Plots

The remaining visualization methods do not restrict plotting to indexed pairs, but instead incorporate all features present in the `propr` object. To exclude features that do not belong to an indexed feature pair, we `simplify` the indexed `propr` object.

```{r, eval = FALSE}
top <- simplify(best)
```

Compositional data are not distributed in Euclidean space, but rather exist in a space known as the Aitchison simplex (Fernandes 2014). The log-ratio transformation (stored in the `@logratio` slot), used when calculating proportionality, transforms the relative count data from the Aitchison simplex to Euclidean space, allowing us to summarize our data through conventional statistics.

First, we look at `pca`, a function for visualizing the samples across the first two dimensions as reduced by principal components analysis (PCA). This is a valid implementation of PCA plotting for relative count data, built using log-ratio transformed data (Gloor 2016). The `group` argument allows us to color the sample IDs by group membership.

```{r, results = "hide", message = FALSE, fig.keep = "last"}
pca(top, group = caneToad.groups)
```

Second, we look at `snapshot`, a function for visualizing the intensity of the log-ratio transformed count data across samples. Heatmap intensity is not scaled.

```{r, results = "hide", message = FALSE}
snapshot(top)
```

Finally, we look at `prism`, `bokeh`, and `bucket` (pronounced *bouquet*), three functions for visualizing the co-clustering of proportional features. We mention these plot functions together because they share some key similiarities. First, these functions are "index-naive" and plot all $\rho$ in the `@matrix` slot of the `propr` object. Second, these functions identify the feature pairs where both of the constituents co-cluster (with the total number of clusters toggled by the argument `k`). Third, these functions return a vector of cluster memberships for the features appear in the `propr` object.

The `prism` function plots the variance of the ratio of the log-ratio transformed feature pair (VLR) versus the sum of the individual variances of each log-ratio transformed feature (VLS). The ratio of the VLR to the VLS equals $1 - \rho$. As such, we use here seven rainbow colored lines to indicate where $\rho = [.01, .05, .50, 0, 1.50, 1.95, 1.99]$, going from red to violet. A low VLR with a high VLS suggests that the feature pair remains in an equilibrium despite high variability among the individual features ($\rho \approx 1$). A very high VLR with a high VLS suggests that the feature pair remains in an inverted equilibrium despite high variability among the individual features ($\rho \approx -1$).

```{r, results = "hide", message = FALSE}
clusts <- prism(top, k = 5)
```

The `bokeh` function plots pairs across the individual variances of the constituent log-ratio transformed features. For clarity of visualization, the figure projects the data on a log-fold scale. Therefore, the highly variable co-clusters appear in the top-right of the figure while the lowly variable co-clusters appear in the bottom-left. Meanwhile, highly proportional pairs tend to aggregate around the $y=x$ diagonal. The user can retrieve the table used to generate the `prism` and `bokeh` plots by passing the `propr` object through the `slate` function.

```{r, results = "hide", message = FALSE}
clusts <- bokeh(top, k = 5)
```

The `bucket` function (pronounced *bouquet*) plots an estimation of the degree to which a feature pair differentiates the experimental groups versus the proportionality between that pair.

```{r, results = "hide", message = FALSE}
clusts <- bucket(top, group = caneToad.groups, k = 5)
```

These figures help us conceptualize high-dimensional data and select a highly proportional module for further analysis. In this example, we have a particular interest in co-cluster 4 because it (1) shows low proportionality in the setting of high individual feature variance, and (2) appears to differentiate the two experimental groups.

## Down-stream

We can extract co-cluster 4 from the `propr` object using the `subset` method.

```{r, results = "hide"}
sub <- subset(top, select = (clusts == 4))
```

Now, we can use the `pca` functions to see how well this cluster differentiates the two experimental groups based on unsupervised methods. We see below in the `pca` plot that projecting our samples across this highly proportional module leads good separation between the experimental groups. This matches the separation achieved in the source publication which used `edgeR` for feature selection (Rollins 2015).

```{r, results = "hide", fig.keep = "last"}
pca(sub, group = caneToad.groups)
```

Having identified a cluster that separates the experimental groups, the next step in this pipeline might involve a gene set enrichment analysis (GSEA) of the transcripts participating in this highly proportional module. The application of GSEA is beyond the scope of this vignette. However, we show below how to extract the names of the transcripts that belong to this cluster.

```{r}
transcripts <- colnames(sub@logratio)
```

## Conclusions

In place of a highly parameterized "black box" algorithm, we introduce a simple tool that is valid for all biological count data, including the most "difficult-to-normalize" datasets. In the example above, we show how this package can uncover a highly proportional transcript module that accurately differentiates the experimental groups. We find this feat particularly impressive considering that, with the exception of the superfluous `bucket` plot, we discovered this module without ever specifying how the experimental groups should guide feature selection.

Although `propr` currently lacks a hypothesis testing framework, we believe that some degree of "certainty" arises from the fact that we discovered a differentially expressed module using unsupervised methods. Still, in light of this limitation, we emphasize here that data analytics is not the end of scientific inquiry, but rather the beginning. No amount of hypothesis testing can substitute experimental validation in the lab.

## References

1. Erb, Ionas, and Cedric Notredame. “How Should We Measure Proportionality on Relative Gene Expression Data?” Theory in Biosciences = Theorie in Den Biowissenschaften 135, no. 1-2 (June 2016): 21-36. http://dx.doi.org/10.1007/s12064-015-0220-8.

2. Fernandes, Andrew D., Jennifer Ns Reid, Jean M. Macklaim, Thomas A. McMurrough, David R. Edgell, and Gregory B. Gloor. “Unifying the Analysis of High-Throughput Sequencing Datasets: Characterizing RNA-Seq, 16S rRNA Gene Sequencing and Selective Growth Experiments by Compositional Data Analysis.” Microbiome 2 (2014): 15. http://dx.doi.org/10.1186/2049-2618-2-15.

3. Gloor, Gregory B., and Gregor Reid. “Compositional Analysis: A Valid Approach to Analyze Microbiome High-Throughput Sequencing Data.” Canadian Journal of Microbiology 62, no. 8 (August 2016): 692-703. http://dx.doi.org/10.1139/cjm-2015-0821.

4. Lovell, David, Vera Pawlowsky-Glahn, Juan José Egozcue, Samuel Marguerat, and Jürg Bähler. “Proportionality: A Valid Alternative to Correlation for Relative Data.” PLoS Computational Biology 11, no. 3 (March 2015): e1004075. http://dx.doi.org/10.1371/journal.pcbi.1004075.

5. Rollins, Lee A., Mark F. Richardson, and Richard Shine. “A Genetic Perspective on Rapid Evolution in Cane Toads (Rhinella Marina).” Molecular Ecology 24, no. 9 (May 2015): 2264-76. http://dx.doi.org/10.1111/mec.13184.
