Abstract
Background
Human heart failure is a complex disease that manifests from multiple genetic and environmental factors. Although ischemic and nonischemic heart disease present clinically with many similar decreases in ventricular function, emerging work suggests that they are distinct diseases with different responses to therapy. The ability to distinguish between ischemic and nonischemic heart failure may be essential to guide appropriate therapy and determine prognosis for successful treatment. In this paper we consider discriminating the etiologies of heart failure using gene expression libraries from two separate institutions.
Results
We apply five new statistical methods, including partial least squares, penalized partial least squares, LASSO, nearest shrunken centroids and random forest, to two real datasets and compare their performance for multiclass classification. It is found that the five statistical methods perform similarly on each of the two datasets: it is difficult to correctly distinguish the etiologies of heart failure in one dataset whereas it is easy for the other one. In a simulation study, it is confirmed that the five methods tend to have close performance, though the random forest seems to have a slight edge.
Conclusions
For some gene expression data, several recently developed discriminant methods may perform similarly. More importantly, one must remain cautious when assessing the discriminating performance using gene expression profiles based on a small dataset; our analysis suggests the importance of utilizing multiple or larger datasets.
Background
Human heart failure is a complex disease diagnosed in over 500,000 American people every year, causing more than 250,000 deaths annually. It may arise from coronary atherosclerosis, exposure to toxins, infection, inflammation, valvular disease leading to volume/pressure overload, or an underlying genetic or idiopathic event [13]. Emerging work suggests the heterogeneity of heart failure. For example, patients with ischemic heart failure have decreased survival compared to the nonischemic heart failure [4,5] and respond differently to therapies [69]. Although benefits can be achieved using ischemic heart failure therapy for idiopathic heart failure (and vice versa), cure rates will be markedly diminished, and unwarranted toxicities problems will be encountered. It may be critical to distinguish these characteristically similar but clinically somewhat distinct heart failures, to better optimize therapy. The ability to distinguish between different etiologies of heart failure may be essential to guide appropriate therapy and determine prognosis for successful treatment.
A new approach to discriminating etiologies of heart failure is gene expression profiling using DNA microarray technology, which has been shown to be promising in the diagnosis of human diseases or subdiseases, especially in cancer [1012]. Recent genomic studies by three separate groups have demonstrated a distinct etiology dependent genomic pattern in the failing heart [1316]. These studies offer hope that the microarray gene expression analysis could potentially add to conventional laboratory approaches to diagnose different underlying etiologies of heart failure while simultaneously enhance prognostic criteria. It was hypothesized that heart failure arising from different underlying etiologies present with different gene expression patterns and that these differences could be used as a diagnostic tool. Here we test the hypothesis with two human heart failure datasets from different institutions.
Sample classification with gene expression data is statistically challenging due to the "small n, large p" problem [17]: the number of samples n is much smaller than the number of genes or predictors p. In our first dataset, we have n = 30 and p > 20000. Many new statistical methods have been developed or adapted to face the challenge. With more and more new methods emerging and existing methods being adapted, it becomes increasingly compelling for practitioners to compare and assess their performance, but there are few such comparative studies [1820]. Huang and Pan [19] compared several methods, including partial least squares (PLS) [21], nearest shrunken centroid (SC) [12], and a penalized PLS (PPLS), for binary classification of gene expression data. They found that these methods are competitive. More recently, some authors [22,20] have shown the promising performance of least absolute shrinkage and selection operator (LASSO) [23] and random forest (RF) [24]. It is our main goal to evaluate and compare these methods using two human heart failure datasets. For this purpose, we also extend the PPLS, originally proposed for binary classification, to multiclass classification. We found that the above five statistical methods perform similarly. Furthermore, our analysis stresses the importance of utilizing multiple datasets for classification purposes.
Results
Minnesota data
Myocardial tissue samples from the left ventricular apex of patients with severe refractory heart failure were collected at the time of the left ventricular assist device (LVAD) placement at the University of Minnesota Medical School. A total of 30 tissue samples were processed for microarray analysis on the Affymetrix Human Genome U133A chip containing ~22,000 genes. The initial data analysis was completed using Affymetrix Microarray Suite (MAS 5.0). A more complete description on the data is provided in [25].
The heart failure patients are divided into three classes according to the underlying etiology. Patients with clinical ECHO and EKG evidence, history of previous myocardial infarctions, and direct observation of the heart for confirmation of infarction at the time of LVAD implantation are defined as ischemic. Patients with an ischemic etiology were further divided into two classes: patients with ischemia but without acute myocardial infarctions (ischemic class) and patients with ischemia that have had an acute myocardial infarctions within ten days of LVAD implant (IM class), the remaining patients were assigned to the idiopathic class. Among the total 30 samples, 10 of them are ischemic, 7 are IM and 13 are idiopathic.
PGA data
The PGA data were obtained in another heart failure study conducted at the CardioGenomics PGA (Programs for Genomic Applications) at the Harvard Medical School. Myocardial samples were collected from patients undergoing heart transplantation whose failure arises from different etiologies (e.g. idiopathic, ischemic, alcoholic, valvular, and hypertrophic) and from normal organ donors whose hearts were not used for transplants. The transcriptional profile of the mRNA in these samples was also measured with Affymetrix oligonucleatide microarray technology. HGU133 plus 2 chips containing 54,675 probe sets were used and data were analyzed in MAS 5.0. In the PGA data set there were 11 normal samples, 11 ischemic samples and 14 idiopathic samples. The PGA dataset is publicly available at Genomics of Cardiovascular Development, Adaptation, and Remodelling. NHLBI Program for Genomic Applications, Harvard Medical School with URL: http://www.cardiogenomics.org webcite.
In order to make the results comparable to those based on the HGU133A chips used in the Minnesota data, we matched the probe sets on a HGU133 plus 2 chip with those on a HGU133A chip. Only six probe sets on the HGU133 plus 2 chip could not be found on a HGU133A chip. Hence we used the remaining 22277 ( = 22283  6) probe sets in the following analyses with the PGA data.
Classification with Minnesota data: Oneagainstothers approach
Table 1 reports the LOOCV misclassification errors of the five classification methods for the Minnesota data with models starting from different numbers of top ranked genes ranging from 50 to all genes. Four kinds of errors are reported for each method: three oneagainstothers twoclass LOOCV misclassification errors (Ischemic vs. others, IM vs. others and Idiopathic vs. others) and one threeclass LOOCV misclassification error. Note that throughout this article, threeclass classification results for SC and RF were obtained by direct applications of SC and RF, rather than by combining multiple binary classifications.
Table 1. LOOCV errors for threeclass classification with Minnesota data: all patients.
Based on Table 1, we first note that the performances of any classifier is sensitive to the number of top ranked genes one starts with. For example, in the threeclass classification, LASSO made 8 errors when only the top 200 genes are considered but made 18 errors when all 22283 genes are used. But there is no obvious relationship between the gene subset size and the five methods' performances.
In the twoclass classification problem of ischemic vs. others, PPLS, PLS, SC, LASSO and RF yielded LOOCV errors range from 7 to 15, 6 to 13, 8 to 13, 9 to 15 and 10 to 12 respectively. The five methods obtained similar numbers of errors in all instances. In the twoclass classification problems of IM vs. others and idiopathic vs. others, all five methods yielded similar numbers of errors. In the overall threeclass classification, again all the five classification methods perform similarly. There is no clear evidence that one method is clearly superior to others.
In order to assess whether a gene expression profile is affected by gender, we classified the 23 samples from male patients only. Among the 23 samples for male patients, 9 of them are ischemic, 7 are IM and 7 are idiopathic.
Table 2 reports the LOOCV misclassification errors of these five methods for the Minnesota data with males only. Again, the models start from different numbers of top ranked genes. The three twoclass LOOCV misclassification errors and one threeclass LOOCV misclassification error are estimated for each method as described before.
Table 2. LOOCV errors for threeclass classification with Minnesota data: males only.
Based on Table 2, again we note that the classification performances of PPLS, PLS, SC, LASSO and RF can be quite sensitive to the number of top ranked genes one uses. For example, in the oneagainstothers twoclass classification problem of ischemic vs. others, PPLS made 10 errors when the top 400 genes are considered but the number of errors suddenly increases to 14 when the top 800 genes are used. Again there is no obvious relationship between the gene subset size and the five methods' performances.
Comparing PPLS, PLS, SC, LASSO and RF, we find that the five methods perform very similarly in almost all instances. There is no evidence that any method is clearly the best. One thing we noticed about LASSO is that when the model contains many genes, say top 12800, then it gives more errors than PPLS, PLS, SC and RF in the threeclass classification. This could happen by chance since we did not observe the same trend in the decomposed oneagainstothers binary classifications.
As in Table 1, the results in Table 2 suggest that discriminating ischemic group from the other two groups was less accurate than distinguishing IM from the other two groups and separating idiopathic group from the other two groups.
If we compare Table 1 (30 patients) with Table 2 (23 male patients), we can see that the misclassification error rates of all the five methods in Table 2 are much higher than those in Table 1. Reduced sample size is likely a factor.
To see whether the high prediction errors are due to the presence of the three classes, we considered a binary classification problem. We applied all the five methods to the 10 ischemic and 13 idiopathic samples. We also assessed the classification accuracy on the male patients with 9 ischemic and 7 idiopathic samples. The classification results were shown in Table 3.
Table 3. LOOCV errors for twoclass classification with Minnesota data: ischemic vs idiopathic.
From Table 3, we can see that all the five methods have very similar performances in classifying ischemic and idiopathic samples. If we compare the classification performances of these five methods with/without female patients, taking the sample size into consideration, we can see that the misclassification error rates with only male patients is much higher than those with all patients. This again is probably because the sample size (16) with only males is smaller. In particular, we noticed that LASSO is more sensitive to the small sample size.
Other models
In the previous classification problems, we only included linear terms of gene expression levels in a model. We also considered expanded models including squared terms of each gene's expression levels. The motivation is to possibly improve model fitting, for instance, to avoid masking in linear models [26]. In this way, the number of variables in the new data is doubled (the original variables plus their squared terms) and we have 44566 variables. We repeated all the previous classification procedures and found that the classification performance did not improve (results not shown).
Classification with Minnesota data: pairwise approach
We repeated the threeclass classification with PPLS via the pairwise approach. The results are included in Table 4. We assessed the classification of PPLS by including all 30 patients and with 23 male only. By comparing Table 4 to Table 1 – 2, we can see that the PPLS via oneagainstothers approach gives much smaller errors. This suggests that for this specific problem, the oneagainstothers approach is probably better. Again, we see the classification with male patients gives much larger LOOCV misclassification error rates.
Table 4. LOOCV errors for threeclass classification: PPLS with pairwise approach.
Classification with PGA data: oneagainstothers approach
Table 5 reports the LOOCV misclassification errors of the five methods for the PGA data. Four kinds of errors are reported for each method: three oneagainstothers twoclass LOOCV misclassification errors (Normal vs. others, Ischemic vs. others, and Idiopathic vs. others) and one threeclass LOOCV misclassification error.
Table 5. LOOCV errors for threeclass classification with PGA data: all patients.
Based on Table 5, we can see that the classification performances of PPLS, PLS, SC, LASSO and RF are quite stable with different numbers of top ranked genes one uses. In the twoclass classification problem of normal vs. others, the five classification methods almost perform perfectly, where PPLS, PLS and RF have 0 errors in all circumstances, SC has 0 errors in all situations except for 3 cases and LASSO has 1 error in one case and 0 errors in all other cases. In the twoclass classification problems of ischemic vs. others and idiopathic vs. others, PPLS, PLS, SC, LASSO and RF yielded 1–5 errors (mostly with 1–3 errors). That the problem of distinguishing normal from the others is the easiest confirms that the normal class is more separable from the other two classes. As for the threeclass classification, the errors range from 1 to 2 and it is almost perfect.
It may be argued that a classification involving normal samples should be much easier because the normal class is very different from the other two classes. Correct diagnosis between ischemic and idiopathic would be much more challenging. Hence we conducted a twoclass classification with 11 ischemic and 14 idiopathic samples. The results are shown in Table 6. The five methods perform almost perfectly: the misclassification errors range from 1 to 3 for all five methods in all models.
Table 6. LOOCV errors for twoclass classification with PGA data: ischemic vs idiopathic.
Genes identified
We consider genes remaining in a final model for each method. To save space, we restrict attention to models starting with all the genes for binary discrimination between ischemic and idiopathic samples. Briefly, LOOCV was first used to select any tuning parameters in a method (e.g. number of components in a PLS model), then a model with the selected parameters was fitted using all the samples. Except that all the genes are used in a final PLS model, for any of the other methods there may be fewer genes remaining in the final model. In particular, LASSO can select at most n genes with n the number of the samples. It turned out that random forest also used many of the genes.
Tables 7 and 8 lists the genes selected by at least four or three methods for the Minnesota data and PGA data respectively. It can be seen that there is no overlap at all between the two gene lists. Although the same genes are not identified from the two datasets, it is clear that the betaadrenergic signalling pathway is likely a discriminatory pathway, given the inclusion of CREM in the Minnesota data and AKAP6 in the PGA data. Furthermore, the inclusion of metabolicrelated genes, such as ATPase and GAPD, is not surprising given the class of ischemic tissue.
Table 7. Genes selected by ≥ 4 methods in twoclass classification (ischemic vs idiopathic) with Minnesota data. The last row gives the total numbers of the genes in the final models.
Table 8. Genes selected by ≥ 3 methods in twoclass classification (ischemic vs idiopathic) with PGA data, The last row gives the total numbers of the genes in the final models.
We also give the univariate ranks of the genes (based on the Fstatistics for the two classes) in the above two tables. It shows that the two sets of genes (or more generally, the genes in a final model) may not include some genes ranked high in the univariate ranking while including some ranked low, highlighting a possible limitation of solely depending on univariate ranks to select important genes.
More numerical results
PLS plots
To further explore why the methods all work much better for the PGA data than for the Minnesota data, we drew some plots using the first and the second PLS components for binary classification of ischemic vs idiopathic groups. We found that, for both datasets, there was a clear separation between the two groups. However, in LOOCV, although again the two groups were separable for both datasets, the leftout sample was more likely to be closer to the other group than to its true group for the MInnesota data, leading to a high LOOCV error rate; Figure 1 gives two examples. In contrast, in the PGA data, a leftout sample tended to be close to its true group, resulting in a low LOOCV error rate. For more details see Supplemental Materials.
Figure 1. PLS plots for two cases in LOOCV for the Minnesota data comparing ischemic vs idiopathic. In both cases, the new sample labeled as "N" (i.e. leftout sample in LOOCV), belonging to class 1 and 2 respectively in the top and the bottom panels, is closer to the other group different from its true class, leading to misclassifications.
Permutation tests
Because of high misclassification error rates with the Minnesota data, it is of interest to investigate whether there is any signal at all in the data. This can be accomplished by a permutation test that compares misclassification errors resulting from using the original data with that from randomly permuted data; a Pvalue is defined as the proportion of permuted datasets with misclassification errors fewer than that of the original data. To generate a randomly permuted data, we randomly permute the class labels of the original data. Because all the methods have similar performance, we consider the nearest shrunken centroid method with the Minnesota data. Tables 9 and 10 summarize the results of misclassification errors for 50 randomly permuted datasets for three and twoclass classifications respectively. It can be seen that the misclassification errors based on the original data are fewer than that based on the permuted data, leading to small Pvalues. This implies that, although there are relatively high misclassification error rates with the Minnesota data, the methods perform significantly better than a random guess.
Table 9. LOOCV threeclass misclassification errors with the original Minnesota data and the percentiles of LOOCV errors with 50 permuted datasets by SC.
Table 10. LOOCV twoclass misclassification errors with the original Minnesota data and the percentiles of LOOCV errors with 50 permuted datasets by SC: ischemic vs idiopathic.
Simulations
We did a simulation study to further evaluate the performance of the five classification methods. To mimic the real data, simulated data were generated from either a fitted PPLS or a fitted LASSO model to the Minnesota data comparing ischemic vs idiopathic, each containing top 50 genes in the initial model. Specifically, suppose that is the fitted response value for sample i based on the original Minnesota data using PPLS or LASSO. Note that is a real number without being dichotomized yet. Suppose that Y_{i }= 1 or 1 is the class label of sample i in the original data, and . To generate a simulated data, we independently draw from Normal distributions for i = 1, ..., 23 and b = 1, ..., 50. Then we apply each method to a simulated dataset , where X_{i }is the gene expression profile of sample i in the original Minnesota data, and obtain fitted values ; the resulting misclassification error number for dataset b is .
Table 11 summarizes the distributions of the misclassification errors of each method based on 50 simulated data with either PPLS or LASSO as the true model. It can be seen that in general all the methods perform similarly, though random forest seems to be most stable and has a slight edge, and the performance of LASSO and nearest shrunken centroid may deteriorate as the number of the genes included in a model is increased. We also did other simulations with the true models starting from various numbers of top genes and various noise levels, and observed similar phenomena: for details see Supplemental Materials.
Table 11. Percentiles of misclassification errors from 50 simulated datasets for twoclass classification. The true model is either PPLS or LASSO fitted with top 400 genes to the Minnesota data to compare ischemic vs idiopathic.
Discussion
With more and more statistical methods being proposed for discriminant analysis for gene expression data, it has become increasing important to compare and evaluate their performances with real data, as it has been done in other contexts [27]. Comparing the five new methods with each other using the two real datasets, we did not find anyone uniformly better than the others. This may be disappointing to someone who wishes to find the best statistical method. However, in the current application, the similar performance of all the five methods on each of the two datasets provides reassurance on the interesting observation that it is not equally easy to distinguish the different etiologies of heart failure using expression profiles in the two datasets.
Both the oneagainstothers approach and the pairwise approach have been widely used in extending a binary classifier to multiclass settings. Our result suggests that, at least for the two datasets used here, the oneagainstothers approach is better, which was found to be true with support vector machines but in general should also depend on which binary classifier is used [28]. We also have observed that any of the five methods may be sensitive to the number of genes being included. This is particularly relevant because, although all the five methods (and many other methods) can handle any large number of genes, this does not dismiss the potential importance of a user's preliminary ranking and screening of genes. Of course, all our observations here are based on the two datasets without consideration of statistical variability, further studies are needed to validate these points.
An interesting finding of this work is that it is difficult to discriminate the different etiologies of human heart failure using one gene expression dataset, and at the same time, it is quite easy for the other dataset. A possible explanation is the different types of the microarray chips used: Affymetrix HGU133A chips were used in the Minnesota study while Affymetrix HGU133 plus 2 chips were used in the PGA study. Because the HGU133 plus 2 chips contain more genes (or ESTs), to minimize the effects of using different genes, we only used the genes present in the Minnesota data and still yielded much better performance for the PGA data. In fact, we used all the genes in the PGA data and obtained similar results for the PGA data. Although we can say that the performance difference in the two datasets is not caused by different genes contained on a chip, we do not know whether the more recent HGU133 plus 2 chips provide more reliable measurements on gene expression. In addition, quality control criteria for the inclusion of a chip were nearly identical between the two datasets. We would suspect that the performance difference may be the result of different patient populations and different study protocols (e.g. lack of clearly prespecified patient inclusion/exclusion criteria). As discussed in [29], a key to validating any prognostic and diagnostic biomarkers is the use of data that can reflect the full range of clinical variability. This highlights the importance of utilizing multiple datasets drawn from multiple subpopulations. Even for the purpose of prediction for one subpopulation, it is possible to improve the performance by borrowing information from other subpopulations [30]. It can be argued that the performance should be weighted on the complexity of the disease. Challenges with the current clinical discrimination of ischemic versus nonischemic heart failure is indeed why defining potential gene expression biomarkers may be a helpful additional approach in this characterization. A recognized limitation of utilizing heart tissue to identify biomarkers is the difficulty of collecting tissue. In summary, the current and other studies stress the importance of collaborating efforts to share tissue/data to strengthen the search for applicable biomarkers.
Conclusions
Many studies have aimed to develop new statistical and machine learning methods for best sample discrimination. Our results suggest that, at least for some gene expression data, several existing methods may work almost equally well. More importantly, because of the quite different performances of the methods on the two datasets, one must remain cautious when assessing the performance of sample discrimination using a small gene expression dataset; it may be necessary to use larger or multiple datasets to draw a more reliable conclusion.
Methods
Binary classifiers: PLS, PPLS and LASSO
We first briefly review the three binary classifiers, which was first designed for regression and can be directly applied to twoclass classification, even when the number of covariates (i.e. genes here) is much larger than the sample size.
We code the response variable (i.e. class label) as Y = 1 for class 1 and Y = 1 for class 2. Suppose that x_{i }is the expression level of gene i, i = 1, ..., p with p as the total number of the genes, and that we have n samples in the training data. A challenge is that we have n <<p.
The main idea of partial least squares (PLS) [21] is to seek a few linear combinations of for j = 1, ..., m, then apply ordinary least squares (OLS) to regress Y on z_{j}'s to obtain
with β's as OLS estimates. The key of course is how to form linear components z_{j}'s. It turns out that
α_{j }= argmax_{α}Corr^{2}(y, Xα) Var(Xα)
with the constraints α = 1, for l = 1,..., j  1, where y is the vector of observed Y's (in the training data), X is the design matrix (i.e. matrix of observed x's), and S is the sample covariance matrix of x's [31]. In practice, the number of linear components m has to be chosen, typically by a form of crossvalidation, such as leaveoneout crossvalidation (LOOCV), to minimize misclassification errors.
PPLS is a penalized regression method in the framework of PLS [19,32]. Suppose that we have built a PLS linear model, which can be rewritten as:
Then we shrink the PLS coefficients by softthresholding [33,34]
where sign(a) = 1 if a ≤ 0 and sign(a) = 1 if a < 0, λ is a shrinkage parameter to be determined, and f_{+ }= max(f, 0). It is common that the shrinkage leads to many , effectively eliminating gene i from the model, thus gene selection is automatically accomplished. Next we construct a linear component . Finally a PPLS model is built by regressing Y on z using OLS
which can be reexpressed as . The parameters involved in building a PPLS model, such as the shrinkage parameter λ and the number of PLS components, are estimated by LOOCV. The goal is to choose the largest shrinkage parameter and the smallest number of PLS components for which the LOOCV misclassification error estimate is minimized.
The LASSO estimates [23] in a linear model
are obtained by
subject to , where Y_{i }is the observed response for sample i and is its LASSO estimate, i = 1, ..., n, and t can be chosen by LOOCV. The constraint can often force many , leading to gene selection.
Note that the class label (1/1) for the response Y is binary, but in any of the above binary classifiers, the response Y is treated as a continuous variable and the estimate could be any real number. To predict the class of a new sample, we use sign(): if the estimated response is greater than or equal to 0, then we classify it into class 1; otherwise, class 2. In particular, this direct use of PLS for binary classification (as in [35]) is different from other approaches [3640]; a distinct advantage of our approach is its simplicity, e.g., avoiding convergence problems when two classes are perfectly separable, which is common in microarray data with a small sample size and a large number of genes.
Multiclass classifiers: nearest shrunken centroids and random forest
Nearest shrunken centroids (SC) is built on a diagonalized linear discriminant analysis (DLDA) [26,41]. Suppose that we have K classes, is the mean expression level of gene i in class k of the training samples, is the pooled sample variance of gene i of the training samples, and π_{k }is the prior probability of a new sample being in class k. The DLDA rule for a new sample is
SC is motivated from the observation that many of the genes will not be predictive of the class membership and should be eliminated from the above DLDA rule. Formally, define
where n_{k }is the number of training samples in class k, and is the overall mean expression level of gene i in all the training samples. Note that by the definition, we have . Let
for all i and k, where Δ is the shrinkage parameter to be chosen by LOOCV. Then substituting in the DLDA rule by , we obtain a SC rule
The new sample x* is assigned to class k_{0 }such that .
Note that, if , then and thus gene i plays no role in classifying for class k. Hence SC effectively accomplishes gene selection by shrinkage.
Random forest (RF) [24] is an ensemble of classification trees [42,43], which have been shown to be useful in tumor classification with microarray data [44]. It is designed to improve over a single classification tree. There are two random aspects that help generate multiple classification trees in RF. First, a bootstrap sample is repeatedly drawn from the original training data and then used to build a classification tree. Second, in building a classification tree, rather than using the best splitting variable (i.e. gene here) from all the available variables at each node, it chooses the best from a small random subset of all the variables. Each tree is grown to the maximum and no pruning is pursued. To predict the class for a new sample, the sample is applied to each tree and each tree votes by giving its prediction, then the majority vote is taken as the final prediction for the sample.
Extending a binary classifier to multiclass classification
Here we describe how a multiclass (K > 2) classification problem can be handled by a binary classifier, such as PLS, PPLS and LASSO. It is achieved by formulating a multiclass classification problem as multiple twoclass classification problems. We consider two most popular approaches: one is to compare each class against all the others, and the other is to compare all possible pairs of classes. Applications of these two approaches can be found, among others, in [4550]. In particular, some have considered the first approach for PLS [50].
The oneagainstothers approach is to reduce a Kclass classification task to K twoclass classification problems. Formally, a new response is defined in the k_{th }binary problem as:
for k = 1, ..., K. Then we build K binary classifiers. To predict a new sample with gene expression profile x*, we apply x* to each binary classifier and yield . Finally, the class of the new sample C(x*) is predicted as
That is, the new sample is classified into the class maximizing .
The pairwise approach reduces a Kclass classification to K(K  1)/2 twoclass classification problems [45]. Specifically, for each of all possible pairs of classes, solve each of the twoclass problems and then, for a new sample, combine all the pairwise decisions to form a Kclass decision. Suppose that the new binary response in a pairwise comparison with classes k_{1 }and k_{2 }(with 1 ≤ k_{1 }<k_{2 }≤ K) is defined as
We build a binary classifier with response using only samples belonging to class k_{1 }or k_{2}, and denote the fitted response value for a new sample with expression profile x* as . As described earlier, we classify the new sample into class k_{1 }or k_{2 }according to the sign of . After this is done for any 1 ≤ k_{1 }<k_{2 }≤ K, the final decision is to assign the new sample to the class that wins the most pairwise comparisons. In the case when there are multiple winning classes, we randomly pick one of the winning classes to be the final winning class. Comparing to the oneagainstothers approach, the pairwise approach is computationally more expensive if K ≥ 4.
Gene ranking
To explore the effect of the number of genes a model starts with on the classification performance, we have a preliminary gene ranking using a usual Fstatistic. This univariate ranking is used throughout, and obviously is by no means to be optimal. For the purpose of the presentation in this section, we only need to consider a given gene. Suppose x_{ik }is the gene expression level of the gene in sample i that is in class k, i = 1, ..., n_{k}, and k = 1, ..., K, where n_{k }is the total number of samples in class k and K is the total number of classes. Let be the mean expression level of classk samples, be the overall mean (across all the samples) and be the total number of samples. We can construct an Fstatistic as the ratio of the mean sums of squares for betweenclass and withinclass variations:
We can rank all the genes based on their corresponding Fstatistics: a gene with a larger Fstatistic indicates a stronger relationship between its expression levels and the class membership in the samples, and therefore has a higher rank as a potential predictor of the class. We started with various models by including different numbers of top ranked genes. We considered models starting from the top 50, 100, 200, 400, 800, 1600, 3200, 6400, 9200, 12800, 16000, 19200, and all (22283 and 22277 for the two datasets respectively) genes respectively.
It is an incorrect practice in microarray experiments to first select genes using all the samples and then perform crossvalidation using the selected genes, which gives downward biased prediction error estimates [51,52]. Hence, it is essential to perform crossvalidation on the entire model building process, including gene selection. In our study, we did honest crossvalidation. In particular, we crossvalidated gene selection (and other aspects of model building, such as parameter selection and estimation). Specifically, in LOOCV, we remove each sample from the data in turn (which is then treated as the test sample), carry out gene selection using Fstatistic based on the remaining samples, build a classifier with the selected genes using the remaining samples, and then test the classifier on the leftout sample.
Data preprocessing
To facilitate the application of penalized regression (i.e. PPLS and LASSO) so that their regression coefficients are in the same unit and thus can be penalized using a global penalty parameter, the expression levels of each gene were scaled to have sample variance 1.
Evaluations
In addition to PLS/PPLS, we will consider the shrunken centroids (SC) method, the LASSO, and the random forest (RF). SC, LASSO and RF have been implemented in R [53], and are easy to use; we applied their R functions using default parameter settings. SC and RF are directly applicable to multiclass classification while LASSO, as PLS/PPLS, is itself a binary classifier. For multiclass classification with PLS and LASSO, we used the same approaches as described for PPLS.
We use the leaveoneout crossvalidation (LOOCV) to estimate the prediction error for each of the methods. Within this firstlevel LOOCV, a secondlevel LOOCV is used to select tuning parameters for each method to minimize cross validation errors. Specifically, in PLS, the smallest number of PLS components is selected among the PLS models that give the minimum LOOCV error. In PPLS, among the models with minimum LOOCV error, we first pick the ones with the smallest number of PLS components, then pick the one with the largest shrinkage parameter. In the SC method, the largest shrinkage is selected among the models that minimize the LOOCV error. The number of candidate threshold values and the number of cross validation folds are both set to be default (i.e. 30 and the smallest class size respectively). In LASSO, the maximum fraction parameter of the models that minimize LOOCV error is selected while the number of the candidate fraction values is set to be 51 (equally spaced from 0 to 1) and the number of cross validation folds is set to be the total sample size. In RF, every parameter is set to be default. For example, the number of trees to grow is set to 500, and the number of candidate splitting variables considered at each split is set to by default, where p is the total number of variables (i.e. genes).
Due to the small sample size (about 10 in each class) in each dataset, it is quite challenging to estimate the prediction error well. Although it is straightforward to apply LOOCV or other crossvalidation methods, their performance may not be optimal. After submitting this work, we became aware of the recent work by Fu et al [54], where a better method than LOOCV was proposed specifically for microarray data. This new method aims to reduce the variability of LOOCV. We reason that with the use of this new method, the main conclusions drawn in this work would not change.
Authors' contributions
XH downloaded the PGA data, did the analysis and simulations. WP conceived of and directed the study. SG cleaned and managed the Minnesota data. XH, YC, SJP, LWM and JH conducted the Minnesota study and generated the data. XH, WP and JH drafted the manuscript.
Additional File 1. PLS plots with the first two PLS components. PLS plot for the Minnesota data comparing ischemic vs idiopathic groups.
Format: PS Size: 6KB Download file
Additional File 2. PLS plots with the first two PLS components. PLS plot for the PGA data comparing ischemic vs idiopathic groups.
Format: PS Size: 6KB Download file
Additional File 3. PLS plots with the first two PLS components. PLS plots for LOOCV for the Minnesota data comparing ischemic vs idiopathic groups.
Format: PS Size: 48KB Download file
Additional File 4. PLS plots with the first two PLS components. PLS plots for LOOCV for the PGA data comparing ischemic vs idiopathic groups.
Format: PS Size: 52KB Download file
Additional File 5. Simulation results with various simulation setups.
Format: PS Size: 39KB Download file
Acknowledgements
XH and WP were supported by NIH grant HL65462. JH was supported by an AHA grant, the Lillehei Heart Institute and the Minnesota Medical Foundation. The authors are grateful to two reviewers for helpful and constructive comments; in particular, section "Other numerical results" was added based on the reviewers' comments.
References

