powerTCR.Rmd

---
title: "powerTCR"
author: "Hillary Koch"
package: powerTCR
abstract: >
    This vignette walks through the functionality of the powerTCR package for analyzing the T cell receptor repertoire. The package's main goals are to fit models to the clone size distribution of the TCR repertoire, and to perform model-based comparative analysis of samples, following the work of Koch et al. (2018). For user simplicity, powerTCR can handle data in the formats MiTCR, MiTCR with UMIs, MiGEC, VDJtools, ImmunoSEQ (both old and new), MiXCR, and IMSEQ. If not using one of these formats, powerTCR only requires that the user be able to turn a sample TCR repertoire into a vector of clone sizes.
output: BiocStyle::html_document
vignette: > 
  %\VignetteIndexEntry{Vignette Title}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}  
---

# Introduction
The powerTCR package allows users to implement the model-based methods discussed in our Koch et al. (2018) paper in PLoS Computational Biology. Specifically, the clone size distribution of the T cell receptor (TCR) repertoire exhibits imperfect power law behavior; powerTCR supports a model that keeps this fact in mind. Additionally, powerTCR contains tools to fit another power law model for the TCR repertoire detailed in Desponds et al. (2016). Given a collection of sampled TCR repertoires, powerTCR equips the user with tools for comparative analysis of the samples, using one of two model-based approaches. This leads to hierarchical clustering of the samples to determine their relatedness based on the clone size distribution alone.

## Summary of features

- Read in and parse TCR sequencing files stored in various formats into the necessary format for powerTCR
- Fit two power law-based models to the clone size distribution of the TCR repertoire
    1. The discrete gamma-GPD spliced threshold model of Koch et al. (2018)
    2. The type-I Pareto model of Desponds et al. (2016)
- Compare sample TCR repertoire samples based on model fit with hierarchical clustering
- Simulate data according from the gamma-GPD spliced threshold distribution

# Fitting a model

In order to fit a model with powerTCR, you only need to be able to supply a __vector of counts__ (that is, a vector of clone sizes). If your data are in a format supported by `parseFile` or `parseFolder`, you can simply read in your file using one of those functions, specify whether or not you want to use only in-frame sequences, and powerTCR will automatically give you a sorted vector of clone sizes for each sample. This functionality is a wrapper for parsing functions found in the `tcR` package.

powerTCR contains a toy data set, called `repertoires`, with two TCR repertoire samples, which we will use throughout this vignette. (In practice, you may have any number of samples.) You can load powerTCR and this data set by typing:

```{r, echo = 4:5, message = FALSE, warning = FALSE}
# This installs packages from BioConductor
# if (!requireNamespace("BiocManager", quietly=TRUE))
    # install.packages("BiocManager")
# BiocManager::install("powerTCR")
library(powerTCR)
data("repertoires")
```

`repertoires` is a list with 2 elements in it, each corresponding to a sample repertoire. Have a look:

```{r}
str(repertoires)
```

These samples are smaller than one might expect a TCR repertoire to be in practice, but for the sake of exploring powerTCR, they permit much faster computation. 

## The discrete gamma-GPD spliced threshold model

The main model that powerTCR focuses on is the discrete gamma-GPD spliced threshold model. This distribution has probability mass function

\[
    f(x) =
    \begin{cases}
        (1-\phi)\frac{h(x|\boldsymbol{\theta}_b)} {H(u-1|\boldsymbol{\theta}_b)} & \text{for $x \leq u-1$} \\
        \phi g(x|\boldsymbol{\theta}_t, u) & \text{for $x \geq u$}
    \end{cases},
\]

where $h$ and $H$ are the density and distribution function of a gamma distribution, and $g$ is the density of a generalized Pareto distribution, or GPD. The gamma distribution has density

\[
    h(x) = \frac{\beta^\alpha}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}
\]

and the GPD has density

\[
    g(x) = \frac{1}{\sigma}\big(1+\xi \frac{x-u}{\sigma}\big)^{-(1/\xi +1)}.
\]

