Augmenting transcriptome assembly by combining de novo and genome-guided tools

Simulating RNA-seq reads

The Arabidopsis thaliana (TAIR10) complete genome sequence, coordinates for genes and transposons were downloaded from ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release and the GFF file was parsed to obtain exonic coordinates. These were used to simulate Illumina-like RNA-seq reads using Flux-simulator (Griebel et al., 2012) (FS-nightly-build_1.1.1-20121119021549) with the following options supplied within a parameter file (Supplemental Methods): NB_MOLECULES (Number of RNA molecules initially in the experiment): 5000000; SIZE_DISTRIBUTION (Gaussian distribution of fragment size as N(Mean, Standard deviation)): N(300,30); READ_NUMBER (Total number of reads to be simulated): 4000000; READ_LENGTH: 76; PAIRED_END: YES; ERR_FILE (Inbuilt Illumina error model for 76nt reads): 76. The resultant Illumina-like reads were split into 2 fastq files corresponding to read1 and read2 using a Python script from the Galaxy tool suite (Supplemental Methods).

Read assembly

The simulated reads were assembled using four de novo (Trinity r2012-06-08; Trans-ABySS v1.3.2; Oases v0.2.08; SOAPdenovo-Trans v1.0) and two genome-guided (TopHat1 v2.0.4 that uses Bowtie1 (Langmead, 2010) v0.12.8 and Cufflinks v2.0.0; genome-guided Trinity r2012-10-05 that uses GSNAP (Wu & Nacu, 2010) r2012-07-20 [V3]) transcriptome assembly pipelines. Trinity and genome-guided Trinity use a fixed k-mer size of 25nt. We used Trans-ABySS on ABySS (Simpson et al., 2009) (v1.3.3) multi-k-mer assemblies (ranging from 20–64nt), Oases on Velvet (Zerbino & Birney, 2008) (v1.2.07) multi-k-mer assemblies (every alternate k-mer ranging from 19–71nt) and SOAPdenovo-Trans with a fixed k-mer size of 23nt. We also tested the TopHat2 (Kim et al., 2013) (v2.0.7 that uses Bowtie v2.0.5)-Cufflinks pipeline on the simulated data but did not observe any difference in the assembly statistics compared to the TopHat1-Cufflinks pipeline that uses Bowtie1. All assemblers were run using default parameters (details in Supplemental Methods). However, we fixed the parameter for minimum length of assembled fragment as 76nt (equal to the length of the read). In the case of SOAPdenovo-Trans, contigs were used instead of scaffolds for all downstream analyses, as the minimum length cutoff could not be set for scaffolds.

Redundancy assessment using CD-HIT-EST

We used CD-HIT-EST (Li & Godzik, 2002) (v4.5.4) to assess redundancy in each assembly. It retained the longest sequence out of a cluster of sequences that share at least 95% sequence similarity (−c 0.95) on either strand (−r 1), based on a word size of 8 (−n 8). The accurate but slow mode (−g 1) was used for clustering. Further details on other optional parameters used to run CD-HIT-EST are provided in Supplemental Methods.

Model assembly

We defined all read-covered transcript regions (TAIR10) as Model Assembly or MA as previously described in the report by Mundry et al. (2012).

Calculation of N50 and N_(MA)50 statistics for an assembly

N50 is defined as the minimum contig length for 50% of the assembly after sorting the contigs in the descending order of their lengths. We calculated N50 values for all assemblies, MA and the TAIR10 simulated transcripts, pre- and post-CD-HIT-EST. We also calculated the N50 values for each assembler, while taking the MA cumulative size (nt) as the denominator instead of the respective assembly sizes. We termed these N50 values as the N_(MA)50.

Mapping assemblies to MA using Megablast

The assembled fragments (query) were mapped against the MA fragments (subject) using Megablast (Altschul et al., 1990) (blast+v2.2.26) with default parameters (Supplemental Methods). The Megablast hits were parsed in order to either maximize query coverage (to compute misassembly statistics) or maximize subject coverage (to compute MA recovery statistics). This was done to choose the best hits and discard the partially overlapping hits when the unique coverage was lower than or equal to 10nt. The scripts used to remove nested and partial overlaps from Megablast hits are provided in Supplemental Methods. For misassembly statistics, each assembled fragment was categorized, based on mapping, to either belong to a single MA source ( ≥ 90% mapping) or be chimeric across multiple sources (misassembled, <90% mapping).

Expression level bins

We estimated the average “per-nucleotide coverage” (pnc) for all MA fragments, based on their read support, as given below. ${\displaystyle \frac{(\text{No. of Reads }\times \text{ Read Length})}{\text{MA fragment length}}.}$ The MA fragments were then categorized into 8 expression level bins, B1 to B8, the pnc for each being: 1 for B1, >1 & ≤ 2 for B2, >2 & ≤ 3 for B3, >3 & ≤ 4 for B4, >4 & ≤ 5 for B5, >5 & ≤ 10 for B6, >10 & ≤ 30 for B7, and >30 for B8. We chose denser sampling for the lower pnc values and sparser sampling for the higher pnc values since we observed the distribution of MA fragments to be denser in the lower pnc categories (Fig. S1).

