Exploiting single-molecule transcript sequencing for eukaryotic gene prediction

Genome assemblies

Genes were predicted based on sugar beet assembly RefBeet-1.2 [13] and spinach assembly Spinach-1.0.1 (assemblies accessible at [21]). Adapted assembly versions compliant with GenBank submission specifications can be accessed at GenBank with accession numbers AYZS02000000 (RefBeet), and AYZV02000000 (spinach).

Sample preparation and SMRT sequencing

Sugar beet seedlings were obtained by incubation of seeds of sugar beet DH line KWS2320 at 20 °C for 48 h. Seedlings were removed from the dish once cotyledons had become fully visible. Seeds that had not germinated were discarded; the success of germination was about 75 %. Total RNA was isolated from the seeds using the Nucleospin Plant RNA kit (Macherey-Nagel, Düren, Germany). Ten nanograms of total RNA was reverse transcribed using the SMARTer PCR cDNA Synthesis Kit (Clontech, Mountain View, CA, USA), and cDNA was amplified using the Advantage 2 PCR kit (Clontech) for 12 cycles. The generated cDNA was then re-amplified in nine independent PCR reactions using the Advantage 2 PCR kit and the IS primer (0.4 μM final concentration) for 30 cycles. A total of 1 μl of generated cDNA was used in each reaction. Re-amplified cDNA was purified using the QIAquick PCR Purification kit (Qiagen, Hilden, Germany) and thereafter size-selected on agarose gels into cDNA fractions of length 1-2 kb, 2-3 kb, and >3 kb. Excised fractions were column-purified using the QIAquick Gel Extraction kit (Qiagen). The fragment size distribution was validated on a Bioanalyzer HS chip (Agilent, Santa Clara, CA, USA) and quantified on a Qubit fluorometer (Life Technologies, Carlsbad, CA, USA). The three cDNA size fractions were submitted to KeyGene N.V. (Wageningen, the Netherlands), where library preparation and SMRT sequencing on the Pacific Biosciences RS sequencing instrument was carried out. MacBead loading and SMRT C2 sequencing chemistry were used together with XL polymerase. Two SMRT cells were run from each of the three fractions, that is, one movie each of 1 × 120 min and one movie each of 2 × 55 min.

Data analysis overview

The analysis steps are summarized in Fig. 6.

Pre-processing of SMRT reads

Raw SMRT sequencing reads consist of one or several subreads representing the same circularized cDNA template separated by adapter sequence. Subreads smaller than 50 bp and reads with a quality less than 0.75 (corresponding to a predicted error rate of >25%) were removed. Remaining subreads were merged into consensus sequence (CCS) reads, whenever the entire cDNA template was covered by at least two subreads, using the Pacific Biosciences SMRT-Analysis pipeline v2.0. Full-insert sequencing reads were identified from single-pass subreads and CCS reads by running hmmer_wrapper.py [22]. Reads with both 5' and 3' primer sequences and a poly(A) tail present were considered to represent full-length transcripts. Primer sequences as well as the poly(A) tail were trimmed off prior to further analysis.

Identification of SMRT reads containing full-length ORFs

The percentage of SMRT reads containing complete ORFs was estimated from SMRT reads overlapping full-ORF BeetSet-2 genes. SMRT CCS and unmerged subreads were aligned to the reference sequence RefBeet-1.2 with GMAP (v2012-07-20, -A -f gff3_match_cdna -B 5). The CDS of 9,422 BeetSet-2 genes was covered (> = 1 base overlap) by one or more SMRT reads. For the BeetSet genes only the transcript with the longest CDS was considered. The completeness of BeetSet-2 genes was inferred from multiple protein alignments to four different eudicot plant species. Protein sequences were downloaded from plants.ensembl.org: A. thaliana (TAIR10.24), V. vinifera (IGGP_12x.24), S. tuberosum (3.0.24), P. trichocarpa (JGI2.0.24). The corresponding protein sequences of the 9,422 BeetSet-2 genes were aligned to the protein sequence of these four other plant species using BLAT (v34, proteins, default). Alignments with less than 50 % matching bases of BeetSet-2 proteins were discarded. Each BeetSet-2 gene and the most similar gene from the other species (at most one per species) were realigned in a multiple alignment step using ClustalW (v1.83). BeetSet-2 genes were regarded as complete if they had either the same start and end coordinates or additional amino acids compared to at least one other species. A total of 7,286 BeetSet-2 genes passed these criteria. 244,542 PacBio reads partially or fully overlapped with these BeetSet genes. The percentage of SMRT reads with full length ORF was calculated from this subset. A SMRT sequence was considered full ORF, when the entire ORF of the underlying sugar beet gene was included in the SMRT sequence alignment. UTR was considered as present if at least 10 additional bases up- and downstream of the BeetSet-2 ORF were present.

