Clustering and differential usage of repertoire CDR3 sequences

In CellaRepertorium: Data structures, clustering and testing for single cell immune receptor repertoires (scRNAseq RepSeq/AIRR-seq)

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.

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>")

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.

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.

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

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.

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

cluster_plot(aa80)

Cluster CDR3 DNA sequences

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

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:

germline_cluster = fine_clustering(germline_cluster, sequence_key = 'cdr3_nt', type = 'DNA')
ggplot(germline_cluster$cluster_tbl %>% filter(chain == 'TRB'), aes(x = v_gene, y = j_gene, fill = n_cluster)) + geom_tile() + theme(axis.text.x = element_text(angle = 90))

Number of pairs

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 pearson 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.

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.

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.

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.

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.

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

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.

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)

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

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

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.

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.

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

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.

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.

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:

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:

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:

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

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.

Colophone

sessionInfo()

CellaRepertorium documentation built on Nov. 8, 2020, 7:48 p.m.