==== Front Bioinformatics Bioinformatics bioinformatics Bioinformatics 1367-4803 1367-4811 Oxford University Press 37387174 10.1093/bioinformatics/btad264 btad264 Genome Sequence Analysis AcademicSubjects/SCI01060 isONform: reference-free transcriptome reconstruction from Oxford Nanopore data Petri Alexander J Department of Mathematics, Science for Life Laboratory, Stockholm University, Stockholm 106 91, Sweden https://orcid.org/0000-0001-7378-2320 Sahlin Kristoffer Department of Mathematics, Science for Life Laboratory, Stockholm University, Stockholm 106 91, Sweden Corresponding author. Department of Mathematics, Science for Life Laboratory, Stockholm University, Stockholm 106 91, Sweden. E-mail: ksahlin@math.su.se 6 2023 30 6 2023 30 6 2023 39 Suppl 1 ISMB/ECCB 2023 Proceedings i222i231 © The Author(s) 2023. Published by Oxford University Press. 2023 https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract Motivation With advances in long-read transcriptome sequencing, we can now fully sequence transcripts, which greatly improves our ability to study transcription processes. A popular long-read transcriptome sequencing technique is Oxford Nanopore Technologies (ONT), which through its cost-effective sequencing and high throughput, has the potential to characterize the transcriptome in a cell. However, due to transcript variability and sequencing errors, long cDNA reads need substantial bioinformatic processing to produce a set of isoform predictions from the reads. Several genome and annotation-based methods exist to produce transcript predictions. However, such methods require high-quality genomes and annotations and are limited by the accuracy of long-read splice aligners. In addition, gene families with high heterogeneity may not be well represented by a reference genome and would benefit from reference-free analysis. Reference-free methods to predict transcripts from ONT, such as RATTLE, exist, but their sensitivity is not comparable to reference-based approaches. Results We present isONform, a high-sensitivity algorithm to construct isoforms from ONT cDNA sequencing data. The algorithm is based on iterative bubble popping on gene graphs built from fuzzy seeds from the reads. Using simulated, synthetic, and biological ONT cDNA data, we show that isONform has substantially higher sensitivity than RATTLE albeit with some loss in precision. On biological data, we show that isONform’s predictions have substantially higher consistency with the annotation-based method StringTie2 compared with RATTLE. We believe isONform can be used both for isoform construction for organisms without well-annotated genomes and as an orthogonal method to verify predictions of reference-based methods. Availability and implementation https://github.com/aljpetri/isONform Swedish Research Council 10.13039/501100004359 2021-04000 ==== Body pmc1 Introduction Long-read sequencing techniques such as Oxford Nanopore Technologies (ONT) and Pacific Biosciences (PacBio) can produce reads covering the majority of transcripts from end-to-end. Such sequencing techniques have proven to enable a better understanding of the transcriptional landscape of cells (Bayega et al. 2018; Byrne et al. 2019; Cole et al. 2020). However, while long-read methods omit the need for transcript assembly as needed with short reads, transcript 3′ and 5′ variability, RNA degradation, long-read sequencing error profiles, and other sequencing artifacts inhibit downstream analysis of long-read data. Particularly, recovering the actual isoforms from long-read transcriptomic datasets has proven difficult, with studies observing thousands of low-quality and potentially spurious transcript predictions (Hoang et al. 2017; Kuo et al. 2020). 1.1 Genome- and annotation-based methods Due to the noise in these datasets, a common approach to predict transcripts is to align the transcriptomic reads to a reference genome using a long-read splice aligner such as minimap2 (Li 2018), deSALT (Liu et al. 2019), or uLTRA (Sahlin and Mäkinen 2021). Several tools for reconstructing transcripts from long-read splice alignments have been proposed (Kovaka et al. 2019; Tung et al. 2019; Tang et al. 2020; Kuo et al. 2020; Wyman et al. 2020; Holmqvist et al. 2021; Orabi et al. 2022; Chen et al. 2022; Volden et al. 2022; Prjibelski et al. 2023). These tools typically predict transcripts by requiring both alignments and transcript annotations (Kovaka et al. 2019; Tung et al. 2019; Tang et al. 2020; Kuo et al. 2020; Wyman et al. 2020; Holmqvist et al. 2021; Chen et al. 2022; Volden et al. 2022), or requiring alignments but with the ability to predict transcripts outside annotations (Orabi et al. 2022; Prjibelski et al. 2023). A thorough community-effort benchmarking of many of these tools using several different sequencing techniques and genomes has recently been performed as part of the LRGASP challenge (Pardo-Palacios et al. 2021), and preliminary results of challenge 1 (LRGASP 2022) show a large discrepancy between methods in isoform detection, even for organisms with high-quality reference genomes such as mouse and human. Larger discrepancies in predicted transcripts between methods can be explained by the different algorithms they use, as well as error levels and artifacts of the sequencing techniques. However, there are also inherent limitations with predicting transcripts from alignments to a linear reference genome. Such limitations include predicting transcripts from fused genes, gene copies, or exons that are not present in the linear genome or shortcomings in splice alignment, which the predictions fundamentally rely on. For example, in the study by Sahlin and Mäkinen (2021), we showed that long-read splice alignment methods have limitations such as (i) small exons not being detected by read aligners, (ii) overfitting to canonical splice sites, and (iii) not aligning over long introns. Additionally, relying on a reference genome and reference annotation limits the methods to organisms and genes for which reference annotations exist. 1.2 Genome- and annotation-agnostic methods By not relying on high quality genome assemblies and annotations, reference-free approaches are important for research of lesser known organisms/gene families. In addition, they can overcome the splice-alignment limitations and the bias introduced by linear reference genomes. So far, relatively few methods have been proposed for reference-free transcript predictions. Some tools have been proposed for PacBio sequencing such as isoseq3 (Gordon et al. 2015) and IsoCon (Sahlin et al. 2018) but are limited to PacBio data and IsoCon is only applicable to targeted sequencing. While there exists recent unpublished work (Nip et al. 2022), the only published tool for reference-free transcriptome reconstruction of ONT reads is RATTLE (de la Rubia et al. 2022). However, as we show in our study, RATTLE has low transcript prediction sensitivity and misses many transcripts present in data. Some challenges of reference-free transcript reconstruction are highly variable abundance, alternative splicing, as well as long-read error profiles. 1.3 Our contribution We have previously demonstrated algorithms for clustering and error correction of long transcriptomic reads through isONclust (Sahlin and Medvedev 2020) and isONcorrect (Sahlin and Medvedev 2021). We here introduce an algorithm and its implementation, isONform, for the final step of this pipeline to produce transcript predictions. IsONform builds a graph based on paired-minimizer seeds, that, similarly to minimizer-based genome assembly approaches (Rautiainen and Marschall 2021; Ekim et al. 2021), aims to sparsely represent the sequence information in the data. IsONform then employs iterative bubble popping on the graph to remove errors while keeping exon differences. When no poppable bubbles exist, isoform predictions can be traced from the graph by following paths of full-length reads. A fundamental difference between our isON pipeline (isONclust, isONcorrect, and isONform) and RATTLE’s pipeline (cluster, error correct, and polish) is that the isON pipeline works with clusters of reads from a gene, while RATTLE clusters by transcript. This allows isON’s error correction and consensus prediction (isONform) to leverage shared exons among splice variants, thus, the potential to recover more transcripts with higher base-level accuracy compared with predictions by RATTLE. However, to benefit from shared exons in different transcripts, it also requires a more sophisticated consensus generation algorithm able to separate isoforms at the consensus calling step, compared with RATTLE, which only does base-level polishing per transcript in the consensus (polishing) step. Using simulated and synthetic reads, we show that isONform recovers substantially more isoforms than RATTLE. On biological data, we show that isONform has much higher consistency with the reference-based tool StringTie2 compared with RATTLE. We believe that isONform (integrated into our isON pipeline) has two primary use cases. First, it can be used for isoform construction for organisms without well-annotated genomes. Second, it can be used as an orthogonal method to verify consistency in reference-based prediction since systematic limitations exist, as described in Section 1.1. 2 Materials and methods In Section 2.1, we give a high level overview of the algorithm. In Section 2.2, we introduce notation and formally define minimizer pairs and the graph structure that isONform uses. Finally, we describe the algorithm in depth in Section 2.3 and implementation details in Section 2.4. 2.1 Algorithm overview The input of our algorithm are clustered and error-corrected reads from isONclust (Sahlin and Medvedev 2020) and isONcorrect (Sahlin and Medvedev 2021). Each cluster represents reads from isoforms from a gene family and can be processed individually and in parallel with isONform. The output is a set of predicted transcripts in fastq format. We describe the algorithm for reads from a single cluster. In step 1, reads are partitioned into intervals using “minimizer pairs” (Sahlin and Medvedev 2021) as seeds. In step 2, a weighted interval scheduling (WIS) problem is solved to find a set of maximally weighted and nonoverlapping minimizer pairs (NOMPs) for each read as done in isONcorrect (Sahlin and Medvedev 2021). In step 3, the set of NOMPs is used to construct a directed acyclic graph (DAG) in which vertices represent the NOMPs and neighboring NOMPs are connected with an edge. In step 4, the graph is simplified by an iterative bubble-popping algorithm designed to remove errors and mutations between reads while keeping exon differences. In step 5, predicted isoforms are generated by finding supported paths through the graph. Finally, step 6 removes any redundant isoform predictions caused mainly by the larger 3′ and 5′ variability using pairwise alignment of predicted isoforms. The isoform predictions remaining after this step constitute the final predictions. An overview of the algorithm is shown in Fig. 1. Figure 1. Overview of the isONform algorithm. IsONform takes as input clustered and error-corrected reads generated with isONclust and isONcorrect, respectively. The example illustrates a cluster consisting of five reads (r1 to r5) from three isoforms, where p denotes positions of minimizers on r1. isONform uses all minimizer pairs from minimizers at distance xmin to xmax. The figure shows an example of 10 minimizer pairs on r2. Each minimizer pair is assigned a weight formed from the length of the minimizer pairs and its abundance in the cluster. The instance is then sent to a weighted interval scheduler, which is used to find the set of NOMPs. In our example, we find five segments for r1 while the scheduler finds four instances for r2. A DAG is built up by the NOMPs, and we add a global source node s as well as a global sink node t (step 3). An iterative bubble-popping approach is employed to merge similar paths in the graph (step 4). This reduces the number of paths supported by reads (step 5). We finally perform postmerging of the transcripts due to 3′ and 5′ variability. For each isoform, a consensus sequence is generated from all reads appointed to it using the SPOA algorithm. 2.2 Preliminaries isONform uses the same generation of minimizers and minimizer-pairs as described in isONcorrect (Sahlin and Medvedev 2021). We describe them here for completeness. 2.2.1 Minimizers Minimizers (Roberts et al. 2004) are subsampled k-mers. Let r∈R denote a string of nucleotides, which we refer to as “read.” We use r[i] to refer to the nucleotide located at the i-th position of r. Given two integers k and w such that 1≤k≤w≤|r|, the minimizer of r at position p is the lexicographically smallest k-mer m that starts inside the interval [p,p+w). We then say that m is a minimizer of the read r, or, alternatively, is a positional minimizer (m,p) of read r. Let M(r) be the set of positional minimizers of read r. In the following example, we use r=ACGGATCAC,k=2,w=4, yielding the set of positional minimizers M(r)=(AC,0),(CG,1),(AT,4),(AC,7). 2.2.2 Minimizer pairs Let xmin and xmax be two positive integer parameters. We let Wr={((mi,pi),(mj,pj))∈M(r)×M(r)|xmin≤pj−pi≤xmax} be the ordered set (according to increasing pi then pj) of paired positional minimizers separated by at least xmin and at most xmax nucleotides in r. Similarly, we denote the sequence of paired minimizers as StrWr={(mi,mj)|((mi,pi),(mj,pj))∈Wr}, i.e. Wr with the positions omitted but duplicates retained. The above example with the following parameters xmin = 2, xmax = 4 therefore yields Wr=((AC,0),(AT,4)),((CG,1),(AT,4)),((AT,4),(AC,7)) and strWr=(AC,AT),(CG,AT),(AT,AC). Given a set of reads R, we let W be the union of all Wr for the reads in R and we let StrW be the union of all StrWr. Besides error correction (Sahlin and Medvedev 2021), minimizer pairs have also been used in genome assembly (Chin and Khalak 2019). 2.2.3 Nonoverlapping minimizer pairs For each read, isONform, similarly to isONcorrect, produces a subset of NOMPs from the set of minimizer pairs. The NOMPs are decided upon by the solution to a WIS problem. Specifically, the minimizer pairs span an “interval” on the read. The input to the WIS problem is a set of minimizer pairs (intervals) I=I1,…In, where Ij∈[aj,bj],aj,bj∈R and aj < bj. The number of reads supporting an interval Ij (provided in StrW) gives a weight wj to Ij. The solution to the WIS problem outputs a subset I′⊂I of nonoverlapping intervals for which the sum of weights is maximized. Our intervals in the solution I′ correspond to NOMPs that we index in the order they appear on the read, e.g. NOMPjr is the jth NOMP in r. The WIS problem can be solved exactly by applying a dynamic programming algorithm that runs with O(n log n) time complexity, with n denoting the number of intervals (Kleinberg and Tardos 2006). 2.2.4 Graphs of nonoverlapping minimizer pairs We represent the NOMPs as vertices, where vir denotes the i-th NOMP on read r. We draw a directed edge eijr between two vertices vir and vjr if vir and vjr represent neighboring NOMPs on a read. This forms a directed graph G=(V,E), where we insert a global source and sink node s and t, respectively. We draw a directed edge (s,v′) from the global source to the first NOMP v′ on each read. Similarly, we draw a directed edge (v″,t) from the last minimizer pair v″ on each read to the global sink. A read spells a path from the global source to the sink. Each node vir also contains information on the start and end position and the length (span) of the NOMP on read r. Furthermore, each edge ei,i+1r has a length l(ei,i+1r) associated with it, representing the distance between the i-th and (i + 1)-th NOMP on r. Therefore, a path p in G that follows a read r has a corresponding length (|pr|), start, and stop position on r. The graph is built up by iterating through the reads in order of appearance in the fastq file. A NOMP that is shared between reads will be represented by a single node. A NOMP is considered shared if it has the same minimizer pair and they do not have a difference in their pairwise distance larger than δ (parameter to isONform). However, isONform stores the length of each NOMP in the vertex. A detailed description of the graph construction is given in Section 2.3.1. 2.2.5 Bubbles and multibubbles We define a bubble as a pair of paths with a shared start node bs and end node be where the internal path nodes are disjoint. Note that nodes that are not part of the bubble can be reached from nodes belonging to the bubble. If the bubble additionally has at least one read supporting each path from bs to be we denote it as “read-supported bubble.” We will detect and resolve what we refer to as read-supported multibubbles. An example of a multibubble structure is shown in Fig. 2a. Multibubbles have similar characteristics as superbubbles (Onodera et al. 2013) except that multibubbles allow nodes not part of the multibubble to be reachable from nodes in the multibubble and for several bubbles to start/end in one node, i.e. they do not fulfil the minimality and the matching criteria (Onodera et al. 2013). Let G=(V,E) a directed graph and (s, t) being an unordered pair of distinct vertices. Then, a multibubble is defined if it satisfies the following criteria: Figure 2. Bubble popping on a graph. (a) Our initial graph. There are three bubbles in this graph ({a, b, c, d},{s’,a, c, d, g, h, i, e}, and {s’,a, b, d, e, f, g, h, i, j, k, l, m, t’}). The bubble popping starts with the smallest bubble (fewest vertices) and continues until no more bubbles can be popped. (b) Bubble {a, b, c, d} was popped and its nodes linearized. (c) Bubble {s’,a, b, c, d, g, h, i, e} was popped and its nodes linearized. (d) The bubble in (c) could also be popped and we now have a single path from s′ to t′, representing a single isoform in the graph. Reachability: t can be reached from s. Acyclicity: the subgraph induced by U is acyclic where U is the set of vertices as described above. then we say that the subgraph in the description of the acyclicity condition forms a multibubble and s, t and U∖{s,t} are the multibubble’s source, sink and interior respectively. 2.3 Algorithm Steps 1 and 2 consist of constructing minimizer pairs and solving the WIS problem to produce NOMPs as input to isONform. These steps are performed identically to isONcorrect. 2.3.1 Step 3: graph construction The graph is built up by iterating through the reads in order of appearance in the fastq file, and adding vertices and edges from the NOMPs produced in step 2. We want to keep the graph acyclic. Repetitive NOMPs within a read may cause a cycle. Also, since NOMPs are fuzzy seeds, they may have different lengths. Therefore, we want to identify larger length differences between NOMPs shared between reads. Finally, there can be gaps between neighboring NOMPs in reads. If two reads share two neighboring NOMPs, the gap between the two NOMPs may differ, indicating a structural difference between the reads. We handle those cases as follows. Preserving acyclicity: For each new NOMPi+1r to be added to G, if NOMPi+1r is already found in the graph it may create a cycle if we add the edge ei,i+1r connecting it to its predecessor. We can prevent G becoming cyclic by looking at the topological order of the graph. If ei,i+1r would introduce a cycle, we instead add a new vertex vi+1r to G representing NOMPi+1r. Differences in NOMP span: For each NOMPi+1r to be added to G, if the NOMP is already found in the graph, we check the length difference of NOMPi+1r and the NOMP of the node existing in the graph. If the NOMP lengths differ by more than δ nucleotides, we add a new node. If a NOMP is already represented by two or more NOMPs, NOMPi+1r is added to the vertex that has the NOMP with the closest distance to NOMPi+1r. Differences between NOMPs: In a similar fashion to checking differences in NOMP span, if we find two reads r and r′ having the same neighboring NOMPs (NOMPir,NOMPi+1r) = (NOMPjr′,NOMPj+1r′) but |l(ei,i+1r)−l(ej,j+1r′)|>δ then we create a separate vertex vi+1r′ for the second NOMP while adding the first NOMP to the vertex already in the graph. 2.3.2 Step 4: iterative bubble popping NOMPs are indel tolerant seeds, but sequencing errors and exon differences may result in local bubbles (Zerbino and Birney 2008) in G. It requires only a few bubbles to lead to many reads having unique paths through G (possibly exponentially increasing with the number of bubbles), which yields redundant consensus sequences. The purpose of the iterative bubble-popping step is to remove local bubbles in the graph caused by read errors and SNP differences while preserving differences at the exon level between reads. Overview: At a high level, G may contain several multibubbles sharing vertices. We find disjoint “read-supported bubbles” that we can pop during each iteration of our iterative bubble popping. We start by analyzing the smallest bubbles (with respect to the number of vertices) and gradually test and pop bubbles until no more new bubbles can be detected in the graph. We test if a bubble is poppable by forming a consensus (using SPOA, Vaser et al. 2017) of the subsequences of reads corresponding to each path in the bubble. If such consensuses do not differ in length by more than δ and have sequence similarity higher than α (parameter to isONform), we pop the bubble. Otherwise, we save it as unpoppable. If a bubble can be popped, we use the relative distance information to reorder (or “linearize”) the vertices of the two paths. At each new iteration, the bubbles we previously investigated could either be popped or were stored as unpoppable. Bubbles that have been deemed to be unpoppable are skipped in the following iterations. During each iteration, the bubble-popping approach reduces the number of edges in our graph. A graphical overview of this process is illustrated in Fig. 2. Details: In detail, the process of iterative bubble popping is performed as follows: We identify all vertices having more than one out-edge S′ as well as all vertices T′ that have more than one in-edge. These vertices are potential source and sink vertices of bubbles. We then generate all potential pairs of (s′,t′), s′∈S′ and t′∈T′ where TOP(s′)