WheatExp: an RNA-seq expression database for polyploid wheat

Data sources and generation

We designed a pipeline specific for polyploid wheat sequence data to analyze previously published RNA-seq datasets using a uniform set of tools and quality controls. The output of our pipeline is a set of expression values for all annotated wheat genes from the IWGSC project. Briefly, raw RNA-seq reads are first trimmed for quality and adapter contamination using two open-source packages, “Sickle” (https://github.com/ucdavis-bioinformatics/sickle) and “Scythe” (https://github.com/vsbuffalo/scythe), respectively, ensuring that only high-quality reads are considered when generating expression profiles. Trimmed reads are mapped to the full set of annotated homoeologue-specific wheat transcripts from the Ensembl genomics platform using BWA [17]. Uniquely-mapped reads are counted using “Htseq-count” [18] and then adjusted to derive RPKM/FPKM (Reads/Fragments per kilobase of transcript per million mapped reads) values for each gene based upon mapping rate, transcript length and library size. This normalization means that within a dataset, expression values are directly comparable across different tissues and developmental time points. Although the same reference was used for each dataset, comparisons across different datasets are less reliable because of differences in the number and length of sequencing reads between different datasets. Our mapping parameters are selected to report only those reads with a mapping quality (MAPQ) score of 40 from the Sequence Alignment/Map (SAM) file, a value which signifies that the read was mapped uniquely. Reads which map ambiguously, either to multiple homoeologues or to other identical sequences, have a lower associated MAPQ score, and are excluded in this step. Table 1 reports the % of reads mapped from each dataset after the application of this selection criterion. Across all six datasets, an average of 50.1 % of reads were mapped uniquely, resulting in homoeologue-specific expression data for each gene. In general, datasets with longer reads (e.g. 101 bp PE reads) resulted in a higher proportion of uniquely mapped reads than those comprised of shorter reads (e.g. 50 bp SE reads).

Web implementation

The web interface was constructed using several different programming packages. The code base for the majority of the project is PHP (https://secure.php.net/) and JavaScript (https://www.javascript.com/). Relational database queries to the backend are performed with the PHP Data Object (PDO) module, allowing for secure queries. An additional advantage of using the PDO module is that the code is compatible with standard database engines such as MySQL, PostGreSQL and SQLite. In order to display dynamic graphs of the data, we implemented the HighCharts JavaScript library (https://github.com/highslide-software/highcharts.com). Specifically, this project uses a PHP module, which implements the HighCharts JavaScript library freely available on github (https://github.com/ghunti/HighchartsPHP). For dynamic text searches in portions of the website, the project implements Asynchronous JavaScript and XML (AJAX) technology using the package JQuery 1.11.3 (https://jquery.com/). Custom PHP and JavaScript code was written to develop a frontend website to enable BLAST [19] searches and to select multiple results for expression display. The site’s frontend was written in HTML and JavaScript with BLAST search [19] and AJAX gene identifier search forms allowing the user to select multiple results for expression display.

Database implementation

The database implementation uses a flexible storage schema to house the data. The storage table has the following MySQL (https://www.mysql.com/) storage datatypes: study_id (varchar), seq_id (varchar), tissue (varchar), mean (float), se (float), se_low (float) and se_high (float). Binary search tree indices (BTREE) were implemented to increase the speed of queries using the study_id and seq_id columns.

System architecture

The WheatExp tool is housed on the GrainGenes server at the following URL http://wheat.pw.usda.gov/WheatExp/. GrainGenes is an internationally recognized database for genomic and genetic resources in Triticeae and Avena species. The web and database is running on a ThinkMate ThinkTank IQ4 system with four Intel Xeon E7-4820 at 2.00 GHz with 52 GB of RAM. It is currently running Linux Kernel 3.13 distribution Ubuntu 14.04 long term support with PHP version 5.5.9 and MySQL version 5.5.43.

Web interface

The WheatExp homepage includes a brief description of the database and project design as well as details of all currently available datasets (Fig. 1a). Our data processing pipeline allows for the rapid incorporation of complementary RNA-seq expression datasets as they are published and we invite suggestions for the addition of new datasets from the user community. We anticipate regular expansions of the database to broaden the range of developmental and temporal expression profiles included. This approach will maximize the utility of the database for researchers studying diverse aspects of wheat development and ensures access to the most relevant high-quality expression datasets.

From this main hub (Fig. 1a), the database can be queried in one of two ways; either by entering the DNA or protein sequence of a gene of interest as a BLAST query, or by a text search for a known gene ID from the Ensembl genomics annotation platform [5] (e.g. Traes_6AS_9E38A95CB.1) or for an annotated functional term associated with the gene’s encoded protein (e.g. “bHLH” or “Cytochrome P450”). For BLAST searches, results are displayed on a new page and include details of each BLAST alignment, sequence and a link to the corresponding gene ID page on the external Ensembl genomics hub for simple cross-referencing (Fig. 1b). A maximum of six matched results may be selected for side-by-side display within the same graph to allow simple comparisons between multiple genes. While this feature was originally implemented to enable comparisons among wheat homoeologues, any set of up to six genes may be selected for comparison, regardless of their relationship.

Likewise, when browsing using the text search function, up to six genes can be selected for addition to the results list, which can subsequently be viewed side-by-side in the results window. For larger-scale analyses, tabular expression data for any number of genes can be downloaded by providing a list of Ensembl gene IDs of interest. The functional terms associated with each gene are obtained through standard gene annotation files in GFF3 format from the IWGSC which are stored within the database for text search function. We chose to adhere to the widely-used standard gene nomenclature format employed by the IWGSC and Ensembl genomics platform [5] and selected the set of annotated cDNA sequences from this platform as our mapping reference. External links to the annotated sequences for each gene are included in the results. This nomenclature format is increasingly becoming the standard for gene annotations within the plant research community, so our use of this reference will allow for the simple translation between projects and will maintain complementarity with the IWGSC project. This will facilitate comparative genomics studies with model plant species and other economically-important crops, such as rice, barley and maize, as the genomic resources contained within the Ensembl platform in each of these species improves. Additionally, comparisons can be made with more distantly related species to analyze functional gene divergence during the course of evolution.

Graphical expression profiles from all datasets are presented on a single results page, displaying mean RPKM/FPKM values +/− Standard Error Mean (SEM) (Fig. 1c). Graphs can be downloaded in one of four image formats and data is also presented in an accompanying table, which can be exported in ‘.csv’ format (Fig. 1c). Gene-level expression data can be downloaded separately, or in bulk as a single tabular file containing all data.

Expression data

All expression profiles in WheatExp are generated from RNA-seq datasets. This approach has several advantages over existing expression studies derived from microarray data, which until recently, was the standard technology used for large-scale expression analysis (e.g. “Plant Expression database (PLEXdb)”, a database of microarray-based expression profiles in different plant species [20]). One of the advantages of RNA-seq is that it is an open platform that does not rely on predetermined sets of probes printed on a gene chip. In addition, this technology provides more reliable expression profiling across a broader dynamic range than is possible with microarrays.

An important advantage of the application of RNA-seq data in polyploid species is that it facilitates the distinction among homoeologues and recently-diverged paralogous genes by allowing the application of stringent read mapping thresholds. Our selection of only uniquely-mapped reads has the dual benefit that the expression data are not only robust, but also homoeologue-specific, since the differences between these genomes (average 97 % identical) are distinguished by the selected mapping parameters. This is illustrated in Fig. 2 by two examples: CIRCADIAN CLOCK ASSOCIATED1, where the expression of the three homoeologous genes is approximately equal (Fig. 2a) and CONSTANS1 where the D-genome homoeologue contributes the majority of transcripts to the overall expression (Fig. 2b).

Simulated RNA-seq data

One drawback of using uniquely-mapped RNA-seq reads for expression analysis is that any read which maps equally well to identical regions in different genes is discarded, potentially resulting in an underestimation of the expression levels of highly similar genes [21]. To determine the extent of this effect in our database, we performed a simulated RNA-seq experiment. We generated 29.4 M synthetic 100bp paired-end reads with random expression levels and Illumina HiSeq2000 error profiles (‘ART’, mode art-illumina, default parameters except -m 500, −s 100 –ss HS20 [22]). All reads were processed using the same pipeline as for all biological RNA-seq data. By comparing the known number of simulated reads with the number of mapped reads, we can determine for each contig the proportion of reads discarded during mapping. Using a set of 3,476 homoeologous triplets (=10,428 genes) identified from a previous study [7], we mapped the subset of reads originating from each homoeologue to a reference comprised only of their genome of origin (i.e. A-genome reads were mapped to A-genome transcripts etc.). For the A, B and D genomes, an average of 98.6, 98.4, and 98.4 % of reads mapped uniquely to their transcripts of origin, respectively, demonstrating that only a small proportion of reads are discarded during mapping when their homoeologous genes are absent from the reference. When we repeated the mapping of all generated reads to the full reference, unique mapping rates were reduced to 82.4, 83.6 and 80.6 % for the A, B and D homoeologous triplets. In each case, this was a slightly lower unique mapping rate than for all remaining transcripts in our dataset (84.4 %). Despite this reduction in the mapping rate, we observed a high level of correlation between the number of generated reads and the observed mapped reads (r = 0.95, 0.96, 0.95 for A, B and D homoeologous triplets, Fig. 3). Therefore, while the estimated expression levels of homoeologous genes in our database are, on average, slightly reduced due to their sequence similarity, the reported expression remains closely correlated with the true expression level. Furthermore, this effect is approximately equal for transcripts originating from the three homoeologous wheat genomes (Fig. 3), demonstrating the absence of bias when comparing homoeologue-specific expression profiles for a gene of interest.

Limitations

The main application of WheatExp is to compare the relative expression levels of the different homoeologues of a single gene across different tissues, developmental stages, environmental conditions and genetic backgrounds. For users interested in comparing the expression of different genes, we have included a statement on the website indicating that comparisons among genes are valid only when the genes being compared have the same number of homoeologues in the reference genome. Based upon results from our simulated RNA-seq experiment, genes where one homoeologue is absent from the reference will exhibit a higher proportion of uniquely-mapped reads and the expression levels of the two remaining homoeologues may also be inflated by the incorrect mapping of reads from the absent homoeologue. Additionally, no expression data will be reported for any genes which lack annotation within the current IWGSC release and any contig assemblies which are duplicated in the reference assembly will exhibit a reduced number of uniquely mapped reads. However, our project design allows for regular updates and refining of the mapping reference as this is expanded through the IWGSC project. As the mapping reference is improved we will re-map and re-process each dataset to generate updated expression sets using new versions of the reference, reducing the incidence and impact of such bias.

Our approach and data analysis pipeline can be applied to other polyploid species for which a homoeologue-specific genomic assembly is available to use as a reference. A critical parameter that must be considered in this application is the average level of identity among homoeologues, since this will affect the selection of the threshold for mapping uniquely mapped reads and thus the ability to discriminate between homoeologues.