A common strategy for partitioning a BCR repertoire dataset into clonal families is to go through all pairs of sequences and identify pairs of clonally related sequences. In the following, we call such related pairs positive, and pairs of sequences belonging to different families negative. Then, the partition is built by single-linkage clustering, which consists of recursively grouping all positive pairs. Two characteristics of the repertoire complicate the search for this partition: large total number of pairs and low proportion of positive pairs. In this section we analyze and model the statistics of pairs of sequences in natural repertoires to inform our choice of the clustering method and parameters. In the next section we will leverage that analysis to design an optimized clustering procedure. To help following notations, a summary of their definitions is provided in Table 1.

A pair of related sequences is expected to share the same V and J genes, as well as the same CDR3 length ${\displaystyle l}$, as determined by alignment to the templates (Figure 1A). The methods developed here begin by partitioning the data into ${\displaystyle \mathrm{V}\mathrm{J}l}$ classes, defined as subsets of sequences with the same V and J gene usage, and CDR3 length ${\displaystyle l}$ (Figure 1B). For a description of the data preprocessing and alignment to the V and J gene templates, see Methods. Clustering will then be performed within each ${\displaystyle \mathrm{V}\mathrm{J}l}$ class independently. While this first step severely limits the number of unnecessary comparisons, some ${\displaystyle \mathrm{V}\mathrm{J}l}$ classes still exceed 10⁵ sequences in large datasets, leading to the order of 10¹⁰ pairs (see Figure 1C for the distribution of the ${\displaystyle \mathrm{V}\mathrm{J}l}$ class sizes ${\displaystyle N}$ for donor 326651 of Briney et al., 2019).

The CDR3 plays an important role in encoding the signature of the VDJ rearrangement. As we will see, the CDR3 length ${\displaystyle l}$ has a strong impact on the difficulty of clonal family reconstruction. The distribution of CDR3 lengths ${\displaystyle l}$ observed in the data is shown in Figure 1D. In what follows we restrict our analysis to sequences with CDR3 lengths a multiple of 3 and between 15 and 105, relying on the common approximation that sequences with no frameshift in the CDR3 come from a productive naive ancestor. The number of sequences with length larger than 105 is too small to reach meaningful conclusions, and sequences of length smaller than 15 are likely nonfunctional (as evidenced by the similar number of in-frame and out-of-frame sequences in Figure 1D).

In each ${\displaystyle \mathrm{V}\mathrm{J}l}$ class, we call prevalence and denote by ${\displaystyle \rho }$ the proportion of positive pairs, i.e., the number of positive pairs divided by the total number of pairs. This quantity is unknown in the absence of the ground-truth partition. However, we can estimate it from the statistics of pairwise distances. We compute the Hamming distance ${\displaystyle n}$ of each pair of CDR3s, defined as the number of positions at which the two nucleotide sequences differ. The distribution of these distances normalized by the CDR3 length, denoted by ${\displaystyle x}$, shows a clear bimodal structure in data (donor 326651 of Briney et al., 2019), with two identifiable components (Figure 2A): the contribution of positive pairs (of proportion ${\displaystyle \rho }$) peaks near ${\displaystyle x=0}$ and decays quickly, whereas the bell-shaped contribution of negative pairs (of proportion ${\displaystyle 1-\rho }$) peaks around ${\displaystyle x=1/2}$.

Figure 2 with 7 supplements see all

Download asset Open asset

The prevalence ${\displaystyle \rho }$ can be formally written as ${\displaystyle {\displaystyle [{\mathrm{\Sigma }}_i{z}_i(z_i-1)/2]/[N(N-1)/2]}}$, where ${\displaystyle z_i}$ denote the sizes of the clonal families in the ${\displaystyle \mathrm{V}\mathrm{J}l}$ class, but we do not know these sizes before the partition into families is found. To overcome this issue, we developed a method to estimate ${\displaystyle \rho }$ a priori, without knowing the family structure (Methods). We do this by fitting the empirical distribution of ${\displaystyle x}$ as a mixture model, ${\displaystyle P(x)=\hat{\rho }{P}_{\mathrm{T}}(x)+(1-\hat{\rho })P_{\mathrm{F}}(x)}$, where ${\displaystyle P_{\mathrm{T}}(x)}$ and ${\displaystyle P_{\mathrm{F}}(x)}$ are the distributions of distances between positive (T as true) and negative (F as false) pairs (Figure 2C and D), estimated as follows. ${\displaystyle P_{\mathrm{F}}(x)=P_{\mathrm{F}}(x|l)}$ is computed for each length ${\displaystyle l}$ by generating a large number of unrelated, same-length sequences with the soNNia model of recombination and selection (Isacchini et al., 2021), and calculating the distribution of their pairwise distances (Methods). ${\displaystyle P_{\mathrm{T}}(x)}$ is approximated by a Poisson distribution, ${\displaystyle P_{\mathrm{T}}(x)=(\mu l{)}^{xl}{e}^{-\mu l}/(xl)!}$, with adjustable parameter ${\displaystyle \mu }$, which is proportional to the average hypermutation rate within the clone. The fit of ${\displaystyle P(x)}$ by the mixture model is performed for each ${\displaystyle \mathrm{V}\mathrm{J}l}$ class with an expectation-maximization algorithm which finds maximum likelihood estimates of the prevalence ${\displaystyle \hat{\rho }}$ and mean intra-family distance ${\displaystyle \hat{\mu }}$, the only free parameters of the mixture model.

