A multiomic characterization of the leukemia cell line REH using short- and long-read sequencing

Overview of the genome

The human pre–B-cell ALL cell line REH was obtained from DSMZ (ACC 22), cultured according to the supplier’s specifications, and authenticated using STR analysis. G-banded karyotyping was then used to verify the presence of the chromosomal features described previously in literature (Uphoff et al, 1997) and by the cell line vendor (DSMZ) (Fig 1A). Karyotyping confirmed most of the large chromosomal aberrations expected in REH: a deletion on the short arm of chromosome 3, trisomy 16, the balanced t(5;12), and monosomy X. It did not, however, completely resolve the four-way t(4;12;21;16). A comparison of the features resolved in the different karyotypes is available in Table S1.

The REH cells were subjected to genomic and transcriptomic sequencing. Of 33 datasets generated (Lysenkova Wiklander et al, 2023), the present study used three whole-genome sequencing (WGS) and two RNA-seq datasets (Table 1). The Illumina WGS dataset had an average depth of coverage of 34 reads, whereas ONT had 18 and PacBio 15; the Illumina WGS was sequenced PE-150 bp, while the ONT and PacBio median read lengths were 13,261 bp and 23,761 bp, respectively (Fig 1B). The IsoSeq median read length was 3,747 bp (Fig S1), whereas the Illumina RNA-seq was sequenced PE-100 bp. All sequencing datasets had an error rate of less than one percent, except for ONT (6.55%). Further sequencing statistics of these datasets are detailed in Table S2.

Three sets of SNV calls were generated by running DeepVariant on the WGS datasets. The variant allele frequency of the SNVs was used to assess uniparental disomies and other potential losses of heterozygosity (LOH) in the cell line, using the SNV callsets generated from Illumina (Fig 1C), ONT (Fig S2A), and PacBio (Fig S2B) data. In addition to confirming the del(3), +16, and −X, a partial trisomy of chromosome 21 was detected, suggesting two copies of the der(16) chromosome resulting from the t(4;12;21;16). Copy neutral loss of heterozygosity (cnLOH) was also detected across much of the short arm of chromosome 9, stretching from p13.2 (between exons 4 and 5 of RNF38) to the p-terminal.

Three SV callsets were created by running TIDDIT on Illumina short-read WGS data and Sniffles on PacBio and ONT long-read WGS data, from which consensus sets were subsequently generated using the SURVIVOR algorithm (Jeffares et al, 2017). To rule out the presence of subclones, the allele fractions of the three-way consensus callset, containing the SVs called in all three WGS datasets, were plotted to ensure that they follow a binomial distribution (Fig S3).

Two haplotype-aware de novo assemblies were generated to provide chromosomal context for SVs. The first assembly was created using PacBio circular consensus sequencing (CCS) reads and hifiasm (Cheng et al, 2021) with no pre-processing or polishing. The second assembly was created using ONT ultralong reads, Flye (Kolmogorov et al, 2019), and Medaka polishing. The hifiasm assembly generated 2,580 contigs, N50 of 3.5 Mb and L50 of 255, which was outperformed by the Flye assembly (2,060 contigs, N50 of 58 Mb and L50 of 20). The Flye assembly was selected for phasing analysis.

The SV landscape stratified by size

In total, TIDDIT generated 5,726 unfiltered SV candidates from the Illumina data, whereas Sniffles generated 36,121 candidates from the PacBio and 36,648 candidates from the ONT data. The distribution of these SVs over chromosomal loci was assessed with heatmap visualizations of the long-read consensus callset (Fig 2A), the TIDDIT callset (Fig S4A), the PacBio and ONT Sniffles callsets (Figs S4B and S4C), and the three-way consensus callset (Fig S4D). Of SVs with length > 100 bp, there were 10,262 variants in the long-read consensus set (called in both ONT and PacBio), whereas the three-way consensus set (intersection of all three WGS callsets) contained 2,072 variants. The non-translocation SVs were binned by size into Small (100 bp–1 kb), Medium (1–10 kb), and Large (>10 kb); in these categories, the long-read consensus sets contained 8,993, 1,420, and 19 variants, respectively, whereas the three-way consensus sets contained 1,485, 414, and 34 variants (Fig 2B). The frequencies of SVs stratified by size were visualized using strip plots (Fig 2C).

