We used two types of cell lines and four different datasets to evaluate the performance of fusion gene detection (see Table 1). We selected the Michigan Cancer Foundation-7 (MCF7) and SKBR3 cell lines because their long-read RNA-seq datasets have been commonly used as a benchmark for fusion-gene detection, and a previously validated list of fusion genes, considered ground truth, is available [15].

The reference genome and transcriptome sequence dataset we used was the GENCODE Release version 19 (https://www.gencodegenes.org/human/release_19.html) [16]. We used JAFFAL version 2.2 (https://github.com/Oshlack/JAFFA) and FusionSeeker version 1.0.1 (https://github.com/Maggi-Chen/FusionSeeker) to compare performance with our tool. We also used minimap2 version 2.17 (https://github.com/lh3/minimap2) to align long reads to the reference genome and transcriptome and BLAT (BLAST-Like Alignment Tool) version 36x2 (https://github.com/djhshih/blat) to align gap regions to the reference genome. JAFFAL and FusionSeeker also used the GENCODE dataset as their references. The annotation of the gene set was carried out on genome assembly GRCh37 (hg19) available at: https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001405.13/.

JAFFAL and FusionSeeker have reported superior accuracy compared to LongGF [14, 15]. Thus, we focused on these two tools for performance comparison. JAFFAL and FusionSeeker face two main challenges: (1) detecting the correct fusion genes, and (2) identifying the correct breakpoints (as shown in Fig. 2).

JAFFAL aligns long reads to both the reference transcript and genome using minimap2 [17] and filters both alignment results by applying constraints. When certain regions of the read near the breakpoint (gap regions as shown in Fig. 2A) fail to align with the genome, possibly due to base-calling errors, it can lead to incorrect alignment with unrelated genes because of partially aligned regions (Fig. 2B), or inaccurate identification of breakpoint positions caused by the gap regions (Fig. 2C). To solve the problem of breakpoint position identification, JAFFAL compares the alignments against annotated genes in the reference genome and realigns breakpoints to exon boundaries [11]. JAFFAL also employs a restriction rule that excludes reads with large gap regions between aligned multiple genes from the analysis. However, the length of the gap region tends to be large in some gene fusion events (such as cis-splicing, as shown in Fig. 1B). In those cases, JAFFAL may not be able to detect fusion genes with large gap regions. Furthermore, either one read of a fusion gene pair can be a short fragment sequence. While minimap2, the alignment tool used in JAFFAL, is robust against sequence noise, it has limitations in accurately aligning short fragment sequences to the genome.

FusionSeeker is a tool that utilizes density-based clustering to accurately identify breakpoints that may contain misalignments. The clustering-based refining process is highly effective, especially given the challenges of achieving precise alignment with the reference genome in long-read RNA-seq. FusionSeeker primarily focuses on identifying isoforms of detected fusion genes and requires a minimum of three supporting reads to generate a consensus sequence.

Previously, we proposed FLBEA (Full-Length and Both-Ends Alignment) to address these challenges [18]. This tool generates fragment sequences from both ends of long reads and aligns them to the genome. The long reads are identified as fusion genes when their fragment sequences are aligned with two different genes. This approach improves precision by incorporating fragment sequences, as relying solely on long reads for fusion gene detection often results in false positives. It is important to note that FLBEA does not directly address the issue of gap regions shown in Fig. 2A.

FLBEA was able to detect fusion genes with higher performance than JAFFAL when the base-calling error rate was low in long reads [18]. However, when the error rate was high (more than 10%), the fragment sequences were not precisely aligned to the reference genome because the fragment sequences were shorter than the original long reads and the accuracy of the alignments was more influenced by the base-calling errors.

We recently proposed a fusion gene detection tool, FUGAREC (Fusion Detection with Gap Re-alignment) version 1.0 (https://github.com/Hideo-Matsuda/FUGAREC) [19]. FUGAREC aims to solve the challenges shown in Fig. 2 by the following three steps: (1) Anchors breakpoints to exon boundaries (Step 1 in Fig. 3A); (2) To align gap regions that cannot be anchored to exon boundaries (Step 2 in Fig. 3A), we use BLAT [20], instead of minimap2; and (3) Reads are clustered every 1000 bp based on the genomic position of their breakpoints. When multiple fusion reads are within a cluster, the read with the highest number of supporting reads is selected (Fig. 3B).

In Step (1), the breakpoint anchoring is an extension of the JAFFAL method [14] to solve the problem of the gap regions (Fig. 2A). While JAFFAL anchors breakpoints to the exon boundaries only when the breakpoints are located within 50 bp of the exon boundaries, FUGAREC does not set any distance limit from the exon boundaries. Instead, it selects only the nearest boundary positions for anchoring breakpoints. Thus, FUGAREC may be able to detect fusion genes with large gap regions. FusionSeeker does not anchor breakpoints to exon boundaries, but it employs a different approach as described in Step (3).

In Step (2), to solve the problem of detecting the incorrect furoin pair (Fig. 2B), FUGAREC adopts BLAT to re-align gap regions to the reference genome. Since BLAT is specifically designed for aligning short sequences to the reference genome and only FUGAREC realigns the gap regions using BLAT, it may more accurately identify the correct fusion pairs compared to JAFFAL and FusionSeeker.

Step (3) aims to solve the problem of incorrect breakpoints (Fig. 2C). JAFFAL also performs clustering of reads on the genomic position of their breakpoints, which will be either the one preserving exon boundaries, or the one with the highest read support with their breakpoints ranging within 50 bp. FUGAREC achieves clustering of reads within wider range intervals, 1000 bp. In addition, FusionSeeker employs a different approach that performs clustering of supporting reads and a partial order alignment of the reads using bsalign [21] to generate a consensus sequence for each fusion gene. Although the approach can contribute to identifying correct breakpoints, it requires a minimum of three supporting reads. As a result, FusionSeeker may not detect fusion genes with low expression levels. FUGAREC does not create consensus sequences and thus can detect fusion genes even with one or two supporting reads.

Since FUGAREC always outperformed FLBEA in performance comparison using simulated data [19], we do not use FLBEA for further analysis in this paper. An earlier version of this paper [22] reported preliminary results of the fusion gene detection from the same long-read RNA-seq datasets shown in Table 1 using FUGAREC. In this paper, we have largely extended the paper [22] by adding the mechanism of gene fusion formation and discussing the possibility of detection of novel fusion genes.

Fusion genes supported by two or more reads were considered for detection. These were evaluated as true positives (TP) if the detected genes appeared in the lists of the previously validated gene fusions [15], and as false positives (FP) otherwise. False negatives (FN) are known fusion genes that were not detected. Subsequently, the performance of fusion gene detection was evaluated based on precision, recall, and F1 score are described as precision = TP/(TP + FP), recall = TP/(TP + FN), and F1 = 2 $\times$ precision $\times$ recall/(precision + recall), respectively. The F1 score, the harmonic mean of precision and recall values, serves as a comprehensive evaluation metric. Therefore, the F1 score was used to evaluate the performance of each tool in detecting fusion genes.

To evaluate the performance of our tool in detecting fusion genes, we compared it with JAFFAL and FusionSeeker using the same data. JAFFAL was executed using the configuration file JAFFAL.groovy. FusionSeeker was run with the -datatype nanopore option for Nanopore sequencing data and the -datatype iso-seq option for PacBio sequencing data. FUGAREC was run using minimap2 version 2.17 with the -x map-ont option and using BLAT version 36x2 with the options: -out=blast8 -minScore=10 -stepSize=5.

We evaluated the performance of FUGAREC compared to those of JAFFAL and FusionSeeker as shown in Tables 2,3. FUGAREC achieved the highest precision among the three tools across all long-read datasets. FusionSeeker had the highest recall in the SKBR3 dataset, but it was lower than FUGAREC in the other datasets. FUGAREC demonstrated the highest F1 score across all datasets.

The overlap in the list of fusion genes detected by each tool is shown in Fig. 4. FUGAREC uniquely detected 4, 9, 67, and 0 fusion genes in the MCF7-ONT-1, MCF7-ONT-2, MCF7-Pacific Bioscience (Pac), and SKBR3-Pac datasets, respectively. Of these, 1, 2, 0, and 0 genes, respectively, matched the previously validated fusion genes [15].

We assume that the gene pairs detected by all three tools (JAFFAL, FusionSeeker, and FUGAREC) are fusion gene candidates. The number of candidates detected by the three tools and not included in the previously validated fusion-gene list is 9, 6, 15, and 0 for the four datasets, respectively. The numbers varied widely across the datasets, and it might not be a good criterion for novel fusion genes due to the difference between the three tools.

Additionally, we counted the number of gene pairs detected in at least two of the four datasets using FUGAREC as fusion gene candidates. The candidates that were not included in the known fusion gene list are presented in Table 4.