The huge diversity of T and B cell receptors is critical for the function of the adaptive immune system: recognizing and eliminating all kinds of antigens, from viruses, to cancer cells. However, this diversity makes it also difficult to understand B and T cell receptor repertoires (or collectively: immune receptor repertoires, IRRs).
You have already learned about some techniques that are useful for IRR analysis when we discussed microbiomes and B cell receptor repertoires (e.g. normalized RADs, Ig gene usage). These techniques can be applied to T cell repertoires, too. Today we will explore further techniques that can deliver a detailed picture of an IRR, namely clustering and graph analysis. These techniques exploit the similarity of TCR sequences:
In our analysis we will ask several questions related to similarities within and between TCR repertoires:
In our analyses we will use a mix of visualization and quantification.
To install ClustIRR please open RStudio and enter [this might take a couple of minutes].
library(devtools)
install_github("snaketron/ClustIRR", upgrade = "never")… after you have installed ClustIRR please load it by entering:
library(ClustIRR)Make sure that ClustIRR version >1.9.39 is loaded. Click on the tab
Packages (bottom right corner in RStudio) search for ClustIRR.
We also need to load other R-packages. Do this next:
library(ggplot2)
library(ggforce)
library(igraph)
library(patchwork)
library(ggrepel)
library(rBLAST)We have two single cell datasets of TCR\(\alpha\beta\) antigen-stimulated repertoires downloaded from the website of Parse Biosciences1 Performance of Evercode™ TCR v4 in Antigen-Stimulated Human T Cells, https://www.parsebiosciences.com/datasets/performance-of-evercode-tcr-v4-in- antigen-stimulated-human-t-cells/; Parse Biosciences; Accessed 2026 June 30:
Load the data:
# load the data
data(D_CMV_FLU)The TCR repertoires are stored in s and the corresponding metadata in
meta. Access them now:
# TCR repertoires
s <- D_CMV_FLU$s
# Meta data
meta <- D_CMV_FLU$metaLet’s peek into the data frame containing both TCR repertoires
str(s)… and also into the metadata
str(meta)Problem: How many clonotypes are in each repertoire? How many CDR3\(\beta\) sequences are in the data frame? How many CDR3\(\alpha\) sequences are in the data frame? How many cells are in each repertoire?
# number of clonotypes
table(s$sample)
# number of unique CDR3b sequences
length(unique(s$CDR3b))
# number of unique CDR3b sequences
length(unique(s$CDR3a))
# count number of cells per samples
print(aggregate(clone_size~sample, data = s, FUN = sum))Problem: Do clone sizes vary significantly within the repertoires?
Create a sina plot to explore the distributions.
ggplot(data = s)+
geom_sina(aes(x = sample, y = clone_size), shape = 21)+
theme_bw()+
scale_y_log10()+
annotation_logticks(base = 10, sides = "l")There are many more clonotypes in the TCR repertoire Flu. Hence we see
many more dots at x=Flu in the sina plot. Overall, we see similar
distributions: most clone sizes fall between 1 and 10 cells. However,
each repertoire also includes a few large clonotypes with more than
1,000 cells (visible at the top). Both sina plots use a log10-scaled
vertical axis.
In addition to single cell TCRs, we also have the gene expression
profiles of individual cells stored in gex.
# Single cell gene expression data
gex <- D_CMV_FLU$gexstr(gex)Problem: Are all ~20,000 human genes included in gex? If not, do
you see common functions among these genes?
Hint: If you do not recognize these gene symbols run gene ontology (GO) enrichment analysis (e.g., by providing these gene symbols as inputs to gProfiler https://biit.cs.ut.ee/gprofiler/gost)
Important remark: Rows in s are TCR clonotypes (groups of clonally
expanded T cells with common ancestors), while rows in gex are T
cells. How can we integrate these two datasets?
Hint: Remember cell barcodes from the lesson on single cell RNA-seq?
Use the cell_id column to access the cell IDs obtained based on these
barcodes.
We have two repertoires, CMV and Flu. We have seen that several TCR clonotypes are expanded in each repertoire. Such clonotype expansions have to be expected as a consequence of the experimental stimulation by CMV and Influenza M1 antigens. With ClustIRR we can detect the effect of expansions (or contraction) on the structure of TCR repertoires, or in general IRR repertoires.
Let’s use ClustIRR to identify how the structure of the TCR repertoire changes from CMV to Flu. ClustIRR operates on CDR3 sequences and clone sizes:
o <- clustirr(s = s, meta = meta,
control = list(blast_gmi = 0.9)) # min. BLAST sequence identitystr(o)The output in o is quite rich. See the help page for explanations.
Check whether you can interpret it!
One of the main outputs generated by clustirr is a joint graph (\(J\)),
stored as an igraph object named graph.
is_igraph(o$graph)Each node in the graph represents a TCR clonotype from either CMV or Flu. Node widths indicate clone sizes, and edge weights represent CDR3\(\alpha\) and CDR3\(\beta\) sequence similarities.
Problem: Visualize \(J\) as an interactive graph using the function
plot_graph! Explore the structure and try coloring the nodes based on
antigen specificity. Do you see any nodes corresponding to CMV-specific
TCR clonotypes? Are these clonotypes typically small or large?
plot_graph(g = o, as_visnet = T, node_opacity = 0.75, show_singletons = F)Important remark: Even in this small case study–where the TCR repertoire contains only about 1,700 TCR clonotypes–the joint graph is already quite complex. This level of complexity makes qualitative analyses impractical. To address this, ClustIRR emphasizes quantitative analyses (see Step 2).
igraph functions to inspect \(J\)We can use the R-package igraph to explore various properties of \(J\).
For example, we can examine the edge weights stored in the edge
attribute nweight. These weights represent alignment scores between
pairs of CDR3\(\beta\) sequences, calculated using BLOSUM62 substitution
scores with a gap opening penalty of -12 and a gap extension penalty of
-2. To ensure comparability between pairs of CDR3 sequences of different
lengths, ClustIRR normalizes the alignment scores by dividing them by
the alignment lengths. Consequently, identical or similar CDR3 sequence
pairs receive large positive nweight scores, while dissimilar pairs
have nweight values close to zero or even negative.
Problem: How many nodes are there in the graph? The igraph
function as_data_frame(x = o$graph, what = "vertices") allows you to
access the node attributes. The output is a data.frame with nodes shown
as rows. How many rows are there?
n <- as_data_frame(x = o$graph, what = "vertices")
nrow(n)Problem: Now access the edges by replacing what = "vertices" with
what = "edges". How many edges are there? Find which pairs of nodes
are connected by the largest/lowest nweight scores.
# get a data.frame 'e' with of edges and their attributes
e <- as_data_frame(x = o$graph, what = "edges")
# max edge weight
e_max <- e[which.max(e$nweight),]
# get CDR3 sequences
n$CDR3b[which(n$name %in% c(e_max$from, e_max$to))]# min edge weight
e_min <- e[which.min(e$nweight),]
# get CDR3 sequences
n$CDR3b[which(n$name %in% c(e_min$from, e_min$to))]ClustIRR employs graph-based community detection (GCD) algorithms–such as Leiden, Louvain, and InfoMap–to identify communities of nodes characterized by a high density of strong edges within the community and a low density of strong edges connecting to nodes outside the community.
Do community detection with the function detect_communities now! Use
similarities based only on the more diverse CDR3\(\beta\) chain for
clustering.
gcd <- detect_communities(graph = o$graph,
algorithm = "leiden",
chains = "CDR3b")The output in gcd is quite rich. See the help page for explanations.
Check whether you can interpret it!
names(gcd)The main element is the community_occupancy_matrix, which contains the
number of T cells in each community (rows) and repertoire (columns).
Here, we have two repertoires (two columns) and approximately 800
communities. This matrix serves as the main input to the function dco
(Step 3), which detects differences in community occupancy between
repertoires.
dim(gcd$community_occupancy_matrix)head(gcd$community_occupancy_matrix)Problem: Visualize the community abundance matrix with the function
get_honeycombs now!
honeycomb <- get_honeycombs(com = gcd$community_occupancy_matrix)
honeycomb[[1]]In the honeycomb plot, we see that most of the dots (communities) do not lie close to the diagonal, indicating that they drastically different number of T cells in CMV (y-axis) and Flu (x-axis). These communities contain more cells from either repertoire relative to the other. The color of the hexagons encodes the density of dots in the 2D scatterplot: dark hexagons represent regions with many communities, while light hexagons indicate sparsely populated areas. The plot uses log10-scaled axes.
Next step: Are there any, and if so which, differentially abundant communities (DACs) between CMV and Flu?
To detect DACs, ClustIRR employs a Bayesian model (remember other models discussed in this lecture series: Sleuth, DADA2, IgGeneUsage).
Using the community occupancy counts stored in the
community_occupancy_matrix as its main input, the
Dirichlet-Multinomial model applied by ClustIRR describes the relative
abundance (occupancy) of community \(i\) in repertoire \(j\) with the
parameter \(\beta_i^j\).
Problem: Run this analysis now!
dco_obj <- dco(community_occupancy_matrix = gcd$community_occupancy_matrix,
mcmc_control = list(mcmc_warmup = 500,
mcmc_iter = 1500,
mcmc_chains = 4,
mcmc_cores = 4))Before interpreting the model results, we need to ensure that the model is valid. A standard approach for this is to check whether the model can retrodict the observed data (i.e., the community occupancy matrix used to fit the model parameters).
General idea of posterior predictive checks:
ClustIRR provides both \(y\) and \(\hat{y}\) for each repertoire in
posterior_summary object within dco_obj. Visualize it with ggplot
now!
ggplot(data = dco_obj$posterior_summary$y_hat)+
geom_abline(intercept = 0, slope = 1, linetype = "dashed", col = "gray")+
geom_errorbar(aes(x = y_obs, y = mean, ymin = L95, ymax = H95),
col = "darkgray", width=0)+
geom_point(aes(x = y_obs, y = mean), size = 0.8)+
xlab(label = "observed y")+
ylab(label = "predicted y (and 95% HDI)")+
theme_bw()Problem: Based on the plot, does the model do a good job of “re-creating” the data it was trained on?
Larger communities are associated with larger, positive \(\beta\)
coefficients. Let’s visualize the mean \(\beta\) values for each community
across the different TCR repertoires using the get_beta_violin_ag
function.
This brings us to a key biological question: Are the communities with high \(\beta\) values specific to antigens from CMV or Flu?
bv_cmv <- get_beta_violin_ag(node_summary = gcd$node_summary,
beta = dco_obj$posterior_summary$beta,
ag_species = TRUE,
ag = "CMV",
db = "vdjdb")
bv_cmv$plotIn the plot, each point represents a community, shown with respect to repertoire (x-axis) and its mean \(\beta\) coefficient (y-axis).
Notice the orange dot with particularly high \(\beta\)s in repertoire CMV? This community contains over 10,000 cells specific for CMV antigens.
Let’s do the same analysis with respect to Flu antigens
bv_flu <- get_beta_violin_ag(node_summary = gcd$node_summary,
beta = dco_obj$posterior_summary$beta,
ag_species = TRUE,
ag = "InfluenzaA",
db = "vdjdb")
bv_flu$plotNotice that now we see the reverse: large orange dots with particularly high \(\beta\)s in repertoire Flu
Problem: what are the community IDs of these large orange dots?
bv_cmv$d[which(bv_cmv$d$sample == "CMV_A0201" &
bv_cmv$d$mean>6.5 &
bv_cmv$d$hit==T),]
# community 9# Large orange dots with high beta in repertoire Flu
bv_flu$d[which(bv_flu$d$sample == "Flu" &
bv_flu$d$mean>3 &
bv_flu$d$hit==T),]
# communities 34, 91, 94Community 9 has 33 clonotypes with paired CDR3\(\alpha\) and CDR3\(\beta\) sequences.
ns <- gcd$node_summary
nrow(ns[ns$community == 9,])We performed clustering based on the CDR3\(\beta\) sequences, hence, we can assume to find similar CDR3\(\beta\)s among the clonotypes. To verify this, we can construct a sequence logo from the CDR3 sequences
get_cdr3_motifs(node_summary = gcd$node_summary, community_id = 9)Based on the logo, community 9 is made of clonotypes with conserved CDR3\(\beta\)s (left logo) and diverse CDR3\(\alpha\)s (right logo).
sg <- subgraph(graph = gcd$graph, vids = V(gcd$graph)$community %in% 9)
V(sg)$Ag <- ifelse(test = regexec(pattern = "CMV", text = V(sg)$Ag_species)!=-1, yes = "+", no = "-")
V(sg)$col <- ifelse(test = V(sg)$Ag=="+", yes = "orange", no = "#A1BFD9")
V(sg)$label <- paste0("Community:", V(sg)$community, " (", V(sg)$clone_size, ")" ,
"\n CDR3a:",V(sg)$CDR3a, "\n CDR3b:",V(sg)$CDR3b)
V(sg)$label[which(ave(V(sg)$clone_size, V(sg)$community, FUN = function(x) x == max(x))==0)] <-
paste0("(", V(sg)$clone_size, ")")[which(ave(V(sg)$clone_size, V(sg)$community, FUN = function(x) x == max(x))==0)]
V(sg)$shape <- ifelse(test = V(sg)$sample=="Flu", yes = "circle", no = "square")
par(mai=c(0,0,0,0))
plot(sg,
vertex.label = V(sg)$label,
layout = igraph::layout_components,
vertex.label.cex = 0.8,
vertex.size = log(V(sg)$clone_size)+3,
vertex.frame.color = NA,
vertex.color = V(sg)$col,
vertex.shape = V(sg)$shape,
edge.width = 0.2)How to interpret this plot? Orange and blue nodes represent TCR clonotypes that are known to be CMV-specific or not, respectively. Most of the clonotypes in this community are known to recognize CMV antigens. Node widths represent the number of T cells within each clonotype (i.e., clonal expansion). Node labels show the community ID, CDR3\(\alpha\) and CDR3\(\beta\) sequences of the largest clonotype and clone sizes (in parentheses for all clonotypes). The nodes are densely connected owing to their similar CDR3\(\beta\) sequences.
d <- igraph::as_data_frame(x = sg, what = "vertices")
d[order(d$clone_size, decreasing = T), ][1:2, c("CDR3a", "CDR3b", "clone_size",
"sample", "Ag_species")]The two largest nodes in this graph have the same CDR3\(\beta\) sequence (CAWSISDLAKNIQYF), but different CDR3\(\alpha\)s. These clonotypes contain 9,754 and 6,052 T cells, respectively.
Problem: do the same analysis for Flu communities (e.g., 34, 91, 94)
# show the CDR3 logos
wrap_plots(lapply(X = c(34, 91, 94),
FUN = get_cdr3_motifs, node_summary = gcd$node_summary), ncol = 1)# Show the networks
sg <- subgraph(graph = gcd$graph, vids = V(gcd$graph)$community %in% c(34, 91, 94))
V(sg)$Ag <- ifelse(test = regexec(pattern = "CMV", text = V(sg)$Ag_species)!=-1, yes = "+", no = "-")
V(sg)$col <- ifelse(test = V(sg)$Ag=="+", yes = "orange", no = "#A1BFD9")
V(sg)$label <- paste0("Community:", V(sg)$community, " (", V(sg)$clone_size, ")" ,
"\n CDR3a:",V(sg)$CDR3a, "\n CDR3b:",V(sg)$CDR3b)
V(sg)$label[which(ave(V(sg)$clone_size, V(sg)$community, FUN = function(x) x == max(x))==0)] <-
paste0("(", V(sg)$clone_size, ")")[which(ave(V(sg)$clone_size, V(sg)$community, FUN = function(x) x == max(x))==0)]
V(sg)$shape <- ifelse(test = V(sg)$sample=="Flu", yes = "circle", no = "square")
par(mai=c(0,0,0,0))
plot(sg,
vertex.label = V(sg)$label,
layout = igraph::layout_components,
vertex.label.cex = 0.8,
vertex.size = log(V(sg)$clone_size)+3,
vertex.frame.color = NA,
vertex.color = V(sg)$col,
vertex.shape = V(sg)$shape,
edge.width = 0.2)Important remark: Notice that community 34 contains many clonotypes, where the largest clonotype contains nearly 8,000 cells. Yet, the largest clonotypes are not known to recognize Influenza antigens, i.e. we see low prevalence of antigen-specific CDR3s in the community. In contrast, community 91 and 94 are large (are associated with large \(\beta\) coefficients) and are also enriched with large clonotypes that are known to recognize Influenza antigens.
We now come to another important biological question: community 9 is expanded in TCR repertoire CMV and contains many clonotypes that are known to recognize CMV antigens. Do we also see higher levels of T cell activation among the corresponding T cells compared to the T cells outside this community?
Answer this question now!
For this we will follow the following procedure:
Extract T-cell activation profiles from gex for cells inside community 9
Compute for each cell the sum of expressions across 12 genes and use this as a signature of T-cell activation (higher expression means more T cell activation)
Do the same analysis for cells outside community 9
Compare the distributions of T-cell activation scores between cells of community 9 and the remaining cells
# this is the data.frame with all clonotypes
ns <- gcd$node_summary
# extract the clonotypes in community 9 from repertoire CMV_A0201
ns_com_9 <- ns[ns$sample == "CMV_A0201" & ns$community == 9,]
# extract the cell_ids (barcodes)
cell_ids <- unlist(strsplit(x = ns_com_9$cell_id, split = "\\,"))# extract the gex profiles of cells in community 9
gex_com_9 <- gex[gex$sample == "CMV_A0201" & gex$cell_id %in% cell_ids,]
# ... and in cells not in community 9
gex_com_not_9 <- gex[gex$sample == "CMV_A0201" & !gex$cell_id %in% cell_ids,]# compute the sum of gene expressions within each cell
t_act_9 <- apply(X = gex_com_9[,1:12], MARGIN = 1, FUN = sum)
t_act_not_9 <- apply(X = gex_com_not_9[,1:12], MARGIN = 1, FUN = sum)
# compare the distributions
# red = 9
# black = not 9
# do we see a clear difference? yes: red distribution shifted towards larger
# values -> up-regulation of genes involved in T-cell activation
plot(density(t_act_9), col = "red",
main = "T cell activation in community 9 (red) vs. rest (black)",
xlab = "T cell activation score")
points(density(t_act_not_9), type = "l", col = "black")Problem: do the same analysis for Flu communities 34, 91, and 94.
Can we draw any conclusion about the T cell activation scores from the specificity structure of these communities?
# this is the data.frame with all clonotypes
ns <- gcd$node_summary
# extract the clonotypes in community 34 from repertoire Flu
ns_com_34 <- ns[ns$sample == "Flu" & ns$community == 34,]
ns_com_91 <- ns[ns$sample == "Flu" & ns$community == 91,]
ns_com_94 <- ns[ns$sample == "Flu" & ns$community == 94,]
# extract the cell_ids (barcodes)
cell_ids_34 <- unlist(strsplit(x = ns_com_34$cell_id, split = "\\,"))
cell_ids_91 <- unlist(strsplit(x = ns_com_91$cell_id, split = "\\,"))
cell_ids_94 <- unlist(strsplit(x = ns_com_94$cell_id, split = "\\,"))# extract the gex profiles of cells in each community
gex_com_34 <- gex[gex$sample == "Flu" & gex$cell_id %in% cell_ids_34,]
gex_com_91 <- gex[gex$sample == "Flu" & gex$cell_id %in% cell_ids_91,]
gex_com_94 <- gex[gex$sample == "Flu" & gex$cell_id %in% cell_ids_94,]
# ... and in cells not in communities 34, 91, and 94
gex_com_not <- gex[gex$sample == "Flu" & !gex$cell_id %in% c(cell_ids_34,
cell_ids_91,
cell_ids_94),]# compute the sum of gene expressions within each cell
t_act_34 <- apply(X = gex_com_34[,1:12], MARGIN = 1, FUN = sum)
t_act_91 <- apply(X = gex_com_91[,1:12], MARGIN = 1, FUN = sum)
t_act_94 <- apply(X = gex_com_94[,1:12], MARGIN = 1, FUN = sum)
t_act_not <- apply(X = gex_com_not[,1:12], MARGIN = 1, FUN = sum)
# compare the distributions
plot(density(t_act_34), col = "pink",
main = "T cell activation in community 34, 91, 94, vs. rest",
xlab = "T cell activation score")
points(density(t_act_91), type = "l", col = "forestgreen")
points(density(t_act_94), type = "l", col = "green")
points(density(t_act_not), type = "l", col = "black")