Generating single-contig haplotype-resolved IGH assemblies

We have established a protocol for comprehensive characterization of the IGH locus using ONT ultra-long sequencing with adaptive sampling (Loose et al. 2016). This method was validated on samples from four donors—two healthy individuals (HD1 and HD2) and two MS patients (MS1 and MS2)—and the Epstein–Barr virus transformed lymphoblastoid cell line (EBV-LCL) HG002. The donors represent diverse ethnic backgrounds: The healthy donors are of African American and Asian ancestry (HD1 and HD2, respectively), and the MS patients are of European and South Asian ancestry (MS1 and MS2, respectively). We established a sequencing protocol with ONT adaptive sampling, in which reads from our region of interest (ROI) were enriched without the need for target amplification during library preparation (Fig. 1A). This approach yielded a mean read depth of approximately 30× over the IGH locus (Fig. 1B), representing a five- to sevenfold increase compared with the mean depth across the genome (Table 1). Furthermore, we obtained ultra-long reads with maximum read lengths exceeding several hundred kilobases, high read N50 values ranging from 65 to 95 kb, and read median Phred scores between 20.3 and 21.6 (Table 1; Supplemental Fig. S1).

Figure 1.

Overview of the experimental setup and bioinformatic pipeline. (A) Schematics of DNA sample preparation and sequencing. High-molecular-weight gDNA was extracted from peripheral blood mononuclear cells (PBMCs), isolated monocytes, or the HG002 Epstein–Barr virus lymphoblastoid cell line (EBV-LCL). The DNA was used as input for the preparation of ultra long-read ONT sequencing libraries. Sequencing was performed using an adaptive sampling strategy to enrich for reads mapping to the IGH locus. (B) Read depth from the four donors on a custom IGH reference. (C) Outline of the bioinformatic pipeline used to generate one-contig haplotype-resolved IGH assemblies.

Open in new tabDownload PowerPointLink to figure

Using these data, we developed a bioinformatic pipeline employing existing tools (Fig. 1C), which enabled the generation of high-quality, contiguous haplotype-resolved assemblies. The pipeline utilizes a custom IGH reference that incorporates several common SVs to enable reference-based phasing of reads. Next, the phased reads are separately assembled de novo. We assessed the accuracy of the initial haplotype-resolved assembly drafts, which allowed us to correct phase shifts or revise potential misassemblies (see below). The final corrected IGH loci in all four donors and in HG002 were assembled in single contigs without any apparent phasing ambiguities.

Assessing the accuracy of IGH assemblies

Large duplications pose a challenge for de novo assembly, especially in cases with high sequence similarity between duplicated segments. These SVs are particularly common at this locus, making it crucial to evaluate the assembly and correct any inaccuracies. To ensure the quality of our IGH assemblies, we employed several strategies for rigorous quality control. First, we used the Flagger pipeline to assess assembly reliability (Liao et al. 2023). In this approach, reads are mapped back to the corresponding IGH draft assemblies, and coverage distribution is analyzed to identify potential misassemblies. Three of the draft diploid assemblies showed indications of collapsed or erroneous sequences (Table 2). We verified true errors by checking for distinct patterns of single-nucleotide variants (SNVs) in the read alignments, or by identifying clear boundaries with soft-clipped reads surrounding the region. If assembly errors were confirmed, we manually corrected them, before reevaluating the new assembly. The size and type of errors detected in the different draft assemblies varied between the donors, which reflects the individual donors’ unique SVs. In the final IGH assemblies, only a few bases remained flagged as misassembled, none of which could be resolved upon manual inspection. These unconfirmed flagged errors represent either assembly errors too large for our reads, or artifacts owing to uneven read coverage or mapping.

