Abstract
Background
Gene set analysis (GSA) methods test the association of sets of genes with a phenotype in gene expression microarray studies. Many GSA methods have been proposed, especially methods for use with a binary phenotype. Equally, if not more importantly however, is the ability to test the enrichment of a gene signature or pathway against the continuous phenotypes which are routinely and commonly observed in, for example, clinicopathological measurements. It is not always easy or meaningful to dichotomize continuous phenotypes into two classes, and attempting to do this may lead to the inaccurate classification of samples, which would affect the downstream enrichment analysis. In the present study, we have build on recent efforts to incorporate correlation structure within gene sets and pathways into the GSA test statistic. To address the issue of continuous phenotypes directly without the need for artificial discrete classification and thus increase the power of the test while ensuring computational efficiency and rigor, new GSA methods that can incorporate a covariance matrix estimator for a continuous phenotype may present an effective approach.
Results
We have designed a new method by extending the GSA approach called Linear Combination Test (LCT) from a binary to a continuous phenotype. Simulation studies and a real microarray dataset were used to compare the proposed LCT for a continuous phenotype, a modification of LCT (referred to as LCT_{2}), and two publicly available GSA methods for continuous phenotypes.
Conclusions
We found that the LCT methods performed better than the other two GSA methods; however, this finding should be understood in the context of our specific simulation studies and the real microarray dataset that were used to compare the methods. Free Rcodes to perform LCT for binary and continuous phenotypes are available at http://www.ualberta.ca/~yyasui/homepage.html webcite. The Rcode to perform LCT for a continuous phenotype is available as Additional file 1.
Background
Gene set enrichment analysis (GSEA) has greatly advanced highthroughput gene expression studies and a number of methods have been proposed to perform this type of analysis (see reviews by Goeman and Buhlmann [1] and Nam and Kim [2]). An important methodological challenge of GSA is the need to deal with large gene sets and small sample sizes. While most GSA methods employ a permutationbased approach to evaluate the significance gene sets, Kim and Volsky [3] gave a parametric view of the test statistic by assuming that the averages of fold changes across the genesets are distributed approximately normally. However, the majority of work in this field has focused on testing the enrichment of gene sets against binary, and sometimes categorical, phenotypes. Equally, if not more importantly, is the ability of the method to test the enrichment of a gene signature or a molecular pathway against a continuous phenotype. Such continuous variables are measured routinely and many important clinicopathological observations such as tumor size or the measurement of marker proteins are continuous. It may not always be technically easy or meaningful to categorize continuous phenotypes into two or more discrete classes. Indeed such artificial categorization may lead to inaccurate classification of the samples, which will eventually affect the downstream enrichment analysis.
We observed an important methodological distinction between the competitive and selfcontained GSA approaches [1,4]. Competitive methods use gene permutation to test whether or not the association of the phenotype with a gene set is similar to its association with the other gene sets (the “Q1 hypothesis”), while selfcontained methods employ sample permutation to test the equality of the two mean vectors of geneset expressions which correspond to the two phenotype groups (the “Q2 hypothesis”). Here, we focused on the selfcontained methods because, unlike the gene permutation approaches, sample permutation preserves correlations within gene sets; a property that we have used to design the proposed method for continuous phenotypes.
Correlations among gene expression measurements have long been observed, especially among measurements for functionally related gene set. Yet in the past, only the multivariate analysis of variance test for gene set analysis (MANOVAGSA) for categorical phenotypes [5] and the Linear Combination Test (LCT) for binary phenotypes [6] have used a covariance matrix estimator of gene expressions to compute the enrichment test statistic. The main challenges in using these methods are the relatively small sample sizes and large gene sets; a situation which is not uncommon in GSA, especially in small microarray studies. To overcome these difficulties, shrinkage methods [7] have been used to estimate the gene expressions covariance matrix. However, GSA has rarely been used for continuous phenotypes, and currently no methods that incorporate a covariance matrix estimator are available. Previously, when we compared the performances of various selfcontained GSA methods for binary phenotypes, we found that LCT was more computationally efficient than MANOVAGSA and approximated its superior power very well. Here, we propose both an extension of LCT to continuous phenotype (hereafter referred to as LCT) and a modified version of LCT (hereafter referred to as LCT_{2}). We compared the proposed methods with two existing GSA methods for continuous phenotype; namely, Significance Analysis of Microarrays for Gene Sets (SAMGS) [8] and Global Test [9]. We used simulations to compare the performances of the GSA methods with small sample sizes and large gene sets.
In addition, we analyzed the performances of the GSA methods using real microarray gene expression data from prostate tumor samples of AfricanAmerican prostate cancer patients [10]. Increased plasma or serum leptin levels have previously been found to be associated with development of prostate cancer [1113]. We, therefore, used the C2 catalog, an extensive collection of metabolic and signaling pathways and gene sets, from the Molecular Signatures Database (MSigDB) of Gene Set Enrichment Analysis application of Broad Institute of MIT and Harvard. The catalog was screened for associations with human leptin gene (LEP) expression, a wellstudied marker of adiposity, and various metabolic and inflammatory conditions, and we identified important molecular pathways that were associated with high expression of this marker in a prostate cancer cohort. In our comparative study, we focused on testing both the power and computational efficiency of the four GSA methods.
Methods
Linear combination test for a continuous phenotype
In this section we derive the LCT for a continuous phenotype. Our derivations follow the binary phenotype framework in that the correlation structure is accommodated in a similar way to the binary phenotype, and the shrinkage covariance matrix estimation is implemented to take care of the small sample size and large gene set problems.
Consider gene expression data consisting of a total of n subjects. The null hypothesis to be tested is, that the expression of a predefined gene set with p genes, {X_{1}, …, X_{p}} is not associated with the phenotype Y. One way of expressing this multivariate hypothesis univariately as a null hypothesis is H_{0}; that is, no linear combination of X_{1},…, X_{p} is associated with the phenotype. Let Z(β) = β_{1}X_{1} + … + β_{p}X_{p} be a linear combination of X_{1},…, X_{p}. Then, for a given vector β of combination coefficients, whether or not the combination Z(β) is associated with the phenotype can be tested in the univariate model as follows: Y_{i} = α_{0} + α_{1}Z_{i}(β) + e_{i}, where i denotes subjects 1, …, n , Y_{i} denotes phenotype measurement of i^{th} subject, α_{0} and α_{1} are the intercept and slope respectively, and, e_{i} is ~ N(0,σ^{2}). This expression describes a classical simple linear regression problem. To test H_{0}, we can consider the mostsignificant linear combination of {X_{1},…, X_{p}}; namely, the linear combination with the maximum sample correlation with the phenotype, among all possible linear combinations. We have
as the coefficients of the mostsignificant linear combination. As the square of the sample correlation, ignoring , we have:
Where is the gene expressions covariance matrix with the hh'– Th entry being
where x_{hl} is the gene expression corresponding to gene h, and subject l. Therefore,
Where Cov_{Y.X}=(Y,X_{1}),…,Cov(Y,X_{p}))^{T} The above optimization problem can be written as
Where A=Cov_{Y,X}Cov_{Y,X}^{T} and . Thus, the solution to this optimization problem is the maximal eigenvector of AB^{1} and is the corresponding eigenvalue [14].
When the size of the gene set is larger than the sample size (i.e., p > n) the matrix B is singular. Similar to the adjustment implemented in MANOVAGSA [5], a possible remedy for the singularity is to employ a shrinkage covariance matrix as proposed previously by Schafer and Strimmer [7]. Thus, the singular covariance matrix can be replaced with a shrinkage covariance matrix given by with shrinkage coefficients if h=h'and , if h≠h' where ρ_{hh '} is the sample correlation between the h– th and h'– th genes, and the optimal shrinkage intensity can be estimated by .
The computational cost of incorporating the covariance matrix estimator into the test statistic in this way is very high. To address the computational efficiency problem, we use a group of normalized orthogonal bases, instead of the original observation vectors. First, we perform an eigenvalue decomposition of the shrinkage covariance matrix (V_{1},..., V_{p}) = (X_{1},..., X_{p})UD^{− 12} The square of the sample correlation can be rewritten as
where γ=D^{1/2}U^{T}β, Cov_{Y.V}=(Y,V_{1}),…,Cov(Y,V_{p}))^{T} According to a matrix algebra calculation [14], the coefficients of the mostsignificant combination are given by γ ∗ ∝ Cov_{Y,V.} This LCT statistic is, therefore, proportional to the sum over the gene set of the square covariance between the phenotype and gene expression measurements, after the orthogonal transformation
where c is a constant. The statistical significance can be evaluated against the null hypothesis with a permutation test (permuting phenotype labels) using this test statistic. The constant c can be ignored in the permutation test. This approach is advantageous computationally because is evaluated only once for the original data, and then there is no need to evaluate it for each permuted version of the data.
A modification of the linear combination test for a continuous phenotype
We also considered an alternative form of LCT (LCT_{2}) which we derived in the linear regression context. A least squares estimate of the regression function is given by where X represents the gene expression matrix and Y represents the vector of phenotype values. In case of singularities, a shrinkage version of the regression function estimate analogous to LCT can be obtained. An alternative version of LCT, can be derived as the square of the L_{2} norm of the shrinkage regression function .
Simulation study design
We carried out a number of simulation studies to compare the performances of the proposed LTC methods with two published selfcontained GSA methods for continuous phenotype; namely, an extension of SAMGS to continuous phenotype via regression analysis [8], and Global Test [9] which uses a random effects model to test the association between gene expression and phenotype. For each gene set of size p, we generated a gene expression matrix X_{nxp} We changed the number of observations n from 10 to 20, 50 and 100, and the gene set size p from 20 to 100, 200 and 400. We focused on scenarios where the gene set size is larger than the sample size, i.e. p > n, because these scenarios are more predominant and are challenging for GSA. We adopted a mixed correlation structure between genes in each set as follows: among the first p_{1} genes, the correlations are constant (ρ_{ij} = ρ); among the next p_{2} genes, the correlation between the ith and jth genes is ρ_{ij} = ρ^{ij} with ρ = 0.0, 0.3, 0.6 and 0.9; and the remaining genes are not correlated. The various simulation scenarios are summarized in Table 1. For each gene set, a continuous phenotype was generated from a normal distribution N(Xμ,I).where μ is a vector of length p. In the null model that we used to compare the size of the tests, we set μ to 0. In the alternative model that we used to check the power of the tests, first, we generated randomly five of the first 20 components of μ from N (ν,ν) with ν ranging from 0 to 2 in increments of 0.1, then, we generated randomly five of the next 20 components of μ from N (−ν,ν) with ν ranging from 0 to 2 in increments of 0.1; the rest remaining components were set at 0. The simulation data were replicated 1,000 times in each model. The pvalues are based on 1,000 permutations.
Table 1. Type I errors for four GSA methods: LCT, LCT_{2}, SAMGS and GlobalTest
The Rpackage to implement Global Test is available for download from http://www.bioconductor.org webcite. The LCT tests and SAMGS for continuous phenotypes were implemented by us using the R statistical software [15].
Results
Simulation study
We found that the type I errors were similar across the four GSA methods (Table 1). As the sample size increased, the type I error moved closer to the nominal level, as is expected when permutation of phenotype labels is used. The empirical power (with n = 20 and p = 100) was calculated using a nominal level of 0.05 for values of the ν parameter ranging from 0 to 5 in increments of 0.25, and correlations between each pair of genes of ρ = 0.0, 0.3, 0.6 and 0.9 (Figure 1). When there was no correlation among genes (ρ = 0.0), the four GSA methods exhibited very similar testing powers. At low correlation values, the LCT_{2} method appeared to be conservative and less powerful; perhaps, because LCT_{2} is based on shrinkage of the regression function, similar to the ridge regression method [16]. However, with increasing correlations among genes (ρ = 0.3, 0.6, 0.9), the differences in power values between the LCT and Global Test methods became increasingly larger. Compared with either the SAMGS or GlobalTest methods, LCT and LCT_{2} both exhibited much better ability to deal with the given correlations among genes.
Figure 1. Power comparison (n = 20 and p = 100) between four GSA approaches: LCT, LCT_{2}, SAMGS and Global Test.
Identifying gene sets associated with human leptin gene expression measurements
Leptin is a wellknown marker protein for human adiposity and the circulating levels of leptin in the blood are directly proportional to the total amount of body fat. Leptin is also associated with various metabolic and inflammatory conditions. We applied all four GSA methods to analyze a real Affymetrix microarray dataset consisting of genomewide transcriptomic measurements of prostate tumor samples from AfricanAmerican prostate cancer patients [10] against the continuous phenotype of the human leptin gene (LEP) expression values. The publicly available data were downloaded from Gene Expression Omnibus [17] [GEO:GSE6956]. The data that we used in the present study are part of a larger microarray study into immunobiological differences in prostate cancer tumors between AfricanAmerican and EuropeanAmerican men. Because LEP expression levels may be different between the two groups, we used only the data from the AfricanAmerican group to test the LCT methods. For our analysis, we used the expression values of 13,233 genes measured in tumor samples from 33 patients. The tumor samples were resected adenocarcinomas from patients who had not received any therapy before prostatectomy and were obtained from the National Cancer Institute Cooperative Prostate Cancer Tissue Resource (CPCTR) and the Department of Pathology at the University of Maryland. According to Wallace et al. [10], the macro dissected CPCTR tumor specimens were reviewed by a CPCTRassociated pathologist who confirmed the presence of tumors in the specimens. The tissues were collected between 2002 and 2004 at four different sites. The median age of prostatectomy was 61 and the median prostatespecific antigen (PSA) at diagnosis was 6.1 ng/mL. Fiftyfive percent of the tumors were stage pT2, and 45% were stage pT3 or more. Detailed RNA extraction, labeling and hybridization protocols were as described previously [10]. The gene expression values were centered and scaled across the samples before the four GSA methods were applied. The need for such standardization was pointed out in an earlier comparative study of GSA methods for a binary phenotype [18].
For comprehensive analysis, we used the C2 catalog from MsigDB [19] consisting of 1,892 gene sets, including metabolic and signaling pathways from major pathway databases, gene signatures from biomedical literature including 340 PubMed articles, as well as other gene sets compiled from published mammalian microarray studies. Following Subramanian et. al. 2005 [19], we restricted the size of gene sets to between 15 and 500 which gave us 1,403 gene sets for use in our analysis. Each gene set was tested for its association with the LEP expression measurements. A limitation of our study is that the findings come from a relatively small observational study and therefore cannot be generalized to other populations.
In terms of computational efficiency, we noted that LCT incorporated the covariance matrix into the estimations for only a small cost (CPU time of 413 seconds) compared with the cost using SAMGS (CPU time of 397 seconds). In contrast, Global Test was very computationally attractive (CPU time of 92 seconds). The CPU times were recorded on our PC (Processor: x86 Family 6 Model 23 Stepping 10 Genuine Intel 3Ghz; 4GB RAM).
We compared the pvalues for the gene sets obtained by the four methods; in particular, the lower pvalues, which we assumed indicted the most interesting gene sets. Table 2 shows the percentages of the gene sets for which the pvalues were less than 0.005, 0.01, 0.05, and 0.10 from the four GSA methods. We found that the performance of LCT and LCT_{2} was similar. The performance of SAMGS and Global Test was also similar but different from the performance of LCT and LCT_{2}, which is consistent with the results of the simulation study. To adjust for multiple comparisons when multiple gene sets are tested, false discovery rate (FDR) could be used instead of Type I error probability; however, the use of adjustment methods would not affect the conclusions of our comparative evaluation study. The FDR values were computed as described by Storey [20].
Table 2. Percentages of gene sets with pvalues less than 0.005, 0.01, 0.05 and 0.10
Gene sets and pathways that were identified, by at least one of the four GSA methods, to be associated with the LEP gene expression measurements (pvalue ≤ 0.05) are listed in Table 3 in ascending order of the pvalues obtained using LCT. The corresponding FDR values were 0.195 for LCT, 0.197 for LCT2, 0.135 for SAMGS, and 0.936 for Global Test. The adipocytokine signaling pathway was predicted to be strongly associated with LEP expression by all four GSA methods. This result was expected, given that adipocytokines are a group of adipose tissuederived hormones that includes leptin. In addition to being linked to obesity and diabetes, adipocytokines may be involved in the regulation of angiogenesis and tumor growth [21]. Regulation of autophagy was found to be associated with LEP expression consistent with previous findings that leptin played a role in the neuroendocrine control of autophagy [22]. Autophagy is a fundamental process in tumorigenesis and treatment response because it can act as a tumorsuppression mechanism, yet it can also enable tumor cell survival under conditions of metabolic stress, including nutrient deficiency [23]. Furthermore, LEP expression was strongly associated with both hypoxiainducible factor1 (HIF1) targets (LCT pvalue = 0.006; LCT_{2}, SAMGS and Global Test pvalue <0.03) and the hypoxia pathway (LCT pvalue = 0.035). Leptin can be activated in response to hypoxia in breast cancer cells where the process is mediated through hypoxiainducible factor1 [24,25].
Table 3. Gene sets and pathways associated with LEP gene expression measurements
Among the gene sets and pathways associated with LEP expression only by the LCT method, we highlight the insulin signaling candidate pathway (LCT pvalue = 0.049). A positive association between serum insulin levels and LEP expression has been reported in obese humans [26]. Furthermore, the association of circulating insulinlike growth factors with increased risk of prostate cancer has been reported in a metaanalysis [27]. Interestingly, the proteasome degradation candidate pathway was found to be significant by both Global Test (pvalue = 0.029) and SAMGS (pvalue = 0.028), but not by LCT (pvalue = 0.12). A small microarray study (N = 10) found that the genes in the proteasome degradation pathway were differentially expressed 72 hours after polyethylene glycolleptin injection [28]. Other gene sets and pathways found to be significantly associated with LEP expression but with less well elucidated roles are shown in Table 3 and may be worthy of future investigation.
Discussion
Many selfcontained GSA methods have been proposed. However, although many of these methods have the potential to be generalized to any design, they have only been illustrated for a binary or categorical outcome. Thorough extension of these methods to a continuous phenotype has rarely been reported, and studies into their implementation, simulation studies to check type I error and power, and their application to real datasets are lacking. Here, we describe the extension of a “selfcontained” GSA method from a binary to a continuous phenotype. The new GSA tests, LCT and LCT_{2}, address several important technical issues. First, they provide a rigorous and computationally efficient approach to extend the enrichment test of a given gene set against a continuous phenotype. This will be of great help in studying a variety of informative measurements that cannot always be easily or meaningfully reduced to binary or categorical phenotypes. Second, because a pathway often consists of genes that are together involved in a biological mechanism or disease, gene expression levels within a pathway are expected to be correlated. Yet most traditional GSA methods fail to accommodate this important characteristic feature of gene expression datasets. While permutation methods using a valid test statistic can result in appropriate Type I error, the incorporation of a covariance matrix estimator into the test statistic is highly desirable because it often results in better power. Furthermore, we noted that when the gene set to be tested is larger than the sample size, the covariance matrix is illconditioned. To address this problem, a shrinkage method for covariance matrix estimation can provide a useful solution; however, shrinkage methods are rarely used in GSA, in spite of their implementation as an Rpackage which is free for download [7]. The computational cost of including a shrinkage covariance matrix estimator, especially for permutationbased hypothesis testing, can be very high. Notably in our LCT algorithm, we overcame this difficulty by using an orthogonal transformation of the gene expression matrix. In the LCT algorithm, therefore, the eigenvalue decomposition of the shrinkage covariance matrix is performed only for the original data, and not for the permuted versions.
We focused here on selfcontained approaches and because competitive and selfcontained methods test different hypotheses, it is important for the user to make an informed choice based on the hypothesis of interest and their understanding of the limitations of the two approaches (see reviews by Nam and Kim [2] and Dinu et. al.[4]). An important limitation of the selfcontained approaches is that only a few genes can drive the association between the gene set and the phenotype. In such cases, posthoc analysis can be used to reduce the gene set to a core subset associated with the phenotype. Such an analysis has been reported in simulations and in a real example for a binary phenotype [4].
The improvements that we have incorporated into our new GSA tests have given these tests a variety of advantages over the existing methods. We hope that they will be used for the rigorous testing of associations between different molecular pathways and gene signatures. At least of the measured clinicpathological phenotypes are continuous. They include tissue features such as tumor size, staining based readouts; cellular characteristics such as the amount of lymphocytic infiltration in a tumor environment; and subjectspecific measurements such as diagnostic or prognostic marker protein or metabolite concentrations. The LCT algorithm can adjust for continuous or categorical covariates following a regression framework. The LEP levels in the prostate tumor example that we considered may also have been influenced by patientspecific covariates such as body mass index (BMI), age, and/or smoking status. Smoking status did not show a significant association with LEP expression levels (pvalue = 0.36), and BMI and age data were not available for our analysis.
To check the linearity assumption, exploratory data analysis should be used prior to running a formal inference. However, we noted that the small sample sizes that are common in microarray studies, would limit a thorough check for nonlinearities. We also noted that the LCT method could be extended to accommodate nonlinearities; however, a larger sample size would be needed. The simulations and real microarray studies which we conducted indicated that the LCT and LCT_{2} methods both performed very well for small sample sizes. The question of how small is small is debatable and depends largely on the study design. In the case of a binary/categorical phenotype, at least five samples per group are desirable. In the case of a continuous phenotype, assessing significance based on less than 10 samples is dangerous; an alternative would be to rely upon representations that are more descriptive/exploratory in nature. While LCT tests only linear associations between sets of genes and a continuous phenotype, SAMGS and Global Test have been extended in a generalized linear model (GLM) framework and can accommodate multiclass, continuous, count, rate, and censored survival phenotypes. SAMGS uses the sum of squares of the Wald statistic for individual genes constituting the pathway as the test statistic. Wald statistics are calculated as the ratio between the regression coefficient for an individual gene and its corresponding standard error. Global Test reduces the GLM to a random effects model, assuming the regression coefficients corresponding to the genes constituting the set are sampled from a common distribution with mean zero and constant variance. A score test statistic is used to test the null hypothesis of no association between the set and the phenotype. The SAMGS and Global Test algorithms can both adjust for covariates, an attractive feature when accounting for other known prognostic factors in the screening of associations between gene sets and a phenotype.
Conclusions
Our proposed LCT method for gene set analysis efficiently incorporates the gene expression covariance matrix into the test statistic. This approach has resulted in a powerful and computationally attractive method for testing the association of a given gene set with a continuous phenotype. Additional file 1.
Additional file 1. R code for the linear combination test (LCT) method for gene set analysis of a continuous phenotype.
Format: DOCX Size: 18KB Download file
Availability and requirements
Project name: Linear Combination Test for GeneSet Analysis of a Continuous Phenotype
Project home page: http://www.ualberta.ca/~yyasui/homepage.html webcite
Operating system: Microsoft Windows XP.
Programming language: R 2.10.1.
Abbreviations
GSA: Gene set analysis; LCT: Linear combination test; SAMGS: Significance analysis of microarray for gene sets.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
ID and XW developed the LCT methodology and designed/conducted the methodological study. SP, LEK and SV provided biological interpretations of the results of the analysis of the realworld dataset. ID, XW and SP drafted the manuscript which was critically reviewed and revised by all authors. All authors read and approved the final manuscript.
Acknowledgements
LEK is supported by Canadian Institutes of Health Research New Investigator award (MSH87734). SP is supported by Ramalingswami Fellowship from DBT, India and grants from MoS&PI and DST, India. We thank the reviewers for their constructive comments which have improved the exposition of this paper substantially.
References