Levy D, Larson MG, Vasan RS, Kannel WB, Ho KK: The progression from hypertension to congestive heart failure.
JAMA 1996, 275:15571562. PubMed Abstract  Publisher Full Text

LloydJones DM: The risk of congestive heart failure: sobering lessons from the Framingham Heart Study.
Curr Cardiol Rep 2001, 3:184190. PubMed Abstract  Publisher Full Text

Nicol RL, Frey N, Olson EN: From the sarcomere to the nucleus: role of genetics and signaling in structural heart disease.
Annu Rev Genomics Hum Genet 2000, 1:179223. PubMed Abstract  Publisher Full Text

Felker GM, Thompson RE, Hare JM, Hruban RH, Clemetson DE, Howard DL, Baughman KL, Kasper EK: Underlying causes and longterm survival in patients with initially unexplained cardiomyopathy.
New Engl J Med 2000, 342:10771084. PubMed Abstract  Publisher Full Text

Dries DL, Sweitzer NK, Drazner MH, Stevenson LW, Gersh BJ: Prognostic impact of diabetes mellitus in patients with heart failure according to the etiology of left ventricular systolic dysfunction.
J Am Coll Cardiol 2001, 38:421428. PubMed Abstract  Publisher Full Text

Kittleson M, Hurwitz S, Shah MR, Nohria A, Lewis E, Givertz M, Fang J, Jarcho J, Mudge G, Stevenson LW: Development of circulatoryrenal limitations to angiotensinconverting enzyme inhibitors identifies patients with severe heart failure and early mortality.
J Am Coll Cardiol 2003, 41:20292035. PubMed Abstract  Publisher Full Text