The results of the fit to real data (donor 326651 of Briney et al., 2019) show that ${\displaystyle \hat{\mu }}$ varies little between ${\displaystyle \mathrm{V}\mathrm{J}l}$ classes, around ${\displaystyle \hat{\mu }\simeq 4\mathrm{\%}}$ (Figure 2—figure supplement 1). In contrast, the prevalence ${\displaystyle \hat{\rho }}$ varies widely across classes, spanning three orders of magnitude (Figure 2B). In addition, when we examine the ${\displaystyle \mathrm{V}\mathrm{J}l}$ classes with increasing CDR3 length ${\displaystyle l}$, we find that the part of the model distribution corresponding to positive pairs, ${\displaystyle P_{\mathrm{T}}(x)}$, varies little, whereas the model distribution over negative pairs ${\displaystyle P_{\mathrm{F}}(x)}$ becomes more and more peaked around 1/2 (Figure 2 and Figure 2—figure supplement 2), making the two categories more easily separable.

We want to build a classifier between positive and negative pairs using the normalized distance ${\displaystyle x}$ alone, by setting a threshold ${\displaystyle t}$ so that pairs are called positive if ${\displaystyle x\le t}$, and negative otherwise. Using our model for ${\displaystyle P(x)}$, for any given ${\displaystyle t}$ we can evaluate the number of true positives (${\displaystyle {\displaystyle \widehat{\mathrm{T}\mathrm{P}}}}$) and false negatives (${\displaystyle {\displaystyle \widehat{\mathrm{F}\mathrm{N}}}}$) among all positive pairs (${\displaystyle {\displaystyle \hat{\text{P}}=\widehat{\mathrm{T}\mathrm{P}}+\widehat{\mathrm{F}\mathrm{N}}}}$), as well as true negatives (${\displaystyle {\displaystyle \widehat{\mathrm{T}\mathrm{N}}}}$) and false positives (${\displaystyle {\displaystyle \widehat{\mathrm{F}\mathrm{P}}}}$) among the negative pairs (${\displaystyle {\displaystyle \hat{\text{N}}=\widehat{\mathrm{T}\mathrm{N}}+\widehat{\mathrm{F}\mathrm{P}}}}$), as schematized in Figure 2—figure supplement 2C and D.

Our goal is to set a threshold ${\displaystyle t}$ that ensures a high precision, ${\displaystyle \hat{\pi }(t)}$, defined as a proportion of true positives among all pairs classified as positive (Figure 2E). In a single-linkage clustering approach, we will join two clusters with at least one pair of positive sequences between them. Therefore, it is critical to limit the number of false positives, which can cause the erroneous merger of large clusters. We can write:

and ${\displaystyle \hat{s}(t)}$ is the estimated sensitivity (Figure 2D), evaluated from the Poisson fit to ${\displaystyle P_{\mathrm{T}}}$ (Methods):

Finally, the estimated prevalence ${\displaystyle {\displaystyle \hat{\rho }\equiv \hat{\mathrm{P}}/(\hat{\mathrm{P}}+\hat{\mathrm{N}})}}$ is inferred from the ${\displaystyle P(x)}$ distribution as explained above.

Figure 2H shows ${\displaystyle \hat{\pi }(t)}$ as a function of ${\displaystyle t}$ for different CDR3 lengths and a fixed value of ${\displaystyle \hat{\rho }}$. For each ${\displaystyle \mathrm{V}\mathrm{J}l}$ class, we define the threshold ${\displaystyle t=t^{\ast }}$ that reaches 99% precision, ${\displaystyle \hat{\pi }(t^{\ast })={\pi }^{\ast }=99\mathrm{\%}}$, by inverting Equation 1. This adaptive threshold depends on the ${\displaystyle VJl}$ class through the CDR3 length ${\displaystyle l}$ and the prevalence ${\displaystyle \rho }$, and it increases with both (Figure 2I): low clonality (small ${\displaystyle \rho }$) means few positive pairs and a smaller adaptive threshold, while short CDR3 means less information and a stricter inclusion criterion.

The predicted sensitivity, ${\displaystyle \hat{s}(t^{\ast })}$, which tells us how much of the positives we are capturing, is shown in Figure 2J. We conclude that for a wide range of parameters, the method is predicted to achieve both high precision and high sensitivity. However, it is expected to fail when the prevalence and the CDR3 length are both low. At the extreme, for small values ${\displaystyle \rho }$ and ${\displaystyle l}$, even joining together identical CDR3s (${\displaystyle t=0}$) results in poor precision because of convergent recombination (reflected by ${\displaystyle t^{\ast }<0}$).

The resulting procedure, which we call HILARy-CDR3, can be applied to Ig repertoire data through the following steps: (1) group sequences by ${\displaystyle \mathrm{V}\mathrm{J}l}$ class; (2) in each class, fit the mixture model to the distribution of pairwise distance to infer ${\displaystyle \hat{\rho }}$ and ${\displaystyle \hat{\mu }}$; (3) invert Equations 1–3 to find the high-precision threshold ${\displaystyle t^{\ast }}$; (4) classify positive and negative pairs according to that threshold; (5) complete the partition by applying single-linkage clustering to positive pairs.

So far we have presented a method to set a high-precision threshold with predictable sensitivity, based on estimates from the distribution of distances ${\displaystyle P(x)}$ only. To verify that these performance predictions hold in a realistic inference task, we designed a method to generate realistic synthetic datasets where the clonal family structure is known. This generative method will also be used in the next sections to create a benchmark for comparing different clustering algorithms.

We first estimated the distribution of clonal family sizes from the data of Briney et al., 2019, by applying HILARy-CDR3 with adaptive threshold as described above to ${\displaystyle \mathrm{V}\mathrm{J}l}$ classes for which the inference was highly reliable, i.e. for which the predicted sensitivity was ${\displaystyle z\in [10,100]}$. In that limit, clusters are clearly separated and the partition should depend only weakly on the choice of clustering method. The resulting distribution of clone sizes follows a power-law with exponent –2.3.