We can fit the model to each of the samples in `repertoires` using the function `fdiscgammagpd`. This function takes a few arguments. The most important are as follows:

First, `fdiscgammagpd` needs to be passed a sample TCR repertoire as a vector of counts. Second, you need to specify a grid of possible thresholds (that is, the parameter $u$) that you are interested in considering. One easy way to do this might be to specify a series of quantiles of the vector of counts. Finally, you also need to specify the shift, which for each sample is ideally the smallest count (at least for TCR repertoire samples). The shift is the minimum value in the support of the distribution, and for clone sizes, should never be smaller than 1.

Let's try fitting the model to the data in `repertoires`.


```{r, warning = FALSE}
# This will loop through our list of sample repertoires,
# and store a fit in each
fits <- list()
for(i in seq_along(repertoires)){
    # Choose a sequence of possible u for your model fit
    # Ideally, you want to search a lot of thresholds, but for quick
    # computation, we are only going to test 4
    thresholds <- unique(round(quantile(repertoires[[i]], c(.75,.8,.85,.9))))
    
    fits[[i]] <- fdiscgammagpd(repertoires[[i]], useq = thresholds,
                               shift = min(repertoires[[i]]))
}
names(fits) <- names(repertoires)
```

The output for a fit looks like this:
```{r}
# You could also look at the first sample by typing fits[[1]]
fits$samp1
```

Each value of the output is described in the `fdiscgammagpd` help file, but the most important are

* nllh: the negative log likelihood of the most likely fit, given the thresholds you've checked
* mle: the maximum likelihood estimates for:

\[\phi, \alpha, \beta, u, \sigma,\text{ and } \xi\text{ respectively}\]

These two important items can be grabbed using convenient accessors contained in powerTCR, called `get_mle` and `get_nllh`.

```{r}
# Grab mles of fits:
get_mle(fits)

# Grab negative log likelihoods of fits
get_nllh(fits)
```

You can also view the likelihoods for every other threshold you checked (in nllhuseq) as well as the output from `optim` for the "bulk"  (truncated gamma) and "tail" (GPD) parts of the distribution.

## The Type-I Pareto model

For reproducibility purposes, powerTCR also provides a means to fit the model of Desponds et al. (2016). This model is investigated and discussed in Koch et al. (2018). The model follows a type-I Pareto distribution, with density:

\[
f(x) = \frac{\alpha u^\alpha}{x^{\alpha+1}}.
\]

For a given threshold $u$, the estimate for parameter $\alpha$ is computed directly as 

\[
\alpha=n\bigg[\sum_{i=1}^n\text{log}\frac{x_i}{u}\bigg]^{-1}+1
\]

where $n$ is the number of clones with size larger than $u$. This value is computed for every possible threshold $u$, and then the parameters that minimize the KS-statistic between empirical and theoretical distributions are chosen.

Let's fit this model to the `repertoires` data, and have a look at the output for the first sample.

```{r}
desponds_fits <- list()
for(i in seq_along(repertoires)){
    desponds_fits[[i]] <- fdesponds(repertoires[[i]])
}
names(desponds_fits) <- names(repertoires)
desponds_fits$samp1
```

Here, min.KS is the minimum KS-statistic of all possible fits. Cmin is the threshold $u$ that corresponds to the best fit. powerlaw.exponent and pareto.alpha are effectively the same -- pareto.alpha = powerlaw.exponent-1. This is just user preference; for the Pareto density given above, pareto.alpha corresponds to the $\alpha$ shown there. However, if the user is more familiar with a "power law" distribution, then powerlaw.exponent is the parameter they should look at.

# Density, distribution, and quantile functions, plus simulating data
powerTCR provides standard functions to compute the density, distribution, and quantile functions of the discrete gamma-GPD threshold model, as well a function to simulate data. These can be very useful for tasks such as visualizing model fit and conducting a simulation study. The functions behave exactly like popularly used functions such as, say, `dnorm, pnorm, qnorm`, and `rnorm`. In order to use these functions, you need to specify all of the model parameters. The one exception is $\phi$, which can go unspecified -- details about how $\phi$ defaults are in the help file for `ddiscgammagpd`.