Felker GM, Benza RL, Chandler AB, Leimberger JD, Cuffe MS, Califf RM, Gheorghiade M, O'Connor CM: Heart failure etiology and response to milrinone in decompensated heart failure: results from the OPTIMECHF study.
J Am Coll Cardiol 2003, 41:9971003. PubMed Abstract  Publisher Full Text

Doval HC, Nul DR, Grancelli HO, Perrone SV, Bortman GR, Curiel R: Randomised trial of lowdose amiodarone in severe congestive heart failure.
Lancet 1994, 344:493498. PubMed Abstract  Publisher Full Text

Singh SN, Fletcher RD, Fisher SG, Singh BN, Lewis HD, Deedwania PC, Massie BM, Colling C, Lazzeri D: Amiodarone in patients with congestive heart failure and asymptomatic ventricular arrhythmia.
New Engl J Med 1995, 333:7782. PubMed Abstract  Publisher Full Text

Golub TR, Slonim DK, Tamayo P, Huard C, Gaasenbeek M, Mesirov JP, Coller H, Loh ML, Downing JR, Caligiuri MA, et al.: Molecular classification of cancer: Class discovery and class prediction by gene expression monitoring.
Science 1999, 285:531537. PubMed Abstract  Publisher Full Text

Hedenfalk I, Duggan D, et al.: Geneexpression profiles in hereditary breast cancer.
New England Journal of Medicine 2001, 344:539548. PubMed Abstract  Publisher Full Text

