Run cluster_irr

Step 1. involves calling the function clust_irr which returns an S4 object of class clust_irr.

# perform clust_irr analysis of repertoire a, b and c
cl_a <- cluster_irr(s = a, meta = meta_a, control = list(gmi = 0.7))

cl_b <- cluster_irr(s = b, meta = meta_b, control = list(gmi = 0.7))

cl_c <- cluster_irr(s = c, meta = meta_c, control = list(gmi = 0.7))

Next, we show the chain-specific similarity scores between CDR3s sequences. Each row is a pair of CDR3 sequences from the repertoire. For each pair we have the following metrics:

• max_len: length of the longer CDR3 sequence in the pair
• max_clen: length of the longer CDR3 core sequence in the pair
• weight: ω = BLOSUM62 score of the complete CDR3 alignment
• cweight= ω_c: BLOSUM62 score of CDR3 cores
• nweight = ω̄: normalized weight by max_len
• ncweight = ω̄_c: normalized cweight by max_clen

The results for CDR3α sequence pairs from repertoire a:

kable(head(cl_a@clust$CDR3a), digits = 2)

… the same table for CDR3β sequence pairs from repertoire a:

kable(head(cl_a@clust$CDR3b), digits = 2)

Notice that very similar CDR3α or CDR3β pairs have high normalized alignment scores (ω̄). We have similar tables for repertoire b and c.

Annotation of CDR3s

The function clust_irr performs automatic annotation of TCR clones based on multiple databases (DBs), including: VDJdb, TCR3d, McPAS-TCR. Each TCR clone recieves two types of annotations (per chain and per database):

• Ag_species: the name of the antigen species recognized by the clone’s CDR3

• Ag_gene: the name of the antigen gene recognized by clone’s CDR3

Annotation process: annotation is performed based on an edit distance threshold, controlled by the user through the parameter control$db_edit (default = 0). If the edit distance between an input CDR3 sequence and a database (DB) CDR3 sequence is less than or equal to control$db_edit, the input CDR3 inherits the annotation data of the matched DB CDR3.

This information will be used for various downstream analyses (see the next steps).

Inspect graphs with plot_graph

We can use the function plot_graph for interactive visualization of relatively small graphs.

For instance, we can highligh TCR clones that recognize certain antigenic species or genes (see dropdown menu), or use the hoovering function to look at the CDR3 sequences of nodes.

Do this now!

plot_graph(g, select_by = "Ag_species", as_visnet = TRUE)

get_graph and get_joint_graph generate igraph objects

We can use igraph functions to inspect various properties of the graph. For instance, below, we extract the edge attributes and visualize the distributions of the edge attributes ncweight and nweight for all CDR3α and CDR3β sequence pairs.

# data.frame of edges and their attributes
e <- igraph::as_data_frame(x = g$graph, what = "edges")

ggplot(data = e)+
  geom_density(aes(ncweight, col = chain))+
  geom_density(aes(nweight, col = chain), linetype = "dashed")+
  xlab(label = "edge weight (solid = ncweight, dashed = nweight)")+
  theme(legend.position = "top")

Can you guess why we observe bimodal distributions? (Hint: CDR3s have different lengths)

Inspecting the outputs of detect_communities

The function detect_communities generates a complex output. Lets investigate its elements:

names(gcd)

[1] "community_occupancy_matrix" "community_summary"         
[3] "node_summary"               "graph"                     
[5] "graph_structure_quality"    "input_config"              

The main element is community_occupancy_matrix, which contains the number of T-cells in each community (row) and repertoire (column). Here we have three repertoires (three columns) and about 260 communities. This matrix is the main input of the function dco (step 5.), to detect differences in the community occupancy between repertoires.

dim(gcd$community_occupancy_matrix)

[1] 261   3

head(gcd$community_occupancy_matrix)

   a  b  c
1 10 10 10
2  2  2  2
3  2  2  2
4  5  7  9
5 28 31 34
6  5  5  5

