Improved annotation with de novo transcriptome assembly in four social amoeba species

Background Annotation of gene models and transcripts is a fundamental step in genome sequencing projects. Often this is performed with automated prediction pipelines, which can miss complex and atypical genes or transcripts. RNA sequencing (RNA-seq) data can aid the annotation with empirical data. Here we present de novo transcriptome assemblies generated from RNA-seq data in four Dictyostelid species: D. discoideum, P. pallidum, D. fasciculatum and D. lacteum. The assemblies were incorporated with existing gene models to determine corrections and improvement on a whole-genome scale. This is the first time this has been performed in these eukaryotic species. Results An initial de novo transcriptome assembly was generated by Trinity for each species and then refined with Program to Assemble Spliced Alignments (PASA). The completeness and quality were assessed with the Benchmarking Universal Single-Copy Orthologs (BUSCO) and Transrate tools at each stage of the assemblies. The final datasets of 11,315-12,849 transcripts contained 5,610-7,712 updates and corrections to >50% of existing gene models including changes to hundreds or thousands of protein products. Putative novel genes are also identified and alternative splice isoforms were observed for the first time in P. pallidum, D. lacteum and D. fasciculatum. Conclusions In taking a whole transcriptome approach to genome annotation with empirical data we have been able to enrich the annotations of four existing genome sequencing projects. In doing so we have identified updates to the majority of the gene annotations across all four species under study and found putative novel genes and transcripts which could be worthy for follow-up. The new transcriptome data we present here will be a valuable resource for genome curators in the Dictyostelia and we propose this effective methodology for use in other genome annotation projects. Electronic supplementary material The online version of this article (doi:10.1186/s12864-017-3505-0) contains supplementary material, which is available to authorized users.


Background
Whole genome sequencing projects are within the scope of single laboratories. The Genomes OnLine Database [1] reports (as of 13 th May 2016) there are 76,606 sequenced organisms, of which 12,582 are eukaryotes. However, only 8,047 are reported as being complete. Annotation of gene models is a requirement for a complete genome [2]. There are several complementary strategies for achieving gene annotation in novel genomes including gene prediction [3,4], expressed sequence tag (EST) libraries [5] and RNA sequencing (RNA-seq) data [6]. Gene prediction methods are limited in the complexity of the gene models they are able to produce; alternative splice sites are unpredictable and untranslated regions (UTRs) have subtle signals [7]. EST libraries, if available, are usually fragmented and incomplete. RNA-seq data is dependent on good alignments to the reference. De novo transcriptome assembly is equally able to fulfil this function, although it can be computationally challenging [8][9][10]. Transcriptome assembly methods can be either reference-guided or reference-free [11,12]. Reference-guided methods have the advantage of simplifying the search space, but are dependent on the relevance, quality and completeness of the reference. Reference-free methods do not have any dependencies, but need to deal with sequencing errors sufficiently well to avoid poor assemblies [11][12][13]. We present the application of a de novo transcriptome assembly to four eukaryotic species: Dictyostelium discoideum, Polysphondylium pallidum, Dictyostelium fasciculatum and Dictyostelium lacteum. The genome of D. discoideum was published in 2005, it is 34 Mb in size and has been assembled into six chromosomes, a mitochondrial chromosome, an extra-chromosomal palindrome encoding ribosomal RNA (rRNA) and three 'floating' chromosomes [14]. The genome was generated via dideoxy sequencing and contigs were ordered into chromosomes by HAPPY mapping [14,15] and still contains 226 assembly gaps.
In contrast, the similar sized genomes of P. pallidum, D. lacteum and D. fasciculatum were sequenced more recently using both dideoxy and Roche 454 sequencing. Their assembly was assisted by a detailed fosmid map and primer walking, leading to only 33 to 54 gaps per genome, but are more fragmented with 41, 54 and 25 supercontigs, respectively [15,16]. The D. discoideum genome has been extensively annotated via the Dictybase project [17], whereas the gene models for P. pallidum, D. fasciculatum and D. lacteum, available in the Social Amoebas Comparative Genome Browser [18], are primarily based on computational predictions.
The social amoeba D. discoideum is a widely-used model organism for studying problems in cell-, developmental and evolutionary biology due to their genetic tractability allowing elucidation of the molecular mechanisms that underpin localized cell movement, vesicle trafficking and cytoskeletal remodeling as well as multicellular development and sociality. The social amoebas form a single clade within the Amoebozoa supergroup and are divided into four major taxon groups according to molecular phylogeny based on SSU rRNA and αtubulin sequences [15]. The four species under study here represent each of the four groups: D. discoideum (group 4), P. pallidum (group 2), D. fasciculatum (group 1) and D. lacteum (group 3). Genome annotations are not static and benefit from the application of additional evidence and new methodologies [7,19]. Therefore we present, for the first time, substantially updated annotations based on a de novo transcriptome assembly for the D. discoideum, P. pallidum, D. fasciculatum and D. lacteum genomes.