Here, we will use `qdiscgammagpd` to compute quantiles from the two theoretical distributions we fit above. Note that we pass `qdiscgammagpd` the quantiles and the fits we obtained. This is just something convenient powerTCR can do if you don't want to manually specify shift, shape, rate, u, sigma, xi, and phi. You don't have to pass a fit, however, as you will see when we simulate data.

```{r}
# The number of clones in each sample
n1 <- length(repertoires[[1]])
n2 <- length(repertoires[[2]])

# Grids of quantiles to check
# (you want the same number of points as were observed in the sample)
q1 <- seq(n1/(n1+1), 1/(n1+1), length.out = n1)
q2 <- seq(n2/(n2+1), 1/(n2+1), length.out = n2)

# Compute the value of fitted distributions at grid of quantiles
theor1 <- qdiscgammagpd(q1, fits[[1]])
theor2 <- qdiscgammagpd(q2, fits[[2]])
```

Now, let's visualize the fitted and empirical distributions by plotting them together. Here, the black represents the original data, with the quantiles of the theoretical distributions plotted on top in color.

```{r, fig.wide=TRUE, echo=2:7}
par(mfrow = c(1,2))
plot(log(repertoires[[1]]), log(seq_len(n1)), pch = 16, cex = 2,
     xlab = "log clone size", ylab = "log rank", main = "samp1")
points(log(theor1), log(seq_len(n1)), pch = 'x', col = "darkcyan")

plot(log(repertoires[[2]]), log(seq_len(n2)), pch = 16, cex = 2,
     xlab = "log clone size", ylab = "log rank", main = "samp2")
points(log(theor2), log(seq_len(n2)), pch = 'x', col = "chocolate")
```

The fits look pretty good!

Let's also try simulating data.

```{r}
# Simulate 3 sampled repertoires
set.seed(123)
s1 <- rdiscgammagpd(1000, shape = 3, rate = .15, u = 25, sigma = 15,
                    xi = .5, shift = 1)
s2 <- rdiscgammagpd(1000, shape = 3.1, rate = .14, u = 26, sigma = 15,
                    xi = .6, shift = 1)
s3 <- rdiscgammagpd(1000, shape = 10, rate = .3, u = 45, sigma = 20,
                    xi = .7, shift = 1)
```

__NB:__ it is possible to simulate data according to a distribution that is totally unrealistic. For example, what if you chose a very light-tailed gamma distribution and a comparatively very high threshold, but insisted (using $\phi$) that data be observed above the threshold? Here is what happens:

```{r}
bad <- rdiscgammagpd(1000, shape = 1, rate = 2, u = 25, sigma = 10,
                     xi = .5, shift = 1, phi = .2)
plot(log(sort(bad, decreasing = TRUE)), log(seq_len(1000)), pch = 16,
     xlab = "log clone size", ylab = "log rank", main = "bad simulation")
```

Fun, but not too realistic for a clone size distribution. There are several ways to go about finding reasonable parameters to simulate. One intuitive and easy technique is to let real data speak for itself -- use parameters similar to those obtained by fitting a distribution to true TCR repertoire data sets.

# Doing comparative analysis

Following the work in Koch et al. (2018), powerTCR provides the tools needed to perform hierarchical clustering of TCR repertoire samples according to their Jensen-Shannon distance. We can test this out on the 3 TCR repertoires we just simulated. First, we need to fit a model to them. For computational efficiency, let's just supply the true thresholds. Then, we can use `JS_dist` to compute the Jensen-Shannon divergence between each pair of theoretical distributions corresponding to each of the TCR samples.

`JS_dist` needs to be supplied two model fits from `fdiscgammagpd` as well as a grid. The grid is important: it is the range over which each distribution gets evaluated. If you are comparing a group of TCR repertoires, the minimum value of your grid should be the smallest clone size across all samples. The upper bound of the grid should be something very large, say 100,000 or more. If you don't select a value large enough, you will not be examining the tail of your fitted distributions sufficiently, and the tail is important! The grid should also contain every integer between its minimum and maximum. For computational efficiency, here the upper bound on our grid is only 10,000.

