---
title: "Clustering and differential usage of repertoire CDR3 sequences"
output: BiocStyle::html_document
author:
- name: Andrew McDavid
  affiliation: University of Rochester, Department of Biostatistics and Computational Biology
  email: Andrew_McDavid@urmc.rochester.edu
vignette: >
  %\VignetteIndexEntry{Clustering and differential usage of repertoire CDR3 sequences}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

In this vignette we demonstrate clustering of 3rd complementary determining region sequence (CDR3) and V-J gene identity of mouse T cells,  ways to visualize and explore clusters that are expanded, pairing of alpha-beta clusters, tests of differential CDR3 usage, and permutation tests for overall clonal properties.

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

```{r setup}
library(CellaRepertorium)
library(dplyr)
library(ggplot2)
library(readr)
library(tidyr)
library(stringr)
library(purrr)
```



# Load filtered contig files

We begin with a `data.frame` of concatenated contig files ('all_contig_annotations.csv'), output from the Cellranger VDJ pipeline.

```{r}
data(contigs_qc)
MIN_CDR3_AA = 6


cdb = ContigCellDB_10XVDJ(contigs_qc, contig_pk = c('barcode', 'pop', 'sample', 'contig_id'), cell_pk = c('barcode', 'pop', 'sample'))
cdb
```

Initially we start with  `r nrow(cdb)` cells and `r nrow(cdb$contig_tbl)` contigs.  We keep contigs that are 

* full - length 
* productive
* high-confidence
* only from T cells
* and with CDR3 sufficiently long.

Then we add a descriptive readable name for each contig.

```{r}
cdb$contig_tbl = dplyr::filter(cdb$contig_tbl, full_length, productive == 'True', high_confidence, chain != 'Multi', str_length(cdr3) > MIN_CDR3_AA) %>% mutate( fancy_name = fancy_name_contigs(., str_c(pop, '_', sample)))

```

After filtering, there are `r nrow(cdb)` cells and `r nrow(cdb$contig_tbl)` contigs.



# Clustering contigs by sequence characteristics

