Data normalization and classification analysis

Before testing the performance, all the data are properly normalized by being first log-transformed, and then scaled to mean 0 and standard deviation 1 (i.e., Z-normalization) in each sample and then each gene. All the normalized datasets are available for free at our Internet server [18]. Also the ME program that runs on Linux OS is available upon request. Many genes have a large number of missing values because heterogeneous data are combined; thus, we adopt a simple imputation method that all the missing values are replaced with the mean value, i.e., 0. Input genes that show high correlation to class labels, or feature genes, are selected by the standard two sample t-statistics [19] in each iteration of the leave-one-out cross-validation (LOOCV) test. The distance constraint matrices (D _ij ) are also generated from the same feature genes. If a sample contains missing values, we again adopt a simple imputation; we replace the one-dimensional Euclidean distance (x _ih - x _jh )² with 2 if x _ih or x _jh is missing. The six kernels are tested with SVMs to analyze their classification performance with various numbers of feature genes and various parameters described in Methods. The maximum accuracy among the tested parameters for each number of feature genes is recorded as the accuracy for each kernel.

Heterogeneous kidney carcinoma data

Classification analysis is performed between normal and carcinoma data. The results of the LOOCV test of 100 samples against various numbers (8–296; increasing 8 genes at each step due to computational limitations) of feature genes are plotted in Figure 3. The figure shows typical prediction curves, namely, accuracy increases with increasing number of feature genes, plateaus at some region, and decreases. Clearly, the ME kernel performs much better in all cases than the other five kernels for small numbers of feature genes (8–192). As regards accuracy, the ME kernel records maximum accuracies as high as 95.0 (89.5/98.4 sensitivity/specificity)% for 152 of feature genes. Statistically, the accuracies of the ME kernel are superior to those of the other five kernels in 64.9% of the tested points (8–296) of feature genes. This percentage increases to 95.8% when accuracies are limited only to the increase and plateau regions (8–192) of the ME kernel.

kNND denoising for AML and ALL data

Acute myeloid leukemia (AML) and acute lymphoblastic leukemia (ALL) data for cancer subtype classification have been reported by Golub et al. [1]. There are 72 samples (47 AML and 25 ALL), all of which are quite homogeneous and of good quality, and are thus suitable for artificial noise experiments. To assess the denoising ability of our ME kernel, we first replace the ν_add × 100% of original data in a gene expression profile with artificial white noise, i.e., the noise is added according to a normal distribution model with a mean of 0 and a standard deviation of twice that of each gene value distribution in the original dataset. Then, we extract 50 feature genes from the training dataset for each iteration of the LOOCV test by standard t-statistics.

As the control experiments using linear and RBF kernels, the standard singular value decomposition (SVD) denoising method is applied to reduce noise immediately after the noise is introduced. In the SVD denoising, three levels of noise removals by dfferent cumulative proportions, 85, 90, and 95%, of eigenvalues are explored. For the ME kernel, the k NND denoising method with the following noise level settings is applied. First, raw noise that is assumed to internally exist in the original data is arbitrarily set at ν_raw = 0.05. Then, we define the total noise level as the sum of the raw noise and the above artificially added noise, ν_add. For example, if 10% noise is added, the total noise level is ν_raw + ν_add = 0.05 + 0.1 = 0.15, and (1 - 0.15)² × 100 = 72.3% of the nearest distance genes out of the feature gene set are considered in calculating the k NNDs between samples (see Methods).

We repeat the above random noise-adding test ten times and average the highest accuracies among various parameter combinations. The results are shown in Figure 4. The artificial noise added is within the range of 0–50%. Since the raw data are quite homogeneous, all kernels except linear and polynomial show the same prediction accuracy of 98.6% when no noise is added. This value decreases gradually with increasing noise levels (10–50%) for the vectorial kernels; for example, the accuracies of the RBF kernel decrease in the order of 96.2, 95.9, 91.0, 82.5, and 79.5%. SVD denoising boosts up these accuracies to 98.0, 96.6, 93.2, 91.0, and 86.5%, respectively. The linear and polynomial kernels also show similar accuracies to the RBF kernel when SVD denoising is used.

