Prediction of Drought-Resistant Genes in Arabidopsis thaliana Using SVM-RFE

Background Identifying genes with essential roles in resisting environmental stress rates high in agronomic importance. Although massive DNA microarray gene expression data have been generated for plants, current computational approaches underutilize these data for studying genotype-trait relationships. Some advanced gene identification methods have been explored for human diseases, but typically these methods have not been converted into publicly available software tools and cannot be applied to plants for identifying genes with agronomic traits. Methodology In this study, we used 22 sets of Arabidopsis thaliana gene expression data from GEO to predict the key genes involved in water tolerance. We applied an SVM-RFE (Support Vector Machine-Recursive Feature Elimination) feature selection method for the prediction. To address small sample sizes, we developed a modified approach for SVM-RFE by using bootstrapping and leave-one-out cross-validation. We also expanded our study to predict genes involved in water susceptibility. Conclusions We analyzed the top 10 genes predicted to be involved in water tolerance. Seven of them are connected to known biological processes in drought resistance. We also analyzed the top 100 genes in terms of their biological functions. Our study shows that the SVM-RFE method is a highly promising method in analyzing plant microarray data for studying genotype-phenotype relationships. The software is freely available with source code at http://ccst.jlu.edu.cn/JCSB/RFET/.


Introduction
Among all kinds of environmental stresses (abiotic and biotic) in worldwide agriculture, drought is a major abiotic stress factor with significant impact on agricultural production. While resistance to biotic stresses is sometimes associated with monogenic traits, abiotic stresses are typically associated with multigenic traits, making it more difficult to study [1]. Hydropenia can trigger a cascade of physiological and metabolic activities in plants so the tolerance and susceptibility to drought are very complex. Among all the plants, Arabidopsis thaliana is the most popular model organism used in studying drought tolerance of plants, as it is a typical glycophyte, and many other xerophytes or desiccationtolerant plants are similar to glycophytes in the drought-resistant mechanism. Hence, we focus on Arabidopsis to explore the droughtresistant genes in this study.
Recently, several research groups have investigated droughtmediated changes in gene expression using microarrays [2][3][4][5]. Microarray datasets typically include several thousands to tens of thousands of genes with relatively a small number of samples, but many genes are irrelevant or redundant for the purpose of this study. Biologically, there are often tens to hundreds of genes significantly associated to a trait like drought resistance. Hence, it is important to develop computational methods to mine these genes based on microarray data.
Prediction of genes associated with a trait can be formulated as a feature selection problem where key features (genes) of microarray data are indicative of a trait. Various feature selection techniques in handling gene expression data have been proposed. In particular, three types of classification-based methods were developed, i.e., filtering methods, wrapper methods, and embedded methods [6]. Feature selection using embedded SVM evaluation criterion to assess feature relevance is a typical and successful method [7]. Several other studies have provided alternative methods. For example, Support Vector Machine-Recursive Feature Elimination (SVM-RFE) was applied to train SVM for obtaining the weight of each feature and removing the one with the smallest weight iteratively [7][8][9][10][11]. This algorithm is superior to the ''naïve'' ranking with only one time RFE [7,12]. However, this study [13] is the only study to apply the SVM-RFE method for identifying genes of an agronomic trait. Instead, current methods typically use a simple t-test for identifying important genes relevant to a trait [14,15]. This could be problematic for integrating data from various sources. Further-more, many researchers currently study genotype-trait relationships in an ad hoc fashion, e.g. by manually tuning various parameters that are relatively poorly understood [16]. These issues may result in vast amounts of plant microarray data being underutilized for studying genotype-phenotype relationships. One of the reasons behind these issues is lack of publicly available advanced tools. For example, while the SVM-RFE method is highly promising in analyzing microarray data, there was no software tool available.
In this study, we developed a systematic tool by improving the SVM-RFE method for identifying trait-specific genes using microarray data. The tool characterizes drought-resistant genes in Arabidopsis thaliana and is generally applicable to study genotypephenotype relationships using gene expression data from microarray or RNA-Seq. Furthermore, the tool is freely available with source code at http://ccst.jlu.edu.cn/JCSB/RFET/.