Sample preparation
Sequencing data were obtained from four RNA-seq experiments. The D. discoideum data were obtained from an experiment comparing gene expression changes between wild-type cells and a diguanylate cyclase (dgcA) null mutant at 22 h of development [20]. The P. pallidum data were obtained at 10 h of development in an experiment comparing wild-type and null mutants in the transcription factor cudA (Du, Q. and Schaap, P. unpublished results). In this experiment P. pallidum cells were grown in HL5 axenic medium (Formedium, UK), starved for 10 h on non-nutrient agar, and harvested for total RNA extraction using the Qiagen RNAeasy kit. The data for D. lacteum and D. fasciculatum were obtained from developmental time series [21]. For these series cells were grown in association with Escherichia coli 281, washed free from bacteria, and plated on non-nutrient agar with 0.5% charcoal to improve synchronous development. Total RNA was isolated using the Qiagen RNAeasy kit at the following stages: growth, mound, first fingers, early-mid culmination, fruiting bodies. D. lacteum RNAs were also sampled at three time points intermediate to these stages.

Illumina paired end sequencing
Total RNA was enriched for messenger RNA (mRNA) using poly-T oligos attached to magnetic beads and converted to a sequencing ready library with the TruSeq mRNA kit (Illumina), according to manufacturer's instructions and 100 basepairs (bp) paired-end sequenced using an Illumina HiSeq instrument. For the D. discoideum and P. pallidum samples, 1 μg of total RNA was used as starting material, with 4 ul of 1:100 dilution External RNA Controls Consortium (ERCC) ExFold RNA Spike-In Mixes (Life Technologies) added as internal controls for quantitation for the RNA-Seq experiment and sequenced at the Genomic Sequencing Unit, Dundee. In total there were 433 M, 413 M, 171 M and 319 M reads respectively for D. discoideum, P. pallidum, D. fasciculatum and D. lacteum.

Data processing and de novo transcriptomics assembly
The quality of the raw reads was checked with FastQC [22] and the reads were found to have high quality scores across their full length. No trimming of the data was performed, as aggressive trimming can negatively impact on the quality of assemblies [23]. All reads for each species were separately combined prior to de novo assembly. Being a more mature genome the D. discoideum data was used to verify the methodology, thereby giving a reference point for the other, less well characterised, species. Figure 1 shows a schematic of the overall workflow.
Trinity version 2013.11.10 [8] was used for de novo assembly, and normalisation of the read data was achieved with a kmer of 25 and aiming for 50x coverage of the genome. Following normalisation there remained 5.3 M, 8.3 M and 16.0 M read pairs in D. discoideum, P. pallidum and D. lacteum, respectively. D. fasciculatum reads were not normalised as there were fewer than the recommended 300 million reads as per the Trinity manual. Trinity was run on the normalised reads using the -jacard-clip parameter and setting -kmin-cov to 4 in an attempt to reduce the number of fused transcripts in P. pallidum only. In the other species the parameter made little difference. For the initial transcript set of D. discoideum and P. pallidum assemblies, any transcripts with BLAT (BLAST-like alignment tool) v35x1 [24] hits to the ERCC spike-in sequences were removed from the D. discoideum and P. pallidum assemblies. D. fasciculatum and D. lacteum were not cultured axenically and thus the samples were contaminated by their bacterial food source. In order to remove the bacterial contamination D. fasciculatum and D. lacteum, transcripts were filtered with the TAGC (taxon-annotated GC-coverage) plot pipeline [25]. TAGC determines for each contiguous sequence (contig) the proportion of GC bases, their read coverage and best phylogenetic match. With this information it is possible to identify which transcripts are mostly likely to be contaminants and removed. In order to remove the contamination, first all the transcripts were aligned to the BLAST 'nt' database using BLAST megablast. Using the trinity assembled transcripts, the BAM file of the reads mapped back to the transcripts and the transcripts to species mapping, non-target related transcripts were removed. The contaminant transcripts were differentiated on the coverage vs GC plots (see Additional file 1: Figure S1).
The normalised set of reads were aligned with bowtie (0.11.3, with parameters applied as per Trinity script alignReads.pl) [26] to the whole transcript set and the total number of reads matching to each transcript were stored (see Additional file 1: Figure S2 for the read distributions for each dataset).

Transcript refinement
Program to Assemble Spliced Alignments (PASA) v2.0.0 [27] was used to refine the Trinity transcripts into more complete gene models including alternatively spliced isoforms. Initially developed for EST data, PASA has been updated to also work with de novo transcriptome data. Using the seqclean tool available with PASA, all the transcripts were screened and trimmed for low complexity regions, poly (A) tails and vector sequences. GMAP (Genome Mapping and Alignment Program) [28] and BLAT [24] were used to align the transcripts to their respective genomes. Trinity transcripts that failed to align to the already existing genome in both GMAP and BLAT were removed as 'failed'. Remaining 'good' transcripts at this stage are termed the PASAaa dataset. Next, PASA takes existing annotations and compares them to the PASAaa dataset. PASA uses a rule-based approach for determining which transcripts are consistent or not with the existing annotation and updates the annotation as appropriate: new genes, new transcript isoforms or modified transcripts. PASAua is the term used for the PASA assembled transcripts after updating with existing annotation.

