A scalable and memory-efficient algorithm for de novo transcriptome assembly of non-model organisms

Background With increased availability of de novo assembly algorithms, it is feasible to study entire transcriptomes of non-model organisms. While algorithms are available that are specifically designed for performing transcriptome assembly from high-throughput sequencing data, they are very memory-intensive, limiting their applications to small data sets with few libraries. Results We develop a transcriptome assembly algorithm that recovers alternatively spliced isoforms and expression levels while utilizing as many RNA-Seq libraries as possible that contain hundreds of gigabases of data. New techniques are developed so that computations can be performed on a computing cluster with moderate amount of physical memory. Conclusions Our strategy minimizes memory consumption while simultaneously obtaining comparable or improved accuracy over existing algorithms. It provides support for incremental updates of assemblies when new libraries become available.


Background
As the advance of high-throughput sequencing makes it possible to sequence billions of bases in a single experiment, this shift in the availability of genomic data allows researchers to focus on biological questions in nonmodel organisms. With the increased availability of de novo assembly algorithms that are designed specifically for assembling millions of short reads [1][2][3][4][5][6][7][8], it becomes possible to study entire genomes or transcriptomes by investigating the assembled sequences.
To obtain a transcriptome, many RNA-Seq libraries are constructed under different experimental conditions or developmental stages, with each library corresponding to one sequencing run of a biological sample. It is preferable to utilize as many libraries as possible to construct one single assembly for each species, as the large amount of data enables simultaneous expression analysis and an increase in coverage support of transcripts that may not be highly expressed under some conditions. In order to obtain the best transcriptomic profile for a given species, there is a need to integrate large amount of accumulated data together from disparate projects and create updated transcriptome assemblies as new data become available. This creates a significant challenge for de novo assembly algorithms, since computational resources are often limited in individual labs while the computational time and memory requirement increase rapidly as the number of reads increases.
These computational challenges have motivated the development of algorithms that are specifically designed for performing transcriptome assembly. While algorithms such as Oases [9] and Trinity [8] aim to extract as much information as possible, they are very memoryintensive, limiting their applications to small data sets with few libraries that biologists produce during experiments. Algorithms such as SOAPdenovo-Trans [10] and Trans-ABySS [11] have high memory requirements for large data sets.
To address these difficulties, our goal is to develop transcriptome assembly algorithms that recover alternatively spliced isoforms while utilizing as many RNA-Seq libraries as possible that contain hundreds of gigabases of data. We subdivide the computations into two stages, in which the first stage collects information from each library independently and in parallel, and the second stage merges these results together while minimizing needed computations. To reduce memory consumption so that computations can be performed on a computing cluster with moderate amount of physical memory, we develop new techniques to enumerate k-mer frequencies in the first stage. We impose appropriate cutoffs in the second stage in order to obtain comparable or improved accuracy over existing algorithms. This strategy supports incremental updates of assemblies when new libraries become available since only the second stage needs to be rerun.
We evaluate the performance of our algorithm by constructing transcriptome assemblies using publicly available libraries from model organisms, and comparing our assemblies to the ones obtained from SOAPdenovo-Trans, Trans-ABySS, Oases and Trinity. We evaluate our performance on non-model organisms both by obtaining publicly available libraries from the naked mole rat Heterocephalus glaber and by constructing new RNA-Seq libraries for the blow fly Chrysomya rufifacies.