Correction of SMRT sequences

Full-insert transcripts determined from CSS and unmerged subreads, remaining CCS, and remaining unmerged subreads were corrected separately with proovread [19] (default settings) using quality-filtered and trimmed Illumina reads [13, 23]. The Illumina data used for correction had been generated from inflorescence tissue of the sugar beet reference genotype KWS2320. In order to accelerate SMRT data correction, Illumina reads were normalized to a maximum coverage of 100-fold per transcript using normalize-by-median.py (parameters: -N 4 -x 48e9 -k 17 -C 100) of the khmer package [24]. Normalization excluded 70 % of the Illumina reads, leaving 21,132,501 reads for the correction. By default, proovread performs quality trimming after error correction, but also provides the error-corrected untrimmed reads. To maintain the full length of the processed SMRT sequences only the corrected but untrimmed SMRT sequences were used for further analysis.

To determine the sequence accuracy of the sugar beet SMRT data before and after error correction, the sequences were aligned against the set of protein coding gene transcripts (v1302) predicted in the sugar beet assembly RefBeet-1.1 using blasr [25, 26] (part of SMRT analysis pipeline v2.0, assembly and gene set available on [21]). For each SMRT sequence only the best alignment was retained and only if at least 90 % of its length matched. The accuracy was determined based on the average alignment identity. RefBeet-1.1 gene predictions have high consensus accuracy [13]. Uncertainties in the accuracy estimation due to potential structural errors in the previously predicted transcripts were avoided by only considering full-length aligning SMRT reads (>= 90% of length).

Stress conditions and mRNA-seq data

KWS2320 sugar beet plants for abiotic stress treatment (salt, heat, light, and control) were grown in hydroculture in Hoagland’s medium with weekly exchange under a 10 h light/14 h dark cycle at constant temperature of 21 °C, 60 % relative humidity, and light intensity of 80–120 μmol m⁻² s⁻¹. For the salt treatment, the plants were transferred into fresh nutrient solution with 50 mM NaCl at day 23 after sowing. The salt concentration was stepwise increased to 300 mM at day 28 and then kept until harvesting. For the light treatment plants were exposed to a light intensity of 800 μmol m⁻² s⁻¹ for 4 h prior to harvest. The heat treatment was performed by exposing the plants to a temperature of 35 °C for 3 h prior to harvest. All plants were harvested 30 days after germination, and the material was immediately frozen in liquid nitrogen. Total RNA was isolated from leaves by phenol:chloroform:isoamylalcohol (25:24:1) extraction and LiCl precipitation. The resulting RNA was treated with DNA-free DNAse (Life Technologies), and subsequently quantified. Success of the stress treatment was validated by detection of differential gene expression for genes expected to be stress-responsive. Among the abiotic stress-induced genes were MYB12, CHS (encoding chalcone synthase), ELIP genes encoding early light-inducible proteins, as well as genes coding for ethylene-responsive transcription factor, heat shock proteins and glutathione S-transferases. Among the 120 genes which showed the highest differential expression (>100-fold) between stressed conditions vs. controls were five genes previously predicted in RefBeet-1.1 with an expression evidence of <1 %.

Non-directional cDNA libraries were sequenced on a HiSeq1500 instrument (Illumina, San Diego, CA, USA) according to the instructions of the manufacturer. Spinach mRNA-seq was performed on an inbred Viroflay variety (Syngenta Seeds, the Netherlands), isogenic to the published reference assembly Spinach-1.0.1 [13]. RNA was isolated from young leaves and apical shoots using a Nucleospin RNA Plant kit (Macherey-Nagel, Düren, Germany). Non-directional mRNA-seq libraries were prepared using Illumina kit RS-100-0801, and were sequenced on the Illumina HiSeq2000 using a 2 × 50 cycle paired-end sequencing protocol.

Identification of accurate sugar beet gene models for parameter training