As a first step to define clonotypes, we will first find equivalence classes of CDR3 sequences with the program [CD-HIT](http://weizhongli-lab.org/cdhit_suite/cgi-bin/index.cgi?cmd=cd-hit). In this case, we use the translated amino acid residues, but often one might prefer to use the DNA sequences, by setting the `sequence_key` accordingly and `type = 'DNA'`. Additionally, a higher identity threshold might be appropriate  (see below).

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

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`.  The `min_length` can be set somewhat smaller, but there is a lower limit for the cdhit algorithm.  `G=1`, the default, specifies a global alignment.  This is almost always what is desired, but local alignment is available if `G=0`.

```{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)`.


```{r}
cluster_plot(aa80)
```



## Cluster CDR3 DNA sequences

```{r, results = 'hide'}
cdb = cdhit_ccdb(cdb, 'cdr3_nt', type = 'DNA', cluster_pk = 'DNA97', identity = .965, min_length = MIN_CDR3_AA*3-1, G = 1)
cdb = fine_clustering(cdb, sequence_key = 'cdr3_nt', type = 'DNA')

cluster_plot(cdb)
```

We can also cluster by DNA identity.

## Cluster by V-J identity

```{r}
germline_cluster = cluster_germline(cdb, segment_keys = c('v_gene', 'j_gene', 'chain'), cluster_pk = 'segment_idx')
```

We can cluster by any other feature of the contigs. Here we cluster each contig based on the chain and V-J genes.  This gives us the set of observed V-J pairings:

```{r}
germline_cluster = fine_clustering(germline_cluster, sequence_key = 'cdr3_nt', type = 'DNA')
filter_cdb(germline_cluster, chain == 'TRB') %>% plot_cluster_factors(factors = c('v_gene','j_gene'), statistic = 'contigs', type = 'heatmap')
```

Number of pairs.  The pearson residual (showing the difference from expected counts given marginals) is probably more informative, set `statistic = 'residual'` for this.

```{r}
ggplot(germline_cluster$cluster_tbl %>% filter(chain == 'TRB'), aes(x = v_gene, y = j_gene, fill = avg_distance)) + geom_tile() + theme(axis.text.x = element_text(angle = 90))
```

Average Levenshtein distance of CDR3 within each pair.  This might be turned into a z-score by fitting a weighted linear model with sum-to-zero contrasts and returning the studentized residuals.  This could determine if a pairing has an unexpected small, or large, within cluster distance.



## Expanded clusters

Next, we will examine the clusters that are found in many contigs.  First we will get a canonical contig to represent each cluster.  This will be the medoid contig, by default.

```{r}
aa80 = canonicalize_cluster(aa80, representative = 'cdr3', contig_fields = c('cdr3', 'cdr3_nt', 'chain', 'v_gene', 'd_gene', 'j_gene'))
```

`aa80` now includes the fields listed in `contig_fields` in the `cluster_tbl`, using the values found in the medoid contig.

```{r}
MIN_OLIGO = 7
oligo_clusters = filter(aa80$cluster_tbl, n_cluster >= MIN_OLIGO)
oligo_contigs = aa80
oligo_contigs$contig_tbl = semi_join(oligo_contigs$contig_tbl, oligo_clusters, by = 'aa80')
oligo_contigs
```

Get contigs/cells/clusters found at least `r MIN_OLIGO` times (across contigs).  Note that replacing `contig_tbl` with the subset selected with the `semi_join` also automatically subsetted the `cell_tbl` and `cluster_tbl`.


```{r}
oligo_clusters = oligo_contigs$contig_tbl %>% group_by(aa80) %>% summarize(`n subjects observed` = length(unique(sample))) %>% left_join(oligo_clusters)

knitr::kable(oligo_clusters %>% select(aa80:cdr3, chain:j_gene, avg_distance, n_cluster))

```

Report some statistics about these expanded clusters, such as how often they are found, how many subjects, etc.

```{r}
oligo_plot = ggplot(oligo_contigs$contig_tbl, aes(x = representative, fill = chain)) + geom_bar() + coord_flip() + scale_fill_brewer(type = 'qual') + theme_minimal()
oligo_plot
```

These always come from a single chain.

```{r}
oligo_plot + aes(fill =   sample) + facet_wrap(~pop)

```

But come from multiple populations and samples.

## Some simple phylogenetic relationships

By using the within-cluster distances, some rudamentory plots attempting to show phylogenetic associations are possible.  (These are most biologically appropriate for B cells that undergo somatic hypermutation.)  

```{r}
library(ggdendro)

dendro_plot = function(ccdb, idx, method = 'complete'){
    h = filter(ccdb$cluster_tbl, !!sym(ccdb$cluster_pk) == idx) %>% pull(fc) %>% .[[1]]
    quer = filter(ccdb$contig_tbl, !!sym(ccdb$cluster_pk) == idx)
    hc = hclust(as.dist(h$distance_mat), method = method) %>% dendro_data(type = "rectangle")
    hc$labels = cbind(hc$labels, quer)
   ggplot(hc$segments, aes(x=x, y=y)) + geom_segment(aes(xend=xend, yend=yend)) + 
  theme_classic() + geom_text(data = hc$labels, aes(color = sample, label = fancy_name), size = 3, angle = 60, hjust =0, vjust = 0) + scale_x_continuous(breaks = NULL) + ylab('AA Distance') + xlab('')
}

to_plot = aa80$cluster_tbl %>% filter(min_rank(-n_cluster) == 1)

map(to_plot$aa80, ~ dendro_plot(aa80, .))

```

A full-blown generative model of clonal generation and selection would be recommended for any actual analysis, but these plots may suffice to get a quick idea of the phylogenetic structure.

## Formal testing for frequency differences

We can test for differential usage of a clone, or cluster with `cluster_logistic_test` and `cluster_test_by`.  The latter splits the `cluster_tbl` by `field = 'chain'`, thereby adjusting the number of cell trials included in the "denominator" of the logistic regression. 
The formula tests for differences between populations, including the sample as a random effect, and only tests clusters that are included in the `oligo_clusters` set.


```{r,  results = 'hide'}
mm_out = cluster_test_by(aa80, fields = 'chain', tbl = 'cluster_tbl', formula = ~ pop + (1|sample), filterset = cluster_filterset(white_list = oligo_clusters)) %>%
  left_join(oligo_clusters)

mm_out = mutate(mm_out, conf.low = estimate-1.96*std.error, 
                conf.high = estimate + 1.96*std.error)

```

```{r per_iso_tests}
mm_outj = filter(ungroup(mm_out), term == 'popbalbc') %>% arrange(desc(representative))

ggplot(mm_outj, aes(x = representative, ymin = conf.low, ymax = conf.high, y = estimate)) + geom_pointrange()  + coord_flip() + theme_minimal() + geom_hline(yintercept = 0, lty = 2) + xlab("Isomorph") + ylab("log odds of isomorph")
```


We test if the binomial rate of clone expression differs between balbc and b6, for the selected clones.  None appear to be different.

## Length of CDR3

```{r}
aa80$contig_tbl = aa80$contig_tbl %>% mutate(cdr3_length = str_length(cdr3_nt))
ggplot(aa80$contig_tbl, aes(fill = pop, x= cdr3_length)) +
  geom_histogram(binwidth = 1, mapping = aes(y = ..density..)) + 
  theme_minimal() + scale_fill_brewer(type = 'qual') + 
  facet_grid(sample ~chain) + theme(strip.text.y = element_text(angle = 0)) + coord_cartesian(xlim = c(25, 55))

```

Some authors have noted that the length of the CDR3 region can be predictive of T cell differentiation. In our study, there doesn't appear to be a noticeable difference between BALB/c and C57BL/6J (b6) mice, but if we needed to make sure, an appropriate procedure would be to run a mixed model with a random `sample` effect (assumed to represent a biological replicate).

```{r cdr3_len, fig.width = 3, fig.height = 3}
cdr_len = aa80$contig_tbl %>% group_by(chain) %>% do(broom::tidy(lme4::lmer(cdr3_length ~ pop + (1|sample), data = .), conf.int = TRUE))
ggplot(cdr_len %>% filter(term == 'popbalbc'), aes(x = interaction(chain, term), y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange() + theme_minimal() + coord_flip() + 
  ylab('Length(CDR3 Nt)') + xlab('Term/Chain') + geom_hline(yintercept = 0, lty = 2)

```

We end up with a (harmless) convergence warning about a singular fit.  This is expected, because the `samples` aren't actually replicates -- they are just subsamples drawn for illustrative purposes.
The Balbc mice  have .5 fewer nucleotides per contig, on average,  and this is not significant.


# Clonal pairs

Next, we can examine the pairing between $\alpha-\beta$ chains and see if any pairs are found more than once.

```{r expanded_clones}
aa80$cluster_pk = 'representative'
aa80 = rank_prevalence_ccdb(aa80)
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'))


```

`pairing_tables` finds all contig combinations of order `table_order` across cells.  Among those combinations that occur at least `min_expansion` times, the expanded combinations and and any other combinations that shared an expanded combo.


```{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

```

## Expanded clones

```{r}
whitelist = oligo_clusters %>% dplyr::select(cluster_idx.1 = representative) %>% unique()
pairing_list = pairing_tables(aa80, table_order = 2, orphan_level = 1, min_expansion = Inf, cluster_whitelist = whitelist,  cluster_keys = c('cdr3', 'representative', 'chain', 'v_gene', 'j_gene', 'avg_distance'))
<<plot_expanded>>
```

By setting `min_expansion = Inf, cluster_whitelist = whitelist` we can examine any pairings for a set of cluster_idx, in this case the ones that were seen multiple times.  Interestingly (and unlike some human samples) the expanded clusters are  $\beta$-chain, and their $\alpha$ chains are sprinkled quite evenly across clusters.

# Permutation tests

Permutation tests allow tests of independence between cluster assignments and other cell-level covariates (such as the sample from which the cell was derived).  The cluster label is permuted to break the link between cell and cluster, and an arbitrary statistic of both cluster label, and cell covariate is evaluated.

```{r}
aa80_chain = split_cdb(aa80, 'chain') %>% lapply(canonicalize_cell, contig_fields = 'aa80')

compare_expanded = function(cluster_idx, grp){
  # cluster_idx contains the permuted cluster assignments
  # grp the cell_covariate_keys.
  # NB: this is always a data.frame even if it is just a single column
  # cross tab by pop
  tab = table(cluster_idx, grp[[1]])
  # count number of times an aa80 class was expanded
  expanded = colSums(tab>=2)
  # compare difference
  expanded['b6'] - expanded['balbc']
}
```

The signature of the statistic should be of a vector `cluster_idx` and `data.frame`.

```{r}
set.seed(1234)
perm1 = cluster_permute_test(aa80_chain$TRB, cell_covariate_keys = 'pop', cell_label_key = 'aa80', n_perm = 100, statistic = compare_expanded)

perm1
```

Although b6 mice had `r perm1$observed` more clones observed to be expanded (occuring >=2 times) than balbc, this is not signficant under a null model where cells were permuted between mouse types (populations), where b6 are expected to have about `r round(perm1$expected)` more expanded clones, just due to the additional number of cells sampled in b6 and the particular spectrum of clonal frequencies in this experiment:

```{r}
knitr::kable(table(pop = aa80_chain$TRB$pop))
```

Indeed if we resample in a way that fixes each group to have the same number of cells:
```{r}
rarify = aa80_chain$TRB$cell_tbl %>% group_by(pop) %>% do(slice_sample(., n = 377))

aa80_chain$TRB$cell_tbl = semi_join(aa80_chain$TRB$cell_tbl, rarify)

cluster_permute_test(aa80_chain$TRB, cell_covariate_keys = 'pop', cell_label_key = 'aa80', n_perm = 500, statistic = compare_expanded)
```

We see that this discrepacy between the number of expanded clones between subpopulations is  mostly explained by a greater number of cells sampled in b6, but also random variability plays a role.

We can also test for oligoclonality, eg, how often is a beta chain expanded in a sample:

```{r}
count_expanded = function(cluster_idx, grp){
  # clusters x sample contigency table
  tab = table(cluster_idx, grp[[1]])
  # number of cluster x samples that occured more than once
  expanded = sum(tab>1)
  expanded
}

perm3 = cluster_permute_test(aa80_chain$TRB,  cell_covariate_keys = 'sample', cell_label_key = 'aa80', n_perm = 500, statistic = count_expanded)
perm3
```

`r perm3$observed` expanded clones were observed in each of the two populations vs `r round(perm3$expected)` expected, and this discrepancy would be significant at $p<$ `r ceiling(perm3$p.value*100)/100`.  This is indicating that there is underdispersion -- fewer clusters are expanded than expected, given the spectrum of clonal frequencies and the number of cells per sample.

To further elucidate this, we can restrict the permutations to maintain certain margins of the table by specifying `cell_stratify_keys.`  This doesn't effect the observed values of the statistics, but will change the expected values (since these are now conditional expectations.)  Here we restrict the permutations within levels of `pop` (eg, only permuting within balbc, and within b6).

```{r}
cluster_permute_test(aa80_chain$TRB,   cell_covariate_keys = 'sample', cell_stratify_keys = 'pop', cell_label_key = 'aa80', n_perm = 500, statistic = count_expanded)
```

In the restricted permutations, the expected number of expanded clusters is even greater.  Both of these effects are due to the fact that the "sample" replicates, within each population actually are not biological replicates, which inflates the `cluster_idx` margin of the table.

## Sequences of permutation tests across cell subpopulations

In many cases, we want test a sequence of contrasts of `cell_covariate_keys` variables vs the `cell_label_key`. For instance, `cell_covariate_keys` might include the cell subpopulation derived from gene expression taken from a cell, and we want to compare the levels of clonal expansion between subpopulations.  Suppose we have four such subpopulation `ident`s:

```{r}
ident = gl(length = length(aa80_chain$TRB$cell_tbl$pop), n = 4, k = 1)
head(ident)
aa80_chain$TRB$cell_tbl$ident = ident
```

We can compare all pairs: 2 vs 1, 3 vs 1, 4 vs 1, 3 vs 2, 4 vs 2, 4 vs 3, and in fact this is the default action.  If `contrasts` (a matrix or list of vectors) is specified, we can perform other sets of comparisons, like 1 vs the average in all other `ident`:
```{r}
contrast_vec = c(1, -1/3, -1/3, -1/3)
```

The key is that the `statistic` should return a vector, as this now does:

```{r}
compare_expanded_vec = function(cluster_idx, grp){
  tab = table(cluster_idx, grp[[1]])
  # count number of times an aa80 class was expanded
  expanded = colSums(tab>=2)
  expanded
}
```


```{r perm-pairwise}
perm4 = cluster_permute_test(aa80_chain$TRB, cell_covariate_keys = 'ident', cell_stratify_keys = 'sample', 
                             cell_label_key = 'aa80', n_perm = 100, statistic = compare_expanded_vec)

plot_permute_test(perm4)
```

The default of all pairwise comparisons.

```{r perm-avg}
perm5 = cluster_permute_test(aa80_chain$TRB, cell_covariate_keys = 'ident', cell_stratify_keys = 'sample', 
                             cell_label_key = 'aa80', n_perm = 100, statistic = compare_expanded_vec, contrasts = contrast_vec)

plot_permute_test(perm5)
```

```{r}
tidy.PermuteTestList(perm5)
```

There is also a `tidy` method.

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