Additionally, we sequenced samples from our four donors using PacBio HiFi technology and used these high-accuracy reads to assess the sequence precision of our ONT assemblies. The HiFi reads were mapped to the final IGH assemblies generated with ONT data, and variants were called using Clair3. No definitive SNVs with high read coverage and allele frequency >0.8 were found (Fig. 2A). However, three SNVs of note were still detected. One SNV with allele frequency of one was called in HD1, but in an area with relatively low HiFi read depth (11×). Two high-coverage SNVs were called with frequency ∼0.5 in HD2 and MS2, indicating either true SNVs or misalignments. Indel errors between assemblies and HiFi reads were more numerous than SNVs. In total, 734 indels of quality greater than one were called (Fig. 2B) in the collective assemblies >10 Mb, giving a rate of 7.2 indels per 100 kb. The majority of these (n = 612) were found within ±1 bp of a homopolymer ≥5 bp.

Figure 2.

ONT assemblies are concordant with HiFi reads Variant caller Clair3 estimated single-nucleotide variants (SNVs; A) and insertions/deletions (indels; B) between the donors’ HiFi reads and their haplotype-resolved IGH assembly from ONT data. (A) Scatter plot of HiFi read depth and variant allele frequency for SNVs from all donors. (B) Dot plot of quality score of indels called, stratified by different filtration. No filter (blue), indels within ±1 bp of homopolymers ≥5 bp (green), and indels within homopolymers or short tandem repeats of dinucleotide of seven or more repeats (orange).

Open in new tabDownload PowerPointLink to figure

Benchmarking with HG002

To further evaluate the performance of our long-read ONT method, we sequenced and assembled the IGH locus for the Genome in a Bottle (GIAB) HG002 reference material. HG002, also known as NA24385, is a thoroughly sequenced reference genome with extensive publicly available high-quality data sets with trio validation. Among these is the recently published telomere-to-telomere (T2T-GIAB Q100) complete diploid HG002 genome (Rautiainen et al. 2023). Using our method, we sequenced and assembled the HG002 IGH locus, resulting in one contig per haplotype. These assemblies were concordant with the Q100 assemblies within the >1.78 Mb encompassing the IGH, as no SNVs were found between the two, and only 102 small indels were detected (indel rate of 5.7 per 100 kb). When we mapped the Q100 Chromosome 14 and our ONT IGH assemblies to the custom IGH reference, we observed, as expected, that large portions of the HG002 IGH locus and corresponding genes were deleted relative to the custom IGH reference (Fig. 3A,B). Notably, our assembly of haplotype 2, corresponding to the paternal haplotype, was missing ∼900 kb of the ∼1.4 Mb IGH sequence in the reference, including all IGHD genes.

Figure 3.

Incompleteness of HG002 IGH assemblies. (A) Schematic representation of the diploid IGH assembly of HG002 aligned to our custom IGH reference. Haplotype 1 (H1) and haplotype 2 (H2) correspond to the maternal and paternal chromosomes, respectively, as derived from the T2T assemblies. The top track displays the IGH gene regions, and the bottom track shows the alignment intervals of HG002 assemblies as blue boxes, with large deletions indicated by solid red lines. (B) Annotation of HG002 IGH assemblies. Empty cells indicate deletion of a gene, and the color of the cell corresponds to the matched allele sequence from the IMGT database. Sequences with no identical match are denoted as “novel.”

Open in new tabDownload PowerPointLink to figure

Annotating assemblies and discovery of novel alleles

A specialized IGH annotation pipeline was applied to the final haplotype-resolved assemblies. Alleles were assigned according to their match identity to known alleles within the ImMunoGeneTics Information System (IMGT) database. Perfect matches to documented IMGT alleles were assigned their corresponding identifiers, whereas unmatched alleles were classified as “novel” and labeled accordingly. Across the four donor assemblies, we identified a total of 650 V, D, J, and C genes, from which we annotated 219 distinct alleles (Fig. 4A–C; Supplemental Table S1). The majority of these alleles were already characterized and cataloged by IMGT, whereas 21 of the remaining 49 alleles were found in VDJbase (as of 26.08.2024), and 28 were not found in either database. Notably, of the 69 constant genes annotated in our donors, about half (∼50%) lacked corresponding entries in either IMGT or VDJbase.