To create a synthetic lineage, we first draw a random progenitor using the soNNia model for IgH generation (Figure 2—figure supplement 4). We then draw the size of the lineage at random, using the power-law distribution above. Mutations are then randomly drawn on each sequence of the lineage in a way that preserves the mutation sharing patterns observed in families of comparable size from the partitioned data (Figure 2—figure supplement 5). We thus generated 10⁴ lineages and 2.5 · 10⁴ sequences. Note that, while that procedure is partially based on real data, in particular the distribution of lineage sizes and mutational co-occurence structure in the lineages, it uses completely random sequences and mutations. In addition, these empirical observables were inferred from ${\displaystyle \mathrm{V}\mathrm{J}l}$ classes that were easy to cluster, ensuring that they are not biased by our inference method, and therefore should not give it an unfair advantage. More details about the procedure are given in the Methods.

We applied the HILARy-CDR3 method to this synthetic dataset. The sensitivity achieved at ${\displaystyle t^{\ast }}$ roughly follows and sometimes even outperforms the predicted one ${\displaystyle \hat{s}(t^{\ast })}$ across different values of ${\displaystyle \rho }$ and ${\displaystyle l}$ (Figure 2J, dashed line), validating the approach and the choice of the adaptive high-precision threshold ${\displaystyle t^{\ast }}$ (the discrepancy is due to the fact that ${\displaystyle \mu }$ is assumed to be constant in the prediction, while it varies in the dataset). These results also confirm the poor performance of the method at low prevalences and short CDR3s.

To improve the performance of HILARy-CDR3, we set out to include the phylogenetic signal encoded in the mutation spectrum of the templated regions of the sequences. Two sequences belonging to the same lineage are expected to share some part of the mutational histories, and therefore sequences with shared mutations are more likely to be in the same lineage.

We focus on the template-aligned region of the sequence outside of the CDR3, where we can reliably identify substitutions with respect to the germline. We denote the length of this alignment by ${\displaystyle L}$, so that the total length of the sequence is ${\displaystyle l+L}$. For each pair of sequences, we define ${\displaystyle n_1,n_2}$ as the number of mutations along the templated alignment in the two sequences, ${\displaystyle n_0}$ the number of mutations shared by the two, and ${\displaystyle n_L=n_1+n_2-2n_0}$ the number of non-shared mutations. Under the hypothesis of shared ancestry, the ${\displaystyle n_0}$ shared mutations fall on the shared part of the phylogeny, and are expected to be more numerous than under the null hypothesis of independent sequences, where they are a result of random co-occurrence (Figure 3A).

Figure 3 with 2 supplements see all

Download asset Open asset

To balance the tradeoff between the information encoded in the templated part of the sequence and the recombination junction, we can compute characteristic scales for the two variables of interest: the number of shared mutations and the CDR3 distance n. Intuitively, in highly mutated sequences, we can expect substantial divergence in the CDR3. At the same time, the number of mutations in the templated regions would increase, possibly leading to more shared mutations. Conversely, sequences with few or no mutations carry no information in the templated region, but we also expect their CDR3 sequences to be nearly identical. To adapt a clustering threshold to the two variables, we compute their expectations under the two assumptions, and define the rescaled variables

where ${\displaystyle ⟨n{⟩}_{\mathrm{T}}=l(n_L+1)/L}$ is the expected value of ${\displaystyle n}$ under the hypothesis that sequences belong to the same lineage (see Methods), and ${\displaystyle ⟨n_0{⟩}_{\mathrm{F}}=n_1{n}_2/L}$ is the expected value of ${\displaystyle n_0}$ under the hypothesis that they do not. The standard deviations are likewise defined as ${\displaystyle {\sigma }_{\mathrm{T}}(n)=\sqrt{⟨n^2{⟩}_{\mathrm{T}}-⟨n{⟩}_{\mathrm{T}}^2}=(1/L)\sqrt{l(l+L)(n_L+1)}}$ and ${\displaystyle {\displaystyle {\sigma }_{\mathrm{F}}(n_0)=\sqrt{⟨n_0^2{⟩}_{\mathrm{F}}-⟨n_0{⟩}_{\mathrm{F}}^2}=\sqrt{n_1{n}_2/L}}}$ (Methods).

For more than 99% of positive pairs, both sequences are mutated, i.e., ${\displaystyle n_1,n_2>0}$ (Figure 3B). Without loss of sensitivity, we focus on the mutated part of the dataset, since we cannot use ${\displaystyle y}$ for non-mutated sequences. The distributions of ${\displaystyle x^{\mathrm{\prime }}}$ and ${\displaystyle y}$ for positive and negative pairs (Figure 3C and D) are well separated, with positive pairs characterized by an overrepresentation of shared mutations. By adding the phylogenetic signal ${\displaystyle y}$ we can identify positive pairs of sequences that have significantly diverged in their CDR3 (${\displaystyle x^{\mathrm{\prime }}>0}$) but share significantly more mutations than expected (large ${\displaystyle y}$).