Goeman JJ, Buhlmann P: Analyzing gene expression data in terms of gene sets: methodological issues.
Bioinformatics 2007, 23:980987. PubMed Abstract  Publisher Full Text

Nam D, Kim SY: Geneset approach for expression pattern analysis.
Brief Bioinform 2008, 9:189197. PubMed Abstract  Publisher Full Text

Kim SY, Volsky DJ, PAGE: Parametric analysis of geneSet enrichment.
Bioinformatics 2005, 6:144. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Dinu I, Potter JD, Mueller T, Liu Q, Adewale AJ, Jhangri GS, Einecke G, Famulsky KS, Halloran PF, Yasui Y: Gene Set analysis and reduction.
Brief Bioinform 2008, 10(1):2434. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Tsai , Chen JJ: Multivariate analysis of variance test for gene set analysis.
Bioinformatics 2009, 25(7):897903. PubMed Abstract  Publisher Full Text

Wang X, Dinu I, Liu W, Yasui Y: Linear combination test for hierarchical gene Set analysis.
Stat Appl Genet Mol Biol 2011., 10(1)
Article 13

Schäfer J, Strimmer K: A shrinkage approach to largescale covariance matrix estimation and implications for functional genomics.

Adeniyi AJ, Dinu I, Potter JD, Liu Q, Yasui Y: Pathway analysis of microarray data via regression.
J Comput Biol 2008, 15(3):269277. PubMed Abstract  Publisher Full Text