The original Sniffles and TIDDIT callsets were then programmatically filtered (see the Materials and Methods section). After automated filtering for large-scale SVs (>100 kb) and interchromosomal translocations, there remained 694 SV candidates from Illumina, 127 from PacBio, and 35 from ONT. All remaining Sniffles candidates were manually inspected in IGV, whereas a pseudorandom sampling from the TIDDIT set selected 29 candidates across 16 chromosomes for further inspection.

Refinement of known REH aberrations at base-pair resolution

The established and expected chromosomal aberrations were detected in each of the three SV callsets (Table 2). These features include a deleted chromosome X, a gain of chromosome 16, a 26-Mb del(3)(p22.3p14.2), a balanced translocation between chromosomes 5 and 12, and, finally, a four-way rearrangement between chromosomes 4, 12, 21, and 16, resulting in the ETV6::RUNX1 fusion gene. The t(5;12) and the four-way translocation were determined to involve two different homologs of chromosome 12, as the two p-arms of chromosome 12 (from p13.2 to the p-terminal) were found to be fused to either chromosome 5 or chromosome 21.

The molecular events of the four-way translocation were reconstructed as follows. The p-arm of chromosome 12 (12 Mb), with a breakpoint at p13.2 in intron 5 of ETV6, was translocated to chr21q22.12, resulting in the canonical subtype-defining fusion gene ETV6::RUNX1. The q-arm of chromosome 21 (12 Mb), with a breakpoint at q22.12 in intron 1 of RUNX1, was translocated to chr16q24.3. The q-arm of chromosome 16 (271 kb), with a breakpoint at q24.3 in intron 4 of PRDM7, was translocated to chr4q32.1. Finally, the q-arm of chromosome 4 (30 Mb), with a breakpoint at q32.1, was translocated to chr12q23.1 (where chromosome 12 had undergone an 85-Mb inversion between p13.2 and q23.1), completing the reciprocal exchange.

Combining the different SV callsets facilitated the resolution of these breakpoints down to base-pair accuracy, thereby resolving three inconsistencies in the previously documented REH karyotypes. First, each of the three callsets resolved the breakpoints of the balanced t(5;12) to 5q23.2-23.3, with a 2.2-Mb deletion on chr5q23 occurring between the breakpoints. The previous karyotypes each erroneously placed the breakpoint on chromosome 5 at q31-q32. Second, the breakpoints of the del(3) were unclear from the karyotypes, ranging from p13-p22 to p21.3-p24, or missing altogether. The SV callsets unanimously resolved this deletion to p14.2-p22.3. Finally, the karyotypes showed a discrepancy in the endpoint of the inv(12)(p13), placing it at either q22 or q23; the callsets resolved its location to q23.1 (Table S1).

Of note, the DSMZ karyotype documents the presence of a sideline containing an extra chromosome 18 and an inversion at chromosome 5. Neither of these aberrations were detected by our STR analysis or G-banding, indicating that the REH cells analyzed by us did not contain this sideline. These aberrations are also absent from the karyotype generated by Uphoff et al (1997).

Known aneuploidies were confirmed in the WGS data: a deletion of chromosome X was detected in all datasets by a marked decrease in the average depth of coverage (DC), whereas an additional chromosome 16 was detected by a corresponding increase (Table S3). The two copies of derived chromosome 16 resulting from the translocation t(16;21) were confirmed through DC and analysis of reads mapping to both chromosomes 16 and 21. The breakpoint was supported by 18 ONT split reads (average ONT DC chr16 = 26.4 versus 18.2 for diploid chromosomes), 19 PB split reads (average PB DC chr16 = 20.8 versus 14.6 for diploid chromosomes), and 32 Illumina discordant read pairs (average Illumina DC chr16 = 51.5 versus 34.7 for diploid chromosomes). A single set of breakends was found for this rearrangement (chr21:34947932, chr16:90067326), suggesting that this aneuploidy was caused by an aberrant mitotic event occurring after the translocation.