Tibshirani R, Hastie R, Narasimhan B, Chu G: Diagnosis of multiple cancer types by shrunken centroids of gene expression.
PNAS 2002, 99:65676572. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Aronow BJ, Toyokawa T, Canning A, Haghighi K, Delling U, Kranias E, Molkentin JD, Dorn GW: Divergent transcriptional responses to independent genetic causes of cardiac hypertrophy.
Physiol Genomics 2001, 6:1928. PubMed Abstract  Publisher Full Text

Tan FL, Moravec CS, Li J, AppersonHansen C, McCarthy PM, Young JB, Bond M: The gene expression fingerprint of human heart failure.
Proc Natl Acad Sci 2002, 99:1138711392. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hwang JJ, Allen PD, Tseng GC, Lam CW, Fananapazir L, Dzau VJ, Liew CC: Microarray gene expression profiles in dilated and hypertrophic cardiomyopathic endstage heart failure.
Physiol Genomics 2002, 10:3144. PubMed Abstract  Publisher Full Text

Kittleson M, Ye SQ, Irizarry RA, Minhas KM, Edness G, Conte JV, Parmigiani G, Miller LW, Chen Y, Hall JL, Garcia JGN, Hare JM: Identification of a gene expression profile that differentiates ischemic and nonischemic cardiomyopathy.
Circulation 2004, 110:344451. PubMed Abstract  Publisher Full Text

