Comparing de novo assemblers for 454 transcriptome data

L. sigmodontis transcriptome data

Comparison of assemblers

Numbers of reads used in assembly

The optimal assembler will use all the reads given, and will deliver assemblies with unambiguous mappings of reads to contigs, and thus putative transcripts. Each assembler has a different way of reporting the number of reads utilised. For example, MIRA reported only 2,170 singletons but classified 115,688 unassembled reads as 'debris'. Similarly, Newbler 2.3 and Newbler 2.5 generate separate lists of singletons, repetitive reads, and 'outliers' (problematic reads such as chimaeras). CAP3 and SeqMan only report assembled and unassembled reads. CLC does not track reads at all and maps reads back to the assembly to estimate where they might belong. Therefore, to compare the assemblies, we mapped all the reads back to each assembly using SSAHA2 [48] and the CLC reference aligner (which is part of the Assembly Cell suite [42]), as recommended in a study comparing short-read aligners [49]. SSAHA2's default settings for Roche 454 reads are clearly tuned for sensitivity, because the number of reads with multiple matches is very high compared to the CLC aligner. For both aligners, the highest numbers of reads were mapped back to the SeqMan assembly, while fewest were mapped to the Newbler 2.3 assembly. CAP3, CLC, and MIRA were comparable in terms of the number of reads used. CLC had the fewest reads with multiple matches by far, indicating that it was the least redundant of the six assemblies.

Number, mean length, and summed length of contigs

The optimal assembler will produce the longest summed length of contigs, while avoiding over-assembly of reads into in silico chimaeras, and avoiding the production of near-identical, largely overlapping contigs from allelic copies or error-rich data. It will also produce a transcriptome estimate with a mean and variance in contig length similar to that expected from the whole transcriptome. Newbler 2.3 generated an assembly with the largest N50 and the longest mean contig length (Table 4), but it also produced the smallest assembly span of the six programmes tested. Inspection of contig assemblies using the next-generation sequence assembly visualisation software Tablet [50] showed that Newbler 2.3 was discarding read overlap information in deriving the assembled contig sequence. The overall assembly sizes of CAP3, CLC and Newbler 2.3 are comparable to each other (between 14.4 and 16.7 Mb), as are the assembly sizes of MIRA, Newbler 2.5 and SeqMan (between 20.0 and 21.3 Mb). Newbler 2.5 has an assembly that is 39% larger than the Newbler 2.3 assembly, and seems to have solved the previous version's problem of discarding sequence data. Of all the six assemblers, Newbler 2.5 produced the highest number of contigs longer than 1 kb, and the most bases in contigs longer than 1 kb (these larger overall contig sizes are represented in Figure 1 as the assembly with the steepest initial slope). MIRA and SeqMan both generated comparable assembly spans, but with at least 8,000 more contigs than Newbler 2.5, indicating that they have shorter contigs overall.

Speed of analyses

The optimal assembler will complete analyses in a short period of time, and the time taken will scale well with increasing data volumes. While speed itself is not an overriding optimality criterion, rapid analyses will permit robust and extensive exploration of parameter space. In addition, some algorithms can make efficient use of multi-threaded processors or cluster computing, effectively reducing their wall-clock run time. The speed metric is not linked to assembly quality. CLC was astonishingly fast (taking only a few minutes to assemble 741,387 reads) compared to MIRA (which took 3 days in 'accurate' mode on the same hardware). This speed comes at a significant cost however, as CLC does not track read placement in the assembly. Tracking of read placement in an assembly is a very valuable feature, as it allows inspection of the data underpinning suspect assemblies, and thus CLC sacrifices speed for verifiability.

Assembly redundancy, shared and novel bases

In comparison to other assemblies, the optimal assembly will include the largest proportion of the unique bases present in the sum of all assemblies. Thus an assembler that produced a large assembly that included many contigs that were redundant could be worse that an assembly that was shorter but only included unique bases. To determine if the differing assembly sizes were due to novel sequences in each assembly, or just due to repetitive and redundant assemblies, we used BLAT [51] with default parameters to pairwise align all six assemblies (Figure 2).