We've wrapped `JS_dist` up into a convenient function `get_distances`, which will create a symmetric matrix of distances (with 0 on the diagonal) for your use.

```{r}
# Fit model to the data at the true thresholds
sim_fits <- list("s1" = fdiscgammagpd(s1, useq = 25),
                 "s2" = fdiscgammagpd(s2, useq = 26),
                 "s3" = fdiscgammagpd(s3, useq = 45))

# Compute the pairwise JS distance between 3 fitted models
grid <- min(c(s1,s2,s3)):10000
distances <- get_distances(sim_fits, grid, modelType="Spliced")
```

Let's have a look at the distance matrix we just computed:
```{r}
distances
```

Note that `get_distances` is just calling `JS_dist` for every pair, and `JS_dist` is simply a wrapper for two functions called `JS_spliced` and `JS_desponds`. If you want to do comparative analysis for data fit using `fdesponds`, then the argument "modelType" must be changed to "Desponds".

We can use this distance matrix to perform hierarchical clustering. This is done easily with the `clusterPlot` function. `clusterPlot` is just a wrapper for `hclust` in R's `stats` package, and takes a matrix of Jensen-Shannon distances like the one we just made, plus a type of linkage. All possible types of linkage are listed in the help file, but we recommend using Ward's method or complete linkage.

```{r, fig.small=TRUE}
clusterPlot(distances, method = "ward.D")
```

The clustering result is exactly what we might expect. Indeed, we simulated s1 and s2 using very similar parameter settings, so we should expect them to be more closely related to each other than to s3. That is exactly what the dendrogram displays.

# Extracting diversity estimators

In Koch et al. (2018), we introduced new measures of diversity, alongise comparisons to often-used estimators of sample diversity borrowed from ecology literature. Providing a list of model fits, powerTCR can easily compute for you sample richness, Shannon entropy, clonality, and our introduced measure -- the proportion of highly stimulated clones. This is done with the function `get_diversity`. Let's demonstrate this on our fits from simulated data:

```{r, echo=2:3, message = FALSE, warning = FALSE}
library(vegan)
get_diversity(sim_fits)
```

# Bootstrapping model fits

Finally, powerTCR allows you to run a parametric bootstrapping procedure on the discrete Gamma-GPD spliced threshold model, and can be executed in parallel if the processors are available to you. Just supply the fits you'd like to bootstrap and the number of resamples you'd like to do (and the number of cores you want to use, if running in parallel). If you want to speed up the bootstrapping by choosing a reduced number of thresholds, you can do so by changing `gridStyle`. Setting `gridStyle` to the default "copy" will just copy the original `useq` parameter in the fit. More details are in the package documentation. For each fit being bootstrapped, this returns a list of fits of length `resamples`. So, for example, if you are bootstrapping 5 fits with 1,000 resamples each, the function will return a list of length 5. Each of those lists will be length 1,000.

Here, we set the number of resamples to something very small to save computational time, but we recommend 1,000 resamples to get reasonable confidence bands.

```{r}
boot <- get_bootstraps(fits, resamples = 5, cores = 1, gridStyle = "copy")
```

Finally, you can get confidence intervals based on your resampling. For example, this is how you can get the 95% confidence interval around the estimate for $\xi$ from the first fit.
```{r, echo = 3:5, message = FALSE, warning = FALSE}
library(magrittr)
library(purrr)
mles <- get_mle(boot[[1]])
xi_CI <- map(mles, 'xi') %>% 
    unlist %>%
    quantile(c(.025,.5,.975))
xi_CI
```



# References

Desponds, J., Mora, T., & Walczak, A. M. (2016). Fluctuating fitness shapes the clone-size distribution of immune repertoires. Proceedings of the National Academy of Sciences, 113(2), 274-279.
Chicago