Computing ${\displaystyle y}$ for each pair of sequences is computationally expensive. To avoid examining all pairs, we first perform two different nested clusterings of each ${\displaystyle \mathrm{V}\mathrm{J}l}$ class using the CDR3-based method: the previously described HILARy-CDR3 ‘fine’ partition with threshold ${\displaystyle t^{\ast }}$ that ensures high precision ${\displaystyle \hat{\pi }=99\mathrm{\%}}$; and a ‘coarser’ clustering with a high threshold ${\displaystyle t=t_{\mathrm{s}\mathrm{e}\mathrm{n}\mathrm{s}}}$ that ensures high estimated sensitivity ${\displaystyle \hat{s}=90\mathrm{\%}}$ (Methods and Figure 3E). When lineages are easily separable (e.g. for sufficiently large prevalence ${\displaystyle \rho }$ and CDR3 length ${\displaystyle l}$), these two partitions coincide, and we do not need to compute ${\displaystyle y}$ at all. When they do not coincide, we can use the phylogenetic signal ${\displaystyle y}$ to refine the coarse high-sensitivity partition. We only need to compute ${\displaystyle y}$ for pairs that belong to the same coarse cluster, but not to the same fine cluster: the phylogenetic signal ${\displaystyle y}$ is used to merge the fine-partition clusters into clonal families (Methods and Figure 3—figure supplement 1). This allows us to considerably reduce the number of pairwise comparisons that we need to make between the templated regions of the sequences.

Using ${\displaystyle x^{\mathrm{\prime }}}$ and ${\displaystyle y}$, we classify pairs of sequences as positive (i.e. belonging to the same family) if ${\displaystyle y\ge x^{\mathrm{\prime }}-t^{\mathrm{\prime }}}$, and as negative otherwise. We can compute the expected sensitivity on the synthetic data, and find that it reaches values ≥90% across the whole range of prevalence ${\displaystyle \rho }$ and CDR3 lengths ${\displaystyle l}$, outperforming HILARy-CDR3 in the low-${\displaystyle \rho }$, low-${\displaystyle l}$ region (Figure 3F and G). This proves that using the phylogenetic signal significantly improves performance over HILARy-CDR3.

The procedure outlined above, which we call HILARy-full, may be summarized as follows: (1) group sequences by ${\displaystyle \mathrm{V}\mathrm{J}l}$ class; (2) apply HILARy-CDR3 twice, once with the high-precision threshold as before to get a fine partition, and once with a high-sensitivity threshold to get a coarse partition, thus obtaining two nested partitions; (3) compute ${\displaystyle x^{\mathrm{\prime }}}$ and ${\displaystyle y}$ using Equation 4 only for pairs that belong to the same coarse cluster but to different fine clusters; (4) merge all fine clusters with at least one pair with ${\displaystyle y\ge x^{\mathrm{\prime }}-t^{\mathrm{\prime }}}$.

We compare our approach to state-of-the-art methods. In addition to our two algorithms—HILARy-CDR3 and HILARy-full—our benchmark includes the alignment-free method of Lindenbaum et al., 2021, partis (Ralph and Matsen, 2016), and the spectral clustering method of SCOPer (Nouri and Kleinstein, 2018). The SCOPer method using V and J gene mutations (Nouri and Kleinstein, 2020) was also tested, but gave worse results (Figure 4—figure supplement 1). Details about the used versions and parameters are referenced in the data availability section. We tested all algorithms on two synthetic datasets: a dataset simulated by the partis package and used in Ralph and Matsen, 2022, to benchmark partis against increasing levels of somatic hypermutations, and the synthetic data described above. That dataset is more realistic in the sense that it represents well the statistics of mutation patterns and, perhaps more importantly, the long-tail distribution of clone sizes observed in the data, with its large impact on the diversity of prevalences, which play an important role in the inference. The partis dataset is generated from a population genetics model. It provides a more independent test since it is not based on data used to develop the method and allows to study performance across different mutation rates.

First, we measure the inference time of each algorithm on our synthetic dataset. We find that the inference time is primarily affected by the size of the largest ${\displaystyle VJl}$ class. Therefore, we measure the inference time using the largest class found in donor 326651 of Briney et al., 2019, with the size of ${\displaystyle N=1.2\times {10}^5}$ unique sequences. We then apply the methods to a series of subsamples of this class to get the computational time as a function of the subsample size (Figure 4A). We only allowed for runtimes below 1 hr. We find that only three methods achieve satisfactory performance (under an hour): the two methods introduced here, and the alignment-free method. The other two methods, SCOPer and partis, are limited to ${\displaystyle \mathrm{V}\mathrm{J}l}$ classes of small size (${\displaystyle <{10}^4}$ and ${\displaystyle <{10}^3}$, respectively).

Figure 4 with 2 supplements see all

Download asset Open asset

To compare the five algorithms in finite time, we test the accuracy of the methods using synthetic datasets with different CDR3 lengths, and with fixed mutation rate of 10% for the partis dataset (the mutation rate is not adjustable in our synthetic dataset as it mimics that of the data). We focus on short CDR3s, ${\displaystyle l\in [15,45]}$, which are the most challenging for lineage inference. Clonal families with longer CDR3s are easy to reconstruct, and simple methods such as single-linkage clustering with a threshold on mutational distance already work very well (Figure 4—figure supplement 2A–C). Each dataset contains 10⁴ unique sequences, so that the dominant ${\displaystyle \mathrm{V}\mathrm{J}l}$ class is typically of size ∼10³ and can be handled by all five algorithms. We measure performance using three metrics applied to the resulting partition: pairwise sensitivity ${\displaystyle s_{\text{post}}}$ (Figure 4B and E), pairwise precision ${\displaystyle {\displaystyle {\pi }_{\text{post}}}}$ (Figure 4C and F), and the variation of information ${\displaystyle v}$ (Figure 4D and G). Performance measures as a function of mutation rate in the partis dataset are presented in Figure 5. Pairwise sensitivity ${\displaystyle s_{\text{post}}}$ and precision ${\displaystyle {\displaystyle {\pi }_{\text{post}}}}$ are a posteriori analogs of the a priori estimates defined before in Equations 1 and 3, now computed after propagating links through the transitivity rule of single-linkage clustering. Their value reflects not only the accuracy of the adaptive threshold but is also affected by the propagation of errors in single-flinkage clustering. Variation of information is a global metric of clustering performance which measures the loss of information from the true partition to the inferred one, and is equal to zero for perfect inference and positive otherwise (Methods).

