A Kernel-Based Multivariate Feature Selection Method for Microarray Data Classification

High dimensionality and small sample sizes, and their inherent risk of overfitting, pose great challenges for constructing efficient classifiers in microarray data classification. Therefore a feature selection technique should be conducted prior to data classification to enhance prediction performance. In general, filter methods can be considered as principal or auxiliary selection mechanism because of their simplicity, scalability, and low computational complexity. However, a series of trivial examples show that filter methods result in less accurate performance because they ignore the dependencies of features. Although few publications have devoted their attention to reveal the relationship of features by multivariate-based methods, these methods describe relationships among features only by linear methods. While simple linear combination relationship restrict the improvement in performance. In this paper, we used kernel method to discover inherent nonlinear correlations among features as well as between feature and target. Moreover, the number of orthogonal components was determined by kernel Fishers linear discriminant analysis (FLDA) in a self-adaptive manner rather than by manual parameter settings. In order to reveal the effectiveness of our method we performed several experiments and compared the results between our method and other competitive multivariate-based features selectors. In our comparison, we used two classifiers (support vector machine, -nearest neighbor) on two group datasets, namely two-class and multi-class datasets. Experimental results demonstrate that the performance of our method is better than others, especially on three hard-classify datasets, namely Wang's Breast Cancer, Gordon's Lung Adenocarcinoma and Pomeroy's Medulloblastoma.