Figure 4.

Annotation of IGH assemblies. Diagrams depicting genes detected in donors’ IGH assemblies. Empty cell indicates deletion of a gene, and a split cell indicates duplication of a gene. The color of the cell corresponds to the matched allele sequence from the IMGT database; sequences with no identical match are denoted as “novel.” Genes are split by groups: V genes (A), D and J genes (B), and C genes (C). In A, a novel structural variant is denoted by an asterisk and represents a large tandem duplication of IGHV3-23 and/or IGHV3-23D, giving a total of seven copies of IGHV3-23. In C, IGHG4 and its duplicate (sometimes termed IGHG4A or IGHG4D) are joined on one row. The red dotted line indicates a duplication extending from IGHE to IGHA1.

Open in new tabDownload PowerPointLink to figure

Characterization of a novel large structural variant in the constant region

A potential large collapse in the IGH constant region was found in one of the draft haplotype assemblies of HD1 (Table 2). Inspecting alignment of the donor's reads to the draft assembly confirmed an increased mean read depth in the flagged region, supporting the presence of a structural anomaly. Further analysis of reads mapping to this region revealed intra-haplotype SNVs with two distinct variant patterns, neither of which completely matched the initial uncorrected assembly. A deeper inspection, involving the alignment of all donor ONT reads to our custom IGH reference, revealed three distinct variant patterns across a ∼100 kb interval in the IGHC region. This finding was corroborated by mapping all HiFi reads to the IGH reference, confirming the presence of three distinct variant patterns spanning the IGHE–IGHA1 interval (Fig. 5A).

Figure 5.

A 120 kb SV in IGHC. (A) Alignment of PacBio HiFi reads from healthy donor 1 (HD1) onto the hg38 reference genome. The top IGV track shows PacBio HiFi reads grouped by origin: centromeric (H1.1) and telomeric (H1.2) side of the duplication in haplotype 1 (H1), and haplotype 2 (H2). The middle track shows a single ultra-long ONT read (>246 kb) from H1 spanning the duplication, with half of it aligning as a supplementary read. The bottom track illustrates constant gene coordinates. (B) Schematic representation of the large duplication involving four constant genes (IGHE–IGHA1) within HD1 H1. (C,D) Gene usage frequency of each subisotype as determined by FLAIRR-seq in HD1. The duplicated genes (C) are stratified by which part of the duplication they originate, except for IGHG4, in which the H2 and H1.1 alleles were identical (ambiguous). (D) Frequency of usage of the nonduplicated genes from each haplotype. (E) Representative phylogenetic trees from clonal families containing both copies of duplicated constant genes. Branch lengths indicate mutational distances between nodes, and node colors denote subisotype and gene location.

Open in new tabDownload PowerPointLink to figure

To further dissect this haplotype, we leveraged our ultra-long ONT reads and identified very long reads (>200 kb) that spanned the reference twice over the IGHE–IGHA1 genes. This strongly suggested that two copies of this array were sequentially arranged on the same chromosome. Collectively, these observations pointed to the existence of a large (∼120 kb) duplication in the IGHC region, which had been collapsed during the de novo assembly process. Because of high sequence similarity between the duplicated segments, the variant had to be manually resolved using ultra-long reads spanning the SV, followed by polishing of the assembly with ONT reads. The corrected HD1 haplotype 1 assembly, now including the resolved duplication, showed no conflicts upon realignment of either the ONT or HiFi reads (Table 2).