Goeman JJ, van de Geer SA, de Kort F, van Houwelingen HC: A global test for groups of genes: testing association with a clinical outcome.
Bioinformatics 2004, 20:9399. PubMed Abstract  Publisher Full Text

Wallace TA, Prueitt RL, Yi MH, Howe TM, Gillespie JW, Yfantis HG, Stephens RM, Caporaso NE, Loffredo CA, Ambs S: Tumor immunobiological differences in prostate cancer between AfricanAmerican and EuropeanAmerican Men.
Cancer Res 2008, 68:927936. PubMed Abstract  Publisher Full Text

Chang S, Hursting SD, Contois JH, Strom SS, Yamamura Y, Babaian RJ, Troncoso P, Scardino PT, Wheeler TM, Amos CI, Spitz MR: Leptin and prostate cancer.
Prostate 2001, 46:6267. PubMed Abstract  Publisher Full Text

Saglam K, Aydur E, Yilmaz M, Goktas S: Leptin influences cellular differentiation and progression in prostate cancer.
J Urol 2003, 169:13081311. PubMed Abstract  Publisher Full Text

Singh SK, Grifson JJ, Mavuduru RS, Agarwal MM, Mandal AK, Jha V: Serum leptin: A marker of prostate cancer irrespective of obesity.
Cancer Biomark 2010, 7(1):1115. PubMed Abstract  Publisher Full Text