West M: Bayesian Factor Regression Models in the "Large p, Small n" Paradigm.

Dudoit S, Fridlyand J, Speed T: Comparison of discrimination methods for the classification of tumors using gene expression data.
J Am Stat Assoc 2002, 97:7787. Publisher Full Text

Huang X, Pan W: Linear regression and twoclass classification with gene expression data.
Bioinformatics 2003, 19:20722078. PubMed Abstract  Publisher Full Text

Wu BL, Abbott T, Fishman D, McMurray W, Mor G, Stone K, Ward D, Williams K, Zhao HY: Comparison of statistical methods for classification of ovarian cancer using mass spectrometry data.
Bioinformatics 2003, 19:16361643. PubMed Abstract  Publisher Full Text

Wold S, Ruhe A, Wold H, Dunn WJ III: The collinearity problem in linear regression. The partial least squares (PLS) approach to generalized inverses.
SIAM J of Scientific and Statistical Computing 1984, 5(3):735742. Publisher Full Text

Segal MR, Dahlquist KD, Conklin BR: Regression approaches for microarray data analysis.
J Comp Biol 2003, 10:961980. PubMed Abstract  Publisher Full Text

Tibshirani R: Regression shrinkage and selection via the lasso.
Journal Royal Statistical Society, Series B 1996, 58:267288.

