To perform solid and comprehensive comparison for currently available NADTs, we employed TIgGER (6), IMPre (9), IgDiscover (5), LymAnalyzer (8) and Partis (7). Their basic information is summarized in Table 1 . As these five NADTs were developed using various programming languages, their installations are subject to various dependencies. With respect to their applications, IMPre and LymAnalyzer work on both T cell receptor (TCR) and B cell receptor (BCR) while the other three only work on BCR. All NADTs support both heavy chain (IGH) and light chain (IGK and IGL) of BCR, while IMPre and LymAnalyzer also support TRB and TRA. TIgGER and Partis only support V genes, IMPre and LymAnalyzer support V and J genes, while IgDiscover supports V, D, and J genes. Except IgDiscover and LymAnalyzer, all other NADTs underwent in silico benchmark during development. Partis developers compared their NADT with others, but no systematic third-party comparison has been performed among them. Therefore, a comprehensive and systematic comparison would benefit the field for novel allele detection using antibody repertoire datasets.

When we compared the supportive features of these NADTs, we found IMPre to be the most versatile and user-friendly NADT before considering its performance for novel allele detection ( Supplementary Table 1 ). To gain more insights into these NADTs, we evaluated their performance with both simulated and real-world Ig-seq datasets ( Figure 1 ). The benchmark result was then summarized and translated into knowledge-based filtration criteria used to obtain credible NACs from collected bulk sequencing dataset.

Generating in silico Ig-seq datasets is a challenging task. An ideal Ig-seq simulating tool should reflect the preferential gene usage, junctional nucleotide insertion and deletion, phylogenetic clonal structure, various allele ratio, and the base errors intrinsic to PCR amplification and next-generation sequencing (NGS). Although several repertoire simulation tools exist (7, 10–13), none of them incorporate the full features of Ig-seq dataset mentioned above. Therefore, we built IMPlAntS (Integrated and Modular Pipeline for Antibody Repertoire Simulation) and it enables both one-stop repertoire simulation and modular calls for adaption to customized pipelines.

Briefly, IMPlAntS consists of three consecutive steps, i) generation of independent V(D)J rearrangements; ii) generation of SHM with phylogenetic structure within clones; and iii) generation of NGS reads incorporating base errors ( Figure 2 ). These steps can be implemented individually or collectively using the corresponding scripts.

In the first step, a series of key parameters can be specified in the configuration files. These parameters include V(D)J gene usage, allele ratio, the distribution of insertion and deletion length, and the percentage of productive rearrangements. In the second step, we generated SHMs in rearranged sequences in a way similar to that reported by Yermanos et al. (13) to create the phylogenetic sequences as in the real repertoire. The resultant repertoire with SHMs comes from several iterations of introducing SHMs to the selected sequences based on the positional mutability and substitutability models. These two models, together with the parameters involved in the first step, derive from our previous large-scale study (14). Finally, we employed a popular NGS simulation tool, ART, to produce NGS reads (15). More details for IMPlAnts can be found in Materials and Methods section.

With this pipeline, we generated four datasets: DEXPR, DSNP, DALLELE, and DSHM ( Table 2 ). Noteworthy is that only DSHM was generated with all three steps mentioned above. In contrast, the other three datasets were generated with only the first and the final step as they contain no SHMs. Each of the four datasets was comprised of 20 repertoires, except for DSHM (n=10). The constituent repertoires within each dataset contained variation only in the studied variable. Except for DSHM, each dataset contained four groups (two groups for DSHM). While each group is represented by five repertoire replicates and has a distinct level as to the studied variables. Other variables were set identically among groups within each dataset and to a level theoretically most favorable to novel allele detection. For each repertoire, we generated 1 million reads to avoid the read number limitation mentioned in the IgDiscover manual (at least 750,000 was recommended). The only exception was with DALLELE, in which repertoires in different groups had varying numbers of reads to make the novel alleles represented by the same number of reads. Lastly, we artificially created “novel” alleles by random selection of the positions and SNPs in germline sequences. The resultant “novel” alleles together with known ones then served as the initial germline database for NADTs’ benchmarking (Materials and Methods).

