Bioinformatics

Bioinformatics Advance Access originally published online on December 13, 2005
Bioinformatics 2006 22(4):472-476; doi:10.1093/bioinformatics/bti827

This Article Abstract FREE Full Text (Print PDF) Supplementary Material A corrigendum has been published All Versions of this Article: 22/4/472    most recent bti827v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Search for citing articles in: ISI Web of Science (7) Request Permissions Google Scholar Articles by Wu, B. Search for Related Content PubMed PubMed Citation Articles by Wu, B. Social Bookmarking          What's this?

© The Author 2005. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oxfordjournals.org

Differential gene expression detection and sample classification using penalized linear regression models

Baolin Wu

Division of Biostatistics, School of Public Health, University of Minnesota A460 Mayo Building, MMC 303, Minneapolis, MN 55455, USA

	ABSTRACT

TOP ABSTRACT 1 INTRODUCTION 2 PENALIZED LINEAR REGRESSION 3 APPLICATION TO MICROARRAY... 4 COMPARISON WITH OTHER... 5 DISCUSSION REFERENCES

 

Differential gene expression detection and sample classification using microarray data have received much research interest recently. Owing to the large number of genes p and small number of samples n (p >> n), microarray data analysis poses big challenges for statistical analysis. An obvious problem owing to the ‘large p small n’ is over-fitting. Just by chance, we are likely to find some non-differentially expressed genes that can classify the samples very well. The idea of shrinkage is to regularize the model parameters to reduce the effects of noise and produce reliable inferences. Shrinkage has been successfully applied in the microarray data analysis. The SAM statistics proposed by Tusher et al. and the ‘nearest shrunken centroid’ proposed by Tibshirani et al. are ad hoc shrinkage methods. Both methods are simple, intuitive and prove to be useful in empirical studies.

Recently Wu proposed the penalized t/F-statistics with shrinkage by formally using the ₁ penalized linear regression models for two-class microarray data, showing good performance. In this paper we systematically discussed the use of penalized regression models for analyzing microarray data. We generalize the two-class penalized t/F-statistics proposed by Wu to multi-class microarray data. We formally derive the ad hoc shrunken centroid used by Tibshirani et al. using the ₁ penalized regression models. And we show that the penalized linear regression models provide a rigorous and unified statistical framework for sample classification and differential gene expression detection.

Availability: For the computer programs, detailed analysis results and R functions for the proposed methods, please see http://www.biostat.umn.edu/~baolin/research/L1C-mc.html

Contact: baolin{at}biostat.umn.edu

Supplementary information: http://www.biostat.umn.edu/~baolin/research/L1C-mc.html

	1 INTRODUCTION

TOP ABSTRACT 1 INTRODUCTION 2 PENALIZED LINEAR REGRESSION 3 APPLICATION TO MICROARRAY... 4 COMPARISON WITH OTHER... 5 DISCUSSION REFERENCES

 
Microarray monitors the expression levels of hundreds of thousands of genes simultaneously in a whole genome level. The spotted arrays (DeRisi et al., 1996) and the high-density oligonucleotide chips (Lockhart et al., 1996) are two commonly used microarray technologies. There are many statistical problems associated with microarray data. The readers are referred to the books by Parmigiani et al. (2003) and Speed (2003) for comprehensive reviews. In this paper we focus on differential gene expression detection and sample classification using microarray data obtained from different groups or conditions, e.g. stage I, II, III, IV breast cancer and normal tissues.

Consider a G-class microarray data, where we have measured the expression levels of p genes for n_g samples from class g, g = 1,..., G. Denote the measured expression values as

and we assume necessary preprocessing has been applied (Yang et al., 2002; Bolstad et al., 2003). We introduce the following class indicators

which satisfy

(1)

In differential gene expression detection, the basic idea is to compare the expression levels across the different classes. We can do the comparison for gene i using the following linear regression model

(2)