A proposed model of the SV is presented in Figure 5B. None of the exonic nucleotide sequences of the duplicated genes were identical when compared pairwise. Between the two copies of IGHG2 and IGHG4, there were no coding differences despite three and five exonic SNVs, respectively. IGHA1 contained one nonsynonymous SNV, leading to an amino acid substitution at position 176. The centromeric copy matched the IGHA1*06 allele, whereas the telomeric copy was a novel allele with a D176E mutation. IGHE exhibited the greatest variation among the duplicated genes, with three synonymous and three nonsynonymous SNVs. The telomeric copy corresponded to the IGHE*02 allele, whereas the centromeric copy was a novel allele with mutations A330V, W334G, and L464I. This novel allele most closely resembled IGHE*05 from our reference set but differed from the IMGT allele owing to the A330V and W334G substitutions.

To determine the extent to which both copies of the duplicated genes were expressed in the BCR repertoire, we performed near-full-length adaptive immune receptor repertoire sequencing (FLAIRR-seq) on peripheral blood mononuclear cells (PBMCs) from the donor (Ford et al. 2023). Notably, we identified productive transcripts from both copies of IGHA1 and IGHG2 within the duplicated haplotype 1 (Fig. 5C). Because the centromeric copy of IGHG4 on haplotype 1 is identical to that on haplotype 2, calling the haplotype for IgG4 transcripts was impossible. The calculated frequency of gene copy usage for each subisotype was evenly distributed between the three copies of IGHA1 (ranging from 29% to 38%). In contrast, IGHG2 transcripts were predominantly derived from the centromeric copy of the gene on haplotype 1 (49%) and the gene on haplotype 2 (46%), whereas the telomeric IGHG2 on haplotype 1 accounted for only 4.9% of the IGHG2 transcript pool. The basis of this biased class‐switching remains unknown (Supplemental Text; Supplemental Fig. S2). For nonduplicated constant region genes, transcripts of IGHA2 and IGHG1 were evenly derived from both haplotypes, whereas IGHG3 from haplotype 2 accounted for 61% of the IGHG3 transcripts (Fig. 5D). Finally, we inferred clonal relationships among transcripts from haplotype 1 using SCOPer (Nouri and Kleinstein 2018) and identified numerous trees in which both copies of a duplicated gene were actively utilized within the same clonal family (example trees are shown in Fig. 5E).

Novel structural variants in the variable region

The donors also harbored several significant novel SVs within the IGHV gene region. In HD2, we characterized the first haplotype with seven copies of the IGHV3-23 gene (Fig. 6A), the longest expansion so far identified. CNV has been noted for this region previously (Sasso et al. 1995), including the recent association with IGHV3-23/D frequency in the expressed Ab repertoire (Rodriguez et al. 2023). We extracted the repetitive ∼10.6 kb segments containing IGHV3-23 and compared them to each other (Fig. 6B). The five centromeric segments had identical V gene sequence, corresponding to the IGHV3-23*03 allele. The two telomeric copies also shared identical nucleotide sequence corresponding to the IGHV3-23*01 allele. Notably, this donor had two additional copies of IGHV3-23*01 on the other haplotype, resulting in nine diploid copies of this V gene.

Figure 6.

Characterization of a seven-copy IGHV3-23 haplotype. (A) Structural variant giving a seven-copy haplotype of IGHV3-23 in HD2. Illustration shows the order of the seven homologous segments. (B) Synteny plot showing sequence identity and structural variation between the seven consecutive IGHV3-23-containing segments. Similarity between segments is represented by lines colored by the percentage identity of matched bases in 250 bp bins. (C) V gene usage in IgG1 transcripts derived from each haplotype in HD2, a “canonical” two-copy and the novel seven-copy IGHV3-23 haplotype. (D) Frequency of allele usage in IgM transcripts containing IGHV3-23 from the seven-copy haplotype.

Open in new tabDownload PowerPointLink to figure