An initial gene set was predicted in the sugar beet reference assembly RefBeet-1.2 with AUGUSTUS v2.7 using the same settings and expression data (mRNA-seq data: SRA accessions SRX287608-SRX287625, Roche/454 data: SRA accession SRX287606, 35,523 ESTs) as in the RefBeet-1.1 annotation v1302 [13]. Error-corrected SMRT sequences were mapped to RefBeet-1.2 using Gmap [27], parameters -A -f gff3_match_cdna -B 5. SMRT sequences aligning to multiple locations or not mapping at their full length were removed. To account for deletion errors, SMRT sequences were considered to map at full length if they aligned with at least 90 % of their lengths, missing at most 50 bases. Gene models of the initial gene set in gff format were validated by SMRT sequences using the custom Perl (v5.10.1) script ‘validate-gene-models-with-PacBio.pl’ (Additional files 2, 3, and 4). A gene model was considered validated if all SMRT sequences aligning to the exons of the gene model confirmed the same number and order of exons as well as exactly the same intron boundaries. The UTR length was allowed to be variable, predicted start and stop codons were required to be covered. Thirty-two genes at scaffold borders (within the first or last 5,000 bases) were removed, since they may represent incomplete genes. In total, 2,267 SMRT gene models were validated by this method.

Due to sequencing errors in SMRT sequences, intron boundaries of the initial gene models did not always coincide with those of the aligned SMRT sequences. For such cases and for cases in which the initial gene model had been erroneously predicted, additional gene models were derived directly from aligned SMRT sequences. Gene models were clustered based on their alignment location and intron boundaries using the custom Perl script ‘derive-gene-models-from-PacBio.pl’ (Additional files 2, 3, and 4). The most abundant isoform per location was selected. Transcript boundaries were derived from the median start and stop positions of all aligned SMRT-sequences representing a selected isoform. Open reading frames (ORFs) were predicted with TransDecoder (Brian Haas et al., unpublished [28]). In this way, 665 additional gene models, non-overlapping with previous SMRT-validated genes, were obtained. The SMRT-based validation was developed as a fast and entirely automatic process without the need of manual curation.

For the manually validated genes, arbitrarily chosen RefBeet-1.1 gene predictions with a hint coverage of at least 85 % were extracted using a customized Perl script ‘parse_AUGUSTUS_gff3.pl’ (Additional files 2, 3, and 4). The exon number of candidate genes was defined based on the intron-exon distribution of the RefBeet-1.1 gene set. Genes were manually inspected using GenomeView [29], visualizing RefBeet-1.1 gene predictions and supporting mRNA-seq evidence. Gene structures were manually inspected and corrected where necessary. In total, 400 validated gene structures were obtained after manual curation.

Training and testing gene prediction parameters

Gene prediction parameter sets specific to Beta vulgaris were trained by the supervised machine learning algorithm implemented in the AUGUSTUS pipeline; we followed the general guidelines indicated in [30]. To select training genes and test genes we combined the set of SMRT-validated and manually validated genes based on initial AUGUSTUS predictions (the 665 SMRT-derived genes were not included at this stage). GenBank files containing the genes together with flanking intergenic regions were generated with gff2gbSmallDNA.pl. As detailed in the instructions on how to train AUGUSTUS, redundant genes and genes with a CDS length not divisible by 3 were removed. A gene was considered redundant if it had a sequence identity of 80 % or larger in an all-versus-all blat [31] protein sequence alignment. The remaining total of 2,542 validated genes were randomly subdivided into 2,000 training genes of which different subsets were used (Table 2) and 542 test genes for the calculation of sensitivity and precision using AUGUSTUS. In two separate approaches we: (1) trained parameters from scratch using the new_species.pl script; and (2) optimized existing A. thaliana parameters included in the AUGUSTUS package using the optimize_AUGUSTUS.pl script. The sensitivity and precision were calculated separately for UTRs, exons, and the entire transcript (referred to as ‘features’). The sensitivity was calculated by dividing the number of correctly predicted features (true positives) by all features in the test gene set (true positives and false negatives) and the precision by dividing the number of correctly predicted features (true positives) by all predicted features within the genomic regions of the test gene set (true positives and false positives).

Masking transposable elements