Interestingly and surprisingly, the three k NND-distance-based methods show high accuracies; for example, the k NND-ME kernel has an accuracy of 97.8% even at 20% noise level and maintains high accuracies of 97.2 and 92.0% at 30–40% noise levels. The EKM and Saigo kernels using k NND-distance also show similar accuracies to the k NND-ME kernel. To verify our results, we extensively analyzed the same data with various parameters including many cumulative proportions in the SVD but obtained similar tendencies, confirming the superior denoising ability of the k NND-based method [21].

Heterogeneous oral cavity carcinoma metastasis data

We further analyze the total performance of the k NND-ME method with a more practical problem-heterogeneous oral cavity carcinoma metastasis data. The data consist of two GEO datasets (GSE2280 and GSE3524) from different authors [22, 23]. One dataset (GSE2280) is derived from primary squamous cell carcinoma dataset of the oral cavity [22], containing 14 metastasis (samples from lymph node tissues are excluded) and eight non-metastasis samples. The other oral squamous cell carcinoma dataset (GSE3524) is comprised of nine metastasis and nine non-metastasis samples (two of stage-unknown samples are excluded) [23]. Both are from the same platform, Affymetrix HG-U133A, where 22,283 genes are analyzed. The size of each dataset is too small and not suitable for SVM classification if analyzed separately. However, combining the two datasets, we obtain as many as 23 metastasis and 17 non-metastasis samples, making it possible to carry out the classification analysis.

The results of the LOOCV test of the 40 samples against various numbers (1–100; increasing one gene at each step) of feature genes with four different kernels, namely, linear, polynomial, RBF, and ME with k NND denoising (k NND-ME), are shown in Figure 5. In the k NND-ME kernel, five different noise levels, ν = 0 (no noise), 0.05, 0.1, 0.15, and 0.2 are evaluated. For comparison, we also classify the two datasets separately and average the accuracies (Figure 5a). The results clearly show that the k NND-ME kernel surpasses the other kernels in both averaged and mixed datasets. Statistically, the accuracies of the k NND-ME kernel are superior to those of the other three kernels in the averaged and the mixed datasets in 98% and 48%, respectively, of all the tested points. The difference in accuracy is much greater in the averaged dataset than in the mixed dataset. The mean differences between the k NND-ME kernel and either of the linear, polynomial or RBF kernel with highest accuracy at each point in the averaged and the mixed datasets are 7.7% and 1.4%, respectively. The accuracies increase and plateau at around 3–30 feature genes in the mixed dataset, while no clear increase or plateau is found for the averaged dataset. The overall maximum accuracy of 87.5 (91.3/82.4 or 82.6/94.1 sensitivity/specificity)% is observed for the k NND-ME kernel at two points, 7 and 15 feature genes, in the mixed dataset. Those accuracies are obtained with ν = 0 and 0.05 denoising parameters. The result also indicates that the k NND-ME kernel shows more stable and higher accuracies than the other kernels for large numbers of feature genes.

Incidentally, the top 15 feature genes that show the highest average ranks by t-statistics in the LOOCV test and that are considered to be associated with oral carcinoma metastasis are: HFE (AF150664), FLJ12529 (NM_024811), CXorf56 (NM_022101), HEATR1 (NM_018072), MGAM (NM_004668), APOL3 (NM_014349), PYY2 (NM_021093), RBP3 (J03912), UBE2V1 (NM_003349, NM_021988, NM_022442, NM_199144, NM_199203), KCNJ15 (U73191), GLS (AB020645), ARHGEF3 (NM_019555), MDM1 (NM_020128), ZC3H13 (AL136745), and C9orf16 (NM_024112).