Machine Learning 2001, 45(1):532. Publisher Full Text

Hall JL, Grindle S, Han X, Fermin D, Park S, Chen Y, Bache RJ, Mariash A, Guan Z, Ormaza S, Thompson J, Graziano J, de Sam Lazaro SE, Pan S, Simari RD, Miller LW: Genomic profiling of the human heart before and after mechanical support with a ventricular assist device reveals alterations in vascular signaling networks.
Physiological Genomics 2004, 17:283291. PubMed Abstract  Publisher Full Text

Hastie T, Tibshirani R, Friedman J: The Elements of Statistical Learning. Data mining, Inference, and Prediction. Springer; 2001.

Titterington DM, Murray GD, Murray LS, Spiegelhalter DJ, Skene AM, Habbema JDF, Gelpke GJ: Comparison of discrimination techniques applied to a complex data set of head injured patients (with discussion).
Journal Royal Statistical Society, Series A 1981, 144:145175.

Rifkin R, Mukherjee S, Tamayo P, Ramaswamy S, Yeang CH, Angelo M, Reich M, Poggio T, Lander ES, Golub TR, Mesirov JP: An analytical method for multiclass molecular cancer classification.
SIAM Review 2003, 45:706723. Publisher Full Text

Simon R: When is a genomic classifier ready for prime time.
Nature Clinical Practice – Oncology 2004, 1:45. Publisher Full Text