Large-scale SVs discovered in the REH genome

In total, 15 intrachromosomal and eight interchromosomal SVs were found (Fig 3), of which one interchromosomal translocation and 14 intrachromosomal SVs (deletions, duplications, and inversions >100 kb) were previously unidentified in the karyotypes (Table 2). A rearrangement was identified involving a 1.4-Mb segment of chromosome 2 inverted at p11.2, 58 kb of which was inserted into chromosome 1, resulting in derived chromosome der(1)inv(2)(p11.2p11.2)ins(1;2)(q21.1;p11.2). Our findings also include a 116-kb duplication on chromosome 1p21.1, a 1.1-Mb inversion on chromosome 16p12.2, and deletions on eight different chromosomes. These include a deletion on chromosome 3 at q26.32 (146 kb), deleting exons 2–5 of the gene TBL1XR1, which has been implicated in GC drug resistance (Jones et al, 2014) (Fig S5). A deletion on chromosome 5q31.3 (205 kb) was found to focally delete exons 2–9 of the GC receptor gene NR3C1 (Xiao et al, 2019) (Fig S6). Meanwhile, LOH on the p-arm of chromosome 9 leading to the homozygosity of a 24.5-Mb deletion at p21.3 was shown to delete both alleles of CDKN2A, which may be associated with inferior outcome (Kathiravan et al, 2019). A deletion on chromosome 12q21.33 (260 kb) was found in the region encoding the ALL-associated gene BTG1 (van Galen et al, 2010) (Fig S7), whereas three deletions were found on chromosome 18, one at q21.1 and two at q23, including the 132-kb deletion affecting exon 10 of the leukemia-associated gene NFATC1 (Medyouf & Ghysdael, 2008) (Fig S8). Finally, two deletions were discovered on chromosome 22q11.22, of which a 214-kb variant was found to entirely delete the gene VPREB1, which encodes the surrogate light chain involved in the formation of the pre–B-cell receptor (pre-BCR) and whose genetic loss contributes to leukemogenesis (Mangum et al, 2014) (Fig S9).

Of the 23 confirmed SVs, a total of 19 were called in all three SV callsets, including the eight aberrations known from the karyotypes and 11 of the novel SVs. An additional two SVs, in repetitive regions of chromosomes 1 and 2, were called in both PacBio and ONT, but not in Illumina data, whereas one SV was called in both Illumina and ONT. Finally, the homozygous deletion in chromosome 9 was only detected in Illumina data. All of the confirmed SVs were supported by at least three split reads from both ONT and PacBio, and at least three discordant read pairs from Illumina, with an average read support of 9.6 from PacBio, 11.2 from ONT, and 15.9 from Illumina (Table 2).

We evaluated the strengths and weaknesses of short- and long-read sequencing technologies in the context of characterizing the REH cell line. Overall, Sniffles run on ultralong ONT reads outperformed in the discovery of large-scale rearrangements, showing high sensitivity and a low false-positive rate relative to the other datasets. Sniffles/ONT called 22 of the 23 confirmed interchromosomal or large-scale SVs > 100 kb (95.65% sensitivity), whereas both Sniffles/PacBio and TIDDIT/Illumina called 21 (91.3%). The Sniffles/ONT callset contained the smallest number of false positives (128; 85.33% FPR), which were fewer than the Sniffles/PacBio callset (606; 96.65% FPR) and TIDDIT/Illumina callset (1,469; 98.59% FPR) (Fig S10A; Table S4).

The TIDDIT callset had a high sensitivity rate, detecting all of the confirmed SVs except for those occurring in highly repetitive regions. It was also the only dataset to detect the homozygous deletion on chromosome 9, despite high coverage of the breakpoints with long reads. However, it suffered from a high false-positive rate. Of the 29 SV candidates randomly sampled from the filtered TIDDIT callset, none were visibly supported by long reads in IGV; these SV candidates were called in regions containing indels >10 bp, which were resolved in the long-read data but mismapped and misidentified as translocations or other SVs in the Illumina data.