where the response has been centered and the intercept is not included in the model. ß_ig is interpreted as the difference between the mean expression value for class g and the mean across all classes.

We can formally test for no difference across all classes by testing ß_i1 = · · · = ß_iG, which can be done using the F-statistics obtained from the least squares fitting (Kutner et al., 2004)

where , and (SSB,SSW) are the between/within-class sum of squares, which are also equal to the regression/error sum of squares (SSR/SSE) of the linear regression model (2). Therefore the F-statistics is proportional to the commonly used ratio of between/within-class sum of squares (Golub et al., 1999).

After we cast differential gene expression detection in the framework of linear regression models, it is natural to consider alternative estimation methods other than least squares. In the following we discuss the differential gene expression detection and sample classification with the penalized linear regression models.

	2 PENALIZED LINEAR REGRESSION

TOP ABSTRACT 1 INTRODUCTION 2 PENALIZED LINEAR REGRESSION 3 APPLICATION TO MICROARRAY... 4 COMPARISON WITH OTHER... 5 DISCUSSION REFERENCES

 
2.1 penalty: LASSO regression
Linear regression with the penalty, known as the LASSO regression, has been shown to have shrinkage and variable selection properties (Tibshirani, 1996; Efron et al., 2004). The penalized regression model corresponds to the following constrained optimization

(3)

And we can show that the parameter estimations are (see Appendix for details)

which is known as soft-thresholding, the subscript plus means positive part (z₊ = z if z > 0 and zero otherwise), and sgn(·) is the sign function: sgn(z) = 1 if z > 0, sgn(z) = –1 if z < 0 and sgn(z) = 0 when z = 0.

With the estimated parameters the predicted mean value for class g is

(4)

Hence for the linear model setup (2), the penalized linear regression (3) is shrinking the difference between the individual class means and the overall mean. As a result the expression value is shrunken towards the global mean. The intuition is that the overall mean is relatively more stable than the individual class means owing to the bigger sample size. So the shrinkage can help to produce more reliable estimates of the true means (Donoho and Johnstone, 1994). Notice that if we make _g the same across classes, then (4) implicitly shrinks more those class means with smaller sample size.

Similar to the least squares estimation we define the total/error sum of squares (SSTO/SSE) for the penalized regression model as

and the penalized F-statistics for testing equal means across groups as

We can show that for gene i (see Appendix for details)

where s_i is the pooled standard error estimation for gene i,

In the ‘nearest shrunken centroid’ method (also known as PAM; Tibshirani et al., 2002, 2003), the following ad hoc soft-shrunken t-statistics and centroid are used

where is the shrinkage parameter and makes m_gs_i equal to the estimated standard error of the mean difference, . Thus d_ig is a t-statistics for gene i comparing class g with the average class. We can equivalently write the previous equations as

(5)

Therefore the shrunken centroid is equivalent to the predicted mean response (4) from the penalized linear model if we take _g = 2n_gm_g s_i. Similar to (4), in (5) the class means with smaller sample size are shrunken more since .

2.2 Differential gene expression detection and sample classification
Most differential gene expression detection methods, being parametric, non-parametric or Bayesian, are based on thresholding. Usually a score is assigned to each gene. Commonly used scores include the P-values and some variants of t/F-statistics. A cutoff value t₀ is chosen and those genes with the score bigger than t₀ are declared as differentially expressed. The false discovery rate (FDR; Benjamini and Hochberg, 1995) has become increasingly adopted in large-scale genomic studies to assess the significance of the identified set of genes. FDR is defined as the expected proportion of false positives and used as a guidance for choosing the cutoff value t₀.

We can also select a set of ‘interesting’ genes using sample classification errors, which can be formulated as the following optimization problem