Applied Multivariate Statistical Analysis, Prentice Hall. 2002, 8081.

R Development Core Team: R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2009.
ISBN 3900051070, URL http://www.Rproject.org webcite

Hoerl AE, Kennard R: Ridge regression: Biased estimation for nonorthogonal problems.
Technometrics 1970, 12:5567. Publisher Full Text

Edgar R, Domrachev M, Lash AE: Gene expression omnibus: NCBI gene expression and hybridization array data repository.
Nucleic Acids Res 2002, 30:207210. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Liu Q, Dinu I, Adewale AJ, Potter JD, Yasui Y: Comparative evaluation of geneset analysis methods.
BMC Bioinforma 2007, 8:431. BioMed Central Full Text

Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP: Gene set enrichment analysis: a knowledgebased approach for interpreting genomewide expression profiles.
Proc Natl Acad Sci USA 2005, 102:1554515550. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Storey JD: A direct approach to false discovery rates.
Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2002, 64:479498. Publisher Full Text

Housa D, Housova J, Vernerova Z, Haluzik M: Adipocytokines and cancer.
Physiol Res 2006, 55:233244. PubMed Abstract  Publisher Full Text

Malik SA, Mariño G, BenYounès A, Shen S, Harper F, Maiuri MC, Kroemer G: Neuroendocrine regulation of autophagy by leptin.
Cell Cycle 2011, 10(17):29172923. PubMed Abstract  Publisher Full Text