PacBio data, which had lower average DC than ONT (15 PB versus 18 ONT), did not detect any SVs that were missed in the other two callsets; however, in both the PB and ONT datasets, Sniffles detected an ins(1;2) and a 116-kb dup(1), both missed in the Illumina data and both occurring in highly repetitive regions indicated by spikes in DC (Fig 3).

The long-read technologies facilitated the confirmation of large-scale aberrations using split-read analysis. The visualization tool Ribbon was used to visualize reads spanning complex breakpoints such as the one between chromosomes 12, 16, and 21, and to explore how specific split reads map to the reference genome (Fig S10B). The ultralong ONT reads in particular helped resolve the rearrangement der(1)inv(2)(p11.2p11.2)ins(1;2)(q21.1;p11.2), where several ultralong reads spanned the entire 58-kb translocated region of chromosome 2 and its flanking regions in the derivative chromosome 1 (Fig S10C).

Phasing of the large-scale SVs was performed by creating a de novo assembly of the ONT reads with Flye and haplotyping with Hapdup, followed by alignment of the haplotypes to the reference genome and SV calling using Hapdiff. One can infer whether a pair of SVs occurs on the same homolog if they occur on the same contig or, alternatively, if a single contig covers the breakpoints of both events and contains neither of the SVs. The three large-scale deletions on the q-arm of chromosome 18 were mapped to a single contig, indicating that they all occur on the same homolog. However, this approach struggled to properly resolve rearrangements with higher complexity; for example, it mapped the large-scale inv(12) to the same homolog as the t(5;12), contradicting cytogenetic evidence. Furthermore, of the 31 SV breakpoints that were assigned to a haplotype, 14 of the contigs containing heterozygous SVs (45%) were erroneously assigned to both maternal and paternal haplotypes, whereas Hapdiff only called 10 of the 24 confirmed large-scale SVs (Table S5).

Seven expressed fusion genes, including two in-frame fusions

Short- and long-read RNA-sequencing data were used together with seven fusion gene callers to detect potential fusion genes in REH. In total, the unfiltered fusion callset consisted of 11,099 fusion gene candidates. After programmatic filtering, 30 candidates remained. Additional manual examination of split and spanning reads, as well as discordant mate pairs, retained seven high-confidence fusion candidates that could be confirmed in IGV with supporting genomic breakpoints in the WGS data (Fig 4).

Of the seven confirmed fusion genes, two were in-frame. First, the expected ETV6::RUNX1 fusion, resulting from the canonical subtype-defining translocation t(12;21), was detected as two splicing variants. Transcripts of both variants retained the sterile alpha motif/pointed domain from ETV6, and the Runt and Runx inhibition domains from RUNX1. ETV6::RUNX1 was detected in both the Illumina RNA-seq and IsoSeq datasets (five and 62 supporting reads, respectively) with support from five fusion callers. The two splice variants of ETV6::RUNX1 were in line with the genomic breakpoints found in intron 5 of ETV6 and intron 1 of RUNX1; however, in the most prevalent ETV6::RUNX1 transcript, exon 5 of ETV6 was spliced to exon 2 of RUNX1 (Fig 4A).

Five splicing variants of a second in-frame fusion gene, RUNX1::PRDM7, resulted from the t(16;21) occurring in two copies of the der(16). In the splicing variant involving exon 5 of PRDM7, the RUNX1::PRDM7 fusion retained the promoter from RUNX1 and the SSXRD motif from PRDM7. RUNX1::PRDM7 was highly expressed, detected by all seven fusion callers, and supported by 144 Illumina reads and 24 IsoSeq reads. The genomic breakpoint in PRDM7 was in intron 4; five splicing variants of RUNX1::PRDM7 were found involving exons 5–11, with the most prevalent fusion transcript taking place between exon 1 of RUNX1 and exon 9 of PRDM7 (Fig 4B).

Two out-of-frame fusion genes were found to arise from the balanced t(5;12): PHAX::AC007450.2 and LRP6::SLC27A6; both were detected in both the short- and long-read RNA-seq datasets. In PHAX::AC007450.2, which was supported by nine Illumina reads and five IsoSeq reads, exon 4 of PHAX was fused with exon 2 of AC007450.2, which lies 85 kb downstream of ETV6 (Fig 4C). Exon 22 of LRP6 was fused with exon 2 of SLC27A6 (Fig 4D). LRP6::SLC27A6 was supported by one short read and two long reads. However, manual inspection of the reads revealed an additional 11 low-quality split long reads that were discarded by the fusion caller and not reported in the supporting read count.

