A Hybrid One-Way ANOVA Approach for the Robust and Efficient Estimation of Differential Gene Expression with Multiple Patterns

Background Identifying genes that are differentially expressed (DE) between two or more conditions with multiple patterns of expression is one of the primary objectives of gene expression data analysis. Several statistical approaches, including one-way analysis of variance (ANOVA), are used to identify DE genes. However, most of these methods provide misleading results for two or more conditions with multiple patterns of expression in the presence of outlying genes. In this paper, an attempt is made to develop a hybrid one-way ANOVA approach that unifies the robustness and efficiency of estimation using the minimum β-divergence method to overcome some problems that arise in the existing robust methods for both small- and large-sample cases with multiple patterns of expression. Results The proposed method relies on a β-weight function, which produces values between 0 and 1. The β-weight function with β = 0.2 is used as a measure of outlier detection. It assigns smaller weights (≥ 0) to outlying expressions and larger weights (≤ 1) to typical expressions. The distribution of the β-weights is used to calculate the cut-off point, which is compared to the observed β-weight of an expression to determine whether that gene expression is an outlier. This weight function plays a key role in unifying the robustness and efficiency of estimation in one-way ANOVA. Conclusion Analyses of simulated gene expression profiles revealed that all eight methods (ANOVA, SAM, LIMMA, EBarrays, eLNN, KW, robust BetaEB and proposed) perform almost identically for m = 2 conditions in the absence of outliers. However, the robust BetaEB method and the proposed method exhibited considerably better performance than the other six methods in the presence of outliers. In this case, the BetaEB method exhibited slightly better performance than the proposed method for the small-sample cases, but the the proposed method exhibited much better performance than the BetaEB method for both the small- and large-sample cases in the presence of more than 50% outlying genes. The proposed method also exhibited better performance than the other methods for m > 2 conditions with multiple patterns of expression, where the BetaEB was not extended for this condition. Therefore, the proposed approach would be more suitable and reliable on average for the identification of DE genes between two or more conditions with multiple patterns of expression.

and the proposed method exhibited considerably better performance than the other six methods in the presence of outliers. In this case, the BetaEB method exhibited slightly better performance than the proposed method for the small-sample cases, but the the proposed method exhibited much better performance than the BetaEB method for both the small-and large-sample cases in the presence of more than 50% outlying genes. The proposed method also exhibited better performance than the other methods for m > 2 conditions with multiple patterns of expression, where the BetaEB was not extended for this condition. Therefore, the proposed approach would be more suitable and reliable on average for the identification of DE genes between two or more conditions with multiple patterns of expression.