Figure 5

Download asset Open asset

Out of the five tested methods, only HILARy-full achieved both high sensitivity and high precision across all CDR3 lengths and for both synthetic datasets. HILARy-full is the only method reaching both high precision and sensitivity for CDR3s shorter than or equal to 30 nucleotides (Figure 4B), which corresponds to ∼10% of a typical repertoire of productive IgGs (Figure 7A, inset).

The HILARy-CDR3 method achieves high precision everywhere by construction, but only reached good sensitivity for CDR3 lengths 24 and above. The alignment-free method also achieves high precision everywhere, but with low sensitivity, meaning that it erroneously breaks up clonal families into smaller subsets. These three methods achieve good precision, thanks to the use of a null model for the negative pairs. On the contrary, SCOPer has excellent sensitivity everywhere but only achieves high precision for large lengths (${\displaystyle l>30}$), suggesting that it erroneously merges short-CDR3 clonal families. Likewise, partis has high sensitivity but loses precision for short CDR3 on our realistic dataset, meaning that many clonal families are erroneously merged again. Note that our definition of precision and sensitivity differs from those used in Ralph and Matsen, 2022, which explains the differences between the performance measures reported here and in Ralph and Matsen, 2022. On the synthetic datasets from Ralph and Matsen, 2022, HILARy-full is the only method achieving high precision and sensitivity across mutation rates (Figure 5).

The variation of information offers a useful summary of the performance (Figure 4D and G and Figure 5C). According to that measure, only HILARy-full performs well across CDR3 lengths, mutation rates, and datasets. In particular, HILARy presents a clear advantage for challenging regions of the parameter space, such as CDR3 lengths below 30 nucleotides, and mutation rates of 20% and more.

For a typical repertoire, performance can be summarized into a global score by averaging over all CDR3 lengths in proportion of their abundance (assuming inference is perfect for CDR3 lengths larger than 45 regardless of the method). In this task, HILARy-full achieved 99.9% precision and 98.5% sensitivity; partis 93% and 96.9%; and SCOPer 96.6% and 99.1%. These scores are high because only a minority of lineages are difficult to infer. Nonetheless, HILARy-full provides a substantial gain in precision, while SCOPer presents a slight advantage in sensitivity.

We conclude that HILARy-CDR3 should be chosen for its consistently high sensitivity, specificity, and speed. In the case of the largest datasets, the faster HILARy-CDR3 is a useful alternative for long enough CDR3s in realistic repertoires.

We added an extension of HILARy to infer lineages from paired-chain repertoires, i.e., with paired light- and heavy-chain sequences. To extend HILARy-CDR3, we generalize the ${\displaystyle \mathrm{V}\mathrm{J}l}$ class to a ${\displaystyle V_H{J}_H{V}_L{J}_L{l}_{H+L}}$ class, using the V and J genes from both the heavy and light chains, and the sum of their CDR3 lengths ${\displaystyle l_{H+L}=l_H+l_L}$. We then apply HILARy-CDR3 using the sum of the Hamming distances between the heavy- and light-chain CDR3s, normalized by ${\displaystyle {\displaystyle l_H+{\displaystyle l_L}}}$, as our new paired-chain ${\displaystyle x}$. The null distribution used is computed with soNNia using a default generation model for paired heavy and light chains. We incorporate the phylogenetic signal of both chains by concatenating their respective template genes, to obtain the total mutation counts ${\displaystyle n_0=n_{0,H}+n_{0,L}}$, and using ${\displaystyle L}$ as the sum of the lengths of the ${\displaystyle V_H}$ and ${\displaystyle {\displaystyle V_L}}$ genes.

In Figure 6A–C we compare our method to SCOPer and partis on the synthetic dataset from Ralph and Matsen, 2022, as a function of CDR3 length, as our method for generating synthetic sequences could not be easily extended to add random light chains. Performance comparison as a function of mutation rate is presented in Figure 6D–F. HILARy performs better than SCOPer and comparably to partis, which was designed and tested against this dataset.

Figure 6

Download asset Open asset

We next use our method to infer the clonal families of the heavy-chain IgG repertoires of healthy donors from Briney et al., 2019. Figure 7 summarizes key properties of the inferred clonal families of donor 326651. We take advantage of the consistency of our method across CDR3 lengths, as evidenced by the benchmark, to study how the lineage structure changes with the CDR3 length. To this end, we divide the dataset into nine quantiles, each containing ∼10% of the total number of sequences (Figure 7A, inset).

Figure 7

Download asset Open asset

We find that across the nine subsets of the data, the statistics of the lineage structure inferred with the mutations-based method are largely universal. The distribution of the clonal family sizes ${\displaystyle z}$ (Figure 7A) follows a power law across all CDR3 lengths under study, with no significant differences between different lengths. This results generalizes an earlier observation used above for generating synthetic datasets, but which was restricted to high-${\displaystyle \hat{\rho }}$, high-${\displaystyle l}$ ${\displaystyle \mathrm{V}\mathrm{J}l}$ classes, and justifies a posteriori the use of a universal power law in the generative model.

For the largest families, of size ${\displaystyle z\ge 100}$, we compute two intra-lineage summary statistics: the site frequency spectrum, which gives the distribution of frequencies of point mutations within lineages, and the distribution of ${\displaystyle d_N/d_S}$ ratios between non-synonymous and synonymous CDR3 polymorphisms within clonal families (estimated by counting). To avoid the bias of the varying family sizes, we subsampled all families to size ${\displaystyle z=100}$.

