Uncertainty-aware quantitative analysis of the structure and dynamics of T cell receptor repertoires

Reviewer #1 (Public review):

In this manuscript, the authors present ClustIRR, a computational tool that analyzes multiple TCR repertoires together, instead of one at a time. It builds one shared similarity graph across all the repertoires, then finds communities (CJs) on that graph that can be compared across samples. Building one shared graph across repertoires, rather than clustering each repertoire on its own, is a real improvement over existing tools. However, there are a few concerns that need to be addressed.

(1) The Introduction motivates ClustIRR by contrasting it with beta-binomial regression, Fisher's exact test, and scCODA/tascCODA (lines 62-77), but none of these are actually run on the same data in the Results. Could the authors include a direct comparison, e.g., applying a standard beta-binomial or Fisher's exact test to the Dataset 1 CJ occupancy matrix, to show the Bayesian model gives a lower false-positive rate or better-calibrated intervals than the alternatives it's positioned against?

(2) Prior predictive and posterior predictive checks are both reported as "(data not shown)" (lines 644, 728). Since the manuscript's central claim is rigorous, uncertainty-aware inference, it would help to include these diagnostic plots, along with Rhat and ESS values, in the supplement rather than stating they were checked.

(3) With 7,505 CJs tested simultaneously for differential occupancy in Dataset 1 alone (Figure 1 legend), what is the expected false discovery rate under the non-overlapping-HDI criterion used throughout? A short discussion of multiple-comparisons correction, or an argument for why it isn't needed under this framework, would strengthen the statistical claims.

(4) Line 115 states the Dataset 1 joint graph produced 10,301 CJs, of which 3,038 were singletons, leaving 7,263 non-singleton CJs. The Figure 1 legend reports 7,505 CJs used for the β modeling. Could the authors clarify how these two numbers relate - whether some singletons were included in the model, or a filtering step was applied that isn't described in Methods?

(5) The Methods section states that archival pretreatment tumor tissue was available for four patients (Pt4, Pt32, Pt36, Pt38; line 531), but the T+/T- DCJ analysis in Fig. 3B-C is only shown for Pt4. Was this analysis attempted in the other three patients? Extending it, even partially, would substantially strengthen the claim that contracting DCJs are enriched for tumor-infiltrating TCRs, which is currently based on a single patient.

(6) Dataset 1 was generated by deliberately stimulating T cells with EBV or MART1 antigen, so recovering EBV/MART1-annotated CJs from VDJdb is closer to a positive control than a blinded validation. Do the authors have, or could they obtain, any independent confirmation (e.g., tetramer data or an orthogonal cohort) for the CJs with large β that lack VDJdb annotation (orange dots, Figures 1B-C)?

Reviewer #2 (Public review):

Summary:

The study confronts a major obstacle in repertoire analysis. Given that individual TCR sequences are diverse and sparsely detected across repertoires, identifying sample-specific enrichment of TCRs based on sample-to-sample comparison of exact clonotype sequences can be intractable. This paper attempts to build on insights that sequence-similar TCRs can share antigen recognition, such that aggregating similar sequences derived from multi-sample joint-graph communities (i.e. clusters of tightly connected nodes) could reduce sparsity and boost signal.

Strengths:

The study is a well-motivated effort to address a need in the field. The paper takes a unique approach. The core method is modeling community occupancy with a hierarchical Dirichlet-Multinomial model that attempts to account for the high level of overdispersion present in repertoire sampling, a technique that has been previously applied to compositional microbiome data.

The authors apply this framework to both single-cell paired-chain and bulk single-chain TCR data. They develop a set of examples from public and synthetic data, with the most promising real-world application shown in reanalysis of longitudinal data during treatment of cancer patients with checkpoint inhibitors.

The manuscript is well structured and cogent. The authors are to be commended for contributing a well-documented, open-source R/Bioconductor package and for providing the underlying analysis datasets in a well-organized repository. In the joint graph construction step, the authors opt to use an existing implementation of the BLAST algorithm on CDR3 sequences, which is a slightly odd choice since it ignores potential contributions of other V-gene germline-encoded CDRs, but the authors also envision that their statistical package could be extended to include community graphs developed with other established TCR clustering tools. This will allow others to potentially explore the utility of Bayesian hierarchical Dirichlet-Multinomial models for differential occupancy analysis of immune receptors under varied clustering criteria.