These pairwise comparisons show that the SeqMan assembly had the most novel bases when compared to the others, and the others five assemblies had the fewest novel bases when compared to SeqMan. At the other end of the spectrum, the Newbler 2.3 assembly had the fewest novel bases and almost all its bases were seen in the other assemblies. For example, the SeqMan assembly had over 5 Mb of novel sequence compared to the assembly from Newbler 2.3, and the latter had only 23,366 novel bases in the reciprocal comparison. CAP3, CLC, MIRA, and Newbler 2.5 were all very similar in that they generated assemblies that had about 1-2 Mb of novel sequence compared to each other.

SeqMan, MIRA and Newbler 2.5 all generated assemblies in excess of 20 Mb, 25-50% larger than the 14-16 Mb assemblies produced by CAP3, CLC, and Newbler 2.3. However, these additional contigs and bases did not represent an excess of novel sequence, except in the case of SeqMan, suggesting that MIRA and Newbler 2.5 produce some redundant contigs (identical except for real polymorphisms or sequencing error). CLC generated the least redundant assembly, a direct consequence of using a de Bruijn graph algorithm with small fragments (21-mers) that collapses most repeats.

Alignments to reference sequences

Another optimality criterion for a novel de novo assembled transcriptome is how well it recapitulates previously determined sequences for the target species, and how well it represents sequences from related organisms. The best assembler will return contigs that match previous data well, and will deliver a high coverage of the conserved proteome of related taxa. We used four comparator datasets:

(i) a set of 2699 Sanger dideoxy sequenced ESTs from L. sigmodontis (Mark Blaxter, unpublished), assembled into 1602 clusters spanning 0.9 Mb (see http://www.nematodes.org/nembase4/species_info.php?species=LSC). We expect the Roche 454 pyrosequencing assemblies to faithfully reflect these longer reads and clusters;

(ii) the proteome predicted from the genome sequence of B. malayi [52]. We expect the B. malayi RefSeq protiens (11,472 GenBank entries, totalling 4.3 M amino acid residues) to be a good model for the proteome of L. sigmodontis, and thus proportional coverage of this proteome should reflect assembly quality;

(iii) gene families conserved across the phylum Nematoda derived from MCL-TRIBE [53] analysis of 62 species assembled in NEMBASE4 database (Mark Blaxter and Ben Elsworth, unpublished; see http://www.nematodes.org/nembase4/tribe_tree.php). We selected 3,681 tribes, containing 120,926 EST clusters, that had representatives in species spanning the inferred base of the phylum. They represent an estimate of the conserved core of genes that may be present in all nematode proteomes, and again proportional coverage should reflect assembly quality.

(iv) a set of 9,782 C. elegans proteins, grouped into 3,731 euKaryotic Ortholog Groups (KOGs), that are known to have orthologs in the proteomes of six other eukaryotic genomes [54]. We expect mRNAs for these conserved proteins to be present in the L. sigmodontis data, and thus proportional coverage should reflect assembly quality.

Merging assemblies to improve credibility

As each assembler uses a different algorithm to derive final contigs, and each of these algorithms may model different portions of the true transcriptome with different accuracies, we reasoned that combining assemblies might identify the subset of the assemblies that had high credibility, as it was found by both (or multiple) approaches. We checked if pairs of assemblies performed better than individual assemblies under the various optimality criteria proposed above, particularly in their recapitulation of the reference data. We combined two assemblies at a time by treating their (first-order) contigs as pseudo-reads and assembled these using a traditional OLC assembler (CAP3 with default settings). Only second-order contigs that contained first-order contigs from both constituent assemblies were considered for further analysis because two independent assemblers had agreed on the consensus sequence in that contig: we term these "robust contigs".

Although there was no consistent pattern indicating a clear superiority of any pair of assemblies (other than the exclusion of the Newbler 2.3 assembly), these merged assemblies show that it is possible to get more long, robust contigs that align better to reference sequences by combining the outputs of two different assemblies. Co-assembling three primary assemblies at a time and considering only the contigs that had 'reads' from all three assemblies gave smaller total contig numbers than pairwise co-assemblies. The number of stringent matches to reference sequences went up by less than 1% in all cases (see Additional File 2, Tables S2-S7), suggesting that the additional effort involved in generating and combining these assemblies may not be worthwhile.

Summary

We compared six assemblers by aligning their assembly contigs to four reference sequence sets. The Newbler 2.3 assembly scored worst in each comparison, probably because it had the smallest span of all the assemblies. Under stricter comparison parameters, the CLC assembly was poorer, suggesting that this assembly is more fragmented than the others. Although no single assembler was optimal in every case, Newbler 2.5 and SeqMan had the best alignments to related reference sequences overall. However, SeqMan generated over 8000 contigs extra for approximately the same total number of assembled bases (21.4 Mb versus 20.1 Mb).

We were able to improve the assembly by merging two assemblies at a time using a traditional OLC assembler (CAP3). This approach generated more high quality alignments to our reference sets.

Need for new assemblers

All the de novo assemblers in this study, with the exception of CLC, use the OLC assembly strategy. CLC uses de Bruijn graphs. As read lengths and throughput increase and sequencing costs come down, an average transcriptome assembly project for a non-model organism may comprise several million 800 base reads. Under the OLC paradigm, the computational time for assembling more reads rises exponentially with data complexity because the number of pairwise comparisons to detect overlaps will increase, and the layout graphs will be harder to resolve into a consensus sequence. The CLC de novo assembler is clearly a step in the right direction because its de Bruijn graph algorithm achieves reasonable results in very little time on large datasets. Other de Bruijn graph assemblers such as Oases [46] and ABySS [55] are promising but are currently not suitable for most non-normalised 454 pyrosequencing de novo transcriptome assembly projects because many transcripts have very low coverage that cannot be assembled reliably. Perhaps the way forward is to use the de Bruijn graph method for transcripts with high coverage and the OLC method for transcripts with low coverage. We have shown that an assembly-merging strategy delivers robust contigs from intermediate assemblies produced by current programs, and this strategy is likely to be of utility in deriving the best assemblies from future programs as well.

Library preparation

The L. sigmodontis samples were prepared by Stella Konstantinou, Eleana Theophilou and Simon Babayan, and sequenced by the GenePool, Edinburgh. Briefly, total RNA from three L. sigmodontis samples (adult male, adult female, microfilaria) was converted to double stranded cDNA using Evrogen's MINT cDNA synthesis kit. First strand cDNA was synthesised using reverse transcriptase (RT) from a 3'-primer with oligo(dT) sequence that annealed to the poly-A stretch of RNA and synthesised cDNA until the 5' end of the mRNA. Finally, double stranded cDNA synthesis was performed using PCR amplification, and the final product contained the same MINT adapter sequence at both 3' and 5' ends. The cDNA was fragmented, size-selected, library-prepped, and sequenced according to standard Roche-454 FLX and Titanium protocols. Sequence files have been deposited into the Sequence Read Archive (SRA) with accession number ERA011678 http://www.ebi.ac.uk/ena/data/view/ERA011678.

Read pre-processing

We chose to remove MINT adapters from the sequences ourselves rather than use the built-in adapter removal tools in MIRA, Newbler, SeqMan, and CLC, because each assembler pre-processes reads slightly differently and CAP3 does not remove adapters at all. The Roche 454 instrument outputs sequences in Standard Flowgram Format (SFF) files that come with quality and clipping information. To ensure that each assembler got exactly the same input, we changed the trim points directly in the SFF file to exclude adapter sequence rather than first convert the SFF files to fasta and quality files. The 29 base MINT adapter sequence was identified using BLASTN [56] from NCBI's new blast+ suite (Version 2.2.23) and by masking all matches with a bit-score greater than an empirically determined threshold of 25. The longest unmasked portion of the sequence was used to set new trim points in the SFF file.

Most EST pipelines remove poly(A/T) regions from reads because aligners would not be able to find good matches back to genomes that do not have the poly(A/T) sequence, and assemblers might try to align on these sequences. The former reason was not a factor for us as this was an unsequenced organism, but we did a trial run both with and without the poly(A/T) sequence and found no misassemblies on these regions by any of the assemblers. Newbler and MIRA account for such sequences when run in the -cdna and EST modes respectively. SeqMan and CAP3 remove low quality ends (homopolymer runs of poly(A/T) sequence have lower qualities in 454 sequencing), and CLC uses a de Bruijn graph which would not allow misassemblies because it would not be able to resolve the branching structure in a poly(A/T) region. Poly(A/T) sequence might also be indistinguishable from genomic sequences in some cases and therefore should not be removed. We found that the total amount of assembled sequence was about 10% lower for assemblies without poly(A/T) sequence. Therefore, we did not trim the poly(A/T) regions in the sequences. An additional advantage of leaving the poly(A/T) sequence in is that assembled contigs with a poly-A at the end (or poly-T at the start) indicate that the assembly end point is almost certainly correct. In our sequences, a poly(A/T) sequence typically had an N about half-way through because a degenerate oligo-dT primer was used to extract the poly-A tails from the cDNA. Any sequences with more than one N were discarded.

Assembly parameters

For each assembler, the version and parameters used are described in this section. CLC, Newbler, and SeqMan used the trimmed SFF files directly. CAP3 and MIRA used fasta and quality files extracted from the SFF files using Roche 454's sffinfo utility. All the assemblers returned a fasta file with contigs, and, except CLC, they all returned an ACE file with read placement information for each contig. All assemblers were run on an 8-core 3.0 GHz Linux workstation with 32 GB of memory, with the exception of SeqMan NGen, which currently runs only on Mac and Windows operating systems. SeqMan NGen was run on a dual core 2.53 GHz Mac mini server with 4 GB of memory.

The CAP3 binary for 64-bit Intel Linux systems was downloaded on June 20, 2010 from http://seq.cs.iastate.edu/cap3.html. The version was not specified in this download, and we did not change any of the default settings.

The CLC Assembly Cell command line application (version 3.02) has only one assembly-specific parameter, the minimum contig length reported, which we set to 100. It converts Roche 454 .sff files to fasta format, strips out all read identifiers, performs the assembly, and returns only a fasta file with contigs. As it does not keep track of reads, we had to use the CLC aligner, also known as the Reference Assembly application, to perform the additional step of mapping all the input reads back to the assembly contigs using default settings.

MIRA version 3.0.0 (production version) was used with the recommended quick switches for a Roche 454 EST assembly: -job = denovo,est,accurate,454. We used the -notraceinfo option because the traceinfo file provides clipping information, whereas we had already used the clipping information to create trimmed fasta and quality files. We also used -GE:not = 8 which is a general option for specifying the number of threads that should be used for steps that can use multiple cores. MIRA's final contig output file includes singletons but we looked for and removed all contigs with the "_s" prefix ("_c" indicates contigs with more than one read).

The latest releases of Roche 454's Newbler (version 2.3, 091027_1459, and version 2.5, newbler v2.5p1-internal-10Jun23-1) provide a "-cdna" option to assemble transcriptomes. In transcriptome assembly, the assembler frequently fails without the -cdna option because it expects approximately even coverage in genome assembly mode. The only other parameters we used were "-ace" (to generate an ACE file at the end) and "-cpu 8" to use all available cores. Like most other OLC assemblers, Newbler stops extending contigs when it cannot resolve branches in the overlap layout graph. Unlike other assemblers, Newbler then tries to create isotigs out of contigs that are consistently connected by a subset of reads. Each isotig corresponds to an alternative transcript, and any contigs or isotigs that share any read overlaps are put into the same isogroup. We noticed that several contig fragments in the Newbler 2.3 ACE file were not reported in the final assembly fasta file (454Isotigs.fna). Initially, the Newbler developers said that these contigs are not parts of actual transcripts and are therefore not reported, but we traced a few of these fragments and found that they were made up of well aligned, overlapping reads that had been split into contigs and then discarded. The developers have since modified Newbler and allowed us to try the pre-release version on our data that addresses this issue, using the -urt flag that allows contigs to be bridged across regions with only single-depth coverage.

DNAStar's SeqMan Ngen version 2 for Mac has a GUI that calls an underlying script. We chose the de novo transcriptome assembly option for 454 unpaired reads. When we tried a run with default settings and no quality or vector trimming, the assembly completed without any error messages but reported only 891 contigs. Subsequently we tried an assembly with quality trimming enabled which completed in 6 hours and reported a 21 Mb assembly.

Assembly read alignment

To get assembly stats on the number of reads used, we used SSAHA2 [48] and CLC's aligner to map all the input reads back to each assembly's contigs (as recommended in a comparison of next-generation aligners [49]). SSAHA2 was used with the -454 option, and the clc_ref_assemble_long binary from CLC's Assembly Cell distribution version 3.00.44070-beta2 was used with default settings.