{'fraction_non_alphanumeric': 0.03129406235312549, 'fraction_numerical': 0.01778160739464202, 'mean_word_length': 4.728741305811083, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mining Gene Expression Profiles: An Integrated Implementation of Kernel Principal Component Analysis and Singular Value Decomposition

The detection of genes that show similar profiles under different experimental conditions is often an initial step in inferring the biological significance of such genes. Visualization tools are used to identify genes with similar profiles in microarray studies. Given the large number of genes recorded in microarray experiments, gene expression data are generally displayed on a low dimensional plot, based on linear methods. However, microarray data show nonlinearity, due to high-order terms of interaction between genes, so alternative approaches, such as kernel methods, may be more appropriate. We introduce a technique that combines kernel principal component analysis (KPCA) and Biplot to visualize gene expression profiles. Our approach relies on the singular value decomposition of the input matrix and incorporates an additional step that involves KPCA. The main properties of our method are the extraction of nonlinear features and the preservation of the input variables (genes) in the output display. We apply this algorithm to colon tumor, leukemia and lymphoma datasets. Our approach reveals the underlying structure of the gene expression profiles and provides a more intuitive understanding of the gene and sample association.


Introduction
Microarray technology has been advanced to the point at which the simultaneous monitoring of gene expression on a genome scale is now possible. Microarray experiments often aim to identify individual genes that are differentially expressed under distinct conditions, such as between two or more phenotypes, cell lines, under different treatment types or diseased and applied to detect genes that have similar expression profiles from microarray experiments. Techniques like clustering (2) and self organization map (3) have been applied to the analysis of gene expression data. We can also find several applications on microarray analysis based on distinct machine learning methods such as Gaussian processes (4,5), Boosting (6) and Random Forest (7). It is useful to find gene/sample clusters with similar gene expression patterns for interpreting the microarray data.
However, due to the large number of genes involved, it might be more effective to display these data on a low dimensional plot. Recently, several authors have explored dimension reduction techniques: Alter et al (8) analyzed microarray data using singular value decomposition (SVD), Fellenberg et al (9) applied correspondence analysis to visualize genes and tissues, Pittelkow and Wilson (10) and Park et al (11) used several variants of biplot methods as a visualization tool for the analysis of microarray data. Visualizing gene expression may facilitate the identification of genes with similar expression patterns.
In this paper we describe a method to visualize gene expression profiles. Our procedure relies on SVD; however, unlike other methods, it incorporates an additional step that involves kernel principal component analysis (KPCA) (12). Kernel representation offers an alternative to nonlinear functions by projecting the data into a high-dimensional feature space, which increases the computational power of linear learning machines (13,14).
Kernel methods enable us to construct different nonlinear versions of any algorithm that can be expressed solely in terms of dot products, known as the kernel trick. Thus, kernel algorithms avoid the explicit usage of the input variables in the statistical learning task. Kernel machines can be used to implement several learning algorithms but they usually act as a black-box with respect to the input variables. This could be a drawback in biplot displays in which we pursue the simultaneous representation of samples and input variables. As we describe in greater detail below, the integrated implementation of SVD-Biplot joint with KPCA enables us to extract the nonlinear features without discarding the simultaneous display of input variables (genes) and samples (microarrays).

Method
We describe an effective procedure that allows the simultaneous display of input variables (genes) and samples (microarrays). The standard SVD-Biplot allows us to represent samples and variables jointly, but our procedure incorporates an additional extraction of nonlinear features that could improve our understanding of the relationships between samples (microarrays) and variables (genes), and enhance the detection of gene expression profiles. Since microarray data are generally nonlinear, finding methods that can handle such data is of great importance if as much information as possible is to be gleaned.
Let X be the preprocessed gene expression data matrix, with n samples (microarrays) in the rows and p variables (genes) in the columns. Preprocessing of the gene expression measurements needs to be considered with caution because preprocessing, such as scaling, normalization and transformation, can have a strong effect on the output visualization (15).
Our procedure represents the gene expression data in a feature space using a kernel method. Initially, we decompose the input data matrix X=GH T into singular values, for an n×r matrix G and a p×r matrix H, representing sample (microarray) and variable (gene) effects, respectively.
Next, we map the variables (genes), the rows of H, into a feature space, and we extract the main nonlinear features of variables (genes) by performing KPCA. Finally, to obtain a simultaneous plot, we project the samples (microarrays) into the subspace spanned by the leading eigenvectors from KPCA.
Our procedure is composed of the following steps: 1. SVD of preprocessed gene expression input matrix X=GH T .
2. Take the rows of H as a set of observations and compute the corresponding kernel matrix K.
3. Extract nonlinear gene expression features by computing KPCA on the kernel matrix K. 4. Project the rows of G onto the subspace expanded by the leading eigenvectors of K.
Step 1 allows us to represent the microarrays (samples) and the gene expressions (input variables) as points on the same r-dimensional space. We obtain the matrices: where vectors g 1 ,…,g n represent microarrays and vectors h 1 ,…,h p represent gene expressions.
Step 2 serves to build the kernel matrix K. To compute K on the sample observations, the choice of the kernel function K(x, y) is crucial because the kernel methods use it to measure the similarity of sample observations. Gene expression data are a set of numerical vectors, so two natural candidates for kernel function K(x, y) are the polynomial kernel: K(x, y)=(xy+c) t where t, c>0 are free parameters, and the radial basis kernel is: K(x, y)=exp(−c||x−y|| 2 ) where c>0 is a free parameter. Once we select a suitable kernel function, we compute the p×p kernel matrix K=(K(h i ,h j )).
Step 3 allows us to extract nonlinear features of gene expressions by solving the dual eigenvalue problem, such as Formula 8 showed below. Notice that in this step we have considered the vectors h 1 ,…,h p , representing gene expressions (Formula 2), as input observations. At the end of this step we find: α 1 ,…,α r , the set of normalized eigenvectors, and λ 1 ≥ λ 2 ≥ … ≥ λ r , the corresponding set of nonzero eigenvalues of K. Formula 9 (see below) allows us to obtain the coordinates of the gene expression h l (l=1,…,p) by projecting onto the leading kernel principal components where j=1,...,r. As is customary, we use a small number of eigenvectors in order to obtain a graphical output in low dimension. The stability of KPCA, that is, the extent to which the projection captures new data sampled according to the same distribution as the training data, increases when the eigenvalues become small at an initial stage in the spectrum of K. Provided we project into a space whose dimension exceeds the index of this stage, we expect to capture most of the variance of the unseen data.
Step 4 is to project the sample (microarrays) points g l , l=1,...,n (see Formula 1), onto the subspace spanned by the eigenvectors of K. The new coordinates are given using Formula 9, which is where j=1,...,r.
We can analyze and better understand the association between genes and samples by superimposing the coordinate Formulas 3 and 4 in the same plot. Furthermore, this integrated implementation of SVD-Biplot and KPCA provides a supplementary tool to detect gene expression profiles in microarray experiments. For brevity we refer to our method as the KPCA-Biplot.

Kernel selection and automatic tuning
We could consider the choice of the kernel and the tuning procedure for the kernel parameters as a model selection problem.
As is well known in the context of supervised learning (e.g., classification tasks), we can handle this selection by including a cross-validation step in the procedure. We therefore guide the kernel selection and the tuning of the parameters automatically using a cross-validation score function based on the minimization of the misclassification rate.
In the context of unsupervised learning (e.g., dimension reduction tasks), an approach to include a cross-validation step for data-driven kernel selection, when the samples are grouped (e.g., different treatment types, cell lines and phenotypes), is to adopt a strategy similar to that followed in the supervised case (16). In summary, assuming a fixed number of dimensions, usually 2 or 3, we perform the dimension reduction for each kernel. Then, using the sample coordinates in the reduced dimension space, a classifier is trained, and assessed with respect to a test dataset. Finally, we select the kernel that minimizes the misclassification rate.

Computational complexity
The speed and the time of calculations have always been important problems in the implementation of algorithm in multivariate data analysis. Microarray datasets typically consist of thousands of variables and less than 100 samples. For such huge datasets, the computation time needed for matrix decomposition using classical SVD algorithm and/or the eigenvalue decomposition for the kernel matrix may be excessive.
Alter et al (8) propose an efficient algorithm for SVD for genome-wide expression data. In summary, the algorithm is as follows. For the dimension of data matrix X (n×p), where n denotes the number of samples and p denotes the number of genes, we assume that the data are centered. For matrices U n×r , D r×r and V r×p , the p×p sample correlation matrix admits the following eigenvalue decomposition A=X T X= VDU T UDV T =VD 2 V T , and the n×n gene correlation matrix admits the following eigenvalue decomposition B=XX T =UDV T VDU T =UD 2 U T . We therefore calculate the SVD of X (step 1 of the proposed method) by diagonalizing B and computing V as V=X T UD −1 .
Furthermore, since the kernel matrix K is p×p (where p is assumed to be from thousands to tens of thousands in real application), reducing the cost of computing might also be appropriate in the kernel matrix eigenvalue decomposition (step 3 of the proposed method). When we need to analyze a large number of genes, we may want to work with an algorithm for computing only the largest eigenvalues, as for instance the power method with deflation (17).

SVD
Let X be the preprocessed gene expression data matrix, with n samples (microarrays) in the rows and p variables (genes) in the columns. Underlying the biplot techniques (18,19) is the SVD: X =UDV T (5) where U and V are matrices of size n×r and p×r, respectively, with orthonormal columns so that U T U=VV T =I r , D is an r×r diagonal matrix with elements λ 1 ≥…≥ λ r in the diagonals, and r is the rank of X, so usually r = min(n, p). Let us define and let G=UD α and H=VD 1−α , where 0 ≤ α ≤ 1, thus X can be factorized as X=UDV T =UD α D 1−α V T =GH T for an n×r matrix G and p×r matrix H. Thus X can be decomposed into two sets of matrices G and H, representing row and column effects, respectively.

KPCA
Given a set of observations x i ∈ℜ p , i=1,…,n, let us consider a feature space F related to the input space by a map φ: ℜ p →F, which may be nonlinear. We assume that we are dealing with centred data In F the covariance matrix takes the We seek eigenvalues λ ≥0 and nonzero eigenvec- Combining Formulas 6 and 7 we obtain the dual representation of the eigenvalue problem for nonzero eigenvalues: K α=nλα (8) where K=(K(x i ,x j )), i, j=1,…,n, is the kernel matrix and K is a kernel function such that the dot product in which is the j-th nonlinear principal component cor-responding to φ.

Validation
In this section we illustrate the application of KPCA-Biplot with data from the colon tumor (20), leukemia (21) and lymphoma (22) datasets. In these examples, the aim of the KPCA-Biplot is to detect genes (variables) that have a similar pattern of up/down-regulation for each sample. By simultaneously displaying both the samples and the genes on the same plot, it is possible both to visually detect genes that have similar profiles and to interpret this pattern by reference to the sample groups. From the position of the genes relative to the samples, it can be deduced that genes, which lie relatively close to any given group, will have higher values (up-regulated) in that group than in the other groups. Genes lying on the opposite side of the origin from a given group will tend to have lower values (down-regulated) in that group. Then gene profiles are useful to reveal differential expression between sample groups.
As an example, we describe the profiles of some illustrative genes that are located away from the central gene cloud in each genomic dataset. In particular, with the aim of detecting different profiles, we explore several directions from the origin of the graphical output and describe the profiles of a set of genes that lie in those directions. By combining the KPCA-Biplot and the plot of the profiles, we can represent all the different kinds of profile on one plot.

Colon tumor dataset
This dataset is composed of 2,000 genes in 40 colon tumor samples and 22 normal colon tissue samples. Gene expression levels were analyzed with Affymetrix oligonucleotide array. The dataset is available on the web at http://www.molbio.princeton.edu/ colondata. We complete the preprocessing of the gene expression data with a microarray standardization and gene centring (10).
We performed the KPCA-Biplot procedure on this dataset as detailed in the previous section. Initially, we factorize the input data matrix using SVD with α=0.5. Next, we compute the kernel matrix in Step 2 using the radial basis kernel with c=0.1. The resulting biplot is given in Figure 1. It shows the projection onto the two leading kernel principal components of genes and microarrays, which were obtained by using Formulas 3 and 4, respectively. Genes are shown as a tilde and microarrays are shown by different marks according to the class of colon tissues to which they belong.
The first axis mainly distinguishes the sample and gene clouds. It also separates those highly expressed genes from the origin. The second axis reveals normal and tumor samples. Nevertheless, normal and tumor samples are not completely separated in this 2D plot. In addition, some genes are labeled with their gene numbers and these have been chosen to illustrate different profiles. Genes that lie towards the left and bottom will have higher expression in tumor samples. In contrast, genes at the top, close to the normal samples, will tend to have higher expression in these samples. Figure S1 shows the gene expression profiles of some genes labeled in Figure 1. The gene expression profiles are computed from the expression values on all arrays for a given gene in different groups of arrays, and they are summarized in a line plot connecting the median of each group. Each row of Figure S1 contains the gene expression profiles of genes lying in the same directions of the KPCA-Biplot. In order to explore different profiles, we have selected three directions as an example. Genes T95018, T61602, T58861, T48804 and T57633 (the top row) have higher expression in tumor samples. In contrast, genes R78934, T92451, M33680 and T42-control (the bottom row), are most strongly expressed in normal cells. Genes T51560, T49703, U05012, T51496 and R44112 (the middle row) are equally expressed in both samples.
These results are consistent with those reported by Alon et al (20), in which the analysis of the same dataset by using two-way data clustering reveals a group of genes whose expression is correlated across tissue types. In particular, it reports a list of 48 ESTs homologous to ribosomal proteins within the set of 2,000 genes used for the clustering. The intensity of the ribosomal protein genes is relatively low in normal colon tissues and high in the colon tumor tissues. We remark that genes T95018, T61602, T58861, T48804 and T57633, which are highlighted by our method, belong to the list of 48 ESTs homologous to ribosomal proteins reported by Alon et al (20). Figure 2 shows the standard biplot representation obtained in analyzing the same dataset. We observed that the sample groups are separated mainly along the second axis but the separation is not completed. In this example, with respect to the sample representation, both the standard biplot and KPCA-Biplot detect the group structure in reduced dimension but the two groups are not fully separated by either method. With respect to the gene representation, both the standard biplot and our method reveal the normal and tumor profiles but we observed that with the KPCA-Biplot the genes are more clearly separated, which can be helpful for mining the microarray data.

Leukemia dataset
The leukemia dataset is composed of 3,051 gene expressions in three classes of leukemia: 19 cases of B-cell acute lymphoblastic leukemia (ALL), 8 cases of T-cell ALL and 11 cases of acute myeloid leukemia (AML). Gene expression levels were measured using Affymetrix high-density oligonucleotide arrays. The data can be downloaded from http://www.genome.wi. mit.edu.
The data were preprocessed according to the protocol described by Dudoit et al (23). In addition, we complete the preprocessing of the gene expression data with a microarray standardization and gene centring.
We performed the KPCA-Biplot procedure on this dataset as detailed in the previous section. Initially, we factorize the input data matrix using SVD with α=0.5. Next, we compute the kernel matrix in Step 2 using the radial basis kernel with c=0.001, which is set heuristically. In this step we also use the polynomial kernel but its behavior, in this case, is similar to that of the radial basis kernel. For brevity we only show results obtained with the radial basis kernel. The resulting biplot is given in Figure 3. It shows the projection onto the two leading kernel principal components of genes and microarrays, which were obtained by using Formulas 3 and 4, respectively. Genes are shown as a tilde and microarrays are shown by different marks according to the class of leukemia to which they belong. The first axis separates ALL (B and T-cells) and AML samples. In particular, the 2D plot detects the three classes of leukemia.
Additionally, in Figure 3 some genes are labeled with their probe IDs and these genes have been chosen to illustrate different profiles. It is expected that genes that lie towards the right will have higher expression values in the AML samples than in the ALL samples. In contrast, genes on the left will have higher expression values in the ALL samples than in the AML samples. Genes placed near the bottom will tend to have higher expression values in T-cells than in B-cells. Figure S2 shows the gene expression profiles of some genes labeled in Figure 3. The gene expression profiles are computed from the expression values on all arrays for a given gene in different groups of arrays and they are summarized in a line plot connecting the median of each group. Each row of Figure S2 contains the gene expression profiles of genes lying in the same directions of the KPCA-Biplot. In order to explore different profiles, we have selected three directions as an example. It is observed that genes lying close to any group of leukemia have higher values in that group than in other groups, and genes towards the opposite group have lower values in that group than in other groups. Genes M57710_at, M28130_rna1_s_at and Y00787_s_at are highly expressed on AML samples and they have a similar profile as shown in Figure S2. In the opposite direction we can find genes L47738_at and Z69881_at, which are down-regulated in AML samples. As expected, these genes have the opposite profile of expression.
The expression of the gene HG3576-HT3779_f_at is highest in AML and ALL-B cell samples but lowest in ALL T-cell samples. In contrast, genes X76223_s_at and X00437_s_at are most strongly expressed in ALL T-cell samples. The shape of this set of genes is opposite to that of HG3576-HT3779_f_at.
Genes X82240_rna1_at and L33930_s_at are highly expressed on ALL-B cell samples. In the opposite direction we can find the gene X05908_at, which is highly expressed in AML and ALL T-cells but expressed most weakly in ALL-B cells.
These results indicate that our method provides a useful tool to find genes with different expression. Our method is complementary to other current methods to detect genes with different expression. For example, genes Z69881_at and M28130_rna1_s_at are 2 of the 50 genes selected in Golub's study (21) to differentiate AML from ALL cells. Genes HG3576-HT3779_f_at, X76223_s_at and X82240_ rna1_at are also selected in Pittelkow and Wilson's study (10).

Lymphoma dataset
The lymphoma dataset comes from a study of gene expression of three prevalent lymphoid malignancies: B-cell chronic lymphocytic leukemia (B-CLL), follicular lymphoma (FL) and diffuse large B-cell lym-phoma (DLCL). Among 96 samples we took 62 samples containging 4,026 genes in three classes: 11 cases of B-CLL, 9 cases of FL and 42 cases of DLCL. Gene expression levels were measured using two-channel cDNA microarrays. The data can be obtained from http://genome-www.stanford.edu/lymphoma.
After preprocessing, all gene expression profiles were base-10 log-transformed and standardized to zero mean and unit variance in order to prevent single arrays from dominating the analysis. Finally, we complete the preprocessing of the gene expression data with gene centring.
To perform the KPCA-Biplot procedure, we factorize the input data matrix using SVD with α=0.5.
Next, we compute the kernel matrix using the radial basis kernel with c=0.01. In this step we use alternatively the polynomial kernel, but it yields the poorest representation of the dataset. The resulting biplot is given in Figure 4. It shows the projection onto the two leading kernel principal components of genes and microarrays, which is obtained by using the steps of kernel selection and automatic tuning. Genes are shown as a tilde and microarrays are shown by different marks according to the type of lymphoma to which they belong. The three classes are clearly separated. The first axis separates DLCL cells and the others, and the second axis separates FL and B-CLL cells.
As we can see in Figure 4, some genes labeled with their clone IDs have been selected to visualize different profiles in this example. It is expected that genes lying towards the right will have higher values for DLCL samples than for FL and B-CLL samples. Genes placed near the top will tend to have higher expression values in FL samples than in B-CLL samples. Figure S3 shows the profiles of some genes that have been labeled in Figure 4. Each row of Figure S3  contains the gene expression profiles of genes lying in the same directions of the KPCA-Biplot. In order to explore different profiles, we have selected four directions as an example. It is observed that genes lying close to any group of lymphoma are up-regulated in that group, and genes lying opposite a group are down-regulated.
Genes 1319066, 1354294 and 428103 show an expression profile with high values in FL samples and low values in DLCL and B-CLL sample cells. In contrast, genes 234357 and 1306024 have opposite profile of expression.
Genes 347751 and 724070 are highly expressed in B-CLL samples and weakly expressed in FL and DLCL samples. In the opposite direction we can find gene 1334260, which is down-regulated in B-CLL samples and up-regulated in the others.
Genes 139009 and 8 are highly expressed in DLCL samples. In contrast, gene 1338881 has the reverse shape. Finally, we observed gene 1352822, which has a linear profile, and genes 685456 and 712829, lying in the opposite direction, also with linear profile but with the opposing slopes.

Conclusion
In this paper we propose an exploratory method, called KPCA-Biplot, that combines KPCA and SVD-Biplot for elucidating relationships between samples (microarrays) and variables (genes). The main properties of KPCA-Biplot are the extraction of nonlinear features together with the preservation of the input variables (genes) in the output display. The method described here is easy to implement and facilitates the identification of genes that have similar or reversed profiles. Our results indicate that KPCA-Biplot is complementary to other tools currently used for finding gene expression profiles, with the advantage that it can capture the usual nonlinear nature of microarray data.