Assembly quality check
At each stage, the transcript datasets were assessed with Benchmarking Universal Single-Copy Orthologs, BUSCO, [29] and Transrate v1.0.0 [13] These methods take complementary approaches in assessing completeness and/or accuracy. Transrate v1.0.0 uses the read data and Fig. 1 The de novo transcriptomics assembly workflow. The reads are input at the top in green, all computational steps are in blue and all data or quality control outputs are shown in grey. PASA is the Program to Assemble Splice Alignments tool [27]. See main text for description of PASAaa and PASAua steps. BUSCO is Benchmarking Universal Single-Copy Orthologs [29] optionally the reference sequence as input. BUSCO defines a set of 429 core eukaryotic genes. These genes are used as a proxy for minimum completeness based on the assumption that a eukaryotic genome or transcriptome assembly should encode a large proportion of the core set of genes. The BUSCO (v1.1b1) tool uses hidden Markov models (HMMs), defined for each of the core genes in the set, returning whether there are complete or partial matches within the de novo transcripts. When run in genome mode, BUSCO additionally uses Augustus [3] to generate a predicted gene set against which the HMMs are tested. Transrate calculates the completeness and accuracy by reporting contig score and assembly score. Contig score measures the quality of the individual contig, whereas assembly score measures the quality of whole assembly.

Orphan RNAs
The full set of Trinity transcripts constitutes the best approximation of the assembly of transcripts expressed in the RNA-seq sequencing data. The transcripts were aligned against the existing genome and coding DNA (cDNA) references (from Dictybase (D. discoideum) and SACGB [18], (D. fasciculatum, D. lacteum and P. pallidum)) using BLAT. Any transcripts not matching the existing references were searched against the NCBI 'nt' database with BLAST [30] and with PSI-BLAST against the NCBI 'nr' database for the longest predicted ORF in any remaining transcripts without a match to 'nr'. This exhaustive search allowed the categorisation of 'annotated' (transcript with match to known genome and/or cDNA), 'known' (match to related species), 'artefact' (match to non-related species (non-Dictyostelid)) and 'putative novel' (remainder) datasets.
PCR and subcloning D. discoideum genomic DNA (gDNA) was extracted using the GenElute mammalian genetic DNA extraction kit (Sigma). Polymerase chain reaction (PCR) reactions were run for 30 cycles with 50 ng of gDNA and 1 μM of primers with 45 s annealing at 55°C, 2 min extension at 70°C and 30 s denaturation at 94°C. The reaction mixtures were size-fractionated by electrophoresis, and prominent bands around the expected size were excised, purified using a DNA gel extraction kit (Qiagen) and subcloned into the PCR4-TOPO vector (Invitrogen). After transformation, DNA minipreps of clones with the expected insert size were sequenced from both ends.