The methods described here were demonstrated on relatively small datasets from 2-5 samples, and future work is likely needed to extend the joint graph differential occupancy concept to larger datasets. The authors are transparent about this and some of the other limitations in their current tool, most notably the computational cost of graph construction based on an all-versus-all sequence alignment to construct a joint sequence graph and the challenge of Bayesian parameter estimation as the number of subgraph entities scales with input data size. Since efficient approximate methods exist to find edges between similar text strings, the underlying idea of applying uncertainty-aware statistical inference to subgraph communities is promising, and the paper advances its primary goal.

Weaknesses:

The paper proposes the utility of the joint-graph community occupancy framework through three examples. I discuss potential weaknesses apparent in each separate example in turn.

(1) Weaknesses in Example 1

A broad weakness of the first results section ("Detecting EBV- and MART1-antigen reactive T cell communities from single cell datasets") is its reliance on a single vendor-generated dataset generated by the company ParseBio with no published experimental protocols and limited, if any, prior peer review. For reasons I will explore in greater detail below, the EBV-sample data may be particularly prone to chimeric pairings that confound the authors' primary analysis goals, and, at the very least, may not reflect physiologically realistic conditions for identifying antigen-reactive TCRs in other contexts.

Let us first consider Dataset 1 in more detail. The authors compare 2 antigen-stimulated repertoires with 3 unstimulated controls. To improve on single clonotype-level comparisons, the authors propose comparing the cell count aggregated across cells within joint graph communities constructed from paired CDR3 sequences across all the samples. Thus, one of the most relevant questions one hopes the authors answer in this section is whether the resolved graph communities are made up of many distinct clonotypes (i.e., are they polyclonal), allowing the method to function as intended by aggregating across multiple clones with putative shared antigen-reactivity.

Supplementary Figure 1B shows the size of all the communities with callouts for the putative EBV-expanded communities. The authors listed communities strongly enriched in the EBV-stimulated sample as e1, e2, e3, and e5. Each contains {greater than or equal to}100 clonotypes, and the authors note they contain at least one clonotype with a CDR3 sequence matching an EBV-annotated clone in VDJdb - a database of TCRs with some experimental evidence of epitope-reactivity. Community "e2" is notable for its remarkable size, including 1,339 unique clonotypes. At first glance, this seems promising for a method attempting to boost signal through community detection. However, it is worth re-investigating the individual clone sizes and sequences within this extraordinary community.

In the EBV-associated community "e2", a look at the data provided by the authors on the paper's GitHub repository suggests a single clonotype (clonotype_9; TRAV12-3 CATQGSNDYKLSF / TRBV9 CASSTGQVATNEKLFF) comprises 26,917 cells. As such, it makes up 29% of the sample with a total of 90,588 cells. If one clone supplies most of a graph community's cell counts, the community-level posterior estimate of β (Figure 1B) effectively tracks a single-clone estimate, and the premise of borrowing statistical power across a polyclonal expansion in this example is hard to assess.

There is also considerable evidence to believe that the apparent mega-polyclonality of cluster e2 may be partially an illusion, stemming from an artifact of this hyper-expanded clone's massive size and the experimental method used to assign TCRα-TCRβ chain pairings. In fact, the same α-chain CDR3 (CATQGSNDYKLSF) appears in ~1,219 clonotypes paired to distinct β chains, generating much of e2's remarkable 1,339-clonotype count. I believe two features warrant caution here. First, a single clonotype making up ~29% of total T cells in the sample is highly unexpected in ex vivo repertoires, suggesting intense non-physiological expansion conditions unlikely to generalize to other settings. That is, one would almost never expect to see a signal this strong.