where is the classification error using genes with some classifiers, e.g. K-nearest-neighbor, support vector machine (SVM), random forest (RF), etc. (Dudoit et al., 2002; Wu et al., 2003). Several difficulties are associated with this approach. First, since it is a combinatorial optimization, it is very hard to get the global optimum. Second, the results will depend on the classifier used. For microarray data it is common to observe p >> n, so overfitting is very likely to occur. Just by chance we can find many combinations of genes to achieve small classification errors.

One way to overcome some of the difficulties is to use some heuristic optimization with regularization. In the PAM method, the shrinkage parameter controls the set of genes selected and the classification errors, i.e.

and the classification error is based on the cross validations using the parameter . Hence we transform the combinatorial optimization into an one-dimensional optimization problem indexed by parameter .

Microarray is often used as an exploratory tool. The goal of differential gene expression detection is usually to identify a set of ‘interesting’ genes for follow-up study. How many genes we select as ‘interesting’ is really a balance of many considerations. The mathematical formulations as minimizing classification error or FDR are really a simplistic approach. Owing to the small sample size, the estimated FDR or classification errors will have large variance. In practice they should be combined with other factors to select genes. For example, owing to cost considerations, only a limited number of genes can be further studied.

2.3 Penalty parameter _g selection
In the PAM method, the tuning parameter is chosen to shrink the majority of the to zero, and those genes with at least one non-zero are selected to build the classifier. Notice that when _g = 2n_gm_gs_i ,

so if we use zero as the cutoff, the penalized F-statistics will select the exact same set of genes as differentially expressed.

For two-class microarray data (G = 2), let ₁ = ₂ = and notice that

hence the penalized F-statistics is

which is equivalent to the two-class penalized F-statistics proposed in Wu (2005), where a different coding for y is used. Wu (2005) proposed two algorithms to select the penalty parameter , which are based on minimizing the dependence between s_i and . This can be applied to select the parameters _g for the penalized F-statistics. In addition, analogous to the sample classification approach, we propose to select both _g and differentially expressed genes simultaneously using the FDR: for specific values of _g, we will declare those genes with as differentially expressed, and the associated FDR can be estimated using permutations.

In the following we first illustrate the selection of penalty parameter using the small round blue cell tumors (SRBCT) microarray data reported at Khan et al. (2001), which has also been analyzed in Tibshirani et al. (2002). We then compare the proposed method with other alternative classification methods using some public microarray data.

	3 APPLICATION TO MICROARRAY DATA

TOP ABSTRACT 1 INTRODUCTION 2 PENALIZED LINEAR REGRESSION 3 APPLICATION TO MICROARRAY... 4 COMPARISON WITH OTHER... 5 DISCUSSION REFERENCES

 
The SRBCT microarray data measured the expression levels of 2308 genes for 63 training samples and 25 testing samples from four tumor types: Burkitt lymphoma (BL, 8 samples), Ewing sarcoma (EWS, 23 samples), neuroblastoma (NB, 12 samples) and rhabdomyosarcoma (RMS, 20 samples). The normalized expression values for these 2308 genes can be downloaded from http://research.nhgri.nih.gov/microarray/Supplement

3.1 Differential gene expression detection and FDR estimation/control
With some selected values of _g, we can shrink the majority of the ₁ penalized F-statistics to zero, and those genes with non-zero can be declared as significant. The number of false positives owing to random chance can be estimated by permutations. Specifically we can randomly permute the class indicator and recalculate the penalized F-statistics. The average number of non-zeros over permutations are used to estimate the number of false positives. Figure 1 shows the false positive estimations for these 2308 genes using _g = . We can see that the number of false positives is extremely small. One possible explanation is that in Khan et al. (2001), these 2308 genes are the quality filtered subset of the original 6567 genes and possibly lots of non-significant genes have been removed. The extremely small proportion of false positives observed in Figure 1 also partially explains the good classification accuracy for the SRBCT microarray data (Khan et al., 2001; Tibshirani et al., 2002).

