RiceFREND: a platform for retrieving coexpressed gene networks in rice.

Similarity of gene expression across a wide range of biological conditions can be efficiently used in characterization of gene function. We have constructed a rice gene coexpression database, RiceFREND (http://ricefrend.dna.affrc.go.jp/), to identify gene modules with similar expression profiles and provide a platform for more accurate prediction of gene functions. Coexpression analysis of 27 201 genes was performed against 815 microarray data derived from expression profiling of various organs and tissues at different developmental stages, mature organs throughout the growth from transplanting until harvesting in the field and plant hormone treatment conditions, using a single microarray platform. The database is provided with two search options, namely, 'single guide gene search' and 'multiple guide gene search' to efficiently retrieve information on coexpressed genes. A user-friendly web interface facilitates visualization and interpretation of gene coexpression networks in HyperTree, Cytoscape Web and Graphviz formats. In addition, analysis tools for identification of enriched Gene Ontology terms and cis-elements provide clue for better prediction of biological functions associated with the coexpressed genes. These features allow users to clarify gene functions and gene regulatory networks that could lead to a more thorough understanding of many complex agronomic traits.

INTRODUCTION

Rice (Oryza sativa L.) has been at the forefront of plant genomics since the completion of the genome sequence of the japonica cultivar Nipponbare. The high-quality map-based sequence and complete annotation of every transcriptional unit in the genome have facilitated remarkable advancements in molecular genetics of rice as well as in many applied aspects of cereal genomics (1–4). Subsequently, elucidating the function of all predicted genes has become an ultimate goal not only for complete understanding of the biology of rice but also in developing novel strategies for crop improvement. Due to worldwide efforts in functional genomics, a wide range of resources for characterizing the function of rice genes have been developed and made available to the scientific community (5–9). However, despite these extensive efforts, more than half of about 32 000 genes that define the rice plant, many of which could be possibly involved in the expression of agronomically important traits, still remains to be fully characterized.

Gene clustering analysis, which is based on the premise that clusters of genes with similar expression patterns across several experimental conditions tend to be functionally related (10), has become a powerful tool for gene function prediction. With the rapid accumulation of microarray data in public repositories, such as NCBI’s Gene Expression Omnibus (GEO) (11,12), European Bioinformatics Institute (EBI)’s Array Express (13) and DNA Data Bank of Japan (DDBJ)’s CIBEX (14), there is also a proliferation of databases on gene expression profiling of various crops in the public domain (15–19). Similarly, an increasing number of gene coexpression databases in various plant species such as Arabidopsis and rice have also been developed (20–26). Coexpression approaches have been widely used in Arabidopsis, particularly in analysis of various biological targets such as transcription factors, enzymes in specific metabolic pathways and protein subunits (27). In rice, gene coexpression data are provided by the RiceArrayNet (25), OryzaExpress (26), ATTED-II (28) and Rice Oligonucleotide Array Database (19). RiceArrayNet was developed based on the Rice 60k Microarray data, and the rest were constructed using rice Affymetrix microarray data.

In the last 5 years, we have collected >800 microarray data from organs and tissues of rice plants at various developmental stages, mature organs throughout the entire growth in the field, cell/tissue types isolated by laser microdissection and rice plants subjected to plant hormone treatment (29,30). These expression profiles generated from the model rice cultivar Nipponbare using a single microarray system (Agilent) are available in our gene expression profile database, RiceXPro (17). Most of the data were collected from rice plants grown in the natural environment, which is much more complex than controlled laboratory conditions. This large-scale field transcriptome profiling can be used as a measure of the corresponding changes in physiological and morphological states of the plant triggered by internal or external stimuli under natural conditions. Therefore, our comprehensive collection of gene expression profiles provides a more accurate landscape of the transcriptome that could be used for understanding the function of genes based on the level, time specificity and position specificity of expression, as well as the overall similarity of gene expression of various organs and tissues in the natural field environment.

Here, we describe the Rice Functionally Related gene Expression Network Database (RiceFREND), a database for retrieving gene coexpression information across a wide range of gene expression profiles which we collected using a single microarray platform.

DATABASE COMPONENTS AND FEATURES

The overall concept of RiceFREND is to provide a platform for efficient prediction of gene functions and gene regulatory networks based on gene expression similarity. A schematic diagram of the major components and various features of the database such as search options, coexpression network viewers, analysis tools and detailed information on coexpressed genes is shown in Supplementary Figure S1.

Calculation for coexpression analysis

All gene expression data described in this database were obtained from microarray analysis using the rice 4×44K microarray RAP-DB platform (Agilent Technologies, G2519F#15241), which contains 35 760 independent probes corresponding to 27 201 annotated loci published in the Rice Annotation Project Database (RAP-DB) (4). A total of 24 data sets corresponding to 815 microarray data from various organs and tissues at different developmental stages and treatment with plant hormones, including redundant data among the data sets, were used as data source for coexpression analysis (Supplementary Table S1). All data have been deposited in GEO (11,12) and accessible through GEO series accession numbers: GSE21396, GSE21397, GSE36040, GSE36042, GSE36043, GSE36044, GSE39423, GSE39424, GSE39425, GSE39426, GSE39427, GSE39429 and GSE39432.

The expression data were obtained either by intensity-based one-color or ratio-based two-color microarray hybridization system. In the 75th percentile normalization for one-color method, the processed raw signal intensities of 35 760 probes were divided by 75th percentile value and transformed to log2 scale. The median expression value within each data set was further subtracted to obtain the relative value for each probe. Mean values were used for redundant probes assigned to the same locus. For data obtained using the two-color microarray system, the expression profile for each gene corresponds to the log2 ratio of cyanine 5 (Cy5) signal intensity for treatment to cyanine 3 (Cy3) signal intensity for mock treatment (log2 Cy5/Cy3). All data sets were combined into one gene expression matrix data for coexpression analysis among genes. The weighted Pearson’s correlation coefficients were calculated to reduce any unsuitable effects that may have resulted from sample redundancies. The mutual rank (MR) was used as an index for coexpression as described in ATTED-II (21,31).

Annotation of genes

The RAP2 representative data in general feature format derived from RAP-DB (4) were basically used as annotation of each gene. The corresponding Michigan State University Osa1 Rice Locus (MSU ID) (32) for each RAP locus ID was obtained using the ID converter tool in RAP-DB. Generic Gene Ontology (GO) annotations in RAP-DB were converted to GO slim using the map2slim.pl script (http://search.cpan.org/∼cmungall/go-perl/). A total of 1308 genes commonly described in PlnTFDB (33) and PlantTFDB (34) were classified as transcription factors. The subcellular localization sites of proteins were predicted by WoLF PSORT (35) based on the deduced amino acid sequences of each gene. The gene symbols/names of 1523 genes and functional categories of 2053 genes were obtained from Oryzabase (36) and KEGG Pathway Database (37), respectively.

Coexpression search with a single guide gene

The coexpression search is based on the ‘guide gene’ approach in which a pre-selected gene or set of genes is used to retrieve coexpressed genes. In the single guide gene search option, coexpression search is initiated by entering a keyword such as RAP locus ID, gene name, MSU ID and other identifiers or by selecting specific chromosome in the panel representing the 12 rice chromosomes. This provides a tabular list of coexpressed genes with information on RAP locus ID, gene description, transcription factor, gene symbol and the corresponding MSU ID. The RAP locus ID opens a pop-up window with direct links to RiceXPro (17) for gene expression profile used in coexpression analysis, RAP-DB (4) for detailed annotations and SALAD database (38) for genome-wide comparative analysis of motifs. The MSU ID provides direct link to the annotation of the gene as described in the MSU Rice Genome Annotation Project database (32). The coexpression data for each gene are shown in two formats, a tabular list of coexpressed genes in ascending order of MR (Figure 1A) and a graphical representation of the coexpression network in HyperTree (http://philogb.github.com/jit/) format (Figure 1B). The MR list also contains information on predicted subcellular localization sites of proteins, KEGG pathways and transcription factors. Coexpressed genes can be selected and used for further analysis, such as GO enrichment test and cis-element analysis as described below. A filter option allows users to retrieve information for gene/genes of interest using RAP locus ID, description, MSU ID, transcription factor and gene symbol/name as keywords.

Figure 1.

Single guide gene search results interfaces. (A) Tabular list of coexpressed genes generated from the search with Os03g0752800 (MADS14) as guide gene. The coexpressed genes are arranged in the order of decreasing MR value. (B) Graphical view of the coexpressed gene network in HyperTree format. Two nodes representing coexpressed genes with MR value below a defined threshold are joined to form an edge. The gene network (MR <7, hierarchy = 3) shown here contains 87 nodes and 114 edges. (C) An enlarged view of the coexpression network using gene symbol as label for nodes. MADS14 is strongly coexpressed with MADS15, MADS18 and MADS34.

The HyperTree graphical viewer shows the relationships of coexpressed genes based on a user-defined hierarchy and MR value or rank based on ascending MR value. In this viewer, the nodes representing the coexpressed genes with MR value below a defined threshold are joined to form edges (Figure 1B and C). The coexpression data can be viewed for a maximum of 100 coexpressed genes and the network size can be set by selecting the MR value and hierarchy. Detailed information on the nodes comprising a network is provided from the ‘Node List’ and via a pop-up window on each node in the HyperTree. The nodes are labeled with RAP locus ID but options for other identifiers such as gene symbol, gene name, transcription factor, MSU ID and KEGG pathway are available. The node for a transcription factor is indicated by a red square and the node with a KEGG pathway annotation is indicated by a tag corresponding to a specific pathway. All genes with KEGG pathway annotation are listed in a table with direct links to the database. In addition to HyperTree, the gene network can also be viewed in Cytoscape Web (39) and Graphviz (http://www.graphviz.org/) formats, thereby providing several options to visually understand the relationship among coexpressed genes. Download options are available for the node list, data for construction of network in Cytoscape (40) and the graph images in Cytoscape Web and Graphviz formats.

The coexpressed gene network generated using MADS14 (Os03g0752800), which encode APETALA1/FRUITFUL-like MADS domain protein, as a guide gene is shown in Figure 1B and C. MADS14 is coexpressed with other MADS-box genes, namely, MADS15 (Os07g0108900), MADS18 (Os07g0605200) and MADS34 (Os03g0753100). It has been reported recently that the four MADS-box genes were coordinately expressed in the shoot apical meristem to specify the identity of the inflorescence meristem during reproductive phase transition (41). Therefore, the gene network would reflect overlapping functions of the four MADS-box genes and provided useful clue that could be used to further elucidate the molecular mechanism of reproductive transition in rice.

Coexpression search with multiple guide genes

The multiple guide gene search uses two or more pre-selected genes as guide genes to expand the coverage of gene coexpression networks. This search is also initiated by entering a keyword such as RAP locus ID, gene name, MSU ID or other identifiers. From the tabular list generated by the search, two or more guide genes are selected, and the hierarchy and MR value or rank to be used for construction of gene network are defined. Information useful for selection of threshold in MR value among query genes can be viewed in the matrix data (Figure 2A). The coexpression network is generated in either a Cytoscape Web format with several options for node label and network layout (Figure 2B) or Graphviz format, both of which provide detailed annotation of nodes, direct links to other databases and the HyperTree and MR List from the single guide gene search against each node in the network. Additionally, the coexpression network can be reconfigured by selecting other guide genes from the node list panel. As in the single guide gene search, download options are also available for note list, data for construction of Cytoscape and the graph images in Cytoscape Web and Graphviz formats.

Figure 2.

Multiple guide gene search results interfaces. (A) Matrix data of MR values among multiple guide genes. Two or more guide genes can be selected to initiate coexpression search. (B) Graphical view of the coexpressed gene network in Cytoscape Web format. The guide genes are shown as large circles with yellow outline. (C) An enlarged view of the coexpression network using gene symbol as label for nodes. Genes involved in momilactone biosynthesis such as CPS4, KS4, MAS, CYP99A2 and CYP99A3 are strongly coexpressed.

The gene network generated by the multiple guide gene search (MR < 7, hierarchy = 2) with ent-kaurene synthase gene Os04g0611800 (OsKS1) and six kaurene synthase-like genes, namely, Os04g0179700 (OsKSL4), Os02g0571300 (OsKSL5), Os02g0571800 (OsKSL6), Os02g0570400 (OsKSL7), Os11g0474800 (OsKSL8) and Os12g0491800 (OsKSL10) as guide genes is shown in Figure 2B. OsKS1 is associated with gibberellin biosynthesis (42), whereas OsKSL4, OsKSL8, OsKSL10 and OsKSL7 are related to biosynthesis of diterpenoid phytoalexin such as momilactones A and B, oryzalexins S, oryzalexins A–F and phytocassanes A–E, respectively (43–45). This multiple guide gene search generated a gene network comprising four clusters. The largest cluster contains OsKSL4, OsKSL7, OsKSL8 and OsKSL10 which have been reported to be induced by UV irradiation and elicitor treatment (42–45). On the other hand, OsKS1, and OsKSL5 and OsKSL6, the expression of which is not affected by such treatment (46), were included in a separate cluster. Interestingly, momilactone biosynthetic genes (47), namely, CPS4, KSL4, MAS, CYP99A2 and CYP99A3 are strongly connected to each other (Figure 2C), and several phytocassane biosynthetic genes (48) are also included in the largest cluster. This coexpression analysis suggests that the resulting gene network can be used to define the biological function of each kaurene synthase genes and may therefore help in exploring novel enzymes related to phytoalexine biosynthesis based on the coexpression signature.

Analysis tools

In order to enhance prediction of biological functions associated with the coexpressed gene network, RiceFREND is provided with analysis tools for identification of enriched GO terms and cis-elements. The GO enrichment test facilitates the identification of overrepresented GO terms among coexpressed genes for better understanding of the functional relationships of these genes. The Fisher’s exact test is applied against GO terms categorized into biological process and molecular function. This generates a tabular list of overrepresented GO terms in ascending order of P-value which is adjusted by the Benjamini and Hochberg’s method (49) for multiple testing collection (Figure 3A). The result can also be viewed in graph format with a filter option (Figure 3B) that can be used to construct a graph against GO terms with case-positive value >5.

Figure 3.

GO enrichment test results interfaces. (A) A tabular list of over-represented GO terms (false discovery rate, FDR < 0.05) generated by GO enrichment test against the top 100 genes coexpressed with LEA3-1. (B) The over-represented terms can be visualized in a graphical format with the enriched GO terms (FDR < 0.05) shown in color scale corresponding to the FDR value. The graph was constructed with filter off setting.

Identification of cis-elements among coexpressed genes could provide additional information for prediction of function. The cis-element analysis tool can be used to clarify whether a specific motif is overrepresented in the 1-kb upstream regions of coexpressed genes based on three statistical tests, namely, χ2 test, Fisher’s exact test and G-test.

We applied the GO enrichment test and cis-element analysis on the top 100 genes coexpressed with LEA3-1 gene (Os05g0542500). It has been reported that over-expression of LEA3-1 improves drought resistance under the natural field conditions (50). Consistent with this observation, the GO terms such as ‘response to water’, ‘response to stress’ and ‘response to abiotic stimulus’ were found to be overrepresented among the coexpressed genes (Figure 3A and B). The expression of LEA3-1 gene is induced not only by drought stress but also by treatment with abscisic acid (ABA). In the cis-element analysis, ABA-responsive element core (ACGT) containing CGMCACGTGB motif (51) is also statistically enriched in the promoter region of the coexpressed genes. Therefore, although the microarray data used for coexpression analysis do not include abiotic stress treatments such as drought and salinity, the results of coexpression search and additional analysis can be used to understand functional relationships among stress-inducible genes based on the coexpression.

FUTURE DIRECTIONS

Coexpression analysis can be used more accurately with a more balanced distribution of expression data from various organs and tissues and a wider coverage in terms of experimental conditions. We are currently undertaking microarray analysis involving various organs and tissues at different developmental stages and rice plants subjected to biotic and abiotic stresses. These data will be added to the database and the annotation information will also be updated to further enhance the functionalities and to provide a more reliable platform for retrieving gene expression networks in rice.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online: Supplementary Table 1 and Supplementary Figure 1.

FUNDING

The Ministry of Agriculture, Forestry and Fisheries (MAFF) of Japan; Genomics for Agricultural Innovation [RTR0002 to Y.N.]. Funding for open access charge: Genomics for Agricultural Innovation [RTR0002].

Conflict of interest statement. None declared.