To compare the sensitivity and specificity of the 5 NADTs in detecting novel alleles (allele level) (LymAnalyzer was excluded as it reported only SNPs) as well as SNPs (SNP level) (Materials and Methods), we used our in silico simulated datasets ( Table 3 ). As expected, lower gene or allele expression and more SNPs or SHMs hampered both sensitivities and specificities for at least one NADT in the detection of novel alleles and SNPs in general ( Table 3 , Figure 3 ). We found TIgGER to work well with respect to both sensitivities and specificities with DEXPR and DSNP, although it did not identify alleles in DALLELE ( Table 3 ). IMPre, though exhibiting lower sensitivities and specificities, identified novel alleles in the datasets with all four variables. IgDiscover manifested very good specificities although it identified fewer alleles than TIgGER. The performance of Partis was less optimal in DSNP than that of TIgGER but excelled in DALLELE and higher SHM datasets. As LymAnalyzer only reports SNPs, it was excluded from allele level comparisons. However, it also showed high sensitivities in all situations in SNP level although the sensitivities were less ideal. The performance of other NADTs was similar in SNP level to that of the allele level.

Taken together, TIgGER, IgDiscover, and Partis showed comparably high specificities and therefore the alleles identified were more reliable. IMPre and LymAnalyzer provided more allele candidates, but none of the NADTs performed well in all situations. However, each of these datasets was simulated with only one variable with a particular quantity to evaluate the effect of these quantitative measures on the performance of NADTs whereas real-world repertoires always consist of combinations of all variables in multiple quantitative measures.

An ideal situation to test NADTs is to genotype all the V alleles in a genome and then compare them with NADTs’ predictions. However, given the high similarities of V alleles and other interspersed tandem sequences among them, sequencing this peculiar region of the genome alone is a challenging task (16). Therefore, we took an alternative approach by acquiring germline V allele sequences from single-cell repertoire sequencing of naïve B cells (scRep-seq) and then conducted novel allele identification on the bulk Ig-seq datasets from the same donor. The naïve state of antibody sequences and the super-high depth of the scRep-seq data ensured the accuracy of acquired germline sequences. Thus, this evaluation represents the real-world situation.

Three healthy donors were included in real-world dataset benchmark. Briefly, we isolated naïve B cells from the peripheral blood of these 3 donors with specific cell surface markers. Then sequencing libraries were constructed according to the protocols of Chromium Single Cell Human BCR Amplification Kit (10X Genomics) and sequenced on an Illumina platform (Materials and Methods). Each single cell sequencing sample contains around 85 million to 91 million paired end reads. With the single naïve B cell sequencing dataset, we identified 4 unique NACs from 3 donors using a customized pipeline ( Table 4 , Materials and Methods). All identified NACs are minor alleles of the involved genes, with expression ratios to the major ones ranging from 0.19 to 0.89. Moreover, each of them only harbors one SNP compared to their nearest known alleles.

We then applied NADTs to the bulk sequencing datasets and compared their novel allele predictions. TIgGER identified three and IgDiscover identified one out of the four NACs ( Table 4 ) while IMPre and Partis missed all of them. Although LymAnalyzer identified two positive SNPs from two NACs, it also falsely predicted 14 and 6 SNPs in these two alleles, respectively. In addition, we found two possible novel germline sequences that harbor a considerable number of mismatches with their nearest known germline sequences ( Supplementary Table 2 ). Notably, the novel germline sequence nearest to IGHV1-NL1*01 was identified in 2 of 3 enrolled donors.

Because of the limited number of bona fide novel alleles (4 NACs in Table 4 were deemed as bona fide novel alleles here), this evaluation was less comprehensive. We thus exploited the genuine Ig-seq dataset (donor1 and donor3) in another way. As many germline sequences were known through the single naïve B cell sequencing dataset, we artificially generated “novel” alleles as those in the simulated dataset mentioned above and evaluated these NADTs in the same way (Materials and Methods). However, because there were not enough genes expressing at around 0.01%, we did not generate novel alleles at this level. Moreover, the allele ratio is hard to precisely infer even with the single naïve B cell sequencing dataset and was thus left unstudied. We denoted the genuine dataset with different initial databases as GD-EXPR, GD-SNP, and GD-SHM.