Huang X, Pan W, Han X, Chen Y, Miller LW, Hall J: Borrowing information from relevant microarray studies for sample classification using weighted partial least squares.
Computational Biology and Chemistry 2005, 29:204211. Publisher Full Text

Frank IE, Friedman JH: A statistical view of some chemometrics regression tools (with discussion).

Huang X, Pan W, Park S, Han X, Miller LW, Hall J: Modeling the relationship between LVAD support time and gene expression changes in the human heart by penalized partial least squares.
Bioinformatics 2004, 20:888894. PubMed Abstract  Publisher Full Text

Donoho DL, Johnstone IM: Ideal spatial adaptation by wavelet shrinkage.

Donoho DL: DeNoising by SoftThresholding.
IEEE Transaction on Information Theory 1995, 41:613627. Publisher Full Text

Hawkins DM, Wolfinger RD, Liu L, Young SS: Exploring blood spectra for signs of ovarian cancer.

Nguyen DV, Rocke DM: Tumor classification by partial least squares using microarray gene expression data.
Bioinformatics 2002, 18:3950. PubMed Abstract  Publisher Full Text

Ghosh D: Singular value decomposition regression models for classification of tumors from microarray experiments.

Ghosh D: Penalized discriminant methods for the classification of tumors from gene expression data.
Biometrics 2003, 59:9921000. PubMed Abstract  Publisher Full Text

