{'fraction_non_alphanumeric': 0.03733429275177815, 'fraction_numerical': 0.008634766003376414, 'mean_word_length': 4.485653560042508, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SPATA: A Seeding and Patching Algorithm for Hybrid Transcriptome Assembly

Transcriptome assembly from RNA-Seq reads is an active area of bioinformatics research. The ever-declining cost and the increasing depth of RNA-Seq have provided unprecedented opportunities to better identify expressed transcripts. However, the nonlinear transcript structures and the ultra-high throughput of RNA-Seq reads pose significant algorithmic and computational challenges to the existing transcriptome assembly approaches, either reference-guided or de novo. While reference-guided approaches offer good sensitivity, they rely on alignment results of the splice-aware aligners and are thus unsuitable for species with incomplete reference genomes. In contrast, de novo approaches do not depend on the reference genome but face a computational daunting task derived from the complexity of the graph built for the whole transcriptome. In response to these challenges, we present a hybrid approach to exploit an incomplete reference genome without relying on splice-aware aligners. We have designed a split-and-align procedure to efficiently localize the reads to individual genomic loci, which is followed by an accurate de novo assembly to assemble reads falling into each locus. Using extensive simulation data, we demonstrate a high accuracy and precision in transcriptome reconstruction by comparing to selected transcriptome assembly tools. Our method is implemented in assemblySAM, a GUI software freely available at http://sammate.sourceforge.net.


Introduction
Human transcriptomes are highly diverse, overlapping, complex, and dynamic. Alternative splicing and structural variations play important roles in enhancing the range of transcriptome complexity [1,2,3,4,5,6,7,8,9]. For example, it is reported that over 90% of human genes are alternatively spliced, and up to 5% of structural variations, such as insertions and deletions are present within exons [10,11,12]. Moreover, all these events are condition-specific that lead to diversity of human transcriptomes. Identifying the expressed transcript sequences is a central task in transcrip-tomics research since it provides critical information for further analysis of the transcriptome. The ever-declining cost and increasing depth of RNA-Seq provide unprecedented opportunities to better identify expressed transcripts [13]. Current efforts to reconstruct the expressed transcripts from short RNA-Seq reads generally follow one of the two strategies: an ab initio strategy, which is also called referenceguided, and a de novo strategy [14,15,16]. Despite initial success, both strategies are facing numerous challenges as described below and there is still a need for new transcriptome assembly approaches.
The reference-guided approaches offer good sensitivity, i.e., they are able to reconstruct transcripts with low abundance. However, since they rely on the alignment to the reference genome, their performance can be compromised by a number of factors. One major factor is that the reference genome being used maybe incomplete and may differ for each individual due to ubiquitous gene structural variations, such as point mutations, insertions, deletions, and gene fusions [4,5,7,29]. Moreover, errors and biases introduced by the splice-aware aligners are carried over to the assembled transcripts. For example, junction reads spanning large introns can be missed because some aligners restrict the intron length due to the computational complexity [23,25,27,15]. Consequently, the assembly programs may fail to reconstruct the corresponding transcripts.
In contrast, the de novo transcriptome assembly approaches generally assemble reads using de Bruijn graph [30,31], whose nodes are all subsequences of length k (k-mers) extracted from the short reads whereas two nodes are connected by an directed edge if they perfectly overlap by k − 1 nucleotides. The graph is then refined to eliminate possible false branches, using the original read sequences, the paired-end information and the coverage levels. The remaining paths are then traversed and reported as transcripts. This elegant computational solution was first introduced for whole genome assembly [32,33,34,35,36], where DNA sequencing depth is expected to be the same across the genome.
The de novo transcriptome assembly programs, such as Rnnotator [37], Trans-Abyss [38,39], Multiple-k [40], Trinity [41], and Oases [42], don't depend on the reference genome and are able to reconstruct transcripts originated from genomes that have undergo rearrangements. However, a straightforward application of de Bruijn graph on transcriptome assembly faces several challenges. First, the coverage depth fluctuates among transcripts, complicating the task of optimizing the trade-off between sensitivity and graph complexity [15,16]. Second, the computing resources required to de novo assemble the large transcriptomes as a whole may be overwhelming. Unlike the reference-guided transcriptome assembly programs, the de novo transcriptome assembly programs do not use a divide-and-conquer strategy to divide the large assembly problem into many smaller assembly problems.
Here we present SPATA, a hybrid approach that combines the strengths of the previously mentioned strategies while avoiding their pitfalls. Comparing with the existing approaches, we claim the following original contributions: (1) We localize as opposed to align reads to individual genomic loci. (2) We develop a novel seeding and patching algorithm to de novo assemble read sequences falling into each genomic locus. (3) We implement the multi-platform software assemblySAM, which takes RNA-Seq data in the FASTA or FASTQ formats as input and outputs transcript sequences.

Methods
SPATA is a novel hybrid approach that seeks to strike a good balance between the de novo and the reference-guided transcriptome assembly approaches. On one hand, the difference between SPATA and the de novo assembly approaches is that SPATA uses the existing reference genome to divide the large assembly problem into many smaller assembly problems and solve them independently. Since the graphs constructed from reads for each genomic locus are much less complex than that of the whole transcriptome, our approach promises a better performance and demands less computational resources. On the other hand, the difference between SPATA and the reference-guided assembly approaches is that SPATA does not rely on the exact mapping of the reads to the reference genome. While the reference guided approaches assemble the transcripts from a graph built from the mapped genomic positions, SPATA as-sembles read sequences falling into each genomic locus via a novel assembly algorithm.
Our approach proceeds into the following three stages: (1) Localization stage: we first align exonic reads to the reference genome. The remaining reads are then split and partially aligned to the reference genome. Reads that are aligned or partially aligned to a genomic locus are then used as input of the second stage. (2) Seeding and growing stage: we fully grow the backbone sequences starting from arbitrarily selected "seeds". The fully grown sequences are then used as input of the third stage. (3) Patching and cutting stage: we connect the backbone sequences to form an isoform graph and then traverse the paths of the graph to report the possible transcript sequences.

Stage I: Localization
We design a split-align-anchor procedure to localize reads into genomic loci. We proceed into two steps as shown in Figure 1.
Step 1: Initial alignment of all the reads. We first align the entire set of reads, single-end or paired-end, to the reference genome using a genome aligner, e.g., Bowtie 2 [43,44]. For single-end reads, this step initially aligns most of the human exonic reads to the reference genome. For paired-end reads, we localize the reads if at least one read of the pair is aligned to the reference genome, including reads with structural variations.
Step 2: Follow-up localization of the split reads. For the remaining reads from step 1, if they are single-end, we equally split each read into three partitions and localize the split reads using anchors. Anchors refer to any partition read(s) that are aligned to the reference genome. For paired-end reads, we equally split each read into two partitions, four partitions for a read pair, and localized them using anchors.
The main difference between our localization procedure and splice-aware aligners is that our procedure does not attemp to map every single nucleotide to the reference genome, which can be computational intensive or even impossible for reads with insertions or deletions. Instead, our goal is to localize reads to genomic loci. Reads falling into each genomic locus will serve as input of our de novo assembly algorithm, which will be described in sections 2.2 and 2.3. Below is the mathematical description of our localization procedure.

Definitions
Definition 1 (Short read). A short read R is a sequence of l letters, (r 1 r 2 ...r l ), where r i ∈ {A, T, C, G} for 1 ≤ i ≤ l. For convenience of the algorithm description, we denote any sequence of letters taken from {A, T, C, G} in the same manner. We denote the length of a sequence as |R|. We also denote an empty sequence as (). A paired-end read of two short reads R 1 and R 2 is denoted as (R 1 , R 2 ).
Definition 2 (Split). Given a short read R, we define the function Split(R, 2) that splits the read into 2 subsequences R 1 and R 2 of approximately equal lengths. If |R| is an even In step 1, RNA-Seq reads are aligned to the reference genome using a genome aligner. In this step, exonic reads are expected to be mapped to the reference genome. The unaligned reads are typically either junction reads or reads with insertions or deletions. In step 2, the remaining unaligned reads are split into three parts and then the split parts are aligned to the reference genome using a genome aligner. If any part of a split read is aligned, then the part is used as an anchor to localize the whole read to the genomic locus.
Similarly, we define the function Split(R, 3) that splits the read R into 3 subsequences of approximately equal lengths.

The Reads Localization Algorithms
Algorithm Reads Localization and Binning (single-end) Input: A set of single-end reads R, a predefined intergenic distance d Output: Binning results L /* L is a set of bins, each bin consists of reads falling into a genomic locus */ /* Align exonic reads */ P ← ∅ /* Initialize the localization result set */ Align R to the reference genome using a genome aligner /* Default is Bowtie 2 */ for each read R i ∈ R do if R i is aligned to the reference genome at position P i then /* Split and localize the rest */ S ← ∅ /* Initialize the split set */ for each leftover read Align S to the reference genome using a genome aligner /* Default is Bowtie 2 */ for each triple (R i1 , R i2 , R i3 ) ∈ S do if R i1 , R i2 , or R i3 is aligned to the reference genome at position P i then P ← P ∪ {(R i , P i )} end if end for /* Bin the reads to genomic loci */ Sort P by increasing order of genomic position. L ← ∅ /* Initialize the binning sets */ B ← {(R 1 , P 1 )} /* Initialize the current bin as the read with smallest genomic position */ for each read-position pair Align S to the reference genome using a genome aligner /* Default is Bowtie 2 */ for each quadruple (R i11 , R i12 , R i21 , R i22 ) ∈ S do if R i11 , R i12 , R i21 , or R i22 is aligned to the reference genome at position P i then P ← P ∪ {(R i , P i )} end if end for /* Bin the reads to genomic loci */ Sort P by increasing order of genomic position. L ← ∅ /* Initialize the binning sets */ B ← {(R 1 , P 1 )} /* Initialize the current bin as the read with smallest genomic position */ for each read-position pair (R i , P i ) ∈ S and i > 1 do if

Stage II: Seeding and Growing
In the seeding and growing stage (stage II), our goal is to construct long backbone sequences from a set of short reads. At the beginning, all the reads are marked as "active". The algorithm then randomly picks a seed from the active set of reads and greedily extends it to both the left and right directions. If there are multiple reads that can extend a seed sequence, the read that most extends the seed sequence is chosen. After the sequence is fully grown, all the reads covered by the seed sequence are marked as inactive. For those remaining active reads, the algorithm repeatedly grows other backbone sequences until no active read left. To ensure that these fully grown backbone sequences are connected to each others, we extend each of them with one inactive read to both the left and right directions if possible. The extended ends will serve as connecting points to construct the isoform graph in the cutting and patching stage (stage III).
An example is displayed in Figure 2a and b. In the seeding and growing stage, the seeds are randomly chosen. In the left panel of Figure 2a, assuming that a green exonic read is chosen to be the first seed, the algorithm first extends the seed to the left until it reaches the leftmost purple read, then it extends the seed to the right until it reaches the rightmost blue read. After the first backbone sequence is constructed, the algorithm marks the following reads as inactive: purple exonic, purple-green junction, green exonic, green-blue junction, and blue exonic reads. The remaining active reads are purple-yellow junction, green-yellow junction, yellow exonic, and yellow-blue junction reads. Assuming a yellow exonic read is chosen to be the second seed, the algorithm grows the seed until it reaches the green-yellow junction read on the left and the yellow-blue junction read on the right. All active reads that are covered by this new sequence are marked as inactive. As described above, the second backbone sequence is then extended by one inactive green read on the left and by one inactive blue read on the right. After this step, only the purple-yellow junction read remains active and is chosen to be the third seed. Since there is no other active read, the seed is extended with one inactive purple read on the left and one inactive yellow read on the right. Similarly, in the right panel of Figure 2a, the first seed is the highlighted yellow exonic read, the second seed is the green-blue junction read, and the third seed is the purple-yellow read. By choosing two different sets of seeds, the seeding and patching stage outputs two different sets of backbone sequences. However, after the cutting and patching stage, both sets of backbone sequences give rise to the same isoform graph as displayed in Figure 2c, from which the same set of transcript sequences are enumerated. Below is the mathematical description of the seeding and growing algorithm.

Definitions
where l o is the length of overlapped letters and k is the minimum overlap cutoff.
In this case we denote the overlapping length between the two short reads as overlap(P, Q) = l o and overlap(Q, P ) = −l o . If short reads P and Q do not overlap with each other, we denote overlap(P, Q) = overlap(Q, P ) = 0. Without loss of generality, we define the overlap between any two sequences in the same manner.
Definition 4 (Extension). Given a short read P and a set of short reads We denote the right extension as T r = ext right(P, T ) and the left extension as T l = ext lef t(P, T ). We also denote ext lef t(P, T ) = () or ext right(P, T ) = () if P cannot be extended.
Definition 5 (Merger). Two short reads can be merged to make a longer sequence if they overlap with each other. Given short reads P = (p 1 p 2 ...p l ) and Q = (q 1 q 2 ...q l ) with overlap(P, Q) = l o , the merged sequence is M = (p 1 p 2 ...p l−lo q 1 q 2 ...q l ) and is denoted as M = merge(P, Q). Without loss of generality, we define the merger between any two sequences in the same manner.
Definition 6 (Cover and Subsequence). Given two sequences S = (s 1 s 2 ...s l ) and Denote it as cover(S ′ , S) = T rue. We also say that S ′ is a subsequence of S and denote it as S ′ = sub(S, i, l ′ ). In the left panel, the first backbone sequence is grown from a green exonic seed and is extended to the leftmost purple read and to the rightmost blue read. The first backbone sequence spans across the purple, green, and blue exons. Likewise, the second sequence is grown from a yellow exonic read and the third sequence is grown from a purple-yellow junction read. The second backbone sequence covers the green-yellow junction, the yellow exon, and the yellow-blue junction whereas the third exon covers the purple-yellow junction. In the right panel, the first sequence is grown from a yellow exonic seed and is extended to the leftmost purple read and to the rightmost blue read. The first backbone sequence spans across all 4 exons. Likewise, the second sequence is grown from a green-blue junction read and the third sequence is grown from a purple-yellow junction read. The panel (b) displays the output backbone sequences of the seeding and growing stage, where the shadings show the overlaps between the backbone sequences. The seeding reads are displayed on the left of their grown sequences. The panel (c) displays the isoform graph after the cutting and patching stage, where the overlapping sequences serve as connecting points to glue the backbone sequences together. After this stage, the two sets of backbone sequences lead to the same isoform graph. The panel (d) displays the output transcript sequences.

The Seeding and Growing Algorithm
Algorithm Seeding and Growing Input: A set of short reads T = {T 1 , T 2 , ..., T t } and the minimum overlap length k Output: A set of sequences S S ← ∅ /* Initialize S */ while T = ∅ do P ← T 1 /* Set P as the first short read in T */ S P ← P /* Initialize S P as the sequence of P */ repeat T l ← ext lef t(P, T ) for P ← T r /* Set P as its right extension */ until ext right(P, T ) = () S ← S ∪ {S P } /* Add sequence S P to set S */ end while

Stage III: Patching and Cutting
Given a set of sequences outputted by the seeding and growing algorithm (stage II), the patching and cutting algorithm (stage III) dynamically constructs an isoform graph. Initially, vertices in the graph are set as sequences from the input and the connection relationship (connected or not connected) between them is "unknown". For every pair of vertices, the algorithm checks their connection relationship flag. If it is "unknown", the two vertices are tested for patching. If patched, then their sequences are cut into segments, the vertices are split, directed edges are added, and the relationship flags between all involved vertices are updated accordingly. Sometimes this patching and cutting process may cause unnecessary cuts, i.e., it may generate some pairs of vertices connected only by one edge in the graph. If this happens, the two vertices are merged into one. After the isoform graph has been constructed, each linear path starting from a vertex with no incoming edges and ending at a vertex with no outgoing edge is put into an isoform structure set, and the joint sequence based on that path is put into a contigs set ( Figure 3). Below is the mathematical description of the patching and cutting algorithm.

Definitions
Definition 7 (Patch). Given two sequences S = (s 1 s 2 ...s l ) and S ′ = (s ′ 1 s ′ 2 ...s ′ l ′ ), S ′ is called a left patch of S if we can find an i and a p ≥ k for which sub(S ′ , 1, p) = sub(S, i, p). We denote it as patch lef t(S ′ , S) = (i, p). Similarly, S ′ is called a right patch of S if we can find an i and a p ≥ k for which sub(S ′ , l ′ − p + 1, p) = sub(S, i, p). We denote it as patch right(S ′ , S) = (i, p).
Definition 10 (Isoform graph, isoform structure and contig). An isoform graph G = (V, E) is a directed graph where each vertex is a sequence of letters from {A, T, C, G} and each edge from vertex V i to vertex V j indicates that there exists a short read R which is a subsequence of the joint sequence of V i and V j , i.e., cover(R, joint(V i , V j )) = T rue. Given an isoform graph G = (V, E), an isoform structure is a linear path in G starting at a vertex with no incoming edge and ending at a vertex with no outgoing edge. For an isoform structure with n vertices as V 1 → V 2 → ... → V n , the corresponding contig is joint(V 1 , V 2 , ..., V n ).

A Patching and Cutting Algorithm
Algorithm Patching and Cutting Input: A set of sequences S and the minimum overlap length k Output: An isoform graph G = (V, E), a set of isoform structures I, and a set of contigs C V ← S /* Initialize V , a set of vertices */ E ← ∅ /* Initialize E, a set of edges */ for each unknown" /* Initialize F , the relationship flags between vertices */ end for tices V i and V j , cutting location variables v and p, and direction variable dir Replace V j with V m /* Update outgoing edges of V j */ end for end if end for Procedure Get Isoform Structures(G) Input: An isoform graph G Output: An isoform structures set I I ← ∅ /* Initialize I */ for each vertex V i ∈ V that has no incoming edge do I ← {V i } end for for each path

Error control
The above algorithms do not take into account the sequencing errors in short reads. In reality, due to those errors, the short reads do not always perfectly overlap with each other. Therefore the algorithms must have error tolerance capability in order to process real sequencing data. For this purpose, we slightly modify the following definitions: Definition 11 (Overlap, updated). A short read P = (p 1 p 2 ...p l ) overlaps with Q = (q 1 q 2 ...q l ) if (p l−lo+1 p l−lo+2 ...p l ) = (q 1 q 2 ...q lo ) for l o ≥ k with at most e 1 l o differences (errors) and no more than e 2 contiguous errors, where l o is the length of overlapped letters, k is the minimum overlap cutoff, e 1 is the maximum error rate in any short read, i.e., there are at most e 1 * l sequencing errors in a short read of length l. Here e 2 is a small constant, for instance, 2. Without loss of generality, we redefine the overlap between any two sequences in the same manner.
Definition 12 (Cover and Subsequence, updated). Given two sequences S = (s 1 s 2 ...s l ) and ..s ′ l ′ ) for some 1 ≤ i ≤ l with at most e 1 * l ′ differences (errors) and no more than e 2 contiguous errors. We also say that S ′ is a subsequence of S.
With these modified definitions, the proposed algorithms can be applied to process data with sequencing errors.