The genuine dataset-based benchmark result exhibited a similar performance spectrum as that based on the simulated dataset ( Table 5 ). These similarities included, i) each of the three studied factors was found to be influential for at least one NADT for novel allele detection ( Supplementary Figure 1 ), ii) TIgGER and IgDiscover were superior to IMPre and Partis in both sensitivity and specificity for detecting novel alleles with multiple SNPs (i. e. 3, 5, and 7) ( Supplementary Figure 2 ), iii) the SNP-level performance spectrum in general resembled that of allele level, iv) IMPre and Partis presented higher sensitivity and specificity for identifying SNPs in DSNP than for alleles, and v) LymAnalyzer remained the most sensitive but least specific NADT in identifying SNPs.

However, several variations were also remarkable and they included, i) for genuine datasets, TIgGER and IgDiscover performed better in identifying novel alleles expressed at a low level (i. e. ~0.1%) than for simulated dataset; both were thus superior to IMPre and Partis, ii) although Partis remained excellent in overcoming SHM noise, it was outperformed by IgDiscover ( Supplementary Figure 2 ), which exhibited a surprisingly high sensitivity of 0.80 and specificity of 1.00 at both SNP and allele levels, iii) IgDiscover manifested significantly higher sensitivities and specificities than TIgGER in three datasets, and iv) LymAnzlyzer displayed low and negligible specificities.

To seek the underlying reasons accounting for these discrepancies, we assessed the output of these NADTs as well as the properties of each input dataset. We found that the inferior performance of IgDiscover and TIgGER on DEXPR in the simulated dataset was caused by low sequence identities to the germline ( Supplementary Figure 3 ), which was caused by sequencing errors simulated by ART, the NGS simulator. Similarly, the failure in DSHM for IgDiscover could also be attributed to the paucity of unmutated sequences as a consequence of the simulation of SHMs and NGS errors ( Supplementary Figure 4 ). In contrast, the success in GD-SHM for IgDiscover indicated that a number of unmutated sequences also exist in IgG dataset. We therefore determined the frequency of such sequences for each simulated novel allele and found that it ranged from 0.07% to 1.11%, with a median of 0.31% ( Supplementary Figure 4 ), which agrees with reported values of previous studies (17, 18). Interestingly, TIgGER failed to detect any novel alleles from both simulated and genuine datasets containing SHMs as its algorithm is expected to be more robust to datasets with SHMs. Our in-depth analyses showed that TIgGER failed to identify novel alleles for DSHM and GD-SHM for different reasons. For DSHM, the simulated SHMs created an overly-diversified repertoire, in which plural sequences for each novel allele were too rare to pass the threshold min_seqs. Whereas in GD-SHM, the diversity of sequences perfectly matching novel alleles failed to meet the default threshold j_max. In addition, we noted a remarkable difference in the diversity filtration criterion between TIgGER and IgDiscover: TIgGER employs a quantitative filtration (j_max) whereas IgDiscovers uses a qualitative filtration (CDR3_exact). When considering only the diversity criterion, TIgGER is stricter than IgDiscover, and this explains the compromised performance of TIgGER. Finally, the lower specificity of LymAnalyzer in the genuine dataset may result from the non-independent mutation events in a genuine dataset that tends to be interpreted as SNPs according to its algorithm.

Together, we concluded that i) TIgGER and IgDiscover outperform all other NADTs considering both sensitivity and specificity in most situations, ii) Partis is characterized by remarkable robustness in overcoming the challenge imposed by SHMs, iii) IMPre is outstanding in detecting minor alleles, and iv) LymAnalyzer is sensitive at the cost of specificity.

With the knowledge obtained above, we designed a scheme to identify reliable NACs using 4 NADTs (excluding LymAnalyzer) from bulk Ig-seq dataset. As intrinsic features (i.e. expression level, allele ratio, and number of SNPs to the nearest allele) of novel alleles were unknown, we took into account the overall performance of each NADT summarized above and gave more credit to TIgGER and IgDiscover. We classified all Ig-seq datasets into two groups with regard to the SHM richness according to the isotypes (Materials and Methods). For IgM datasets, NACs found by at least 2 NADTs with at least one being either TIgGER or IgDiscover were retained. For datasets in which SHMs were expected to be enriched, only NACs called by two out of three NADTs, namely TIgGER, IgDiscover and Partis, were retained.