Second, one TCR-α chain paired to ~1,220 distinct and diverse β chains in one sample is also biologically unexpected given what we know of the best-characterized epitope-specific responses for EBV, including to the well-known HLA-A*02 EBV BMLF-1 epitope, which recruits a tetramer-stained repertoire with conserved CDR3 motifs in both chains and at least some constraint on favored V-gene/α-β pairing (See Extended Data Figure 5 in Dash et al., Nature 2017). This raises the concerning possibility of barcode collision and mispairing against a hyperexpanded clone in the ParseBio split-pool method, unlikely to be robust to a clone occupying a third of the sample. Most of the ~1,219 β chains paired to the dominant α have a cell count of 1, further raising the concern of artifactual pairing versus genuine convergence. The second largest community "e1" also seems to suffer from the same issue, with a TCRβ sequence from one super clone making up 7% of the sample potentially being artifactually over-paired to >500 rare single-cell-count TCR α chains.

Taken together, these observations suggest that the authors' first positive control example passes but probably for the wrong reason since at least some of the "antigen-specific communities" are strongly anchored by a single hyper-clone. This could be remedied by repeating the same type of analysis on an ex vivo single-cell repertoire following more modest stimulation or natural infection (e.g., yellow-fever vaccine, influenza, or SARS-CoV-2 single-cell TCR datasets). I would advise future work using an alternative data source with better-documented experimental protocols, given the concerns above.

(2) Weaknesses in Example 2

Example 2 explores the application of Bayesian methods for identifying sample-enriched joint graph communities found in longitudinal data from many participants at two time points and longitudinal data from 1 person (Pt4) at 5 time points. A potential weakness of Case Study 2 is that it yields limited additional biological insight compared to what was previously shown by the authors of the underlying input data. Previously, Formenti et al. 2018 showed that the number of expanded clones after treatment strongly reflects responder status in this cohort, greater in patients with CR/PR versus SD or PD (See Figure 2b of the study). This 2018 primary analysis showed that tracking individual clonotypes was sufficient to reveal biological insight without the need for the computational demands of constructing a massive sequence similarity graph, finding communities on joint graphs, or Bayesian statistical inference. Thus, the impact of Case Study 2 in proving the unique utility of ClustIRR is somewhat diminished.

It is not clear how this study's result is more "robust" than the original analysis. Perhaps the authors could further clarify what is learned from the uncertainty-aware approach that could not be learned from exact clone tracking in time series.

Thus, example 2 shows that a complex method recapitulated the finding of a much simpler method for analyzing longitudinal TCR data where a strong signal of expansion was already present at the single-clonotype level. Since the Dirichlet method is sensitive to absolute counts, the large expanding clone in each community at 22 days may alone have carried most of the signal, which is not fully explored.

In this section, the authors make an interesting observation that CDR3β detected in both PBMC and patient-matched tumor samples were enriched in contracting communities (7/14) versus expanded communities (1/41). The authors do not indicate whether a similar enrichment of tumor-infiltrating lymphocytes (TILs) matched Day 0-22 contracting communities in the other 4 patients with tumor-matched samples, which, if consistent, would strengthen their finding.

(3) Weaknesses in Example 3

The final example explores "convergent repertoire differences between species." This is intriguing in principle but less informative due to the use of synthetic data with somewhat predictable properties that some may reasonably consider baked in by the data-generating process. That is, some of the results might be guaranteed by the way the data is constructed using the OLGA/IGoR generative model. For instance, the authors observe a positive correlation between CJ community size and Pgen of constituent clonotypes, stating: "This indicates that CDR3 sequences with high Pgen are statistically more likely to be generated, leading to convergence of similar sequences into public CJs." I may be mistaken, but this conclusion is almost guaranteed by the way the data-generating OLGA model outputs more similar high-Pgen sequences and fewer lower-Pgen sequences.

An interesting finding in this section is shown in Figure 4B, where community-level aggregation allowed for discrete clustering of human samples away from mouse samples that was not possible by comparing cosine similarity of a sparse clonotype occurrence matrix, a result that would be higher impact if it could be shown to separate real repertoire samples from humans with differential serology, vaccination status, or HLA backgrounds.

(4) Weaknesses in General

More generally, one aspect of the method that seems under-emphasized is the fact that many of the nodes in a multi-sample joint sequence similarity graph may have no edges. These zero-degree nodes would probably frequently occur in only one sample but be absent in other samples. It is not strongly emphasized in the paper how the model would infer whether such a singleton found in only one sample in the graph could be reliably inferred to be sample-specific enriched (see, for example, the large single node in Figure 3D, the orange node labeled "GQYF" in the far-right position of the lowest row in panel D). Presumably the number of cell counts represented in this single-sequence node is so great at sampled timepoint Day 22 as to yield a statistically strong signal in the multinomial model; however, the authors may wish to comment on how, for such singleton sequences, the power to detect sample-specific enrichment differs from prior single-clonotype-based methods.