To investigate whether the IGHV3-23 gene expansion translated to proportionally higher usage in the expressed BCR repertoire, we performed FLAIRR-seq on PBMCs from the donor. HD2 is heterozygous for IGHG1 and IGHM, which allowed us to differentiate IgG1 and IgM transcripts derived from the two-copy and the seven-copy IGHV3-23 haplotype. Our analysis revealed no difference in the proportional usage of IGHV3-23 between the haplotypes, in neither the IgG1 nor IgM repertoires (Fig. 6C,D), despite the extra five copies on one haplotype. In the seven-copy haplotype, IgM transcripts utilizing the IGHV3-23 gene were more frequently derived from the *03 allele than the *01 allele (Fig. 6D), aligning with the greater number of *03 allele copies.

In addition, HD2 carried single-copy IGHV1-69/IGHV2-70 haplotypes on each of their chromosomes. However, analysis of these haplotypes in the context of the GRCh38 two-copy IGHV1-69/IGHV2-70 haplotype (which harbors IGHV1-69D, IGHV1-69-2, and IGHV2-70D) indicated that they represent distinct deletion breakpoints. Breakpoint analysis of haplotype 2 was consistent with the deletion of IGHV1-69, IGHV2-70D, and IGHV1-69-2 (Supplemental Fig. S3D). However, determination of breakpoint in haplotype 1 was less conclusive as the haplotype is less congruent with the reference (Supplemental Fig. S3C). The haplotype lacks IGHV1-69-2 and IGHV1-69D, but fine mapping of the breakpoint indicates it occurred close to the boundary of, or within, IGHV2-70D or IGHV2-70. This means the deletion may not technically be a deletion of one or the other gene but rather a hybrid of the two. Critically, this is in contrast to the single-copy IGHV1-69/IGHV2-70 haplotype in MS1, which represents the deletion of IGHV1-69D, IGHV1-69-2, and IGHV2-70D, similar to the gene deletion profile previously characterized in GRCh37 (Rodriguez et al. 2023).

Hemizygous deletion or identical IGH haplotypes

One of the MS patients (MS2) lacked heterozygous variants necessary for accurate phasing of the IGH locus. Consequently, a nonphased single-haplotype assembly was generated from all reads in the locus. Although a small number of heterozygous variants were detected (Supplemental Fig. S4A), closer examination revealed that these variants had low read support and were located in highly repetitive regions of the locus, consistent with known sequencing errors associated with the ONT platform (Sereika et al. 2022). The quality distribution of the few heterozygous variants mirrored this observation, in that most variants were of low quality (Supplemental Fig. S4B). Furthermore, when corresponding HiFi reads were aligned on the haploid assembly, only one candidate SNV was called (Fig. 2A). These findings suggest that this individual either is hemizygous for the IGH locus or possesses two near-identical haplotypes. Additional coverage analysis of HiFi reads from the patient supports the hypothesis of a diploid locus, as no reduction in read depth was seen in IGH compared with other loci, and their coverage pattern mirrored those of the other donors (Supplemental Fig. S5). Accordingly, further bioinformatic separation of the two haplotypes was not warranted.

Ethics statement

The study and sample collection was approved by the regional ethical committee of SouthEast Norway (2009/23) and the University of Louisville institutional review board (IRB 14.0661). All donors provided written informed consent before study inclusion.

Sample collection and preparation for ONT sequencing

PBMCs were collected from two MS patients using Vacutainer CPT (Becton Dickinson Biosciences), before monocytes were negatively selected using the Pan Monocyte Isolation Kit, human (Miltenyi Biotec). Cells were frozen and kept in nitrogen until DNA extraction and sequencing as stated below.

Commercial frozen human PBMCs were purchased from StemCell Technologies. HG002 EBV-LCL were purchased from the Coriell Institute of Medical Research and expanded for 14 days in Gibco RPMI 1640 (Thermo Fisher Scientific) supplemented with 15% heat-inactivated fetal calf serum and penicillin/streptomycin at 100 U/mL and 100 µg/mL, respectively, before aliquots were freezed.

ONT DNA library preparation