View larger version (5K): [in this window] [in a new window]   Fig. 1 False positive estimation for 2308 genes in SRBCT microarray data: R is the number of significant genes for different and FP is the number of false positives averaged over 1000 permutations.

 
3.2 Sample classification
Similar to the PAM method, we can apply the nearest centroid classification using the shrunken centroid estimated from the penalty models (4) to select genes and classify samples simultaneously. Considering the small sample sizes for different classes of tumors (8, 23, 12, 20), we choose 5-fold CV to estimate the classification errors.

There are two levels in terms of the penalty parameter selections: the relation of _g across different classes within the same gene and the relation of penalty parameters across different genes. Notice that in the penalized linear regression model (3), the penalty parameter _g implicitly depends on the sample size n.

For within-gene penalty parameter selection, we can use common _g across classes: = _g/n; or use standardized penalty regression: each predictor is standardized first and then the same penalties are applied to all regression coefficients, which corresponds to using

where is the standard deviation for y_gj. Let ß_ig = _ig/_g and we have

For across gene penalty parameter selection, we can use the same for different genes. Or we can scale by the standard error of individual gene expressions, i.e. using penalty parameter s_i for gene i.

Table 1 summarizes the parameter estimations for the combination of standardization within-gene and scaling across-gene.

View this table: [in this window] [in a new window]   Table 1 Penalized parameter estimation

 
In 5-fold CV, we randomly partition the 63 samples into five similar groups. Each time one group is held out as a testing set. Notice that in the CV, n and n_g will be different for different combination of training groups. Figure 2 shows the classification error estimations from the 5-fold CV. We can see that the scaling helps to decrease the classification error and reduce the number of genes used. Using both standardization and scaling can substantially improve the classification.

View larger version (13K): [in this window] [in a new window]   Fig. 2 Average classification errors over 100 random 5-fold CVs: the y-axis is the average number of misclassified samples, the x-axis labels the number of genes with non-zero . The four plots represent the combination of standardization (first row: unstd versus second row: std) and scaling (first column: unscaled versus second column: scaled).

 
Next we compare the proposed method (using both standardization and scaling) with some popular classification methods by application to some public microarray data.

	4 COMPARISON WITH OTHER CLASSIFICATION METHODS

TOP ABSTRACT 1 INTRODUCTION 2 PENALIZED LINEAR REGRESSION 3 APPLICATION TO MICROARRAY... 4 COMPARISON WITH OTHER... 5 DISCUSSION REFERENCES

 
We evaluate classification performance using the following three public microarray data listed at Table 2.

View this table: [in this window] [in a new window]   Table 2 Three public microarray data

 
We randomly choose one-third of the samples as testing set and the other samples as training set. For each training set, K-fold CV is used to train the classifier, which is then used to predict the testing set. We use K = 5 for the breast cancer data and K = 3 for the other two data due to their relatively small sample sizes in each group. We repeat the process 50 times to get the average misclassification errors. The same procedures are applied to SVM (Burges, 1998), RF (Breiman, 2001), PAM (Tibshirani et al., 2002) and the proposed method (called RegCL in the following discussion). For SVM multiclass-classification, we employ the ‘one-against-one’-approach, and the appropriate class is found by a voting scheme (Chang and Lin, 2001).

Table 3 summarizes the average classification errors over 50 random CVs. We can see that the proposed method has favorable performance compared with other methods (for detailed results please see the Supplementary website).

View this table: [in this window] [in a new window]   Table 3 Average classification errors (%) for the microarray data

 

	5 DISCUSSION

TOP ABSTRACT 1 INTRODUCTION 2 PENALIZED LINEAR REGRESSION 3 APPLICATION TO MICROARRAY... 4 COMPARISON WITH OTHER... 5 DISCUSSION REFERENCES

 
Due to the ‘large p small n’ commonly observed for microarray data, regularization can help to prevent over-fitting and produce more reliable estimations. Some ad-hoc shrinkage methods have been proposed to utilize the shrinkage ideas and prove to be useful in empirical studies. In this paper we show that the ad hoc shrinkage methods can be rigorously derived from the penalized linear regression models, which provide a unified statistical framework for the penalized differential gene expression detection (Wu, 2005) and the ‘nearest shrunken centroid’ classification methods (Tibshirani et al., 2002). Through applications to public microarray data, we illustrate the importance of the regularization parameter selection. By comparing with some other popular classification methods, we showed the favorable performance of the proposed methods.