Although a sufficient number of genes (a total of 22,283 genes) are used for SVD denoising, the denoised dataset does not significantly improve the raw accuracies in small numbers (≤ 30) of feature genes, where the overall maximum range accuracy (84.0%) exists. SVD denoising affects only large numbers (≥ 31) of feature genes. This is probably related to the property of SVD denoising that affects the ratio of information to noise content. Further analysis is needed to understand the full property of the SVD denoising method. In summary, the maximum accuracy (87.5%) of the k NND-ME kernel in raw data is not improved by SVD denoising (data not shown).

Properties of kernel methods

Kernels are numerical expressions of similarity metrics between two samples. The basic form of kernels for two sample vectors with M dimensions x _i = (x_{i 1},..., x _iM ) and x _j = (x_{j 1},..., x _jM ) is represented by an inner product function such as K(x _i , x _j ) = ϕ(x _i )·ϕ(x _j ). ϕ(·) means an arbitrary mapping of vectors to another space with generally different dimensions called 'reproducing kernel Hilbert space (RKHS)' [7], which has many properties common to those of the Euclidean space. For technical convenience, rather than defining the mapping functions ϕ(·) for x _i and x _j , the inner product forms, i.e., K(x _i , x _j ) = ϕ(x _i )·ϕ(x _j ) of the mappings of x _i and x _j are preferably used in practical calculation [7]. The function K(x _i , x _j ) is called a kernel.

Three standard kernels popularly used in microarray studies are the linear kernel: K(x _i , x _j ) = x _i ·x _j , the polynomial kernel: K(x _i , x _j ) = (x _i ·x _j + 1) ^D , and the RBF kernel: $K(x_i,x_j)=\exp (-{\displaystyle {\sum }_{k=1}^M{(x_{ik}-x_{jk})}^2/{\sigma }^2)}$. All these kernels belong to the vectorial data kernel family that takes vectorial data for N samples as input, and we can fill all the (i, j) elements in the N × N kernel matrix with the specified kernel function.

Any kernel matrix generated from such kernel functions possesses a necessary property for SVM learning called positive semidefiniteness (see Appendix for details). If a kernel matrix is positive semidefinite, the mapped vectors ϕ(x) exist in the RKHS where the triangle inequalities among the mapped vectors are conserved. Our aim is to develop kernels that are robust to heterogeneous and noisy gene expression data. To this end, we first devise a new distance metric called k NND (detailed later) that can fulfill our requirements. However, unfortunately, the triangle inequalities are not conserved in the metric. To construct valid kernels from such distances, we introduce the following ME kernel algorithm.

ME kernel with kNND denoising

kNND denoising method for ME kernel

The main issue addressed herein is how to handle missing or noisy values that exist in a large portion of a gene expression profile consisting of heterogeneous data. To effectively eliminate such spurious values without removing the entire gene, we devised the following simple method. Assume that we have a gene expression table (i.e., M genes × N samples matrix) where a sample contains ν (×100)% of noisy genes on average. In such a case, only 1 - ν of genes in that sample contain no noise. Therefore, for any pair of samples, the ratio of common genes not containing noise is expected to be (1 - ν)². Based on this observation, we compute the distance between two samples x _i and x _j as follows: First, we compute the one-dimensional Euclidean distances d _h = (x _ih - x _jh )² for h = 1, ... , M genes. Then, we select k = (1 - ν)² × M of one-dimensional Euclidean distances d _h from the nearest (smallest) ones. Finally, we take the sum of the selected d _h s as the distance between x _i and x _j . We refer to this method as k-nearest neighbor gene distance denoising (k NND denoising) method hereafter. For instance, if a sample with M = 100 feature genes contains ν = 15% of noisy values, k = (1 - 0.15)² × 100 ≃ 72 of the nearest distance genes out of the 100 feature genes are only considered in calculating k NNDs between samples.

