BACKGROUND
with the recent and rapid advance of sequencing techniques, transcriptome sequencing  has emerged as a powerful tool for obtaining large amount of functional genomic data in non-model organisms.
this has encouraged efforts such as the one thousand plants project, or 1kp  <cit> , and many other transcriptome projects.
each of these data sets contains sequence information for thousands of genes, showing huge potential for phylogenomic inference.
however, there are many analytical and computational challenges that come with analyzing these data sets.
one of the biggest challenges is to accurately assemble the short reads from non-model organisms that do not have any reference genome, or de novo transcriptome assembly.
because this is the first step in any phylogenomic analysis, problems at this stage  cause difficulties for downstream analyses including ortholog and paralog identification, alignment, and matrix construction.
these problems increase the amount of missing data in the final aligned matrix, ultimately limiting the amount of useful transcriptomic data for phylogenomics.
currently, the most popular software packages for short-read rna-seq assembly include oases  <cit> , trinity  <cit> , trans-abyss  <cit>  and soapdenovo .
all four packages are based on constructing, simplifying, and resolving de bruijn graphs to extract likely transcripts .
two of these, oases  <cit>  and trans-abyss  <cit>  start with constructing de bruijn graphs directly from sequencing reads, remove potential errors, and then resolve each de bruijn graph to extract transcripts for each connected component  in the graph.
both packages use a range of k-mer sizes to accommodate variation in read coverages among genes.
trinity  <cit> , on the other hand, uses a single k-mer with size fixed at  <dig> bp.
trinity first carries out a greedy extension step starting from the most abundant k-mer to build linear contigs, groups overlapping contigs into connected components, and constructs a de bruijn graph for each component.
sequencing reads are then mapped to the graphs, the graphs are simplified, errors are removed , and finally likely isoforms are extracted for each component or subcomponent.
all four also use the information from mate pairs to assemble contigs into scaffolds when paired end reads are available.
each “locus” from the oases output  consists of one or more “transcripts”   <cit> .
biologically a locus or a component/subcomponent can each contain one gene or several paralogs, and a single gene can have fragments distributed among multiple loci or components/subcomponents  <cit> .
trans-abyss does not explicitly output sequences in hierarchical groups.
soapdenovo-trans is available only as precompiled executables without formal publication or source code  <cit> .
therefore we are unable to evaluate its assembly method in detail, and only include it here for completeness.
all the published de novo transcriptome assemblers are optimized for building references for comparing gene expression levels, identifying splice variants, and determining gene fusion events  <cit> .
for phylogenomic purposes, however, only one representative transcript for each gene is required.
splice variants are not only unutilized but also complicate the detection of true paralogs for phylogenomic matrix construction.
one method for picking the highest covered transcript from a locus from oases assemblies includes choosing the isoform with the highest geometric mean read coverage across nodes  <cit> .
by using geometric mean, regions of very low expression that may be associated with assembly error were penalized.
another way of selecting isoforms, as suggested by the trinity documentation  <cit> , is by mapping sequencing reads back to the assembled transcripts, and either picking the isoform with the highest coverage, or removing the isoform with the lowest coverage in a subcomponent.
however, strategies for picking the dominant and correct transcript have not been extensively explored.
a second problem that has been largely ignored by previous de novo transcriptome analyses is the creation of chimeras.
only “fusion transcripts” formed from multiple genes have been discussed in previous studies without further discussion of their nature and sources  <cit> .
chimeras can come from misassembly of short reads or pcr-induced recombination during library preparation.
chimeras may also be real biological products from gene fusion or trans-splicing.
in the case of short read assemblies, it is reasonable to assume that mis-assembly is the predominant cause, and by comparing different de novo assembly strategies, the difference is mainly caused by the assembly process.
in this study, by “chimera” we are referring to structural anomalies caused by all three potential sources.
if a chimera contains the same gene repeating itself, it is a self chimera.
if a chimera is the result of multiple genes, it is a multi-gene chimera.
if the sequences are assembled together in the same direction, it is considered to be a cis chimera.
if the sequences are assembled in opposite directions, it is a trans chimera.
considering both gene composition and direction, we recognized four types of chimerisms: cis-self, cis-multi-gene, trans-self and trans-multi-gene .
recognizing these four types of chimeras helps to identify the potential causes of misassemblies, and enables us to design effective post-assembly filters for removing chimeras.
here, we examine the extent of problems with chimeras as well as redundancy in de novo transcriptome assembly with the goal of optimizing transcript choice for phylogenomic analyses.
the castor bean  provides a unique opportunity for addressing these problems because it allows for direct comparison of a sanger sequencing-derived draft genome sequence  <cit>  with a typical de novo assembly of  <dig> bp paired end rna-seq data, recently generated as part of the 1kp project  <cit> .
here we compare assemblies from oases, trinity, trans-abyss and soapdenovo-trans to the genome-derived gene annotations of r. communis, to evaluate assembly strategies and post-assembly filters that best reduce redundancy while simultaneously maximizing total gene coverage and minimizing assembly chimeras.
methods
data set
short-read rna-seq results were obtained from the 1kp project database.
sampling and voucher information is available from the 1kp website   <cit> .
procedures of rna extraction, quality control, and library construction are detailed in johnson et al.
<cit> .
in summary, rna was obtained from mixed sample of leaves and flower buds.
sequencing was carried out on the illumina hiseq <dig> platform, with paired end  <dig> bp reads and insert size of around  <dig> bp.
non-redundant transcripts from genome annotation of r. communis were downloaded from phytozome v.  <dig>   <cit> .
r. communis chloroplast  and mitochondrial  genome sequences were downloaded from the genbank, and together they constitute the organellar genome sequences.
sequence cleanup and de novo assembly
downloaded reads were subjected to the following cleaning steps:  read pairs were removed if either of the reads in the read pair had average quality score ≤ 32;  read 3’ end nucleotides with quality scores <  <dig> were trimmed, and post-trim read pairs with either of the reads ≤  <dig> bp were removed; and  reads with adaptor contamination were removed.
possible chloroplast operons were observed when attempting de novo assembly with all reads, creating enormously complex graphs with thousands of assembled transcripts per locus in oases.
no mitochondrial operons were observed.
in order to reduce complexity and separate naturally occurring multi-gene transcripts from chimeras induced by assembly errors, cleaned read pairs that concordantly mapped to r. communis organellar genomes were removed using bowtie  <dig> v.  <dig> .0-beta <dig>  <cit> .
cleaned and filtered nuclear reads were assembled de novo using four software packages.
trinity v. 2012-02- <dig>  <cit>  was used with the default settings including a fixed k-mer size of  <dig> as suggested by the authors.
oases v.  <dig> .
<dig>  <cit>  was used with k-mer sizes from  <dig> to  <dig>  with an increment of  <dig> and an average insert length  <dig> bp.
the oases-m pipeline was used to merge assemblies from single k-mer oases assemblies either with all transcripts from k =  <dig> to  <dig> , or with a narrower range of k =  <dig>   <dig>   <dig>   <dig> and  <dig> .
abyss v.  <dig> .
<dig>  <cit>  was used with k-mer sizes  <dig>   <dig>   <dig>   <dig> and  <dig>  and contigs were filtered, extended and merged using trans-abyss v.  <dig> .
<dig>  <cit> .
soapdenovo-trans v.  <dig>   <cit>  was used with k-mer sizes  <dig>   <dig>   <dig>   <dig>   <dig> and  <dig> with default settings.
only assembled transcripts longer than  <dig> bp were kept.
evaluating de novo assembly by comparing to the reference transcripts
downloaded non-redundant reference transcripts from r. communis genome annotation were compared to the r. communis organellar genome sequences using default blat settings  <cit> .
transcripts with ≥ 80% combined region that mapped with ≥ 95% similarity to organellar genome sequences were removed.
the remaining nuclear transcripts were used as reference for accessing the quality of de novo assemblies of nuclear reads.
results from de novo assembly were mapped to the reference nuclear transcripts using default blat settings  <cit> .
blat hits were processed in three steps.
first, we ignored hits with lower than 95% similarity or less than  <dig> matching base pairs.
next, if an assembled transcript had only one blat hit, it was considered to be the best hit, and the assembled transcript was considered to be non-chimeric.
finally, if an assembled transcript had more than one blat hits, the best hit with the highest number of matching nucleotides was found.
the assembled transcript was identified as a chimera if any of the rest blat hits had no, or minimal, overlap in query coverage to the best hit .
if none of these comparisons identified chimeras, the longest hit was considered the best hit, and the assembled transcript was identified as non-chimeric.
selecting highly-covered isoforms for reducing redundancy
all four de novo assembly packages output large numbers of sequence isoforms, some of which are real splice variants, but many are chimeras.
for phylogenomic analysis, only a single representative, correct transcript is desired.
we investigated, in more detail, strategies to choose the representative isoforms in a multiple-k package oases, and in a single-k package trinity, since these two packages group assembled sequences into either locus or component/subcomponent explicitly in their output.
both also have relatively detailed documentations for their algorithms  <cit> .
in oases, very short sequence motifs in a locus tend to have high read coverage and need to be filtered out first.
we tested only considering transcripts longer than  <dig> ,  <dig> ,  <dig>  or  <dig>  of the longest transcripts in the same locus.
after filtering transcripts by their relative lengths, the remaining transcript with the highest read coverage as measured by the geometric mean of read coverages across nodes was chosen as the representative transcript.
we then calculated the percentage of the representative transcripts being chimeric vs. k-mer size, and number of transcripts per locus.
by examining percent chimera vs. relative length filter, k-mer size, and number of transcripts per locus plots, we applied a relative length filter of  <dig> , and chose the representative transcripts from k =  <dig>   <dig>   <dig>   <dig> only when there was  <dig> transcripts per locus, with no limit on number of transcripts per locus at k =  <dig> .
with this parameterization, gene families would only be represented at relatively high values of k. in addition of filter <dig>  we also tested retaining the representative transcripts from k =  <dig>   <dig>   <dig>   <dig> when there was either  <dig> or  <dig> transcripts per locus with no limit on number of transcripts per locus at k =  <dig> .
the program cap <dig>  <cit>  was then used to merge exemplars with overlap length cutoff of  <dig> and overlap percent identity cutoff of  <dig> .
we also tested using cd-hit-est  <cit>  instead of cap <dig> for removing redundancy with sequence identity cutoff set to  <dig>  .
to identify the representative isoforms in trinity, we mapped reads to the assembled transcripts using rsem  <cit> , and either only retained the isoforms with the highest isoform percentage  within each subcomponent , or removed the isoforms with the lowest isopct if there were more than one isoforms per subcomponent .
trans-abyss does not explicitly group assembled transcripts by graph component.
therefore we were unable to investigate strategies of choosing the dominant isoforms for each component other than by merging sequences by similarity.
trans-abyss multiple-k results were subject to cap <dig>  to combine similar sequences.
no multiple-k pipeline was provided with the current release  <dig>  of soapdenovo-trans.
to explore the effect of combining multiple k-mers in soapdenovo-trans, we combined all scaffolds from soapdenovo-trans k =  <dig>   <dig>   <dig>   <dig> and  <dig> using cap <dig> .
alternatively, we tested combining all contigs of k =  <dig>   <dig>   <dig>   <dig> and  <dig> instead of scaffolds, since the ns inserted in the scaffolds interfere with setting a similarity score cutoff for combining.
post-assembly trans chimera removal without reference
assembled sequences were blasted against non-redundant peptide sequences from  <dig> model eudicot species downloaded from the phytozome v.  <dig>  database  <cit> , excluding r. communis.
the closest relative to r. communis in the database was manihot esculenta, which split from the lineage leading to ricinus approximately  <dig> million years ago  <cit> .
in order to facilitate chimera identification, blastx was carried out with a relatively high e-value cutoff of  <dig> , and max_target_seqs set to  <dig>
trans chimeras were detected from blastx results with three steps.
first, high-scoring segment pairs  with less than 30% identity or with query coverages shorter than  <dig> bp were ignored.
second, trans-self chimeras were detected from hsps between a query-target pair.
detection of trans-self chimeras  included the following steps:
1) since hsps of the same direction were most often separated in the blastx results due to indels from sequencing errors, we considered all the query coverages in the same direction represented one single orf.
if all hsps from a query-target pair were in the same direction, we ignored all subsequent hsps and continued to the next query-target pair.
2) if hsps between a query-target pair were in the opposite directions, we merged all query coverages of plus direction, and merged all query coverages of minus direction.
3) if the overlap between the merged plus and minus query coverages were less than  <dig> bp and less than 20% of the shorter one of the two, the sequence was labeled as trans-self chimera, and the transcript was cut leaving only the longer of the merged query coverages.
all subsequent hsps from the same query were ignored.
finally, if no trans-self chimera was detected, we continued to check for trans-multi-gene chimeras  from all hsps of the same query against different targets.
we did this with the following steps:
1) all query coverages of plus direction were merged.
similarly, all query coverages of minus direction were merged.
2) if the overlap between the merged plus and minus query coverages from the same query were less than 20% of the shorter one of the two and less than  <dig> bp, the query sequence was cut leaving only the longer query coverage.
chimeras detected at this step are of the trans-multi-gene type.
this is because those that have different regions hitting the same target were detected as trans-self chimeras at the previous step and the query range being cut out.
in this case the trans-multi-gene detection steps were skipped.
only sequences that were longer than  <dig> bp after cutting were retained.
no cis chimeras were removed.
this is because when blasting to distantly related model species, tandem duplication, rearrangement, and heterogeneity in evolutionary rates among segments within a gene could result in the false detection of cis chimeras.
calculating total reference coverage and redundancy
for each assembly strategy tested, only the best blat hit from non-chimeric sequences was used for calculating total reference coverage.
for each nuclear reference transcript, only the longest non-chimeric reference coverage from blat hits was used to calculate the total coverage, as downstream phylogenomic analyses will only use one sequence per gene.
final reference coverage from each assembly strategy was calculated by percentage of total reference base pairs assembled to longer than  <dig> bp, percentage of reference genes assembled to longer than  <dig> bp, and percentage of reference genes assembled to 80% or more.
we measured redundancy as the number of assembled sequences with blat hits divided by the number of reference transcripts with blat hits.
all scripts used in the methods are available from bitbucket .
RESULTS
de novo assembly of short rna-seq reads recovered up to half of total genes
the rna-seq data set consisted of  <dig> , <dig> read pairs.
of these,  <dig> , <dig>  met our quality criteria and were free of adaptor contamination.
organellar reads accounted for a third of the cleaned reads, leaving  <dig> , <dig>  nuclear reads.
of all the cleaned and filtered nuclear reads,  <dig> % retained their original read length of  <dig> bp.
the remaining ranged from  <dig> to  <dig> bp in length .
out of the total  <dig>  annotated nuclear and organellar genes in r. communis,  <dig>   were nuclear.
comparison of de novo assembly of nuclear reads to reference nuclear transcripts shows that up to  <dig>   of these were assembled to longer than  <dig> bp with de novo methods.
this confirms that de novo assembly of short read rna-seq data is capable of recovering close to half of the genes in a nuclear genome with mixed tissue types of leaves and flower buds.
this number could be further increased with deeper sequencing depth or a higher diversity of tissue types.
such high gene coverage demonstrates the huge potential of rna-seq data in obtaining exome sequences in non-model organisms.
however, it also raises the question of why many phylogenomic analyses that use short-read rna-seq data only include hundreds of genes  <cit> , instead of thousands or even tens of thousands of genes, as in similar studies that incorporate longer reads from sanger or  <dig> sequencing  <cit> .
close scrutiny of assembly and post-assembly cleanups are the first steps towards increasing the matrix occupancy in phylogenomic analyses.
types and frequencies of chimeras vary dramatically among assembly strategies
our chimera detection criteria differ from previous studies in taking self chimeras into account and allowing short overlap between query coverages.
these detection criteria enabled us to thoroughly examine a wider range of chimeras.
among the four assembly packages, overall percentages of chimeras  were lowest among trans-abyss assemblies .
oases assemblies, both from single k-mer  and multiple k-mer , produced some of the highest percentages of chimeras.
the relationship between chimera rate and k-mer size varied among the three multiple-k assembly packages .
among single k-mer assemblies in oases, the percentage of chimeras was lowest when k =  <dig> , and increased towards both lower and higher k. the trend among trans-abyss single k-mer assemblies was opposite, with the highest percentage of chimera at k =  <dig> , and decreased towards both lower and higher k. as for soapdenovo-trans, percentages of chimeras decreased from  <dig> % when k =  <dig> to  <dig> % when k =  <dig>  both default multiple-k assembly pipelines by oases and trans-abyss produced higher percentages of chimeras as compared to single k-mer assemblies using the same software packages.
the composition of chimeras also varied dramatically among different assembly strategies.
the majority of chimeras produced by oases and trinity were trans-self , except in oases single-k assemblies when k was very large .
vast majority of all chimeras produced by soapdenovo-trans were either trans- or cis-multi-gene chimeras.
all four types of chimeras except cis-self had roughly equal share among assemblies using trans-abyss.
the difference in overall percentage of chimeras among assembly methods was largely contributed by the drastic difference in trans-self chimeras .
such drastic differences may relate to how contigs are assembled into scaffolds and how hairpin loops in the de bruijn graph are resolved  <cit> .
however, so far very little discussion is available on how loops in de bruijn graphs should be resolved, and whether different approaches should be taken between genome  vs. transcriptome assemblies .
the analyses presented here suggest a few general conclusions about transcriptomic assemblies.
1) the number and composition of chimeras differ dramatically among different assembly strategies, and the difference is largely due to the number of trans-self chimeras.
2) higher k-mer size does not necessarily lead to a smaller percentage of chimeras.
3) merging multiple single k-mer assemblies increases the percentage of chimeras.
with a chimera rate of at least 16%, oases-m, including both narrow and wide range of k, is not suitable for downstream analyses.
strategies for choosing highly covered isoforms to reduce redundancy in oases
another important consideration, especially for phylogenomics, is redundancy in assembly results.
each component/subcomponent in the trinity output, or each locus in oases, includes fragments of a gene or clusters of paralogous genes  <cit> , and often contains splice variants and/or chimeras.
in oases, very short transcripts that represent conserved motifs often have much higher read coverage compared to longer transcripts in the same locus.
these very short transcripts need to be filtered out before selecting exemplars.
we examined subsets of transcripts with lengths longer than  <dig> ,  <dig> ,  <dig> , and  <dig>  of the longest transcript within the same locus respectively.
we found that with a higher proportion cutoff, it was more likely that the dominant transcript was a chimera .
regardless of the value of k, loci that had high numbers of transcripts were likely to contain many chimeras.
the probability of the highest covered transcript being chimeric generally increased with size of the locus to as high as  <dig>  .
however, the probability remained relatively low  when size of locus was equal to  <dig> or equal to  <dig>  the spike of chimeras at two transcripts per locus  was almost exclusively caused by change in trans-self chimeras.
in oases, chimera rate decreased with increasing k until k =  <dig>  this was largely a result of the production of smaller loci with fewer numbers of transcripts per locus as k increased.
therefore, by only including loci with either one transcript , or one and three transcripts , the vast majority of chimeras were excluded.
to preserve members of gene families and alternatively-spliced genes, we included all transcripts at k =  <dig> regardless of number of transcripts per locus.
in this way, we only retained assemblies of these genes at a relatively high k-mer size when each locus had fewer transcripts, and therefore less likely to produce chimeras.
although we only investigated filtering loci by numbers of transcripts per locus and k-mer sizes in oases, similar methods are potentially useful for other software packages using multiple-k assembly strategies.
post-assembly trans chimera removal by blastx
after assembly, a blastx analysis of the transcripts against known protein sequences is a routine step for quality checking and downstream analyses such as homolog clustering and functional annotation.
here we take advantage of results from blastx to detect chimeras as well.
cis chimeras cannot be reliably detected when compared to sequences in a related species.
tandem duplication and rearrangement of gene segments can cause a false identification of cis-self chimera, and heterogeneity in base pair substitution rate within a gene can produce blastx hits similar to cis-multi-gene chimera.
trans chimeras, on the other hand, are much easier to detect from blastx results.
in the majority of eukaryotic nuclear genomes, a transcript is unlikely to have two different orfs of the opposite direction, especially if each of these orfs is highly similar to known coding sequences, of sufficient length, and there is no substantial overlap between these orfs.
our results regarding chimera removal were most successful in oases and trinity, both having trans chimeras as dominant chimera types.
the percentage of chimeras in the oases-mn assembly were reduced from  <dig> % to  <dig> % after using cap <dig> to merge similar reads and subsequent blastx and cutting of trans chimeras .
similarly, chimeras were reduced from  <dig> % to  <dig> % in oases k =  <dig>  and from  <dig> % to  <dig> % in trinity.
chimera rates were reduced from  <dig> % to  <dig> % in trans-abyss mk, and from  <dig> % to  <dig> % in soapdenovo-trans k =  <dig>  where trans chimeras account for a relatively small percentage of total chimeras.
in addition to using results from blast against other model species, one can also detect chimeras by mapping reads back to assembled transcripts to identify areas with a sudden change in read coverage.
however, this may not work for self chimeras.
methods have been developed for detecting fusion transcripts by mapping reads to a conspecific reference genome.
how effective it is doing so using transcripts assembled de novo has not yet been explored.
overall comparison of gene coverage, chimera and redundancy
composition of assembled transcripts, reference gene coverages, and redundancy are summarized in figure  <dig> and additional file  <dig>  the best assembly strategies are those that maximize number and length of reference covered, measured as percent gene covered, percent base pairs covered, and percent genes assembled to >80% , yet with a low chimera rate and redundancy .
overall, we found a trade-off between maximizing reference coverage vs. minimizing chimera rate and redundancy within each package before post processing using cap <dig> and blastx.
with proper post processing, however, results from seven assembly strategies among all four assembly packages converged to very similar values: all seven successfully assembled  <dig> – <dig> % of reference genes to more than  <dig> bp, with a chimera rate of  <dig> – <dig> %, and  <dig> – <dig>  transcripts per reference gene assembled .
the best score for each column, measured as highest in reference coverage and lowest in chimera rate and redundancy, is in bold.
among oases assemblies , while single k  assemblies suffered from low reference coverages, oases single k  and oases mn, mw assemblies all suffered from relatively high chimera rates and redundancies.
the problems with high chimera rates and redundancies persisted even after using cap <dig> to reduce redundancy and blastx to remove trans chimeras in two of the strategies with relatively high reference coverages .
we also tested using cd-hit-est instead of cap <dig> to remove redundancy, and the outcome was very similar to cap <dig> results .
the default oases-m pipeline takes transcripts from single-k assemblies as input and assembles them using k =  <dig>  this strategy is one step removed from the original read coverage information, and keeps most isoforms and errors from the original single k-mer assemblies.
our two novel oases mk pipelines “oases filter1-cap3-blast” and “oases filter1&3-cap3-blast” utilize the read coverage and number of transcripts per locus information in single k-mer assemblies, and achieve both low redundancy and low chimera rates with high reference coverages .
all trans-abyss single k-mer assemblies suffered from relatively low reference coverages , whereas trans-abyss mk, especially after post processing, ranks among one of the highest performing assembly strategies tested .
however, trans-abyss, as a package, carries out many functions beyond sequence assembly without detailed documentation.
the current version  <dig> .
<dig> is not as user-friendly as the other packages tested in this study.
trinity, trinity blast, and trinity cap3-blast are the top three assembly strategies in terms of total reference genes assembled .
between trinity blast and trinity cap3-blast, adding a cap <dig> step reduced redundancy from  <dig>  to  <dig> , while slightly reduced reference coverage.
an alternative approach of reducing redundancy by removing the lowest covered transcript in a subcomponent  reduced both redundancy and chimera rate only slightly, at a cost of slightly reducing reference coverage as well.
this suggests that the lowest covered transcripts can be either chimeras or paralogous genes that have relatively low expression levels.
similarly, only keeping the highest covered transcripts for each subcomponent  further reduced reference coverage, redundancy, and chimera rate.
this strategy is a more aggressive way of reducing redundancy than cap <dig> alone or removel.
however, since trinity started with the highest reference coverage by the number of genes, after pickh combined with cap <dig> and blastx, this very aggressive post processing strategy produced one of the cleanest, lease redundant assemblies, and the overall reference coverage is not very far below the other high performing assembly strategies .
soapdenovo-trans k=21-cap3-blast produced the highest reference coverage measured by base pairs only after oases mw .
it does this by aggressively scaffolding from contigs, and this process more than doubled the percentages of chimeras .
the ns inserted in many of the scaffolds also made it difficult to find a suitable similarity cutoff for merging scaffolds from multiple single k-mer soapdenovo-trans assemblies using cap <dig>  instead of using scaffolds, we merged all contigs from soapdenovo-trans k =  <dig>   <dig>   <dig>   <dig> and  <dig> using cap <dig>  this strategy, soapdenovo-trans contigs cap3-blast, turns out to be among the highest performing strategies tested , comparable to or outperforms most strategies incorporating mate pair information.
this suggests that with more accurate scaffolding, soapdenovo-trans has the potential to further improve its performance.
in summary, in our particular data set we identified seven strategies among four assembly packages that produced results that are most suitable for downstream phylogenomic analysis.
all seven scored very close to each other in reference coverage, chimera rate, and redundancy, suggesting that the most suitable assembly strategy can vary with different data sets, and that there is no single best strategy for all assembly tasks.
CONCLUSIONS
de novo assembly of short read rna-seq data is capable of recovering up to half of the total expressed genes to more than  <dig> bp.
however, oases, trinity, and to a lesser extent trans-abyss, all produce large amounts of trans-self chimeras.
we find that although oases-m produces the highest gene coverage among popular assembly packages, its high redundancy and chimera rate make it unsuitable for phylogenomic purposes, even with extensive post processing.
trinity achieves high reference coverage that is similar to oases-m, but with a much lower redundancy and chimera rate, especially after trans chimera removal using blastx and effective redundancy reduction.
soapdenovo-trans, although yet unpublished, is highly promising in producing some of the cleanest assemblies with the highest reference coverage.
many packages used in next generation sequence analyses were initially optimized for purposes other than phylogenomics, and care should be taken when utilizing these packages for phylogenomics.
with the rapid development of both sequencing techniques and software packages, one needs to examine the specific types and structures of assembly problems in order to minimize them.
future evaluation of de novo assembly tools for phylogenomics should focus on completeness  <cit> , chimerism, and redundancy, instead of ambiguous measures borrowed from genome assembly, such as n <dig>  lastly, since a significant number of the chimeras are trans-chimeras, strand-specific library preparation and assembly, as demonstrated by garg et al.
<cit> , can eliminate a major source of assembly errors.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
yy and sas designed the study, carried out the analyses, drafted the manuscripts, and read and approved the final manuscript.
supplementary material
additional file 1
length distribution of cleaned and filtered nuclear reads.
click here for file
additional file 2
chimera compositions among assembled transcripts before post-processing.
click here for file
additional file 3
reference coverages, chimera percentages and redundancies among assembly strategies.
click here for file
acknowledgements
the authors would like to thank mark chase, gane ka-shu wong, pamela s. soltis and douglas e. soltis for granting us permission to use the 1kp data; michael j. moore, samuel f. brockington, joseph w. brown and cody hinchliff for helpful discussions and comments on the draft; and two anonymous reviews for help improving the manuscript.
the department of ecology and evolutionary biology at university of michigan, ann arbor, provided funding.

 with the recent and rapid advance of sequencing techniques, transcriptome sequencing has emerged as a powerful tool for obtaining large amount of functional genomic data in non-model organisms. this has encouraged efforts such as the one thousand plants project, or 1kp, and many other transcriptome projects. each of these data sets contains sequence information for thousands of genes, showing huge potential for phylogenomic inference. however, there are many analytical and computational challenges that come with analyzing these data sets. one of the biggest challenges is to accurately assemble the short reads from non-model organisms that do not have any reference genome, or de novo transcriptome assembly. because this is the first step in any phylogenomic analysis, problems at this stage cause difficulties for downstream analyses including ortholog and paralog identification, alignment, and matrix construction. these problems increase the amount of missing data in the final aligned matrix, ultimately limiting the amount of useful transcriptomic data for phylogenomics. currently, the most popular software packages for short-read rna-seq assembly include oases, trinity, trans-abyss and soapdenovo. all four packages are based on constructing, simplifying, and resolving de bruijn graphs to extract likely transcripts. two of these, oases and trans-abyss start with constructing de bruijn graphs directly from sequencing reads, remove potential errors, and then resolve each de bruijn graph to extract transcripts for each connected component in the graph. both packages use a range of k-mer sizes to accommodate variation in read coverages among genes. trinity, on the other hand, uses a single k-mer with size fixed at bp. trinity first carries out a greedy extension step starting from the most abundant k-mer to build linear contigs, groups overlapping contigs into connected components, and constructs a de bruijn graph for each component. sequencing reads are then mapped to the graphs, the graphs are simplified, errors are removed, and finally likely isoforms are extracted for each component or subcomponent. all four also use the information from mate pairs to assemble contigs into scaffolds when paired end reads are available. each “locus” from the oases output consists of one or more “transcripts”. biologically a locus or a component/subcomponent can each contain one gene or several paralogs, and a single gene can have fragments distributed among multiple loci or components/subcomponents. trans-abyss does not explicitly output sequences in hierarchical groups. soapdenovo-trans is available only as precompiled executables without formal publication or source code. therefore we are unable to evaluate its assembly method in detail, and only include it here for completeness. all the published de novo transcriptome assemblers are optimized for building references for comparing gene expression levels, identifying splice variants, and determining gene fusion events. for phylogenomic purposes, however, only one representative transcript for each gene is required. splice variants are not only unutilized but also complicate the detection of true paralogs for phylogenomic matrix construction. one method for picking the highest covered transcript from a locus from oases assemblies includes choosing the isoform with the highest geometric mean read coverage across nodes. by using geometric mean, regions of very low expression that may be associated with assembly error were penalized. another way of selecting isoforms, as suggested by the trinity documentation, is by mapping sequencing reads back to the assembled transcripts, and either picking the isoform with the highest coverage, or removing the isoform with the lowest coverage in a subcomponent. however, strategies for picking the dominant and correct transcript have not been extensively explored. a second problem that has been largely ignored by previous de novo transcriptome analyses is the creation of chimeras. only “fusion transcripts” formed from multiple genes have been discussed in previous studies without further discussion of their nature and sources. chimeras can come from misassembly of short reads or pcr-induced recombination during library preparation. chimeras may also be real biological products from gene fusion or trans-splicing. in the case of short read assemblies, it is reasonable to assume that mis-assembly is the predominant cause, and by comparing different de novo assembly strategies, the difference is mainly caused by the assembly process. in this study, by “chimera” we are referring to structural anomalies caused by all three potential sources. if a chimera contains the same gene repeating itself, it is a self chimera. if a chimera is the result of multiple genes, it is a multi-gene chimera. if the sequences are assembled together in the same direction, it is considered to be a cis chimera. if the sequences are assembled in opposite directions, it is a trans chimera. considering both gene composition and direction, we recognized four types of chimerisms: cis-self, cis-multi-gene, trans-self and trans-multi-gene. recognizing these four types of chimeras helps to identify the potential causes of misassemblies, and enables us to design effective post-assembly filters for removing chimeras. here, we examine the extent of problems with chimeras as well as redundancy in de novo transcriptome assembly with the goal of optimizing transcript choice for phylogenomic analyses. the castor bean provides a unique opportunity for addressing these problems because it allows for direct comparison of a sanger sequencing-derived draft genome sequence with a typical de novo assembly of bp paired end rna-seq data, recently generated as part of the 1kp project. here we compare assemblies from oases, trinity, trans-abyss and soapdenovo-trans to the genome-derived gene annotations of r. communis, to evaluate assembly strategies and post-assembly filters that best reduce redundancy while simultaneously maximizing total gene coverage and minimizing assembly chimeras. methods data set short-read rna-seq results were obtained from the 1kp project database. sampling and voucher information is available from the 1kp website. procedures of rna extraction, quality control, and library construction are detailed in johnson et al.. in summary, rna was obtained from mixed sample of leaves and flower buds. sequencing was carried out on the illumina hiseq platform, with paired end bp reads and insert size of around bp. non-redundant transcripts from genome annotation of r. communis were downloaded from phytozome v.. r. communis chloroplast and mitochondrial genome sequences were downloaded from the genbank, and together they constitute the organellar genome sequences. sequence cleanup and de novo assembly downloaded reads were subjected to the following cleaning steps: read pairs were removed if either of the reads in the read pair had average quality score ≤ 32; read 3’ end nucleotides with quality scores < were trimmed, and post-trim read pairs with either of the reads ≤ bp were removed; and reads with adaptor contamination were removed. possible chloroplast operons were observed when attempting de novo assembly with all reads, creating enormously complex graphs with thousands of assembled transcripts per locus in oases. no mitochondrial operons were observed. in order to reduce complexity and separate naturally occurring multi-gene transcripts from chimeras induced by assembly errors, cleaned read pairs that concordantly mapped to r. communis organellar genomes were removed using bowtie v..0-beta. cleaned and filtered nuclear reads were assembled de novo using four software packages. trinity v. 2012-02- was used with the default settings including a fixed k-mer size of as suggested by the authors. oases v.. was used with k-mer sizes from to with an increment of and an average insert length bp. the oases-m pipeline was used to merge assemblies from single k-mer oases assemblies either with all transcripts from k = to, or with a narrower range of k = and. abyss v.. was used with k-mer sizes and and contigs were filtered, extended and merged using trans-abyss v... soapdenovo-trans v. was used with k-mer sizes and with default settings. only assembled transcripts longer than bp were kept. evaluating de novo assembly by comparing to the reference transcripts downloaded non-redundant reference transcripts from r. communis genome annotation were compared to the r. communis organellar genome sequences using default blat settings. transcripts with ≥ 80% combined region that mapped with ≥ 95% similarity to organellar genome sequences were removed. the remaining nuclear transcripts were used as reference for accessing the quality of de novo assemblies of nuclear reads. results from de novo assembly were mapped to the reference nuclear transcripts using default blat settings. blat hits were processed in three steps. first, we ignored hits with lower than 95% similarity or less than matching base pairs. next, if an assembled transcript had only one blat hit, it was considered to be the best hit, and the assembled transcript was considered to be non-chimeric. finally, if an assembled transcript had more than one blat hits, the best hit with the highest number of matching nucleotides was found. the assembled transcript was identified as a chimera if any of the rest blat hits had no, or minimal, overlap in query coverage to the best hit. if none of these comparisons identified chimeras, the longest hit was considered the best hit, and the assembled transcript was identified as non-chimeric. selecting highly-covered isoforms for reducing redundancy all four de novo assembly packages output large numbers of sequence isoforms, some of which are real splice variants, but many are chimeras. for phylogenomic analysis, only a single representative, correct transcript is desired. we investigated, in more detail, strategies to choose the representative isoforms in a multiple-k package oases, and in a single-k package trinity, since these two packages group assembled sequences into either locus or component/subcomponent explicitly in their output. both also have relatively detailed documentations for their algorithms. in oases, very short sequence motifs in a locus tend to have high read coverage and need to be filtered out first. we tested only considering transcripts longer than,, or of the longest transcripts in the same locus. after filtering transcripts by their relative lengths, the remaining transcript with the highest read coverage as measured by the geometric mean of read coverages across nodes was chosen as the representative transcript. we then calculated the percentage of the representative transcripts being chimeric vs. k-mer size, and number of transcripts per locus. by examining percent chimera vs. relative length filter, k-mer size, and number of transcripts per locus plots, we applied a relative length filter of, and chose the representative transcripts from k = only when there was transcripts per locus, with no limit on number of transcripts per locus at k =. with this parameterization, gene families would only be represented at relatively high values of k. in addition of filter we also tested retaining the representative transcripts from k = when there was either or transcripts per locus with no limit on number of transcripts per locus at k =. the program cap was then used to merge exemplars with overlap length cutoff of and overlap percent identity cutoff of. we also tested using cd-hit-est instead of cap for removing redundancy with sequence identity cutoff set to. to identify the representative isoforms in trinity, we mapped reads to the assembled transcripts using rsem, and either only retained the isoforms with the highest isoform percentage within each subcomponent, or removed the isoforms with the lowest isopct if there were more than one isoforms per subcomponent. trans-abyss does not explicitly group assembled transcripts by graph component. therefore we were unable to investigate strategies of choosing the dominant isoforms for each component other than by merging sequences by similarity. trans-abyss multiple-k results were subject to cap to combine similar sequences. no multiple-k pipeline was provided with the current release of soapdenovo-trans. to explore the effect of combining multiple k-mers in soapdenovo-trans, we combined all scaffolds from soapdenovo-trans k = and using cap. alternatively, we tested combining all contigs of k = and instead of scaffolds, since the ns inserted in the scaffolds interfere with setting a similarity score cutoff for combining. post-assembly trans chimera removal without reference assembled sequences were blasted against non-redundant peptide sequences from model eudicot species downloaded from the phytozome v. database, excluding r. communis. the closest relative to r. communis in the database was manihot esculenta, which split from the lineage leading to ricinus approximately million years ago. in order to facilitate chimera identification, blastx was carried out with a relatively high e-value cutoff of, and max_target_seqs set to trans chimeras were detected from blastx results with three steps. first, high-scoring segment pairs with less than 30% identity or with query coverages shorter than bp were ignored. second, trans-self chimeras were detected from hsps between a query-target pair. detection of trans-self chimeras included the following steps: 1) since hsps of the same direction were most often separated in the blastx results due to indels from sequencing errors, we considered all the query coverages in the same direction represented one single orf. if all hsps from a query-target pair were in the same direction, we ignored all subsequent hsps and continued to the next query-target pair. 2) if hsps between a query-target pair were in the opposite directions, we merged all query coverages of plus direction, and merged all query coverages of minus direction. 3) if the overlap between the merged plus and minus query coverages were less than bp and less than 20% of the shorter one of the two, the sequence was labeled as trans-self chimera, and the transcript was cut leaving only the longer of the merged query coverages. all subsequent hsps from the same query were ignored. finally, if no trans-self chimera was detected, we continued to check for trans-multi-gene chimeras from all hsps of the same query against different targets. we did this with the following steps: 1) all query coverages of plus direction were merged. similarly, all query coverages of minus direction were merged. 2) if the overlap between the merged plus and minus query coverages from the same query were less than 20% of the shorter one of the two and less than bp, the query sequence was cut leaving only the longer query coverage. chimeras detected at this step are of the trans-multi-gene type. this is because those that have different regions hitting the same target were detected as trans-self chimeras at the previous step and the query range being cut out. in this case the trans-multi-gene detection steps were skipped. only sequences that were longer than bp after cutting were retained. no cis chimeras were removed. this is because when blasting to distantly related model species, tandem duplication, rearrangement, and heterogeneity in evolutionary rates among segments within a gene could result in the false detection of cis chimeras. calculating total reference coverage and redundancy for each assembly strategy tested, only the best blat hit from non-chimeric sequences was used for calculating total reference coverage. for each nuclear reference transcript, only the longest non-chimeric reference coverage from blat hits was used to calculate the total coverage, as downstream phylogenomic analyses will only use one sequence per gene. final reference coverage from each assembly strategy was calculated by percentage of total reference base pairs assembled to longer than bp, percentage of reference genes assembled to longer than bp, and percentage of reference genes assembled to 80% or more. we measured redundancy as the number of assembled sequences with blat hits divided by the number of reference transcripts with blat hits. all scripts used in the methods are available from bitbucket.  de novo assembly of short rna-seq reads recovered up to half of total genes the rna-seq data set consisted of, read pairs. of these,, met our quality criteria and were free of adaptor contamination. organellar reads accounted for a third of the cleaned reads, leaving, nuclear reads. of all the cleaned and filtered nuclear reads, % retained their original read length of bp. the remaining ranged from to bp in length. out of the total annotated nuclear and organellar genes in r. communis, were nuclear. comparison of de novo assembly of nuclear reads to reference nuclear transcripts shows that up to of these were assembled to longer than bp with de novo methods. this confirms that de novo assembly of short read rna-seq data is capable of recovering close to half of the genes in a nuclear genome with mixed tissue types of leaves and flower buds. this number could be further increased with deeper sequencing depth or a higher diversity of tissue types. such high gene coverage demonstrates the huge potential of rna-seq data in obtaining exome sequences in non-model organisms. however, it also raises the question of why many phylogenomic analyses that use short-read rna-seq data only include hundreds of genes, instead of thousands or even tens of thousands of genes, as in similar studies that incorporate longer reads from sanger or sequencing. close scrutiny of assembly and post-assembly cleanups are the first steps towards increasing the matrix occupancy in phylogenomic analyses. types and frequencies of chimeras vary dramatically among assembly strategies our chimera detection criteria differ from previous studies in taking self chimeras into account and allowing short overlap between query coverages. these detection criteria enabled us to thoroughly examine a wider range of chimeras. among the four assembly packages, overall percentages of chimeras were lowest among trans-abyss assemblies. oases assemblies, both from single k-mer and multiple k-mer, produced some of the highest percentages of chimeras. the relationship between chimera rate and k-mer size varied among the three multiple-k assembly packages. among single k-mer assemblies in oases, the percentage of chimeras was lowest when k =, and increased towards both lower and higher k. the trend among trans-abyss single k-mer assemblies was opposite, with the highest percentage of chimera at k =, and decreased towards both lower and higher k. as for soapdenovo-trans, percentages of chimeras decreased from % when k = to % when k = both default multiple-k assembly pipelines by oases and trans-abyss produced higher percentages of chimeras as compared to single k-mer assemblies using the same software packages. the composition of chimeras also varied dramatically among different assembly strategies. the majority of chimeras produced by oases and trinity were trans-self, except in oases single-k assemblies when k was very large. vast majority of all chimeras produced by soapdenovo-trans were either trans- or cis-multi-gene chimeras. all four types of chimeras except cis-self had roughly equal share among assemblies using trans-abyss. the difference in overall percentage of chimeras among assembly methods was largely contributed by the drastic difference in trans-self chimeras. such drastic differences may relate to how contigs are assembled into scaffolds and how hairpin loops in the de bruijn graph are resolved. however, so far very little discussion is available on how loops in de bruijn graphs should be resolved, and whether different approaches should be taken between genome vs. transcriptome assemblies. the analyses presented here suggest a few general conclusions about transcriptomic assemblies. 1) the number and composition of chimeras differ dramatically among different assembly strategies, and the difference is largely due to the number of trans-self chimeras. 2) higher k-mer size does not necessarily lead to a smaller percentage of chimeras. 3) merging multiple single k-mer assemblies increases the percentage of chimeras. with a chimera rate of at least 16%, oases-m, including both narrow and wide range of k, is not suitable for downstream analyses. strategies for choosing highly covered isoforms to reduce redundancy in oases another important consideration, especially for phylogenomics, is redundancy in assembly results. each component/subcomponent in the trinity output, or each locus in oases, includes fragments of a gene or clusters of paralogous genes, and often contains splice variants and/or chimeras. in oases, very short transcripts that represent conserved motifs often have much higher read coverage compared to longer transcripts in the same locus. these very short transcripts need to be filtered out before selecting exemplars. we examined subsets of transcripts with lengths longer than,,, and of the longest transcript within the same locus respectively. we found that with a higher proportion cutoff, it was more likely that the dominant transcript was a chimera. regardless of the value of k, loci that had high numbers of transcripts were likely to contain many chimeras. the probability of the highest covered transcript being chimeric generally increased with size of the locus to as high as. however, the probability remained relatively low when size of locus was equal to or equal to the spike of chimeras at two transcripts per locus was almost exclusively caused by change in trans-self chimeras. in oases, chimera rate decreased with increasing k until k = this was largely a result of the production of smaller loci with fewer numbers of transcripts per locus as k increased. therefore, by only including loci with either one transcript, or one and three transcripts, the vast majority of chimeras were excluded. to preserve members of gene families and alternatively-spliced genes, we included all transcripts at k = regardless of number of transcripts per locus. in this way, we only retained assemblies of these genes at a relatively high k-mer size when each locus had fewer transcripts, and therefore less likely to produce chimeras. although we only investigated filtering loci by numbers of transcripts per locus and k-mer sizes in oases, similar methods are potentially useful for other software packages using multiple-k assembly strategies. post-assembly trans chimera removal by blastx after assembly, a blastx analysis of the transcripts against known protein sequences is a routine step for quality checking and downstream analyses such as homolog clustering and functional annotation. here we take advantage of results from blastx to detect chimeras as well. cis chimeras cannot be reliably detected when compared to sequences in a related species. tandem duplication and rearrangement of gene segments can cause a false identification of cis-self chimera, and heterogeneity in base pair substitution rate within a gene can produce blastx hits similar to cis-multi-gene chimera. trans chimeras, on the other hand, are much easier to detect from blastx results. in the majority of eukaryotic nuclear genomes, a transcript is unlikely to have two different orfs of the opposite direction, especially if each of these orfs is highly similar to known coding sequences, of sufficient length, and there is no substantial overlap between these orfs. our results regarding chimera removal were most successful in oases and trinity, both having trans chimeras as dominant chimera types. the percentage of chimeras in the oases-mn assembly were reduced from % to % after using cap to merge similar reads and subsequent blastx and cutting of trans chimeras. similarly, chimeras were reduced from % to % in oases k = and from % to % in trinity. chimera rates were reduced from % to % in trans-abyss mk, and from % to % in soapdenovo-trans k = where trans chimeras account for a relatively small percentage of total chimeras. in addition to using results from blast against other model species, one can also detect chimeras by mapping reads back to assembled transcripts to identify areas with a sudden change in read coverage. however, this may not work for self chimeras. methods have been developed for detecting fusion transcripts by mapping reads to a conspecific reference genome. how effective it is doing so using transcripts assembled de novo has not yet been explored. overall comparison of gene coverage, chimera and redundancy composition of assembled transcripts, reference gene coverages, and redundancy are summarized in figure  and additional file the best assembly strategies are those that maximize number and length of reference covered, measured as percent gene covered, percent base pairs covered, and percent genes assembled to >80%, yet with a low chimera rate and redundancy. overall, we found a trade-off between maximizing reference coverage vs. minimizing chimera rate and redundancy within each package before post processing using cap and blastx. with proper post processing, however, results from seven assembly strategies among all four assembly packages converged to very similar values: all seven successfully assembled – % of reference genes to more than bp, with a chimera rate of – %, and – transcripts per reference gene assembled. the best score for each column, measured as highest in reference coverage and lowest in chimera rate and redundancy, is in bold. among oases assemblies, while single k assemblies suffered from low reference coverages, oases single k and oases mn, mw assemblies all suffered from relatively high chimera rates and redundancies. the problems with high chimera rates and redundancies persisted even after using cap to reduce redundancy and blastx to remove trans chimeras in two of the strategies with relatively high reference coverages. we also tested using cd-hit-est instead of cap to remove redundancy, and the outcome was very similar to cap results. the default oases-m pipeline takes transcripts from single-k assemblies as input and assembles them using k = this strategy is one step removed from the original read coverage information, and keeps most isoforms and errors from the original single k-mer assemblies. our two novel oases mk pipelines “oases filter1-cap3-blast” and “oases filter1&3-cap3-blast” utilize the read coverage and number of transcripts per locus information in single k-mer assemblies, and achieve both low redundancy and low chimera rates with high reference coverages. all trans-abyss single k-mer assemblies suffered from relatively low reference coverages, whereas trans-abyss mk, especially after post processing, ranks among one of the highest performing assembly strategies tested. however, trans-abyss, as a package, carries out many functions beyond sequence assembly without detailed documentation. the current version. is not as user-friendly as the other packages tested in this study. trinity, trinity blast, and trinity cap3-blast are the top three assembly strategies in terms of total reference genes assembled. between trinity blast and trinity cap3-blast, adding a cap step reduced redundancy from to, while slightly reduced reference coverage. an alternative approach of reducing redundancy by removing the lowest covered transcript in a subcomponent reduced both redundancy and chimera rate only slightly, at a cost of slightly reducing reference coverage as well. this suggests that the lowest covered transcripts can be either chimeras or paralogous genes that have relatively low expression levels. similarly, only keeping the highest covered transcripts for each subcomponent further reduced reference coverage, redundancy, and chimera rate. this strategy is a more aggressive way of reducing redundancy than cap alone or removel. however, since trinity started with the highest reference coverage by the number of genes, after pickh combined with cap and blastx, this very aggressive post processing strategy produced one of the cleanest, lease redundant assemblies, and the overall reference coverage is not very far below the other high performing assembly strategies. soapdenovo-trans k=21-cap3-blast produced the highest reference coverage measured by base pairs only after oases mw. it does this by aggressively scaffolding from contigs, and this process more than doubled the percentages of chimeras. the ns inserted in many of the scaffolds also made it difficult to find a suitable similarity cutoff for merging scaffolds from multiple single k-mer soapdenovo-trans assemblies using cap instead of using scaffolds, we merged all contigs from soapdenovo-trans k = and using cap this strategy, soapdenovo-trans contigs cap3-blast, turns out to be among the highest performing strategies tested, comparable to or outperforms most strategies incorporating mate pair information. this suggests that with more accurate scaffolding, soapdenovo-trans has the potential to further improve its performance. in summary, in our particular data set we identified seven strategies among four assembly packages that produced results that are most suitable for downstream phylogenomic analysis. all seven scored very close to each other in reference coverage, chimera rate, and redundancy, suggesting that the most suitable assembly strategy can vary with different data sets, and that there is no single best strategy for all assembly tasks.  de novo assembly of short read rna-seq data is capable of recovering up to half of the total expressed genes to more than bp. however, oases, trinity, and to a lesser extent trans-abyss, all produce large amounts of trans-self chimeras. we find that although oases-m produces the highest gene coverage among popular assembly packages, its high redundancy and chimera rate make it unsuitable for phylogenomic purposes, even with extensive post processing. trinity achieves high reference coverage that is similar to oases-m, but with a much lower redundancy and chimera rate, especially after trans chimera removal using blastx and effective redundancy reduction. soapdenovo-trans, although yet unpublished, is highly promising in producing some of the cleanest assemblies with the highest reference coverage. many packages used in next generation sequence analyses were initially optimized for purposes other than phylogenomics, and care should be taken when utilizing these packages for phylogenomics. with the rapid development of both sequencing techniques and software packages, one needs to examine the specific types and structures of assembly problems in order to minimize them. future evaluation of de novo assembly tools for phylogenomics should focus on completeness, chimerism, and redundancy, instead of ambiguous measures borrowed from genome assembly, such as n lastly, since a significant number of the chimeras are trans-chimeras, strand-specific library preparation and assembly, as demonstrated by garg et al., can eliminate a major source of assembly errors. competing interests the authors declare that they have no competing interests. authors’ contributions yy and sas designed the study, carried out the analyses, drafted the manuscripts, and read and approved the final manuscript. supplementary material additional file 1 length distribution of cleaned and filtered nuclear reads. click here for file additional file 2 chimera compositions among assembled transcripts before post-processing. click here for file additional file 3 reference coverages, chimera percentages and redundancies among assembly strategies. click here for file acknowledgements the authors would like to thank mark chase, gane ka-shu wong, pamela s. soltis and douglas e. soltis for granting us permission to use the 1kp data; michael j. moore, samuel f. brockington, joseph w. brown and cody hinchliff for helpful discussions and comments on the draft; and two anonymous reviews for help improving the manuscript. the department of ecology and evolutionary biology at university of michigan, ann arbor, provided funding.