Approximately 6 × 10⁶ cells, either monocytes (MS1 and MS2), PBMCs (HD1 and HD2), or EBV-LCLs (HG002), were used as input for DNA extraction and subsequent ultra-long sequencing library preparation. DNA was extracted using Monarch HMW DNA Extraction from Cell & Blood (New England Biolabs) following the protocol “High Molecular Weight DNA (HMW DNA) Extraction from Cells” (NEB T3050) with the following alterations: (1) the agitation speed during cell lysis was set to 1400 rpm for 10 min; (2) the incubation time with beads during gDNA binding to beads was increased to 8 min; (3) the elution buffer was replaced with 560 µL extraction EB from SQK-ULK114 (ONT); (4) an extra incubation step for 1 h at room temperature was added after elution with agitation at 56°C; (5) the centrifugation of beads with DNA was increased to 1 min at 16,000g; and (6) an additional 200 µL extraction EB buffer was added after complete elution to get correct input volume for ONT library preparation. Sequencing libraries were generated using an Ultra-Long DNA Sequencing Kit from ONT (SQK-ULK114) following the manufacturer's guidelines. In line with recent literature and ONT documentation, we use the term “ultra-long” for libraries with N50 > 50 kb that also include a substantial proportion of bases in reads >100 kb.

ONT sequencing

Sequencing was performed on a PromethION 2 solo (P2 solo) device connected to a high-performance computer running MinKNOW software version 23.07.12. The loading and washing of PromethION R10.4.1 flow cells were performed as per the manufacturer's guidelines. Each donor was sequenced using one flow cell, and the sequencing ran for a maximum of 72 h with washing and reloading performed around every 24 h using the Flow Cell Wash Kit (EXP-WSH004) to increase output.

To serve as an IGH reference for adaptive sampling during sequencing, T2T-CHM13v2.0 was altered at the telomere end of Chromosome 14. Here, the variable part of IGH was replaced with a custom IGH reference as described elsewhere (Rodriguez et al. 2020), whereas the rest of Chr 14 including the IGH constant region remained. Adaptive sampling was enabled in MinKNOW to enrich for the ROI, specified as the IGH region with an additional 200 kb lead-in from the custom CHM13 reference. During sequencing, the high accuracy basecalling model was enabled to accommodate adaptive sampling. The resulting POD5 files were rebasecalled post-run using the superaccurate model, filtering out reads with a Phred score below 10, and were subsequently used for further downstream processing.

Processing ONT data to IGH assemblies