Visualizing community abundance matrices with get_honeycombs

honeycomb <- get_honeycombs(com = gcd$community_occupancy_matrix)

In the honeycomb plots shown below, several communities (black dots) appear far from the diagonal. This indicates that these communities contain more cells in repertoires b and c (y-axes in panels A and B) compared to repertoire a (x-axes in panels A and B). Meanwhile, the same points are generally close to the diagonal in panel C but remain slightly more abundant in repertoire c (y-axis) compared to b (x-axis).

In step 5, ClustIRR will provide a quantitative answer to the question: Which communities are differentially abundant between pairs of repertoires?

Importantly, the color of the hexagons encodes the density of communities in the 2D scatterplots: dark hexagons indicate high a frequency of communities, while light hexagons represent sparsely populated regions.

wrap_plots(honeycomb, nrow = 1)+
    plot_annotation(tag_levels = 'A')

Also see community_summary. In the data.frame wide we provide community summaries in each row across all samples, including:

• clones_a, clone_b, clone_c, clones_n: the frequency of clones in the community coming from repertoire a, b, c and in total (n)
• cells_a, cells_b, cells_c, cells_n: the frequency of cell in the community coming from repertoires a, b, c and in total (n)
• w: average inter-clone similarity
• ncweight_CDR3a/b, nweight_CDR3a/b: raw and normalized similarity for CDR3α and CDR3β sequences
• n_CDR3a, n_CDR3b: number of edges between CDR3α and CDR3β sequences

kable(head(gcd$community_summary$wide), digits = 1)

In the data.frame tall we provide community and sample/repertoire summaries in each row.

kable(head(gcd$community_summary$tall), digits = 2)

Node-specific (TCR clone-specific) summaries are provided in node_summary

kable(head(gcd$node_summary), digits = 2)

Decoding communities with decode_community

Community with ID=8 contains 42 TCR clones. These are connected based on their CDR3β sequences.

ns_com_8 <- gcd$node_summary[gcd$node_summary$community == 8, ]

table(ns_com_8$sample)

 a  b  c 
14 14 14 

We can extract it and visualize it using igraph:

par(mai = c(0,0,0,0))
plot(subgraph(graph = gcd$graph, vids = which(V(gcd$graph)$community == 8)))

Furthermore, we can “decompose” this graph using decode_community based on the attributes of the edges (edge_filter) and nodes (node_filter).

# we can pick from these edge attributes
edge_attr_names(graph = gcd$graph)

[1] "nweight"  "ncweight" "chain"    "w"       

The following edge-filter will instruct ClustIRR to keep edges with with edge attributes: nweight>=3 AND ncweight>=3

edge_filter <- rbind(data.frame(name = "nweight", value = 3, operation = ">="),
                     data.frame(name = "ncweight", value = 3, operation = ">="))

# we can pick from these node attributes
vertex_attr_names(graph = gcd$graph)

 [1] "name"             "clone_id"         "Ag_gene"          "Ag_species"      
 [5] "info_tcr3d_CDR3a" "db_tcr3d_CDR3a"   "info_tcr3d_CDR3b" "db_tcr3d_CDR3b"  
 [9] "info_mcpas_CDR3a" "db_mcpas_CDR3a"   "info_mcpas_CDR3b" "db_mcpas_CDR3b"  
[13] "info_vdjdb_CDR3a" "db_vdjdb_CDR3a"   "info_vdjdb_CDR3b" "db_vdjdb_CDR3b"  
[17] "clone_size"       "sample"           "CDR3b"            "CDR3a"           
[21] "cell"             "HLA_A"            "age"              "TRAV"            
[25] "TRAJ"             "TRBV"             "TRBJ"             "community"       

The following node-filter will instruct ClustIRR to retain edges between nodes that have shared node attributed with respect to ALL of the following node attributes:

node_filter <- data.frame(name = c("TRBV", "TRBJ", "HLA_A"))