we then utilized results from blastx against model protein sequences to effectively remove trans chimeras. to address these problems, we compare de novo assemblies of paired end bp rna-seq reads using oases, trinity, trans-abyss and soapdenovo-trans to transcripts from genome annotation of the model plant ricinus communis. in order to reduce redundancy, we investigated three methods: 1) using cap and cd-hit-est to combine highly similar transcripts, 2) only retaining the transcript with the highest read coverage, or removing the transcript with the lowest read coverage for each subcomponent in trinity, and 3) filtering oases single k-mer assemblies by number of transcripts per locus and relative transcript length, and then finding the transcript with the highest read coverage. after optimization, seven assembly strategies among all four packages successfully assembled – % of reference genes to more than bp, with a chimera rate of – %, and on average – transcripts per reference gene assembled. rna-seq has shown huge potential for phylogenomic inferences in non-model organisms. the differences were largely due to the number of trans-self chimeras that contain repeats in the opposite direction. more than half of the total chimeras in oases and trinity were trans-self chimeras. with rapidly improving sequencing and assembly tools, our study provides a framework to benchmark and optimize performance before choosing tools or parameter combinations for analyzing short-read rna-seq data. our study demonstrates that choice of assembly package, k-mer sizes, post-assembly redundancy-reduction and chimera cleanup, and strand-specific rna-seq library preparation and assembly dramatically improves gene coverage by non-redundant and non-chimeric transcripts that are optimized for downstream phylogenomic analyses. within each package, we found a trade-off between maximizing reference coverage and minimizing redundancy and chimera rate. however, error, incompleteness, and redundant assembled transcripts for each gene in de novo assembly of short reads cause noise in analyses and a large amount of missing data in the aligned matrix. by doing so we evaluate strategies for optimizing total gene coverage and minimizing assembly chimeras and redundancy.
 rna-seq has shown huge potential for phylogenomic inferences in non-model organisms. however, error, incompleteness, and redundant assembled transcripts for each gene in de novo assembly of short reads cause noise in analyses and a large amount of missing data in the aligned matrix. to address these problems, we compare de novo assemblies of paired end bp rna-seq reads using oases, trinity, trans-abyss and soapdenovo-trans to transcripts from genome annotation of the model plant ricinus communis. by doing so we evaluate strategies for optimizing total gene coverage and minimizing assembly chimeras and redundancy.  we found that the frequency and structure of chimeras vary dramatically among different software packages. the differences were largely due to the number of trans-self chimeras that contain repeats in the opposite direction. more than half of the total chimeras in oases and trinity were trans-self chimeras. within each package, we found a trade-off between maximizing reference coverage and minimizing redundancy and chimera rate. in order to reduce redundancy, we investigated three methods: 1) using cap and cd-hit-est to combine highly similar transcripts, 2) only retaining the transcript with the highest read coverage, or removing the transcript with the lowest read coverage for each subcomponent in trinity, and 3) filtering oases single k-mer assemblies by number of transcripts per locus and relative transcript length, and then finding the transcript with the highest read coverage. we then utilized results from blastx against model protein sequences to effectively remove trans chimeras. after optimization, seven assembly strategies among all four packages successfully assembled – % of reference genes to more than bp, with a chimera rate of – %, and on average – transcripts per reference gene assembled.  with rapidly improving sequencing and assembly tools, our study provides a framework to benchmark and optimize performance before choosing tools or parameter combinations for analyzing short-read rna-seq data. our study demonstrates that choice of assembly package, k-mer sizes, post-assembly redundancy-reduction and chimera cleanup, and strand-specific rna-seq library preparation and assembly dramatically improves gene coverage by non-redundant and non-chimeric transcripts that are optimized for downstream phylogenomic analyses. 1kpchimerade novo assemblyredundancyrna-seqsoapdenovo-transtrans-abysstranscriptometrinityoases 