Results and discussion
De novo transcript assembly Table 1 shows a summary of the Trinity output for the D. discoideum, P. pallidum, D. lacteum and D. fasciculatum de novo transcriptome assemblies. Overall, the raw assemblies are similar in terms of total transcripts, GC content, and contig N50 or E90N50 (N50 for the top 90% expressed transcripts). D. discoideum is slightly anomalous in N50, E90N50, mean length and transcripts ≥ 1,000 bp with all features being smaller than the other three assemblies. The mean length over all the annotated Dictybase coding sequences is 1,685 bp which is substantially larger than in the assembled transcripts (867 bp) suggesting that the D. discoideum transcripts are fragmented. Figure 2 shows the distribution of transcript lengths for D. discoideum, P. pallidum, D. fasciculatum, D. lacteum (cyan) when compared to the available cDNA datasets (magenta). The D. discoideum cDNAs are manually curated, whereas the others are predicted. The transcript sets are enriched in short transcripts (<1000 bp) as compared to their cDNAs with the effect being most marked in D. discoideum, D. fasciculatum and D. lacteum (Fig. 2a,c and d). The P. pallidum assembly is more similar to its cDNA reference dataset (Fig. 2b). Interestingly, the longest assembled transcript in D. discoideum (21,679 bp) was found to be approximately half of the mitochondrial chromosome. We speculate that as the mitochondrion is gene rich and highly expressed, Trinity was unable to resolve overlapping reads from adjacent genes thereby joining them all into one 'supercontig'. The subsequent steps in the assembly were performed with PASA [27] which uses reference genome and transcript datasets to generate a refined and updated transcriptome assembly. The first stage takes the transcriptome assemblies, aligns them against the genome and clusters them into gene structures according to their genome alignments. Any transcripts, which do not align adequately to the genome are filtered out by PASA, under the assumption that they are misassemblies. This dataset will be referred to as 'PASA annotated assemblies' , PASAaa. In unfinished and complex genomes, it is possible that there are missing gene loci in the genome reference. The missing loci may appear in a de novo transcriptome assembly and would be filtered out by PASA. The second stage uses the aggregated and filtered set of transcripts to refine the existing annotations for each of the species. At this stage, the gene models are updated with new or extended UTRs, new alternatively spliced isoforms are added, and introns are added or removed. New genes are identified, and existing genes are split or merged as required by the de novo assembly data. This dataset is referred to as 'PASA updated annotations' , PASAua. Table 2 shows the results of each stage of the assembly workflow from Trinity to each of the PASA steps and compared to the existing set of gene models from DictyBase (D. discoideum) or Augustus predictions (P. pallidum, D. lacteum and D. fasciculatum). It is clear that at each stage the assemblies become more similar to the existing gene models ( Table 2). For example, in all the species the total number of transcripts was 3-4-fold larger in the Trinity data than in the existing annotations. Although de novo assembly has the potential to identify novel genes and transcripts, a 3-fold increase is unlikely. By the end of PASAua, the transcript counts were within 1,500 of the existing models, with D. fasciculatum, D. lacteum and P. pallidum having more genes than in their Augustus-predicted models, and D. discoideum having 760 fewer genes than in the DictyBase-curated models. This is to be expected as the gene prediction algorithms are unlikely to have found all transcripts, whereas the D. discoideum curated set will include genes expressed under certain conditions only (e.g. developmental time points) that were not part of the experiment included here. Mean transcript lengths increased through the workflow. In particular, for D. discoideum, the mean Trinity transcript length was 871 bp and the final PASAua length was 1,787 bp indicating that the high fragmentation observable by an excess of short transcripts (Fig. 2) has been reduced. Similarly, the total number of identified exons was reduced from the initial Trinity dataset. Overall, the initial Trinity assemblies have been refined from a fragmentary and redundant dataset to a more full-length and less redundant set of transcripts, which are more similar to the existing reference datasets in terms of total transcript counts, mean length, number of exons and exons per transcript (Table 2).

