Briefings in Bioinformatics, 2022, 23(5), 1–11 https://doi.org/10.1093/bib/bbac349 Problem Solving Protocol IgIDivA: immunoglobulin intraclonal diversification analysis Laura Zaragoza-Infante , Valentin Junet , Nikos Pechlivanis, Styliani-Christina Fragkouli , Serovpe Amprachamian, Triantafyllia Koletsa, Anastasia Chatzidimitriou, Maria Papaioannou, Kostas Stamatopoulos, Andreas Agathangelidis and Fotis Psomopoulos Corresponding author: Fotis E. Psomopoulos, Institute of Applied Biosciences (INAB), Centre for Research and Technology Hellas (CERTH), Thessaloniki, Greece 6th km Charilaou - Thermi, 6th km Charilaou - Thermi, Thessaloniki 570 01, Greece. Tel.: +30-2310-498-478; Fax: +30 2310 498 270; E-mail: fpsom@certh.gr Abstract Intraclonal diversification (ID) within the immunoglobulin (IG) genes expressed by B cell clones arises due to ongoing somatic hypermutation (SHM) in a context of continuous interactions with antigen(s). Defining the nature and order of appearance of SHMs in the IG genes can assist in improved understanding of the ID process, shedding light into the ontogeny and evolution of B cell clones in health and disease. Such endeavor is empowered thanks to the introduction of high-throughput sequencing in the study of IG gene repertoires. However, few existing tools allow the identification, quantification and characterization of SHMs related to ID, all of which have limitations in their analysis, highlighting the need for developing a purpose-built tool for the comprehensive analysis of the ID process. In this work, we present the immunoglobulin intraclonal diversification analysis (IgIDivA) tool, a novel methodology for the in-depth qualitative and quantitative analysis of the ID process from high-throughput sequencing data. IgIDivA identifies and characterizes SHMs that occur within the variable domain of the rearranged IG genes and studies in detail the connections between identified SHMs, establishing mutational pathways. Moreover, it combines established and new graph-based metrics for the objective determination of ID level, combined with statistical analysis for the comparison of ID level features for different groups of samples. Of importance, IgIDivA also provides detailed visualizations of ID through the generation of purpose-built graph networks. Beyond the method design, IgIDivA has been also implemented as an R Shiny web application. IgIDivA is freely available at https://bio.tools/igidiva Keywords: intraclonal diversification, B cell receptor immunoglobulin, high-throughput sequencing, graph networks, graph metrics Introduction with diverse BcR IG is warranted to endow the host with immune Recognizing antigens is B cells’ ‘raison d’être’. This is accom- competence. The extraordinary diversity of the human BcR IG plished through the immunoglobulin (IG), which forms the part of repertoire relies largely on V(D)J recombination, a combinatorial the B cell receptor (BcR) that mediates antigen recognition [1, 2]. association of distinct IG heavy and light chain variable (V), Considering the enormous antigen diversity in nature, it is obvious diversity (D; for heavy chains only) and joining (J) genes occurring that a correspondingly vast repertoire of antigen-specific B cells in developing B cells. Moreover, the variable regions of the IG Laura Zaragoza-Infante is a PhD student of immunogenetics at the Institute of Applied Biosciences (INAB), at the Centre for Research and Technology Hellas (CERTH) and is associated with the Faculty of Medicine at the Aristotle University of Thessaloniki. Her research interests include the immunogenetic study of B cell lymphomas and the application of bioinformatics in immunogenetics. Dr Valentin Junet was, at the time of this study, employed by Anaxomics Biotech S.L. and was a PhD candidate at the Autonomous University of Barcelona. His research interests include bioinformatics and the application of machine learning to biomedical data. Nikos Pechlivanis works as a bioinformatics engineer at INAB|CERTH. He is pursuing a PhD in the field of computation biology and bioinformatics. His research focuses on studying machine learning algorithms and computational techniques with respect to their corresponding applications on genomic data. Styliani-Christina Fragkouli obtained her BSc in physics from the Aristotle University of Thessaloniki in 2019 and is working as research assistant at INAB|CERTH. Her interests include machine learning and bioinformatics. Serovpe Amprachamian is currently an undergraduate student of electrical and computer engineering at the Aristotle University of Thessaloniki. His research interests include software development, machine learning, deep learning and their for data analysis and mining in life and medical sciences. Dr Triantafyllia Koletsa is an associate professor of pathology at the School of Medicine, Aristotle University of Thessaloniki. Her research interests focus on tumor biomarkers, the immunogenicity of a tumor and its microenvironment. Dr Anastasia Chatzidimitriou is a senior researcher at INAB|CERTH. Her research activities focus on the field of immunobiology and genomics applying NGS methodologies for the detection of new biomarkers in cancer. Dr Maria Papaioannou is a professor of hematology at the Aristotle University of Thessaloniki and is a collaborating researcher at INAB|CERTH. Beyond her clinical work, her research interest focuses on immunogenetics of hematologic diseases (ALL, HIV Lymphomas, CLL and PNH). Dr Kostas Stamatopoulos is a specialist in hematology, director of INAB|CERTH. His research interests include precision medicine in cancer, immunobiology of lymphoid malignancies and patient empowerment in care. Dr Andreas Agathangelidis is an assistant professor of molecular genetics at the Department of Biology of the National and Kapodistrian University of Athens. His research interests concern the characterization of genes, signaling pathways and cellular mechanisms associated with the ontogeny and pathogenesis of lymphoid malignancies. Dr Fotis Psomopoulos is a researcher at INAB|CERTH. His research interests lie at the intersection of bioinformatics and machine learning, primarily working on the design and implementation of novel algorithms for knowledge extraction from large datasets in life sciences. Received: April 20, 2022. Revised: July 25, 2022. Accepted: July 26, 2022 © The Author(s) 2022. Published by Oxford University Press. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. 2 | Zaragoza-Infante et al. heavy and light chains, representing the antigen binding sites, identity between the different sequences, often differing in only comprise four framework regions (FR) and three hypervariable a few mutations, and given that traditional phylogenetic analysis complementarity determining regions [3]. Successful completion are not suitable for these cases [41, 42]; (iii) they do not provide any of V(D)J recombination leads to the expression of functional BcR analysis of graph metrics that could greatly assist in the quantifi- IG of both IgM and IgD isotypes on the surface of naive B cells, cation of ID levels and, subsequently, in performing comparisons rendering them competent to effectively recognize antigens [2, 4]. between samples or groups of samples; (iv) some of them are not Once this happens, B cells mature further in specific microenvi- publicly available; (v) they are restricted to the analysis of the ronments within the secondary lymphoid organs, called germinal IG heavy chain gene repertoire; (vi) they do not provide a user- centers, through two distinct molecular processes: somatic hyper- friendly interface that would enable their application without mutation (SHM) and class-switch recombination (CSR) [1, 5, 6]. prior programming knowledge and (vii) they are visualization- Both processes are catalyzed by the enzyme activation-induced only tools. deaminase (AID) [7, 8]. SHM mostly entails the introduction of Here, we present the immunoglobulin intraclonal diversi- point mutations in the IG variable domain. These mutations can fication analysis (IgIDivA) tool, a purpose-built tool for the alter the affinity of the antibody for its cognate antigen, with detailed assessment of the ID process through the analysis mutations that lead to an increase in affinity being promoted [9]. of high-throughput NGS data. IgIDivA provides a detailed The introduction of mutations within rearranged genes occurs at characterization of all mutations occurring in the context of rates of 10−5 –10−3 mutations per base pair per generation, 106 - the ID process as well as their connections. Subsequently, fold higher than spontaneous mutations occurring elsewhere in connections are used for ‘building’ mutational pathways, while the genome [10, 11]. On the other hand, CSR is responsible for a series of graph metrics are calculated toward a multifaceted the replacement of the IG heavy chain constant gene from IGH- characterization of ID. Finally, statistical analysis is provided M/IGHD to IGHG or IGHE or IGHA, switching antibody production for comparisons of different features of ID among samples or from IgM/IgD to a different class, such as IgG, IgE or IgA, without groups of samples. Overall, IgIDivA is a user-friendly pipeline, altering the antigen specificity of the antibody [12]. also available as an R Shiny app, freely available at https://bio. The aforementioned BcR IG diversity of the immune system tools/igidiva. in a healthy individual is ref lected in the polyclonality of the respective repertoire. Human diseases implicating B cells may vary in terms of BcR IG gene repertoire diversity: some are poly- Materials and methods clonal (for instance, systemic lupus erythematosus is associated Data pre-processing with intense polyclonal B cell activation) [13], whereas others The main pre-processing workf low, depicted in Figure 1, consists are characterized by oligoclonal (e.g. rheumatoid arthritis and of (i) quality filtering and synthesis of raw reads (in the case multiple sclerosis) [14, 15] or even monoclonal B cell expansions of paired-end sequencing) in order to obtain high-quality, full- (B lymphoid malignancies) [16]. length BcR IG sequences; (ii) data annotation with IMGT/HighV- An additional level of complexity may arise when focusing on QUEST [43] and (iii) meta-data analysis with tripr (T cell receptor/ specific, relevant B cell clonal expansions. In such a context, BcR immunoglobulin profiler in R) for the grouping of sequences into IG repertoire diversity may increase through a process known as clonotypes and their alignment [44]. IgIDivA requires as input two intraclonal diversification (ID), which entails the introduction of tripr output files per sample. First, the ‘clonotypes computation’ ongoing SHMs due to continuous antigenic pressure [6, 17, 18]. file, i.e. grouping of BcR IG gene rearrangement sequences into Studies of the ID process have provided valuable insight into the clonotypes. Second, the ‘grouped alignment’ file, which concerns ontogeny and evolution of B cell clones in health and disease the alignment of all BcR IG gene rearrangement sequences against [19–23]. However, most relevant studies were performed using the respective germline V genes at the nucleotide (nt) level. Each low-throughput, Sanger sequencing; hence, they were inherently unique sequence or group of identical sequences is considered as limited with regard to analytical depth and breadth [22–28]. This a nucleotide variant (nt var). The combination of these two files limitation was recently surpassed due to the advent of next- allows for the subsequent selection of specific clonotypes. generation sequencing (NGS), allowing a deeper and, thus, more accurate capture of the diversity of the BcR IG gene repertoire, Processing of the input both at the clonal and the subclonal levels, the latter being directly The IgIDivA approach is based on the assumption that, during the associated with ID [29–31]. However, in order to fully under- process of ID, SHMs are acquired sequentially due to continuous stand the complex immunogenetic ‘mechanics’ of ID, purpose- interactions with antigen(s). Hence, identifying the ‘timeline’ of built bioinformatic tools are required. SHM is a critical aspect of the process. Currently, several different bioinformatic approaches exist for The first step in our approach concerns clonotype selection: the analysis and visualization of SHM within the BcR IG gene depending on the context, the simplest option would be to focus rearrangement sequences and their classification in the context on the most expanded (i.e. the most frequent) clonotype. However, of ID, such as ClonalTREE [32], GCTree [33], GLaMST [34], IgTree any number of clonotypes can be selected per sample and each [35], MTree [36], ViCloD [37], Alakazam [38, 39] and AncesTree [40]. clonotype will be processed independently. While ClonalTREE was developed for a different purpose (bacterial Subsequently, the ‘grouped alignment’ file should be filtered in composition evolution), the rest are dedicated to the study of order to limit the analysis to the selected clonotype(s). BcR IG repertoires. That notwithstanding, most existing solutions display one or more of the following limitations: (i) the analytical Network generation process in all of them is based on the inference of mutational Network generation comprises the identification of all different variants, relying on statistical analysis for network generation, nt vars, their SHM patterns and their connections. which can lead to arbitrary results and overcomplicated graphs; First, the main nt var (i.e. the largest group of identical BcR IG nt (ii) BcR IG sequence classification is based on phylogenetic tree- sequences within a clonotype) of a given clonotype is identified. constructing programs, which is not ideal, given the high level of This represents the central point of the network since it comprises IgIDivA | 3 Figure 1. Workf low of the pre-processing of the data before the use of IgIDivA. The main steps are represented—from the obtention of the NGS raw data to the obtention of the two files (clonotypes computation and grouped alignment) that are used as input for IgIDivA. The original input is represented as an oval. The processes are represented as rectangles and the files are represented as parallelograms. The output files of the processes are the input files of the following ones. the largest fraction of sequences and, in turn, B cells of the respec- mentioned, all nt vars with fewer SHMs than the main nt var are tive clone. By default, SHMs of the main nt var are calculated from selected. In regard to nt vars with additional SHMs, only the nt the beginning of FR1 to the end of FR3. However, this analysis can vars belonging to the following pathways are selected for further be performed from a more downstream position within the BcR analysis (Figure 2): IG, and suggestions are given when selecting the parameters. (i) Most relevant pathway. Also named as the main block of Second, nt vars that have the same SHMs as the main nt var as SHMs, it is defined as the group of nt vars leading to a specific well as additional SHMs are selected. Whenever the main nt var end of pathway nt var (end node) and including the highest does not carry any SHM (i.e. 100% germline identity), all other nt number of BcR IG sequences within a given clonotype. vars will be assigned to that category. Samples that do not contain (ii) Longest pathway(s). They contain the highest number of nt vars with additional SHMs or if they contain only nt vars with a nodes. single additional SHM (one level of extra SHMs) will not be further (iii) Longest mutational pathway(s). They include connected analyzed. nodes with the highest cumulative number of mutations. Third, nt vars are connected according to their SHMs. This It is identical to the ‘longest pathway’ if ‘jumps’ are not process takes place only if the connection is consistent with allowed. progressive acquisition of SHMs. The tool also allows for the inclusion of ‘jumps’ (i.e. connection of nt vars with common SHMs differing by two or more SHMs). Additionally, a filter regarding the Graph metrics calculations minimum size of analyzed nt vars can be applied (default option is 10 sequences). In complex cases where a given nt var could High levels of complexity in the ID ref lect a strong selective be connected either to a nt var with only one SHM difference or antigenic pressure, leading to concrete, specific and non-random one with jumps (difference of two or more SHMs), only the first mutational patterns. To obtain an overview of such process, we connection would be shown for simplification purposes. introduced a series of graph network metrics. They can be used The last part of this process concerns the selection of nt vars to quantify different features related to the complexity of the ID with common yet fewer SHMs than the main nt var. These nt process, such as the level of convergence or the length of the vars are considered as having emerged prior to the main nt var. mutational pathways (Table 1). In brief, a sample with a high level Then, they are ordered according to their number of SHMs. At of complexity in the ID process would be characterized by high this step, all nt vars with fewer SHMs compared to the main nt ‘average degree’, ‘average distance’, ‘relative reads convergence’ var while also passing the threshold of the minimum number of score, ‘maximal path length’ and ‘maximal mutational length’ sequences are considered. The reason for this additional criterion as well as by a low ‘end nodes density’ score. The metrics ‘aver- is that these nt vars are extremely rare, given that they represent age degree’ and ‘average distance’ are established Graph Theory early steps of the ongoing SHM process. metrics [45]. The rest (‘relative reads convergence’, ‘end nodes density’, ‘maximal path length’ and ‘maximal mutational length’) have been specifically designed for ID analysis. Mutational pathways selection In the case of the relative reads convergence, a low value After the generation of the graph network, only the nt vars with would be close to 0, with the possibility of having cases with values SHMs considered as critical parts of the ID process within a >1 (if the nt vars with additional mutations accumulate more BcR given BcR IG clonotype are selected. In more detail, as previously IG sequences than the main nt var). For the end nodes density, 4 | Zaragoza-Infante et al. Figure 2. ID network and pathway selection. Representation of the complete networks and its components. The nodes represent nt variants (nt vars) and the edges represent the connections among them. The letters below the nt vars represent the mutations at the nt level that each nt var has. Note that the nt mutations are shown in this figure for the purpose of this example; they do not appear in the actual graphs, only the mutations at the aa level are shown, if any. The germline is represented in this figure for explanation purposes, but it does not appear in the actual graphs. The main nt is indicated at the center of the networks in a dark circle with three mutations (A–C). Toward it, nt vars lacking some of those mutations appear (the darker the circle, the less mutations they have and the closer to the germline they are). From the main nt var, nt vars having the same mutations (A–C) and additional ones are shown in circles orange and green (the colors are intercalated in order to distinguish better the different level of mutations; the darker the color, the higher is the number of additional mutations). The size of the nodes is proportional to the number of reads that constitutes the nt var. Panels A and B represent complete networks, the first allowing jumps (from +2 mutations to +4) and the second, without jumps allowed. As it can be seen in panel A, if jumps are allowed and a nt var (A + B + C + D + E + F) can be connected to the previous one (A + B + C + D + E) but also to another one with jumps (A + B + C + F), it will only be connected to the one with less difference in terms of number of mutations to simplify the graph. Panels C–F represent the pathways that are taken into account for the network construction. C represents the ‘less mutations pathway’ (the pathway that connects the nt vars between the main nt var and the germline). D represents the most relevant pathway (the one that gathers the maximal number of reads). E represents the longest pathway (length = 3). F represents the longest mutational pathway (+4 mutations). For simplicity, the three types of pathways displayed only contain one pathway each. It is, however, possible for each type to have more than one pathway (i.e. many pathways with the same maximal pathway length or with the same maximal mutational pathway length). the closer to 0, the lower the randomness of the ID process. For Graph network visualization the rest of the proposed metrics, there cannot be a single defined For graph network visualization, the main nt var is selected and range of values that can be readily considered as ‘low’ or ‘high’. depicted as the central point. Nodes that represent nt vars with In principle, the further a value is away from 0, the greater the fewer SHMs are depicted with directionality toward the main nt level of the complexity is for the respective metric. Ultimately, the var according to the connectivity pattern; in specific, the greater comparison between samples across multiple metrics can allow its SHM difference, the further away it is depicted. For the nt vars for better estimation of the ID complexity evaluation. with additional SHMs, only nodes and edges that fall into the pre- viously mentioned pathways are selected and represented with directionality from the main nt var. Depending on the selected Identification of replacement SHMs parameters, the network representation may also contain: (i) aa For each node, SHMs at the nt level are classified into silent and substitutions (this does not include the SHMs of the main nt var replacement SHMs, which lead to amino acid (aa) substitutions. in order to better visualize the acquisition of new replacement This process is performed for the main nt var as well as for all nt SHMs); (ii) the relative size of each individual nt var can be vars with additional SHMs. proportional to its number of sequences using a log10 scale and IgIDivA | 5 Table 1. Description of the graph metrics used for the characterization and comparison of ID Metric Calculation Meaning High value = Question that is answered by the metric RELATIVE READS Number of sequences of Tendency for the BcR IG High convergence Will the clone evolve toward the CONVERGENCE the most relevant pathways sequences to accumulate in accumulation of more SHMs or will it Number of sequences of the main nt var or to acquire remain in its current state? the main nt var additional convergent SHMs Number of end nodes END NODES Randomness or specificity of High randomness Is the acquisition of SHMs specific or Number of nt vars with DENSITY the mutational path random? additional SHMs MAXIMAL PATH Number of levels of additional Complexity of the High complexity Is this an ongoing SHM process LENGTH SHMs mutational path happening in many steps? MAXIMAL Maximum level of additional Complexity of the mutational High complexity Are there many different SHMs MUTATIONAL SHMs pathway, allowing happening progressively? LENGTH non-consecutive SHMs AVERAGE DEGREE Average total number of Complexity and connectivity High connectivity Are the mutations strongly connections of each nt var of the mutational pathways interconnected? AVERAGE DISTANCE Average number of steps along Complexity of the High complexity Is the ID process the shortest pathways between mutational path complex? each pair of nt vars (iii) SHMs jumps (if allowed, nodes with non-consecutive SHMs In more specific, the dataset consisted of 72 analyzed PCR will be shown). amplicons from 44 patients with CLL: 32 cases belonged to stereo- The main nt var is always depicted in dark petrol blue, whereas typed subset #2 (32 IGHV-IGHD-IGHJ, 21 IGLV-IGLJ gene rearrange- the nt var with fewer SHMs are represented using different shades ments), 7 cases belonged to stereotyped subset #169 (7 IGHV- of lighter blue (the darker the shade, the greater the distances IGHD-IGHJ, 6 IGLV-IGLJ gene rearrangements), while 6 cases were from the main nt var). For the nt vars with additional SHMs, not stereotyped (6 IGLV-IGLJ gene rearrangements). different shades of orange and green are intercalated, allowing for a rapid distinction of mutational levels and the identification of non-consecutive SHMs (nt vars with ‘jumps’). Jumps can be Results and discussion visually identified whenever two connected nodes have the same General output color and/or there is a gap in intensity between their respective The output of IgIDivA consists of a series of graphs and summary shades of colors (the darker the code, the greater the distance tables in order to provide a systematic analysis of the level of ID. from the main nt var). In more specific terms, the following types of output are provided for each given sample: Group comparisons (i) Summary calculations table. It contains information about Information regarding different groups and their respective sam- the number of clonotypes that are considered for analysis, ples should be provided by the user. Then, all groups are compared the number of nt vars, the number of total sequences, the for each individual metric. First, a Kruskal-Wallis test is applied number of singletons, the number of expanded sequences in order to determine whether the comparisons are overall sta- and the number of sequences of the main nt var. tistically significant [46]. Subsequently, pairwise comparisons are (ii) Extra mutations calculations. It provides the number of nt performed for all samples and the statistical significance of each vars with additional SHMs for each given number of SHMs as comparison is calculated with the Wilcoxon test [47]. The user can well as the total number of sequences. It includes the total select whether the p-values will be adjusted or not. For a given number of nt vars and sequences (i.e. not only the ones that metric, the corresponding values of each group are represented are plotted and used for the calculation of graph metrics). with side-by-side boxplots. The x-axis corresponds to the groups, (iii) Fewer SHMs calculations. This table shows the number of while the y-axis corresponds to the metrics values. Furthermore, sequences lacking SHMs of the main nt var for each different the plot contains the mean/median value as a red dot/horizontal number of SHMs. line within the boxplot of each group. Whenever a pairwise com- (iv) Evolution. It provides information for all unique SHMs or parison between two groups is statistically significant, the plot combinations of SHMs of all the nt vars that are part of will also contain a bar connecting these two groups as well as the the connected graph network. It also shows the number corresponding p-value directly above. of SHMs in comparison to the germline, the number of sequences with those SHMs and the mutational level to The dataset which they belong. In more detail, the mutational level is The dataset used for the validation of IgIDivA was extracted from ‘less’ if they have fewer SHMs than the main nt var, ‘main’ a previous NGS study focusing on the ID analysis in chronic for the SHMs of the main nt var and ‘additional’ for the lymphocytic leukemia (CLL) samples from stereotyped subsets #2 cases with more SHMs than the main nt var. This information and #169 [31, 48]. The dataset is publicly available in the repository may assist in the identification of SHMs that are represented European Nucleotide Archive (ENA) with the accession number in a given network at the nt level and their connectivity PRJEB36589. pattern. 6 | Zaragoza-Infante et al. (v) Block table. This table contains the following information subset #169 was 5101 variants (p-value = 1.000). The non-subset for each end node: (i) the block score (i.e. the ratio of the cases displayed lower ID levels for the dominant clonotype (p- total number of sequences of the nodes forming the block value < 0.05) compared to subset #2, but no statistically signif- of pathways that leads to that particular end node to the icant differences were found against subset #169. Despite these total number of sequences of all the nodes of the network findings, the analysis was not able to capture the full extent of with more SHMs than the main nt var); (ii) convergence score the ID process. Thus, we decided to analyze the ID process in these (i.e. relative reads convergence metric, Table 1); (iii) number two stereotyped subsets in CLL in a more systematic way. of sequences of the block of pathways; (iv) number of nodes We followed the same pre-processing of data as outlined in of the block of pathways; (v) id of nodes forming the block of the original study by Gemenetzi et al. [31] and then proceeded to pathways; (vi) pathways length and (vii) maximal mutational the dedicated ID analysis with IgIDivA. Because of differences in length of the block. From this table the group of nodes with the experimental approach for the PCR amplification and NGS the highest block score is selected as the ‘main block’ or the of the heavy and light chain sequences, the starting point of ‘most relevant pathways’. the alignment had to be different for the samples with IGHV- (vi) Graph network. It contains all features described in the ‘Net- IGHD-IGHJ gene rearrangements and the ones with IGLV-IGLJ work generation’ subsection of the Materials and Methods gene rearrangements (column start = 5 and 23, respectively). Thus, section. samples were divided in two different batches and analyzed with (vii) Graph info. It summarizes the values of the graph metric. IgIDivA (Table 2). It contains the germline identity %, the values of the graph Some of the samples were discarded from the analysis, as metrics (Table 1) as well as information related to those described in Table 3. metrics. Finally, comparisons were undertaken between the different (viii) Replacement mutations. It consists of two different tables groups for all metrics (running time = 6.5 s). In total, all the analy- with the information of the replacement SHMs in the main sis was performed in 7 min 26 s, with an average of 6.05 s/sample. nt var and the rest of the nt vars, respectively. Information The analysis was performed using a standard mid-range laptop, for all SHMs at the aa level and the number of sequences with two processor cores (Intel Core i3) and 4 GB RAM. carrying each mutation is provided. These changes might Evidence of ID was found in all three groups (Figure 3). Overall, have implications at the functional level. subset #2 cases displayed a higher level of complexity in the ID process, followed by #169 cases and non-subset cases. In When analyzing more than one sample, additional output is most cases, the most prominent differences concerned the produced: IGLV-IGLJ data. For example, in terms of convergence score, (i) Discarded samples table. It provides the names of samples the IGL sequences of subset #2 had mean/median values of that have been discarded from the analysis (e.g. samples 0.079/0.0105, followed by IGH sequences of subset #2 cases with no connections among nt vars). (mean/median 0.038/0.01), IGH subset #169 sequences (mean/- (ii) aa mutation tables. It contains all identified replacement median 0.011/0.003), non-subset IGL sequences (mean/median SHMs, including those present in more than one sample. 0.003/0.003) and, finally, subset #169 IGL sequences (mean/me- (iii) Metrics table. It includes the graph metrics information for dian 0.002/0.002). Nevertheless, differences between the groups all samples. If a sample has been discarded, the cause is were small and only a few pairwise comparisons displayed provided. statistically significant differences (Figure 4): (iv) Comparisons. If samples are classified into groups, the tool (i) Convergence score: Subset #2 IGH sequences versus non- performs pairwise comparisons for all groups. This is per- subset IGL sequences (mean/median 0.038/0.010 and formed independently for each of the graph metrics. 0.003/0.003, respectively; p-value = 0.045); subset #2 IGH sequences versus subset #169 IGL sequences (mean/median Application: ID assessment in CLL stereotyped 0.038/0.01 and 0.002/0.002, respectively; p-value = 0.012). subsets #2 and #169 (ii) Maximal mutational length: Subset #2 IGH sequences versus To assess the capacity of IgIDivA to provide a systematic and subset #169 IGL sequences (mean/median 6.35/4 and 2.2/2, multifaceted analysis of ID, we processed the data from the study respectively; p-value = 0.021). by Gemenetzi et al. [31, 48], concerning the analysis of samples (iii) Average degree: Subset #2 IGL sequences versus subset #169 from patients with CLL belonging to stereotyped subsets #2 and IGL sequences (mean/median 2.881/2.850 and 2.268/2.250, #169. respectively; p-value = 0.003); subset #2 IGL sequences In CLL, ‘stereotyped subsets’ are defined as groups of patients versus subset #2 IGH sequences (mean/median 2.881/2.850 presenting common immunogenetic features as well as common and 2.45/2.31, respectively; p-value = 0.008); subset #2 IGL clinical presentation and outcome [49–52]. Previous studies have sequences versus non-subset IGL sequences (mean/median shown that stereotyped subset #2, the largest in CLL, displays 2.881/2.850 and 2.502/2.500, respectively; p-value = 0.038). pronounced ID within the IG heavy and (lambda) light chain (iv) Average distance: Subset #2 IGH sequences versus subset genes [31, 53]. Stereotyped subset #169 is considered a satellite #169 IGH sequences (mean/median 1.94/2 and 1.286/1, to subset #2, i.e. it is closely related at the immunogenetic level respectively; p-value = 0.047). [54]. Gemenetzi et al. explored this relationship through the study of ID within the BcR IG utilizing NGS [31]. The authors reported All p-values were calculated using the Wilcoxon test [47] in a median of 3662 nt variants for the dominant heavy chain gene order to evaluate the statistical significance of changes in the clonotype of subset #2 cases and 1803 variants/dominant clono- metrics’ distribution medians between the provided groups. The type for subset #169 cases (p-value = 0.479). Regarding the light other metrics and pairwise comparisons did not show any sta- chain clonotypes, a median of 7015 variants/dominant clonotype tistically significant differences. This highlights the overall high was reported for subset #2 cases, while the respective value for level of similarity between the analyzed groups and, perhaps IgIDivA | 7 Table 2. Number of samples and running time of the two different batches of samples analyzed with IgIDivA IG gene rearrangements Number of samples Running time Running time/sample IGHV-IGHD-IGHJ 39 3 min 20 s 5.13 s IGLV-IGLJ 33 3 min 49 s 6.9 s Table 3. Number of samples discarded from the analysis and number of samples remaining for the analysis Subset Rearrangements Discarded (maximal Discarded (no connected Remaining for mutational length = 1) nodes with extra SHMs) the analysis #2 IGHV-IGHD-IGHJ 12 3 17 #2 IGLV-IGLJ 5 0 16 #169 IGHV-IGHD-IGHJ 0 0 7 #169 IGLV-IGLJ 1 0 5 Non-subset IGLV-IGLJ 1 0 5 Figure 3. Different levels of complexity and convergence in the ID process can be depicted with the IgIDivA graph networks. Panels A–C are subset #2 CLL samples (A and B are samples with IGHV-IGHD-IGHJ gene rearrangements, while C is a sample with IGLV-IGLJ rearrangement), D and E are subset #169 samples (IGHV-IGHD-IGHJ and IGLV-IGLJ gene rearrangements, respectively) and F is a non-subset sample (IGLV-IGLJ gene rearrangement). All networks have different levels of complexity and convergence of the ID (e.g. A has the longest pathways, B has only three nodes with additional mutations, C has many (121) end nodes, while D has only one end node and E and F have three and two end nodes, respectively). The figure shows the need for graph metrics to measure and quantify these characteristics and the differences or similarities among groups. For a better understanding, the values of the metrics convergence score and end nodes density of each sample has been added next to their respective identifier. The sample with the highest convergence score (C, CS = 1.017, indicating a higher accumulation of the BcR IG sequences in the nt vars with additional mutations compared to the main nt var), is also the sample with the highest END (indicating a higher level of end nodes and, thus, a higher level of randomness). This is an example of how it is necessary to have several metrics in order to correctly characterize one sample. A has the second highest convergence score, with a low END, showing a high level of convergence and complexity of the mutational pathways in the ID process. 8 | Zaragoza-Infante et al. Figure 4. Some of the pairwise comparisons between the sample groups displayed statistically significant results. Each set of boxplots represents a different graph metric (A, convergence score; B, maximal mutational length; C, average degree; D, average distance). Pairwise comparisons were performed for each metric among the groups (non-subset IGLV, subset #169 IGHV, subset #169 IGLV; subset #2 IGHV and subset #2 IGLV). The values are shown as boxplots, with each boxplot being a group, the y-axis representing the value of the metric, the horizontal line in the middle of the boxplot representing the median and the dot inside the boxplot representing the mean. If a pairwise comparison is statistically significant, the p-value is shown. For each metric, the p-value of the overall Kruskal Wallis comparison is shown as well. ‘Light’ = light chain, IGLV-IGLJ rearrangement; heavy = heavy chain, IGHV-IGHD-IGHJ rearrangement. most importantly, the ability of IgIDivA to identify even small cell receptor) gene repertoire and is publicly available. Alakazam differences in terms of ID that otherwise could remain unseen. has functionalities related to aa substitutions and graph metrics calculation and has been already used for the assessment of ID [30]. In terms of functionality, the tool offers the capacity IgIDivA versus other bioinformatic tools to study clonal lineages, diversity analysis, gene usage and lin- As previously mentioned, IgIDivA was designed in order to over- eage reconstruction [38]. This comparison between Alakazam and come limitations of previous bioinformatic tools [32–40]. First, IgIDivA may provide some insights regarding the corresponding IgIDivA does not perform any kind of inference for missing nt var strengths and challenges. To this end, we selected one sample connections and non-consecutive related mutations are allowed from the previously analyzed dataset (H33) and proceeded to a due to their potential relevance. This is particularly important similar type of analysis with both tools. Given the fact that the considering that inferred nt variants could lead to an artificially individual functionalities of both tools are quite different, a direct, higher level of complexity, especially in cases with extensive ID. quantitative comparison was not feasible; yet, some qualitative Also, IgIDivA selects certain pathways to display in the graph differences could be identified. (see Mutational pathways selection) leading to the generation of First, the input was modified in order to enable processing ‘cleaner’ graphs containing only relevant information. Moreover, it with Alakazam. The sample was analyzed with IMGT/HighV- IgIDivA is f lexible since a variable number of clonotypes can be QUEST and the output was extracted in AIRR format. Second, analyzed, depending on the user needs. Also, IgIDivA is publicly a script with a combination of some of the functionalities of available and easy-to-use without requiring a certain level of Alakazam had to be written due to the fact that this tool is not programming skills from the user. specialized for ID analysis and it has a wider range of uses: (i) For the comparison of IgIDivA with an already existing bioin- makeChangeoClone and (ii) buildPhylipLineage. Then, a lineage formatic tool, we chose Alakazam since it is a complete solu- tree was obtained, which was compared to the network obtained tion for the study of the adaptive immune receptor (IG and T with IgIDivA (Figure 5). IgIDivA | 9 Figure 5. Graphs obtained with Alakazam and IgIDivA. Sample H33 analyzed with (A) Alakazam and (B) IgIDivA. IgIDivA provides a more detailed characterization of the nt vars that are part of the ID process, including information of their size, their SHM content and the mutations at the aa level. Alakazam displays inferred nt vars (uncoloured circles) and depicts all the nt vars in the sample (coloured circles). IgIDivA does not infer intermediate nt vars but allows the presence of ‘jumps’ (two nt vars can be connected while having more than one SHM of difference) and the nt vars that are represented are the ones that are part of the ‘less mutational pathway’, the ‘most relevant pathway’, the ‘longest pathway’ and the ‘longest mutational pathway’ (i.e. the pathways that contribute the most to the ID process). Consequently, the results are more focused on the ID analysis and provide more meaningful results in that regard. Some of the differences found between the two methods were metrics are specific to the ID analysis, and comparisons the following: between groups of samples can be performed. (v) Purpose-built end-to-end tool: Unlike Alakazam, IgIDivA has (i) Processing time: Alakazam took approximately 3 min to been designed specifically for the analysis of ID; thus, it provide the output, while IgIDivA’s output was produced in offers a complete workf low that requires a minimal effort 10 s. from the user’s side. (ii) Presence of inferred nt vars: Alakazam includes in its lineage (vi) User-friendly tool: IgIDivA is also an easy-to-use Shiny app tree inferred nt vars together with all the nt vars present in without any need for programming knowledge. Moreover, the the data. In contrast, IgIDivA focuses only on real nt vars scripts are available and can be modified if necessary. included in pathways directly related to the ID process. (iii) aa mutations: Alakazam offers the possibility to study the Overall, Alakazam offers many useful functionalities for a aa content of sequences and its properties. However, it is complete repertoire analysis. However, IgIDivA is more specific in not possible to integrate this information into the network. the context of the ID process and provides a complete view in that IgIDivA plots the aa substitutions in the networks and pro- regard. vides tables with all the aa mutations present in all the samples. (iv) Graph metrics: Alakazam provides the option of calculating some metrics, however, most of them are calculated Conclusions individually for each nt var. For example, the function The study of ID in different contexts, extending from infections ‘getPathLengths’ calculates the distance from the germline and vaccination to autoimmunity, allergy and B cell malignancy, for each nt var; the ‘summarizeSubtrees’ function gives as can assist in gaining insight into the implicated processes, with result a table with some properties of each node (such as the diagnostic, prognostic and therapeutic relevance. Hence, dedi- size, the pathway length or the outdegree). While Alakazam cated tools can meaningfully assist in thorough immunogenetic provides comparisons between samples of the distributions analysis; against that, however, already available tools are not of normalized subtree statistics for a population of trees always specific for this purpose, have some drawbacks or do not with the ‘plotSubtrees’ option, it is not possible to get provide the necessary detailed characterization. For these rea- the metrics results per sample (instead of per node) or to sons, we developed IgIDivA, a purpose-built and publicly available make comparisons between groups of samples. IgIDivA’s tool for the complete and detailed analysis of the ID process in 10 | Zaragoza-Infante et al. high-throughput sequencing data. IgIDivA offers several function- References alities: (i) retrieval of different types of immunogenetic informa- 1. Meffre E, Casellas R, Nussenzweig MC. Antibody regulation of B tion, (ii) identification of SHMs and establishment of mutational cell development. Nat Immunol 2000;1:379–85. pathways, (iii) visualization of graph networks representing SHM 2. Rajewsky K. Clonal selection and learning in the antibody sys- connections involved in the ID process, (iv) calculation of new tem. Nature 1996;381:751–8. graph-based metrics that were developed for the evaluation of ID 3. Stewart AK, Schwartz RS. Immunoglobulin V regions and the B levels and (v) application of statistical analysis for the assessment cell. Blood 1994;83:1717–30. of comparisons between groups of samples. Moreover, IgIDivA 4. Nussenzweig MC. Immune receptor editing: revise and select. provides summary tables throughout the analysis that the user Cell 1998;95:875–8. can use for downstream study, while, at the same time, it is fast, 5. Sagaert X, Sprangers B, Wolf-Peeters D. The dynamics of the B open-source and easy-to-use for both experienced and inexperi- follicle: understanding the normal counterpart of B-cell-derived enced users. To our knowledge, IgIDivA is the first tool developed malignancies. Leukemia 2007;21(7):1378–86. specifically for the detailed study and characterization of the ID 6. De Silva NS, Klein U. Dynamics of B cells in germinal cen- process, offering relevant information, visual representations and tres. Nat Rev Immunol [Internet] 2015;15(3):137–48. Available from. useful metrics for quantification together with statistical analysis https://doi.org/10.1038/nri3804. for the comparison among samples. It can be applied for the study 7. Cattoretti G, Büttner M, Shaknovich R, et al. Nuclear and cyto- of ID in many different contexts and samples. plasmic AID in extrafollicular and germinal center B cells. Blood Future steps regarding the development of IgIDivA include 2006;107(10):3967–75. compatibility to the AIRR format as well as support of addi- 8. Crouch EE, Li Z, Takizawa M, et al. Regulation of AID tional visualizations of the output. Moreover, the f lexibility of expression in the immune response. J Exp Med 2007;204(5): IgIDivA offers the possibility of further future extensions, includ- 1145–56. ing more graph metrics, additional parameters or continuous 9. Methot SP, Di Noia JM. Molecular mechanisms of somatic hyper- customization. mutation and class switch recombination. Advances in immunol- ogy 2017;133:37–87. Key Points 10. Peled JU, Kuang FL, Iglesias-Ussel MD, et al. The biochemistry of • The availability of bioinformatic tools tailored for the somatic hypermutation. Annu Rev Immunol 2008;26:481–511. analysis of ID is limited. 11. Klien U, Goasens T, Fischer M, et al. Somatic hypermuta- • We present IgIDivA, a purpose-built tool that provides tion in normal and transformed human B cells. Immunol Rev a detailed characterization and visualization of the ID 1998;162(1):261–80. process, allowing for comparisons between samples or 12. Stavnezer J, Guikema JEJ, Schrader CE. Mechanism and regula- groups of samples accompanied by advanced statistical tion of class switch recombination. Annu Rev Immunol 2008;26: analysis. 261–92. • IgIDivA can provide relevant information about the ID 13. Dörner T, Giesecke C, Lipsky PE. Mechanisms of B cell autoim- process in infections, vaccination’s reactions, autoim- munity in SLE. Arthritis Res Ther 2011;13:243–12. munity, allergy and B cell malignancies, among others, 14. Itoh K, Patki V, Furie RA, et al. Clonal expansion is a characteristic possibly having implications in their diagnosis, prognosis feature of the B-cell repertoire of patients with rheumatoid and treatment. arthritis. Arthritis Res 2000;2(1):50–8. • IgIDivA is fast, publicly available and user-friendly and 15. Baranzini SE, Jeong MC, Butunoi C, et al. B cell repertoire diver- available also as a Shiny app web application. sity and clonal expansion in multiple sclerosis brain lesions. J Immunol 1999;163(9):5133–44. 16. Agathangelidis A, Chatzidimitriou A, Gemenetzi K, et al. Higher- order connections between stereotyped subsets: implications Data availability for improved patient classification in CLL. Blood. 2020. The data underlying this article are available in the repository ENA 17. Küppers R, Klein U, Hansmann ML, et al. Cellular origin of human with the accession number PRJEB36589, at https://www.ebi.ac.uk/ B-cell lymphomas. N Engl J Med 1999;341:1520–9. ena/browser/view/PRJEB36589. 18. MacLennan ICM. Germinal centers. Immunol Rev 1994;12(1): 117–39. 19. Jacob J, Kelsoe G, Rajewsky K, et al. Intraclonal generation of Acknowledgements antibody mutants in germinal centres. Nature 1991;354(6352): The authors would like to thank Giorgos Karakatsoulis, from the 389–92. Institute of Applied Biosciences, for his relevant insights in the 20. Clarke SH, Huppi K, Ruezinsky D, et al. Inter-and intraclonal statistical analysis. diversity in the antibody response to influenza hemagglutinin. J Exp Med 1985;161(4):687–704. 21. Schröder AE, Greiner A, Seyfert C, et al. Differentiation of B Funding cells in the nonlymphoid tissue of the synovial membrane of This work was supported in part by COSMIC (Combatting Disor- patients with rheumatoid arthritis. Proc Natl Acad Sci 1996;93(1): ders of Adaptive Immunity with Systems Medicine), a Marie Curie 221–5. European Training Network funded from the European Union’s 22. Bikos V, Karypidou M, Stalika E, et al. An immunogenetic signa- Horizon 2020 research and innovation program under grant agree- ture of ongoing antigen interactions in splenic marginal zone ment no. 765158. This work was also supported by ELIXIR, the lymphoma expressing IGHV1-2∗04 receptors. Clin Cancer Res research infrastructure for life-science data. 2016;22(8):2032–40. IgIDivA | 11 23. Sutton L, Kostareli E, Hadzidimitriou A, et al. Extensive intr- 38. Stern JNH, Yaari G, Vander Heiden JA, et al. B cells populating the aclonal diversification in a subgroup of chronic lympho- multiple sclerosis brain mature in the draining cervical lymph cytic leukemia patients with stereotyped IGHV4-34 recep- nodes. Sci Transl Med 2014;6(248):248ra107–7. tors: implications for ongoing interactions with antigen. Blood 39. Gupta NT, Vander Heiden JA, Uduman M, et al. Change-O: a 2009;114(20):4460–8. toolkit for analyzing large-scale B cell immunoglobulin reper- 24. Sanger F, Nicklen S, Coulson AR. DNA sequencing with toire sequencing data. Bioinformatics 2015;31(20):3356–8. chain-terminating inhibitors. Proc Natl Acad Sci 1977;74(12): 40. Foglierini M, Pappas L, Lanzavecchia A, et al. AncesTree: 5463–7. an interactive immunoglobulin lineage tree visualizer. PLoS 25. Gurrieri C, McGuire P, Zan H, et al. Chronic lymphocytic Comput Biol [Internet] 2020;16(7):e1007731–10. Available from. leukemia B cells can undergo somatic hypermutation and intr- https://doi.org/10.1371/journal.pcbi.1007731. aclonal immunoglobulin VHDJH gene diversification. J Exp Med 41. Morel B, Barbera P, Czech L, et al. Phylogenetic analysis of SARS- 2002;196(5):629–39. CoV-2 data is difficult. Mol Biol Evol 2021;38(5):1777–91. 26. Bagnara D, Callea V, Stelitano C, et al. IgV gene intraclonal diver- 42. Philippe H, Brinkmann H, Lavrov DV, et al. Resolving difficult sification and clonal evolution in B-cell chronic lymphocytic phylogenetic questions: why more sequences are not enough. leukaemia. British journal of haematology. 2006;133(1):50–8. PLoS Biol 2011;9(3):e1000602. 27. Jain R, Roncella S, Hashimoto S, et al. A potential role for anti- 43. Li S, Lefranc M-P, Miles JJ, et al. IMGT/HighV QUEST paradigm gen selection in the clonal evolution of Burkitt’s lymphoma. J for T cell receptor IMGT clonotype diversity and next generation Immunol 1994;153(1):45–52. repertoire immunoprofiling. Nat Commun 2013;4(1):1–13. 28. Chapman CJ, Zhou JX, Gregory C, et al. V(H) and V(L) gene analy- 44. Kotouza MT, Gemenetzi K, Galigalidou C, et al. TRIP - T cell recep- sis in sporadic Burkitt’s lymphoma shows somatic hypermuta- tor/immunoglobulin profiler. BMC Bioinform 2020;21(1):1–21. tion, intraclonal heterogeneity, and a role for antigen selection. 45. Koutrouli M, Karatzas E, Paez-Espino D, et al. A guide to con- Blood 1996;88(9):3562–8. quer the biological network era using graph theory. Front Bioeng 29. Agathangelidis A, Vlachonikola E, Davi F, et al. High-throughput Biotechnol 2020;8(January):1–26. immunogenetics for precision medicine in cancer. Semin 46. Kruskal WH, Wallis WA. Use of ranks in one-criterion variance Cancer Biol [Internet] 2021;84(September):80–8. Available from. analysis. J Am Stat Assoc 1952;47(260):583–621. https://doi.org/10.1016/j.semcancer.2021.10.009. 47. Wilcoxon F. Individual comparisons by ranking methods. Biomet- 30. Bagnara D, Tang C, Brown JR, et al. Post-transformation IGHV- rics Bull [Internet] 1945;1(6):80–3. Available from http://www.jstor. IGHD-IGHJ mutations in chronic lymphocytic leukemia B cells: org/stable/3001968. implications for mutational mechanisms and impact on clinical 48. Gemenetzi K, Psomopoulos F, Carriles AA, et al. Higher-order course. Front Oncol 2021;11(May):1–14. immunoglobulin repertoire restrictions in CLL: the illustrative 31. Gemenetzi K, Psomopoulos F, Carriles AA, et al. Higher-order case of stereotyped subsets 2 and 169. [Internet]. European immunoglobulin repertoire restrictions in CLL: the illustrative Nucleotide Archive (ENA), PRJEB36589. 2021. Available from case of stereotyped subsets 2 and 169. Blood 2021;137(14): https://www.ebi.ac.uk/ena/browser/view/PRJEB36589. 1895–904. 49. Stamatopoulos K, Agathangelidis A, Rosenquist R, et 32. Ismail WM, Tang H. Clonal reconstruction from al. Antigen receptor stereotypy in chronic lymphocytic time course genomic sequencing data. BMC Genom leukemia. Leukemia [Internet] 2017;31(2):282–91Available from. [Internet] 2019;20(Suppl 12):1002–11Available from. https://doi.org/10.1038/leu.2016.322. https://doi.org/10.1186/s12864-019-6328-3. 50. Messmer BT, Albesiano E, Efremov DG, et al. Multiple distinct sets 33. DeWitt WS, Mesin L, Victora GD, et al. Using genotype abun- of stereotyped antigen receptors indicate a role for antigen in dance to improve phylogenetic inference. Mol Biol Evol 2018;35(5): promoting chronic lymphocytic leukemia. J Exp Med 2004;200(4): 1253–65. 519–25. 34. Yang X, Tipton CM, Woodruff MC, et al. GLaMST: grow lineages 51. Widhopf GF, Rassenti LZ, Toy TL, et al. Chronic lymphocytic along minimum spanning tree for B cell receptor sequencing leukemia B cells of more than 1% of patients express virtually data. BMC Genom 2020;21(9):1–11. identical immunoglobulins. Blood 2004;104(8):2499–504. 35. Barak M, Zuckerman NS, Edelman H, et al. IgTree ©: creating 52. Agathangelidis A, Darzentas N, Hadzidimitriou A, et al. Stereo- immunoglobulin variable region gene lineage trees. J Immunol typed B-cell receptors in one-third of chronic lymphocytic Methods 2008;338:67–74. leukemia: a molecular classification with implications for tar- 36. Dunn-walters DK, Belelovsky A, Edelman H, et al. The dynam- geted therapies. Blood 2012;119(19):4467–75. ics of germinal centre selection as measured by graph- 53. Kostareli E, Sutton LA, Hadzidimitriou A, et al. Intraclonal diver- theoretical analysis of mutational lineage trees∗ . Dev Immunol sification of immunoglobulin light chains in a subset of chronic 2002;9(December):233–43. lymphocytic leukemia alludes to antigen-driven clonal evolu- 37. Abdollahi N. B cell receptor repertoire analysis in clinical con- tion. Leukemia 2010;24(7):1317–24. text: new approaches for clonal grouping, intra-clonal diver- 54. Agathangelidis A, Chatzidimitriou A, Gemenetzi K, et al. Higher- sity studies, and repertoire visualization. Immunol Sorbonne Univ order connections between stereotyped subsets: implications [Internet] 2021. Available from; (Doctoral dissertation, Sorbonne for improved patient classification in CLL. Blood, J Am Soc Hematol Université). https://tel.archives-ouvertes.fr/tel-03482030. 2021;137(10):1365–76.