De Bruijn graph
Given a set of reads and a parameter k, a de Bruijn graph is defined by taking each k-mer that appears within the reads as a vertex. Two k-mers s 1 s 2 · · · s k and s 2 · · · s k s k+1 are connected by a directed edge if the (k + 1)-mer s 1 s 2 · · · s k s k+1 appears in the reads and the (k −1)-suffix of the first k-mer is the same as the (k − 1)-prefix of the second k-mer, where s 1 and s k+1 can be arbitrary letters. By linking together the same k-mer that appears in different reads, the de Bruijn graph can be used to implicitly assemble these reads [12,13]. Since the size of the de Bruijn graph depends on the number of distinct k-mers from the reads that is often much smaller than the total size of reads, this strategy is especially suitable for assembling high-throughput sequencing data [2,3,[5][6][7].

Independent computation of k-mer frequencies
While most existing short read assembly algorithms use hashing techniques [5,8,14] or suffix arrays [15] to enumerate k-mer frequencies, the memory requirement per k-mer can be high with large multiplicative constants. While techniques such as sparse hashing (http:// code.google.com/p/google-sparsehash) or entropy-based compression [16] can be used to reduce the memory overhead per k-mer, there is a need to handle collisions and the memory requirement can still be high. Recently, two algorithms DSK [17] and KMC [18] were developed based on disk-based partitioning of the kmer space, which allow the user to specify a memory consumption limit.
We consider the following iterative algorithm to enumerate k-mer frequencies independently for each library (see Fig. 1). At the start of each iteration, we assume that a list of all k -mers that appear in the library in either the forward or the reverse complementary direction are given in sorted order for k < k. This list can be represented by encoding each nucleotide by two bits and using a 64bit or 128-bit integer to encode each k -mer. Suppose that there are n such k -mers. We create an array of size 4n that contains four slots for each k -mer. We use it to count the number of each (k + 1)-mer that appears in the library by using binary search to find the location of its k -prefix within the array and updating one of the four slots that corresponds to its last nucleotide. At the end of the iteration, we remove slots with zero counts to obtain a list of all (k + 1)-mers.
To make sure that each edge in the de Bruijn graph corresponds to a (k + 1)-mer that appears in the library, we repeat this procedure until the frequencies of (k + 1)-mers are obtained and store all k-mer frequencies along with edge information. To reduce computational time, we start the process with a moderate value of k (between 10 to 15) by assuming that all k -mers appear Fig. 1 Illustration of the iterative algorithm to enumerate k-mer frequencies. For the k -mer a 1 · · · a k , its two frequency slots with zero counts for nucleotides c and t are removed to obtain (k + 1)-mers a 1 · · · a k a and a 1 · · · a k g in the library. One advantage of this procedure is that the memory requirement per k-mer is low with a multiplicative constant of four. When multiple assemblies with different values of k are needed, frequencies with smaller values of k for one assembly can be used to obtain frequencies with larger values of k for another assembly. This significantly reduces the computational time over all values of k.

Construction of de Bruijn graph
Given a list of all k-mer frequencies in sorted order for each library, we combine these lists by performing a merge sort and adding the corresponding frequencies for each k-mer. We exclude a k-mer from the de Bruijn graph if its frequency is less than a given k-mer coverage cutoff c. This strategy is different from the one employed by other short read assembly algorithms such as Velvet [5] or ABySS [6], which apply the cutoff after constructing the de Bruijn graph and removing redundant paths that are of less support. Since our results show that there is a need to increase the values of k and c to obtain good performance as the size of the data set increases, our strategy keeps the size of the de Bruijn graph manageable. Edges in the de Bruijn graph are constructed by linking together adjacent k-mers through binary search. After the initial construction of the de Bruijn graph, each linear path that contains a maximal succession of vertices with no branches is collapsed into a single node. While the entire merging process needs to be run sequentially for each given setting of k and c, different settings of k and c can be run in parallel.

Construction of splicing graphs
To simplify the de Bruijn graph, we remove short tips that may correspond to sequencing errors by iteratively removing end nodes with sequence length less than 2k. We follow the strategy in [19] to obtain splicing graphs from the de Bruijn graph so that each splicing graph mostly represents alternatively spliced variants of only one gene, with new strategies to handle paired-end reads. Note that this strategy is different from algorithms that construct a set of predicted transcripts from the de Bruijn graph [8][9][10][11].
In order to remove obvious SNPs that create branches in the graph, we search for split-then-merge branching structures in which all the branches from a node merge immediately into a single node and the sequences associated with each branch are of the same length with very few mismatches. Successive split-then-merge structures are merged into a single node.
For each paired-end read, we identify the node u in the collapsed de Bruijn graph in which the last k-mer of the forward read resides and the node v in which the first k-mer of the reverse read resides, and increase the frequency count of the paired edge u → v by one. We repeat this procedure over all paired-end reads and retain all paired edges that have frequency counts above a given cutoff c 2 that is proportional to the total number of bases in the data set. The resulting de Bruijn graph contains two types of edges, including normal edges and paired edges.
Our results show that there is always a big tangle in the de Bruijn graph that contains a large number of nodes within a single connected component, while most of the other tangles are much smaller. In order to address these tangles, we decompose each connected component into strongly connected components, in which each strongly connected component is either just a single edge or a maximal subgraph with each node reachable from all other nodes. The regions within a strongly connected component that are not just a single edge represent the complicated regions that always contain a cycle, while the other regions represent the simpler regions in which each connected region is likely to belong to the same gene. With the assumption that the reads are not strandspecific, it is also possible to have forward-backward tangles in which a connected component contains both a forward node and its corresponding backward node Table 1 Data sets used in the evaluation of transcriptome assembly, with organism denoting the organism, type denoting whether the organism is model or non-model, libraries denoting the total number of libraries, size denoting the total number of bases in all the reads after quality trimming, reference denoting the publication that describes the libraries, and tick marks within assembly on 32 GB machines denoting the algorithms that can be used for assembly on machines with 32 GB physical memory We extract the strongly connected components that are not just a single edge and all the forward-backward tangles. We ignore the junction information within these subgraphs, and treat each node as an individual splicing graph that consists only of a single node. We remove these nodes along with their adjacent edges, and extract each connected component in the remaining graph as a splicing graph that does not contain cycles. We remove the overlapping sequence fragments within the nodes that arise according to the definition of a de Bruijn graph, and make the junction locations precise in obvious cases when a node does not have multiple incoming edges and multiple outgoing edges at the same time. We retain one of the two possible orientations for each splicing graph. Only splicing graphs with length (in nucleotides) of the longest path of at least 100 are retained.
In order to study expression of nodes in a splicing graph, we incorporate the de novo expression measure of number of reads per kilobase of node per million reads (RPKM) developed in [19]. This measure is similar to the number of reads per kilobase of exon per million mapped The predicted units are locus for SOAPdenovo-Trans that is represented as a splicing graph containing nodes and edges, transcript (trans) for Trans-ABySS that is a linear concatenation of constituent nodes, and splicing graph for ASplice. For SOAPdenovo-Trans and ASplice, N50 denotes the N50 value of the length (in nucleotides) of the longest path in each splicing graph. For Trans-ABySS, N50 denotes the N50 value of the length of a predicted transcript, and only predicted transcripts of length at least 100 are retained. Total hits denotes the total number of hits from nucleotide BLAST search of nodes to the transcriptome of the same organism. Isoforms are considered to be the same gene. Only the top hit with E-value below 10 −7 is considered. Hits from nodes within the same predicted unit to the same gene are counted only once. Unique hits denotes the number of unique hits to different genes. Mem (GB) denotes the physical memory requirement in gigabytes over all stages of each algorithm reads used by [20] and [21], except that reads that appear in the assembly are used instead of mapped reads, and each node in a splicing graph is evaluated instead of each exon, with each read that contains a k-mer within a node contributing to that node. Within each node, a RPKM estimate is computed independently for each library. Alternatively, measures similar to transcripts per million (TPM) [22] can be used, which are more comparable across libraries.
In order to make the results directly applicable to downstream analysis, we represent each assembly in an annotated FASTA format, in which each splicing graph is given as a collection of nodes, with connecting normal and paired edges and RPKM values for each library embedded within the name of each node. Since it is possible to have empty nodes that do not contain any nucleotide after adjustment of junction locations, RPKM values are computed before junction adjustment to reflect the original coverage values across a branch. Such empty nodes correspond to additional isoforms that skip nodes (e.g., exons) within a branch.

Model organisms
To compare the performance of our algorithm ASplice to other algorithms, we extracted reads from publicly available RNA-Seq libraries in model organisms Schizosaccharomyces pombe, Arabidopsis thaliana and Drosophila melanogaster (see Table 1). We trimmed each read by removing all positions including and after the first position that has a quality score of less than 15. We applied our algorithm to obtain a de Bruijn graph for a given k-mer length and a given k-mer coverage cutoff c. We compare the performance of our  algorithm ASplice to SOAPdenovo-Trans and Trans-ABySS on machines with 32 GB physical memory (except for D. melanogaster, in which there is not sufficient memory to run SOAPdenovo-Trans). Nucleotide BLAST search is applied to the transcriptome of the same organism to evaluate the performance. Since each algorithm returns different structures, note that the results are not completely comparable. Table 2 shows that while there were performance tradeoffs among different values of k and the k-mer coverage cutoff c, SOAPdenovo-Trans and Trans-ABySS generally recovered longer transcript structures and more genes in the transcriptome (except for S. pombe, in which ASplice recovered more genes). Figure 2 shows that ASplice had higher specificity with respect to the percentage of predicted positions that are included in the top BLAST alignments. ASplice also often had comparable sensitivity with respect to the percentage of nucleotide positions in the transcriptome that are recovered through the top BLAST alignments, especially for larger values of the k-mer coverage cutoff c when the assembly conditions are more stringent, and had higher sensitivity for S. pombe. Figure 3 further shows that, with respect to alternative splicing junctions that are derived from the splicing graphs and annotated positions of the gene transcripts, ASplice was often more conservative, and could have higher sensitivity when the assembly conditions are more stringent (no comparisons were made to SOAPdenovo-Trans since its assemblies often contain many gap positions around junctions, making comparisons difficult). The poor performance with respect to alternative splicing junctions in A. thaliana is due to the relatively small size of the data set. Figure 4 shows that our de novo expression estimates were highly correlated to the ones obtained from applying RSEM [22] to map the reads in each library to the transcriptome of the same organism, with 0.49 ≤ R 2 ≤ 0.84 for gene transcripts without alternative splicing and 0.28 ≤ R 2 ≤ 0.40 for gene transcripts with alternative splicing, where R 2 denotes the coefficient of determination that measures how well the data fit a regression line. The lower R 2 in the latter case is probably due to higher assembly difficulties.
In terms of memory requirement, Table 2 shows that while SOAPdenovo-Trans and Trans-ABySS had large increases as the size of the data set increases (compare to   to reduce memory requirement during the parallel stage of ASplice, the maximum memory consumption is obtained during the sequential stage, in which the main structure to store is the de Bruijn graph after the k-mer coverage cutoff c is applied.

Non-model Heterocephalus glaber
We also assess the performance of our algorithm in the non-model naked mole rat Heterocephalus glaber (there is not sufficient memory to run SOAPdenovo-Trans on machines with 32 GB physical memory for this data set). Table 3 shows that large values of the k-mer coverage cutoff c were needed to obtain reasonable assemblies due to the large size of the data set (see also Table 1). Similar trends were observed as before when translated BLAST search to the rat Rattus norvegicus is applied, with ASplice recovering more genes when the assembly conditions are more stringent.

Non-model Chrysomya rufifacies
We applied our algorithm to assemble the transcriptome of the blow fly Chrysomya rufifacies from a set of RNA-Seq libraries that we have constructed, which includes the following developmental stages: embryos, first instar larvae, second instar larvae, predator and non-predator third instar larvae, early pupae, mid pupae, late pupae, thelygenic and arrhenogenic adult females, and adult males. There are totally 66 libraries with 6.8 G reads and average read length 86 after quality trimming.
The blow fly C. rufifacies has monogenic sex determination in which a female either produces only female offspring (thelygenic) or produces only male offspring (arrhenogenic) based on the genotype [23], which is a distinct mechanism among flies. Sex determination in flies is typically achieved (in part) by alternative splicing, in which sex-specific isoforms of genes like doublesex and transformer lead to female or male development [24,25]. Within the genus, there is also a human-associated male-eye phenotype that is hypothesized to have evolved multiple times in concert with human civilization [26]. Genomic tools for this blow fly enable the study of the evolution of sex determination and co-evolution with humans.
Since there is not sufficient memory to run either SOAPdenovo-Trans or Trans-ABySS on machines with 32 GB physical memory for this large data set with 590 G bases, we only run ASplice, which allows computations to be performed on a computing cluster due to its low memory requirements. We considered larger values of k and further subdivided large libraries into multiple parts during the parallel stage. Table 4 shows that the assemblies were of high quality, with long splicing graphs, moderate amount of branches that may represent alternative splicing, and over 60% of the D. melanogaster genes recovered. The ratio of the total number of BLAST hits from different splicing graphs to the number of unique BLAST hits to different genes was between 1.5 to 2, indicating a small amount of sequence fragmentation of the same gene into different splicing graphs. There were only a small number of splicing graphs that have BLAST hits to more than one gene, and the maximum number of different genes that have BLAST hits to a splicing graph was small, thus each splicing graph specifies the alternatively spliced variants of one gene in most cases.
By comparing to the D. melanogaster homologs of assembled nodes (see Fig. 5), we found expected alternative splicing in the doublesex gene and consistent bias of expression within female-specific (upper right node of length 117) and male-specific (lower right node of length 876) segments and between thelygenic and arrhenogenic females (with generally higher expression within thelygenic females in the upper right node).

Small Drosophila melanogaster libraries
Since Oases and Trinity are very memory-intensive, we assess their performance by extracting reads from a small set of three D. melanogaster RNA-Seq libraries in [27] at the sequence read archive [28] that includes the following developmental stages: 2-16 hours embryos (SRR058885), third instar larvae (SRR059066), and mixed pupae (SRR042298). These libraries have 1.8 G bases after quality trimming, and k is fixed to 25. Table 5 shows that Oases had the longest assemblies and recovered the largest number of genes. While ASplice had that is related to the doublesex gene in D. melanogaster. Each node is represented by a rectangle of width that is proportional to (and labeled by) the length of its sequence. The histogram adjacent to each node gives the RPKM values for each adult library, with each light bar and the dark bar to its right denoting an adult female and its mate (an adult male) respectively. The left half of each histogram (marked by solid bars) denotes thelygenic females and their mates, while the right half of each histogram (marked by dashed bars) denotes arrhenogenic females and their mates longer assemblies than Trinity, Trinity recovered more genes. Figure 6 shows that ASplice had slightly higher specificity at the expense of lower sensitivity with respect to recovered mRNA positions, and had higher specificity with respect to alternative splicing junctions when the assembly conditions are more stringent. Sensitivity was low for all algorithms due to the small size of the data set. When compared to the results on the large D. melanogaster data set, this shows that as many libraries as possible should be utilized to obtain a more complete transcriptome, justifying the use of our scalable algorithm when possible.

Conclusions
We have developed an algorithm for de novo transcriptome assembly of non-model organisms that utilizes a large amount of RNA-Seq libraries in order to obtain a transcriptome that is as complete as possible, while simultaneously extracting alternative splicing information and expression levels in different libraries. When compared to existing algorithms, our algorithm is more conservative and generally has higher specificity at the expense of lower sensitivity, but is able to utilize larger amount of data to obtain more complete assemblies. As the size of the data set increases, larger values of k and usually much larger values of the k-mer coverage cutoff c are needed to obtain reasonable assemblies.
Since large libraries can be further subdivided into multiple parts during the parallel stage, our algorithm is scalable and the parallel stage can be run on different computing nodes. Since our k-mer counting technique requires iterating over every read and performing binary search repeatedly over increasing values of k, our algorithm is generally much more computationally intensive than existing algorithms, although our memory requirement is much lower for large data sets. Our strategy is especially suitable when multiple assemblies with different values of k are desired, as our algorithm is based on processing iteratively larger values of k. Incremental updates of assemblies are easy to perform, as it is only necessary to run the parallel stage on the new libraries before running the sequential stage on all libraries. Such a strategy is especially important to iteratively obtain more complete transcriptome assemblies over time through collaboration across research communities.
For small to medium data sets, it takes a few hours to a few days to complete the parallel stage over all values of k as long as the libraries are subdivided into small enough parts. The time to finish the sequential stage ranges from a few hours for a small data set to one to two days for a larger data set. For our largest C. rufifacies data set with 590 G bases, it takes a few days to finish the parallel stage over all values of k when the libraries are subdivided into parts with about 10 G bases each, with each computing  25_3  32277  748  46371  11311  26  50656  264  48079  9987  12  31949  305  30423  9560  9   25_5  21246  880  28366  9930  20  35797  254  34487  8412  12  20533  335  20061  8043  9   25_10  11509  842  15489  7075  18  20281  245  19759  5946  12  11243  338  11149  5722  9 For Oases and Trinity, the predicted unit is locus that contains a set of predicted transcripts, N50 denotes the N50 value of the longest transcript length in a predicted locus, and predicted transcripts of length at least 100 are retained. Other notations are the same as in Table 2  process consuming less than 14 GB memory. For each given setting of k and the k-mer coverage cutoff c, it takes a few days to finish the sequential stage. With values of c between 50 and 200 in our largest C. rufifacies data set, the memory consumption is less than 14 GB. Thus, our algorithm can assemble large data sets on a computing cluster with moderate resources.