---
title: "An Introduction to CellaRepertorium"
output: BiocStyle::html_document
package: CellaRepertorium
vignette: >
  %\VignetteIndexEntry{An Introduction to CellaRepertorium}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
author:
- name: Andrew McDavid
  affiliation: University of Rochester, Department of Biostatistics and Computational Biology
  email: Andrew_McDavid@urmc.rochester.edu
bibliography: citations.bib
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>")
```

Congratulations on your installation of CellaRepertorium.
This package contains methods for manipulating, clustering, pairing, testing and conducting multimodal analysis single cell RepSeq data, especially as generated by [10X Genomics Chromium Immune Profiling](https://support.10xgenomics.com/single-cell-vdj/software/overview/welcome).

# Ethos and data structure

The fundamental unit this package operates on is the **contig**, which
is a section of contiguously stitched reads from a **single cell**. Each
contig belongs to one (and only one) cell, however, cells may generate 
multiple contigs. 


```{r, echo = FALSE}
knitr::include_graphics(system.file('help/figures/contig_schematic.png', package = 'CellaRepertorium'))
```

Contigs belong to cells, and can also belong to a **cluster**.  A
`ContigCellDB()` object tracks these two types of membership by using a sequence of
three `data.frames` (`dplyr::tibble()`, actually). `ContigCellDB()` also
tracks columns (the primary keys) that uniquely identify each row in
each of these tables. The `contig_tbl` is the `tibble` containing
**contigs**, the `cell_tbl` contains the **cells**, and the
`cluster_tbl` contains the **clusters**.


```{r,  echo = FALSE}
knitr::include_graphics(system.file('help/figures/table_schematic.png', package = 'CellaRepertorium'))
```


The `contig_pk`, `cell_pk` and `cluster_pk` specify the columns that
identify a contig, cell and cluster, respectively. These will serve as
foreign keys that link the three tables together.

# Manipulation

We'll start *in media res* with an example of minimially-processed annotated contigs from 
Cellranger.  For details on how to import your own CellRanger data, as well QC steps that could be performed, see `vignette('mouse_tcell_qc')`.

```{r}
library(CellaRepertorium)
library(dplyr)
data("contigs_qc")
cdb = ContigCellDB_10XVDJ(contigs_qc, 
                          contig_pk = c('barcode', 'pop', 'sample', 'contig_id'), 
                          cell_pk = c('barcode', 'pop', 'sample'))
```

This constructs a `ContigCellDB` object, specifying that the columns `barcode`, `pop`, `sample`, and `contig_id` unique identify a contig, so are its **primary keys**, and that columns `barcode`, `pop`, `sample` are the cell primary keys.

We can manipulate the `contig_tbl` with the `$` operator.
```{r}
cdb$contig_tbl$cdr_nt_len = nchar(cdb$contig_tbl$cdr3_nt)
```

Or with the `mutate_cdb` function, which saves a few keystrokes.
```{r}
suppressPackageStartupMessages(library(Biostrings))
cdb = cdb %>% mutate_cdb(cdr3_g_content = alphabetFrequency(DNAStringSet(cdr3_nt))[,'G'], tbl = 'contig_tbl')

head(cdb$contig_tbl, n = 4) %>% 
    select(contig_id, cdr3_nt, cdr_nt_len, cdr3_g_content)
```

Other functionality, some of which is depicted in  `vignette('mouse_tcell_qc')` and 
`vignette('cdr3_clustering')` includes:

*  Exploring pairing (classical alpha-beta or heavy-light, single chain, or >1 chain type) with `enumerate_pairing`.
*  Splitting by a factor with `split_cdb`.
*  Filtering with `filter_cdb`.
*  Canonicalization with `canonicalize_cell`.  This chooses a representative contig for each cell and copies various fields into the `cell_tbl` so that the contig-cell relationship in these fields is now one-to-one.  This is useful for any analysis of contigs that requires the cell as its base.  Likewise  `canonicalize_cluster` choses a representative contig for each cluster.

# Clustering

We provide an R port of the fast biostring clustering algorithm [CD-HIT](http://weizhongli-lab.org/cdhit_suite/cgi-bin/index.cgi?cmd=cd-hit) [@Fu2012].

```{r}
aa80 = cdhit_ccdb(cdb, 'cdr3', type = 'AA', cluster_pk = 'aa80', 
                  identity = .8, min_length = 5)