Repeats were identified and classified using the RepeatModeler software [32], v1.0.7, downloaded from [33]) applied on the sugar beet assembly RefBeet-1.2 and the spinach assembly Spinach-1.0.1. Repeat sequences classified as transposable elements were selected in both species (searching for ‘DNA, LTR, LINE, SINE, or Helitron’ in the RepeatModeler output). The repeat catalog of the sugar beet assembly RefBeet-1.1 had been manually curated [13], resulting in a number of repeats identified as transposons or retrotransposons which were ‘Unknown’ or ‘Simple repeat’ according to RepeatModeler. Those sequences were included in the combined sugar beet and spinach collection of transposable elements. The collection was used as input for RepeatMasker (further parameters: -gff -no_is -norna -nolow) to mask the sugar beet assembly RefBeet-1.2. Positions of transposable elements in RefBeet-1.2 were converted into repeat hint annotation interpretable by AUGUSTUS (source = RM in gff hint file). In one approach, genes were predicted on the repeat-masked genome, and in a second approach, the unmasked genome was used along with the positions as repeat hint information. To enforce the exclusion of repeat regions from gene models, the bonus factor for the prediction of repeat hints as non-exonic regions was increased from 1.01 to 1e + 10, and the repeat hint priority level was set to 6 (priority=6 in gff), which is above the default priority level of expression evidence of 4. The effect on the gene prediction accuracy of both masking approaches was evaluated by applying AUGUSTUS on genomic regions containing SMRT-derived genes (GenBank format) and by manual inspection of 200 genes of the genome-wide prediction in sugar beet (see ‘Manual quality assessment of gene predictions’). Both evaluations showed a slightly higher number of correctly predicted genes for hint masking when combined with additional expression hints and higher intron hint weighting (sensitivity hint masking 85.0 % versus sequencing masking 84.5 %). For gene annotation in Spinach-1.0.1, transposons were only provided to AUGUSTUS as repeat hint annotation.

Processing of expression evidence

KWS2320 sugar beet mRNA-seq data from stress conditions and their controls were quality-filtered and trimmed according to criteria described elsewhere [23], resulting in 526.9 million reads. mRNA-seq data that had been employed in predicting RefBeet-1.1 genes consisted of 616.3 million reads of which 396.9 were derived from genotype KWS2320. For this work, only mRNA-seq data of the KWS2320 genotype were utilized, in total 923.8 million quality-filtered reads. The reads were aligned to RefBeet-1.2 with blat, and processed and converted into hint information for AUGUSTUS as described previously [13].

mRNA-seq expression data were filtered using the custom Perl script ‘RNA-seq-noise-reduction.pl’ (Additional files 2, 3, and 4) to reduce the background noise and to facilitate the prediction of the most abundant isoform per locus correctly taking into account that some isoforms were not reported. The mRNA-seq coverage was reduced by 10 % of the local peak coverage (95 percentile, in 1-kbp windows). If the coverage difference between two overlapping intron hints (gff format) was greater than 90 % the intron hint with the lower coverage was removed. The mRNA-seq coverage (wig format) was reduced within boundaries of introns by 50 % of the adjacent exon coverage. Introns were considered if they were smaller or equal than 50 kbp in size and showed a coverage drop of at least 50 % at their exon-intron junction when comparing the coverage 10 bases upstream and downstream of each junction. To increase the weight of intron hints derived from mRNA-seq data (bonus factor 1), anchors were added for introns that were supported by at least 50 aligned mRNA-seq reads (source = M, bonus factor 1e + 100). The intron malus factor was adjusted from 0.34 to 0.001.

Aligned full-insert SMRT sequences (gff format) were converted into hints using the custom Perl script ‘derive-gene-models-from-PacBio.pl’ (Additional files 2, 3, and 4). In brief SMRT sequences were clustered based on their location and intron boundaries. Per location, the most abundant isoform was converted to exon hints, exon part hints (for terminal exons) and intron part hints, grouped together by using the group tag (gff column 9). The source tag of these hints was set to E when employed as EST hints and to M when employed as anchors. It was required that an isoform was represented by at least two SMRT sequences or by one SMRT sequence for which intron boundaries were confirmed by aligned mRNA-seq reads. The terminal exon positions were set to the median alignment start and stop positions of all SMRT sequences representing a selected isoform. In order to prevent AUGUSTUS from appending exons to SMRT hint groups, flanking intergenic hints of length one were added at a distance of 50 bases.

Calculation of an optimized sugar beet gene set

For the final parameter training, SMRT-validated genes, SMRT-derived genes, and manually validated genes were combined. In order to exclude additional genes from the flanking regions of the training genes, non-overlapping genes of the initial prediction were temporarily added to generate the GenBank file using gff2gbSmallDNA.pl. Redundant and problematic genes were removed as described above, and 2,794 training genes and 349 test genes were randomly selected using randomSplit.pl. Parameters were trained with optimized settings as evaluated above and using optimize_AUGUSTUS.pl. The extended training gene set resulted in further improvement of the ab initio performance (Table 4 and 7). The sugar beet reference gene set ‘BeetSet-2’ was generated on the assembly RefBeet-1.2 using AUGUSTUS version 2.7 with the newly generated B. vulgaris parameters, filtered sugar beet KWS2320 mRNA-seq hints, KWS2320 SMRT hints as well as Sanger and Roche/454 EST hints (Table 4; Additional files 2, 3, and 4). The options to predict and print UTRs were switched on, the gene model was set to ‘complete’ and no in-frame stop codons were allowed.