Materials and Methods
Data source GEO (http://www.ncbi.nlm.nih.gov/geo/) currently contains 1664 datasets with 90 GDSs, 299 platforms and 1275 series for Arabidopsis thaliana. Among these datasets, we used GSE10670 concerning global expression profiling of wild type and transgenic Arabidopsis plants in response to water stress published on September 1, 2008, and last updated on March 15, 2009. From GEO this is the largest gene expression dataset available up to date for studying plant drought resistance. The data was generated from the Affymetrix platform GPL198 Arabidopsis ATH1 Genome Array. It consists of 22,810 probes and each probe corresponds to one gene. In all we identified 22 samples from GSM269812 to GSM269833, including the wild type (WT), two independent transgenic lines (T6 and T8), and the vector control line (C2), which was used as an additional control. In all related experiments, the relative water content (RWC) of wild type and transgenic leaves during a period of dehydration was monitored. At day 7 when the transgenic plants were still at an RWC of .85%, the wild type and vector control plants were at an RWC of ,50-60% [17]. Hence, we first chose experimental samples just on T6 and T8, which were expected to reveal drought tolerant genes. Table 1 gives the specific description of the data used, and the class label (0/1) is according to the different stress conditions. Table 2 presents the susceptibility genotype (WT and C2) data, which were used to discard the tuning genes as described later.

Data preprocessing stage
We applied a quantile-based [18] RMA (Robust Multi-chip Averages) method for normalizing microarray data. The RMA feeds probes data stored in Affymetrix CEL into a stochastic model to estimate gene expression and converts the probe data to gene expression data. We conducted the RMA analysis using Bioconductor (http://www.bioconductor.org/), which is an open-source tool for bioinformatics using the R statistical programming language.

T-test method for preliminary selection
For tens of thousands of genes, it would be of high complexity to use the SVM-RFE directly. Hence, we first employed a t-test [19] to filter out unlikely genes involved in drought tolerance. In our preliminary selection we assigned 0.001 as the p-value threshold, resulting in 736 genes, which are still too many for agronomic studies.
Using SVM-RFE method for gene selection RFE is an iterative procedure for SVM classifier. A cost function J computed on training samples is used as an objective function. Expanding J in Taylor series to the second order using the OBD algorithm [20], and neglecting the first order term at the optimum of J, yield: Here (w i ) 2 was used as the ranking criterion and we used LIBSVM (a library for Support Vector Machines) [21] with a linear kernel. We present below an outline of the SVM-RFE in the linear kernel. For more details about this method, see Guyon et al. [7].

Inputs: Training samples (microarray datasets)
X 0~½ x 1 ,x 2 ,:::,x 12 T Class labels (1 for well watered or 0 for drought) Repeat until s~½ Output: A gene-ranking list r. We trained the classifier, computed the ranking for the 736 genes obtained and then removed the gene with the lowest ranking. We repeated the process until all the genes were removed. This iterative process is a sequence backward selection (SBS) procedure and at last the method produces a gene-ranking list with weights from high to low.
To improve prediction accuracy we conducted several rounds of bootstrapping in the SVM-RFE procedure and in each round one sorted list was produced. However, there is a shortage of experimental samples that are needed to train the SVM. This brings up two issues: one is how to generate training and test sets, and the other is how to combine the weights of each gene in each sorted list. To address these two issues, we developed the following two solutions: 1) In order to make good use of limited data for predicting drought-resistant genes in Arabidopsis, the generation of training set is a key factor. The dataset was randomly split into n subsets of approximately equal size, then one subset was removed, and the remaining samples formed the training set. Each time a different subset was selected in such a way that all the samples had an equal chance to be selected as the training data. We call it n-CV, where n could be equal to 12, 6, 4, or 3, respectively. In the following examples, there are 12 samples from transgenic lines described in Table 1, and a 6-CV was used for each subset with two samples. Then the Leave One Out Cross Validation (LOOCV) was used for training SVM (see Figure 1). After that we performed 100 times of 6-CV and each 6-CV ran 6 times of the SVM-RFE procedure. The computational complexity of each n-CV is n*(g*O(s 3 )), where s is the number of samples and g is the number of genes [22]. The computational time on the training stage primarily depends on the number of genes given the small number of training samples. Hence, the computational complexity can be easily handled by filtering out unlikely genes involved. The 100 times 6-CV training took around 90 minutes to get the results on a desktop computer with Intel Core 2 Duo E6750 and DDR2 2GB memory. 2) We acquired 100 sorted lists, and designed a re-ranking measure to take occurrence and ranking of every gene into account to form a final ranking list of all genes: where p is the number of times for CV, W f is the sum of one gene's weights in p experiments, and W L k is the weight in the k-th occurrence in sorted lists. L k shows the ranking of one feature in the k-th list.
Here, we have top 10 (N) genes and 100 (p) SVM-RFE trainings. Thus, we obtain a final list containing the optimal genes sorted by W f in decreasing order. The algorithm flowchart is shown in Figure 2.