Introduction
Microarray gene expression based cancer classification is one of the most important tasks in bioinformatics. A typical classification task is to separate healthy patients from cancer patients, based on their gene expression ''profile''. However, because cancers are usually marked by changing in the expression levels of certain genes [1], therefore it is obvious that not all measured features are discriminative features for target. Hence, feature selection problem is ubiquitous in cancer classification.
Feature selection techniques for microarray data can be broadly grouped into three categories that are wrapper (classifierdependent) methods [2,3], embedded (classifier-dependent) methods [4,5] and filter (classifier-independent) methods [6,7]. The primary distinguishing factors among them are computational complexity and the chance of overfitting [8]. Generally, in terms of computational cost, filters are faster than embedded methods, which are in turn faster than wrappers. In terms of overfitting, wrappers have higher learning capacity so are more likely to overfit than embedded methods, which in turn are more likely to overfit than filter methods [9]. Filter methods can be divided into two classes, univariate-based filters and multivariate-based filters [8]. Univariate filter methods have attracted much attention because of their low complexity and fast performance for high dimensionality of microarray data analyses. However, some valuable genes discarded by univariate methods may have great contribution for classification [10]. Therefore, the major reason of their less accurate performance is that they disregard the effects of feature-feature(we use without distinction the term ''feature'' and ''gene'' in the paper) interactions. The applications of multivariate filter methods are simple bivariate-based methods which are almost based on entropy(or conditional entropy) and mutual information [9,11], such as mRMR [7,12], CFS [13] and several variants of the Markov blanket filter method [14]. However, they also abandon presumably redundant variables that can result in a performance loss [15].
Partial least squares(denoted as PLS), which shares the characteristics of other regression and feature transformation techniques(such as canonical correlation analysis and principal component analysis), has proven to be useful in situations when the number of observed variables(D) are significantly greater than the number of observations(N) (e.g.N%D). In other words, PLS is a popular approach to solve problems when there is high multicollinearity among features [16]. SlimPLS [17],PLSRFE [18,19] and TotalPLS [20] are multivariate-based feature selection methods that were proposed by Gutkin et al. and You et al., respectively. Unfortunately, classical PLS technique is essentially a linear regression method that only can capture the linear relationships between genes in original space. In real biological applications, linear relationship often fails to fully capture all the information among genes. Kernel method, which approaches the problem by projecting the data into a high dimensional feature space, is commonly used for revealing the intrinsic relationships that are hidden in the raw data.
Motivated by mentioned above, in this paper, we develop a feature selection method based on the partial least squares(abbreviated PLS) [21] and theory of Reproducing Kernel Hilbert Space [22], we called it kernelPLS(publicly available at https://github. com/sqsun/kernelPLS). Determining the number of components is a thorny problem in PLS(also in kernelPLS) method. In order to obtain a reasonable number of components, we make use of the relationship between PLS and linear discriminant analysis to determine the number of components in kernel space based on kernel linear discriminant analysis. We find that the two classifiers combined with our feature selection method obtained promising classification accuracy on eleven microarray gene expression datasets.
The rest of this paper is organized as follows. In section 2 we proposed a filter method based on PLS and kernel method. Then we proceed in section 3 to determine the optimal parameters for our method. In section 4 we compared our approach with several competitive filters. The conclusion can be found in section 5.

Methods
In the following, let X [ N|D represents a data matrix of N inputs (N samples) and Y [ N|C stands for corresponding response matrix of C-dimensional(C classes). Further we assume columns of X and Y are zero-mean.

Kernel partial least squares
PLS is one of the widespread use of a class of multivariate statistical analysis technique introduced by [21], and a popular regression technique in Chemometrics [23]. It differs from other methods in constructing the fundamental relations between two matrices (X and Y ) by means of latent variables called components, leading to a parsimonious model which shared characteristics with other regression and feature transformation techniques [16]. The goal of PLS is to calculate vectors of its X -weight (v), Y -weight (c), X -score (t) and Y -score (u) by an iterative method for the optimization problem: arg max EvE~1,EcE~1 cov(t,u)~cov(Xv,Yc). Where t~Xv and u~Yc, are called components of X and Y , respectively.
When the first two components t 1 and u 1 are obtained, the second pair t 2 and u 2 is extracted from their residuals E X~X {t 1 p T and E Y~Y {t 1 q T , respectively. Here p and q are called the loadings of t with respect to X and Y , respectively. This process can be repeated until the required halt condition is satisfied. The detail description of the algorithm can be found in [17]. The geometric representation of PLS can be found in Figure 1(a).
The kernel version of PLS uses a nonlinear transformation W( : ) to map the gene expression data into a higher-dimensional(even infinite dimensional) kernel space ; i.e. mapping W : x i [ D ?W(x i )[ . However, we do not need to know the specific mathematical expression of nonlinear mapping, we only need to state the entire algorithm in terms of dot products between pairs of inputs and substitute kernel function K( : , : ) for it. This is so-called the ''kernel trick''.
In order to state dot product operation in the algorithm, we can restrict v to belong to the linear spans of the points. They can therefore be expressed as: be an element of the Gram matrix K X in feature space and h is the desired number of components. Deflating Y will, however, be needed for kernel partial least squares.
The first component for kernel PLS can be determined as eigenvector of the following square kernel matrix for where l is an eigenvalue. The size of the kernel matrix K Y K X is N|N. Hence, no matter how many variables there are in the original matrices X and Y , the size of these kernel matrices will not be get affected by it. Therefore, the combination of PLS with kernel produces a powerful algorithm that will solve this problem rapidly and effectively. The geometric representation of kernel PLS can be found in Figure 1(b). The kernel PLS algorithm procedure and the number of determined components can be found in Table 1 (https://github.com/sqsun/kernelPLS).

The importance of each feature
In original space, let T is a set of components, T~ft 1 ,t 2 , Á Á Á ,t h g. The accumulation of variation explanation of T to Y is given by [24,25] where h is the number of components and v il is the weight of the ith feature for the lth component. Y(Y ,t l )~X C j~1 Y(y j ,t l ) is the correlation between t l and Y , where Y( : , : ) is correlation function. The larger value of w i , the more explanatory power of the ith feature to Y .
It is worth noting that the above equation can also be used in kernel space. The reason is holding of equation W(y j )~y j , because here y j is class label. So the expression Y(W(y j ),t W l ) can be expressed as

Model selection
Two issues are still unresolved before applying kernel PLS for feature selection. The number of components and the number of features are unknow.

The number of components
In order to determine the number of components h, there are two widely used methods in the previous works, one is setting a fixed number, such as h~3, and another is by cross validation (CV). Different datasets contain various data structures, therefore, a fixed number is not suitable for all datasets. Although the CV combined with various classifiers lead to good performance, it suffers from huge computational burden.
To fully circumvent these difficulties, [26] has given an implication of close relationship between PLS and Fisher's linear discriminant analysis (FLDA) in original space. FLDA can be considered as an optimization problem a T S 1 a=a T S 2 a È É , e.g. finding an appropriate projection vector a. Where S 1 presents the inter-class scatter matrix, S 2 denotes the intra-class scatter matrix.
In kernel space, the FLDA turns out to be an optimization , where S W 1 and S W 2 are the inter-class scatter matrix and the intra-class scatter matrix in kernel space, respectively. We consider It denotes the contribution of the lth component for classification. Where N i indicates the number of samples in the ith class, here m W il represents the mean vector of the i th class with respect to l th component in projection space and the c l represents segmentation threshold of classification, the larger c l corresponds to the more significant in classification. Figure 2 shows how classification performance varies with the change in number of features which were selected. The average classification error rate was calculated by two classifiers on all test datasets. An improvement in performance could be evident if the number of related features increase from 1 to 25, but after increasing number of features beyond 25, no significant improvement was obvious. In order to find optimum results for all the datasets, we extend the range from 20 to 50 features configurations in our study.

Test datasets
To assess the performance of our method, we have conducted several experiments on a number of publicly available datasets. Summary of all data sets we used in our experiments can be found in table 2 and following is the brief description of each data set.
N AMLALL(A)( [27]). There are two parts containing the initial (train), 38 bone marrow samples from two classes: 27 cases of acute lymhoblastic leukemia(ALL) and 11 cases of acute myeloid leukemia(AML); independent (test), 34 samples from two classes: 20 cases of ALL and 14 cases of AML. Each case is described by expression levels of 7129 probes from 6817  Table 1. Algorithm 1: kernelPLS.
Output: w -the weight of each feature N Prostate(P) ( [30]). This dataset contains 52 prostate tumor samples and 50 normal samples with 12600 genes. An independent set of testing samples is generated from the training data, 25 tumor and 9 normal samples are extracted

Comparison of selected genes
In our first experiment, we used two datasets, namely the Leukemia data (two-class) of [27] and the Lymphoma data(threeclass) of [33], to compare our method with previous works with respect to the selected genes.
For the Leukemia data, we collected several most important genes (in table 3) that were published in several papers. It can readily be seen that three probes, X95735_at, M27891_at and M23197_at were reported by five published papers, and their ranking by our method are 4th, 17st and 8st, respectively. We notice that there are many overlapping of genes among the list of papers.
For Leukemia data, the top-ranked 40 features obtained by our procedure are shown in table 4 in which genes are in columns from 1 to 40. There is a worthwhile result achieved by our method, that is, it obtained the genes with the highest weight. Many of these genes are known as differentially expressed genes by many foregoing studies. 24 out of 40 genes are listed in this table that were also selected by [27], which shows the effectiveness of our method.
For the Lymphoma data of [33], the missing values are imputed by KNN-imputed method(k~10). The top 40 genes ranked by our procedure are listed in table 5. From the table, We can see that important genes can be captured easily by our method. There are many genes that are also chosen by [38]. Figure 3 illustrates the differentially expressed genes for two datasets, namely the Leukemia data and the Lymphoma data. No single gene is uniformly expressed across the class, all these genes as a group appear correlated with class which is illustrating the effectiveness of the Kernel PLS method. In Figure 3(a) the top panel is consist of three genes GENE1622X, GENE2402X and GENE1648X that are highly expressed in DLCL, middle panel is comprised of GENE1606X, GENE896X and GENE1617X that are highly expressed in DLCL but moderately expressed in FL. Bottom panel compose of three genes, namely GENE1602X,-GENE681X and GENE1618X, are more highly expressed in CLL. In Figure 3(b) the top panel shows three probes highly express in AML and the bottom panel shows three probes more highly expression in ALL. The probe U377055_rna1_s_at was found by our method to distinguish AML from ALL. Figure 3

Comparison of several multivariate-based feature selectors
In our second experiment, we compared several feature selectors with our procedure based on two classifiers, SVM and KNN. In our experiments, we choose the RBF kernel for each dataset to perform classification. To determine the best values of C(-c) and c(-g), we conducted particle swarm optimization algorithm to pick the pair (C,c) with best accuracy in the range of C[f10 {3 , Á Á Á ,10 2 g and c[f10 {3 , Á Á Á ,10 4 g. We set the parameter to k~5 for k-nearest neighbor. To obtain a statistically reliable predictive measurement, we performed 10-fold cross validation for two-class datasets and 5-fold cross validation for multi-class datasets. The results are evaluated by classification accuracy(Acc), area under receiver operating characteristic curve (AUC) for two-class problems and classification accuracy(Acc), Cohen's Kappa coefficient(Kappa) for multi-class problems. The reason of using 5-fold cross validation for multi-class datasets is that there is just a few number of samples in some groups (classes) of these datasets. Therefore to ensure the presence of samples of each class in training and also in test datasets we need to perform 5-fold cross validation for multi-class datasets.
In this paper, the comparison was conducted with four competitive algorithms, PLS, ReliefF, SVMrfe and mRMR. The Table 4. Top-ranked 40 features selected using kernelPLS for the Leukemia dataset. The boldfaced probes were selected by [27]. doi:10.1371/journal.pone.0102541.t004 Table 5. Top-ranked 40 features selected using kernelPLS for the Lymphoma dataset. parameter setting of them are as follows: for the PLS-based feature selection, we used the SIMPLS method and the number of components determined by self-adaptive manner which is the same as the kernelPLS (the proposed method). The parameter k of ReliefF is equal to the number of sample according to the published paper [39]. For SVMrfe, in order to ensure acceptable running time, we use SVM with RBF kernel and its parameter settings are same as in LIBSVM. Without loss of generality, we used two datasets, Breast(twoclass) and Lymphoma(three-class) to show the performance of our method. Figure 4 shows the comparison of error rate between our method and four other methods. One can see that when number of selected features are 30, error rate of our method is less than other methods for both classifiers and both datasets. Table 6 and 7 summarized the comparison of results generated by our method and other methods with respect to Acc and AUC for two-class datasets. From the results, we can see that the performance of our method is better than others. Refers to table 6 we can see that for Breast(B) and Prostate(P) datasets, accuracy of our method is considerably higher as compare to other methods, which shows the effectiveness of our method.
Similarly in table 7 for datasets Breast, Lung, DLBCL, Medulloblastoma, Prostate and Stjude, kernelPLS shown better accuracy rate for SVM classifier wrather than KNN. Both Acc and AUC values of our method have higher values among others and finally the average results likewise are best. Although for few datasets our results are similar to their results but in these cases time taken by our method is significantly smaller than other methods. For example in table 7 for AMLALL dataset, including our method, the AUC is 100% for many methods but time consumed by our method is only 0.0891 s while the time taken by other methods, ReliefF, mRMR, SVMrfe and PLS, are about 5 s, 52 s, 210 s and 12 s, respectively. So time consumption by our algorithm is many times less than others which depicts overall well performance of our method.
It is worth noting that our method outperforms others on three hard-classify datasets, Wang's Breast cancer, Gordon's Lung adenocarcinoma and Pomeroy's Medulloblastoma. We also make a comparison with the results of other feature selectors in published papers. Fox example, the reference [40] reported that the accuracies of k-TSP+SVM on these datasets were 67.1%, 72.2% and 64.2%, respectively. The reference [41] combined multiple feature selection (or feature transform) approaches for Medulloblastoma dataset and the obtained highest Acc was 70%.
To estimate the performance of our method we did not limit our evaluation to only two-class datasets we also used 5 multi-class datasets in our experiments. Tables 8 and 9 demonstrate the comparison of kernelPLS with other methods for multi-class datasets on the bases of results obtained for two evaluation measures, namely Acc and Kappa. Results shown in table 8 and  table 9 are for two classifiers KNN and SVM, respectively. In table 8 results obtained by kernelPLS are better than Relief, SVMrfe and PLS and highly competitive to mRMR method for several multi-class datasets. For example in case of Stjude dataset for Acc and Kappa values by kernelPLS are 96.4% and 0.956 respectively which are highest among all values achieved by other methods. Likewise table 9 authenticates the high performance by kernelPLS over other methods for SVM classifier. Here one can see that kernelPLS give outperforming results for all datasets by achieving accuracies and Kappa coefficients values superior than all other methods. As a conclusion the overall high average Acc and Kappa values in both tables show the effectiveness and significance of our method as compare to other popular methods. Table 10 shows the comparison between running time taken by our method and other methods. There is no single method among these that can perform faster than our method. It clearly shows that kernelPLS is faster than the other algorithms. For example for

Discussion
In this article, we proposed an effective multivariate-based feature filter method for cancer classification, namely, kernelPLSbased filter method. We showed that gene-gene interactions cannot be ignored in feature selection techniques to improve classification performance. In other words the nonlinear relationship of gene-gene interactions is a vital concept that can be taken into account to enhance accuracy. To capture these nonlinear relations of interaction between genes we used kernel method because kernel method can be used to reveal the intrinsic relationships that are hidden in the raw data. In order to capture the reasonable number of components, we make use of the relationship between PLS and linear discriminant analysis to determine the number of components in kernel space based on kernel linear discriminant analysis. To verify the importance of gene-gene interactions we compared our feature selector with other multivariate-based feature selection methods by using two classifiers SVM and KNN. Experimental results, expressed as both accuracy(Acc) and area under the ROC curve(AUC), showed that our method leads to promising improvement in ACC and AUC. We can conclude that the gene-gene interactions whats more, nonlinear relationships of gene-gene interactions are core interactions that can improve classification accuracy, efficiently. We can summarize the characteristics of proposed approach as follows: (1)Fast and efficient. The time complexity of deflation procedure used after the extraction of each component scale is O (N 2 ), where N is the number of sample. In most cases, the number of sample in microarray data is less than 150, therefore, the running speed of kernelPLS procedure(feature selection time) is faster than others, which are summarized in table 10. (2)Modelfree, e.g. no need the distributional assumptions. Because of small sample size, it is difficult to validate distributional assumptions, such as Gaussian distribution, Gamma distribution etc. (3)Applicable to both two-class as well as multi-class classification problems.
In our method, the choice of kernel functions can affect the results. When high dimensionality exist(such as microarray datasets), the performance of linear kernel is better than Gauss   Table 8.  Table 9. Comparison of kernelPLS with four other methods. For 5-fold cross validation classification accuracy(%) and Cohen's kappa coefficient of SVM on multi-class datasets. kernel for our method. What's more, in case of linear kernel there is no noticeable effect on the results while adjusting its parameters.