Quality assessment
Benchmarking Universal Single-Copy Orthologs (BUSCO) and Transrate are tools which allow the assessment of completeness and accuracy of transcriptome assemblies. A set of 429 core eukaryotic genes (CEGs) was defined by BUSCO, for the purpose of assessing completeness in eukaryotic genomes [29]. CEGs are conserved across taxa and the majority should be present in the majority of eukaryotic species. A large fraction of missing BUSCO genes could be indicative of an incomplete assembly. Figure 3 shows the comparison of complete and partial BUSCO matches in all four species for the genome reference, Trinity assembly, PASAaa refined transcripts and PASAua updated annotations. In the ideal situation all the BUSCOs would be detected in an assembly, however high sequence divergence or absence in the species will give a lower maximum detection level. The whole genome BUSCO score represents the upper limit for any of the assemblies. All except the PASAaa datasets have >80% complete or fragmented BUSCOs and are close to the whole genome count, suggesting the transcriptome assemblies are nearly complete. It is noticeable, that the number of identified BUSCOs is consistently lower in the PASAaa data for all four species (Fig. 3). This drop is due to the strict PASA filtering during transcript assembly. PASAaa only retains transcripts, which align to the reference with 95% identity and 90% length coverage. Manual checking of the BUSCOs that are identified in the Trinity data, but not in PASAaa reveals that they all are labelled as failed alignments. This suggests that either BUSCO is overly permissive in defining the orthologues or that PASAaa is overly aggressive in filtering transcripts. PASAua appears to 'rescue' this behaviour, presumably by including good annotations for genes that are poorly assembled in the Trinity data.
Transrate assesses transcript quality by calculating several contig-level metrics based on the input RNA-seq data, and measures how well the read data support the contigs. Contigs are scored individually and then combined into an overall assembly score which ranges from 0 to 1. An optimal score is also reported, which predicts the best potential assembly score achievable by removing the worst scoring contigs in the dataset. An assembly score of 0.22 and optimised score of 0.35 were found to be better than 50% of 155 published de novo transcriptome assemblies [13]. A high Transrate score with a small improvement in the optimal score indicates a good de novo assembly, which is unlikely to be improved without further data or information. Figure 4 compares the distribution of Transrate contig scores from the Trinity assembly, PASAaa refinement, PASAua update and reference transcript/coding sequence (CDS) datasets for each of the four species. In contrast to the BUSCO data, the PASAaa data shows an improvement in Transrate contig scores when compared to the raw Trinity output meaning that the PASAaa transcripts are more consistent with the data, confirming that perhaps BUSCO is too permissive when assigning orthologues rather than PASAaa being too aggressive with its filtering. Notably the reference sequence datasets ('CDS' Fig. 4) for D. discoideum and P. pallidum, show a lower median score than the PASAua data, indicating that PASAua is working well in combining the data with the existing annotations. There is little difference in D. lacteum. In D. fasciculatum the CDS data shows the best Transrate score of any of the assemblies. Figure 5 compares the Transrate assembly scores and optimal scores between PASAua and the annotated CDS over the four species. The assembly scores range from 0.31 (D. fasciculatum) to 0.42 (D. discoideum) and the optimal scores range from 0.32 (D. fasciculatum) to 0.53 (D. discoideum). It is clear that PASAua has better Transrate scores (Fig. 5a filled circles) than the annotated CDS (Fig. 5a filled triangles), except for D. fasciculatum, with all the PASAua assemblies scoring better than 50% of published transcriptome assemblies (Fig. 5a dotted black line). The optimal scores for PASAua are also all better than 50% of published transcriptome assembly data (Fig. 5a dotted cyan  (0.32) is small (Fig. 5a green filled circles), suggesting that there is little improvement to the assembly possible given the read data for this species. Using the optimal score, Transrate defines a set of 'good' contigs which best fit the data. The proportion of PASAua good contigs ranges from 79.9% (D. discoideum) to 97.2% (D. fasciculatum) which, for all species, is a higher proportion than the annotated CDS (Additional file 1: Table S3). Transrate additionally has a reference-based measure, which aligns the transcripts to the reference protein sequences and the results are shown in Fig. 5b. The y-axis in Fig. 5b shows the proportion of reference protein sequences covered with transcript sequences at several thresholds (25,

Interpretation
What does an RNA-seq-based de novo assembly achieve when there is an already existing annotation either manually curated or generated via prediction? Is it worth it? Table 3 details the results following PASA refinement of the existing gene models. Despite being a manually curated genome, the D. discoideum gene models where extensively modified by PASA with 7,182 being updated. Most of the updates in D. discoideum (6,750, 94%) are the result of UTR additions at 5′ and 3′ ends of genes, which were mostly missing in the existing models. The assemblies in the other species have a similar number of updates, but UTR-only updates to transcripts are a smaller fraction of the total. 187 new alternatively spliced transcripts, in 170 genes, were identified in D. discoideum (Table 3). There are currently 70 alternatively spliced transcripts, in 34 genes, annotated in Dictybase so this new data represents a 2.7-fold increase in the number alternatively splice transcripts and a 5-fold increase in genes. This number in D. discoideum could be an underestimate as the D. fasciculatum, D. lacteum and P. pallidum assemblies all have~1000 alternate splice isoforms. Figure 6 provides examples, in each of the four species, of changes to the transcript models determined by PASA that are well supported by all the data. Each panel highlights a different type of change to the reference model. Gene DDB_G0295823 has a single transcript (DDB0266642 Fig. 6a) with two exons and a single intron. The RNA-seq data (brown), Trinity assembly (purple) and PASAaa refinement (red) identifies extensions to the model, adding 5′ and 3′ UTRs to the annotation (green, narrow bars). The Trinity transcript (purple) is on the opposite strand to the reference transcript (black) and is corrected by PASA (red & green). The example in P. pallidum (Fig. 6b) shows three new alternatively spliced products of the gene (Fig. 6b green bars 1, 2, 3 labels). The three new models have the same coding region, but differ in their 5′-UTRs: two with differently sized introns and one without an intron. The new models also include a longer second coding exon (Fig. 6b arrow) Reference Coverage Proportion of reference proteins Fig. 5 Transrate assembly scores and reference coverage metric. a Compares the Transrate [13] assembly score and the optimised score between the CDS and PASAua [27] datasets in the four species (Ddis: D. discoideum, Ppal: P. pallidum, Dfas: D. fasciculatum, Dlac: D. lacteum). The dotted lines represent the Transrate scores that would be better than 50% of 155 published de novo transcriptomes as found by Smith-Unna and co-workers [13]: 0.22 overall score (black horizontal dotted line) and 0.35 optimal score (cyan horizontal dotted line). b The proportion of reference protein sequences covered by transcripts in the CDS and PASAua datasets by at least 25%, 50%, 75%, 85% and 95% of the reference sequence length  which increased the sequence of the protein product by 9 amino acids. Figure 6c shows an example, in D. fasciculatum, where an alternatively spliced transcript alters the protein product. The alternatively spliced isoform (Fig. 6c, labelled 1) removes the first intron and extends the 5′-UTR when compared to the updated gene model (labelled 2). The CDS is shortened by 45 amino acids with the use of alternate start site, but the rest of the protein is identical. In the RNA-seq data it appears that this new alternative transcript is not the dominantly expressed isoform in the context of the whole organism. The final example is the merging of two D. lacteum genes into one (Fig. 6d). The black bars show two distinct genes (DLA_11596 and DLA_04629), but the RNA-seq data (brown) and the Trinity assembly (purple bars) show uninterrupted expression across the intergenic region between the two genes (arrow). The PASA refinement and re-annotation (red and green bars) encapsulate the expression as a contiguous region with the coding region being in-frame over the two existing gene models. The annotation for the upstream DLA_11596 gene in SACGB [17] gives its best bi-directional hit in Uniprot/ TrEMBL as gxcN in D. discoideum (DDB0232429, Q55 0V3_DICDI). gxcN codes for a 1,094 amino acid protein where DLA_11596 codes for a 762 amino acid protein and the pairwise alignment of DLA_11596 with DDB0232429 shows no overlap over the C-terminal 300 residues. The PASAua gene fusion of DLA_11596/DLA_04629 (Fig. 6d) codes for a longer, 1,029 protein which aligns across the full length of DDB0232429 in a pairwise alignment. We suggest that the existing gene model, DLA_11596, is a truncated form of a D. discoideum gxcN orthologue and that the fusion with the downstream DLA_04629 gene represents the more accurate gene model. Given that D. discoideum has been extensively studied and the annotation curated by Dictybase, it is of note that our pipeline identified putative changes which altered the protein sequence of 554 genes (4.5% of total reference models) ( Table 3). D. discoideum has been the focus of many functional studies including about 400 deletions in genes that are required for normal multicellular development [16]. Comparing the 554 D. discoideum genes with modified proteins to the developmentally essential genes, we found 16 genes (2.9%) that overlapped (see Additional file 2: Figure S3 for domain diagrams). Out of the 16, nine are either truncated or extended at the N-or C-terminal. In the remaining seven proteins, there is loss or gain of exons. Five proteins were updated with additional exons: DDB_G0268920, DDB_G0269160, DDB_G0274577, DDB_G0275445 and DDB_G0277719, and two proteins have an exon deletion: DDB_G0271502 and DDB_G0278639.
Investigating these protein changes in more detail revealed some errors in the underlying genome sequence, which resulted in some unusual gene models. Figure 7 shows clcD (chloride channel protein, DDB_G0278639) as an example. In the domain architecture of clcD, there are two CBS (cystathionine beta-synthase) domains present at positions 827-876 and 929-977 in the transcript sequence. In the updated sequence the protein is truncated and these two domains have been removed. This is likely to be incorrect since all eukaryotic CLC proteins require the two C-terminal CBS domains to be functional [31]. How did this change occur in the de novo transcript assembly? In the existing annotation, there is an impossibly short two-base intron between the CLC domain and first CBS domain. Splicing requires a two-base donor and a two-base acceptor at either end of the splice site meaning at least four bases are required, not including any insert sequence. Careful investigation of the RNA-seq genome aligned reads reveals a singlebase insertion immediately after the intron in 22 out of 23 reads overlapping the region (yellow inset, Fig. 7). The RNA-seq data turns the two-base intron into a three-base, in-frame codon inserting an isoleucine into the protein sequence and retaining the CBS domains. By implication there is a missing base in the genome reference, which interrupts the open reading frame with a premature stop upstream of the CBS domains (arrow, Fig. 7). PASA cannot deal with missing bases in the reference and erroneously truncates the, now out-of-frame, coding region four codons downstream of the missing base at a TGA stop codon. It also cannot create an impossible intron, which a human annotator presumably added in order to keep the transcript in-frame and retain the conserved CBS domains. PASA did make an error updating this gene, but it does not seem possible for it to have dealt with the missing base any other way.
Inspection of all the D. discoideum gene models identified 119 sites in 102 genes with introns shorter than 5 bp (see Additional file 3: Table S1). Of these genes, five have three tiny introns each. Four of them are either in poorly expressed genes or in poorly expressed regions within genes. One gene (DDB_G0279477), however, is well expressed across the full length. The gene contains two 3 bp introns and one 1 bp intron. The two 3 bp introns contain a TAA sequence encoding a stop codon, but according to the RNA-seq data the codons should be TTA (Leu) with evidence from 56 and 33 reads in the two sites, respectively, 100% of which contain the TTA codon. The 1 bp intron region is covered by 38 reads and one would not expect to see introns in RNA-seq data, by definition, it does seem highly unlikely for a 1 bp intron to exist given our current knowledge of mRNA splicing: canonical GU-AG dinucleotides and a branch point >18 bp upstream from the 3′ splice site. For this gene, there are clear errors in the genome sequence, which have lead to the creation of an erroneous gene model to compensate for them. It is arguable that none of the 119 < 5 bp introns are genuine but are artificial constructs to fix problems with the gene models. We recommend that gene annotators revisit these genes and consider updating the models [7,32] and the underlying genome using RNA-seq data as evidence [33,34].
In addition to what we have shown here, it would be possible to use the RNA-seq data to directly improve the genome assembly of the four dictyostelid species mentioned herein. Xue et al. [35] have shown with their 'L_RNA_Scaffolder' tool that improved scaffolding of complex genomes such as human and zebrafish is possible with RNA-seq indicating the feasibility in more gene dense species.
The protein changes in D. fasciculatum, D. lacteum and P. pallidum number in the thousands (Table 3) highlighting that computational gene prediction is only a first step in annotating a genome. A reliable genome annotation requires evidence from many sources of information [19]. The types of protein changes seen in these three species range from inappropriately fused or split genes (see Fig. 6 bottom panel for an example) via insertions/deletions to changes in protein coding start/ stop codons positions resulting in extended or truncated coding sequences. All the PASAua outputs are in the form of GFF files viewable within any genome browser. We have made an IGB Quickload server available for easy browsing of the data (http://www.compbio.dundee.ac.uk/ Quickload/Dictyostelid_assemblies).
In the D. discoideum, D. fasciculatum, D. lacteum and P. pallidum datasets 44, 19, 21 and 175 novel putative genes were identified by PASA respectively (Table 3). These novel genes are in genomic loci with no current annotated gene model or where an existing model is substantially modified. The 44 D. discoideum novel genes, defined by 47 transcripts, were examined by eye in IGB [36] against all known D. discoideum reference datasets, including predicted gene models (see Additional file 4: Table S2). Of the 47 transcripts, 8 are novel alternate splice transcripts (Additional file 4: Table S2). Although 'novel' suggests there is no existing annotation at the locus of interest, if a gene update is sufficiently different from the reference gene model, PASA may consider that locus as a novel gene. In most of these cases the new transcript represents a corrected model for a previously computationally predicted gene. Many of the predicted gene models were annotated in Dictybase as pseudogenes and were originally ignored by PASAua, which only considers protein coding genes. Fragments of the pseudogenes do encode ORFs and PASA has reported them as being novel genes (Additional file 4: Table S2), but it is not possible to be sure whether the protein products are expressed in vivo with this data. Out of the 47, it appears only 6 are truly novel as they do not overlap any previously annotated transcripts: novel_model_13, novel_model_23, novel_model_30, novel_model_31, novel_model_38 and novel_model_39. All except novel_model_23 have a sequence match to existing genes, suggesting that they are paralogues. The longest novel unannotated model is 510 AA in length (novel_model_31) and appears to be a duplicate copy of the leucine rich repeat protein lrrA present on the chromosome 2.
Notwithstanding the large number of updates to the existing D. discoideum annotations it is clear from Table 3 that there are substantially more changes in the other three species. In particular, the numbers of modified protein sequences are 4, 9 and 10-fold larger in P. pallidum (2,252), D. lacteum (4,741) and D. fasciculatum (5,393), respectively. Similarly, there are 7, 5 and 6-fold more novel alternate splice isoforms in the three species, respectively. For P. pallidum (1,321), D. lacteum (1,088) and D. fasciculatum (842), the gene models were predicted with Augustus (G. Glöckner, personal communication) which, given the updates found with PASA, suggests that although the predicted gene models are in the correct locus, many are inconsistent with empirical RNA-seq evidence. With respect to novel genes annotated by PASA, it is notable that D. fasciculatum and D. lacteum have fewer than either D. discoideum or P. pallidum. It is unclear why this would be. Many genes were inspected by eye with IGB [36] and overall the annotations appear appropriate, but there are many occasions where human intervention would make further improvements.

Orphan RNAs
As mentioned above, PASA requires that transcripts align to the genome before it can consider them for further analysis. It makes sense to use the genome as a filter for valid transcripts, however this makes the assumption that the genome is complete. Any gaps in the genome that include genes will result in filtering out perfectly valid transcripts.
To determine whether this has happened here, we isolated the transcripts that did not align to the genome and used a process of elimination to identify those transcripts that could be genuine. Table 4 breaks down the number of orphan RNAs and whether they match nondictyostelid genes ('artefact'), genes in other dictyostelids ('known') or neither ('novel'). D. fasciculatum and D. lacteum have far more 'novel' non-genome transcripts (6,559 and 6,465, respectively) than D. discoideum (69) or P. pallidum (26). This is likely due to the fact that these species, which were cultured on bacteria, contain chimeric misassemblies of bacterial and dictyostelid transcripts. Despite this, they still have 525 and 945 'known' transcripts which have sequence matches to other Dictyostelids, higher than seen in D. discoideum (14) and P. pallidum (82). These transcripts are probably the best candidates for experimental assessment as genuinely non-genome transcripts.
We further investigated the 69 D. discoideum 'novel' transcripts with a more sensitive PSI-BLAST search on their longest ORFs and queried their cognate proteins for functional domains using SMART. Table 5 shows the 11 most interesting hits based on the sequence match, read count and ORF length. They are all well expressed and have ORF lengths consistent with functional proteins. Three novel transcripts (comp4660_c0_seq1, comp4660_c4_ seq1 and comp5569_c2_seq1) show similar sequence matches to DDB_G0292950 via PSI-BLAST searching, in spite of very low sequence similarity between them. DDB_G0292950 codes for a hypothetical protein which is not conserved in other dictyostelids and is poorly expressed (RPKM <1) at all time points in dictyExpress [37]. The three transcripts match across different parts of DDB_G0292950 indicating that they are different parts of the same larger gene. All transcripts identified in Table 5 were selected for experimental validation via PCR amplification, 8/11 were confirmed.
The comp5787_c28_seq1 transcript has putative homologues in D. fasciculatum, D. purpureum and P. pallidum as shown in Fig. 8. Sequence conservation is high as well as conservation of the Importin-beta N-terminal domain (IBN_N) and HEAT-like repeat (HEAT_EZ) domain architecture although the D. discoideum sequences appears to have an additional HEAT repeat domain (Fig. 8).

Genomic cloning of orphan Dictyostelium discoideum mRNAs
The newly assembled transcripts that could not be mapped onto the genome are either contaminants or genuine mRNAs for which the genomic counterpart is in an assembly gap of the genome. To investigate the latter option, we used PCR to attempt to amplify the genes from D. discoideum genomic DNA (gDNA). Oligonucleotide primers were designed to amplify regions of about 0.5 -1.4 kb of 11 transcripts (Additional fie 1: Table S4). The amplified size can however be larger due to the presence of introns. For eight transcripts, corresponding gDNAs could be amplified, but for two genes two transcripts were part of the same gene (Additional file 2: Figure S3). The six genes in total were all protein coding genes. For three transcripts, comp470_c0_seq1, comp4678_c1_seq1 and comp2066_c1_seq1 no PCR products were obtained, but the first two transcripts contained multiple stop codons in all reading frames and are likely assembly errors. The amplified PCR products were sub-cloned and sequenced from both ends. Sequences were assembled and aligned with the transcript sequence. Apart from just a few mismatches, the transcript and gDNA sequences were identical (Additional fie 1: Table S5). Only one amplified fragment contained introns (Additional file Table 1: S5). Six out of seven of the protein coding orphan transcripts therefore had a counterpart in the genome. Overall, deciphering the genomes of organisms is a key step in being able to probe their biology. With the advent of high-throughput sequencing technologies this has become a simpler problem to solve. Yet it is still not trivial to finish a genome assembly without any gaps [38]. The genome sequence on its own, however, imparts very little functional information and requires annotation of genes, transcripts and regulatory regions to be scientifically useful [7]. Many gene annotation methods are dependent on either homology to related species [30,39] or via gene finding prediction algorithms [40,41] or ideally both. However, the first method will miss all unusual or species-specific genes, while both methods fall short of accurately predicting intron-rich genes, genes with alternative or non-canonical splice sites or genes with very short exons. The ability to generate a whole transcriptome for a given species and use it to empirically annotate the genome has the power to confirm and correct any errors introduced with other methods. This has been achieved with expressed sequence tags (ESTs) in the past [42], but now can be performed with RNA-Seq short read data [32]. This evidence-based methodology is non-trivial and is not perfect. There are examples where the data is not adequately represented in the final transcript set when interpreted by the human eye. In addition, PASA only defines protein-coding genes meaning that all non-coding RNAs (ncRNAs) will be ignored and will not be in the final annotation unless already identified in the reference. Identifying ncRNAs is difficult as they have no obvious products and well-defined sequence features [43]. This does not negate their importance or relevance to the Dictyostelia.

Conclusion
In this study, we present a de novo transcriptome assembly in four social amoeba species for the first time and with these data we have: Created a final set of of 11,523 (D. discoideum), 12,849 (P. pallidum), 12,714 (D. fasciculatum) and 11,315 (D. lacteum) transcripts. Substantially updated the existing transcript annotations by altering models for more than half of all the annotated transcripts.
Identified changes to thousands of transcripts in the predicted gene models of P. pallidum, D. lacteum and D. fasciculatum many of which affect the protein coding sequence. Identified and validated six novel transcripts in D. discoideum.
Putatively identified dozens to hundreds of novel genes in all four species. Identified errors in the genome sequence of at least two D. discoideum genes (clcD and DDB_G0279477). With the possibility of, at least, another 104 genes having sequence errors. Found hundreds of putatively alternatively spliced transcripts in all species, something which has not been identified before in P. pallidum, D. lacteum or D. fasciculatum.
By combining methodologies we now have a better and more complete description of the transcriptome for these four species. This is not an end-point, however, but a further step towards fully finished genomes. More data and more manual refinement will be required to improve the annotations further. Fig. 8 Protein sequence for comp5787_c28_seq1 alignment with homologues from D. fasciculatum, D. purpureum and P. pallidum. a Jalview [45] multiple sequence alignment together with Jpred secondary structure prediction and its associated confidence, 'JNETCONF' [46]. Green arrows represent extended strands and red bars represent helical regions. In the alignment IBN_N (purple) and HEAT_EZ (red) domains are highlighted. b MrBayes [47] phylogenetic tree annotated with SMART [48] domain architectures determined. Each amino acid in the multiple alignment is coloured according to the clustalx [49] colour scheme