dbWFA: a web-based database for functional annotation of Triticum aestivum transcripts

Abstract

The functional annotation of genes based on sequence homology with genes from model species genomes is time-consuming because it is necessary to mine several unrelated databases. The aim of the present work was to develop a functional annotation database for common wheat Triticum aestivum (L.). The database, named dbWFA, is based on the reference NCBI UniGene set, an expressed gene catalogue built by expressed sequence tag clustering, and on full-length coding sequences retrieved from the TriFLDB database. Information from good-quality heterogeneous sources, including annotations for model plant species Arabidopsis thaliana (L.) Heynh. and Oryza sativa L., was gathered and linked to T. aestivum sequences through BLAST-based homology searches. Even though the complexity of the transcriptome cannot yet be fully appreciated, we developed a tool to easily and promptly obtain information from multiple functional annotation systems (Gene Ontology, MapMan bin codes, MIPS Functional Categories, PlantCyc pathway reactions and TAIR gene families). The use of dbWFA is illustrated here with several query examples. We were able to assign a putative function to 45% of the UniGenes and 81% of the full-length coding sequences from TriFLDB. Moreover, comparison of the annotation of the whole T. aestivum UniGene set along with curated annotations of the two model species assessed the accuracy of the annotation provided by dbWFA. To further illustrate the use of dbWFA, genes specifically expressed during the early cell division or late storage polymer accumulation phases of T. aestivum grain development were identified using a clustering analysis and then annotated using dbWFA. The annotation of these two sets of genes was consistent with previous analyses of T. aestivum grain transcriptomes and proteomes.

Database URL:urgi.versailles.inra.fr/dbWFA/

Introduction