The Localization and Binning Algorithm
Given a read aligner with running time T (g, l) per short read, where g is the length of the reference genome and l is the short read length, the localization stage runs in O(T (g, l)r) time, where r is the total number of short reads. Running time of the binning step is O(r log r + r) = O(r log r) because sorting genomic position takes O(r log r) time and the binning time after sorting is O(r). Therefore, the total time for localization and binning is O(T (g, l)r + r log r). Since T (g, l) is Ω(g) and g ≫ log r, we have O(T (g, l)r + r log r) = O(T (g, l)r), i.e., the time complexity of the localization and binning algorithm is mainly determined by the time complexity of the genome aligner.  The total time complexity of the algorithms for short reads localization, transcriptome assembly and isoforms reconstruction is thus O(T (g, l)r + sr). Since T (g, l) is Ω(g) and g ≫ s, the complexity is mainly determined by the read mapping algorithm.

Simulation Experiments
We used FluxSimulator [http://flux.sammeth.net], an open source software package simulating whole transcriptome sequencing experiments, to generate RNA-Seq data sets. FluxSimulator first randomly generates integer copies of each splicing isoform according to the user-provided annotation file. It then constructs an amplified, size-selected library and sequences the library in silico. The resulting cDNA fragments are then sampled randomly for simulated sequencing, where the initial and terminal ends of each selected fragment are reported as reads. The current version of FluxSimulator (1.2) provides two sequencing error models, one is for read length 76 and another is for read length 36. The default error model was used to generate six data sets consisting of 15 million, 20 million, 30 million, 50 million, 75 million, and 100 million paired-end reads with length 76 from the reference protein-coding transcripts available from the Ensembl database (version GRCh37.69). More details of the simulated data sets and the error model can be found in the Appendix. Both simulated RNA-Seq data sets and implementation of SPATA are available for download at http://sammate.sourceforge.net/.
We compared the performance of SPATA with Cufflinks [19] and Trinity [41]. Cufflinks is a reference-guided transcriptome assembly whereas Trinity is a de novo transcriptome assembly [15,16]. As mentioned before, SPATA is a hybridization of the reference-guided and the de novo strategies, i.e., it uses the reference genome to localize reads to genomic loci and then assembles the read sequences falling into each genomic locus using a new assembly algorithm. Simulated FASTQ files were used as input for SPATA and Trinity. Since Cufflinks does not process raw reads, we used TopHat [23] to align reads from the FASTQ files to the human reference genome. We then used TopHat's BAM files as the input for Cufflinks. We used default parameters to run all the tools.

Precision and Accuracy Analyses
For each dataset, the ground truth is a set of expressed transcript sequences whereas the output of each assembly software is a set of contigs, i.e., a continuous sequence of bases constructed from the short reads in the FASTQ files. In the rest of the manuscript, we use the term reference transcripts to refer to the expressed transcript sequences. Output contigs of each assembly program were compared against the reference transcripts using the pairwise sequence alignment  Figure 5: Accuracy analysis of Cufflinks, Trinity, and SPATA using simulated data. The figure displays the fraction of reference transcripts fully reconstructed by different coverage quantiles. The six panels are arranged in an increasing order of read numbers from 15 million to 100 million of paired-end reads. From the ground truth, we know the reference transcripts and their expression abundances. In each panel, the horizontal axis shows the expression percentile (in 10% increments) and the corresponding mean coverage of the reference transcripts. The left vertical axis shows the fraction of reference transcripts that can be fully reconstructed by each method. The right vertical axis displays the number of reference transcripts that were fully reconstructed by each method. A transcript is considered to be fully reconstructed by a method if it has at least 90% base-identity with a contig from the method's output. We examine all reference transcripts having mean coverage of 5 or higher.

Cufflinks
Trinity SPATA Figure 6: The number of reference transcripts reconstructed by Cufflinks, Trinity, and SPATA using simulated data. The horizontal axis shows the number of reads for each data set in an increasing order whereas the vertical axis shows the total number of reference transcripts that were fully reconstructed. A transcript is consider to be fully reconstructed by a method if it has at least 90% base-identity with a contig from the method's output. We examine all reference transcripts having mean coverage of 5 or higher.
program SSAHA2 [45]. We evaluated and compared the performance of the assembly programs according to "precision" and "accuracy". When output contigs were aligned to the reference transcripts, we define the precision as the fraction of contigs that were properly mapped to the reference transcripts. A contig is considered to be mapped if at least 90% of its sequence was covered by a reference transcript ( M M+N +G ≥ 90%, where M is the number of matched bases, N is the number of mismatched bases and G is the cumulative insertion and deletion length). When the reference transcripts were aligned to the constructed contigs, we define the accuracy as the fraction of reference transcripts that can be fully reconstructed. A transcript is defined to be fully reconstructed if it has at least 90% base-identity with one of the output contigs.
When evaluating the performance of different assembly software, high values for both precision and accuracy are important. On one hand, a high precision indicates that most of the method's contigs are subsequences of the reference transcripts. On the other hand, a high accuracy indicates that a method reconstructed a large number of the reference transcripts.
Regarding precision, we calculated the fraction of output contigs of Cufflinks, Trinity, and SPATA that were properly mapped to the reference transcripts. Figure 4 displays the fraction of contigs that were mapped to reference transcripts. The performance comparison in terms of precision is overall comparable with slightly higher values for Cufflinks. It might due to the following reasons: First, Cufflinks only processes reads that were aligned to the reference genome and discards reads with sequencing errors. Second, Cufflinks attempts to identify the minimum number of isoforms that cover the majority of reads.
Regarding accuracy, we calculated the fraction of reference transcripts that can be fully reconstructed by Cufflinks, Trinity, and SPATA. Figure 5 displays the number of reference transcripts that were fully reconstructed across a broad ranges of expression levels and sequencing depths. The expression of the reference transcripts is calculated from the simulated ground truth origin of the paired-end reads, i.e., the mean coverage of a transcript is calculated as the number of reads multiplied with the read length divided by the transcript length. We examined all reference transcripts having mean coverage of 5 or higher.
In general, Cufflinks performed better at low expression levels whereas SPATA and Trinity performed better at high expression levels ( Figure 5). Cufflinks may perform better at lower expression levels since the small gaps within a transcript caused by sequencing errors or lack of read coverage can be corrected using the reference genome [19,15,16]. Trinity and SPATA, however, outperform Cufflinks in assembling highly reference transcripts. A possible reason is that the parsimony strategy used by Cufflinks attempts to identify the minimum number of isoforms to cover the majority of reads while filtering out uncertain reads. The latter may result in missing some alternative splicing events [19,15,16,22].
Another notable phenomenon is that SPATA reconstructs more reference transcripts than Cufflinks and Trinity for any of the six data sets. Figure 6 displays the total number of reference transcripts reconstructed by each method across all of the six data sets. Interestingly, Trinity and SPATA reconstruct more reference transcripts as the number of reads increases whereas the number of reference transcripts reconstructed by Cufflinks fluctuates.

Robust Analysis Using Incomplete Reference Genomes
In the above analyses, Cufflinks and SPATA used the human reference genome version GRCh37.69 downloaded from Ensembl. The same genome sequence was used to generate the simulated data sets by FluxSimulator. We consider this latest version of human genome as the "complete" reference genome. In order to access the robustness of Cufflinks and SPATA against incomplete reference genomes, we also examined their performance using different versions of the human reference genome.  Figure 7 displays the precision and accuracy of Cufflinks, Trinity, and SPATA for the data set consisting of 15 million reads. Since Trinity did not use any reference genome, its constant performance was used as a benchmark to compare the performance of Cufflinks and SPATA. In Figure 7, the left panel displays the precision of the assembly tools whereas the right panel displays their accuracy. The number of reference transcripts reconstructed by Cufflinks, Trinity, and SPATA is displayed in Figure 8 in the Appendix. In terms of precision, the performance of Cufflinks and SPATA are comparable and are consistently higher than that of Trinity. This translates into a robust precision of Cufflinks and SPATA, regardless of quality of the reference genome being used. In terms of accuracy, as the quality of the reference genome declines, the performance of Cufflinks and SPATA also declines and a sharper drop is observed for Cufflinks. For the reference genome sequence hg10, a low-quality reference genome, the accuracy of Trinity and SPATA is almost the same and is higher than that of Cufflinks. For all other reference genome sequences, SPATA achieves a better accuracy than Cufflinks and Trinity.

Conclusion
In this paper, we presented SPATA, a hybrid transcriptome assembly approach, to reconstruct transcriptomes via a novel three-stage algorithm: localization, seeding and growing, and patching and cutting. We first assessed the performance of SPATA by reconstructing the transcritomes using multiple simulated RNA-Seq data sets for which SPATA achieves high precision and accuracy. We also showed that SPATA reconstructs more reference transcripts than the selected transcriptome assembly tools for all of the simulated data sets. We further demonstrated the robustness of SPATA against incomplete reference genomes by using different versions of human reference genome over the past 12 years. SPATA consistently achieved a high precision and re-constructed more reference transcripts than the competing methods. Overall, the analyses favor the use of SPATA in transcriptome reconstruction wherever a reference genome is available, regardless of quality.
SPATA is also expected to be a valuable tool to assemble transcriptomes of non-model organisms, where the reference genome sequences are likely to be incomplete. Despite advances in sequencing technologies, assembling a complete reference genome is still costly and difficult for many species. As a result, most of the existing genomes are available only as unfinished drafts with gaps and excessive assembly errors [46,47,48]. In these cases, SPATA is expected to reliably capture both known transcript structures and novel variations due to its high precision and accuracy. SPATA is also conveniently accessible for both informatics and life science researchers via an easy-to-use GUI software. Figure 8 displays the number of reference transcripts that were fully reconstructed by Cufflinks, Trinity, and SPATA using different reference genomes for the data set consisting of 15 million reads. Since the total number of reference transcripts is a constant, the number of the reference transcripts reconstructed by the assembly tools also reflects their accuracy. For the low-quality reference genome genome hg10, SPATA and Trinity reconstruct almost the same number of reference transcripts and reconstruct more than Cufflinks. For all other reference genome sequences, SPATA reconstructs more reference transcripts than Cufflinks and Trinity.

Appendix
The statistics of the simulated data are displayed in Ta- Reference Genomes Figure 8: The number of reference transcripts reconstructed by Cufflinks, Trinity, and SPATA using different reference genomes for the data set consisting of 15 million reads. The horizonal axis displays the reference genomes being used in chronological order. The vertical axis displays the total number of reference transcripts that were fully reconstructed. A transcript is consider to be fully reconstructed by a method if it has at least 90% base-identity with a contig from the method's output. We examine all reference transcripts having mean coverage of 5 or higher. ble 1, Figure 9 and Figure 10. Table 1 shows the number of reference transcripts in each data set. Figure 9 displays the quality score distribution of sequenced bases across 76 read positions whereas Figure 9 displays the expression abundance distribution of the reference transcripts.