Characteristic Gene Selection via Weighting Principal Components by Singular Values

Conventional gene selection methods based on principal component analysis (PCA) use only the first principal component (PC) of PCA or sparse PCA to select characteristic genes. These methods indeed assume that the first PC plays a dominant role in gene selection. However, in a number of cases this assumption is not satisfied, so the conventional PCA-based methods usually provide poor selection results. In order to improve the performance of the PCA-based gene selection method, we put forward the gene selection method via weighting PCs by singular values (WPCS). Because different PCs have different importance, the singular values are exploited as the weights to represent the influence on gene selection of different PCs. The ROC curves and AUC statistics on artificial data show that our method outperforms the state-of-the-art methods. Moreover, experimental results on real gene expression data sets show that our method can extract more characteristic genes in response to abiotic stresses than conventional gene selection methods.


Introduction
The growth of plants is greatly affected by a variety of abiotic stresses, such as cold, drought, salt, heat, UV-B light, osmotic press, and so on. In response to these abiotic stresses, plants have evolved a number of defense mechanisms that can increase tolerance to the adverse conditions. The underlying concept is that there exists a specific set of interacting genes responding to the abiotic stresses. Therefore, understanding abiotic stress responses is now thought to be one of the most important topics in plant science [1].
In order to obtain the characteristic genes responding to these stresses, many conventional experimental methods were proposed, such as RT-PCR [2,3] and Northern blotting [4,5], etc. RT-PCR can accurately position the genes in tissue or cell and Northern blotting can display the information of detected genes. However, these two methods have the following fatal flaw: only a limited amount of genes can be simultaneously studied. To overcome this flaw, people have developed the gene microarray technology [6][7][8], which makes it possible to monitor gene expression levels on a genomic scale [9,10].
With the rapidly development of gene microarray technology, how to efficiently analyze gene expression data becomes a matter of great urgency. During the last decade, feature selection from gene expression data has been extensively studied. The most commonly used methods of feature selection first calculate a score for each gene, respectively, then select the genes with high scores [11]. These methods are often denoted as univariate feature selections (UFS). The main virtues of UFS are: (a) intuitive and easy to understand; (b) computationally simple; and (c) fast [12]. A common disadvantage of UFS based methods is that each feature is separately considered, thereby ignoring feature dependencies. In order to handle the problem, the method of multivariate feature selection (MFS), also denoted as dimension reduction, was introduced [13]. MFS uses all the gene expression data simultaneously to select the genes. Until now, many mathematical methods for MFS have been used for gene expression data analysis. For example, Park et al. gave the theoretical analysis on feature extraction capability of class-augmented PCA [14]. Ma et al. used PCA to identify differential gene pathways in [15]. De Haan et al. used PCA to analyze microarray data [16]. Musumarra et al. used PLS to identify genes for new diagnostics [17]. Boulesteix et al. provided a systematic comparison of the PLS methods for the analysis of gene data [18].
However, the classical methods, such as PCA and PLS, still have some drawbacks, e.g. the principal components (PCs) of PCA or the latent components (LCs) of PLS are usually dense, which makes it difficult to interpret PCs or LCs without subjective judgment. To overcome these drawbacks, many mathematical tools have been devised to reduce the complexity of the data. Among them, sparse methods have significant advantage, while giving up little statistical efficiency. For example, Zou et al. proposed sparse PCA using the lasso [19]. In [20], Journée et al. described a general power method for sparse PCA. Lai et al. used sparse local discriminant projections for feature extraction [21]. Moreover, many sparse methods have been widely used for gene expression data analysis. Lass et al. used the SPCA for clustering and feature selection [22]. In [23], Witten et al. proposed a penalized matrix decomposition, which was used to analyze plant gene expression data by Liu et al. [24]. Cao et al. used sparse PLS discriminant analysis for biologically relevant feature selection [25].
Though sparse methods are useful, yet these methods used the first LC of PLS or PC of SPCA to select feature. For example, Boulesteix et al. used the first LC of PLS to select the important genes [18]. Liu et al. used the first PC of SPCA for characteristic genes selection [26]. These methods indeed assumed that the first component of PLS or SPCA plays a dominant role in gene selection. However, we identify that in a number of cases this assumption is not satisfied and conventional PCA-based gene selection methods usually provide poor results. Actually, not only the first LC or PC but also the remaining LCs or PCs include the important information for gene selection [27]. So if only the first LC or PC is used in selecting genes, the poor results may have been obtained with loss of some important information.
In this paper, in order to select characteristic genes, a novel method is proposed that is based on the weighted PCs of SPCA by singular values (WPCS). First, the PCs of SPCA and singular values are calculated. Second, using the singular values as the weights of PCs, the weighted PCs (WPC) are gotten. Then, as the absolute value of the i-th entry of every WPC somewhat denotes the importance of the i-th gene, we take the sum of all the i-th entries of all the WPC as the extent of importance of the i-th gene and use it to select features. The genes corresponding to the first m largest sums are selected as characteristic genes. The experimental results show that our method is efficient and powerful for gene selection. Our work has the following contributions: first, it proposes, for the first time, the method of weighting PCs by singular values for gene selection. Second, from the viewpoint of minimizing reconstitution error, it clearly presents the idea of the proposed method. Third, it conducts a large number of gene selection experiments.