c8 <- decode_communities(community_id = 8, 
                         graph = gcd$graph,
                         edge_filter = edge_filter,
                         node_filter = node_filter)

par(mar = c(0, 0, 0, 0))
plot(c8, vertex.size = 10)

as_data_frame(x = c8, what = "vertices")[,c("name", "component_id", 
                                            "CDR3b", "TRBV", "TRBJ",
                                            "CDR3a", "TRAV", "TRAJ")]

       name component_id           CDR3b     TRBV    TRBJ              CDR3a
a|8     a|8            1  CASSGRGGDNEQFF  TRBV7-9 TRBJ2-1      CASGTLSYNEQFF
a|271 a|271            1  CASSQAGGRNEQFF  TRBV7-9 TRBJ2-1   CTSSQVGVGMNTEAFF
a|22   a|22            1 CASSQAGYSYNEQFF    TRBV9 TRBJ2-1     CASSPPLAEETQYF
a|23   a|23            1  CASSYSGHYNEQFF  TRBV6-5 TRBJ2-1     CSVVGQLTGYEQYF
a|124 a|124            1  CASAPGTSYNEQFF  TRBV7-6 TRBJ2-1   CASSLYGTERADEQYF
a|166 a|166            1  CASSLGTPYNEQFF   TRBV13 TRBJ2-1      CASGQTDPNGYTF
a|340 a|340            1  CASSGRSSYNEQFF TRBV12-4 TRBJ2-1   CASSLPAERPDYGYTF
a|379 a|379            1  CASSQGRSYNEQFF  TRBV4-1 TRBJ2-1      CSASGTEDNEQFF
a|408 a|408            1  CASSPLTGRNEQFF   TRBV27 TRBJ2-1     CASSLGASEYEQYF
a|429 a|429            1  CASSGAGQFNEQFF   TRBV19 TRBJ2-1       CSVRDLGETQYF
a|456 a|456            1  CASSELQSYNEQFF  TRBV6-1 TRBJ2-1 CSGAGWIGFPDQYNEQFF
a|464 a|464            1  CASSLAGGYNEQFF  TRBV6-8 TRBJ2-1    CASSIYMGKNTEAFF
a|482 a|482            1  CASSLLGSYNEQFF  TRBV6-2 TRBJ2-1  CASMGRQGRDGNEKLFF
a|500 a|500            1  CASSSSPGDNEQFF  TRBV7-3 TRBJ2-1     CASSYLAAQETQYF
b|8     b|8            1  CASSGRGGDNEQFF  TRBV7-9 TRBJ2-1      CASGTLSYNEQFF
b|271 b|271            1  CASSQAGGRNEQFF  TRBV7-9 TRBJ2-1   CTSSQVGVGMNTEAFF
b|22   b|22            1 CASSQAGYSYNEQFF    TRBV9 TRBJ2-1     CASSPPLAEETQYF
b|23   b|23            1  CASSYSGHYNEQFF  TRBV6-5 TRBJ2-1     CSVVGQLTGYEQYF
b|124 b|124            1  CASAPGTSYNEQFF  TRBV7-6 TRBJ2-1   CASSLYGTERADEQYF
b|166 b|166            1  CASSLGTPYNEQFF   TRBV13 TRBJ2-1      CASGQTDPNGYTF
b|340 b|340            1  CASSGRSSYNEQFF TRBV12-4 TRBJ2-1   CASSLPAERPDYGYTF
b|379 b|379            1  CASSQGRSYNEQFF  TRBV4-1 TRBJ2-1      CSASGTEDNEQFF
b|408 b|408            1  CASSPLTGRNEQFF   TRBV27 TRBJ2-1     CASSLGASEYEQYF
b|429 b|429            1  CASSGAGQFNEQFF   TRBV19 TRBJ2-1       CSVRDLGETQYF
b|456 b|456            1  CASSELQSYNEQFF  TRBV6-1 TRBJ2-1 CSGAGWIGFPDQYNEQFF
b|464 b|464            1  CASSLAGGYNEQFF  TRBV6-8 TRBJ2-1    CASSIYMGKNTEAFF
b|482 b|482            1  CASSLLGSYNEQFF  TRBV6-2 TRBJ2-1  CASMGRQGRDGNEKLFF
b|500 b|500            1  CASSSSPGDNEQFF  TRBV7-3 TRBJ2-1     CASSYLAAQETQYF
c|8     c|8            1  CASSGRGGDNEQFF  TRBV7-9 TRBJ2-1      CASGTLSYNEQFF
c|271 c|271            1  CASSQAGGRNEQFF  TRBV7-9 TRBJ2-1   CTSSQVGVGMNTEAFF
c|22   c|22            1 CASSQAGYSYNEQFF    TRBV9 TRBJ2-1     CASSPPLAEETQYF
c|23   c|23            1  CASSYSGHYNEQFF  TRBV6-5 TRBJ2-1     CSVVGQLTGYEQYF
c|124 c|124            1  CASAPGTSYNEQFF  TRBV7-6 TRBJ2-1   CASSLYGTERADEQYF
c|166 c|166            1  CASSLGTPYNEQFF   TRBV13 TRBJ2-1      CASGQTDPNGYTF
c|340 c|340            1  CASSGRSSYNEQFF TRBV12-4 TRBJ2-1   CASSLPAERPDYGYTF
c|379 c|379            1  CASSQGRSYNEQFF  TRBV4-1 TRBJ2-1      CSASGTEDNEQFF
c|408 c|408            1  CASSPLTGRNEQFF   TRBV27 TRBJ2-1     CASSLGASEYEQYF
c|429 c|429            1  CASSGAGQFNEQFF   TRBV19 TRBJ2-1       CSVRDLGETQYF
c|456 c|456            1  CASSELQSYNEQFF  TRBV6-1 TRBJ2-1 CSGAGWIGFPDQYNEQFF
c|464 c|464            1  CASSLAGGYNEQFF  TRBV6-8 TRBJ2-1    CASSIYMGKNTEAFF
c|482 c|482            1  CASSLLGSYNEQFF  TRBV6-2 TRBJ2-1  CASMGRQGRDGNEKLFF
c|500 c|500            1  CASSSSPGDNEQFF  TRBV7-3 TRBJ2-1     CASSYLAAQETQYF
          TRAV    TRAJ