Under models of neutral evolution with fixed population size, the distribution of point-mutation frequencies ${\displaystyle \nu }$ goes as ${\displaystyle {\nu }^{-1}}$. Here, we observe a non-neutral profile of the spectrum, with an upturn at large allele frequencies ${\displaystyle \nu >0.5}$ (Figure 7B). It is a known signature of selection or of rapid clonal expansion (Horns et al., 2019; Nourmohammad et al., 2019). We find that site frequency spectra are universal for all CDR3 lengths, suggesting that the dynamics that give rise to the structure of lineages and the subsequent dynamics that influence the sampling of family members do not depend on the CDR3 length.

The lineage ${\displaystyle d_N/d_S}$ ratio is also largely consistent across CDR3 lengths (Figure 7C), while spanning two orders of magnitude, suggesting a wide gamut of selection forces. We could have expected longer loops to be under stronger purifying selection (lower ${\displaystyle d_N/d_S}$) to maintain their specificity and folding. Instead, we observe that short CDR3s have more lineages with low ${\displaystyle d_N/d_S}$. This may be due to different sequence context and codon composition in short versus long CDR3s. Short junctions are largely templated, whereas long junctions have long, non-templated insertions, and it was shown that templated regions have evolved their codons to minimize the possibility of non-synonymous mutations (Saini and Hershberg, 2015), which would lead to a lower ${\displaystyle d_N/d_S}$, regardless of selection.

We focus the analysis high-throughput RNA sequencing data of IgH-coding genes (Briney et al., 2019). The sequences were barcoded with unique molecular identifiers (UMIs) to correct for the PCR amplification bias and correct sequencing errors. We aligned raw sequences using presto of the Immcantation pipeline (Vander Heiden et al., 2014) with tools allowing for correcting errors in UMIs and deal with insufficient UMI diversity. Reads were filtered for quality and paired using default presto parameters. We selected only sequences aligned with the IgG primer and therefore the lineage analysis is limited to the IgG subset of the repertoire. Preprocessed data was then aligned to V, D, and J templates from IMGT (Giudicelli et al., 2006) database using IgBlast (Ye et al., 2013). After processing, all UMI count information is discarded and only unique nucleotide sequences are kept for further analysis.

Pairs of sequences stemming from the same VDJ recombination are expected to have the same CDR3 length ${\displaystyle l}$ and align to the same V and J templates. An exception could be caused by a insertion or deletion within the CDR3 that would alter its length as a result of the somatic hypermutation process. Such indel events are rare and generally selected against (Lupo et al., 2022), therefore in what follows we shall assume the effect of these events is negligible. The inference could also be affected by the misalignment of either V or J templates but we previously found the effect of alignment errors to be insignificant for identifying VJ classes (Spisak et al., 2020) (the alignment of the D template is error-prone and unreliable, hence not used in the inference procedure). Importantly, the two simplifications described here would result in decreased sensitivity of inference but are not expected to affect its precision.

The extraordinary diversity of VDJ rearrangements can be efficiently described and quantified using probabilistic models of the recombination process as well as subsequent purifying selection. Sequence-based models can assign to each receptor sequence ${\displaystyle s}$, its total probability of generation, ${\displaystyle P_{\text{gen}}(s)}$ (Murugan et al., 2012; Elhanati et al., 2015; Marcou et al., 2018) as well as a selection factor ${\displaystyle Q(s)}$, inferred so as to match frequencies ${\displaystyle P_{\text{data}}(s)}$ of the sequences with a model-based distribution (Elhanati et al., 2014; Sethna et al., 2020; Isacchini et al., 2021)

The ${\displaystyle P_{\text{gen}}}$ model was inferred using unmutated out-of-frame sequences from Briney et al., 2019, using the IGoR software (Marcou et al., 2018). The selection function ${\displaystyle Q}$ model was learned using unmutated productive IgM sequences from Briney et al., 2019, using the soNNia software (Isacchini et al., 2021).

The post-selection distribution ${\displaystyle P_{\text{post}}}$ describes the diversity of the CDR3 regions and in doing so provides an expectation of pairwise distances between unrelated, independently generated sequences of same length ${\displaystyle l}$ (Isacchini et al., 2021). As the soNNia software does not include somatic hypermutations, the underlying assumption is that additional diversity on the CDR3 caused by hypermutations doesn’t affect the distribution of pairwise distances. This assumption is justified by the quality of the fit. We can define

where ${\displaystyle |s_1-s_2|}$ stands for (Hamming) distance between sequences ${\displaystyle s_1}$ and ${\displaystyle s_2}$. This definition of the null distribution is a straightforward recipe for its estimation using (Monte Carlo) samples from ${\displaystyle P_{\text{post}}}$.

Should ${\displaystyle P_{\text{post}}}$ differ significantly from the empirical frequencies ${\displaystyle P_{\text{data}}}$ one can resolve to the following alternative

the equivalent of the negation distribution as defined in Lindenbaum et al., 2021, and used in our evaluation of the alignment-free method (Lindenbaum et al., 2021) in the method benchmark analysis.

Pairwise prevalence is defined as the ratio of pairs of related sequences to the total number of pairs of sequences in a given set. Related sequences share an ancestor and have diverged by independent somatic mutations, post-recombination. Low prevalence can be a major difficulty for any inference procedure as any misassignment (or fallout) will result in a drastic loss of sensitivity or precision. It is instrumental to have an a priori estimate of pairwise prevalence before the families are identified.

To estimate the prevalence from the distribution of distances ${\displaystyle P(n)}$ for a given set of sequences (typically a ${\displaystyle \mathrm{V}\mathrm{J}l}$ class or ${\displaystyle l}$ class), we propose the following expectation-maximization procedure. We stipulate the distribution in question is a mixture distribution of two components, ${\displaystyle P_{\mathrm{F}}(n)}$, the expectation for unrelated sequences defined as above, and ${\displaystyle P_{\mathrm{T}}(n)}$, describing related sequences, modeled using a Poisson distribution