Recovery of isoforms

Using simulated data, we obtained the MA-equivalent for the exonic regions of isoform-bearing genes, and performed a Megablast search of the assemblies against it (Supplemental Methods). The Megablast hits were parsed to maximize subject coverage, after removing nested and partial overlaps (same as described earlier). For all assemblers, we calculated the number of exons recovered per isoform and the length recovery of each exon.

Augmenting Trinity assembly with TopHat1-Cufflinks assembly

We mapped the Trinity assemblies against the TopHat1-Cufflinks transcripts using Megablast. Using the subtractBed option from BEDTools (Quinlan & Hall, 2010), the regions unique to TopHat1-Cufflinks transcripts were identified. These unique transcript regions were used to augment the Trinity assembly.

Receiver Operating Characteristic (ROC) Curve

The ROC curve was plotted for each assembler using the simulated data. The sensitivity (True Positive Rate, TPR) was estimated as % total length recovered by each assembler out of the total MA size. The FPR (False Positive Rate, 100-specificity) was estimated as the % assembled fragments that did not map to the MA fragments. For TopHat1-Cufflinks assembly, and the unique regions to TopHat1-Cufflinks used to augment Trinity, mapping was performed against TAIR10 transcripts.

Zebrafish transcriptome data analyses

To exemplify the recovery trends of isoforms and non-isoforms, we assembled the transcriptome of zebrafish, Danio rerio, embryo (2dpf) using an RNA-seq dataset from SRA (ERR003998, containing 22,286,504 paired reads, 76nt in length with a 200bp insert size). These reads were assembled using de novo and genome-guided transcriptome assemblers in the same manner as described for the simulated data (Trinity r2012-06-08; Trans-ABySS v1.3.2; Oases v0.2.08; SOAPdenovo-Trans v1.02; TopHat1 v2.0.4 that uses Bowtie1 v0.12.8 and Cufflinks v2.0.0; genome-guided Trinity r2012-10-05 that uses GSNAP r2013-03-31 [V5]). The ENSEMBL zebrafish genome, Zv9, was used as the reference for genome-guided assemblies. Further, 85 transcripts from the zebrafish hox gene cluster were downloaded from ENSEMBL. These correspond to 49 genes, 23 of which harbour 59 isoforms. We focused on the recovery of these transcripts by all assemblers.

Since isoforms carry greater sub-sequence similarity than non-isoforms, we pooled all the 85 transcripts and identified their shared and unique regions using Megablast. RNA-seq reads were then mapped to these, again using Megablast, in order to identify their MA-equivalents and to estimate their average pnc. This served as the reference for recovery comparison across all assemblers. Finally, the recovery of MAs in the shared and unique categories was estimated also using Megablast, as described earlier, while maximising subject coverage.

Assembly statistics

We simulated Illumina-like paired-end 76nt RNA-seq reads for the exonic regions of the Arabidopsis thaliana TAIR10 genome (see Methods for details), covering 15,532 TAIR10 transcripts. The transcript regions contiguously covered by reads were termed as Model Assembly (MA) fragments and were used as a valid reference for mapping the assemblies. A given transcript, therefore, comprised of one or more MA fragments, the shortest being 76nt in length. We obtained 70,382 MA fragments from the 15,532 TAIR10 transcripts.

The reads were assembled using Trinity (TNY), Trans-ABySS (TA), Oases (OS), SOAPdenovo-Trans (SDT), TopHat1-Cufflinks (TC) and genome-guided Trinity (GGT) (see Methods for details) with default parameters and a minimum assembled fragment length of 76nt. We compared the numbers of assembled fragments, the minimum and maximum assembled fragment lengths, their length frequency distributions, the N50 and N_(MA)50 statistics for the six assemblers, pre- and post-redundancy assessment in the assembly using CD-HIT-EST.

The total number and lengths of assembled fragments varied widely across the six assemblers (Table 1). The shortest fragment reported by all assemblers was 76nt (as fixed by the minimum reporting length threshold). The range of fragments at the long end varied across assemblers with the longest for TopHat1-Cufflinks at 10,502 nt (same as the maximum MA fragment length, Table 1). This was expected since TopHat1-Cufflinks is a genome-guided assembler that allows recovery of full-length or near full-length transcripts.