a|8     TRAV19 TRAJ2-1
a|271  TRAV4-2 TRAJ1-1
a|22   TRAV7-9 TRAJ2-5
a|23  TRAV29-1 TRAJ2-7
a|124 TRAV12-4 TRAJ2-7
a|166 TRAV25-1 TRAJ1-2
a|340  TRAV7-2 TRAJ1-2
a|379 TRAV20-1 TRAJ2-1
a|408 TRAV12-4 TRAJ2-7
a|429 TRAV29-1 TRAJ2-5
a|456 TRAV29-1 TRAJ2-1
a|464   TRAV19 TRAJ1-1
a|482   TRAV19 TRAJ1-4
a|500  TRAV6-6 TRAJ2-5
b|8     TRAV19 TRAJ2-1
b|271  TRAV4-2 TRAJ1-1
b|22   TRAV7-9 TRAJ2-5
b|23  TRAV29-1 TRAJ2-7
b|124 TRAV12-4 TRAJ2-7
b|166 TRAV25-1 TRAJ1-2
b|340  TRAV7-2 TRAJ1-2
b|379 TRAV20-1 TRAJ2-1
b|408 TRAV12-4 TRAJ2-7
b|429 TRAV29-1 TRAJ2-5
b|456 TRAV29-1 TRAJ2-1
b|464   TRAV19 TRAJ1-1
b|482   TRAV19 TRAJ1-4
b|500  TRAV6-6 TRAJ2-5
c|8     TRAV19 TRAJ2-1
c|271  TRAV4-2 TRAJ1-1
c|22   TRAV7-9 TRAJ2-5
c|23  TRAV29-1 TRAJ2-7
c|124 TRAV12-4 TRAJ2-7
c|166 TRAV25-1 TRAJ1-2
c|340  TRAV7-2 TRAJ1-2
c|379 TRAV20-1 TRAJ2-1
c|408 TRAV12-4 TRAJ2-7
c|429 TRAV29-1 TRAJ2-5
c|456 TRAV29-1 TRAJ2-1
c|464   TRAV19 TRAJ1-1
c|482   TRAV19 TRAJ1-4
c|500  TRAV6-6 TRAJ2-5