Transferring stable sugar beet gene identifiers

From the previous sugar beet gene set [13], the longest CDS per gene was mapped against RefBeet-1.2 using gmap (version 2012-07-20, -A -f gff3_match_cdna -B 5 -t 20). A stable identifier of a previously annotated gene was transferred to the BeetSet-2 gene with the longest summed CDS match length using the custom Perl script ‘transfer-stable-stable-identifier.pl’ (Additional files 2, 3, and 4). Whenever the CDS overlap was equally long to multiple genes, the gene order was tried to be preserved. Considering both evidence and non-evidence genes, in 558 cases one previously annotated gene overlapped multiple adjacent BeetSet-2 genes. In these cases the identifier was transferred to the BeetSet-2 gene with the longest partial CDS overlap. In 2,020 cases, multiple adjacent previously annotated genes matched best to a single BeetSet-2 gene. In these cases the identifier of the best reciprocal matching previously annotated gene was transferred. New stable identifiers were assigned to the remaining BeetSet-2 genes.

Manual quality assessment of gene predictions

Of the previous sugar beet gene set calculated with default A. thaliana parameters [13], 100 genes with 100 % evidence and 100 genes with 1-99 % evidence were randomly selected. The evidence for a gene was considered to be 100 % if all CDS, UTRs, and introns were supported by expression data. Genes and expression evidence were visualized in Gbrowse [34], and the correct structure was inferred from the combined Illumina mRNA-seq, SMRT, Sanger, and Roche/454 expression data. The number of correct structures was compared between the annotation in RefBeet-1.1 and BeetSet-2. Spinach genes predicted with A. thaliana or B. vulgaris parameters were inspected in the same way, except that genes at scaffold borders were kept.

Calculating number of stress-condition specific genes

Genes specifically supported by mRNA-seq reads from plants under stress conditions were determined by selecting those genes with no expression hint coverage under non-stress conditions and at least 90 % coverage of the CDS length under stress conditions (salt, heat, or light).

Applying Beta vulgaris parameters on spinach gene prediction

Newly generated B. vulgaris parameters were used together with 108.0 million quality filtered spinach mRNA-seq reads (RNA isolated from spinach leaves) to predict genes in the Spinach-1.0.1 genome assembly. Hint information was produced from these mRNA-seq data (isogenic to assembly) and processed in the same way as sugar beet mRNA-seq data.

Prediction of 1:1 orthologous genes between sugar beet and spinach

Evidence and non-evidence genes were combined. Per gene, the transcript encoding the longest protein was selected. Protein sequences of spinach and sugar beet were aligned to each other using blastp [35] (expect cut-off 1e-5, minimum alignment length 50 % of protein length). 1:1 orthologous genes were predicted using the reciprocal best blast hit approach [20].

Data access

Raw transcript sequencing data from sugar beet genotype KWS2320 and spinach genotype Viroflay were deposited in the NCBI Sequence Read Archive (SRA) with the following accession numbers: SRX674050 (sugar beet Pacific Biosciences SMRT sequences); SRR1508751, SRR1508753, SRR1508755, SRR1508756, and SRR1508758 (Illumina sugar beet mRNA-seq data from stressed and unstressed plants); SRX674044 (spinach mRNA-seq). Sugar beet and spinach gene models, assemblies, sugar beet SMRT validated genes used for training corrected SMRT sequences, and hints and anchors for Augustus gene prediction can be downloaded from [21]. This site further includes a genome browser for visualizing sugar beet and spinach annotation. Gene models were deposited on GenBank and have been assigned stable locus identifiers BVRB_1g000010 - BVRB_1g023380, BVRB_2g023390 - BVRB_2g047930, BVRB_3g047940 - BVRB_3g070810, BVRB_4g070820 - BVRB_4g097640, BVRB_5g097650 - BVRB_5g127170, BVRB_6g127180 - BVRB_6g156500, BVRB_7g156510 - BVRB_7g180820, BVRB_8g180830 - BVRB_8g202200, BVRB_9g202210 - BVRB_9g226140, and BVRB_000010 - BVRB_043090 for sugar beet, and SOVF_000010 - SOVF_217030 for spinach. Genome assemblies have been assigned accession numbers AYZV02000000 (spinach) and AYZS02000000 (sugar beet). GenBank gene sets and assemblies were adapted to be compliant with GenBank submission specifications.