Recently Storey et al. (2005) proposed the optimal discovery procedure to fully utilize the information from all the genes to improve differentially gene expression detection. It would be interesting to study how to incorporate dependence information among genes into the proposed penalized F-statistics for differential gene expression detection.

	Acknowledgments

 
This research was partially supported by a startup fund from the Division of Biostatistics and a research grant from the Graduate School of University of Minnesota. I would like to thank two referees for their constructive comments, which have greatly improved the presentation of the paper.

Conflict of Interest: none declared.

	FOOTNOTES

 
Associate Editor: John Quackenbush

Received on October 11, 2005; revised on December 6, 2005; accepted on December 7, 2005

	REFERENCES

TOP ABSTRACT 1 INTRODUCTION 2 PENALIZED LINEAR REGRESSION 3 APPLICATION TO MICROARRAY... 4 COMPARISON WITH OTHER... 5 DISCUSSION REFERENCES

 

Alizadeh, A., et al. (2000) Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature, 403, 503–511[CrossRef][Medline].

Benjamini, Y. and Hochberg, Y. (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B, . 57, 289–300.

Bolstad, B., et al. (2003) A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics, 19, 185–193[Abstract/Free Full Text].

Breiman, L. (2001) Random forests. Mach. Learning, 45, 5–32[CrossRef].

Burges, C.J.C. (1998) A tutorial on support vector machines for pattern recognition. Data Min. Knowl. Disc, . 2, 121–167[CrossRef].

Chang, C.C. and Lin, C.J. (2001) Libsvm : a library for support vector machines. Technical Report, available at http://www.csie.ntu.ed.tw/~cjlin/libsvm.

DeRisi, J., et al. (1996) Use of a cDNA microarray to analyse gene expression patterns in human cancer. Nat. Genet, . 14, 457–460[CrossRef][ISI][Medline].

Donoho, D.L. and Johnstone, I.M. (1994) Ideal spatial adaptation by wavelet shrinkage. Biometrika, 81, 425–455[Abstract/Free Full Text].

Dudoit, S., et al. (2002) Statistical methods for identifying differentially expressed genes in replicated cDNA microarray experiments. Stat. Sinica, 12, 111–139.

Efron, B., et al. (2004) Least angle regression. Ann. Stat, . 32, 407–499[CrossRef].

Golub, T.R., et al. (1999) Molecular classification of cancer: class discovery and class prediction by gene expression monitoring. Science, 286, 531–537[Abstract/Free Full Text].

Khan, J., et al. (2001) Classification and diagnostic prediction of cancers using gene expression profiling and artificial neural networks. Nat. Med, . 7, 673–679[CrossRef][ISI][Medline].

Kutner, M., Nachtsheim, C., Wasserman, W. Applied Linear Regression Models, (2004) 4th edn. , Mcgraw-Hill/Irwin, NY.

Lockhart, D., et al. (1996) Expression monitoring by hybridization to high-density oligonucleotide arrays. Nat. Biotechnol, . 14, 1675–1680[CrossRef][ISI][Medline].

In Parmigiani, G., Garrett, E.S., Irizarry, R.A., Zeger, S.L. (Eds.). The Analysis of Gene Expression Data: Methods and Software, . (2003) , Springer, NY.

Pomeroy, S.L., et al. (2002) Prediction of central nervous system embryonal tumor outcome based on gene expression. Nature, 415, 436–442[CrossRef][Medline].

In Speed, T. (Ed.). Statistical Analysis of Gene Expression Microarray Data, (2003) Chapman & Hall/CRC, USA.

Storey, J.D., Dai, J.Y., Leek, J.T. (2005) The optimal discovery procedure II: applications to comparative microarray experiments. UW Biostatistics Working Paper Series, , pp. 260 available at http://www.bepress.com/uwbiostat/paper260.

Tibshirani, R. (1996) Regression shrinkage and selection via the lasso. J. R. Statistical Soc. B, Methodological, 58, 267–288.

Tibshirani, R., et al. (2002) Diagnosis of multiple cancer types by shrunken centroids of gene expression. Proc. Natl. Acad. Sci. USA, 99, 6567–6572[Abstract/Free Full Text].

Tibshirani, R., et al. (2003) Class prediction by nearest shrunken centroids, with application to DNA microarrays. Stat. Sci, . 18, 104–117[CrossRef][ISI].

Tusher, V.G., et al. (2001) Significance analysis of microarrays applied to the ionizing radiation response. Proc. Natl. Acad. Sci. USA, 98, 5116–5121[Abstract/Free Full Text].

West, M., et al. (2001) Predicting the clinical status of human breast cancer by using gene expression profiles. Proc. Natl. Acad. Sci. USA, 98, 11462–11467[Abstract/Free Full Text].

Wu, B. (2005) Differential gene expression detection using penalized linear regression models: the improved SAM statistics. Bioinformatics, 21, 1565–1571[Abstract/Free Full Text].

Wu, B., et al. (2003) Comparison of statistical methods for classification of ovarian cancer using mass spectrometry data. Bioinformatics, 19, 1636–1643[Abstract/Free Full Text].

Yang, Y.H., et al. (2002) Normalization for cDNA microarray data: a robust composite method addressing single and multiple slide systematic variation. Nucleic Acids Res, . 30, e15[Abstract/Free Full Text].

CiteULike    Connotea    Del.icio.us    What's this?

This article has been cited by other articles:

				Y. Lai Genome-wide co-expression based prediction of differential expressions Bioinformatics, March 1, 2008; 24(5): 666 - 673. [Abstract] [Full Text] [PDF]

				T. Abeel, Y. Saeys, E. Bonnet, P. Rouze, and Y. Van de Peer Generic eukaryotic core promoter prediction using structural features of DNA Genome Res., February 1, 2008; 18(2): 310 - 323. [Abstract] [Full Text] [PDF]

				F. Tai and W. Pan Incorporating prior knowledge of gene functional groups into regularized discriminant analysis of microarray data Bioinformatics, December 1, 2007; 23(23): 3170 - 3177. [Abstract] [Full Text] [PDF]

				S. Wang and J. Zhu Improved centroids estimation for the nearest shrunken centroid classifier Bioinformatics, April 15, 2007; 23(8): 972 - 979. [Abstract] [Full Text] [PDF]

This Article Abstract FREE Full Text (Print PDF) Supplementary Material A corrigendum has been published All Versions of this Article: 22/4/472    most recent bti827v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Search for citing articles in: ISI Web of Science (7) Request Permissions Google Scholar Articles by Wu, B. Search for Related Content PubMed PubMed Citation Articles by Wu, B. Social Bookmarking          What's this?

Online ISSN 1460-2059 - Print ISSN 1367-4803

Copyright © 2008 Oxford University Press

Oxford Journals Oxford University Press

• Site Map
• Privacy Policy
• Frequently Asked Questions

Other Oxford University Press sites: Oxford University Press Oxford Journals Japan American National Biography Booksellers' Information Service Children's Fiction and Poetry Children's Reference Corporate & Special Sales Dictionaries Dictionary of National Biography Digital Reference English Language Teaching Higher Education Textbooks Humanities International Education Unit Law Medicine Music Online Products Oxford English Dictionary Reference Rights and Permissions Science School Books Social Sciences Very Short Introductions World's Classics