A fusion transcript arose from the 260-kb deletion del(12)(q21.33q21.33), which resulted in a truncated BTG1 gene, with a breakpoint in exon 2. In this fusion, the truncated BTG1 transcript was fused with an expressed non-coding region originating between two long non-coding RNA genes LINC02404 and AC090049.1. The BTG1::LINC02404/AC090049.1 fusion was only called in the short-read RNA-seq dataset using Arriba, but by no long-read callers. The resulting out-of-frame truncated transcript was highly expressed with 120 supporting Illumina reads and nine IsoSeq reads (Fig 4E).

The 205-kb del(5)(q31.3q31.3) resulted in the antisense transcript NR3C1::ARHGAP26, fusing exon 1 of NR3C1 with exon 20 of ARHGAP26. This fusion was only detected by one short-read fusion caller, with 10 supporting Illumina reads, but also found support in the long-read data, with five supporting IsoSeq reads (Fig 4F).

Finally, the 137-kb del(6)(q21q21) resulted in the fusion gene TRAF3IP2::REV3L, fusing exon 8 of TRAF3IP2 to exon 3 of REV3L. This fusion was not detected by any long-read callers but was found by three short-read fusion callers and was supported by six Illumina reads and 13 IsoSeq reads (Fig 5G).

Further splicing and read support details can be found in Table S6.

Of note, three of the seven confirmed fusion genes involve the partner genes ETV6 and RUNX1 or regions in their immediate vicinity (ETV6::RUNX1, PHAX::AC007450.2, and RUNX1::PRDM7). One homolog of chromosome 21 was found to be involved in both RUNX1 fusions, whereas both homologs of ETV6 were involved in different chromosomal aberrations occurring on chromosome 12, with one forming the ETV6::RUNX1 fusion and the other undergoing a 584-kb deletion between the breaking and fusion events of the t(5;12), deleting the entire gene and leaving no WT ETV6 in the genome.

The performance of the seven fusion gene calling tools varied widely. The number of fusion gene candidates found by the short-read callers had a wide range: STAR-fusion found only seven, Arriba found 31, and FusionCatcher found 78, while at the higher end of the range were Squid (n = 393) and pizzly (n = 5,520). The two long-read callers returned a large number of candidates: 336 from Cupcake and 4,927 from JAFFAL. The percentage of fusion gene candidates passing automated filtering were as follows: JAFFAL, 0.24%; Squid, 0.25%; pizzly, 0.36%; Cupcake, 1.79%; FusionCatcher, 7.69%; Arriba, 29.03%; and STAR-fusion, 57.14%.

The overall performance of the fusion callers also varied widely, as measured by sensitivity (percentage of the manually confirmed fusion genes detected) and FPR, calculated after discounting fusion gene candidates also found in the GM12878 RNA-seq dataset. Arriba detected each of the seven manually confirmed fusion genes (100.0% sensitivity, but 76.67% FPR), and STAR-fusion had the lowest FPR (50.0%) but only detected three fusion genes (42.86% sensitivity). Squid had the lowest sensitivity rate (14.29%; 1 confirmed gene fusion found), with an FPR of 99.65%. The remaining tools each called three of the confirmed fusion genes (42.86% sensitivity rate) and had an FPR between 94.23% and 99.94% (Table S7).

Key SNVs conclude the detailed characterization of the REH cell line

Finally, we characterized the SNVs in the REH genome. For SNV discovery, only the Illumina data were used, because of the higher coverage and lower error rate of this dataset in comparison with the long-read datasets (Table S2). Key SNVs, detailed below, were subsequently corroborated in the long-read data.