Results and Discussion
In this section, our proposed WPCS method is compared with the following existing methods. (a) SPCA-1 method uses only the first one PC of SPCA (proposed by Journée et al. [20]) to identify the characteristic genes; (b) SPCA-2 method uses the first two PCs of SPCA to identify the characteristic genes; (c) PLS method uses partial least squares regression (PLS) (proposed by Boulesteix et al. [18]) to identify the characteristic genes.
First, these methods are carried on the artificial data. Then, these methods are used to extract the characteristic genes responding to abiotic stresses from real gene expression data.

Simulation on Artificial Data
To investigate the performance of the methods, the average receiver operator characteristic (ROC) curves are shown in Figure 1 with six different SNR. Figure 1(A) and 1(B) show that our WPCS and the competitive methods can identify the patterns even with very low SNR. As shown in Figure 1, with different SNR, WPCS outperforms other methods. For example, with SNR = 1.0, PLS method are dominated by SPCA-2 and SPCA-1; and our WPCS can achieve the best results.
The area under curve (AUC) statistics are listed in Table 1, from which we can conclude that under the same SNR, the ascending order of accuracy given by these methods is: PLS, SPCA-1, SPCA-2 and WPCS.
From the experiments on artificial data, a conclusion can be drawn that WPCS method outperforms other methods for feature selection.

