Identification and Analysis of Rice Yield-Related Candidate Genes by Walking on the Functional Network

Rice (Oryza sativa L.) is one of the most important staple foods in the world. It is possible to identify candidate genes associated with rice yield using the model of random walk with restart on a functional similarity network. We demonstrated the high performance of this approach by a five-fold cross-validation experiment, as well as the robustness of the parameter r. We also assessed the strength of associations between known seeds and candidate genes in the light of the results scores. The candidates ranking at the top of the results list were considered to be the most relevant rice yield-related genes. This study provides a valuable alternative for rice breeding and biology research. The relevant dataset and script can be downloaded at the website: http://lab.malab.cn/jj/rice.htm.


INTRODUCTION
Rice (Oryza sativa L.) is one of the most important food crops worldwide, being used as the main food source by more than half of the global population (Mahender et al., 2016;Li et al., 2017). In the developing world, rice provides 27% of dietary energy and 20% of dietary protein (Huang et al., 2013). However, despite genetic improvements in grain yield delivered by the exploitation of semi-dwarfism and heterosis over the past 50 years, a substantial increase in grain productivity of the major crops is still required to feed a growing world population (Abe et al., 2018). The prime breeding target is to increase both grain size and grain number, because they impact both on yield potential and its end-use quality (Okada et al., 2018). However, the simultaneous improvement of grain quality and grain yield is a major challenge because of the well-established negative correlation between these two traits which is controlled by quantitative trait loci and influenced by environmental changes. Additionally, determining which genes in quantitative trait loci regulate grain size and number has not been clarified (Borzee et al., 2018;Li et al., 2018). Therefore, the identification genetic variants associated with improvements in grain yield would facilitate the breeding of new high-yielding rice varieties and may also be applicable to other crops (You et al., 2017).
Vast numbers of genetic variants have been detected by traditional genome-wide association studies and recent sequencing studies, and connecting the functional implications of these results to known genes has become a standard task (Li et al., 2015;Dehury et al., 2017;Torres and Henry, 2018;Wu et al., 2018). We previously developed a database, RicyerDB, to collect all known rice yieldrelated genes by integrating multiple omics data, information from the literature, and associated databases (Jiang et al., 2018). This work also established a search tool to query a particular gene, and to provide insights into gene functions and locations. Any rice yield-related gene can therefore be easily queried and the findings downloaded through the webpage, while candidate genes can be screened and prioritized to identify those most likely to be associated with known genes.
To achieve this goal, several approaches have been proposed from the perspective of computational systems biology (Behroozi-Khazaei and Nasirahmadi, 2017;He et al., 2017;Liu E. et al., 2017;Liu Y. et al., 2017;Xiong et al., 2017;Maione and Barbosa, 2018;Zhou et al., 2018). For example, the Endeavor tool uses the guilt-by-association principle to rank candidate genes according to their functional similarities to a set of predefined seed genes (Aerts et al., 2006;Tranchevent et al., 2008Tranchevent et al., , 2016. In recent years, a protein-protein interaction (PPI) network has been developed to achieve a global inference of entire genes (Liu et al., 2010;Lee, 2011;Rezadoost et al., 2016;Wang et al., 2016;Zeng et al., 2016;Luo and Liu, 2017;Holland and Johnson, 2018;Vlaic et al., 2018). PPI networks have also been used to provide a simplified yet systematic measure of functional similarities between genes (Chen et al., 2017a(Chen et al., , 2018a. Some methods for identifying yield-related genes have linked profile and sequence technology to facilitate the prediction of related genes. For example, Odilbekov et al. (2018) used machine learning and integrated this analysis with data obtained from spectroradiometer, infrared thermometer, and chlorophyll fluorescence measurements to identify the most predictive proxy measurements for studying Septoria tritici blotch disease of wheat.
Hybrid breeding is an effective tool to improve yield in rice, although parental selection remains a difficult issue. Xu et al. (2018) compared six genomic selection methods, such as least absolute shrinkage and selection operation and support vector machine, to evaluate predictabilities for different methods, and demonstrated their implementation to predict the hybrid performance of rice. Although good results have been achieved by these studies, the techniques of microarray and sequencing are nevertheless expensive.
The main target of this research was to use current knowledge to identify rice yield-related genes with network prediction methods. We proposed a computational systems biology approach for the identification of candidate genes via a random walk model on a PPI network with functional similarities (Kohler et al., 2008). Starting from known nodes, our method simulates the process in which a random walker travels to its neighbors or jumps to itself in the network, scores a gene using the probability that the walker stays in the gene at a steady state, and then ranks candidate genes according to their scores. Using a series of cross-validation experiments, we systematically demonstrated the robustness of our method, and applied our approach to predict a landscape of associations between known genes and candidates.

Flowchart Overview
We modeled the problem of identifying candidate genes associated with a set of known genes as a prioritization problem, and proposed to solve this problem using a three-step approach. As shown in Figure 1, taking the set of known genes as input, we first standardized the genes between STRING (Szklarczyk et al., 2015) and RicyerDB (Jiang et al., 2018). Then, we constructed a protein-protein network that scores the edges through functional similarities. This procedure applied a RWR algorithm to the network to calculate a score for each candidate gene, and then ranked the candidates to obtain a ranking list as the output (Chen et al., 2012a,b;Chen, 2016;Chen X. et al., 2016;Li et al., 2016;Peng et al., 2016;Zhu et al., 2018). Finally, the top candidate gene was verified according to its function and by the published literature.

Construction of the Functional Similarity Network
The functional similarity network is described as a graph G = (V, E), where V represents the nodes of the network and E stands for the edges of the network. The background network comes from the STRING database because of existing potential associated interactions among the proteins. The known rice yield-related genes were identified from our previous work with RicyerDB (Jiang et al., 2018). To standardize gene names between STRING and RicyerDB, genes were retrieved by reference to National Center for Biotechnology Information gene names. Functional similarities among genes in the background network were considered by scoring E for GO annotations. Using the latest release of the GO database (Ashburner et al., 2000;Chen L. et al., 2016;Raza, 2016; The Gene Ontology, 2017), edges were scored for a shared functional significance score of genes in the network that were annotated with GO terms.
The shared functional significance score F(i,j) between gene i and j was measured by the Weighted Shared Functions approach, which considered a gene's functions as a set of functional categories in GO. The functions shared by a small number of genes are taken to be far more significant than ones shared by a large number of genes. Each function had its own significance, which was defined as the inverse number of genes sharing the function. When two genes, i and j, have m functions in common, i.e., F(i)∩F(j) = ( f 1 , f 2 , . . ., f m ), F(i,j) was given as the total sum of the significance of the functions shared between them as follows: Here sig(f n ) denotes the significance of a function f n (n = 1,2,..., m) shared between genes i and j, | Genes (f n )| is the number of genes sharing a function f n . We calculated the ranking score, p, for each gene in the disease-related network and ranked these genes in the descending order of p.
FIGURE 1 | Illustration of the proposed method. Our method takes a set of seed genes as the input, and gives a ranking list of the candidates as the output. A functional similarity network was constructed by applying a random walk with restart algorithm to the network to obtain scores for candidate genes, and then the candidates were ranked according to their scores.

Random Walking on the Functional Similarity Network
We achieved the goal of identifying candidates related to known seeds by calculating a score for each candidate and then ranking the candidates to obtain a ranking list. The higher the rank, the more likely the gene was to be related to the given source nodes.
For this purpose, we adapted the RWR method in the functional similarity network. At the beginning, the walker chooses the seeds as the starting point. In each step of the walking process, the walker may start on a new journey with probability r or move on with probability 1−r. When moving on, the walker may move at random to one of its direct neighbors.
In our application, the initial probability vector P 0 was constructed such that equal probabilities were assigned to the nodes representing members of the disease, with the sum of the probabilities equal to 1. This is equivalent to letting the random walker begin from each of the known disease genes with equal probability. The transition matrix W is the column-normalized adjacency matrix of the graph, and P t is a vector in which the ith element holds the probability of being at node i at time step t. Formally, the RWR is defined as: Candidate genes were ranked according to the values in the steady-state probability vector P. P vector changes with time t, while it is possible to obtain it by explicitly calculating Equation (1) until convergence. The iteration is finished when the change between P t and P t+1 falls below 10 −10 . In this paper, we set default values for parameters r = 0.3 (see Results section for details).

Validation Method
We adopted a five-fold cross-validation experiment to assess the capability of RWR to identify the left seeds. All seed genes were divided equally into five parts, then one part was removed as a test set, and added to the candidate genes. All candidate genes were ranked by RWR to determine the ranking of the test gene. This procedure was repeated until all seed genes were used up as test genes.
In the context of the functional similarity network, the above validation procedure was equivalent to removing one part of the seed genes to candidate genes and determining whether candidates containing these seeds could receive a high rank. The r parameter of RWR ranged from [0,1] and was used to identify the ranking of the five parts. ROC curves were plotted, and areas under the ROC curve (AUC) values were used to evaluate the performance of r.

Data Sources
We obtained the rice background protein-protein network from the STRING database. In the network, protein associations were either directly derived from physical interactions or functional links from experimental evidence and computational methods (Jensen et al., 2009). The network composes of 6561 nodes and 567034 edges, which represent proteins and interactions between them, respectively. In our study, 136 known genes were selected as seed genes and other genes as candidate genes. We downloaded O. sativa Japonica protein network data through STRING version 10.5 (Szklarczyk et al., 2015).
Proteins with accurate functional annotations are vital to biological research. We obtained functional annotation information from the GO Consortium (Ashburner et al., 2000), and downloaded GO annotations of O. sativa from the most recent GO version. GO enrichment analysis is used to interpret high-throughput molecular data. GO annotation is the list of all annotated genes linked to ontological terms describing those genes.
The RicyerDB database integrates publicly available resources to construct a public platform for browsing and the interactive visualization of yield-related genes. The first release of RicyerDB contained more than 400 manually curated gene information entries which were all associated with rice yield.

Performance of the Proposed Method
The score vector P (the probability of being at the current node) for all genes in the network was calculated based on the ranking of corresponding r coefficients. Candidate genes were then ranked in the descending order of P score.
For optimal parameters, genes were also ranked according to the calculated p scores with nine different r-values (r = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 0.9). The matching numbers of the five-part seed genes were applied to assess the effectiveness of RWR. In Figure 2 listed the five cases of all, the number of matched seeds among the top 500 (every 100 is a measurement cutoff) in the ranking list of r = 0.3 was higher than other r-values in most cases.
The sum of the numbers of matched seed nodes in all ranking results was determined, and r = 0.3 was shown to have the maximum match in general. Finally, the parameter r = 0.3 was selected to calculate vector P to obtain the ranking results.  Further to detect the robustness of parameter r, we repeated the five-fold cross validation 100 times. Then we applying statistical analysis to compare the ranking of all seeds at different r-values in our model, the results were shown as Figure 3.

Prioritization of Candidate Genes and Validation by Literature Review
In the functional similarity network, all candidate genes were prioritized by RWR according to vector P at the final status. We manually searched the 100 top candidate genes ( Table 1) in PubMed 1 for their association with yield. This verified eight candidate genes associated with rice production. The LOC_Os11g40150 (rank 39) alias is OsRad51A1, which is a key component of homologous recombination in DNA repair. Direct interaction with OsNAC14 recruits factors involved in DNA damage repair and defense response, resulting in an improved tolerance to drought (Shim et al., 2018). LOC_Os04g37619 (rank 11) named ZEP, which is one of the key genes that involved hormone abscisic acid biosynthesis in rice by ion beam. Irritation can enhance the expression of genes involved in ABA biosynthesis, resulting in increasing content of endogenous plant hormone abscisic acid in rice (Chen et al., 2014).
Taken together, of the top 100 candidate genes in the ranking list, 46 candidate genes predicted by our method had been confirmed to be correlated with rice yield in PubMed literature ( Table 1). Top-ranked candidates were found to have a high confirmation rate in terms of their association with rice yield, especially top 20 candidates ( Table 2).
We conducted GO analysis to assess the functional enrichment of the top 100 candidate genes (Figure 4). The GO term having the most candidates annotated to was GO: 0005524 ∼ ATP binding, which is a binding motif within the primary structure of an ATP binding protein. A recently identified rice ATP binding cassette plays multiple roles in plant 1 http://www.ncbi.nlm.nih.gov/pubmed growth, development and environmental stress responses (Zhang X.D. et al., 2018). ATP binding has also been shown to play an important role in rice development (Coneva et al., 2014;Zhao et al., 2015;Chang et al., 2016;Lei et al., 2018).

DISCUSSION
In the present study, we identified genes associated with rice yield using the RWR method on a functional similarity network. We demonstrated the high performance of the RWR approach via a five-fold cross-validation experiment and showed the robustness of the parameter r. As an application of the RWR approach, we predicted a landscape of associations between known seeds and candidate genes. Our work has the following advantages. First, the RWR method can predict associations among known seed genes and candidate genes with the ability to spread the information that known seeds carried via their neighbors. Second, the interaction network provides a systematic view of functional similarities between genes by calculating GO terms. Finally, the robustness of the parameter r leads to a high level of accuracy in making predictions, and the method that achieving parameter can be adapted to other dataset.
Rice is the most important food crop worldwide. Use of the RWR method in the function similarity network can identify candidate genes associated with known rice yield-related genes, while gene ranking saves experimental time in the exploitation of rice as a major crop. Future development of our research will include the collection of more rice yield-related genes via online databases and the analysis of literature. Subsequent accurate analysis involving an effective prediction algorithm will enable The confirmation rate was calculated by dividing the confirmation number by the corresponding number of top n. It represented the effectiveness of the confirmation.
FIGURE 4 | GO terms in which the top 100 candidate genes are enriched. The abscissa shows GO terms, and the ordinate represents the number of GO terms.
the prediction of novel genes that can boost rice yield. In the future, we would further develop computational models for the identification and analysis of rice yield-related microRNAs/Long non-coding RNAs based on Chen et al.'s researches (Chen and Yan, 2013;Chen and Huang, 2017;Chen et al., 2017bChen et al., , 2018b.

AUTHOR CONTRIBUTIONS
CW designed the research. XZ performed the research. FX analyzed the data. JJ wrote the manuscript. All authors read and approved the manuscript.

FUNDING
The work was supported by the Natural Science Foundation of China (Nos. 91735306, 61872114, and 61872309).