RNA sequencing read depth requirement for optimal transcriptome coverage in Hevea brasiliensis

Background One of the concerns of assembling de novo transcriptomes is determining the amount of read sequences required to ensure a comprehensive coverage of genes expressed in a particular sample. In this report, we describe the use of Illumina paired-end RNA-Seq (PE RNA-Seq) reads from Hevea brasiliensis (rubber tree) bark to devise a transcript mapping approach for the estimation of the read amount needed for deep transcriptome coverage. Findings We optimized the assembly of a Hevea bark transcriptome based on 16 Gb Illumina PE RNA-Seq reads using the Oases assembler across a range of k-mer sizes. We then assessed assembly quality based on transcript N50 length and transcript mapping statistics in relation to (a) known Hevea cDNAs with complete open reading frames, (b) a set of core eukaryotic genes and (c) Hevea genome scaffolds. This was followed by a systematic transcript mapping process where sub-assemblies from a series of incremental amounts of bark transcripts were aligned to transcripts from the entire bark transcriptome assembly. The exercise served to relate read amounts to the degree of transcript mapping level, the latter being an indicator of the coverage of gene transcripts expressed in the sample. As read amounts or datasize increased toward 16 Gb, the number of transcripts mapped to the entire bark assembly approached saturation. A colour matrix was subsequently generated to illustrate sequencing depth requirement in relation to the degree of coverage of total sample transcripts. Conclusions We devised a procedure, the “transcript mapping saturation test”, to estimate the amount of RNA-Seq reads needed for deep coverage of transcriptomes. For Hevea de novo assembly, we propose generating between 5–8 Gb reads, whereby around 90% transcript coverage could be achieved with optimized k-mers and transcript N50 length. The principle behind this methodology may also be applied to other non-model plants, or with reads from other second generation sequencing platforms.