White E, DiPaola RS: The doubleedged sword of autophagy modulation in cancer.
Clin Cancer Res 2009, 15:53085316. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Cascio S, Bartella V, Auriemma A, Johannes GJ, Russo A, Giordano A, Surmacz E: Mechanism of leptin expression in breast cancer cells: role of hypoxiainducible factor1α.
Oncogene 2008, 27:540547. PubMed Abstract  Publisher Full Text

Terrasi M, Fiorio E, Mercanti A, Koda M, Moncada CA, Sulkowski S, Merali S, Russo A, Surmacz E: Functional analysis of the 2548G/A leptin gene polymorphism in breast cancer cells.
Int J Cancer 2009, 125:10381044. PubMed Abstract  Publisher Full Text

Russell CD, Petersen RN, Rao SP, Ricci MR, Prasad A, Zhang Y, et al.: Leptin expression in adipose tissue from obese humans: depotspecific regulation by insulin and dexamethasone.
Am J Physiol 1998, 275:E507E515. PubMed Abstract  Publisher Full Text

Renehan AG, Zwahlen M, Minder C, O‘Dwyer ST, Shalet SM, Egger M: Insulinlike growth factor (IGF)I, IGF binding protein3, and cancer risk: systematic review and metaregression analysis.
Lancet 2004, 363(9418):13461353. PubMed Abstract  Publisher Full Text

Taleb S, Haafen RV, Henegar C, Hukshorn C, et al.: Microarray profiling of human white adipose tissue after exogenous leptin injection.
Eur J Clin Invest 2006, 36:153163. PubMed Abstract  Publisher Full Text