Mutect2 was run on the Illumina WGS dataset, filtered against a panel of normals from the 1000 Genomes Project and allele frequencies from gnomAD, to produce a set of SNVs that were further filtered to include only variants affected by loss of function or nonsense-mediated mRNA decay (NMD). The resulting set contained 131 variants (Table S8), including a nonsense mutation in exon 5 of GC receptor NR3C1 (p.Gln528Ter; Fig S11) and a frameshift mutation in exon 8 of PAX5 (p.A322Rfs*19; Fig S12). Phasing analysis placed the NR3C1 mutation on the same homolog of chromosome 5 (contig 56) as the t(5;12).

Combining cytogenetics, SV analysis, and SNV detection, our integrative short- and long-read genomic and transcriptomic analysis allows us to provide a detailed characterization of the REH cell line, shown in Fig 5.

REH cell line

The cell line REH was obtained from DSMZ (ACC 22) and cultured according to the supplier’s specifications. The cell line was verified by Eurofins human STR profiling cell line authentication services. DNA and RNA were extracted from 5 M cells using QIAGEN’s AllPrep DNA/RNA Micro Kit. For the long-read WGS libraries, PacBio’s Circulomics Nanobind UHMW DNA Extraction was used with an input of 15 M cells.

G-banded karyotyping was performed by the Department of Clinical Genetics, Uppsala University Hospital, Uppsala, Sweden. The cells from the cell culture were harvested following standard cytogenetic methods (KaryoMAX Colcemid Solution in HBSS, Gibco), and chromosome preparations were stained by G-banding trypsin–Giemsa standard procedures (0.25% trypsin, Gibco; Wright’s stain, Sigma-Aldrich; Giemsa’s stain, VWR Chemicals and Buffer tablets, pH 6.8, Merck). Twenty-five metaphases were analyzed using a Carl Zeiss scanner microscope and Metafer and Ikaros software from MetaSystems. The karyotype was carried out according to the International System for Human Cytogenetic Nomenclature (ISCN) 2020.

Sequencing

An Illumina PCR-free WGS library was generated with 1,000 ng input DNA and sequenced with PE-150 on a HiSeq X instrument. An Illumina short-read RNA-seq library was generated with 300 ng input total RNA and the TruSeq Stranded Total RNA kit, and sequenced with PE-100 on a NovaSeq 6000 instrument. One ONT WGS library was generated from 40 μg input DNA prepared with the Ultra-Long Sequencing Kit (SQK-ULK001) and sequenced on PromethION 24 using a R9.4.1 flow cell. Two PacBio IsoSeq RNA-seq datasets and one PacBio SMRTbell HiFi WGS dataset were prepared with SMRTbell Express Template Prep Kit 2.0, Sequel II Binding Kit 2.0, and Sequel II Sequencing Plate 2.0. For the IsoSeq library, 300 ng input RNA was used for each of the two libraries; a varying bead ratio was used to generate one library with standard-length transcripts and one library with full-length transcripts. Full-length non-chimeric (FLNC) reads were generated using PacBio tools ccs and lima. For the PacBio WGS, 10 μg input DNA was used to create a continuous long-read (CLR) library and 15 μg to create a HiFi library; the resulting CCS reads from both libraries were merged for analysis.

One RNA-seq dataset was also generated from GM12878 (B lymphoblastoid cell line) by merging four Illumina short-read RNA-seq libraries, each generated using 500 ng of input total RNA with the TruSeq Stranded Total RNA with Ribo-Zero depletion, and sequenced on a NovaSeq 6000 instrument.

Data analysis

Reads were aligned to the reference genome version GRCh38, followed by analysis with tools for SNV, SV, and fusion gene detection (Table 1). For gene annotation, GENCODE V37 (Ensembl 103) was used. Sequencing statistics (Table S2) were generated using Alfred. Versions for all software packages are listed in Table S9.

For visualization, we ran Copycat on the short- and long-read BAM files to bin the read coverage in preparation for visualization. Visualization was performed with SplitThreader (Nattestad et al, 2016 Preprint), Ribbon (Nattestad et al, 2021), and IGV. We used SplitThreader to visualize interchromosomal SVs across the Illumina, PacBio, and ONT callsets, whereas Ribbon was used for inspection of intrachromosomal features and split-read analysis. Circos plots were drawn using Circa. Heatmaps were generated with the chromoMap R package.

WGS data analysis