Introduction
Microarray technology has enabled the expression levels of thousands of genes to be investigated simultaneously. However, this technology poses statistical challenges by virtue of the large number of transcripts surveyed with small sample sizes. The identification of transcripts that are differentially expressed (DE) between two or more conditions is a common task that is undertaken to reduce the dimensionality of the transcripts, as important genes belong to the reduced set of DE transcripts. Useful information regarding the regulatory network can be obtained by associating differential expressions with the genotypes of molecular markers [1]. By assigning DE genes to the list of gene sets, it is possible to obtain useful biological interpretations [2,3]. Furthermore, the number of DE genes that influence a certain phenotype may be large, whereas their relative proportions are typically small; therefore, identifying these DE genes from among the large number of recorded genes is challenging [4][5][6][7][8].
In general, four types of statistical procedures are used to identify DE genes: (i) classical parametric approaches, such as t-test, Ftest (ANOVA) and likelihood ratio test (LRT)-based asymptotic χ 2 -test; (ii) classical nonparametric approaches [8][9][10][11]; (iii) empirical Bayes (EB) parametric approaches [5-7, 12, 13]; and (iv) EB nonparametric approaches [14,15]. In the classical procedures, DE genes are generally detected based on p-values (significance levels) that are estimated either via permutation or based on the distribution of a test statistic, whereas in EB procedures, the posterior probability (pp) of differential expression is used to identify DE genes. However, most of the aforementioned algorithms are not robust against outliers [12,16]. Thus, they may produce misleading results in the presence of outlying transcripts or irregular patterns of expression. Several recent studies have reported that the assumption of normality does not hold for some existing microarray datasets [17]. One of the causes for this breakdown of the normality assumption may be related to the presence of outliers in the data. cDNA microarray data are often corrupted by outliers that arise because of the many steps that are involved in the experimental process, from hybridization to image analysis [12,16].
Some nonparametric approaches, such as KW [9], are somewhat robust against outliers between two or more conditions with multiple patterns of expression in the case of large sample sizes; however, these approaches are sensitive to outliers in the case of small sample sizes. To overcome this problem, the β-divergence-based empirical Bayes (BetaEB) approach [16] was developed for the robust identification of DE genes. This approach performs well in the presence of outlying expressions with up to 50% genes for both small-and large-sample cases. The parameters of the BetaEB approach are estimated based on the expressions of all genes. It can tolerate up to 50% outlying genes if the mean vector is initialized by the median vector. In the presence of more than 50% outlying genes, it is difficult to initialize the shifting parameters (mean vector) of this approach with the good part of the dataset. Therefore, the BetaEB approach occasionally produces misleading results. Note that more than 50% outlying genes may also occasionally occur with at least one patient/tissue sample. Moreover, this approach was not extended for the detection of DE genes in the case of more than two conditions/groups with multiple patterns of expression due to the computational complexity. Therefore, in this paper, an attempt is made to develop a hybrid approach that unifies the robustness and efficiency of estimation in one-way ANOVA using the minimum β-divergence method [18][19][20][21] to overcome all of the aforementioned problems that arise in the BetaEB approach. The advantage of the proposed algorithm compared to BetaEB is that it performs considerably better than the BetaEB approach for both the small-and large-sample cases in the presence of more than 50% outlying genes with 5% outlier samples for each outlying gene. The major advantage of the proposed algorithm is that it performs well in the case of multiple conditions/groups (m > 2) for identifying DE genes with multiple patterns of expression, whereas BetaEB was not extended for multiple conditions/groups (m > 2) due to the computational complexity. The proposed method introduces a weight function, which plays a key role in its performance. This weight function assigns smaller weights to outlying observations, thereby ensuring the robustness of the inference. Appropriate initialization of the parameters also improves the performance of the proposed method, as is also discussed in this paper.