A pipeline to process ONT reads into haplotype-resolved IGH assemblies was established using existing tools and custom scripts. First, catfishq (version 1.4.0) (https://github.com/philres/catfishq) was used to filter out adaptive sampling rejected reads from the rebasecalled FASTQ files. Next, the reads were aligned to our CHM13v2.0 reference with the custom IGH described above using minimap2 (v 2.26) with the option “-ax map-ont” (Li 2018). Variant calling was performed on the resulting BAM file with Clair3 (v 1.0.0) (Zheng et al. 2022) with the options “‐‐enable_phasing”, ‘‐‐platform=’ont”, and “‐‐bed_fn” enabled, specifying a BED file with IGH coordinates (chr14:99630000-101325184). Reads were haplotagged using Clair3 <phased_merged_output.vcf> with WhatsHap (v 1.4) (Martin et al. 2016) haplotag and the options “‐‐output-haplotag-list ‐‐ignore-read-groups ‐‐tag-supplementary ‐‐region=chr14:99630000-101325184”. Continuity of haplotype tags between phase sets was evaluated by manual inspection of the haplotagged BAM in the Integrative Genomics Viewer (IGV) (Robinson et al. 2017). If switch errors were detected, the Clair3 VCF was edited before running “whatshap split” with default parameters to produce hap1 and hap2 phased FASTQs. Lastly, de novo assembly was produced from the haplotype-separated reads using Flye (v 2.9.2) (Kolmogorov et al. 2019) with the arguments “‐‐genome-size 1200K ‐‐nano-hq <phased.fastq>”.

Manual inspection in IGV and Flagger (v 0.4.0) (Liao et al. 2023) was used to evaluate haplotype-resolved assemblies. ONT reads from the ROI were mapped back to the diploid assembly with minimap2- and Flagger-highlighted misassembled regions. Problematic regions of the assemblies had to be manually resolved when collapsed by the assembler. In the example of the collapsed large duplication in the IGHC region, we isolated reads that spanned the duplication. Subsequently, one of the longer isolated reads was extracted to a FASTA and polished with the rest of the isolated reads using “flye ‐‐nano-hq <isolated_reads.fastq> ‐‐polish-target <long_read.fasta> ‐‐iterations 3”. The polished long read was mapped back to the first-draft assembly with minimap2 with the “-ax asm20” option, and coordinates for cutting and pasting of the resolved duplication were determined. The entire collapsed sequence was removed, and the polished long read was trimmed to fit this gap in the assembly and inserted here. The resulting new draft assembly was finally verified again as outlined above with Flagger.

PacBio HiFi sequencing of IGH using capture probes

Parallel sequencing of samples from all donors was performed according to a previously published method (Rodriguez et al. 2020). In brief, probe-based targeted capture of the IGH region and long-read single-molecule, real-time (SMRT)-seq was used to sequence IGH reads. Libraries were sequenced on PacBio instruments SEQUEL II (HD1, MS1, and MS2) or REVIO (HD2).

IgA and IgG FLAIRR-seq

RNA was extracted from HD1 and HD2 PBMCs using the AllPrep DNA/RNA Mini Kit (Qiagen). IgG and IgA transcripts were resolved by targeted amplification using FLAIRR-seq as described previously (Ford et al. 2023). FLAIRR-seq cDNA libraries were prepared using the SMRTbell prep kit 3.0 (PacBio) according to the manufacturer's specifications. The prepared library was sequenced on the Revio long-read system (PacBio). FLAIRR-seq data were processed via the Immcantation tool suite using pRESTO and Change-O as previously described (Vander Heiden et al. 2014; Gupta et al. 2015; Ford et al. 2023). For the identification of clonal relationships, we used SCOPer hierarchical clustering with a threshold of 0.15 (Gupta et al. 2017). Dowser and Igphyml were used to build and visualize the clonal trees (Hoehn et al. 2019, 2022).

Validation of IGH assemblies

Personalized alignment of HiFi long reads were used to infer orthogonal validation, and base-level concordance was obtained to ensure the highest quality of haplotype-resolved assemblies. The HiFi reads were aligned with minimap2 to the corresponding ONT IGH assemblies with the option “-ax map-hifi”. The variant caller Clair3 was used to assess base accuracy as a measure of agreement between the high-accuracy reads and the polished assemblies. Clair3 image was ran with singularity and the options “‐‐platform=‘hifi’ ‐‐model_path=‘hifi_sequel_2’ ‐‐no_phasing_for_fa”.

The alignment of HiFi reads were also used to infer orthogonal validation, and base-level concordance was obtained with other metrics. Custom Python scripts were utilized to assess the proportion of assembly bases with at least 80% concordance from the reads, ensuring a minimum read depth of 5× and 10× for estimating coverage and accuracy.

Annotating IGH assemblies

IGH assemblies were annotated in a reference-guided approach using custom Python and bash scripts available at GitHub (https://github.com/Watson-IG/wasp) and as Supplemental Code. Finished ONT assemblies and reads were aligned to a human reference genome (Rodriguez et al. 2020) using minimap2 with the presets “-ax asm20” and “-ax map-ont”, respectively. Sequences for IGH exons were extracted from assembly BAM using gene coordinates mapped to the reference from GRCh38 and searched against available human IGH alleles in IMGT database (IGHV, -D and -J sets were downloaded February 13, 2024; IGHC alleles were collected November 24, 2022). HiFi reads mapping over the same intervals were used to validate support for assigned and novel alleles, using the same Python scripts and metrics that were used to assess read coverage in the assembly as a whole.