Large-scale data mining of four billion human antibody variable regions reveals convergence between therapeutic and natural antibodies that constrains search space for biologics drug discovery

ABSTRACT

The naïve human antibody repertoire has theoretical access to an estimated > 10¹⁵ antibodies. Identifying subsets of this prohibitively large space where therapeutically relevant antibodies may be found is useful for development of these agents. It was previously demonstrated that, despite the immense sequence space, different individuals can produce the same antibodies. It was also shown that therapeutic antibodies, which typically follow seemingly unnatural development processes, can arise independently naturally. To check for biases in how the sequence space is explored, we data mined public repositories to identify 220 bioprojects with a combined seven billion reads. Of these, we created a subset of human bioprojects that we make available as the AbNGS database (https://naturalantibody.com/ngs/). AbNGS contains 135 bioprojects with four billion productive human heavy variable region sequences and 385 million unique complementarity-determining region (CDR)-H3s. We find that 270,000 (0.07% of 385 million) unique CDR-H3s are highly public in that they occur in at least five of 135 bioprojects. Of 700 unique therapeutic CDR-H3, a total of 6% has direct matches in the small set of 270,000. This observation extends to a match between CDR-H3 and V-gene call as well. Thus, the subspace of shared (‘public’) CDR-H3s shows utility for serving as a starting point for therapeutic antibody design.

KEYWORDS:

• CDR-H3
• database
• repertoire

Introduction

Antibodies are the largest group of biotherapeutics, with nearly 200 approved in at least one country¹ (https://www.antibodysociety.org/antibody-therapeutics-product-data/) and comprising more than a quarter of new approvals by the US Food and Drug Administration in 2022.² These therapeutics derive from the successful harnessing of the natural process wherein an organism rapidly develops an antibody binder against a pathogen. The organism achieves this by a priori having a large number of naïve antibodies (~10¹¹),³ at least one of which ideally binds the pathogen and starts an immune response. Thereafter, in germinal centers, during the somatic hypermutation process, the organism explores an even larger (>10¹⁵) theoretically possible space of these molecules to develop an ever better binder. Understanding how this immense diversity is explored is crucial to our understanding of the immune system and developing better monoclonal antibody therapies.

The great diversity was proposed to be the result of accruing variability in response to pathogens, accumulated on an evolutionary scale by vertebrates.⁴ However, this number is not plausible when taking into account the total number of epitopic sites that an organism would need to address within a reasonable time to mount an immune response.⁵ It was noted that, despite theoretically great diversity, amino acid preferences in the complementarity-determining region (CDR)-H3s (the most variable loop) are biased.⁶ The sequence space that is theoretically available to antibodies is ‘astronomical’, with some estimates placing it at as many as 10¹⁸ possible molecules.^5,7,8 With such large numbers, it appears implausible that an organism would explore all the possibilities, suggesting that a biased sampling of this space is more likely than purely random sampling.

The theoretically accessible sequence repertoire might indeed be theoretically ‘astronomical’ in size, but for practical reasons, there seems instead to be another, smaller, more manageable dimension that is being explored.⁵ However, productive and nonproductive sequences appear to have similar profiles of somatic hypermutation, suggesting a stronger role in the random biochemical process rather than strategic selection.⁹ To address this conundrum, theories of ‘minimal’ repertoires were proposed decades ago.¹⁰ The structural dimension was proposed as one explanation, as there appears to be more conservation in structural configurations than there are sequences of antibodies.^11–13 Nevertheless, even at the most basic, sequence level, certain commonalities between independent repertoires are observed.

Despite the great variability available to repertoires, different individuals can evolve the same antibodies toward the same antigen challenge.¹⁴ Commonality exists between repertoires of individuals as well, as documented by deep sequencing studies of three⁷ and ten individuals.⁸ In both cases, the prevalence of clonotypes (same V-J assignment and CDR-H3 amino acid sequence) between different individuals was studied. In the Soto et al. study, 0.3% of heavy chain clonotypes were shared by all three individuals in the study (1–6% were shared between any two individuals). In the Briney et al. study, 0.022% of clonotypes were shared by all ten subjects in the study (on average 0.95% between any two individuals). Therefore, despite great diversity available to antibody repertoires, independent natural samples draw same/similar sequences, suggesting that the space is explored not randomly, but in a directed manner.

While antibodies are used as biotherapeutics, developing a biologic diverges from natural antibody processes. Monoclonal antibody harnessing originated in 1975 with Kohler and Milstein’s hybridoma technology, which involved fusing immunized mouse lymphocytes with myeloma cells. Yet, using mouse antibodies in therapy may trigger an immune response, generating anti-murine antibodies and reducing efficacy. Strategies such as chimerization, humanization, and use of fully human sequences aim to minimize this response.

Throughput limitations of hybridoma technology prompted the development of other selection methods, such as display systems producing fully human immunoglobulins that resemble adaptive immune response in vitro. Typically in this approach, antibody-binding domains such as antigen-binding fragments or single-chain variable fragments are presented on the outside of the cell wall or membrane. Both processes need to be followed by lead optimization.¹⁵

Despite divergent development processes, therapeutic sequences can be found in naturally sourced repertoires.^16,17 Position-specific scoring matrices developed on the basis of natural repertoires aligned with mutations selected for in therapeutics.¹⁸ Therefore, despite the great diversity available to antibodies, a set of these developed by independent artificial methods appears to follow similar selection pressures.

Convergence between natural and therapeutic antibodies may facilitate drug discovery. Next-generation sequencing (NGS) repertoires are growing in numbers, providing an unprecedented sampling of the vast theoretical space available to antibodies. Drawing therapeutic molecules from this space is desirable, as human-origin antibodies can produce fewer downstream issues than non-human or display platforms.¹⁹ Mapping the characteristics of the observed convergence between natural and therapeutic antibodies could potentially be used to limit the space of natural sequences that are used for drug discovery, accelerating the process.

Here, to interrogate the sampled natural antibody space, we automatically data-mined the public NGS repositories to create a heterogeneous set of naturally sourced B-cell receptors that is an order of magnitude larger than previously published datasets.^20–22 We used this dataset to identify commonalities between it and therapeutic antibodies.

Results

Large-scale data mining public next-generation sequencing repositories for B-cell receptor repertoires

We automatically mined the SRA repositories, to identify bioprojects associated with B-cell receptors (BCRs). For context, when conceived in 2018, the Observed Antibody Space (OAS), which has been a widely used resource for BCR NGS data,^23,24 was created using manual identification of BCR-containing sequences. This resulted in 55 bioprojects, which was recently updated to 90 bioprojects.²⁰ Here, using an automatic identification pipeline requiring processing of more than 500,000 bioprojects, we identified 287 bioprojects with BCR data. Of these 287, 220 were earmarked for processing because they did not include single-cell data and had more than 10% of reads identifiable as immunoglobulin genes and associated metadata indicated that immunoglobulins were the main focus of the study.

The 220 bioprojects contributed a total of 7 billion sequences that were processed into IMGT-aligned amino acid sequences, allowing for CDR delineation. We processed the data using a single uniform pipeline, creating the AB-NGS database (Table 1). By contrast, previously in OAS, the raw sequences were processed using IgBLAST²⁵ to convert nucleotides to amino acid sequences. The amino acid sequences were IMGT numbered and germlines identified using ANARCI.²⁶ Both pieces of software used IMGT²⁷ as the underlying germline database. It is known, however, that annotation results depend on the underlying germline database and piece of annotation software used in the process.²⁸ For this reason, we adapted best practices from IgBLAST, ANARCI, and MMSeqs²⁹ programs to create a single pipeline that filters and numbers nucleotide sequences in a single pass, using the VDJBase germline database³⁰ (see Methods).

Table 1. Statistics of the number of heavy-chain sequences processed to make the AB-NGS database. We stratified our datasets by the organism and whether the sequence was productive (no stop codons, no truncation, see Methods for more details). The multiple repeats within the bioproject were taken as identical variable region sequences on the nucleotide level to increase the stringency of including only the confident sequences.

	All	Heavy	Light
Raw reads	7,025,951,034	5,203,593,751	1,822,357,283
Productive sequences	6,004,004,237	4,443,879,409	1,560,124,828
Productive human	5,352,603,279	3,948,237,955	1,404,365,324
Unique productive human	2,497,759,415	2,121,082,795	376,676,620
Unique, with multiple (>1) repeats within bioproject	316,064,499	264,597,059	51,467,440

The processing of 7 billion sequences resulted in 6 billion productive sequences (Table 1). The processing allowed us to identify the CDR regions, closest germlines and, wherever available, the isotype (Figure 1). The gene usage across our datasets places some of the most common genes as having the largest coverage,^7,8 so our annotation pipeline recapitulates the data distribution from its constituent experiments.^7,8,14 In some projects, isotype was provided in metadata, but we performed sequence level annotation by aligning constant region of each sequence to C genes database, obtained from NCBI. For a considerable number of sequences, it was impossible to identify the isotype. The most common constant region identification was IGHM. In second place is IGK, with IGHG2 only in third position. Altogether, IGG isotype is typically associated with an active immune response whereas IGM is associated with the naïve repertoire. The constant region identification suggests a predominantly naïve makeup of our database, though with the caveat that in multiple cases identification of the constant region was not possible. This observation, however, is further confirmed by metadata analysis.

Figure 1. Annotated C regions as cumulative plot. We automatically identified the constant region, wherever it was possible. On the basis of the heavy C region call, the isotype can be inferred. For most of the sequences it was not possible to infer the C region (the blank leftmost column).

For each bioproject, we collected and cleaned its associated metadata. The 220 bioprojects were associated with more than 11,000 biosamples, that in turn are linked to author-submitted metadata, such as organism, B-cell type, tissue, and disease. This data format is not uniform, requiring substantial manual processing to achieve normalization across all 11,000 samples. For each sample, where present, we extracted the organism, B-cell type, tissue, and disease and their per-sequence distribution is given in Figure 2.

Figure 2. Metadata landscape of our database. a. Author-annotated organism distribution of sequences in our dataset. b. Source of B-cells in our dataset. c. Type of B-cells in our dataset. d. Author-annotated disease present in the dataset.

The BCR sequences in our database cover a range of organisms, with the most common being human, followed by mouse (Figure 2a). This reflects the therapeutic focus of sequencing efforts that address humans, wherein mice are a typical model organism in immunization campaigns. Other model organisms such as macaques, rats and rabbits are also represented. The presence of camelids (Vicugna pacos, Camelus bactrianus) indicates some studies that target the naturally occurring heavy-chain only antibodies, or VHHs, marketed therapeutically as nanobodies.^31,32

Most of the sequences in our database come from peripheral blood mononuclear cells, followed by spleen and bone marrow (Figure 2b). Most of the sequences come from naïve B-cells, with memory B-cells in the second place (Figure 2c). The majority of sequences in our database come from individuals who were not mounting an immune response to any disease (Figure 2d). Nevertheless, in a large proportion of cases, the disease status information was not available.

Altogether, sequence and metadata analysis indicates that the contents of our database are largely human, naïve, and not responding to any immune challenge. We believe that it is the biggest compilation of publicly available human BCR sequencing data to date. By contrast, the current version of OAS (Oct 2023), held 2.2 billion redundant human antibody sequences from 70 bioprojects, contrasted to our 5.3 billion from 135 bioprojects (Table 1). We thus consider the dataset to be well placed to study commonalities to constrain the great antibody diversity for therapeutic discovery.

Highly public CDR-H3s account for 0.07% of all CDR-H3s and 0.31% of confident CDR-H3s

Drawing from observations that natural antibodies can be sampled by different individuals and therapeutic antibodies can arise naturally, we hypothesized that there is a reduced set of ‘shared’ antibodies within which therapeutic varieties are more likely to be found.

Taking advantage of the large set of independent human studies, we quantified the extent to which the human natural antibody sequence space is sampled in the same way, by checking how many of the sequences are ‘shared’ or ‘public’. To avoid confounding factors such as gene usage and disease states, we focused only on the CDR-H3 region, as the most diverse and instrumental element of an antibody for ease of analysis and interpretation.

We took bioproject as the quantum of an independent sample, as many studies differ in variables such as number of individuals and sequencing regimens. Therefore, any commonalities that result do so despite the diversity in our datasets. To further stratify the diversity, we divided the dataset into the Confident and All sequences. The Confident set was composed of sequences that had to occur at least twice in the source bioproject. This further ensured that the same CDR-H3s observed across different bioprojects did not occur as a result of technical artifacts, such as sequencing errors.

To test for ‘public sequences’, we counted how often the same CDR-H3 occurred among different naturally sourced NGS studies. Our dataset includes data from studies that are known for their depth of sequencing (e.g., Soto et al., Briney et al.,^7,8 which could potentially skew any statistics. All but one bioproject of 135 contributed CDR-H3s to the public datasets, suggesting that no single dataset is a biasing factor. Altogether, we note that when bioproject redundancy of just two is imposed, it captures a greater number of CDR-H3s (1.58% for All and 4.17% for Confident) than for more stringent redundancies (Table 2). Furthermore, the difference between public CDR-H3s from All or Confident datasets becomes blurred with higher redundancies (Table 2). This suggests that, at bioproject redundancy of just two, chance occurrences can be observed. The percentage of public CDR-H3s appears to drop asymptotically with the higher stringencies, appearing to stabilize in the range of one promile (0.07% for redundancy 5 or 0.01% at redundancy 9). The corresponding figure for the confident CDR-H3s is 0.31 at redundancy 5 and 0.05 at redundancy 9. Altogether, the public CDR-H3 levels are within the same range that the public clones (that also included the V call) reported between individuals within the Briney et al. and Soto et al. studies^7,8 (though direct comparison between the numbers is tenuous). Nevertheless, here, we report public sequence occurrences that are drawn from much larger and more heterogeneous samples than either of the two studies.

Table 2. Number of public CDR-H3s by the number of independent bioprojects they occurred in. The ‘Confident’ CDR-H3 dataset is a collection of CDR-H3s from sequences that occurred at least once within each bioproject. By contrast the ‘All’ dataset includes all the sequences without such constraints. Within Confident and All datasets, we noted how many separate bioprojects each CDR-H3 occurred in. The table shows the number of CDR-H3s in each dataset with CDR-H3s occurrences equal to or greater than specified bioproject redundancy. The percentages in both columns refer to different totals, in case of All CDR-H3s it is 385,878,857 sequences and in case of Confident 85,688,290.

Bioproject redundancy	All CDR-H3%)	Confident CDR-H3s (%)
2	6,126,871 (1.58%)	3,577,323 (4.17%)
3	1,410,087 (0.36%)	1,104,768 (1.28%)
4	570,412 (0.14%)	504,856 (0.58%)
5	288,311 (0.07%)	271,147 (0.31%)
6	164,861 (0.04%)	159,848 (0.18%)
7	102,684 (0.02%)	101,010 (0.11%)
8	68,275 (0.01%)	67,709 (0.07%)
9	47,611 (0.01%)	47,390 (0.05%)

Any occurrences of shared CDR-H3s between samples could either be attributed to biological convergence or the presence of common features that increase likelihood of multiple occurrences. The biological convergence is hard to define and difficult to check. Features driving higher probability of occurring between different studies are likely CDR length and sequence composition. To check for both, we plotted the CDR-H3 lengths of All CDR-H3s and the public datasets (Figure 3). The results shows unequivocally that the public CDR-H3s are generally shorter. One would expect increasingly more common CDR-H3s to be increasingly shorter. Nonetheless, increasing the stringency of the number of CDR-H3 occurrences does not reduce the CDR-H3 length. To assess CDR-H3 composition, which we expressed as positional diversity, we plotted the positional entropy for a range of most common CDR-H3 lengths (Figure 4). The results show that All CDR-H3s have slightly higher diversity (entropy) than the public datasets. Nevertheless, the diversity is not reduced radically if samples of, comparatively small, 700 are drawn from the public datasets.

Figure 3. Length distribution of the CDR-H3s across the entire dataset and public datasets. The solid black line shows the length distribution in the entire dataset of 385 m sequences. The public datasets are divided into the Confident ones (dashed lines, subscript ‘r’ in legend) and the ones not requiring redundancy (solid lines, subscript ‘all’ in legend). The numeric subscripts in public datasets indicate the across bioproject redundancy required, e.g., public_all_3, has all the CDR-H3s that were observed across at least three bioprojects.

Figure 4. Diversity of the CDR-H3s in the public datasets. Here we show the Confident (redundant) dataset diversity against the background of all the CDR-H3s (‘All’ dataset, black line). Diversity was calculated as amino acid frequency entropy calculated for each position on length-matched CDR-H3s. The entropy is lower in the in the N-terminal region because of higher conservation (e.g., the AR motif). Each panel shows specific CDR-H3 lengths according to the IMGT definition. As expected, downsampling each public dataset to 700 (the number of therapeutic CDR-H3s we have) and recalculating the entropy reduces the diversity further.

Therefore, our results indicate that public CDR-H3s bear features making them more likely to be observed across multiple datasets: they are shorter and less diverse. We checked whether such a collection of CDR-H3s with common features indeed has any strategic immunological diversity bias that can be used therapeutically or whether it is purely a statistical coincidence.

Public CDR-H3s define a reduced set of clonotypes closely reflecting antibodies derived from therapeutic programs

We have previously noted that very good sequence identity matches within the NGS space can be found for some therapeutic CDR-H3s.^16,17 Though our previous studies indicated that the public NGS for therapeutic sequences can be mined, they did not shed light on patterns in the convergence across naturally sourced and artificially developed molecules. Any immune repertoire selection biases, especially those shared with the therapeutic space, could be made clinical actionable by biasing selection in that direction. To check whether the set of public CDR-H3s, that is a priori ‘shared’ and radically constrains the diversity space, bears any therapeutic actionables, we checked to what extent it overlaps with therapeutic CDR-H3s.

We identified 700 CDR-H3s that are unique within the set of INN-drawn 827 therapeutics, at various stages of development.¹⁶ We next checked the length matched Hamming distance from each of the 700 therapeutic CDR-H3s to one of our ca. 385 m CDR-H3s. For instance, a Hamming distance between sequences POTATOS and TOMATOS is two as they differ in that many length-matched aligned characters. We performed the CDR-H3 matching for the Confident and All datasets and for five CDR-H3 public stringencies (e.g., redundancy 3 indicates that CDR-H3s occurred in at least three different bioprojects).

We plotted the best Hamming distance from our 700 therapeutics to all our varieties of public datasets and the results are shown in Table 3. For a non-trivial proportion of therapeutic CDR-H3s, we find perfect matches at all levels of public CDR-H3 stringency. For Confident public sets at bioproject redundancy of 2, we find perfect matches for a total of 9% of all CDR-H3s. Increasing the redundancy to 5, decreases the amount of public CDR-H3s from 3,577,323 to 271,147 (7% of 3,577,323). However, perfect therapeutic matches at bioproject stringency 2 at 9% decreases to only 6%. Moreover, the 271,147 CDR-H3s from Confident dataset at redundancy 5 have 1 Hamming distance and 2 Hamming distance matches for 10% and 15% of therapeutic CDR-H3s, respectively. Therefore, the relatively small number of 271,147 public CDR-H3s cover ca 30% of all therapeutic CDR-H3s within two substitutions. Similar good correspondences between public and therapeutic CDR-H3s are noted for all bioproject redundancies regardless of whether CDR-H3s are from All or Confident datasets (Table 3). Comparatively, if we take the entire dataset of 385 m CDR-H3s, it covers 67% of all therapeutics up to two substitutions, despite being ca. 1420 larger in size than the public dataset with 271,147 CDR-H3s. Increasing the bioproject redundancy, radically reduces the number of CDR-H3s in a public dataset, and does not proportionally reduce the number of therapeutic matches. This suggests that there might be a non-trivial convergence between public and therapeutic datasets.

Table 3. Number of therapeutic CDR-H3s within a specific Hamming distance from a dataset. We created a range of reference datasets of CDR-H3s that we compare to the 700 therapeutic CDR-H3s using Hamming distance. As references, we created variants of the Confident and All datasets akin to Table 2, where the CDR-H3s were only constrained to those that occurred in a sufficient number of bioprojects (given by bioproject redundancy column). We also created randomizations of each such dataset either by randomizing the CDRs in the NGS dataset (the -D subscript) or the 700 therapeutic CDR-H3s we compared to (−T subscript). For instance, Public All-2 is the subset of CDR-H3s from the ‘All’ dataset, constrained to CDR-H3s occurring in at least two bioprojects and compared to the original unmodified 700 therapeutic CDR-H3s. Random All T-2 is the same dataset, but the 700 therapeutics CDR-H3s are randomized. In Random All D-2, only the original NGS dataset is randomized and therapeutic CDR-H3s are left unmodified. The CDR-H3s are compared by length-matching and hamming distance at five levels, with 0 Hamming distance indicating identity and 4, indicating four differences. Additionally we compared the therapeutic contrast with germline CDR-H3s, in original form from IMGT and applying our randomization scheme to the same germline CDR-H3s.

			Best length-matched Hamming distance
Bioproject redundancy	Sequence Redundancy	DS	0	1	2	3	4
none	none	All CDR-H3s	104 (14%)	193 (27%)	182 (26%)	99 (14%)	58 (8%)
2	All	Public All −2	70 (10%)	114 (16%)	145 (20%)	132 (18%)	97 (13%)
		Random All T-2	19.3 (2%)	44.0 (6%)	72.4 (10%)	106.8 (15%)	117.9 (16%)
		Random All D-2	15.9 (2%)	74.0 (10%)	100.5 (14%)	140.5 (20%)	109.4 (15%)
	Confident	Public Confident-2	63 (9%)	106 (15%)	144 (20%)	138 (19%)	95 (13%)
		Random Confident T-2	16.9 (2%)	41.8 (5%)	72.4 (10%)	98.2 (14%)	119.0 (17%)
		Random Confident D-2	16.0 (2%)	68.1 (9%)	95.4 (13%)	137.5 (19%)	113.4 (16%)
3	All	Public All −3	59 (8%)	85 (12%)	138 (19%)	137 (19%)	99 (14%)
		Random All T-3	16.1 (2%)	37.0 (5%)	58.3 (8%)	89.9 (12%)	112.4 (16%)
		Random All D-3	13.5 (1%)	61.0 (8%)	84.9 (12%)	131.0 (18%)	119.5 (17%)
	Confident	Public Confident-3	58 (8%)	83 (11%)	133 (19%)	140 (20%)	102 (14%)
		Random Confident T-3	14.6 (2%)	37.8 (5%)	58.6 (8%)	88.3 (12%)	112.2 (16%)
		Random Confident D-3	12.9 (1%)	57.0 (8%)	86.1 (12%)	127.3 (18%)	122.7 (17%)
4	All	Public All −4	49 (7%)	81 (11%)	121 (17%)	144 (20%)	105 (15%)
		Random All T-4	12.9 (1%)	32.2 (4%)	54.2 (7%)	84.9 (12%)	110.9 (15%)
		Random All D-4	10.8 (1%)	56.4 (8%)	75.5 (10%)	119.9 (17%)	124.8 (17%)
	Confident	Public Confident-4	49 (7%)	79 (11%)	120 (17%)	144 (20%)	108 (15%)
		Random Confident T-4	10.9 (1%)	32.4 (4%)	54.4 (7%)	80.9 (11%)	110.9 (15%)
		Random Confident D-4	10.8 (1%)	53.0 (7%)	76.9 (10%)	121.5 (17%)	124.2 (17%)
5	All	Public All −5	46 (6%)	75 (10%)	111 (15%)	144 (20%)	109 (15%)
		Random All T-5	11.3 (1%)	28.8 (4%)	51.7 (7%)	74.4 (10%)	107.7 (15%)
		Random All D-5	10.6 (1%)	47.0 (6%)	74.3 (10%)	110.3 (15%)	125.6 (17%)
	Confident	Public Confident-5	46 (6%)	75 (10%)	111 (15%)	144 (20%)	108 (15%)
		Random Confident T-5	12.1 (1%)	28.6 (4%)	51.5 (7%)	75.3 (10%)	108.5 (15%)
		Random Confident D-5	10.1 (1%)	48.9 (6%)	72.7 (10%)	112.6 (16%)	122.4 (17%)
		Germline	0 (0%)	0 (0%)	1 (0%)	4 (0%)	7 (1%)
		Germline random	0.0 (0%)	0.0 (0%)	0.0 (0%)	0.0 (0%)	0.1 (0%)

A non-trivial element that could be contributing to the overlap we observe is the germline nature of the public and therapeutic datasets. It could be the case that the public CDR-H3s are in fact not very diversified away from the germline. To carry out the germline test, we created a set of artificial germline CDR-H3s by combining all possible human V, D and J amino acid segments found in the IMGT database.⁶ We did not account for V(DD)J fusions³³ because of their low frequency measured in peripheral blood. An alternative approach would be to use a simulation framework such as immuneSIM, SoNNIa,^34,35 though we opted for IMGT CDR-H3s because we consider these unencumbered by mutations, representing a germline ground truth.

To check germline-therapeutic overlap, we compared the germline IMGT CDR-H3s and for each therapeutic CDR-H3, calculated the smallest length-matched Hamming distance. The results are shown in Table 3. Only at a Hamming distance of 2 do we start to observe any overlap, suggesting that trivial germline nature is not the driving factor in the convergence we see between public and therapeutic CDR-H3s.

Another factor that could be contributing to the therapeutic overlap with the public database is a bias inherent to the discovery procedure of a given therapeutic. We hypothesized that therapeutics derived from murine sources are less likely to be found as good matches in our highly public dataset of 271,147 CDR-H3s, as opposed to human platforms. Of our 700 therapeutic CDR-H3s, we identified 226 CDR-H3s whose origin can be labeled as human and 312 that can be labeled as humanized. There were some chimeric, murine, and antibodies of other sequence source, but their numbers were much smaller with respect to these two groups so we considered a contrast between the human/humanized as the best way forward. A total of 20 human therapeutic CDR-H3s found direct matches in our highly public dataset as opposed to 19 humanized ones. We used Fisher’s exact test to check whether there is a statistically significant difference between observing direct match within a human or humanized dataset. The test statistic was not significant at a strict significance level of 0.01. Therefore, the development platform or organism source did not appear to bestow features that separated the CDR-H3s between those that are frequently found in our highly public dataset.

Since we established previously that public CDR-H3 datasets have features that increase the likelihood of occurrences among datasets, we checked whether the good correspondence between therapeutic and public datasets is due to chance. We devised a conservative randomization scheme for each CDR-H3 that kept two N-terminal residues and three C-terminal residues constant (as these positions are typically less diverse,³⁶ see Figure 4), whilst uniformly randomizing the rest. Using this scheme, the shorter CDR-length distribution is already factored in. Since we compare the 700-member CDR-H3 therapeutic set versus public datasets of varying sizes, we introduced two dataset randomization schemes, therapeutic-sided and public-sided. In the therapeutic-sided test, we randomize the therapeutic dataset and compare it against the undisturbed reference dataset (e.g., public dataset). In the public-sided randomization, we modify only the reference dataset leaving the therapeutic CDR-H3s undisturbed. We performed ten separate runs in each case of therapeutic-sided and public-sided randomizations, with the results shown in Table 3. Despite quite lenient length-matched randomization scheme, it is clear that randomizations consistently produce much worse matches to therapeutics than undisturbed public datasets. For instance, whereas perfect matches for public datasets are in the range of 6–10% for all therapeutics, they do not go beyond 1–2% for the randomization experiments. It appears that randomizing the public dataset rather than therapeutics produces slightly better results. Though public CDR-H3s have features that are common among different datasets, there is a nontrivial signal driving therapeutic-public dataset convergence that cannot be explained by shorter lengths and random diversity alone.

Though most of the study focused on CDR-H3 alone to avoid the many confounding factors such as isotypes and germlines, selection of antibody-based therapeutics does not hinge on this loop alone. Whilst constructing phage display libraries, a known binder or a framework with known good biophysical properties, such as trastuzumab, might be used.³⁷ For this reason, we also checked whether the public-therapeutic correspondence disappears when information beyond CDR-H3 is introduced, such as the clonotype,^7,8 that encompasses information of the entire variable region. For each public CDR-H3 noticed in a public dataset, we also noted its gene call, either V only or in combination with the J-gene. A match between a therapeutic clonotype was called when the CDR-H3 matched to a specific extent (e.g., Hamming distance 2) and there was a match in gene calls.

The best matches to therapeutic clonotype (V or V+J call and CDR-H3 Hamming distance) within a public dataset are shown in Table 4. For the Confident dataset at bioproject stringency 5, we found perfect CDR-H3 matches for 6% of therapeutic CDR-H3s compared to 5% when clonotype information is required. At bioproject stringency 5, public CDR-H3s were within 2 substitutions of ca. 30% of all therapeutic CDR-H3s, but this decreases to only ca 27% when clonotype information is required. A similar high correspondence between clonotype (Table 4) and CDR-H3 (Table 3) only matches is seen for other bioproject stringencies. This indicates that correspondence between public and therapeutic molecules is not only constrained to CDR-H3s but also to clonotype, which is already information that can be used for drug discovery, e.g., for library design or virtual screening. The set of public CDR-H3s with their assigned clonotypes is much smaller than the number of possible antibodies, but appears to closely reflect diversity of known antibodies discovered via therapeutic programs.

Table 4. Clonotype-level matches of public H3s to therapeutics. As in Table 2 & Table 3, for All and Confident datasets, we constrained them to only those CDR-H3s that occurred in a number of bioprojects specified by the Bioproject redundancy column. The CDR-H3 Hamming distance between a reference and therapeutic is calculated only if their V or V+J calls match between the NGS sequence and a therapeutic sequence.

	Hamming distance	0		1		2		3		4
Bioproject redundancy	Dataset redundancy/V call type	V	V+J	V	V+J	V	V+J	V	V+J	V	V+J
2	All	42 (6%)	38 (5%)	82 (11%)	77 (11%)	130 (18%)	109 (15%)	146 (20%)	123 (17%)	104 (14%)	105 (15%)
	Confident	42 (6%)	38 (5%)	76 (10%)	72 (10%)	131 (18%)	107 (15%)	140 (20%)	119 (17%)	108 (15%)	113 (16%)
3	All	40 (5%)	37 (5%)	70 (10%)	66 (9%)	112 (16%)	89 (12%)	149 (21%)	122 (17%)	111 (15%)	113 (16%)
	Confident	40 (5%)	37 (5%)	67 (9%)	63 (9%)	111 (15%)	87 (12%)	145 (20%)	120 (17%)	116 (16%)	118 (16%)
4	All	38 (5%)	36 (5%)	66 (9%)	59 (8%)	103 (14%)	82 (11%)	144 (20%)	123 (17%)	109 (15%)	109 (15%)
	Confident	38 (5%)	36 (5%)	63 (9%)	56 (8%)	106 (15%)	85 (12%)	142 (20%)	121 (17%)	111 (15%)	111 (15%)
5	All	38 (5%)	36 (5%)	60 (8%)	51 (7%)	100 (14%)	83 (11%)	136 (19%)	117 (16%)	114 (16%)	112 (16%)
	Confident	38 (5%)	36 (5%)	60 (8%)	51 (7%)	100 (14%)	83 (11%)	135 (19%)	115 (16%)	115 (16%)	114 (16%)

Discussion

Therapeutic antibody discovery hinges on our ability to translate the binding diversity of these molecules into developable therapeutics. The task is daunting due to the vast sequence space of approximately 10¹⁸ possible (natural) molecules. In that space, a molecule needs to be found with the right combination of suitable target-binding and manufacturability properties. Antibodies are naturally bestowed with target-binding properties, but not manufacturability. ‘Developable’ antibodies appear to nevertheless form a subset of all available natural antibodies.³⁸

Identifying potentially therapeutic antibodies in this vast space can be facilitated by a combination of an ‘oracle’³⁹ calling binding (or specificity) and developability properties or by reducing the space to more viable molecules. Traditional discovery methods such as display protocols have the benefit of an oracle in the form of a process that filters out binders.^40,41 Display methods have no way of calling developable antibodies, which can lead to multiple downstream issues.¹⁹ A workaround these methods are using is in fact constraining the initial libraries (or search space) to sequences that are more developable, e.g., without any liabilities.⁴¹

Computational binding oracles have been proposed.^37,42,43 Opposed to display approaches that ‘constrain’ the space a priori because of their low-throughput nature, computational methods promise to enable exploration of an almost unconstrained set of antibodies. In theory, by being able to perform binding and developability checks quickly computationally, all the combinatorial combinations of sequences could be checked. In practice, there is currently no universal ‘binding oracle’ available, with published approaches relying on a priori deep mutational scans against individual antigens to train the oracles. Developability filtering approaches, though becoming standard, still do not fully reflect experimental performance.⁴⁴

Both experimental and computational approaches could benefit from ‘constraining’ the search space to only antibodies that a priori have a higher chance of becoming a therapeutic. Here, we tackled this task by setting out an approach to mapping the ‘common diversity’ shared across natural and therapeutic antibodies.

Studying this diversity poses a multi-dimensional challenge even on the sequence level. For example, different elements, such as frameworks and CDRs, all associated with different functions (frameworks as scaffolds⁴⁵ or intrinsic features (e.g., canonical classes of CDRs,⁴⁶ are included within a single variable region. Such considerations become exacerbated with other dimensions, such as structural⁴⁷ or embeddings from transformers.^48,49 To allow straight-forward analysis and interpretation of our results, we used only the IMGT-derived CDR-H3 sequence. We thus did not account for all the higher-dimensional considerations, even the clonotype, which is a limitation of the study.

Our analysis shows substantial overlap between public CDR-H3s and the therapeutic ones that is not due to statistical prevalence of shared features alone. As little as ca. 270k CDR-H3s are within two substitutions of 30% of all known therapeutic antibodies and 6% are perfect matches. The finding was possible via gross oversimplification of the antibody space to just CDR-H3, but is actually reflected by extending it to clonotypes where public CDR-H3s are observed. Therefore, combining the 270,000 CDR-H3s and the common therapeutic frameworks IGHV3–23, for instance, in the form of a display library, should statistically cover a great number of possible therapeutics.

The biological significance/cause of this observation is beyond the scope of this study (e.g., similar to the GRAB motifs.⁵⁰ It is plausible to assume that despite great theoretical diversity, antibodies in fact sample a smaller space,⁵ just by more constrained means, such as conserved structure.⁵¹ Our study appears to show that there are biases in which repertoires are constructed, which should be reflected in somatic hypermutation trajectories whenever these can be reasonably reconstructed.^52–54 We showed that there does not appear to be statistical significance between humanized or human CDR-H3s being found in the highly public dataset. Therefore, the convergences across immunoglobulins could go beyond therapeutic/natural, but between different organisms as well. Thorough investigation of organism-specific diversity, including murine as well as single-domain antibodies,³¹ is needed to reveal more actionable observations that constrain the search space for developable therapeutics.

The public-therapeutic overlap we note has an immediate actionable for ‘constraining’ the search space for either experimental or computational approaches. Naive display libraries could be biased toward only the public clones/CDR-H3s to increase the chances of molecules hitting the therapeutic, though unlikely in practice as display library providers need to guarantee ‘uniqueness of their sequences’. With poor performance of generalistic virtual antibody screening methods, constraining their search space to sequences with higher probability of being ‘therapeutics’, could fill the predictive inadequacy gap and make them usable in practice.

We hope that the findings and the data associated with the study will help in antibody engineering efforts.

Materials and methods

Sequence processing

We implemented an internal antibody numbering pipeline to avoid issues with versioning of germline databases and those caused by using different pieces of software at nucleotide and amino acid annotation level. As the human germline reference, we used the freely available germlines from VDJBase.³⁰

The software operates similarly to IgBLAST,²⁵ with a major difference of using the MMseqs2²⁹ software instead of BLAST,⁵⁵ for efficient sequence alignment.

The processing pipeline consists of following steps:

1. Prepare V, J, and C genes alignments. V gene alignments are the most time consuming, as they are performed on a full germline database as well as both forward and reverse complemented queries. From this result, we can determine a locus as well as a reading frame of the query sequence, which is used to narrow the number of possible germlines during J and C alignments. Next, the already aligned query sequence fragment (for V gene) is masked, and only the rest is taken into consideration during J alignment. This process is repeated for the C gene.

2. Join alignment results into AIRR⁵⁶ rearrangement format. Besides germline ids and aligned query sequence fragments, information such as V and J frames, region offsets as well as several useful flags (presence of stop codon, presence of frameshift or other defects that might render the immunoglobulin nonfunctional) are calculated and added to the results.

3. Number sequences with IMGT schema. We again use MMSeqs2 to calculate pairwise alignments between translated query sequence and gapped germline sequence. Resulting CIGAR string is converted to a state vector and then used directly as input to a modified version of ANARCI,²⁶ without the IMGT database and HMMER. This allows for bypassing the default HMMER alignment, which is computationally expensive.

4. The annotated alignment uses custom validation rules. This produces a collection of flags, part of which are meant to extend the alignment quality checks (absence of empty regions, conserved amino acid presence) and various technical sanity checks, such as are region offsets correct, are there no negative offsets, is CDR3 in junction.

5. Add IMGT numbering and validations to alignment data. Results in extended AIRR format.

Summarizing, there are nine extra columns in our extended AIRR format:

1. study_accession

2. run_accession

3. cluster_size

4. v_frame

5. j_frame

6. imgt_residue_mapping

7. positional_imgt_mapping

8. exc – optional, exception from IMGT numbering

9. additional_validation_flags – from step 4 described above

As one can notice, the alignment is performed twice: 1) By precisely aligning query nucleotide sequences to multiple ungapped germline sequences; and 2) By aligning translated query to gapped (using IMGT schema) amino acids sequences of corresponding genes. The second alignment, however, is much less time consuming, as we only need to align to a single, previously selected germline.

MMseqs2 performs this alignment in the following steps:

1. Frame Generation – generating all possible reading frames (forward and reverse complemented).

2. Prefiltering – the secret to MMseqs2 performance lies in this step. Using k-mers generated for query and target sequences, the software finds a matching subset and rearranges it into a diagonal. Next, top n (5 in our case) genes sharing the longest diagonal with query sequence are selected. Thanks to that, computational-heavy pairwise alignment is performed on a significantly reduced gene dataset.

3. Pairwise alignment – a fast implementation (SIMD) of Smith-Waterman algorithm is executed for query-target pairs that were selected during prefiltering.

4. Choosing best – assign the germline with the lowest E-value. It is defined as the number of expected alignments with similar scores, between a target and a set of random sequences. This way we are taking into consideration also the quality of alignment, which raw Smith-Waterman scores do not provide.

Automatic identification of B-cell receptor-containing bioprojects

The sequencing datasets are abstracted into bioprojects, which consolidate sequencing runs within one study (e.g., scientific publication) into SRA files. To identify novel bioprojects containing antibodies, we analyzed SRA files from ca. 500,000 bioprojects. For each SRA, we applied our antibody numbering pipeline, which returns an empty value for a sequence if a sequence cannot be reliably matched to a germline sequence. Because of the data volume, the pipeline was implemented in Amazon AWS Batch to perform the calculations faster, via parallelization.

We only considered SRA files where more than 10% of sequences were identified as antibodies. This is a very low threshold that was designed to increase the coverage, as SRA files with only 10% antibodies are not considered to be of good quality. A total of 287 bioprojects were identified with more than 10% antibody content. Each of the bioprojects was further manually checked that the subject of the study was BCR sequencing and were not single-cell. This brought the number to 220 high-quality bioprojects.

Metadata processing pipeline

Bioprojects are associated with biosample data that might contain a range of biological parameters relating to sequencing runs such as age and disease. The attributes are not consistent between bioprojects and thus we created a normalized attribute set where attributes available in each biosample were transferred. In the first instance, biosample metadata were automatically downloaded for our 220 bioprojects, resulting in 11,000 biosamples. Each of the biosamples was then inspected manually to fix the correspondence between attributes available in the biosample and our normalized attribute set.

Disclosure statement

PD, DC, JK, TS, BJ, SW, TG, IJ, WB, and KK are current or former employees of NaturalAntibody that commercializes computational solutions for the therapeutic antibody discovery industry.

Data availability statement

The AB-NGS dataset is available for nonprofit organizations for noncommercial use via following instructions on our site: https://naturalantibody.com/ngs/.

Additional information

Funding

The author(s) reported there is no funding associated with the work featured in this article.