We then explored the efficiency of this scheme in identifying NACs from a total number of 424 Ig-seq datasets either generated in-house or from the public resource (Materials and Methods). The selected datasets stemmed from 382 donors and were all obtained from RNA samples amplified with RACE (rapid amplification of cDNA ends) protocols. Detailed metadata for these datasets are outlined in Supplementary Table 3 . According to the dataset classification criteria (Materials and Methods), we obtained 336 (79.2%) SHM-rich datasets (enriched for IgG sequences) and 88 (20.8%) SHM-sparse datasets (enriched for IgM sequences) ( Figure 4A ). Despite the lower fraction in overall datasets, IgM datasets contain more reads than IgG datasets ( Figure 4B ). Applying the selected 4 NADTs to these datasets, we found clear differences between the four NADTs in both the number of samples identified with NACs and the number of unique NACs ( Table 6 ). IMPre discerned NACs for 71.0% of the datasets, whereas the other three NADTs found NACs for only 16.3% to 18.2% datasets. Moreover, the other three NADTs reported NACs for a sharply lower (over 10-fold) percentage of SHM-rich datasets than SHM-sparse datasets, which was likely due to more SHMs and low number of input reads that had reduced the confidence for these NADTs to make novel calls. In contrast, IMPre reported NACs for a large fraction of IgG datasets (63.4%) and also more NACs overall for individual samples ( Table 6 and Figure 4C ), which likely reflects its higher sensitivity to those underrepresented sequences ( Table 3 ). However, the genuine sensitivity and specificity for the NADTs were elusive through these bulk sequencing datasets, for which we have no access to the genotype information.

Applying this scheme to 424 Ig-seq datasets, we identified 23 and 2 reliable NACs from SHM-sparse and SHM-rich group, respectively ( Supplementary Table 4 ). One NAC, IGHV3-33*01_G72C, was identified in both groups. Three of the 24 unique NACs were found to harbor more than one SNPs to their corresponding nearest alleles, while eleven were found in more than one donor ( Figure 5 ). The most frequent NAC was found in 29 donors. Notably, 17 of the 24 NACs can also be identified from public databases or independent reports in the literature (Materials and Methods) ( Supplementary Table 4 ), which demonstrated the high efficiency of our scheme. To enlarge the knowledge database of NACs, we also included 263 multiplex datasets we collected in a previous study into our analysis (14). These latter datasets were derived from 71 donors and consisted of 186 SHM-rich datasets and 77 SHM-sparse datasets. Considering the degenerate primers designed against framework region 1 (FR1) of V genes, we considered only the sequence downstream of FR1 for each NAC for these multiplex datasets. Applying the same scheme to these datasets as to RACE datasets, we identified in total 22 NACs ( Supplementary Table 5 ) and found that 21 of them were from SHM-sparse datasets. Eleven of the 22 NACs were cross-validated in previous public or published resources. Combining the two NAC sets, we identified 43 unique NAC sequences from a total number of 687 Ig-seq datasets (3 NACs were found in both RACE and multiplex dataset).

We then characterized these 43 NACs and found that all NACs derive from core V genes we defined in a previous study according to their prevalence in antibody repertoires (14) ( Figure 6A ). This result further suggested that gene usage is critical in novel allele identification through Ig-seq dataset. Furthermore, the number of NACs did not correlate with the number of known alleles for a typical gene (Pearson correlation coefficient: 0.43) ( Figure 6B ). However, IGHV1-69, the gene with the second largest known polymorphisms in IMGT, was found with up to 10 additional NACs. Since germline V gene polymorphisms have been implicated in immune response capability (1–3), these NACs will facilitate the elucidation of the role of germline variants in disease susceptibility. Finally, we classified all identified SNPs (n=75) into two categories, replacement (R) SNPs and silent (S) SNPs, according to the variation of encoded amino acids. Overall, the R/S ratio for these SNPs was around 2 (1.88) ( Figure 6C ). Nevertheless, the R/S ratio was larger for complementarity-determining region (CDR) SNPs (2.78) than framework region (FR) SNPs (1.41), which indicated a varied selection pressure between FRs and CDRs.