Gene Ontology (GO) Analysis
The Gene Ontology (GO) Term Enrichment tool can be used to help discover what those genes may have in common [28]. GOTermFinder is a web-based tool that finds the significant GO terms shared among a list of genes. The analysis of GOTermFinder provides significant information for the biological interpretation of high-throughput experiments. In this paper, our proposed method will be evaluated by GOTermFinder [29], which is publicly available at http://go.princeton.edu/cgi-bin/ GOTermFinder. Its threshold parameters are set as following: maximum p-value = 0.01 and minimum number of gene products = 2.
Here, only the main results of GO are given. Figure 2 shows the sample frequency of response to stimulus (GO: 0050896) given by the four methods. From Figure 2(A), WPCS method outperforms the others in all the data sets of shoot samples with six different stresses. Figure 2(B) shows that only in droughtstress data set of root samples, our method is dominated by SPCA-1 and SPCA-2 methods. In other data sets, our method is superior to the others. Figure 3(A) shows the sample frequency of response to stress (GO: 0006950) in shoot samples. From Figure 3(A), it can be seen that only in drought-stress data set, the PLS method is slightly superior to our method. In other data sets, our method is superior to the other methods. Figure 3(B) shows that only in drought-stress data set of root samples, our WPCS gives a similar result to that of SPCA-1 and SPCA-2 methods, and exceeds that of PLS method. In other data sets, our WPCS method surpasses the others.
The remarkable results are listed in Tables 2-5. The number of genes responding to stimulus (GO: 0050896) selected by the four methods in shoot and root samples are listed in Table 2 and  Table 3, respectively.
As Table 2 listed, in shoot samples, WPCS method outperforms the others in all the data sets with six different stresses. As Table 3 listed, in root samples, only in drought-stress data set, WPCS method is dominated by SPCA-2. For other stresses data sets, WPCS outperforms our competitive methods. Table 4 and Table 5 give the gene numbers and P-value of response to stress (GO: 0006950) selected by the four methods in shoot and root samples, respectively.
To sum up, for all the data sets except drought-stress data set, our method is superior to other methods. For the drought-stress Table 1. AUC statistics for artificial data. data set in shoot samples, only the PLS method slightly suppresses our WPCS method. To further study the characteristic genes closely related to the stresses, the cold stress in shoot samples and UV-B stress in root samples are analyzed. Table 6 lists the numbers of response to cold (GO: 0009410) in shoot samples selected by these methods. The background sample frequency of response to cold (GO: 0009410) is 0.9% (276/29887). As Table 6 listed, our method can select more genes than others.
In detail, we compare the genes selected by WPCS with the genes selected using others. Different genes selected using WPCS and neglected by other methods are listed in Table 7. As Table 7 listed, the functions of genes selected using WPCS are closely related with cold stress. Table 8 gives the numbers of response to light stimulus (GO: 0009416) in root samples selected using these methods. The background sample frequency of response to light stimulus (GO: 0009416) is 1.8% (547/29887).
As Table 8 listed, WPCS can select more genes than others. Moreover, we compare the genes selected by WPCS with the ones by other methods. The genes selected using WPCS and neglected by others are listed in Table 9. As Table 9 listed, the functions of genes selected using WPCS are closely related with UV-B stress.
From the experiments and analyses on gene expression data, a conclusion can be drawn that WPCS method is very efficient and powerful for gene selection.

Conclusion
In this paper, a novel method of gene selection, WPCS, is proposed, that uses the weighted PCs by SVs as the basis of selection. The idea of WPCS is clearly shown. WPCS works as follows. First, it obtains the PCs of SPCA and SVs. Second, using the SVs as the weights of PCs, it obtains the WPC. Then, it sums the absolute value of the WPC in row, and sorts the sum in descending order. Finally, it selects the genes corresponding to the top part of the sum as the characteristic genes. A large number of experiments on artificial data and gene expression data demonstrate that the proposed WPCS method outperforms the state-ofthe-art gene selection methods. For gene expression data, WPCS can extract more characteristic genes in response to abiotic stresses than the other methods.

Artificial Data
The artificial data are in R p with p~2000 and generated as N(0,1). Then the noise matrix is added toṽ v with different Signal-to-Noise Ratios (SNR). The first four eigenvectors of S 4 are chosen to be v k~ṽ v k =ṽ v k k k, k~1,2,3,4 . To make these four eigenvectors dominate, we let the eigenvalues bec 1~4 00,c 2~3 00, c 3~2 00, c 4~1 00 and c k~1 for k~5, Á Á Á ,2000. Then the simulation scheme in [30] is used to generate the artificial data, which include ten samples in each test.