Ding B, Gentleman R: Classification using generalized partial least squares. [http://www.bepress.com/bioconductor/paper5] webcite
Technical Report 5, Bioconductor Project Working Papers 2004.

Fort G, LambertLacroix S: Classification using partial least squares with penalized logistic regression.
Bioinformatics 2005, 21:11041111. PubMed Abstract  Publisher Full Text

McLachlan GJ: Discriminant Analysis and Statistical Pattern Recognition. New York: Wiley; 1992.

Breiman L, Friedman JH, Olshen RA, Stone CJ: Classification and Regression Trees. Belmont: Wadsworth; 1984.

Zhang H, Singer B: Recursive Partitioning in the Health Sciences. SpringerVerlag: New York; 1999.

Zhang H, Yu CY, Singer B, Xiong M: Recursive partitioning for tumor classification with gene expression microarray data.
PNAS 2001, 98:67306735. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Friedman J: Another approach to polychotomous classification.

Hastie T, Tibshirani R: Classification by pairwise coupling.
Annals of Statistics 1998, 26:451471. Publisher Full Text

Allwein E, Schapire R, Singer Y: Reducing multiclass to binary: A unifying approach for margin classifiers.
Journal of Machine Learning Research 2000, 1:113141. Publisher Full Text

Ramaswamy S, Tamayo P, Rifkin R, Mukherjee S, Yeang CH, Angelo M, Ladd C, Reich M, Latulippe E, Mesirov JP, Poggio T, Gerald W, Loda M, Lander ES, Golub TR: Multiclass cancer diagnosis using tumor gene expression signatures.
PNAS 2001, 98:1514915154. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Dettling M, Buhlmann P: Boosting for tumor classification with gene expression data.
Bioinformatics 2003, 19:10631069. Publisher Full Text

Tan Y, Shi L, Tong W, Hwang GTG, Wang C: Multiclass tumor classification by discriminant partial least squares using microarray gene expression data and assessment of classification models.
Computational Biology and Chemistry 2004, 28:235243. Publisher Full Text

Ambroise C, McLachlan GJ: Selection bias in gene extraction on the basis of microarray geneexpression data.
PNAS 2002, 99:65626566. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Simon R, Radmacher MD, Dobbin K, McShane LM: Pitfalls in the use of DNA microarray data for diagnostic and prognostic classification.
J Natl Cancer Inst 2003, 95:1418. PubMed Abstract  Publisher Full Text

Ihaka R, Gentleman R: R: A language for data analysis and graphics.
Journal of Computational and Graphical Statistics 1996, 5:299314.

Fu WJ, Carroll RJ, Wang S: Estimating misclassification error with small samples via bootstrap crossvalidation.
Bioinformatics 2005, 21:19791986. PubMed Abstract  Publisher Full Text