We observed a 4- and ∼1.5-fold reduction in assembly size post-CD-HIT-EST for Trans-ABySS and Oases, respectively (Table 1). For Trinity, genome-guided Trinity and SOAPdenovo-Trans, we observed no difference pre- and post-CD-HIT-EST across the entire distribution of assembled fragment lengths (Fig. 1). However, Trans-ABySS resulted in a higher number of longer, redundant assembled fragments (Fig. 1) and a low number of shorter assembled fragments.

Like the results from the frequency distribution statistics for long- and short-assembled fragments, Trinity, genome-guided Trinity and SOAPdenovo-Trans yielded N50 values closer to that of MA. However, Trans-ABySS, Oases and TopHat1-Cufflinks yielded N50 numbers that were much higher than that of the MA. Post-CD-HIT-EST, the N50 values were not affected for most assemblers except for Trans-ABySS (Table 1).

In contrast, the N_(MA)50 values, which are not dependent on the total size of any assembly, were lower than that of the MA for Trinity, genome-guided Trinity and SOAPdenovo-Trans, both pre- and post-CD-HIT-EST. For Trans-ABySS and Oases, the N_(MA)50 values were relatively higher, and were reduced ∼3 and ∼1.5-fold respectively, post-CD-HIT-EST. The N_(MA)50 values were higher than that of the MA for TopHat1-Cufflinks, both pre- and post-CD-HIT-EST.

Mapping-based statistics

MA fragments, from TAIR10 transcripts corresponding to transposable elements and transcript isoforms as per the GFF annotation, were not included in the expression-based binning. We excluded these elements since they have stretches of identical sequences, which posed a problem in assigning the pnc index due to unreliable mapping of reads and assemblies to the correct isoform and/or transposable element. Out of a total of 70,382 MA fragments, we were thus left with 46,805 MA fragments, which still had certain level of sub-sequence similarity, presumably arising from gene paralogs and SSRs. The numbers of MA fragments finally obtained in different expression bins were: 15,706 in B1 (the lowest expression bin), 9,722 in B2, 5,205 in B3, 2,797 in B4, 1,820 in B5, 4,522 in B6, 4,862 in B7 and 2,171 in B8 (the highest expression bin).

MA recovery statistics

The assembled transcript fragments from all assemblers were mapped to the MA using Megablast to assign common identifiers for comparisons. We compared the numbers and lengths of MA fragments recovered within all expression bins, B1 to B8, using the best hits from Megablast.

Trinity detected the highest number of MA fragments in B1, followed by genome-guided Trinity and SOAPdenovo-Trans (Fig. 2). However, the number of unique transcripts detected was highest for TopHat1-Cufflinks, followed by Trinity, for the lowest expression bin B1 (Fig. 2). The extent of overlapping (common detections) between Trinity and either genome-guided Trinity or SOAPdenovo-Trans was 30% for the lowest expression bin B1. The number of MAs detected was highest for Trinity for all expression bins (B1–B8), with other assemblers showing better recovery with increase in expression. In comparison to Trinity, the detection sensitivity for Oases and Trans-ABySS increased with increased expression, as stated above, but the overall number of recovered MA fragments remained lower in all expression bins (Fig. 2). TopHat1-Cufflinks showed a drop in its detection sensitivity from B1 to B2, and a steady rise thereafter. For all assemblers, with the exception of Trans-ABySS, we observed an increase in higher-order intersections (common detections by 2, 3, 4, 5 or all 6 assemblers) across bins B2–B8, proportional to the decrease in unique detections and lower order intersections (Fig. 2).

We observed that longer MA fragments tend to have higher read support, thus justifying their assignment only to the higher expression bins, B3 onwards (Fig. 3, Fig. S1). The shorter MA fragments displayed all levels of read support, and, therefore, were assigned to all expression bins, B1 to B8 (Fig. 3, Fig. S1). Given this non-uniform distribution of lengths in each expression bin, we compared the length recovery across assemblers, for each expression bin, in different MA fragment length categories (Fig. 3).

For the shortest length category, 76nt, the median length recovery was close to 100% across all expression bins for all assemblers, with the exception of Trans-ABySS (Fig. 3). The outlier (datapoints falling outside the interquartile range of % MA length recovery) trend was highly variable across assemblers, and clustered closer to the median for SOAPdenovo-Trans, followed by Trinity and genome-guided Trinity. For Trinity, Trans-ABySS, Oases and genome-guided Trinity, we observed the outliers clustering around ∼40% length recovery (∼30nt) for the MA fragments in the lowest expression bin (B1) that matched the word size for Megablast (28nt). The outliers in the lowest expression bin for TopHat1-Cufflinks, and those in the highest expression bin for Oases, spanned all the way from ∼40% to the median.