Gene Expression Data
The raw data include two classes: roots and shoots in each stress, which were downloaded from NASCArrays [http://affy. arabidopsis.info/] [31], reference numbers are: control, NAS-CArrays-137; cold stress, NASCArrays-138; osmotic stress, NASCArrays-139; salt stress, NASCArrays-140; drought stress, NASCArrays-141; UV-B light stress, NASCArrays-144; and heat stress, NASCArrays-146. The sample numbers of each stress type are listed in Table 10. There are 22810 genes in each sample. The data are adjusted for background of optical noise using the GC-RMA software by Wu et al. [32] and normalized using quantile normalization. The results of GC-RMA are gathered in a matrix for further processing.

Selection of the Parameters
In SPCA, we take l 0 -norm penalty and set c~0:1. In PLS, only the first component is used. For the sake of comparison, on gene expression data, 300 genes are roughly selected by all the methods.

Singular Value Decomposition (SVD)
In this subsection, the details of the WPCS method are presented. Let A denote an n|p matrix of real-valued gene expression data, which consists of p genes in n samples. In the case of gene expression data, a ij is the expression level of the j-th gene in the i-th sample. The elements of the j-th column in A form the n-dimensional vector r j , which is referred to as the transcriptional response of the j-th gene. Correspondingly, the elements of the i-th row in A form the p-dimensional vectors i , referred to as the expression profile of the i-th sample. Usuallypwwn, so it is a classical small-sample-size problem. To integrate SPCA and SV, the singular value decomposition (SVD) and the SPCA via cardinality penalty (l 0 -penalty) are introduced as follows.
If the variables contained in the columns of A with rank rƒ min (n,p) are centered, the equation for SVD of A is as follows:  where U is an n|r matrix of left singular vectors with U T U~I, V is a p|r matrix of right singular vectors with V T V~I, and D is a diagonal matrix of singular values. Let u k denote column k-th of U, let v k denote column k-th of V, and note that d k denotes the kth diagonal element of the matrix D. According to Eckart et al. [33], which is the closest rank-l matrix to A. The term ''closest'' means that A (l) minimizes the square error sum between the elements of A and A (l) :

Sparse Principal Component Analysis (SPCA)
The results given with l 0 -and l 1 -norm penalty in SPCA are similar, which is also shown in [20]. Since the l 0 -norm is faster than l 1 -norm, l 0 -norm penalty is taken on SPCA. Let F~UD and Z~V, eq.(1) can be written as follows: where Z T Z~I and F T F~S~A T A. It is the classical PCA formation. Extracting one principal component (PCs) amounts to computing the dominant eigenvector of S (or, equivalently, dominant right singular vector of A). That is, PCA seeks to project the data onto the linear combination of variables that maximizes the sample variance. It is well-known that the solution to this problem is given by the right singular vector of A. In general, PCs is not expected to have many zero coefficients. So, to makes it easy to interpret PCs without subjective judgment, Sparse PCA proposed by Journée et al. in [20] is used to generate the sparse PCs.
Let us consider the optimization problem.
with sparsity-controlling parameter c §0, : k k 0 denotes the l 0norm, that is, the number of non-zero components (cardinality).
According to [20], eq.(4) can be rewritten as follows: where the maximization with respect to z[B p for a fixed v[B n has the closed form solution According to [20], eq.(5) can be cast in the following form:  Here, fort[R, sign(t) denotes the sign of the argument and t z~m axf0,tg.
we get.
for all nonzero vector z, when c . a i Ã k k 2 2 .So, this derivation assumes that cv a i Ã k k 2 2 , then.
Otherwisec § a i k k 2 2 , so z Ã i (c)~0, i~1,:::,p: The Idea of the Proposed Method The residual sum of squares (RESS) can be used for evaluating the quality of the reconstitution of A with M PCs of SPCA. According to [34], it can be computed as follows: whereÂ A ½M is a reconstitution of A, Y denotes the sum of all the squared elements, and l l is the singular value of the l-th component. The smaller the value of RESS is, the better the SPCA model is. From eq.(12), we can see that a larger M may give a better estimation ofÂ A ½M . If M takes value r (rank of A),  Table 7. Different genes of response to cold (GO: 0009410) in shoot samples.

At1g21910
Participates in plant developmental processes as well as biotic and/or abiotic stress signaling.

At1g22770
Regulates several developmental processes, such as circadian clock, carbohydrate metabolism, and cold stress response.

At1g29395
Expression is induced by short-term cold-treatment, water deprivation, and abscisic acid (ABA) treatment.

At2g19450
Role in senescence and seed development induced by cold-stress.

At2g25930
Temperature stress reduced the pyk20 transcript level.

At2g28900
Predominantly expressed in leaves and is also inducible by cold treatment.

At2g33380
Plays a role as a peroxygenase involved in oxylipin metabolism during biotic and abiotic stress.

At2g38470
Involved in response to various abiotic stresses

At2g47180
Increases tolerance to chilling stress At3g05880 Induced by low temperatures, dehydration and salt stress.

At3g48360
Mediates multiple responses to nutrients, stresses, and hormones.

At3g53990
Low temperature and salt responsive protein family.

At4g30650
Low temperature and salt expression protein homologous.

At4g30660
Putative low temperature and salt responsive protein.

At4g37610
Under cold stress indicates increased expression.

At5g52300
Induced by low temperature, exogenous abscisic acid (ABA) and drought.

At5g57560
Controlling tolerance to cold stress doi:10.1371/journal.pone.0038873.t007 where the matrix A can be perfectly reconstituted. Let M~1, the RESS 1 can be expressed as follows: Substituting eq.(13) into eq. (14), the RESS 1 can be obtained as follows: As eq. (15) shown, if only one PC is used to reconstitute the matrix A, RESS 1 may be larger P r l~2 l l than RESS r . So, if only the first PC is used for characteristic gene selection, some important information may be lost, especially when the second one or two SVs are approximately equal to the first one. In order to obtain a better reconstitution, all the PCs of SPCA need to be utilized. In SPCA, F is the matrix of factor scores and Z is a loading matrix of the principal components (PCs), which transforms the original data matrix into factor scores. The data matrix A, factor scores matrix F and PCs Z are shown in Figure 4.
As Figure 4 shown, the PCs Z give the coefficients of the linear combinations used to compute the factors scores F. So the bigger the absolute value of the elements in PCs Z is, the more contribution it gives for the factor scores matrix, the more important the corresponding gene in A is. So the characteristic genes can be selected according to the PCs Z. Let.
denote the i-th PC, and then the PCs can be given as the follows: Substituting F~UD into eq.(3), it can be reformed as follows: where D is the diagonal matrix of singular values. Let eq.(18) can be reformed as follows: Substituting eq.(17) into eq.(19),  Table 9. Different genes of response to light stimulus (GO:0009416) in root samples.

Gene No. Function of Gene
At2g29500 HSP20-like chaperones superfamily protein.

At3g54890
Encodes a component of the light harvesting complex associated with photosystem I.

At3g55120
Catalyzes the conversion of chalcones into flavanones.

At5g02810
Acts as transcriptional repressor of CCA1 and LHY.

At5g12030
Encodes a cytosolic small heat shock protein with chaperone activity that is induced by heat and high light intensity stress.

At5g24470
Encodes a pseudo-response regulator whose mutation affects various circadian-associated biological events such as red light sensitivity of seedlings during early photomorphogenesis.
The matrixẐ Z is referred to as weighted PCs (WPC), which can be obtained via weighting PCs by the diagonal matrix of SVs (WPCS).
Substituting eq.(16) into eq.(21), the WPC can be reformed as follows:Ẑ As the absolute value of the i-th row of WPCẐ Z somewhat denotes the importance of the i-th gene, the absolute value sum of all the entries in the i-th row as the evaluating vector EV, which can be expressed as follows: : ð23Þ In particular, if the dimensionality of the gene data is p, the EV has p entries. After sorting the evaluating vector EV, the genes corresponding to the first m largest entries can be selected as characteristic genes.
In summary, the main steps of WPCS method are shown as follows.
(2) To obtain the PCs Z, execute SPCA on the A.    Figure 4. The graphical depiction of SPCA of a matrix A with factor scores F and PCs Z. In this figure, with factor scores F and PCs Z.r r j is the row vector of PCs Z the j-th gene, which transforms the original data vector r j into factor scoresr r j . Correspondingly,s s i is the column vector of PCs Z, which transforms the original data vector s i into factor scoresŝ s i . doi:10.1371/journal.pone.0038873.g004 (7) Sort the EV in descending order. (7) Select the genes corresponding to the first m largest entries as characteristic genes.
The workflow diagram of our method is shown in Figure 5.