With any large effort to find statistically significant features from a large candidate set examined all at once, a reader might be concerned with the potential for false discovery. Throughout, the authors seem to assign statistical significance when the 95% high-density interval (HDI) of the posterior estimate excludes zero or when the 95% HDIs of two features do not overlap. There is little discussion of how this implicitly handles multiplicity adjustment via shrinkage, which the authors could address more directly and explain more clearly to a broad audience, including many non-statisticians, who will read this paper.

Reviewer #3 (Public review):

Summary:

Analysis of immune receptor repertoires (IRR) needs to take into account the underlying diversity of the repertoires analysed, and the limitations inherent to the technologies used to measure IRRs: limited sampling depth relative to total number of cells and clonotypes, and experimental noise. In this work, Kitanovski and colleagues present ClustIRR. ClustIRR proposes to improve the analysis of immune receptor repertoires, specifically TCRs in the presented applications, by performing two steps: (1) consistent and comparable sequence clustering across repertoires, to account for sparsity of sampling, and (2) estimation of sequence clusters of interest using a Bayesian approach. They showcase ClustIRR performance in 3 scenarios: detection of antigen-specific T cells in a peptide-stimulation, detection of T cells responding to immunotherapy in the context of lung cancer, and analysis of mouse and human T cell repertoires.

Strengths:

(1) The cluster occupancy calculation presents an important conceptual framework which would be of interest and useful to the TCR repertoire field as it smoothly integrates information over a set of experimental conditions or time points. The application to longitudinal TCR sequencing datasets is particularly interesting, and could easily be extended to BCR sequencing datasets. Moreover, the calculation can be performed with any user-defined grouping of TCR clones, which allows for usage of other existing methods as the user wishes.

(2) The results of the human and mouse repertoires provide a very insightful argument for the use of metaclones as opposed to single clones for analysis of repertoires compared to single clones.

Weaknesses:

(1) While ClustIRR provides an interesting framework to analyse TCR sequencing datasets, it is not clear whether ClustIRR can perform more informative sequence clustering than state-of-the-art methods. A comparison of obtained clusters with existing methods would provide a useful benchmark. Moreover, computation time scales quite fast with the number of sequences included. This is a major limitation, as the authors state that a time of 2.5 hours is required for clustering of ~100,000 sequences, a number of clones that can easily be reached when analysing multiple samples together.

(2) The authors claim in the abstract that ClustIRR is integrated with gene expression data. However, in the results presented, the gene expression and TCR sequencing data are analysed separately, and the results are simply correlated. No real integration in the analysis exists for these two data types. The claim should be removed from the abstract. Moreover, the differential gene expression section, while it presents interesting results, lacks clarity and transparency. Presentation of the data in more transparent ways (such as showing violin plots or clustering on the UMAP) would increase the strength of the claims.

(3) The score calculation does not seem to be normalized to take into account the underlying diversities of CDR3a and CDR3b. While still a useful metric for sequence clustering, I worry about the impact of the lack of normalisation on the conclusion that the "alpha chain drives functional convergence through germline bias". The observation that clustering is mostly driven by the J gene is known and expected (https://pmc.ncbi.nlm.nih.gov/articles/PMC5553937/), as the J gene has lower diversity and greater overlap with the definition of CDR3. Because CDR3a is lower diversity generally from CDR3b, it will likely dominate the sequence similarity signal. Thus, it will appear that Ja drives the signal. The authors do try to address this by looking at clusters driven by CDR3b similarity. However, a low number of CDR3b-driven clusters is consistent with the similarity definition. I wonder if instead the appropriate control for this analysis would be to run the same analysis disregarding CDR3a sequence altogether, and quantify whether similar species-specific DCJ are identified when only CDR3b similarity is used? Absence or reduction of species-specific clustering would confirm that the effect is driven by the CDR3a sequence.