The short-read WGS data were mapped and processed with the nf-core/sarek pipeline (Garcia et al, 2020). PacBio CCS reads were mapped and sorted using pbmm2. After adapter sequences were trimmed using Porechop, the ONT reads were aligned using minimap2.

For the variant allele frequency plots (Figs 1 and S2), three sets of SNV calls were generated by running DeepVariant on the WGS datasets (in conjunction with PEPPER, for ONT).

For the filtered SNV callset (Table S8), Mutect2 was run on the Illumina alignments in tumor-only mode, filtered against a panel of normals generated from the 1000 Genomes Project and allele frequencies from gnomAD, both provided by the GATK Best Practices bundle. Calls with a population frequency of over 0.0001 were excluded. SNV calls were subsequently annotated with SnpEff and filtered for loss of function or NMD affecting at least 25% of the gene’s transcripts.

TIDDIT (Eisfeldt et al, 2017) was used for SV detection in the short-read data. In both long-read WGS datasets, Sniffles2 (Sedlazeck et al, 2018) was run with non-germline mode enabled. SURVIVOR was used to filter and merge VCF files using a minimum size of 100 bp and excluding variants found in the ENCODE Blacklist. For consensus callsets, SVs were required to agree on type and strand, with breakpoints at a maximum allowed distance of 1 kb.

For the evaluation of large-scale variants and translocations, the original callsets were pre-filtered, retaining SV candidates with either length > 100,000, or of a breakend type and mapping to two different chromosomes; discarded candidates were excluded from accuracy statistics. The remaining candidates were further filtered, requiring them to meet all the following criteria: (1) passing all quality filters set by SV calling software (2) a minimum of five supporting reads, with the number of supporting reads >20% of the average depth of coverage (DC) for the dataset and (3) position coverage no greater than 150% of the average DC for the dataset.

Of the remaining candidates, the following were selected for manual inspection: (1) all candidates with support in at least one Sniffles callset and (2) TIDDIT candidates with exactly 15 supporting reads, providing a pseudorandom sampling from this callset. SV candidates were confirmed visually in IGV, requiring the presence of split long reads and/or discordant mates in at least two of the three SV callsets for confirmation. Final candidates were screened against gnomAD SVs v2.1 to exclude possible germline variants.

The de novo assembly used for SV phasing was created with Flye and Medaka on ONT ultralong reads. Adapter sequences were trimmed from the ONT data before assembly using Porechop. We in addition trimmed the first 50 bp from all reads and removed all reads shorter than 500 bp using NanoFilt (De Coster et al, 2018). After assembly using Flye, polishing was performed with Medaka, which uses neural networks and is trained on specific configurations of ONT pore type, instrument, and basecaller. This assembly was polished using the Medaka model r941_prom_hac_g507. The quality of the assemblies was assessed with QUAST (Mikheenko et al, 2018). The assembly was phased using Hapdup, then realigned to the reference genome, and subjected to SV calling using Hapdiff.

RNA-seq analysis

The nf-core/rnafusion pipeline was run on short-read REH RNA-seq data to call and visualize putative fusion genes. This pipeline uses five different fusion tools: Arriba, FusionCatcher, pizzly, Squid, and STAR-fusion. The nf-core/rnafusion pipeline was also run on the GM12878 short-read RNA-seq dataset to generate a set of non-leukemic fusion gene candidates for filtering.

Putative fusion genes were called in the long-read IsoSeq data with two different workflows. First, cDNA_Cupcake was used in conjunction with clustering using isoseq3 and fusion gene classification using SQANTI3. Second, the JAFFAL pipeline was run with default settings.

Candidate fusion genes were programmatically pre-filtered as follows. All candidate fusion events called by rnafusion in the GM12878 cell line were removed and excluded from accuracy statistics. The remaining candidate fusion genes were required to have different partner genes and to fulfill any of the following criteria: (1) containing one or more known ALL-associated genes (Marincevic-Zuniga et al, 2017), and supported by at least five reads (2) called by at least one short-read fusion detection tool and one long-read fusion detection tool and (3) called by at least three short-read fusion detection tools, or (4) supported by at least 10 long reads. The remaining candidates were inspected manually in IGV.