where μ is the mean divergence per base pair. If a particular CDR3 length ${\displaystyle l}$ is represented by unusually large number of ${\displaystyle \mathrm{V}\mathrm{J}l}$ classes, the resultant shape of the positive distribution is often closer to a geometric profile, and is then modeled using ${\displaystyle P_{\mathrm{T}}(n)=(1-M)M^n}$, where ${\displaystyle M=\frac{1}{1+\mu l}}$. In sum

In a standard fashion, we proceed iteratively by calculating the expected value of the log-likelihood (pairs of sequences indexed by i)

where the membership probabilities are defined as

We then find the maximum

and iterate the expectation and maximization steps until convergence, ${\displaystyle |{\rho }_{t+1}-{\rho }_t|<ϵ}$, to obtain ${\displaystyle \hat{\rho }={\rho }_{t+1}}$.

Results for largest ${\displaystyle \mathrm{V}\mathrm{J}l}$ class within each ${\displaystyle l}$ class can be found in Figure 2—figure supplement 3 and results for ${\displaystyle l}$ classes using a geometric distribution can be found in Figure 2—figure supplement 6. Dependence of maximum likelihood prevalence estimates ${\displaystyle \hat{\rho }}$ on class size ${\displaystyle N}$ is plotted in Figure 2—figure supplement 7.

The standard method for CDR3-based inference of lineages proceeds through single-linkage clustering with a fixed threshold on normalized Hamming distance divergence (fraction of differing nucleotides) (Kepler, 2013; Uduman et al., 2014; Yaari and Kleinstein, 2015; Nourmohammad et al., 2019). This crude method suffers from inaccuracy as it loses precision in the case of highly mutated sequences and junctions of short length (see Figure 4—figure supplement 2). If junctions are stored in a prefix tree data structure (Knuth, 2013) single-linkage clustering can be performed without comparing all pairs and hence is typically orders of magnitude faster than alternatives. The prefix tree is a search tree constructed such that all children of a given node have a common prefix, the root of the tree corresponding to an empty string, and leaves corresponding to unique sequences to be clustered. To find neighbors of a given sequence it suffices to traverse the prefix tree from the corresponding leaf upward and compute the Hamming distance at branchings. This method limits the number of unnecessary comparisons and greatly improves the speed of Hamming distance-based clustering (Boytsov, 2011). We implement the prefix tree structure to accommodate CDR3 sequences. Briefly, all the CDR3 sequences of identical length are stored in the leaves of a prefix tree (Navarro, 2001; Boytsov, 2011), implemented as a quaternary tree where each edge is labeled by a nucleobase (A, T, C, or G). The neighbors of a specific sequence are found by traversing the tree from top to bottom, exploring only the branches that are under a given Hamming distance from the sequence. Clusters are obtained by iterating this procedure and removing all the neighbors from the prefix tree until no sequences remain. The package is coded in C++ with a Python interface and is available independently. The time performance of this method for high-sensitivity and high-specificity partitions is studied as a part of the method benchmark analysis.

We take advantage of the speed of a prefix tree-based clustering to perform single-linkage clustering. Besides the algorithmic speed-up afforded by the prefix tree, the difference with previous methods is that we use an adaptive threshold. For any dataset, we define two CDR3-based partitions, high-sensitivity and high-precision clustering, corresponding to two choices of threshold.

The high-precision partition is obtained by setting the threshold ${\displaystyle t}$ to ${\displaystyle t_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}}^{\ast }}$ as the largest ${\displaystyle t}$ such ${\displaystyle \hat{\pi }(t)\le {\pi }^{\ast }}$, with ${\displaystyle {\pi }^{\ast }=0.99}$ (99% precision), where ${\displaystyle \hat{\pi }(t)}$ is given by Equation 1–3. To get the high-sensitivity partition, we set the threshold to ${\displaystyle t_{\mathrm{s}\mathrm{e}\mathrm{n}\mathrm{s}}^{\ast }}$, the smallest ${\displaystyle t}$ such that ${\displaystyle \hat{s}(t)\ge s^{\ast }}$, where ${\displaystyle s^{\ast }=0.9}$ (90% sensitivity), where ${\displaystyle \hat{s}(t)}$ is given by Equation 3.

We apply these thresholds to the single-linkage clustering described above to generate the precise and sensitive partitions, which are then used by the mutations-based method to find an optimal partition that merges the fine clusters within the coarse clusters (Methods and Figure 3—figure supplement 1). We refer to the high-precision partition from the CDR3 alone as HILARy-CDR3, and the mutation-based method as HILARy-full.

Finally, the structure of families leads to propagation of errors that lowers the precision with respect to the a priori estimate ${\displaystyle \hat{\pi }}$. Denoting family size as ${\displaystyle z}$, one error accounted for in ${\displaystyle {\displaystyle \widehat{\mathrm{F}\mathrm{P}}}}$ causes, on average, ${\displaystyle ⟨z{⟩}^2-1}$ extra errors by merging two families. If the a priori precision ${\displaystyle \hat{\pi }}$ is high, we can neglect the second order effect of these two families simultaneously affected by other ${\displaystyle {\displaystyle \widehat{\mathrm{F}\mathrm{P}}}}$ pairs. Therefore the expected precision (Equation 1) of the resulting partition reads

where we assumed ${\displaystyle \hat{s}\simeq 1}$. For ${\displaystyle \hat{\pi }=99\mathrm{\%}}$ and ${\displaystyle ⟨z⟩\simeq 2}$ this formula gives ${\displaystyle ⟨{\pi }_{\text{post}}⟩\simeq 97\mathrm{\%}}$.