Triticum aestivum (L.), common wheat or bread wheat, is one of the most important staple crops in the world. It is cultivated worldwide and provides >20% of the calories and proteins in the human diet (http://faostat.fao.org). Although ongoing sequencing efforts have already produced important genomic resources (1–4), the complete sequencing and annotation of the hexaploid (2n = 6× = 42, AABBDD) T. aestivum genome has yet to be achieved. A first version of the genome of the bread wheat cv. Chinese Spring has recently been published (4), providing the scientific community with highly valuable genomic and evolutionary information, which will facilitate genome-wide analysis of bread wheat. However, because of the low-coverage (5-fold) shotgun sequencing method used in this project, this resource does not represent a high-quality draft of the wheat genome in terms of sequence completion, quality and annotation. Genome-wide analysis of gene expression by means of expression microarrays or transcriptome sequencing (RNA-Seq) is now being adopted in T. aestivum (5, 6), but analysing such large datasets requires extensive annotation efforts. Data fragmentation and technical and semantic heterogeneity can severely limit the efficient extraction and interpretation of biological data (7, 8).

More and more genomic information is becoming available for T. aestivum research. Various resources and associated tools grant the user a structural overview of expressed sequence tag (EST) (ITEC, http://avena.pw.usda.gov/genome/) (9, 10) or bacterial artificial chromosome clone libraries (11, 12) for instance (13–15). Important initiatives are underway to facilitate the breeding of improved Triticeae varieties. The TriticeaeGenome project (www.triticeaegenome.eu) grants access to comprehensive information extracted from experimental data to provide a better understanding of Triticeae genomes (16). The global database GrainGenes (http://wheat.pw.usda.gov/) provides a variety of services and bioinformatics tools for the Triticeae and Avena sativa research communities. The HarvEST database (http://harvest.ucr.edu/) (17), dedicated to several crop species, including T. aestivum and Hordeum vulgare, provides access to curated EST assemblies, comparative analysis tools and links to orthologues in related model plant species. Together, these resources compile and cross-reference a great deal of information on physical and genetic mapping, markers, sequence variations and quantitative trait loci. To some extent, they also provide information leading indirectly to predicted gene product functions, but none of them is focused on functional gene annotation, and it is necessary to navigate through numerous unlinked resources to extract functional information.

Recently, pipelines for the automated annotation of genomic sequences of T. aestivum and related species have been developed (3, 18). These pipelines are based on the prediction of gene models within genome sequences, so they are not able to functionally annotate sequences originating from transcriptome sequencing like ESTs. Because no reference genome sequence is yet available for T. aestivum, a massive sequencing effort has produced more than a million ESTs (http://wheat.pw.usda.gov/genome/). To deal with the high level of redundancy of this resource, these sequences were clustered (i.e. overlapping and partial polyA-tailed expressed sequences are grouped) to provide a reference set of unique expressed genes, NCBI UniGenes (http://www.ncbi.nlm.nih.gov/UniGene). The assembly conditions used to build NCBI UniGenes make it the most comprehensive coding DNA (cDNA) assembly available to date, and UniGene assemblies have been used as a reference set of sequences for many species. An additional effort was made to construct full-length cDNA sequences that were included in the TriFLDB database (19). Full-length cDNA sequences are most commonly used in genome annotation as a resource for cross-species comparative analyses. Currently, TriFLDB is the most reliable source of full-length cDNA sequences in T. aestivum. TriFLDB includes annotations based on homologies found by searching protein databases, extensive Gene Ontology (GO) annotations and InterProScan results. Recently, a new collection of nearly 1 million ESTs, assembled into contigs and singlets, was annotated with GO terms (20), but meaningful prediction of gene function requires more than one system of annotation.

After the sequencing of the first plant genome, Arabidopsis thaliana in 2000 (21), several plant sequencing projects have been successful. The sequenced genomes most closely related to the T. aestivum genome are those of Oryza sativa ssp. indica (22), Oryza sativa ssp. japonica (23), Zea mays (24), Glycine max (25), Sorghum bicolor (26), Brachypodium distachyon (27) and Hordeum vulgare (28). Both structural and functional annotation resources for these species are developing steadily. One of the most effective methods to annotate a transcript is to find its orthologous counterparts in well-annotated closely related genomes (29, 30). Although H. vulgare (L.) would be expected to be the most useful reference because it is more closely related to T. aestivum, comprehensive and high-quality annotations of gene function are only available for O. sativa and A. thaliana (2), essentially owing to the lengthy and accurate annotation efforts undertaken.

To annotate ESTs or transcripts using sequence homology, it is necessary to navigate through unrelated databases. Some tools using this homology approach have been developed. For instance, Blast2GO (31) can be queried using T. aestivum sequences to give GO results. ONDEX (7), developed with the challenges of functional annotation of T. aestivum genome in mind, combines data integration from various sources and various mining methods, including graph-based analyses, to annotate wheat gene functions according to a wisely chosen set of annotation standards. However, ONDEX does not provide easy access to workable static results that are often required in research. To fill this gap, we developed dbWFA, an open-access database relating the T. aestivum UniGene set and the full-length cDNA sequences from TriFLDB to A. thaliana (TAIR10) (32) and O. sativa (pseudomolecules version 7.0) (33) annotation through BLAST (34) results. dbWFA also includes the inventory database of T. aestivum transcription factors (wDBTF) (35) and hand-curated gene families (36). As an all-in-one interface for the annotation of T. aestivum sequences, dbWFA will be useful to the researchers working on T. aestivum and more generally on cereals, particularly for comparative cereal genomics and functional genomics. The web implementation of dbWFA provides an easy-to-use interface to annotate transcript sequences from T. aestivum, with functional information from multiple pervasive annotation systems. Here, the use of dbWFA is illustrated with several query examples, and the quality of the annotation method is assessed by comparing the MapMan bin annotation of all the transcripts of the T. aestivum NimbleGen 40 k microarray (37) with that of A. thaliana and O. sativa. The use of dbWFA is further illustrated by analysing the annotation of 433 genes specifically expressed during either the early cell division or the late storage polymer accumulation (SPA) phases of grain development.

Data Content, Database Architecture and Web Interface

Five functional classification/annotation systems were integrated (Figure 1) in dbWFA to offer a fast and efficient functional annotation tool for T. aestivum UniGenes:

• GO (http://www.geneontology.org) (38), a non-redundant structured hierarchy of ontologies, which is the most widely used functional annotation system in bioinformatics. The GO project provides an efficient annotation standard that can be applied to numerous species. It is built on a controlled vocabulary of terms for describing gene function. dbWFA includes GO annotation data (OBO version 1.2) for A. thaliana and O. sativa.

• Plant Metabolic Network (PMN; http://www.plantcyc.org) (39), which provides a broad network of curated databases on primary and secondary plant metabolism, including pathways, enzymes, genes, compounds and reactions from several plant species. dbWFA contains data from AraCyc (version 9.0) for A. thaliana and RiceCyc (version 3.2) for O. sativa.

• MapMan (http://mapman.gabipd.org) (40), which is a user-driven tool for large datasets (e.g. gene expression data from microarrays) visualized in the context of diagrams of metabolic pathways or other processes. MapMan annotation data (bin tree version 1.1) for both A. thaliana and O. sativa are stored in dbWFA. The dbWFA database also provides a function to automatically generate MapMan T. aestivum mapping files.

• Munich Information Center for Protein Sequences Functional Catalogue (MIPS FunCat; http://www.helmholtz-muenchen.de/en/mips/projects/funcat) (41), which provides a hierarchical scheme for the functional description of proteins of prokaryotic and eukaryotic origin. MIPS FunCat annotations for A. thaliana (MAtDB version 2.1) are stored in dbWFA.

• A. thaliana Gene Family Information (TAIR version 10; http://www.arabidopsis.org/browse/genefamily) (42), which provides gene family information for the plant model species A. thaliana.

The 17 541 full-length cDNA sequences from the TriFLDB and other public databases and all the transcript sequences from the T. aestivum UniGene set (builds #55, #58, #59 and #60) were processed using the BLASTx algorithm against A. thaliana and O. sativa predicted cDNA sequences (Figure 1). Build #55 (the one used to develop the T. aestivum NimbleGen 40 k microarray) (37) and the following major releases were retained, as users may have developed resources based on different builds of the UniGene even though NCBI only stores the most recent build. BLAST results with an e-value >10⁻³ were not stored in the database, as we considered this would be too poor a match for most research. No other filter was applied to the BLAST results before their insertion into the database. All the parameters from the BLAST tabular results were kept, and >30 × 10⁶ BLAST results for the UniGenes and 95 × 10⁶ BLAST results for the ESTs shaping the UniGene clusters were stored, so they could be rapidly screened when querying the database.

The database also contains curated information on T. aestivum transcription factors (2891 transcripts), E3 ubiquitin ligases of the ubiquitin-proteasome system (876 transcripts), hormone-responsive genes (467 transcripts) and seed storage proteins (55 transcripts; Figure 1). Transcription factor UniGenes were retrieved from the wDBTF database (34). E3 ligase and hormone-responsive UniGenes were recovered from the NCBI and TAIR databases using all A. thaliana and O. sativa E3 ligase and hormone-responsive sequences as the query in homology searches using the BLASTn, BLASTx and tBLASTx programs (36). The BLAST hits were filtered using an e-value threshold of 10⁻⁵ and an alignment length exceeding 80 bp. All sequences were checked for consistency and for the presence of specific protein signatures using the InterProScan program (http://www.ebi.ac.uk/Tools/pfa/iprscan/). For seed storage protein UniGenes, homology searches were performed on the whole UniGene build #55, using BLASTx and T. aestivum seed storage protein sequences as reference. No preliminary filter was applied to BLASTx results. Instead, all the alignments were carefully examined, and similarity in known conserved critical regions of seed storage proteins was given priority over e-value and BLAST score alone. In dbWFA, curated UniGene annotations are assigned to T. aestivum transcripts without any intermediate BLAST result.

Following the recommendations of the International Wheat Genome Sequencing Consortium (IWGSC) for annotating T. aestivum genomic sequences (3), the percentages of coverage (with respect to the length of the orthologous proteins) and identity are used to assign functional annotations to a transcript. In dbWFA, users can define the value of these two parameters, but we strongly recommend using the cutoff values suggested by the IWGSC, where BLAST results with an identity >45% and coverage >50% are assigned a ‘putative function’ and BLAST results with identity and coverage >90% are assigned a ‘known function’.

All data are stored in a MySQL database. The integration of the database allows one to assign the functional annotation from any of the systems described above to the transcripts of interest and vice versa. The dbWFA database thus provides a very powerful resource for the annotation of T. aestivum UniGenes. To find the most commonly sought types of information from dbWFA, simple yet pertinent queries with their parameters can be sent through a web-based interface (Figure 2). The results are delivered as html pages, and an export procedure is available to retrieve data in spreadsheet. The html result pages provide links redirecting the user to websites of the different annotation systems, allowing a global analysis of the annotation results. The web interface can also be used to automatically create MapMan mapping files for the search results. Although the dbWFA web interface only allows data mining of common queries, specific queries can be performed using the SQL database, which can be downloaded from the dbWFA website. The modularity of the database will facilitate the integration of new T. aestivum data as transcripts are sequenced and annotated through different pipelines.

Using dbWFA: Percentage of Annotated UniGenes, Comparison of T. aestivum UniGene and A. thaliana and O. sativa Whole-Genome Annotation and Query Examples

Thirty-four percent (13 713 transcript sequences), 40% (14 843), 35% (20 016) and 35% (20 034) of the transcript sequences of the UniGene builds #55, #58, #59 and #60, respectively, have a putative functional annotation in at least one of the annotation resources. Eighty-one percent of the 17 541 full-length cDNA sequences from TriFLDB have a putative functional annotation in at least one of the annotation resources. The number of transcripts and full-length cDNA sequences annotated in the different resources are given in Table 1. BLASTn analysis revealed that 12 478 full-length cDNA sequences matched a sequence in the UniGene set (build #60) with a coverage and identity threshold value >50 and 90%, respectively. Among these 12 478 correspondences, 10 996 and 5 932 full-length cDNA sequences and UniGene sequences, respectively, have a putative functional annotation in at least one of the annotation resources. This result highlights the additional information brought by the full-length cDNA sequences.

The quality of the annotation method is illustrated by comparing the MapMan bin annotation of all the transcripts of the T. aestivum NimbleGen 40 k microarray (developed with UniGene build #55) and the full-length cDNA sequences from TriFLDB with the annotation of A. thaliana and O. sativa imported from MapMan and recorded in the database. The MapMan bins were used here because this annotation system is available for the three species. Overall, there was no clear bias between the three species (Figure 3), and the percentages of genes in the 26 categories for the three species were well correlated (T. aestivum versus A. thaliana: r = 0.96, P < 0.001; T. aestivum versus O. sativa: r = 0.69, P < 0.001), with no significant bias (P < 0.001). The higher correlation found with A. thaliana compared with O. sativa is mainly because there are fewer annotated transcripts in the DNA bin for O. sativa than for A. thaliana and T. aestivum (r = 0.90 for T. aestivum versus O. sativa when this bin is not considered). For the full-length coding sequences retrieved from TriFLDB and other public databases, the correlations between T. aestivum and A. thaliana and between T. aestivum and O. sativa were the same (r = 0.90, P < 0.001). The pairwise correlations between the four MapMan bin annotations presented were remarkably high, all >0.9 when the DNA bin was omitted. Similar results were obtained for the PlantCyc pathway reactions and GO (data not shown).

Unlike many annotation tools, dbWFA makes it possible to query multiple annotation systems simultaneously. To demonstrate various features of the dbWFA database, some query examples are presented in Box 1, using either the website or the database installed on a local machine.

Identification and Annotation of UniGenes Specifically Expressed During Either the Early or Late Stage of Grain Development

A total of 39 029 transcripts from the UniGene set (build #55) and 1613 transcription factors from the wDBTF database not present in the UniGene set are spotted on the custom T. aestivum NimbleGen 40 k microarray (36). Previous studies have shown that 18 140 (44.6%) of these transcripts are expressed during T. aestivum grain development (47). In dbWFA, 34–40% (depending on the build) of these transcripts have a putative functional annotation.

T. aestivum grain development comprises several distinct phases, starting with a syncytial then a cellularization phase (ca. 0–100°Cdays after anthesis), followed by a first differentiation phase of active endosperm cell division (ECD), expansion and differentiation (ca. 100–250°Cdays after anthesis), a second differentiation phase when storage polymers rapidly accumulate (ca. 250–750°Cdays after anthesis) and a maturation phase when grain rapidly desiccates (ca. 750–900°Cdays after anthesis) (48, 49). The transitions between these phases are associated with major changes in the grain transcriptome (5, 36, 50, 51) and proteome (52, 53).

To validate the principle underpinning the database and provide another example of the usefulness of dbWFA, we analysed the functional annotation of the transcripts specifically expressed during either the ECD or SPA phase of grain development. We used transcriptome data obtained with the custom T. aestivum NimbleGen 40 k microarray for the T. aestivum cultivar Recital grown under standard conditions in a greenhouse and sampled every 34–117°Cdays between 132 and 686°Cdays after anthesis (35). Transcripts with different patterns of expression were classified with J-Express 2012 software package (54) using Euclidean distance–based k-means clustering. For this analysis, the number of clusters was empirically set at 25 because it allowed us to clearly discriminate between gene expression clusters specific to the ECD and SPA phases of grain development. One cluster of 238 genes contained genes expressed exclusively during ECD stages (Figure 4A and B). Two other clusters contained genes expressed exclusively during SPA stages. The latter two clusters were merged to form a single SPA cluster of 195 genes. dbWFA was then used to retrieve the functional classification of the transcripts from both clusters. The MIPS Functional Classification was used, as it was the most informative and straightforward for comparisons with previous studies.

Using IWGSC-recommended coverage (50%) and identity (45%) percentages, 68 (29%) ECD-specific transcripts and 129 (66%) SPA-specific transcripts were assigned to an MIPS functional category, respectively (Figure 4C). The annotation results were consistent with previous transcriptome (5, 51) and proteome (52, 53) studies of developing T. aestivum and H. vulgare (55) grain. The functional classifications of the gene clusters were different. Not surprisingly, 12 transcripts involved in cell fate and cell type, tissue differentiation and organ differentiation were specifically expressed during the ECD phase, whereas no SPA phase–specific transcripts were annotated as being in these MIPS functional categories. Also in good agreement with our knowledge of grain development, 55 seed storage protein transcripts were specifically expressed during the SPA phase, while none was found among the annotated ECD-specific genes.

Quantitative differences in the annotation of these two clusters of transcripts were also observed. Several transcripts in the SPA-specific cluster were involved in cell rescue, defence and virulence and in the interaction with the environment. In particular, transcripts involved in plant hormonal regulation were overrepresented in the SPA-specific gene cluster. Transcripts coding for proteins involved in protein synthesis and proteins with metabolic functions were overrepresented in the ECD cluster. These results coincide with previous transcriptome (5, 36, 56) and proteome analyses (53). Finally, we note a substantial difference in the MIPS functional category ‘protein with binding function or co-factor requirement’, with more transcripts involved in DNA binding in the ECD cluster and more transcripts involved in RNA binding in the SPA cluster.

All these data show great similarity to published results for T. aestivum and H. vulgare, reflecting the accuracy of the automatic annotation provided by dbWFA. Unlike several other T. aestivum transcriptome analyses where a complex process had to be carried out to assign a functional annotation to selected transcripts and/or proteins (5, 57), here only a single request was made to the dbWFA database to retrieve the functional classification of 45% of the transcripts of interest, and 40% of the 40 642 transcripts of the T. aestivum NimbleGen 40 k microarray. This percentage of annotated transcripts is similar to that previously reported (38%) for the T. aestivum Affymetrix GeneChip® microarray (4).

Outlook

The dbWFA database was created by integrating numerous data sources. As a result, it is a practical source of heterogeneous data for functional annotation of T. aestivum transcripts. The website grants access to the most common queries that can be applied to the database, and the freely available MySQL database is a powerful tool for more specific requests. Although further analyses are required to confirm the dbWFA annotation results, the database provides an efficient and fast solution for acquiring a wide range of functional information. cDNA resources are useful to predict exonic regions from genomic sequences; thus, efforts to annotate the UniGene resources will significantly contribute to the analysis of sequence data produced by ongoing T. aestivum initiatives and other genome sequencing projects.

The version of dbWFA presented here is operational, but the aim is not to restrict the database to storing O. sativa and A. thaliana annotations but to expand it to include data from other plant species genomes as their functional annotation becomes more consistent. Integration of InterProScan (58) in the workflow could be a valuable way of augmenting the process. Also, the integration of AFAWE (59) would provide an annotation workflow with different function prediction tools. However, the current version of AFAWE cannot be used independently from its web interface and would thus have to be implemented using the tools called in its workflow, which are available as web services. Finally, the upcoming integration of a BLAST program in the workflow will allow users to annotate their own sequences and will also make dbWFA applicable to other species.

Acknowledgements

The authors thank Dr Etienne Paux and Dr Catherine Feuillet (INRA, UMR1095 GDEC, Clermont-Ferrand, France) for useful discussions and advice, and Mr Sébastien Reboux, Ms Claire Viseux and Mr Michael Alaux (INRA, URGI, Versailles, France) for installing and maintaining the database on the URGI server.

Funding

This work was supported by a PhD grant from the French Ministry for Higher Education and Research to J.V.

Conflict of interest. None declared.

References

Author notes