For the subsequent length categories, we observed a gradual increase in the median length recovery from as low as 20% to all the way up to 100%, with an increase in the expression levels. This gradual increase was also seen, to a minor extent, for the same bin across MA length categories for most assemblers. TopHat1-Cufflinks showed a median recovery of 100% across all expression bins and all MA fragment length categories. The pattern of length recovery was similar for Trinity, SOAPdenovo-Trans and genome-guided Trinity across all the expression bins and MA fragment length categories. Overall, these three assemblers outperformed Trans-ABySS and Oases in terms of higher median MA fragment lengh recovery and tighter distribution around the median (Fig. 3).

Misassembly statistics

We used Megablast-based mapping to evaluate the accuracy of assembled fragments, and determine whether each assembled fragment belonged to a single MA source or was a chimera across multiple sources (misassembled). We classified the assembled fragments into three categories, ≥ 90%, between 60%–90%, and <60%, based on the extent of their lengths mapped to any single MA fragment. They were further classified into various assembled fragment length categories to check their relationship with assembly quality, pre- and post-CD-HIT-EST. Trans-ABySS and Oases showed relatively higher numbers of misassembled fragments (between 60%–90%, and <60% mapping) in the 200–400nt and 500–600nt ranges, respectively (Fig. 4). In the case of Trans-ABySS, the extent of misassembly decreased post-CD-HIT-EST, whereas in the case of Oases, it went up. We did not observe any difference in the extent of misassembly, pre- and post-CD-HIT-EST, for any other assembler (Fig. 4). The extent of misassembly was overestimated in TopHat1-Cufflinks as it additionally assembled transcript regions not represented by reads, and therefore the MA, using the genomic information (Table 1). For Trinity, SOAPdenovo-Trans and genome-guided Trinity, there was a clear trend of decreasing misassembly with increasing assembled fragment length (Fig. 4). The degree of misassembly was less for SOAPdenovo-Trans than for Trinity or genome-guided Trinity for the shorter fragments.

Recovery of isoforms using simulated and zebrafish data

Next, we compared the extent of recovery of isoforms from the transcripts assembled with various assemblers. We had a total of 24,846 exons in the MA fragments corresponding to 2,970 isoforms. We mapped assembled fragments from all assemblers to these exons, using Megablast, and maximized the exon coverage. Trinity, followed by genome-guided Trinity and SOAPdenovo-Trans recovered the highest numbers of isoforms to 80%–100% of their length (Fig. 5). We observed a correlation between the median pnc for each length recovery category and the numbers of exons recovered per isoform. The median pnc was in the range of 2–3 for unrecovered isoforms and 17–23 for those which were fully recovered or close to fully recovered by all assemblers. Trans-ABySS and Oases, which showed a relatively lower recovery, correlated with a higher median pnc, suggesting a higher threshold of pnc needed for good recovery by those assemblers (Fig. 5).

In order to test whether our understanding of isoform recovery from simulated reads holds true for a real dataset, we measured the recovery of shared and unique transcript regions of the hox gene cluster in zebrafish (see Methods for details). We observed that across assemblers, Trinity, followed by genome-guided Trinity and SOAPdenovo-Trans recovered most and closer to full-length transcripts (Figs. 6–8). TopHat1-Cufflinks recovered mostly full-length, but fewer, transcripts. Trans-ABySS and Oases recovered fewer and truncated transcripts. The length recovery (% of MA length) by all assemblers, except TopHat1-Cufflinks, displayed a dependency on pnc which was more obvious in the recovery of shared regions of transcripts, as expected due to a wider range of read depth (Figs. 6–7).

Augmenting transcriptome assembly

Since each assembler produced a set of unique transcripts or fragments, we proceeded towards augmenting the assemblies, one with another. We ruled out most combinations of assemblers, and chose Trinity and TopHat1-Cufflinks, as they produced the maximum number of valid transcript fragments and full-length transcripts, respectively. After aligning the Trinity transcript fragments to the TopHat1-Cufflinks transcripts and augmenting the unique regions from TopHat1-Cufflinks to the Trinity assembly, we obtained a cumulative size increase of 1.37 Mb in the assembly size (∼20% of original Trinity assembly). The augmented assembly detected 1,377 MAs more than the original Trinity assembly, and the MA recovery increased by 5,23,127nt after augmentation. Relaxing the stringency of Megablast word size, from 28 to lower, would have improved the MA length recovery further, though at the expense of losing isoform reconstruction capability and sensitivity to variation in read-depth maintained by the fragmented structure of Trinity transcriptome assembly. We further observed that the median length recovery for the augmented Trinity (AGT) was overall better than Trinity alone (Fig. 9). For MAs longer than 1500nt, this augmentation yielded full-length transcripts across all expression levels (Fig. 9). Even for MAs around 500nt in length, we saw the advantage of augmentation, at least for recovering more number of full-length transcripts. The outliers in augmented Trinity were fewer than TopHat1-Cufflinks in all length categories. The improved recovery with augmented Trinity proves that an integrative approach is useful, particularly when one has access to genome in addition to the RNA-seq reads.