To generate synthetic data we make use of the statistics of tree topologies of the lineages identified in the high-sensitivity and high-precision regime of CDR3-based inference from the data (yellow region above the black line in Figure 3F). We denote the set of these lineages by ${\displaystyle \mathcal{L}}$. We assume that to good approximation the mutational process and the selection forces that shaped the mutational landscape in these lineages do not depend on the CDR3 length.

To test the performance of different inference methods across CDR3 lengths, we build synthetic datasets of fixed length.

In the first step, we choose the number of families ${\displaystyle N}$. We then draw ${\displaystyle N}$ independent family sizes from the family size distribution of the form observed in healthy datasets

where ${\displaystyle Z_{\alpha }=\sum _{z\ge 1}{z}^{-\alpha }=\zeta (\alpha ,1)}$. In the next step, we assign a naive progenitor to each lineage by sampling from the ${\displaystyle P_{\text{post}}}$ distribution, selecting sequences with a prescribed length ${\displaystyle l}$ (Figure 2—figure supplement 4). We then choose a lineage in the set of reconstructed lineages ${\displaystyle \mathcal{L}}$ at random among lineages of size ${\displaystyle z}$ (or, for large sizes, the lineage of the closest size smaller than ${\displaystyle z}$). To create a lineage with the same mutation patterns as the real data, we then identify all unique mutations in the lineage from ${\displaystyle \mathcal{L}}$ using standard alignment and tree recontruction methods described in Spisak et al., 2020, and for each mutation denote the labels of members of the lineage that carry it. For each mutation, this defines a configuration of labels, one of ${\displaystyle 2^z-1}$ possible. We subsequently loop through observed configurations and choose new positions for all mutations to apply them to the synthetic progenitors of the ancestor, using the position- and context-dependent model of Spisak et al., 2020. The number of mutations assigned to a given configuration is rescaled by a factor ${\displaystyle \frac{L+l}{L_0}}$ where ${\displaystyle L}$ is the templated length of the synthetic ancestral sequence and ${\displaystyle L_0}$ is the templated length of the model lineage from ${\displaystyle \mathcal{L}}$.

This way a synthetic lineage preserves all properties of the lineages of long CDR3s found in the data, particularly the mutational spectra (Figure 2—figure supplement 5) except for the ancestral sequences and the identity of mutations.

We compute the expected distributions of the CDR3 Hamming distance ${\displaystyle n}$, and the number of shared mutations ${\displaystyle n_0}$, under a uniform mutation rate assumption. In other words, we assume that the probability that a given position was mutated, given a mutation happened somewhere in a sequence of length ${\displaystyle L}$, equals ${\displaystyle L^{-1}}$ (we know this not to be true, see, e.g., Spisak et al., 2020, but it allows for simple computations). It follows that the probability that a given position has not mutated once in a series of ${\displaystyle n}$ mutations is ${\displaystyle (1-L^{-1}{)}^n}$.

In this section, we introduce the variation of information ${\displaystyle v}$, used for evaluating alternative methods for clonal family inference in the benchmark analysis. It is a useful summary statistic to quantify the performance of inference as it is affected by its precision as well as sensitivity (Brown et al., 2007). Variation of information ${\displaystyle v(r,r^{\ast })}$ measures the information loss from the true partition ${\displaystyle r^{\ast }}$ to the inference result ${\displaystyle r}$ (Zurek, 1989; Meilă, 2003). To define the variation of information we first introduce the entropy ${\displaystyle S(r)}$ of a partition ${\displaystyle r}$ of ${\displaystyle N}$ sequences into clusters ${\displaystyle c}$ as

where ${\displaystyle n(c)}$ denotes the number of sequences in cluster ${\displaystyle c}$. The mutual information between two partitions ${\displaystyle r}$ and ${\displaystyle r^{\ast }}$ can then be computed as

where ${\displaystyle n(c,c^{\ast })}$ denotes the number of overlapping elements between cluster ${\displaystyle c}$ in partition ${\displaystyle r}$ and cluster ${\displaystyle c^{\ast }}$ in partition ${\displaystyle r^{\ast }}$. Finally, variation of information is given by

Variation of information is a metric in the space of possible partitions since it is non-negative, ${\displaystyle v(r,r^{\ast })\ge 0}$, symmetric, ${\displaystyle v(r,r^{\ast })=v(r^{\ast },r)}$, and obeys the triangle inequality, ${\displaystyle v(r_1,r_3)\le v(r_1,r_2)+v(r_2,r_3)}$ for any three partitions (Zurek, 1989).

We used version 1.2.0 for spectral SCOPer, 1.3.0 for SCOper using the V and J mutation presented in Figure 4—figure supplement 1, version 1.2.0 for HILARy, version 0.16.0 for partis, and the code from this repository https://bitbucket.org/kleinstein/projects/src/master/Lindenbaum2020/Example.ipynb for the alignment-free method. The HILARy tool with Python implementations of the CDR3 and mutation-based methods introduced above can be found at https://github.com/statbiophys/HILARy (copy archived at Athènes, 2024). The standalone prefix tree implementation can be found at https://github.com/statbiophys/ATrieGC (copy archived at Dupic, 2024). A complete guide to our benchmark procedure can be found in the README of the folder https://github.com/statbiophys/HILARy/tree/main/data_with_scripts where we make available scripts to infer lineages and reproduce the benchmark figures of this article. We also upload this folder with all input and output data at https://zenodo.org/records/10676371.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

The study was supported by European Research Council COG 724208 and ANR-19-CE45-0018 ‘RESP- REP’ from the Agence Nationale de la Recherche grants and DFG grant CRC 1310 ‘Predictability in Evolution’.

© 2024, Spisak et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.