aa80 = fine_clustering(aa80, sequence_key = 'cdr3', type = 'AA')
```

This partitions sequences into sets with >80% mutual similarity in the amino acid sequence,  adds some additional information about the clustering, and returns it as a `ContigCellDB` object named `aa80`.  The primary key for the clusters is `r aa80$cluster_pk`.

```{r}
head(aa80$cluster_tbl)
head(aa80$contig_tbl) %>% select(contig_id, aa80, is_medoid, `d(medoid)`)
```

The `cluster_tbl`  lists the `r nrow(aa80$cluster_tbl)` 80% identity groups found, including the number of contigs in the cluster, and the average distance between elements in the group.
In the `contig_tbl`, there are two columns specifying if the contig `is_medoid`, that is, is the most representative element of the set and the distance to the medoid element `d(medoid)`.

Other functionality to operate on clustering include:

* `cluster_germline()` defines clusters using combinations of factors in the `contig_tbl`, such as the V- and J-gene identities.

# Pairing

One of the main benefits of single cell repertoire sequencing is the ability to recover both light and heavy chains, or alpha and beta chains of B cells and T cells.  Pairing is a property of the `cell_tbl`.  We provide a number of tools to analyze and visualize the pairing.

```{r}
library(ggplot2)

paired_chain = enumerate_pairing(cdb, chain_recode_fun = 'guess')

ggplot(paired_chain, aes(x = interaction(sample, pop), fill = pairing)) + 
    geom_bar() + facet_wrap(~canonical, scale = 'free_x') + 
    coord_flip() + theme_minimal()
```

We first determine how often cells were paired, and how often *non-canonical* multi-alpha or multi-beta cells are found.  
T cells with two alpha chains, and to a lessor extent, two beta chains are an established phenomena [@Padovan1993; @Padovan1995; @He2002].
However, if the rate of these so-called dual TCR cells is too high, then multiplets or excess ambient RNA may be suspected.  (B plasma cells seem to be particularly sticky, and are laden with immunoglobulin RNA).

## Enumerating clonotypically expanded pairs

Besides this high-level description of the rate and characteristics of the pairing, we want to discover pairs that occur repeatedly.  For that, we can use the `pairing_tables` function.  

First we copy some info into the `cluster_tbl` from the medoid contig:
```{r}
aa80 = canonicalize_cluster(aa80, representative = 'cdr3', 
                            contig_fields = c('cdr3', 'cdr3_nt', 'chain', 'v_gene', 'd_gene', 'j_gene'))
aa80$cluster_pk = 'representative'

```

The `representative` just gives the clusters a more useful unique name (the CDR3 animo acid sequence).  The other information would be helpful with visualizing and understanding the results of the pairing.  

Next we provide an ordering for each contig in a cell:
```{r}
aa80 = rank_chain_ccdb(aa80)
```

In this case the contigs will be picked by the beta/alpha chain of the cluster.  Other options are possible, for instance the prevalence of the cluster with `rank_prevalence_ccdb`, which would allow detection to see expanded, dual-TCR pairings (alpha-alpha or beta-beta).

Finally, we generate the pairings.
```{r}
pairing_list = pairing_tables(aa80, table_order = 2, orphan_level = 1, min_expansion = 3, cluster_keys = c('cdr3', 'representative', 'chain', 'v_gene', 'j_gene', 'avg_distance'))

```

By default, this is subset to a list suitable for plotting in a heatmap-like format, so includes only expanded pairings.  You can get all pairings by setting `min_expansion = 1`, and force inclusion or exclusion of particular clusters with `cluster_whitelist` or `cluster_blacklist`.

```{r plot_expanded}
pairs_plt = ggplot(pairing_list$cell_tbl, aes(x = cluster_idx.1_fct, y = cluster_idx.2_fct)) + 
  geom_jitter(aes(color = sample, shape = pop), width = .2, height = .2) + 
  theme_minimal() + xlab('TRB') + ylab('TRA') + 
  theme(axis.text.x = element_text(angle = 45))

pairs_plt = map_axis_labels(pairs_plt, pairing_list$idx1_tbl, pairing_list$idx2_tbl, aes_label  = 'chain')
pairs_plt

```

# Testing

If the experiment that generated the data was a designed experiment, it might be of interest to test clusters for differential abundance.  We implement ordinary and mixed-effect logistic and binomial tests.  See  `cluster_logistic_test` for details and `vignette('cdr3_clustering')` for an example.

We also implement a permutation test which may be suitable for testing for various phylogenetic properties of clonotypes, such as the clonotype diversity or polarization.  Cluster assignments are permutated by cells, possibly conditioning on less granular covariates.  See `cluster_permute_test` and `vignette('cdr3_clustering')`.

# Multimodal Analysis

`ContigCellDB` objects can be included as a field on the `colData` of a `SingleCellExperiment`.  This permits various multimodal analyses, while maintaining the correspondence between the cells in the `SinglCellExperiment` and the cells in the `ContigCellDB`.  See `vignette('repertoire_and_expression')`.

# Acknowledgments

Development of CellaRepertorium was funded in part by UL1 TR002001 (PI Bennet/Zand) pilot to Andrew McDavid.

# Colophone
```{r}
sessionInfo()
```

# References