Results and Discussion
By applying the SVM-RFE method in analyzing the droughtresistant genotype based on the microarray data, we selected 10 genes in the final list according to W f in Eq. 2, as summarized in Table 3. Figure 3 shows the occurrence of these genes in top-10 list and also in top-30 list when conducting 100 times of 6-CV. It shows that the occurrence of these genes in the top list is very high and consistent. We checked the functional annotations using GO (http://www.geneontology.org/), KEGG pathways (http://www. genome.jp/kegg/) and the literatures [2,[23][24][25][26][27][28], and obtained the related information shown in Table 3. From Table 3 it can be seen that most of the genes have no specific molecular function annotations. This is not surprising as the water stress tolerance is a complex trait, which resulted from sophisticated coordination of physiological and biochemical alterations at the cellular and molecular levels [24]. Table 3 reveals some interesting biological features. Both the first and second genes directly respond to water deprivation. The cellular component of the second gene is chloroplast where photosynthesis of plants occurs. As photosynthesis uses carbon dioxide and water releasing oxygen, it is likely that this gene plays some role in water utilization in chloroplast. The third and fifth genes are both related to defense response, whose pathway is known to have cross-talk with drought resistance pathway [29].
The sixth gene in the table responds to wounding and it may help repair the damage of cell caused by water loss. Some other genes in the list may also be related to drought resistance. The eighth-ranked gene responds to osmotic stress and salt stress. Osmotic adjustment is an important physiological mechanism adapting to water stress [24]. Osmotic adjustment can maintain a dynamic balance between damage and repair of cellular components to relieve plants injury and improve plants' ability of stress resistance.
In the column GO: Component of Table 3, all of the 3 rd , 5 th , 6 th , 7 th , and 8 th genes belong to membrane systems, like endomembrane, plasma membrane, and so on. Membrane system is the key part damaged by drought stress and it is the most sensitive original reaction site against adversity [25]. Membrane, together with associated proteins, provides cells with not only a relatively stable internal environment, but also provides a switch to material transportation, energy exchange and information transmission between cells and the environment. Therefore, these five genes may help adjust osmotic membranes to boost drought resistance. In all we have demonstrated that seven genes may be closely related to water tolerance for Arabidopsis, i.e., the 1 st , 2 nd , 3 rd , 5 th , 6 th , 7 th , and 8 th genes in Table 3.
We also repeated the computational process with the susceptibility genotype samples (WT and C2), and obtained the top-10 gene list as described in Table 4. From the relative function annotations, it appears that most of these genes have little relationship with the ability of drought resistance, which may explain why these genotypes do not have capacities of drought resistance.
Our understanding of the functions of these seven genes is far from complete. Compared the top-10 gene list obtained from the resistant genotype with the susceptibility genotype list, 3 genes are the same, which are the 1 st , 2 nd , and 8 th in Table 3, and the 10 th , 8 th , and 5 th in Table 4. We call the same ones tuning genes [16]. Maybe their adaptability to hydropenia is the result of irritable reactions to environmental changes. So by removing the tuning genes from the top-10 gene list with transgenic genotype, the inference is that the real drought-resistant genes should be in the 7 ones in Table 5 (a subset of Table 3). And with the same thought, we analyzed the two top 100 gene lists using the resistant genotype and the susceptibility genotype, respectively (see Table S1 and Table S2 for the detailed information). Comparing the top-100 gene list obtained from the resistant genotype with the susceptibility genotype list, it can be seen that 37 genes overlap (tuning genes). In Table S1, the ''Overlap'' column indicates the tuning genes. We compared our result against the result published by Huang et al. [30] to look for the overlap of the genes identified to be involved in drought resistance. Our comparison shows that 5 out of the 10 genes from resistant genotype and 6 out of the 10 genes from susceptible genotype are identified to be the same. When we expanded our analysis to include the top 100 genes, 50 genes of the resistant genotype and 42 of the susceptible genotype overlap with the gene lists in the published results. GO term enrichment analysis was also performed for the annotations of the top 100 genes in the resistant and susceptibility genotypes, respectively using Amigo GO Term Enrichment Tool (http:// amigo.geneontology.org/cgi-bin/amigo/term_enrichment) as shown in Table 6. All the GO terms identified for both genotypes with a significant p-value less than 0.001 have functional categories related to the drought stress.  Table 3. doi:10.1371/journal.pone.0021750.g003 The SVM-RFE algorithm was used in cancer marker gene prediction with abundant datasets. However, there are much fewer microarray samples for finding the drought-resistant genes of Arabidopsis. Hence, we modified the original SVM-RFE method to address this drawback by using bootstrapping and leave-one-out cross-validation. Our study shows that the improved method is effective for identifying drought-resistance genes. Since the sparseness of gene expression data for studying genotype-trait relationships is a common issue, the method provides a framework for handling this issue. The framework has some advantages over some other feature selection methods that require extensive training data, such as random forest [31]. This is by no means a replacement of additional experimental data, but it can effectively utilize the sparse data available to generate useful hypotheses and guide further targeted experimental work.
Although Arabidopsis is a model organism for plant gene function analysis and gene expression studies, few genes related to drought resistance mechanisms are annotated with direct experimental evidences. There is a need to predict additional drought-resistance genes based on gene expression data. The predictions of drought-resistance genes generated by the JU/MU development of the SVM-RFE method-based software provides useful hypothesis for experimentalists to verify. For example, the 3 rd , 8 th and 10 th genes in Table 3 are not annotated as droughtresistance genes, but they are highly likely involved in drought resistance. Perhaps a major challenge in the future is to inquire into the relative contribution of each gene to water tolerance. The JU/MU approach is applicable to the study of plant genes related to other stress resistance and genes associated with any agronomic trait in general.  Supporting Information