Multiplying the above k NNDs by constant G for N samples, an N × N distance constraint matrix (D _ij ) is generated. Note that since the samples use different gene sets in k NND metric, positive semidefiniteness will not hold when directly imported to the RBF kernel function. Instead, however, when our ME kernel algorithm is applied and those k NNDs are used as constraint, the resulting kernel matrix holds positive semidefiniteness as well as reflects similarities between samples. Subsequently, we will train SVMs for the optimized 'ME' kernel with sample labels.

Other distance-based kernels

For comparison, we tested two approaches to obtain a kernel matrix from distance matrix D _ij . Both approaches are originally devised for conversion of a non-positive-semidefinite similarity matrix S into a kernel matrix. The first approach is to take S^⊤S as a new kernel matrix. The kernel is sometimes called empirical kernel mapping (EKM) [16]. The second approach is to subtract the smallest negative eigenvalue of the similarity matrix S from its diagonal. We call it the Saigo kernel [17]. We obtain a similarity matrix from distance matrix D _ij via S _ij = exp(-D _ij /σ ²).

Singular value decomposition of vectorial kernels

where U_{M × N}is M × N column-orthogonal matrix (columns are called left singular vectors) and ${\mathbf{V}}_{N\times N}^{\top }$ is N × N orthogonal matrix (rows are called right singular vectors). The matrix W_{N × N}is N × N diagonal matrix and 1/(N - 1)W² equals eigenvalues of the uncentered covariance matrix of A. We choose the largest q eigenvalues and replace other diagonal elements of W with 0, creating the W _q matrix. Finally, we obtain the denoised matrix, A _q , with A _q = U·W _q ·V^⊤. The vectorial data kernels are subsequently computed from the denoised matrix, A _q .

For actual analysis, software called SVDMAN developed by Wall et al. [27] is used.

Support vector machines

subject to

y _i (w·ϕ (x _i ) + b) ≥ 1 - ξ _i , ξ _i ≥ 0

for all N samples where y _i and C are the class label of i-th sample and a constant parameter, respectively. Optimization of this problem yields the hyper-plane that maximizes the margin between the two classes. This is called the SVM learning algorithm. The above optimization problem can be transformed into an equivalent form of the other equations in which only the kernel function K(x _i , x _j ) = ϕ(x _i )·ϕ(x _j ) appears and mapped vectors ϕ(x _j ) are not explicitly described [7]. Hence, the SVM learning algorithm needs only kernel matrix K and mapped vectors ϕ(x _j ) are not necessary.

Leave-one-out cross validation

where TP, FP, TN and FN are true positive, false positive, true negative, and false negative frequencies, respectively, in the classification.

Parameter selection

Since classification accuracies are dependent on parameters in the kernel-SVM method, we tested various parameter values to obtain the best performance possible. For all the six (linear, polynomial, RBF, EKM, Saigo, and ME) kernels tested here, seven SVM parameters, C = {10^-3, 10^-2, 10^-1, 1,10, 10², 10³}, are tested. For the polynomial kernel, D = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10} are tested. For the RBF, EKM, and Saigo kernels, σ = {10^-10, 10^-9, 10^-8,10^-7, 10^-6, 10^-5,10^-4, 10^-3, 10^-2,10^-1, 1} are tested. In the ME kernel, we used only one parameter G that magnifies the distance constraints D _ij to adjust the trade-off between over-learning and generalization of classification models (for details, see [21]). The parameter G has to be chosen carefully. When G → 0, typically K → 11^⊤/N. When D _ij > 2/N for ∀i, ∀j, K → I/N. The two are somewhat extreme cases. However, if the value of G is positive but too small, SVM cannot find the hyper-plane separating the positive class from the negative one clearly. If the value of G is too large, it leads to the so-called diagonal dominant problem [16]. We tested the parameter in the range of G = {2^-5, 2^-4, 2^-3, 2^-2, 2^-1, 1, 2, 2², 2³, 2⁴, 2⁵}. Note that the number of parameter combinations in the ME kernel is equal to those in the RBF, EKM and Saigo kernels in this study.