A total of 28 samples from peripheral blood, tumor and normal tissues, and bone marrow were collected. Of these, 7 peripheral blood samples were derived from healthy individuals (without recent infection events), 6 peripheral blood samples were from hepatitis B virus-infected patients, 1 bone marrow sample and 2 peripheral blood samples were from graft-versus-host disease (GvHD) patients, 4 peripheral blood samples, 1 normal intestine sample and 2 intestine tumor samples were from colorectal cancer (CRC) patients, 2 peripheral blood samples were from individuals involved in traffic accidents, and 3 peripheral blood samples were from patients with adolescent idiopathic scoliosis, sore throat, and chronic pharyngitis, respectively. Peripheral blood mononuclear cells (PBMCs) and bone marrow mononuclear cells were isolated using Ficoll (TBD Science) density-gradient centrifugation. The tissues were cut into small pieces and grind with liquid nitrogen. These experiments were handled under the guidelines of the Ethics Committee of Southern Medical University. For human naïve B cells isolation, PBMCs were counted and washed with DPBS supplemented with 1% bovine serum albumin (BSA), and then were stained with a cocktail of fluorescent conjugated antibodies, including ECD-CD19 (Beckman Coulter, A07770), FITC-IgD (Beckman Coulter, B30652), APC-CD27 (BD Bioscience, 561400), and 7-AAD (BD Bioscience, 559925). Human naïve B cells (CD19+IgD+CD27-7-AAD-) were sorted using a cell sorter (MoFlo XDP, Beckman Coulter) and collected for single-cell V(D)J sequencing.

RNA purification was carried out using the RNeasy Mini Kit (Qiagen, 74106) according to the manufacturer’s instructions. Total RNA was used as a template to synthesize cDNA with a SMARTer RACE (Rapid Amplification of cDNA Ends) cDNA Amplification Kit (Clontech, 634928) according to the manufacturer’s protocol. Heavy chain variable regions were amplified using 1 μl of RT reaction product and 10 pmol of each primer in a 50 μl total reaction volume (KAPA HiFi HotStart ReadyMix, Roche) using the following thermal cycling program: 95°C for 3 min; 30 cycles of 98°C for 20 s, 60°C for 15 s, and 72°C for 15 s; 72°C for 5 min. PCR products were purified using the Nucleospin Gel & PCR Clean-up kit (Macherey-Nagel, 704609.25) and subjected to library preparation using VAHTS Universal DNA Library Prep Kit (Vazyme, ND607-01). Libraries were quantified by capillary electrophoresis (Bio-Fragment analyzer, Bioptic). After quantification, libraries were pooled and sequenced on an Illumina platform (MiSeq PE300). All primers are listed in Supplementary Table 6 .

The concentration of the single cell suspension was counted and adjusted to 1000 cells/μl. The single cell suspensions were loaded onto the Chromium Controller microfluidics device (10X Genomics) and processed using Chromium Next GEM Single Cell 5’ Kits v2 according to manufacturer’s protocol. The remaining procedures, including library construction, were performed according to the protocols of the Chromium Single Cell Human BCR Amplification Kit (10X Genomics). Following library construction, the BCR libraries were sequenced on an Illumina platform (NovaSeq 6000) using 2×150bp kit.

IMPlAntS (Integrated and Modular Pipeline for Antibody Repertoire Simulation) was developed to as much as possible mimic real-world antibody repertoires and meet the requirements (i.e. minor allele frequency control and NGS data simulation) in this study.

