Characterization of Chronic Lymphocytic Leukemia Immunoglobulin Rearrangements from Partial Read Sequencing

Abstract

The determination of the mutational status in the immunoglobulin variable region is an established prognostic biomarker for chronic lymphocytic leukemia (CLL). The length and inner variability of the variable, diversity, and joining (VDJ) rearranged sequences compromise B-cell clone characterization using next-generation sequencing (NGS), and a standardization is needed to adapt the procedure to the current clinical guidelines. Here, we develop a complete strategy for sequencing the variable domain of the immunoglobulin heavy chain (IGH) locus with a simple, low-cost, and efficient method that enables sequencing using shorter reads (MiSeq 150 × 2), allowing for faster results. Clonality and mutational status determination are performed within the same analysis pipeline. We tested and validated the method using 319 CLL patients previously diagnosed with IGH locus characterized using Sanger sequencing, along with 47 healthy donor samples. The analysis method follows a clone-centered consensus sequence strategy to identify B-cell clones and establish a clonal threshold specific for each patient's clonality profile, thereby overcoming the limitations of Sanger sequencing which is the gold standard used for determining immunoglobulin heavy variable (IGHV) mutational status.

Keywords: IGH locus; B cell; Chronic lymphocytic leukemia; Immune repertoire; NGS.

Graphical Abstract

Figure 1

Workflow for IGH locus characterization in CLL patients A. Sequencing with multiplexed FR primer sets and Illumina MiSeq 2 × 150 bp kit. B. Reconstruction of the region of interest by overlapping FR reads using the in-house pipeline B-MyRepCLL (https://github.com/afuentri/B-MyRepCLL). The pipeline generates a consensus sequence for each B-cell rearrangement followed by filtering steps to minimize artifacts. C. Repertoire structure determination. Automatic KNN classification distinguishes between CLL and healthy repertoires, followed by prioritization of the predominant rearrangements to identify clonal rearrangements (clonal/subclonal cut-off). VDJ, variable, diversity, and joining; KNN, K-nearest neighbors; FR, framework region; CLL, chronic lymphocytic leukemia; IGHV, immunoglobulin heavy chain variable; IGHD, immunoglobulin heavy chain diversity; IGHJ, immunoglobulin heavy chain joining.

Figure 2

Bioinformatics pipeline basis A. Reads from theoretical clonal IG rearrangements in a single sample (different colors represent different B-cell clones) are initially assigned independently to different IGHV@ and IGHJ@ alleles. Dotted lines indicate reads aligned to both V and J alleles. B. In the next step, these reads are used to infer the V–J allele combinations. Reads corresponding to each IGHV–IGHJ pairing are isolated and mapped against a joined reference of the specific IGHV–IGHJ pair. A consensus sequence is generated for each combination, representing an individual IG rearrangement. These consensus sequences are used to calculate the percentage identity against germline IGHV@ alleles. Asterisks represent the gap for IGHD sequence. Red and gray boxes indicate somatic hypermutation events and the junction sequence which is a gap in the reference alleles combination, respectively. C. CDR3 and IGHD sequence extraction is performed. CDR3 amino acid sequence is retrieved by searching for the conserved amino acid motifs (Cys104, Typ118, and WGXG) in different open reading frames. IGHD@ is detected as an insertion considering the combined sequences of IGHV–IGHJ alleles as reference. IG, immunoglobulin; CDR3, complementary determining region 3.

Figure 3

Optimization of the procedure A. Accuracy of KNN classification (k = 1–10) in a random split of the dataset between training and test sets. B. Scatter plot of the KNN test classification with clonal and polyclonal labels. C. F1-micro and F1-macro average accuracy scores for 10-fold cross-validation for KNN classification (k = 1–10). D. Box plot for MAX_DIFF values per sample grouped by polyclonal, 1CLONE, and >1CLONE. After Mann-Whitney U test, Bonferroni-corrected P values are annotated to show differences between group distributions (polyclonal vs. 1CLONE: P = 1.303E−27; polyclonal vs. >1CLONE: P = 5.562E−13; 1CLONE vs. >1CLONE: P = 3.213E−03). ns, not significant (P > 0.05); *, 0.01 < P ≤ 0.05; **, 0.001 < P ≤ 0.01; ***, 0.0001 < P ≤ 0.001; ****, P ≤ 0.0001. The scipy.stats Python module was used to perform the statistical test. MAX_DIFF maximum clonal difference within a sample.

Figure 4

Representation of clonal ratios in polyclonal (healthy), 1CLONE, and 2CLONE groups from the test samples Bars on the positive Y-axis represent the clonal percentage of the different IG rearrangements detected per individual, with different colors per donor. Bars on the negative Y-axis represent clonal percentage ratios between consecutive clones in a sample, ordered by abundance. The maximum clonal ratio is highlighted in red and the remaining clonal ratios in black. A. Polyclonal repertoire (five healthy donors). B. 1CLONE repertoire (24 samples with a single predominant clone). C. 2CLONE repertoire (10 samples with double predominant rearrangements).