Statistical Approaches for Gene Selection, Hub Gene Identification and Module Interaction in Gene Co-Expression Network Analysis: An Application to Aluminum Stress in Soybean (Glycine max L.)

Selection of informative genes is an important problem in gene expression studies. The small sample size and the large number of genes in gene expression data make the selection process complex. Further, the selected informative genes may act as a vital input for gene co-expression network analysis. Moreover, the identification of hub genes and module interactions in gene co-expression networks is yet to be fully explored. This paper presents a statistically sound gene selection technique based on support vector machine algorithm for selecting informative genes from high dimensional gene expression data. Also, an attempt has been made to develop a statistical approach for identification of hub genes in the gene co-expression network. Besides, a differential hub gene analysis approach has also been developed to group the identified hub genes into various groups based on their gene connectivity in a case vs. control study. Based on this proposed approach, an R package, i.e., dhga (https://cran.r-project.org/web/packages/dhga) has been developed. The comparative performance of the proposed gene selection technique as well as hub gene identification approach was evaluated on three different crop microarray datasets. The proposed gene selection technique outperformed most of the existing techniques for selecting robust set of informative genes. Based on the proposed hub gene identification approach, a few number of hub genes were identified as compared to the existing approach, which is in accordance with the principle of scale free property of real networks. In this study, some key genes along with their Arabidopsis orthologs has been reported, which can be used for Aluminum toxic stress response engineering in soybean. The functional analysis of various selected key genes revealed the underlying molecular mechanisms of Aluminum toxic stress response in soybean.


Introduction
With the advent of fast and cheaper genome sequencing technologies, huge genomic data is being generated and deposited in public domain databases over the years by different research organizations across the globe [1,2]. Most of these datasets are related to expression of genes from various experiments conducted to understand behavior of biological mechanism of species under biotic and abiotic stresses. In due course of time huge gene expression data is generated through microarray experiments under these stresses. Integration and analysis of data generated by microarray experiments for the same stress or related conditions is essential to enhance the sensitivity of the hypothesis under consideration for drawing valid conclusions [3]. For instance, meta-analysis of microarray data pertaining to different experiments in rice and Arabidopsis revealed the presence of highly connected key genes that are central to the plant defense system under various biotic and abiotic stresses [4,5].
Usually, microarray data are used for gene selection and modules detection in genetic network analysis, which suffers from the inherent limitation of its high dimensionality, i.e. the number of genes is much larger than the number of subjects/samples [6]. Therefore, it is important to select most relevant genes related to stresses/conditions from thousand(s) of genes with the help of appropriate computational approach(s). In this regard, Volcano plot method [7] is quite popular among the researchers in which genes are selected by considering their relevance with their classes. However, such method may not be sufficient to discover some complex relationships among genes for a certain trait or condition [8]. Besides, several statistical and machine learning methods, viz. t-score, F-score, Information Gain (IG) measure, Random Forest (RF) and Support Vector Machine-Recursive Feature Elimination (SVM-RFE) [7,[9][10][11][12][13] have also been used for gene selection. However, in these methods genes are selected by considering only their relevance with classes. In such case, there is a possibility that genes which are spuriously associated with the classes may also get selected.
In order to understand the interrelationship among the selected genes, identification of gene modules and key genes responsible for a particular stress/condition, analysis of gene coexpression networks need to be carried out. Weighted Gene Co-expression Network Analysis (WGCNA) [14] is a latest and popular technique used to decipher co-expression patterns among genes. The WGCNA approach typically deals with the identification of gene modules by using the gene expression levels that are highly correlated across samples [14]. This technique has been successfully utilized to detect gene modules in Arabidopsis, rice, maize and poplar for various biotic and abiotic stresses [5,[15][16][17][18]. Further, this approach also leads to construction of Gene Co-expression Network (GCN), a scale free network, where, genes are represented as nodes and edges depict associations among genes [14,19]. In such network, highly connected genes are called hub genes, which are expected to play an important role in understanding the biological mechanism of response under stresses/conditions [20][21][22][23][24]. Identification of hub genes will also help in mitigating the stress in plants through genetic engineering. The existing approaches [21][22][23][24] have mainly focused on hub gene identification, based only on gene connection degrees in the GCN. Moreover, these techniques select such genes empirically without any statistical criteria. Besides, few approaches can be found in the literature for the identification of hub nodes in a scale free network [22][23][24].
Aluminum (Al) toxic stress is a major impediment to the crop production on acidic soils that affects about 30-40% of the world's arable lands [25]. Soybean (Glycine max L.), that provides major source of proteins, unsaturated fats, carbohydrate and fibers, is one of the most important legume crop, capable of providing nutritional security to the global population. Soybean is preferably grown on acidic soil and its productivity is significantly reduced by Al toxic stress. In acidic soil, Al stress causes rapid inhibition in root growth and subsequently inhibits water and nutrient uptake by plants. This increases the susceptibility of plants to other environmental stresses and results in reduction of crop productivity [26]. Under heavy pressure of population explosion and global warming, achieving nutritional security in general and protein security in particular through enhancing the productivity of soybean is of paramount importance. However, the underlying mechanisms for Al toxic stress response in plants in general and soybean in particular are not so clearly deciphered till now [27].
In this study, a statistical technique i.e. Bootstrap SVM-RFE (Boot-SVM-RFE) is proposed for selection of informative genes. In this technique, genes are selected after reducing the effect of spurious associations between genes and classes. The performance of the proposed gene selection technique is found to be better than the existing techniques, while compared by using three different datasets. Further, a statistical approach for identification of hub genes in the GCN was also proposed. Again, this approach is evaluated on the genes selected from the above datasets and found to be superior in terms of scale free property of biological network. Besides, an R package has been developed based on the proposed hub gene identification approach. Further, an attempt has been made to integrate and analyze the gene expression datasets generated by different experiments for the identification of Al toxic stress responsive genes in soybean, by using the proposed techniques. Hub genes responsible for Al toxic stress have been identified and their functional analysis has been done.

Materials and Methods
The soybean microarray experimental datasets under Al stress were collected from Gene Expression Omnibus with platform GPL4592 (http://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GPL4592). This platform contains 3855 experimental samples on 37,593 probes generated using Affymetrix Soybean Genome Array. Out of these samples, 80 samples related to Al stress were collected for further study. Initially, raw CEL files of these collected samples were processed using Robust Multichip Average (RMA) algorithm available in affy Bioconductor package of R [28][29][30]. This includes background correction, quantile normalization and summarization by the median polish approach [31]. Then the microarray experimental samples with mean ! 7.1 and standard deviation 2.5 were selected, as uniformity of colors in the correlation plot was observed for these parameters setting (S1 Fig). Through this process, 78 samples (generated over 3 different experiments) were selected. The descriptions about the selected samples are given in S1 Document. The log2 scale transformed expression data from the RMA for these selected experimental samples were used for further statistical analysis. J was given by where, Δw i is change in weight due to removing i-th gene from existing dataset and ΔJ(i) is used as weight pruning criterion. It may be noted here that, the cost function J is a quadratic function of w i and both are directly proportional to each other. Hence, measurement of either w i or J provides equivalent information. Keeping this in view w 2 i is used as the ranking criterion for evaluating impact of ith gene on this classification [6]. In this process, genes are eliminated with the smallest w 2 i iteratively in a backward elimination manner and ranked gene list is prepared at the end. Moreover, most of the gene selection methods are sensitive to small permutation of experimental conditions [13]. The ranking of genes based on high dimensional expression data may also lead to the selection of spurious genes and make the selection process unreliable [33]. Therefore, it is essential to select genes based on statistical testing instead of their ranks. Keeping in view the above fact, a test statistic has been proposed for selecting informative genes.
In this testing procedure, n bootstrap samples each of size m are selected randomly with replacement to construct a training set for SVM from available M samples in a dataset. Then SVM-RFE procedure was applied to each of these n bootstrap samples to get n list of genes along with their ranks. Therefore, each of genes will have n number of ranks (one for each bootstrap). Let a score function i.e. Rank Score (R ij ) is defined to convert these ranks of each gene into corresponding score in each bootstrap sample, as where, N represents total number of genes considered in the dataset and p ij (1 p ij N) is the ranked position of i-th gene in j-th bootstrap sample. After getting the rank scores of all genes over n bootstrap samples, following proposed hypothesis needs to be tested for selection of informative genes. 1) is the rank score for j-th bootstrap sample (j = 1, 2, . . ., n). Further, R j (i) is a random variable (rv). Since, R j (i) is a function of rank, therefore, its empirical distribution is symmetric about the second quartile. So, without loss of generality, we defined another variable r j (for fixed i) as: In order to test the statistical significance (H 0 vs. H 1 ) for gene i, the r j 's are arranged in ascending order of their magnitude and subsequently, the ranks 1, 2, . . ., n are assigned, keeping in mind their original signs. Let T + be the sum of the ranks of positive r j 's and T − be the sum of the ranks of negative r j 's. Now for finding distribution of test statistic T + , another variable Z (k) is defined as: Here, k = {1, 2, . . .n}. Now, the variables Z (k) are independent Bernouli variates and its order of moments can be obtained as: Var ðZ ðkÞ Þ ¼ EfZ ðkÞ ð1 À EðZ ðkÞ ÞÞg ð7Þ Then, the first two moments of the statistic (T + ) can be written as: VarðT þ Þ ¼ Let R i be expected rank score for i-th gene over all bootstrap samples. Under the simple null hypothesis H 0 : R i = Q against H 1 : R i > Q the expressions in Eqs 8 and 9 can be written as As the number of bootstrap samples are quite large, then under Linberg's cental limit theorem [34,35], T + follows normal distribution asymptotically, i.e.
Based on test statistic (Eq 12) under H 0 , it can be tested that whether a gene is informative or not. This procedure was repeated for other genes. The performance of the proposed technique for selection of informative genes was compared with respect to existing techniques viz. SVM-RFE, t-score, F-score, RF and IG on gene expression data of Al stress in soybean, salinity and cold stress in rice (Table 1). In order to assess the performance of genes selection techniques, the fixed number of top ranked genes selected through these techniques were used in SVM classifier to discriminate the class labels of samples between stress (+1) and control (-1). In this case, SVM-Classification Accuracy (CA) was computed through a sliding window size technique. Here, the window sizes (number of sorted genes) were taken as 50, 100, . . ., 500 with sliding length of 50. The performance of these techniques was adjudged on the basis of CA and it's co-efficient of variation (CV).

Gene co-expression network analysis
GCNs were constructed by using gene co-expression measure that depicts association among genes [23,36]. Let x i be the expression profile of i-th gene, i.e. the expression values of i-th gene across all the microarray samples. Then, gene co-expression similarity measure s ij between i-th and j-th gene is computed as the absolute value of Pearson's Correlation Co-efficient (PCC) [14,23], which is given by: The adjacency score (a ij ) between i-th gene and j-th gene is defined in terms of s ij [18] as: where, β (!1) is soft threshold power, determined by using the concept of scale free property of biological networks [23]. The detail methodology for determination of the soft threshold power has been discussed elaborately by Zhang and Horvath (2005) [14]. This soft threshold approach leads to a weighted GCN that satisfies the scale free property of biological networks. For both Al stress and control conditions, the value of β was taken as 8 for calculation of adjacency score (S2 Fig), with best approximation to scale free criteria [36] using R 2 > 0.80 through fitting of Power law model. In order to identify the gene modules (i.e. group of tightly co-expressed genes) within the selected informative genes, the Topological Overlap Matrix was constructed based on the adjacency scores [14]. The BlockWiseModules function available in WGCNA package [19] of R was executed to identify these modules. For this purpose, various parameters like module size, deep split level and tree merge cut height was set at 20-30, 4 and 0.15-0.25 respectively. In order to find the consensus modules showing co-expression patterns of genes across stress and control conditions, the function blockwiseConsensusModules was used with parameter settings 8, 30 and 0.15 as power, minimum module size and merge cut height respectively.

Proposed statistical approach for identification of hub genes
In network theory, a node is defined as hub node [20][21][22][23][24], if its connection degree is greater than average connection degree of the network [21]. In the existing approach, a gene is declared as hub gene based on an indicator function [21,23], i.e. Hub i = [I(ki > τ)] and number of hub genes (NHub) in the genetic network is calculated as NHub ¼ where, Hub i : hub status of i-th gene (i.e.1 or 0); k i : connection degree of i-th gene; τ: threshold value i.e. average connection degree of the network. This technique selects hub genes empirically based on only observed gene connectivity without taking into account any statistical consideration. Therefore, an alternate statistical approach based on statistical significance of gene connectivity was proposed for detection of hub genes in the GCN. The proposed statistical approach is described as follows: The Weighted Gene Score (WGS) for i-th gene in terms of weighted gene connectivity (a ij ) can be written as: where, WGS i represents the relative importance of i-th gene based on its connections to the remaining genes in GCN. For the purpose of hub gene identification, following hypotheses are constructed.
H 0 : WGS i m i:e: i À th gene in the GCN is not a hub gene where, μ is average connection degree of the complete network model. Here in order to get the distribution of the test statistic under H 0 , a resampling procedure was used. In this procedure, m microarray samples were selected randomly with equal probability from M microarray samples to construct one subsample (for one GCN) (m M). Then statistical measures (Eqs 13-15) were applied to get WGS for each gene in that GCN. This procedure was repeated large number of times say S to get S sets of WGS. In this study, S = 500 was taken to get 500 random GCNs under stress and control conditions separately. For testing H 0 vs. H 1 , a NP test statistic was proposed to test significance of the WGS for each gene, i.e. for testing whether WGS of a gene is greater than the average connection degree of the complete network or not. The proposed procedure for testing the hypothesis is as follows: Let for a particular gene (i), WGS k (i) be the WGS for k-th subsample (k = 1, 2,. . ., S). Here WGS k (i) 's are rvs. So, without loss of generality, another variable X k can be defined as: In order to test the statistical significance of connectivity for gene i, the X k 's are arranged in ascending order of their magnitude and subsequently, the ranks 1, 2, . . ., S are assigned, keeping in mind the original signs of X k . Let W + be the sum of the ranks of positive X k s and Wbe the sum of the ranks of negative X k s. The distribution of the test statistic (W + ) under H 0 can be obtained by following the above approach to get the distribution of T + in Boot-SVM-RFE. Further, under large number of subsamples (S = 500), following the central limit theorem the distribution of W + is approximately standard normal, i.e.
The procedure was repeated for each gene in the GCN and the statistical test was applied to identify hub gene based on significance values for both control and stress conditions separately. Algorithm: Step 1: Begin with all genes (nodes) in the GCN Step 2: Construct a data set say T k with m samples randomly taken from M microarray samples Step 3: Calculate WGS for all genes Step 4: Repeat Step 2 and 3 S times to get S sets of WGS for each gene Step 5: Take a particular gene (i-th gene) along with its WGS Step 6: Test the hypothesis for i-th gene and obtain its p-value Step 7: Repeat the Step 5-6 for all genes (i = 1, 2, . . ., G) Step 8: Rank the p-values and select the hub genes The proposed hub gene identification approach for the GCNs constructed under two contrasting conditions (stress vs. control) can be called as Differential Hub Gene Analysis (DHGA). By this approach, the identification of hub genes is possible in both these GCNs based on statistical test of significance. On the basis of p-values, genes in the GCNs under either condition can be grouped into various groups, viz. Housekeeping Hub Genes (HHG), Unique Hub Genes (UHG) for stress, UHG for control, Non-hub genes based on a decision matrix ( Table 2).

Modeling of Module Interaction Network
A Module Interaction Network (MIN) can be defined as a digraph, where modules are nodes or units and edges depict regulatory relationships among the modules. Most of the approaches available in literature for modelling of genetic networks are not applicable to the gene expression data due to the problem of high dimensionality [37]. It is biologically required to consider gene module as a functional unit [18,20]. Further, these modules can be taken as unit to study the interaction among the gene modules. The expression levels of the modules can be calculated as: For Al stress data, all the selected microarray samples obtained through meta-analysis belong to time series experiments (the gene expression values were measured over 5 time points e.g. 0, 2, 12, 48 and 72 hours) (S1 Document). Further, the bspline method of data interpolation was used to interpolate the module expression values up to 50 time points in the interval of [0, 72 hours]. Then, we modeled the expression level of module d at time t as a linear regression with the expression levels of other modules at time (t−1). The model which depicts the interaction between these modules can be written as: where, the β h s are regression coefficients and ε is the random noise with mean 0 and variance σ 2 . To compute the regression coefficients, which depicts the interaction among the modules, Bayesian Model Averaging (BMA) [38] algorithm was executed. Further, the posterior probabilities for each module interaction were calculated by using iterative BMA algorithm [39]. The module interactions were arranged in descending order by the value of posterior probabilities and significant module interactions were selected by fixing proper threshold. The MIN was constructed by using RCytoscape package [40].

Performance analysis of Boot-SVM-RFE
In order to study the performance of the proposed Boot-SVM-RFE technique, the top 1000 genes obtained based on ranking from each of the gene selection techniques were used for classification of crop microarray samples into control and stress classes through SVM classifier. The CAs were measured for each sliding window size over 5 fold cross validation. The CAs for different sliding window sizes is given in Table 3. It is observed that for Al stress data, the CAs of Boot-SVM-RFE are higher than that of other techniques viz. SVM-RFE, t-score, F-score, RF and IG for the sliding window sizes 50, 100, 150, 200, 250 and 300. However, for higher window sizes i.e. 350, 400, 450 and 500, the CAs for Boot-SVM-RFE are at par with that of SVM-RFE but higher than that of t-score, F-score, RF and IG (Table 3). In general the performance of Boot-SVM-RFE is highest followed by SVM-RFE, RF, IG, F-score and t-score with respect to CA values for Al stress. However, in case of salinity and cold stress, the performance of Boot-SVM-RFE is observed to be better than other gene selection techniques irrespective of sizes of sliding window ( Table 3). The order of the performance of different gene selection techniques in case of salinity stress is Boot-SVM-RFE > SVM-RFE > F score > RF > F score > t score > IG, whereas, in case of cold stress the order of performance is Boot-SVM-RFE > RF > IG > SVM-RFE > F score > t score (Table 3). From this performance analysis, it can be seen that the performance of the proposed Boot-SVM-RFE is consistently better over other contemporary techniques across different datasets related to abiotic stresses. Further, it is also observed that Boot-SVM-RFE has less CV for most of the sliding window sizes with respect to other five techniques when applied to the datasets of these three stresses ( Table 3).

Selection of informative genes for Al stress in soybean
Since the Boot-SVM-RFE was found to be superior as compared to other gene selection techniques, it was further employed to select informative genes for Al stress in soybean. In order to get a robust and minimal set of informative genes, the fold change in Volcano plot was replaced with-log10 (p-values) obtained from Boot-SVM-RFE and then a gene selection plot was constructed. The threshold values for Y and X-axis of the gene selection plot were fixed as 4 and 2.5 respectively, which lead to selection of 981 genes (Fig 1). The consensus sequences of these 981 genes obtained from GeneChip Soybean Genome Array of Affymetrix were then used to identify the Arabidopsis orthologs [41] and it was found that 554 genes have unique orthologs in Arabidopsis (Fig 1). Further, the annotations of these selected genes were obtained from SoyBase (http://soybase.org) [42]. A brief description about these selected genes is given in S1 Table. Functional analysis of selected genes for Al stress in soybean The Gene Ontology (GO) enrichment analysis of the 981 selected genes was performed by using AgriGO [43], a plant-specific GO term enrichment analysis tool. It is observed that most of the selected genes are responsible for transition metal ion binding, metal ion binding, cation binding, ion binding, etc. (Fig 2A). These molecular functions (MF) might be activated due to high concentration of Al ions in water or soil. Two other MF i.e. oxido-reductase (redox) and kinase activities are also present in these selected genes (Fig 2A). The significant behavior of the genes in redox activity might be related to electron transport in complex chemical reactions that balances the charges during ion transport. The redox activity might also be related to Reactive Oxygen Species (ROS) that are produced in response to oxidative stress due to water deficit during abiotic stress like Al toxic stress [44]. In biological process categories, such as cellular nitrogen compound metabolic process, amine metabolic process, cellular amino acid and derivative metabolic process, oxoacid metabolic process, organic acid metabolic process, carboxylic acid metabolic process, cellular ketone metabolic process and ion transport activity, the number of selected genes is more as compared to other biological processes (Fig 2A). It may be inferred that some of these chosen genes are involved in ion transport activities, i.e. involved in transporting the ions outside the cell to maintain the proper pH in the cell [45]. In case of cellular components, chosen genes are related to transcription factor complex, cytoplasmic membrane-bounded vesicle, membrane-bounded vesicle, cytoplasmic vesicle, vesicle and nucleoplasm part (Fig 2A). It can be seen that the maximum number of the genes is related to vesicle and membrane, which is consistent with the detoxifying mechanism of metal ions available in Al stress condition, especially in sequestration by vacuole [46,47]. Some of the selected genes present on membrane are found to be involved in transporting of metal ions outside the cell or to the vacuole to maintain pH and transmembrane proton gradient [48].

Gene co-expression network analysis for Al stress in soybean
Using WGCNA, the selected 981 genes were divided into 19 and 18 modules (including grey colour module, which is the module of the non-modular genes) for Al stress and control conditions respectively (Fig 3). In both the cases, module represented by turquoise colour contains maximum number of genes, hence designated as the largest module for either condition. Based on the expression profiles of these selected genes for both Al stress and control conditions, 23 consensus modules (set of genes with similar co-expression patterns) were obtained.   Table 4. It is observed that, every module is significantly annotated with GO terms, except gene modules represented by green-yellow and grey colour (Table 4). So, it can be inferred that functions of genes present within these two modules are still largely unknown. Furthermore, 19 identified modules (including grey module) were used as functional unit to model the MIN. The posterior probabilities of each module interactions were computed (using iBMA algorithm) and used to construct the MIN for Al stress as shown in Fig 4. It is seen that modules 1 and 18 are interacting with most of the modules in the MIN (Fig 4). Thus, these modules along with their members may play a significant role in Al stress response in soybean.

Performance analysis of proposed hub gene detection approach
Based on the existing approach (i.e. WGS alone), 39.05% and 36.91% genes in the GCN are detected as hub genes in soybean for stress and control conditions respectively (Table 5). Thus, Statistical Approaches for Gene Co-Expression Network Analysis large proportions of genes are identified as hub genes in the GCN based on the existing approach, which contradicts the scale free property of biological networks (as GCN is a scale free network) [20,36]. Similar findings are observed for salinity and cold stresses in rice (Table 5). However, in case of proposed approach (computing p-values) only 23.24% and 19.14% of genes in the GCN are found to be hub genes in soybean (for p value < 1E-10) for stress and control situations respectively (Table 5). Moreover, the number of hubs can be further reduced by decreasing the level of significance in the proposed approach. This indicates that the proposed approach will be able to identify relatively small subset of genes as hubs in the GCN i.e. fewer WGS are statistically significant. Similar results were obtained for the salinity and cold stresses in rice ( Table 5). The distributions of WGS (i.e. are heavy right tailed distributions) contain lower and upper values, which are not much discriminated between low and high connection degree of genes (Fig 5). On the contrary, from the distribution of p-values, the genes with high connection degrees are well separated from that of low connection degrees in the GCN (Fig 6). In other words, the distinction between statistically strongly and weakly connected genes in the GCN can be better seen from Fig 6 than Fig 5. Further, the results obtained from the proposed DHGA approach for two contrasting conditions (stress vs. control) are given in Table 6. The DHGA approach allowed grouping of 187 (stress) and 208 (control) hub genes of rice in the GCN into HHG (141), UHG to salinity stress (46) and UHG to control (67) ( Table 6). Similar interpretations can be made for Al stress in soybean and cold stress in rice.

Differential hub gene analysis for Al stress condition in soybean
Following the above approach, 228 and 187 genes were identified as hub genes whose p-values were 1E-10 for Al stress and control conditions of soybean respectively (Table 6). From the DHGA result, it is seen that 98 hub genes are common whereas 130 and 89 hub genes are unique for Al stress and control conditions respectively (Fig 7C). The mapping of the HHG and UHGs in soybean to Arabidopsis genome leads to the identification of corresponding Arabidopsis orthologs genes (Fig 7C). The GCNs constructed for these two differential conditions (Al stress vs. control) in soybean along with the positions of hub genes and UHGs are shown in Fig 7. The functional analysis of the selected hub genes under Al stress revealed their associated cellular mechanisms. From the GO analysis (in cellular components), it is observed that most of the hub genes are present in plastid, vacoule, membrane bounded vacuole, cytoplasmic vacuole (Fig 2B), mainly responsible for pumping out ions from the cell. In MF category, majority Statistical Approaches for Gene Co-Expression Network Analysis of the hub genes were found to be involved in nucleic acid, cation ion, metal ion and zinc ion binding activities (Fig 2B), which may be responsible for fixing metal ions. Further, a large portion of the hub genes were found to be responsible for nitrogen compound metabolic process, response to stress and chemical stimulus, defense response and signal transduction under biological process category. The module membership of the hub genes as well as UHG under Al stress showed that most of the hub genes under Al stress condition belong to turquoise (86), yellow (45) and brown (38) modules (Table 4). Similarly, out of 130 UHG, mainly 45, 32 and 26 are the members of turquoise, brown and yellow colour modules respectively. Interestingly, it can be seen that the blue colour module contains the second highest number of selected genes (137), but have no hub genes, while the brown colour module is the third largest module (100 genes) contains 38 hub genes out of which 32 are UHG for Al stress. Further, brown colour module is found to be associated with various important functions like ion binding, redox activity, kinase Statistical Approaches for Gene Co-Expression Network Analysis activity, phosphortransferase activity (Table 4), which are important for the abiotic stress response in plants [45]. From the molecular biology point of view, the brown colour module along with its members seems to be very important for breeding Al stress resistant varieties. A brief description about hub genes and UHG for Al stress condition is provided in S2 and S3 Tables respectively.

Discussion
The proposed Boot-SVM-RFE technique was found to be superior for selection of informative genes from high dimensional gene expression data. This approach is also advantageous over classical gene selection techniques like t-test and F-score, as it does not require any distributional assumptions about the data. In this technique, a p-value was assigned to each gene and genes with lower p-values were considered as informative for the particular condition/trait under investigation. The selection of informative genes based on p-values is scientific as well statistically meaningful to experimental biologists as compared to other techniques. Further, the bootstrap procedure used in this technique was expected to remove the spurious associations of the genes with their classes. The comparative analysis showed that the Boot-SVM-RFE performed better than existing techniques i.e. SVM-RFE, t-score, F-score, RF and IG in terms of CA. Besides, its performance can be considered as robust due to the lower CV values in CA for all window sizes. The proposed statistical approach for hub gene identification allowed the ranking and selection of candidate hub genes in the GCN, based on an assessment of the statistical significance of the gene connections. This was done with a randomized resampling based procedure where statistical significance values were calculated based on the NP test, which does not require Gaussian assumptions of data. Further, genes with lower p-values represent highly connected genes in the GCN and thus designated as hub genes. Moreover, the randomisation procedure used in this approach allows one to test, whether the observed gene connectivity is greater than expected gene connectivity value by chance (i.e. rejection of null hypothesis of random association). This was also able to remove the spurious association among genes, as these associations are measured on the basis of PCC. It seems to be more statistically convincing to select hub genes based on p-values rather than WGS alone, because in comparison to WGS, the p-values provides a reliable measure of gene connectivity based on a statistical criterion (lower p-value indicates high gene connectivity and vice-versa). Further, the detected hub genes tend to have Statistical Approaches for Gene Co-Expression Network Analysis higher connection degrees and are widely separated from the genes with low connection degrees in the GCN. Moreover, based on this approach a few as well as important genes were identified as hubs in the GCN as compared to existing approach, which is in accordance with the scale free property of biological networks.
Using DHGA approach, genes in the GCN were grouped into various categories like HHG, non-hubs, UHG for stress and control based on the computed p-values for these two contrasting conditions. These identified hub genes may be considered as biomarkers for further studies, including analysis of their involvement in diverse cellular mechanisms. Further, the HHG can be used for the maintenance of basal cellular functions that are essential for the existence of a cell [49], whereas, UHG can be used in stress response engineering in crops for developing stress tolerant cultivars.
Understanding Al stress response mechanism in soybean is of paramount importance for plant breeders to develop Al stress tolerant cultivars. In public domain databases, there are few samples available related to Al stress in soybean, which have been generated over varying experimental conditions by multiple studies. Thus, meta-analysis was performed to combine these datasets and the meta-data was used for further statistical analysis. Then, developed techniques were applied to identify the responsible genes to understand stress response mechanism in this crop. It has been reported that there are two main processes involved in Al stress response in plants (i) exclusion of Al ions from root cells and (ii) detoxification of Al ions in the plant cells [50]. Some selected genes were found to be involved in transporting of metal ions outside the cell, which might be associated with the first process. The function like redox activity related to electron transport under chemical reactions that balances the charges during metallic ions transport [45], might be associated with the second process. The redox activity might also be related to ROS generation that is produced in plants in response to Al stress. Further, ROS also seriously disrupts normal metabolism of cell through peroxidation of lipids [51], proteins and nucleic acids [52]. The increased redox activity is consistent with activation of the anti-oxidative enzymes such as catalase, ascorbate peroxidase and guaiacol peroxidase under abiotic stress condition [53]. The activities under BP taxonomy like cellular glucan and cellular amino acid metabolic processes are known to increase in plants in response to various abiotic stresses [54] and other reported biological processes need to be studied in the context of Al toxic stress. The role of phosphotranferase activity in conferring tolerance against abiotic stresses like drought and salt in rice and Arabidopsis are well established [55]. The role of stress induced organic acid synthesis in conferring Al tolerance in higher plants are also well reported [56]. These processes might be related to detoxification of Al ions, which occurs rapidly after exposure to Al stress in plants [57].

dhga R Software package
In order to facilitate the use of proposed hub gene detection and DHGA approaches, we have developed an R software package which includes dhga R package accompanying documentation and model real data examples. This package can be freely downloaded from https:// cran.r-project.org/web/packages/dhga. This software is capable of providing weighted adjacency matrix, edge list, node list for constructing GCN, p-values for gene connections along with weighted gene scores, etc. It also able to identify hub genes and perform differential hub analysis in the GCN based on the proposed approaches. Also the outputs provided by this can be directly used as inputs for gene co-expression networks construction software like Cytoscape, Visant, etc. However, it is difficult to construct GCN using this software package on R-platform.

Conclusions
This investigation has some contributions to gene co-expression analysis. First, a statistically sound Boot-SVM-RFE gene selection technique was proposed for the selection of informative genes from high dimensional gene expression data. Second, a new statistical approach was proposed for the identification of hub genes in a GCN. Third, the proposed DHGA approach may be used to group genes in the GCN into various categories based on their gene connectivity. Fourth, a statistical modelling approach was employed to find the interaction among the gene modules. Moreover, the proposed Boot-SVM-RFE and DHGA approach can be used for other case vs. control genomic studies including NGS expression study. This study also throws some light to understand the mechanism of Al stress response in soybean and some key important genes were reported. Moreover, functional enrichment analysis of these key genes revealed their associated intracellular functions under Al stress. This information revealed in this study on various molecular mechanisms like biosynthesis of secondary metabolites and stress specific roles of certain plant products may be useful for mitigation of Al stress in plants, particularly in soybean. These identified genes can act as potential targets for bio-engineering of Al toxic stress response in soybean.  Table corresponds to modules under individual condition (labeled by color names as well as text along with the number of genes in the modules), and column corresponds to consensus modules. Numbers in the Table indicate