Materials and Methods
Let x jk be the kth observed random expression of a gene in the jth condition (j = 1, 2, . . ., m; k = 1, 2, . . ., n j ), which follows the one-way ANOVA model as expressed below: where μ j is the mean of all expressions of a gene in the jth condition and jk is the random error term that follows Nð0; s 2 j Þ. We wish to test the null hypothesis (H 0 ) : μ 1 = μ 2 = . . . = μ m = μ against the alternative hypothesis (H 1 ) : H 0 is not true, assuming that s 2 1 ¼ s 2 2 ¼ ::: ¼ s 2 m ¼ s 2 . Thus, the generalized likelihood-ratio test (LRT) criterion yields the following F-statistic to test H 0 against H 1 : which follows the F-distribution with (m-1) and (n-m) degrees of freedom under H 0 [22], where n = n 1 + n 2 +, . . ., + n m andm ¼ P m j¼1 n jmj =n. Here,m j andŝ 2 j are the maximum likelihood estimates (MLEs) of μ j and s 2 j , respectively, for the jth condition/group. The critical region (CR) for testing H 0 against H 1 at the (1-α)100% level of significance is defined by is the upper 100α% points of the Fdistribution with m − 1 and n − m degrees of freedom. F 0 is also known as the cut-off point or critical value of the test. However, it is well known that the MLEsθ j ¼ ðm j ;ŝ 2 j Þ of y j ¼ ðm j ; s 2 j Þ for j = 1, 2, . . ., m in Eq (2) are highly sensitive to outliers. Therefore, the identification of DE genes using classical ANOVA may produce misleading results because gene expression data are often corrupted by outliers, as previously discussed. Thus, in this paper, we consider the minimum β-divergence method [20,21] to improve the robustness and efficiency of estimation in one-way ANOVA. The minimum β-divergence estimatorsθ j;b ¼ ðm j;b ;ŝ 2 j;b Þ of the parameters θ j ¼ ðm j ; s 2 j Þ are computed iteratively as follows: and where which we call the β-weight function [20,21]. The notation θ t+1 represents the update to θ t in the (t+1)th iteration. The robustness of these estimators is discussed in the context of influence functions in [20], and their consistency is discussed in [21]. Note that the minimum β-divergence estimatorsθ j;b ¼ ðm j;b ;ŝ 2 j;b Þ reduce to the classical MLEsθ j ¼ ðm j ;ŝ 2 j Þ for β = 0. Note that MLEs of a Gaussian distribution are consistent and asymptotically efficient in the absence of outliers [23]. Therefore, in this paper, an attempt is made to develop a hybrid approach in which the classical MLEsθ j are used in the absence of outliers and the minimum β-divergence estimatorsθ j;b (Eqs (3) and (4) are used in the presence of outliers for the estimation of θ j in one-way ANOVA. The minimum β-divergence method offers two approaches for unifying the robustness and efficiency of estimation in ANOVA. One method is to select the tuning parameter β via cross-validation (CV), as discussed in detail in a previous publication [20]. In the absence of outliers, the CV method produces β = 0 for the minimum β-divergence estimators and is thus equivalent to the classical estimators, as discussed above. In the presence of outliers, it produces β > 0 for the minimum β-divergence estimators. To develop the alternative approach, we consider the β-weight function (Eq (5)) with β = 0.2 for outlier detection. This weight function assigns smaller weights (! 0) to outlying/unusual observations and larger weights ( 1) to uncorrupted/usual observations. An outlying gene expression x jk is defined based on the β-weight function as follows: where the threshold value δ j is the pth quantile value of the distribution of W b ðx jk jθ j;b Þ. The predicted distribution of the β-weight function is discussed in S1 File.
Thus, we can unify the minimum β-divergence estimator with MLEs for θ j in the jth condition as follows:θ which leads to the formulation of the hybrid one-way ANOVA approach for the robust and efficient estimation of differential gene expression with multiple patterns with a modified F-statistic, denoted by F β , as given by To test the null hypothesis (H 0 ) against the alternative hypothesis (H 1 ) from the robustness perspective, we can compute p-values under the assumption that F β approximately follows the F-distribution. Note that this modified F-statistic (F β ) reduces to the classical F-statistic (Eq (2)) for β = 0. However, we can also compute permutation-based p-values to test whether H 0 is true or false, as discussed in S1 File.

Robust Multiple Comparison Test
Multiple comparison procedures are commonly used in ANOVA after a significant omnibus test result is obtained via the F-test. A significant ANOVA result suggests the rejection of the global null hypothesis H 0 , which states that the means are identical across all groups being compared. Multiple comparison procedures are then used to determine which of the means differ significantly. In a one-way ANOVA involving m group means, there are m(m − 1)/2 pairwise comparisons. There are several methods for performing such a multiple comparison test.
In this paper, we consider the robustification of Fisher's least significant difference (LSD) test for multiple comparison tests. We reject H 0 : μ i = μ j at the 100(1-α)% level of significance for all i 6 ¼ j when jm i;b Àm j;b j> LSD b . Here, LSD β is computed as follows: where MSE b ¼ ½n 1ŝ

Results and Discussion
We investigated the performance of the proposed method in comparison with other popular methods using both simulated and real gene expression data.

Performance Evaluation Based on Simulated Gene Expression Profiles
To investigate the performance of the proposed method in comparison with several popular existing methods, such as the classical parametric approach ANOVA (F-test), the nonparametric approaches SAM [10] and KW (Kruskal-Wallis test) and the empirical Bayes (EB) approaches LIMMA [13], EBarrays [5,24], eLNN [25] and BetaEB [16], in the detection of DE (important) genes, we considered gene expression profiles simulated based on classical (ANOVA) and Bayesian (EBarrays LNN) data generation models for both small-and largesample cases with two and multiple groups/conditions in both the absence and the presence of outlying expression, as discussed in subsections 3.1.1 and 3.1.2.
3.1.1 Performance Evaluation Based on Simulated Gene Expression Profiles with m = 2 Conditions. To investigate the performance of the proposed method in comparison with seven popular methods (ANOVA, SAM, eLNN, LIMMA, KW, EBarrays and BetaEB), as mentioned above, for both small-and large-sample cases with m = 2 groups/conditions, we considered 100 datasets for both cases with sample sizes of n 1 = n 2 = 3,4 and 15, respectively. Each dataset for each case represented gene expression profiles for 20,000 genes, each with n = (n 1 + n 2 ) samples, where the expression of each gene was generated using the one-way ANOVA model described by Eq (1) for arbitrary values of (μ 1 , μ 2 ) 2 (2,5) and σ 2 = 0.05. Among the 20,000 genes represented in each dataset for each case, we generated 300 DE genes (μ 1 6 ¼ μ 2 ) and 19,700 EE (equally expressed) genes (μ 1 = μ 2 ). To investigate the robustness performance, we generated three types of outlying datasets from each of the original datasets by replacing at most 5% of the expression values from (n 1 + n 2 ) samples by outliers with each 5%, 10% and 75% genes, respectively. The kth expression in the jth group for a gene was corrupted by an outlier generated as follows: x Ã jk ¼ dþ max(x jk ;k = 1, 2, . . ., n j ;j = 1, 2),, where d 2 (5,10) is an arbitrary value.
We first investigated the goodness of fit of the predicted distribution of the β weights (Eq 12 in S1 File) compared with the observed distribution of β weights (W b ðÁ jθ j;b Þ). We produced both the predicted distribution (solid curve) and the observed distribution based on simulated gene expression data, assuming that m j ¼m j;b and s 2 j ¼ŝ 2 j;b (Eq (1)) in both the absence and the presence of 5% or 10% outlying expressions, as shown in Fig 1a and 1b, respectively. It is clear that the predicted distributions (solid curves) in both figures reflect the corresponding observed distributions (histograms). In Fig 1b, smaller β weights or insignificant weights with p-values < 10 −5 are evident, which correspond to outlying expressions. Note that the predicted distributions (solid curves) are identical in both figures. Therefore, we can use either of these two distributions to select the cut-off point (Eq (6)) to determine whether an expression is an outlier. Table 1 and Table A in S1 File summarize the average performance results of the eight investigated methods (ANOVA, SAM, LIMMA, eLNN, EBarrays, BetaEB, KW and Proposed) based on 100 datasets generated using a one-way ANOVA model with m = 2 groups/conditions and σ 2 = 0.05 for cases of both small (n 1 = n 2 = 3, 4) and large (n 1 = n 2 = 15) samples. The performance indices/measures, including the true positive rate (TPR), the false positive rate (FPR), the true negative rate (TNR), the false negative rate (FNR), the false discovery rate (FDR), the misclassification error rate (MER), the area under the ROC curve (AUC) and the partial AUC (pAUC), were calculated for each method based on their estimated top 300 DE genes, under the assumption that the other estimated genes for each method were the EE genes for each dataset (recall that each dataset contained 300 true DE genes and that the remainder were true EE genes). The values of each of these performance measures lie between 0 and 1. The method that produces the smallest values of FPR, FNR, FDR and MER and the largest values of TPR, TNR, pAUC and AUC is considered to be the best performer. The results presented in Table 1, Table A in S1 File and S1, S2, and S3 Figs clearly show the performance of all eight methods for the identification of DE genes in the absence and presence of 5% or 10% outlying expressions with 5%, 10% and 75% genes for both the small-and large-sample cases. We clearly observe that all eight methods performed similarly in the absence of outlying genes for both the small-and large-sample cases. An interesting coincidence occurs between FNR and FDR for all methods in all cases because we declared the top (1 − p 0 ) × 100% tp be DE genes, where p 0 = 0.985 is the proportion of EE genes. A similar situation was also observed in [26]. For the same reason, we also obtained almost identical results for FPR and TNR in all cases for all methods.
However, in the presence of 1 or 2 outlying expressions with 5%, 10% and 75% genes for the small-or large-sample or both sample cases, four methods (eLNN, KW, BetaEB and Proposed) exhibited better performance on average than the other four methods (ANOVA, SAM, EBarrays and LIMMA) because the former four methods produce larger values of TPR, AUC and pAUC and smaller values of FNR, FDR and MER compared to the other four methods. Two methods (KW and eLNN) showed slightly good performance for the large-sample case, but they were sensitive to outliers for the small-sample case. The empirical Bayes approach (BetaEB) exhibited good performance for both the small-and large-sample cases in the presence of outlying expressions with 5% and 10% genes, but it exhibited weak performance in the presence of outlying expressions in 75% genes. Thus, the robust BetaEB method and the proposed method exhibited better performance than the other six methods in presence of at most 50% outlying genes (although we provided results with 5%, 10% and 75% outlying genes only) with 5% outliers for each outlying gene. The BetaEB approach appears to be slightly better than the proposed method for the small-sample case, whereas the proposed method exhibited considerably better performance than the BetaEB method for both the small-and large-sample cases in the presence of more than 50% outlying genes. We also observed similar results when we generated gene expression profiles based on a Bayesian (EBarrays LNN) data generating model considering all aspects of the ANOVA model (Table B in S1 File).
In this study, we investigated the performance of the proposed method in comparison with only four popular methods (ANOVA, SAM, LIMMA and KW) because the robust BetaEB approach is not applicable for multiple patterns of expression because of its computational complexity. Table 2 shows that all the methods, except KW, produce almost identical values of FDR, AUC and pAUC in the absence of outlying genes for the small-sample case. KW was found to have considerably lower values of AUC and pAUC and a larger value of FDR for multiple groups/conditions (m = 4) in the small-sample case. Additionally, it exhibits greater sensitivity than the other methods (ANOVA, SAM, LIMMA and proposed) in the presence of outlying genes. However, all of the methods, including KW, perform similarly in the absence of outlying genes for the large-sample case. In this case, the performance of KW is better than that of the other three methods (ANOVA, SAM and LIMMA) in the presence of outlying genes, but it is slightly weaker than that of the proposed method. Thus, in all cases, the proposed method exhibited performance similar to that of the other approaches in the absence of outliers and better performance in the presence of outliers with lower FDR values and higher AUC and pAUC values. We also investigated the performance of the proposed method in       Interestingly, we observed that all methods, including SAM and KW, perform well for the large-sample case in the absence of outlying genes, whereas we observed a loss of efficiency of SAM and KW for the small-sample case. The nonparametric approach KW is more efficient than the other three methods (ANOVA, SAM and LIMMA) in the presence of outliers in different proportions of genes for the large-sample case (Table 4). Therefore, we observed that the proposed method performed similarly in the absence of outliers, but it exhibited the best  performance in the presence of outliers compared to all other methods in all cases of both sample sizes (Tables 3 and 4). The results presented in Tables 2, 3 and 4 were calculated under the global null hypothesis (H 0 ) of no differential gene expressions in the m = 4 conditions. However, after rejection of the global H 0 , a multiple comparison test was required to identify the pattern of differential gene expression. Tables 5 and 6 present the adjusted p-values of the multiple comparison tests for ANOVA, LIMMA, KW and the proposed method based on the datasets corresponding to the true patterns μ 1 = μ 2 6 ¼ μ 3 = μ 4 and μ 1 = μ 2 6 ¼ μ 3 6 ¼ μ 4 , respectively. Both tables indicate that all four methods (ANOVA, LIMMA, KW and Proposed) were successful in identifying the correct patterns of differential gene expression in the absence of outlying expressions, whereas in the  presence of a few outliers (5%), the proposed method exhibited superior performance in identifying the correct patterns. For example, in the absence of outlying genes, all four methods (ANOVA, LIMMA, KW and Proposed) provided larger p-values, such as 0.9379, 1.000, 0.8808 and 0.5537, respectively (Table 5), because we considered no difference between group 1 and group 2 in the pattern of a gene. Among them, the LIMMA performed better (p-value = 1) than the other approaches in obtaining the true difference between group 1 and group 2 of the pattern. In the case of group 1 and group 3, all four methods detected the true difference between group 1 and group 3 of the pattern, but among them, the ANOVA, LIMMA and proposed methods performed better than the KW method in obtaining the true pattern. However, two (ANOVA and LIMMA) of the four methods failed to detect the true patterns in the  presence of 5% outliers in the previous expression, as shown in Tables 5 and 6, respectively. The nonparametric approach KW performed more poorly than these two methods in detecting the true patterns with larger p-values, even in the absence of outlying genes for the small-sample case, but better for the large-sample case with multiple groups in both the absence and presence of outlying genes, whereas the proposed method remained similar, as previously mentioned, to the case where outliers were absent in both sample sizes.
To investigate the pattern-detection performance of the proposed method in comparison to the others, we generated 300 DE genes among the 20,000 genes in the dataset for m = 4 conditions with different patterns for the sample size (n1 = n2 = n3 = n4 = 6) with a 2-fold change in expression between the groups. Then, we first applied four methods (ANOVA, LIMMA, KW and Proposed) to test the null hypothesis of no differential expression among four groups for each gene. If the test was rejected, then we applied the multiple comparison test for pattern detection of gene expression. Table 7 presents the multiple comparison results, where the values presented in the form {x, x, x, x} indicate the numbers of downregulated (DR) or upregulated (UR) differentially expressed (DE) genes estimated using ANOVA, LIMMA, KW and the proposed methods, respectively. This table also describes the true patterns of the DE genes in the various groups. We observe that all four methods exhibited nearly identical performance in the absence of outlying genes; however, in the presence of a single outlying expression in each of 10% genes, the proposed method exhibited far superior performance compared to the other methods.

Performance Evaluation Based on Real Gene Expression Profiles
We considered three publicly available microarray gene expression datasets in the analyses presented in subsections 3.2.1, 3.2.2 and 3.2.3 to evaluate the performance of the proposed method in comparison with several popular methods, as discussed above. These three datasets are (i) the platinum spike gene expression dataset [27], (ii) the colon cancer gene expression dataset [28] and (iii) the pancreatic cancer gene expression dataset [29]. All three datasets were generated using Affymetrix technology. 3.2.1 Analysis of the Platinum Spike Gene Expression Dataset with m = 2 Conditions. This dataset was previously analyzed in [4,27]. It consists of 18 spike-in samples (9 controls versus 9 test cases). We downloaded this dataset from the GEO website under the accession number GSE21344. We also downloaded the designated fold change (FC) dataset associated with the probes from www.biomedcentral.com/content/supplementary/1471-2105-11-285-s5. txt. After pre-processing (using RMA) and filtering the dataset, we obtained gene expressions with 18707 probes, among which 1944 probes are known as the designated DE genes under spiked-in fold changes of 0.25, 0.28, 0.40, 0.66, 0.83, 1.5, 1.7, 2, 3 and 3.5. In our analysis, we consider these designated DE genes as designated 'DE gene-set' and the rest of the genes as designated 'EE gene-set' for performance evaluation of the proposed method in a comparison with the other seven methods (ANOVA, SAM, eLNN, LIMMA, KW, EBarrays, BetaEB). We applied all eight methods to the dataset to identify the DE genes. We considered the estimated top 1944 genes for each method and crossed with the designated 'DE gene-set' to calculate the same summary statistics (TPR, TNR, FPR, FNR, FDR, MER, AUC and pAUC) for performance evaluation as used for the simulation studies. Table 8 shows that all the methods produce similar results with the original spike dataset. Although the performance of all methods appears to be Table 8. Performance evaluation based on Spike gene expression profiles with 2 conditions for the sample case (n 1 = n 2 = 9).  In the M-A plot, the marker (Á) with a gray color represents all genes, the marker (°) with a green color represents designated DE genes, the marker ( Ã ) with a blue color represents designated EE outlying genes, and the marker ( Ã ) with a red color represents designated DE outlying genes. We clearly observed that some designated DE genes belong to the zero line (EE line); however, no designated EE outlying genes belong to the zero line (EE line). Therefore, the outlying 233 designated EE genes detected by the proposed method as DE genes with p-values of less than 0.05 should belong to the designated DE geneset. Thus, the robust methods (KW, BetaEB and Proposed) did not perform better than the classical methods (ANOVA, SAM, LIMMA, eLNN and EBarrays), even for the dataset containing outlying genes. To show the outlier effects on the methods, we corrupted low label expressions of 200 designated DE genes by a single outlier larger than the maximum value of expressions. Then, we applied all methods again and evaluated the performance, as shown in Table 8. We found that the proposed and BetaEB methods remained almost unchanged for all the summary statistics, as previously calculated in the case of the original expression, whereas all of the other methods decreased TPR, AUC, and pAUC and increased FNR, FDR and MER significantly. Thus, the application of robust methods would be better than classical methods for real gene expression data analysis.

Analysis of the Colon Cancer Gene Expression
Dataset with m = 2 Conditions. This dataset has been analyzed in a previous study [28] and consists of 22 control and 40 colon cancer samples with 2000 genes. We downloaded this dataset from the website http:// microarray.princeton.edu/oncology/affydata/index.html. We then directly applied three methods (KW, BetaEB and the proposed method) to this dataset, as before, to identify the DE genes. To obtain some insights into the possible mechanisms that may be important in the development of resistance to drugs, the WebGestalt2 software package (available from http://bioinfo.vanderbilt.edu/webgestalt) [30] was used to query two pathway databases, KEGG (Kyoto Encyclopedia of Genes and Genomes) and GO (Gene Ontology) Pathways, for these 8 DE genes, yielding many important pathways to be enriched. We identified three important genes, which were described as L21993 (ADCY2; adenylate cyclase 2 (brain)), U21090 (POLD2; polymerase (DNA directed) delta 2, accessory subunit) and X74330 (PRIM1; primase DNA, polypeptide 1 (49 kDa)), which are involved in purine metabolism, DNA replication, pyrimidine metabolism and other metabolic pathways.
Using the GO database, we determined that POLD2 and PRIM1 are involved in biological processes and that the genes POLD2, PRIM1, U29092 (ubiquitin-conjugating enzyme E2I  Table C in S1 File, various independent studies have reported significant gene under-/over-expression in colon cancer tissues compared with corresponding normal tissues. We found that four independent studies have reported significant over-expression of PRIM1 (p-value 0.0002) in colon cancer tissues compared with corresponding normal tissues. Moreover, in one of these four studies, the rank of PRIM1 overexpression was within the top 4%, whereas in the remaining studies, this over-expression was ranked in the top 6-21%, similar to other genes (Table C in S1 File). These studies strongly suggest that these 8 genes may also be associated with colon cancer.
In our analysis, we also found some DE genes detected by other methods, such as only 5 genes by KW and 10 genes by BetaEB (Fig 2a). We applied all of these 15 genes to WebGestalt2 software for KEGG pathway analysis. Only five genes, X53743 (FBLN1: fibulin 1), X67699 (CD52: CD52 molecule), D38551 (RAD21: RAD21 homolog (S. pombe)), X16356 (CEA-CAM1: carcinoembryonic antigen-related cell adhesion molecule 1 (biliary glycoprotein)) and X07290 (ZNF3: zinc finger protein 3), were mapped using the KEGG Map in WebGestalt. We did not obtain any specific pathway for these genes, but some genes were very crucial from the previous investigations by different authors. As an example, the gene FBLN1 (fibulin 1) (detected by the KW approach only) was identified as a tumor suppressor gene whose inactivation may contribute to carcinogenesis [31]. It plays a tumor suppressive role in colorectal cancer (CRC) [32] and gastric cancer [33], among others. The gene RAD21 (detected only by the BetaEB method) is a component of the cohesion complex and is integral to chromosome segregation and error-free DNA repair. This gene expression in CRC is associated with aggressive disease, particularly in KRAS mutant tumors, and resistance to chemoradiotherapy. The gene RAD21 may be an important novel therapeutic target [34]. The gene CEACAM1 (detected only by the BetaEB method) is a tumor suppressor whose expression is known to be lost in the great majority of early adenomas and carcinomas. The loss of CEACAM1 expression is more common in neoplastic tumors than the adenomatous polyposis coli (APC) mutations [35]. Notably, these genes were also statistically significant by the proposed method with p-values (p) of 0.0003 < p < 0.0020, but they were not included in the top 100 genes. On the other hand, 8 DE genes of the proposed method were found by the KW and BetaEB approaches with p-values of 0.004 < p < 0.02 and posterior probabilities (pp) of 0.50 < pp < 1, respectively.
3.2.3 Analysis of the Pancreatic Cancer Gene Expression Dataset with m = 4 Conditions. This dataset was used in the research reported in [29]. It consists of 24 samples, with 6 replicates in each of m = 4 conditions -circulation tumor samples (CTC), hematological cells (G), original tumor tissue (T), and non-tumor pancreatic control tissue (P) from patients and represents 8152 genes. We downloaded this dataset from the website http://www.ncbi.nlm.nih. gov/geo/query/acc.cgi?token = zbotduky//wssgujm/&acc=GSE18670. This dataset had been pre-processed by previous users [29]. Therefore, we directly applied four methods (ANOVA, LIMMA, KW and Proposed) to this dataset to test the null hypothesis (H 0 ) of no differential expressions among four groups for each gene. If the H 0 was rejected, we performed a pairwise comparison test for all 4 methods to identify the pattern of gene expressions . Fig 3a presents a Venn diagram of the DE genes estimated by the four methods based on pairwise comparisons of CTC vs T, CTC vs P, CTC vs G, T vs P, T vs G and G vs P, which was also the approach used in [29]. We used the Bonferroni method to adjust the p-values for ANOVA and KW and the Benjamini-Hochberg (BH) method to adjust the p-values for LIMMA and the proposed method. We used adjusted p-values at the 5% level of significance with an absolute fold change of 2 for each pairwise comparison in identifying DE genes. From this Venn diagram, it is be involved in cancer cell migration. In the p38 MAPK pathway, PP2CB and MEF2C were significantly upregulated. The simulation study results in Table 7 also increase our confidence in the results of the proposed method for this real dataset because that simulated dataset was generated in accordance with the parametric properties of this real dataset.

Conclusion
We proposed a hybrid one-way ANOVA approach that unifies the robustness and efficiency of estimation of model parameters for the discovery of differential gene expressions with two or more conditions with multiple patterns. The proposed approach is controlled by the β-weight function of the minimum β-divergence method such that MLEs are used in the absence of outliers and minimum β-divergence estimators are used in the presence of outliers for estimating the group parameters in the ANOVA model. The proposed method produces robust and efficient results because MLEs are consistent and asymptotically efficient under a Gaussian distribution in the absence of outliers and because the minimum β-divergence estimators are highly robust and asymptotically efficient in the presence of outliers. It overcomes the problems that arise in the existing robust methods for both small-and large-sample cases with multiple patterns of gene expression. The β-weight function plays the key role in the performance of the proposed method. It can accurately and significantly detect outlying expressions. The simulation results showed that all eight methods (ANOVA, SAM, LIMMA, EBarrays, eLNN, KW, robust BetaEB and proposed) perform almost identically for both the small-and large-sample cases with m = 2 conditions in the absence of outliers. The robust BetaEB method and the proposed method exhibited better performance than the other six methods in the presence of at most 50% outlying genes (although we provided results with 5%, 10% and 75% outlying genes only) with one or two outliers for each outlying gene. However, the BetaEB approach appears to be slightly better than the proposed method for the small-sample case, whereas the proposed method exhibited considerably better performance than the BetaEB method for both the small-and large-sample cases in the presence of more than 50% outlying genes with 5% outlier samples for each outlying gene. To investigate the performance of the proposed method in the case of multiple (m > 2) conditions with multiple patterns of expression, four popular methods (ANOVA, SAM, LIMMA and KW) were considered based on the availability of software/extended-version of the algorithms for testing the equality multiple means. The nonparametric approach KW exhibited weak performance in a comparison of four methods (ANOVA, SAM, LIMMA and proposed) for the small-sample case with m > 2 conditions in the absence of outliers. However, the proposed method exhibited improved performance compared to the other methods for the small-sample case in the presence of outliers. The nonparametric KW method and the proposed method showed better performance in a comparison of the other three methods (ANOVA, SAM and LIMMA) for the large-sample case with m > 2 conditions in the presence of few outlying genes, where the proposed method appears to be slightly better than the KW method. The results on real gene expression datasets demonstrated that the proposed method can identify some additional outlying genes as differential expression compared to the other methods. Moreover, we found that the identified genes are reported as being important genes in different studies by other researchers. Therefore, the proposed approach would be more suitable and reliable on average for the identification of DE genes between two or more conditions with multiple patterns of expression.  (8) genes. This directed acyclic graph (DAG) shows the gene ontology categories of eight (8) genes, detected by the proposed method only in the colon cancer dataset, obtained using the WebGestalt database. These enriched GO categories were hierarchically organized into a DAG tree; each box in the tree lists the name of the GO category, the number of genes in that category, and the FDR-adjusted p-value (adjP) if the enrichment is significant. The categories shown in red are enriched (adjusted p-value of < 0.05), whereas those in black are non-enriched. (TIF) S1 File. Results for simulated and real gene expression datasets. Performance evaluations based on simulated gene expression profiles are presented in Tables A and B and real gene expression colon cancer dataset results are presented in Table C