Visualizing the distribution of β with get_beta_violins

beta_violins <- get_beta_violins(node_summary = gcd$node_summary,
                                 beta = d$posterior_summary$beta,
                                 ag_species = NULL,
                                 db = "vdjdb",
                                 db_dist = 0,
                                 chain = "both")

beta_violins$violins

[[1]]

Compare βs of clones specific for CMV, EBV, flu or MLANA?

beta_violins <- get_beta_violins(node_summary = gcd$node_summary,
                                 beta = d$posterior_summary$beta,
                                 ag_species = c("CMV", "EBV", "Influenza"),
                                 ag_genes = c("MLANA"),
                                 db = "vdjdb",
                                 db_dist = 0,
                                 chain = "both")

patchwork::wrap_plots(beta_violins$violins, ncol = 2)

Compare β between repertoires with get_beta_scatterplot

beta_scatterplots <- get_beta_scatterplot(node_summary = gcd$node_summary,
                                          beta = d$posterior_summary$beta,
                                          ag_species = c("CMV"),
                                          db = "vdjdb",
                                          db_dist = 0,
                                          chain = "both")

patchwork::wrap_plots(beta_scatterplots$scatterplots[[1]], ncol = 3)

Posterior predictive checks

Before we can start interpreting the model results, we have to make sure that the model is valid. One standard approach is to check whether our model can retrodict the observed data (community occupancy matrix) which was used to fit model parameters.

General idea of posterior predictive checks:

1. fit model based on data y
2. simulate new data ŷ
3. compare y and ŷ

ClustIRR provides y and ŷ of each repertoire, which we can visualize with ggplot:

ggplot(data = d$posterior_summary$y_hat)+
  facet_wrap(facets = ~sample, nrow = 1, scales = "free")+
  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)")

Differential community abundance results → par. δ

Given two repertoires, a and b, we can quantify the differential community occupancy (DCO):

Importantly, ClustIRR always computes both contrasts (a vs. b and b vs. a): δ_i^a − b and δ_i^b − a.

ggplot(data = d$posterior_summary$delta)+
    facet_wrap(facets = ~contrast, ncol = 2)+
    geom_errorbar(aes(x = community, y = mean, ymin = L95, ymax = H95), 
                  col = "lightgray", width = 0)+
    geom_point(aes(x = community, y = mean), shape = 21, fill = "black", 
               stroke = 0.4, col = "white", size = 1.25)+
    theme(legend.position = "top")+
    ylab(label = expression(delta))+
    scale_x_continuous(expand = c(0,0))

Given two repertoires, a and b, ClustIRR also quantifies absolute differences in community probabilities:

Again, both contrasts are computed: ϵ_i^a − b and ϵ_i^b − a.

ggplot(data = d$posterior_summary$epsilon)+
    facet_wrap(facets = ~contrast, ncol = 2)+
    geom_errorbar(aes(x = community, y = mean, ymin = L95, ymax = H95), 
                  col = "lightgray", width = 0)+
    geom_point(aes(x = community, y = mean), shape = 21, fill = "black", 
               stroke = 0.4, col = "white", size = 1.25)+
    theme(legend.position = "top")+
    ylab(label = expression(epsilon))+
    scale_x_continuous(expand = c(0,0))