As mentioned in Results, IMPlAntS consists of three consecutive steps, i) generation of independent V(D)J rearrangements; ii) generation of BCRs with SHMs of proper phylogenetic structure within clones; and iii) generation of NGS reads incorporating base errors ( Figure 2 ). These steps can be implemented individually or collectively using the corresponding scripts hosted on github (https://github.com/Xiujia-Yang/IMPlAntS).

In the first step, a customizable number of independent rearranged sequences are in silico simulated by considering two major features of the real-world rearrangement repertoire: preferential gene usage and junctional nucleotide modification (P and N nucleotide insertions and deletions). To investigate the influence of allelic diversity on novel allele identification, we equipped IMPlAntS with the ability to simulate alleles of a certain gene with varied ratios (only two alleles are supported), which can be customized by modifying the gene usage configuration file. Moreover, IMPlAntS also allows simulation of nonproductive rearrangements and their percentages in antibody repertoire can be specified by users for specific aims. Notably, the four simulated datasets in this study include only productive rearrangements.

The second step can be further divided into two stages: generation of clonally related sequences with proper phylogenetic structure and various numbers for each sequence. Clonally related sequences are created by a certain number of iterations (to mimic the affinity maturation of real-world antibody sequence) where SHMs are induced for randomly selected sequences across the variable region based on the positional mutability and substitution models similar to Yermanos et al. (13). In each iteration, a fraction of sequences in the current sequence pool are randomly selected for SHM simulation and new sequences with simulated SHMs will be added into the current sequence pool that will be subjected to random selection in the next iteration. Independent rearranged sequences serve as the input in the first iteration. Because the positional mutability model stores mutation probabilities for different positions observed in end repertoires (repertoires containing sequences have already undergone multiple rounds of maturations), a parameter named ‘—mut_ability_fold’ (less than 1) is introduced here to prevent the generation of hyper-mutated sequences after a number of iterations. Iterations above produce nonredundant clonally related sequences. Then selective sequences will be populated according to the power law (20) to mimic the clonal expansion of B cells with a various number of replicates. The key parameters in this step, including the number of iterations, the maximum number of sequences, the alpha value of the power law, and the largest size of sequences, are customizable. ART is employed in the last step to produce NGS data with Illumina MiSeq system settings.

For the above steps, parameters of gene usage, junctional modification, positional mutability and substitution models, were obtained from a population-level antibody repertoire study (14) and are set as defaults of IMPlAntS. Gene usage is calculated as the percentage of clones (sets of sequences sharing the same V and J gene and CDR3 nucleotide sequence) in a repertoire recombined from a certain gene. In this study, V, D and J gene usages are taken from normalized medians of gene usages from 2152 antibody repertoires of 582 donors. Junctional modification parameters consist of 10 entities (i.e. V3D, V3P, N1, D5D, D5P, D3D, D3P, N2, J5D and J5P (D, deleted nucleotide; P, palindromic nucleotide; N, nontemplated nucleotide), as demonstrated also in Figure 2 ). The probabilities of modification lengths for each of these entities are derived from the observation of a combination of 2152 antibody repertoires of 582 donors. The positional mutability and substitution models were obtained from IgG repertoires of PBMC from 353 healthy donors. All parameters above can be found on the github and are set as defaults by IMPlAntS. Supplementary Figure 5 and Supplementary Figure 3B show the approximation of the real-world repertoire for repertoires in the four simulated datasets in this study.

In this study, the germline reference sequences for V, D, and J genes were obtained from IMGT GENE-DB and provided as Supplementary Table 7 . The artificially “novel” alleles for V genes were created for both simulated dataset and real Ig-seq dataset. Only germline reference sequences used in the simulation were extracted to serve as the initial reference sequences for the simulated dataset. The set of alleles subject to the artificial SNP generation for each dataset was selected according to the criteria defined as Table 2 . We randomly created SNPs in the sequence of selected alleles. These artificial SNPs were set to locate in the first 280 bp of V genes at the 5’ ends to avoid the possible failure in novel allele detection caused by junctional modification. A pitfall here is that there exists a possibility that the rearranged sequences fail to be best aligned against the artificially novel sequences, and this brings challenges in the evaluation of novel allele identification for NADTs. Therefore, we performed pairwise alignment between customized reference sequences and the germline sequences contained in each dataset and removed those unaltered allele sequences that were found to be more similar to the germline sequences than the “novel” allele sequences. The novel alleles identified by NADTs were in fact the real-world germline sequences, while “novel” is just a concept relative to the altered germline reference sequences.

The pair-end simulated dataset and bulk sequencing dataset were firstly assembled using PEAR (v0.9.6). The successfully assembled sequences were then taken as the input for IgDiscover and LymAnalyzer. As TIgGER can only accept a formatted database of well-annotated sequences as input, we further annotated and formatted the assembled sequences with IgBLAST (v2.8.0+) and Change-O toolkits (v0.4.4), respectively (IgBLAST was selected for its excellent performance (20) and easy output format conversion through Change-O toolkits). For IMPre and Partis, the input assembled sequences were corrected in a forward orientation at first. The script employed by IMPre (‘IMPre.pl’) was modified to enable germline reference customization. The revised script, ‘IMPre_revised.pl’, can be found on the github (https://github.com/Xiujia-Yang/IMPlAntS). All parameters used by the five NADTs were set in default or as suggested. We provided the detailed commandline parameters as below,

Sensitivity is defined as the proportion of true positives that are correctly identified among all true positives, whereas specificity is defined as the proportion of true positives among all the identified positives. For individual SNPs (SNP level), a hit is considered as a true positive only when its nearest allele (same as the allele selected for artificial SNP generation), loci and nucleotide variant are correct at the same time. For individual sequences (allele level), a hit is considered a true positive only when it covers all the genuine SNPs and contains no mismatches with the genuine novel sequence in all other reported loci. A schematic diagram is provided here to demonstrate the cases of true positive and false positive in identifying individual sequences for novel alleles ( Scheme 1 ).

Cell Ranger (v3.1.0) was exploited to preprocess the raw single naïve B cell sequencing dataset. Contig assembly, annotation, and clonotype analysis were performed using “cellranger vdj” with the Cell Ranger V(D)J compatible reference (refdata-cellranger-vdj-GRCh38-alts-ensembl-3.1.0). Then the assembled contig sequences (“all_contig.fasta”) of the two replicates for each donor were pooled and then annotated using IgBLAST (v2.8.0+) with the germline references obtained from IMGT/GENE-DB (refer to above). Afterwards, the V segment (or allele) sequence was extracted from each annotated sequence and then each unique V segment sequences was counted. It is worth mentioning here that those short V segment sequences were merged into the longer ones provided they are with the same V allele annotation as the longer ones and were included in them. The counts for the short V segment sequences were also added to the longer ones. We discarded those with a length less than 290 bp or with a count less than 10 and determined the most frequent V segment sequence for each gene as the most confidential germline sequence for a gene. Apart from that, we also retained the second most frequent V segment sequence for a gene provided that its abundance was at least one tenth of that of the most frequent one (21).

All enrolled Ig-seq datasets (i.e. 424 RACE datasets and 263 multiplex datasets mentioned in the Results section) were analyzed using MiXCR (v3.0.7) per the method in our previous study (14). After clonotype assembly, a constant gene will be assigned for each clone if antibody sequences from this clone cover constant region. The isotype (i.e. IgM, IgD, IgG, IgA, and IgE) was extracted for each clone and the clone-level isotype frequency was calculated for each dataset. IgM and IgD are deemed as SHM-sparse isotypes while IgG, IgA and IgE are deemed as SHM-enrich isotypes (22). Datasets will be classified as IgM datasets if they contain more SHM-sparse isotypes than SHM-enrich isotypes, otherwise IgG datasets. Constant genes were required to be assigned for more than a half number of clones in each dataset. All 687 Ig-seq datasets we enrolled in this study met this requirement.

To double-check NACs we identified through NADTs, we collected antibody heavy chain V allele sequences from five public databases (IMGT/GENE-DB, http://www.imgt.org/genedb/; IgPdb, https://cgi.cse.unsw.edu.au/~ihmmune/IgPdb/information.php; VBASE2 (23), http://www.vbase2.org/; Lym1K (24), http://maths.nuigalway.ie/biocluster/database/; OGRDB (4), https://ogrdb.airr-community.org/) and nine independent reports (5, 6, 19, 25–30) and compared them with identified NACs. Before the sequence comparison, degenerate bases or N nucleotides in collected allele sequences were substituted with ‘A’, ‘C’, ‘G’, or ‘T’, accordingly. NACs whose sequences were identical to any of the sequences from a source were considered cross-validated NACs. As the set of V allele sequences used as germline reference to identify novel alleles is not as complete as those in the later release of IMGT/GENE-DB, several NAC sequences were included in the later release of IMGT/GENE-DB and thus were also cross-validated in it. All collected V allele sequences are outlined in Supplementary Table 8 .