Background
Transcriptome analysis has become increasingly powerful through advances in second generation sequencing technologies from companies such as Illumina, Roche and Life Technologies. Improvements in sequencing chemistry and read length have enabled unprecedented depth of sequencing, limited only by cost and availability of biological material. In particular, RNA sequencing (RNA-Seq), also known as whole transcriptome shotgun sequencing, has emerged as a valuable tool for profiling expressed genes in plants and other organisms [1][2][3]. The depth of transcriptome sequencing provided by RNA-Seq has thus provided a cost-effective means of qualitative and quantitative analyses of gene transcripts in many non-model plant species including the rubber tree, Hevea brasiliensis, the subject of the present study.
Challenges in de novo transcriptome assembly in higher plants lie in the immense number of gene transcripts, large variations in transcript expression levels, presence of alternatively spliced transcript variants and issues in strand directionality [11,12]. Owing to such problems, de novo assembly requires significantly greater sequencing depth as compared with reference-based assembly. Especially in the case of non-model plants, few guidelines are available for determining the amount of reads to generate to enable deep coverage of transcripts expressed in a particular sample. The general practice commonly adopted, especially for new entrants in second generation sequencing, is to piggy-back on ballpark estimates adopted for the model species. From a survey of recent publications on de novo transcriptome analysis in non-model plants, the read generation per sample could fall below 100 Mb or it could be high as 7 Gb, with 2-5 Gb being the most common sequencing depths [13,18,19,21,23,. Therefore, there is need for a practical procedure to estimate the reads needed for deep coverage of gene transcripts in de novo assembly where such information is unavailable.
Hevea brasiliensis is the commercial source of natural rubber (cis-polyisoprene). The tree, a diploid species (2n = 36) from the Euphorbiaceae family, has a C value of about 2 pg as estimated by flow cytometery [50]. Compared to model plants, transcript resources in this economically important species have only been initiated relatively recently, with the analysis of Expressed Sequence Tags (ESTs) from latex [29,[51][52][53][54]. The first crop of publications on the application of second generation sequencing in Hevea rubber genomics emerged in 2011 (Table 1). In these papers, functional profiling of gene transcripts was reported in latex, leaf and bark tissues, focussing on rubber biosynthesis, molecular markers and transcription factors [29,[55][56][57][58][59][60]. Overall, the depth of sequencing reported in these papers was between 0.2-5 Gb per sample based on Illumina PE-RNA-Seq or 454 pyrosequencing platforms (Table 1). Hevea being a tree, analysis of its gene expression is often in RNAs prepared from distinct cells, tissues or organs, including RNAs from the same sample types but under different physiological conditions [61][62][63][64]. Having a means of assessing the degree of coverage of genes that are expressed in various samples would be important to achieve a normalized set of Hevea gene transcripts. In this paper, we describe a method to estimate datasize requirement for high transcriptome coverage based on an analysis of assembly statistics of Illumina PE RNA-Seq reads generated from Hevea bark. We first optimized and validated a 16 Gb bark assembly across a range of datasizes and k-mers. Subsequently, we applied a transcript mappability method to illustrate the trend of gene transcript coverage by different read amounts that had been assembled using a range of k-mers.

Findings
Generation of Illumina PE RNA-Seq Hevea tissue libraries and de novo assembly The development of an approach to determine the datasize required for deep coverage of a de novo transcriptome arose from the generation of three Hevea tissue read libraries as part of a programme in developing genomic resources for the rubber tree. Considerably more PE RNA-Seq raw reads were available for bark (16.9 Gb) as compared with latex and leaf tissues (4.9-5 Gb each) ( Table 2; see Materials and methods). Hitherto, similar work on the rubber tree involving second generation sequencing (Table 1) has far lower transcriptome coverage per sample compared to what we have generated. Taking cognizance of the immensity of the bark read library, we asked the following questions: (a) What is the trend of assembly characteristics of a significantly larger read library? and (b) Can we devise a method to relate datasize requirement to the degree of coverage of transcripts expressed in the sample?
We first processed the bark, latex and leaf raw reads for quality, and found a very high percentage of clean reads (approximately 96-98%) ( Table 2; see Materials and methods). The distribution of read length categories also indicated high PE RNA-Seq quality based on the fact that a vast majority of clean paired reads was in the largest size category (100 nt) ( Figure 1). Henceforth in this report, the clean paired reads from bark are referred to as the 16 Gb read set while those from latex or leaf as the 5 Gb read set ( Table 2). A number of plant de novo transcriptome projects have reported the use of multiple k-mers in assembly optimization [13,21,27,65]. Therefore, we analyzed the effect of two parameters on transcriptome assembly, namely read amount and k-mer size.
To do this, the 16 Gb bark read set was assembled in incremental quantums of 1, 3, 5, 8, 10, 13 and 16 Gb using the Oases assembler and a k-mer range of 51-77 (see Materials and methods). As shown in Table 3, the higher the k-mer, the lower was the number of transcripts generated by a particular datasize. At the same time, the larger the dataset size, the greater the number of transcripts generated using a particular k-mer. On the other hand, the transcript N50 length (see Materials and methods) showed a bell shape distribution from k-mer 51 to 77 in all datasize sub-assemblies with the exception of the 1 Gb datasize (N50 peak values are highlighted in Table 3). The k-mer at which the N50 value peaked in this bell shape distribution was thus referred to as the "optimized k-mer" of assembly for a particular datasize (and accordingly the "optimized N50"). As seen in Table 3, the optimized N50 became larger with increasing read depth, and this corresponded also with increment in the size of optimized k-mers. Hence, larger data sizes facilitated not only an increase in the number of transcripts assembled, but also an increase in transcript length. We therefore propose that by a judicious combination of datasize and kmer range, the ensuing N50 trend may serve as a criterion for determining optimal de novo assembly.

Validation of the bark transcriptome
In this study, a k-mer of 73 assembled 87,612 transcripts from the 16 Gb read set with a transcript N50 value of 2,068 bp (Table 3). At this stage, this assembly was referred to as the optimized bark assembly for 16 Gb reads. Due to the importance of quality de novo assembly, we performed three mapping analyses to validate the 16 Gb bark transcriptome (see Materials and methods). First of all, 255 publicly available Hevea cDNAs which had been verified to contain complete open reading frames or ORFs were mapped to 87,612 bark transcripts using the Megablast software. In the absence of more cDNAs containing complete Hevea proteins, rubber-specific ORF quality of bark transcripts could only be evaluated using these sequences which also included isoforms for several gene families (see Additional file 1: Table S1). The results showed that 250 Hevea ORFs (of 255) had hits to bark transcripts with sequence identity match ranging from 86-100% (Table 4 and Additional file 1: Table S1). Also, 80% of the 250 ORFs showed a minimum of 70% utilization of ORF sequence length in their match alignments with bark transcripts (Table 4 and Additional file 1: Table S1). These observations indicated that a high proportion of transcripts within the optimized bark assembly encoded complete or near-complete Hevea proteins. Secondly, completeness of gene representation in the 16 Gb transcriptome was assessed by mapping 87,612 bark transcripts to a set of 248 core eukaryotic genes (CEGs) that had been shown to be a reliable indicator of completeness of gene space in eukaryotic species [66]. Although initially used to assess gene space in newly sequenced genomes, the approach was recently applied to transcriptomes and it also complements other transcript metrics such as N50 length. Using BlastX, 87,612 bark transcripts detected 247 out of 248 CEG proteins (98%) (e-value ≤ 1e -10 ). Thus, these results support the completeness or depth of bark gene representation in this transcriptome.
Thirdly, the 87,612 bark transcripts were validated by using the Exonerate software to map them to rubber genome scaffolds (BioProject ID: PRJNA80191). Figure 2 shows the number of bark transcripts which were mapped to categories of transcript-to-scaffold coverage. The higher the percentage of transcript-to-scaffold coverage, the greater the significance of match alignment. A large proportion of bark transcripts (84,471 or 96.41% of total) could be mapped to the scaffolds and of these, almost 70% showed transcript-to-scaffold coverage of 90-100% ( Figure 2). Therefore, this indicated the presence of a sizeable proportion of high quality bark transcripts.
As a whole, results of the three mapping analyses carried out provided sufficient validation for the quality of 87,612 bark transcripts assembled from the 16 Gb read set. Fragmented or erroneous transcripts could still be present to some extent in any assembly but we think that the proportion of bark transcripts which did not show meaningful mapping or alignment with Hevea transcripts or genome scaffolds could also be explained by reasons such as inherent variations between sequences derived from different tree clonal varieties.

Mapping saturation test for bark transcript accumulation
Next, we addressed the question of read depth requirement by mapping the series of incremental bark transcripts to total transcripts from the optimized 16 Gb bark assembly. The principle behind mapping subsets of bark transcripts to total transcripts from the full 16 Gb bark transcriptome is that the number of the former aligning to the latter should follow a saturation pattern. As outlined in Figure 3, transcript sets from 1, 3, 5, 8, 10 and 13 Gb bark reads that had been assembled independently across k-mers 51-77 were mapped to 87,612 bark transcripts. The extent of transcript mapping, expressed as a percentage of total bark transcripts (or the full transcriptome) was taken as a measure of transcript representation by each sub-assembly.
Results of the number of mappings to 87,612 bark transcripts were displayed as a colour matrix as shown in Figure 4 and Additional file 2: Table S2. Overall, transcript mapping by incremental subsets of bark transcripts approached full saturation (i.e. 100%) as datasize increased (for any k-mer) or as k-mer decreased (for any datasize). Transcript representation was high in all cases, being more than 80% transcript mapping level with the exception of three higher k-mer assemblies (73-77) of 3 Gb reads (Figure 4 and Additional file 2: Table S2). At the lower end of the saturation spectrum, 3-8 Gb reads  The peak N50 length from each datasize assembly (other than for 1 Gb) is highlighted in bold italics.
would already be sufficient to obtain even 80-85% mapping level or transcript coverage. This datasize range could also readily generate transcript coverage exceeding 90% if the k-mer range for assembly was decreased. At the upper end of the saturation spectrum, transcripts from nearly all of the 10-13 Gb assemblies across all kmers showed mapping levels greater than 90%. However, although the colour matrix indicated generally high transcript coverage by assemblies of 3-13 Gb reads, it is important to select a datasize that would also produce the optimal transcript N50 length at the desired transcript coverage level. As determined previously, the optimized N50 increased with datasize (  Table S2). Therefore, based on the colour matrix, the optimized 3-5 Gb assemblies would fall within the 85-90% transcript coverage bracket and the optimized 8-13 Gb assemblies within the 90-95% coverage bracket (Figure 4 and Additional file 2: Table S2). This also indicated that based on the mapping saturation test in this study, a shift in bark transcript coverage bracket by the incremental assemblies occurred between 5-8 Gb.
In essence, the amount of reads for optimal coverage of a tissue transcriptome should take into consideration the requirements for transcript coverage level (reflected by percentage of mapping saturation) and for N50 length. In order to attain both high transcript representation and best N50 length, our general recommendation for Hevea is to generate between 5-8 Gb reads for de novo assembly. Firstly, as observed in the shift in transcript coverage bracket, a minimum of nearly 90% (i.e. 89.48%) transcript coverage could already be achieved by an optimized 5 Gb assembly (Figure 4 and Additional file 2: Table S2), a level that is within range of the following coverage bracket (90-95%). Secondly, it is noted from Table 3 that the improvement in optimized N50 length was more rapid from 3-8 Gb than from 8-16 Gb assemblies. Therefore, generating less than 5 Gb reads may lead to reduction in complete transcripts, and sequencing beyond 8 Gb reads may not yield significantly more new or complete transcripts other than the rarely expressed ones.

Analysis of latex and leaf transcriptome assembly
RNAs from different types of tissues and growth conditions are of interest in Hevea transcriptome profiling studies. Application of the recommended 5-8 Gb sequencing depth for Hevea would assume the assembly trends to be the same between reads generated from bark and from other tissues. To validate this, we performed Oases assembly of 5 Gb latex and 5 Gb leaf read sets using the same parameters as those for bark assembly (see Materials and methods). Assemblies were performed in incremental read amounts of 1, 3 and 5 Gb across a k-mer range of 51-77. Table 5 shows the statistics of latex and leaf assemblies which correspond to the same quantums of read assemblies in Table 3. As with bark assembly, the transcript N50 length showed a peaking trend only in the 3 and 5 Gb leaf read assemblies across the k-mer range but not in the 1 Gb assembly. On the other hand, the transcript N50 length showed a peak for all three latex assemblies across the k-mers. For the 5 Gb datasize assemblies, the optimized N50 was highest for bark (1,900 bp) followed by latex (1,281 bp) and leaf (1,086 bp) (Tables 3 and 5). Differences in statistics such as total assembled transcripts and optimized N50 values could be due to the presence of tissue-specific genes and variation in dynamic range of transcript levels in tissues, all of which have effect on the outcome of assembly. However, similar to that for bark, the optimized N50 of latex and leaf assemblies increased with read amounts. In conclusion, assembly of reads generated from different tissues did not display major or unexplainable differences between one another, and therefore, the mapping saturation test should be applicable to other Hevea tissues.

Discussion and conclusions
Knowing whether transcriptome sequencing and assembly have substantially captured all the genes expressed in a sample is an important consideration for plants having limited genomic resources as reference. Generally, the amount of reads for comprehensive coverage of a de novo transcriptome is often determined by a balance of budget, capacity of sequencing platform and guesstimates or "best practices" based on other species. In this work, we report a systematic approach which we name the "transcript mapping saturation test" to assess the amount of reads required for optimal transcriptome coverage in the Hevea rubber tree. This was made possible by the availability of 16 Gb Illumina PE RNA-Seq reads from Hevea bark which enabled us to map transcripts from incremental sub-assemblies to transcripts from the entire assembly (or the full transcriptome) in order to detect the mapping saturation point. The workflow of this methodology is outlined in Figure 3, beginning with assembly optimization and validation of the full transcriptome, followed by the mapping saturation test. Because sequencing has become increasingly affordable, obtaining as much as 16 Gb reads per sample as a starting point is not insurmountable. Using our approach, sequencing to this extent has to be done only once in the beginning, after which the user is equipped with a guide (the colour matrix) to estimate optimal coverage of expressed genes in the plant species of interest. In developing this approach, we used the Oases assembler because Velvet, for which Oases is an extension, had previously been found to be suitable for producing quality transcripts from Hevea short reads [29]. Thus, we progressed to Oases, which additionally has the ability to resolve alternatively spliced transcripts [5]. We would suggest that a de novo project intending to adopt our approach should first test if the assembler of choice is suited to their transcriptome. Even though our method development is based on Illumina PE RNA-Seq reads, the principle behind this approach should also be applicable to other plant species and to reads from other sequencing platforms.
Although the transcript mapping approach is based on mapping saturation, this does not reduce the need to validate the assembly quality of the full transcriptome. In this work, the N50 trend was used in the initial selection of best k-mer for assembly. Subsequently, the completeness and correctness of assembled transcripts were supported by results of mapping to rubber genome scaffolds [60] whereby a significant proportion showed transcriptto-scaffold coverage of 90% and above. This was also supported by detection of all but one of 248 core genes expressed in eukaryotes [66] and significant alignments with known Hevea protein coding frames. However, we should point out that what this paper proposes is essentially a methodology; we do not specifically assert that a 5-8 Gb read depth would be sufficient for optimal transcriptome coverage universally. The optimal read depth  Between 90-95% Between 95-99% 99% and above Figure 4 Colour matrix representing the transcript mapping saturation test results. An asterisk corresponds to the assembly with the optimized k-mer and transcript N50 length for a particular datasize (1 Gb assembly transcripts were not included since the transcript N50 length was not optimized by any k-mer for this datasize). See Additional file 2: Table S2 for full details of BlastN matches by subsets of bark transcripts to the optimized 16 Gb bark transcriptome.

Plant material and RNA isolation
Latex and bark shavings were obtained from 15-year old RRIM 928 Hevea brasiliensis trees growing in the Rubber Research Institute of Malaysia Research Station, Sungai Buloh. Equal volumes of latex were tapped from three trees and collected directly into 2× RNA extraction buffer [67]. The bark just below the tapping cut of the trees was scraped to remove surface matter before bark shavings (approximately 1 cm depth) were excised with a tapping knife. Young leaves of RRIM 928 trees were collected from the source bush nursery in the Rubber Research Institute of Malaysia Research Station, Sungai Buloh. Total RNA was isolated from latex and leaf tissues using the phenol-chloroform method [67]. Bark total RNA was isolated using a modified procedure of the Qiagen RNeasy Plant Mini Kit [68]. RNA samples were assessed for quality and quantity using the Nanodrop spectrophotometer (Thermo Scientific).

Sequence generation and quality assessment
Bark, latex and leaf total RNAs (20 μg each) were sent to the Illumina Fast Track sequencing service in San Diego, USA where 200 bp fragment size libraries were produced for paired-end RNA sequencing (PE RNA-Seq). Each RNA-Seq sample was sequenced 100 nucleotides at each end (2 × 100 nt), resulting in about 50 million raw reads each from latex and leaf and nearly 170 million raw reads from the bark ( Table 2). Raw reads from bark, latex and leaf are available from the NCBI Sequence Read Archive (accession nos. SRX278513-5).
Clean reads were obtained by trimming raw reads at a minimum phred score of Q = 20, followed by removal of reads below 30 bp and subsequently reads which contained 'N' nucleotides. Clean paired reads from the bark (163,316,702 reads; see Table 2) were referred to as the 16 Gb read set and clean paired read sets from the latex and leaves (48,650,932 and 46,062,766 reads respectively; see Table 2) as the 5 Gb read sets. These read sets were used for subsequent transcriptome assembly. Clean paired reads were classified into arbitrary nucleotide size categories to confirm good PE RNA-Seq data quality (Figure 1).

Transcriptome assembly and transcript mapping
Clean paired reads from the bark, latex and leaf (Table 2) were assembled with the Velvet (Version 1.1.05) [4] and Oases assembler (Version 0.1.22) [5] using default parameters and selection of a minimum transcript length of 100 bp. A range of hash lengths (k-mers 51-77) was used for assembly of a read set to determine the k-mer which produced the highest transcript N50 length. This best N50 value was termed as the "optimized N50" while the hash length which produced it was the "optimized kmer". Note: N50 length is the length of the shortest transcript whereby the sum of transcripts of equal length or longer is at least 50% of the total length of all transcripts.
Incremental quantums of bark reads (1, 3, 5, 8, 10, 13 and 16 Gb) were obtained by partitioning the subsets from the 16 Gb read set. Each subset was random as the 16 Gb read set was already fully randomized. (Similarly, 5 Gb read sets from latex and leaf were also fully randomized). Serial mapping of bark transcripts was performed using BlastN [69] and the top hit by any query transcript with e-value ≤ 1.0e -5 was counted as a match. The complete methodology for the transcript mapping saturation test is shown in Figure 3.

Bark transcript validation
For evaluation of rubber-specific ORF quality of 87,612 bark transcripts from the 16 Gb assembly (k-mer 73), 255 Hevea sequences which were confirmed to encode complete ORFs were selected from the NCBI GenBank non-redundant database (www.ncbi.nlm.nih.gov). These Hevea cDNAs, which were isolated by traditional gene cloning approaches such as cDNA library hybridization and PCR, were generally of high quality as they had mainly been obtained by Sanger sequencing (see Additional file 1: Table S1). Megablast [69] was used to map the 255 Hevea ORFs to 87,612 transcripts from the 16 Gb bark assembly. Top hits from this analysis (with 86-100% sequence identity match) were screened for high quality matches based on a minimum of 70% coverage of Hevea ORFs (or query coverage) in their alignments to bark transcripts (the subject) (see Table 4 and Additional file 1: Table S1).
For evaluation of completeness of assembled bark transcripts, 248 core eukaryotic genes (CEGs) of Arabidopsis thaliana were downloaded from the CEGMA resource at korflab.ucdavis.edu/Datasets/genome_completeness/. This approach was based on a list of 248 highly conserved but low copy number genes that had been shown to be a reliable indicator of completeness of gene space in eukaryotic species [66]. Using BlastX [69], 87,612 bark transcripts were mapped to the CEGs with any hit of e-value ≤ 1.0e -10 counted as a match.
Rubber genome scaffolds from the BioProject ID: PRJNA80191 (www.ncbi.nlm.nih.gov/nuccore/448814761) were used for validating bark transcripts. Bark transcripts were mapped to genome scaffolds by Exonerate (Version 2.2.0) [70] using default settings with the exception of the following parameters: heuristic mode, est2genome model and alignment score of at least 10 percent of the maximal score for each query. The significance of mapped transcripts was evaluated by calculating query coverage which is expressed as percentage of the transcript sequence (query) that overlaps with the scaffold sequence (subject).