An extended data mining method for identifying differentially expressed assayspecific signatures in functional genomic studies
 Derrick K Rollins^{1, 2}Email author and
 AiLing Teh^{2}
DOI: 10.1186/17560381311
© Rollins and Teh; licensee BioMed Central Ltd. 2010
Received: 30 June 2010
Accepted: 17 December 2010
Published: 17 December 2010
Abstract
Background
Microarray data sets provide relative expression levels for thousands of genes for a small number, in comparison, of different experimental conditions called assays. Data mining techniques are used to extract specific information of genes as they relate to the assays. The multivariate statistical technique of principal component analysis (PCA) has proven useful in providing effective data mining methods. This article extends the PCA approach of Rollins et al. to the development of ranking genes of microarray data sets that express most differently between two biologically different grouping of assays. This method is evaluated on real and simulated data and compared to a current approach on the basis of false discovery rate (FDR) and statistical power (SP) which is the ability to correctly identify important genes.
Results
This work developed and evaluated two new test statistics based on PCA and compared them to a popular method that is not PCA based. Both test statistics were found to be effective as evaluated in three case studies: (i) exposing E. coli cells to two different ethanol levels; (ii) application of myostatin to two groups of mice; and (iii) a simulated data study derived from the properties of (ii). The proposed method (PM) effectively identified critical genes in these studies based on comparison with the current method (CM). The simulation study supports higher identification accuracy for PM over CM for both proposed test statistics when the gene variance is constant and for one of the test statistics when the gene variance is nonconstant.
Conclusions
PM compares quite favorably to CM in terms of lower FDR and much higher SP. Thus, PM can be quite effective in producing accurate signatures from large microarray data sets for differential expression between assays groups identified in a preliminary step of the PCA procedure and is, therefore, recommended for use in these applications.
Introduction
It is well known that living organisms have complicated gene structures. However, while major advancements have been made in recent years, understanding of the biological functions of each individual gene is still quite limited. Active research is strongly focused on understanding the behavior of genes and as well as the highly complex metabolism and regulatory network inside living cells [1]. This effort falls under a molecular biological field called functional genomics (FG). There are at least three areas in which experimental techniques are widely applied in FG: transcriptomics, proteomics, and metabolomics [2]. A combination of leading scientific techniques as well as powerful mathematical and statistical tools for data analysis makes the task of identifying important transcriptome, proteome, and metabolome corresponding to a biological effect promising. Typical studies in these areas involve the identification of possible behavior and responses of species under various genetic backgrounds as well as environmental factors (i.e. assay).
There are different high technology techniques applied in FG field to advance understanding of the transcriptional genetic response of many organisms in various environmental perturbations [1]. One of the techniques that have been adopted in this field is a multiplex technology called DNA microarray [3]. A new technique that is becoming popular and will probably displace arraybased measurement in FG is nextgeneration sequencing (RNAseq) [4, 5]. These techniques have the ability to generate data sets that consist of expression levels of thousands of genes, providing a wealth of information that is hidden by high noise levels, low signal levels, and a relatively small number of experimental units to the number of genes studied. More specifically, since the data set containing the gene expression measurements consists of a lot more genes than assays, analytical techniques are needed to provide accurate gene identification under a large number of gene candidates that is much greater than the number of experimental runs.
To achieve this objective, traditional statistical methods, such as principal component analysis (PCA) [2–8], the focus of this article, are being retrofitted to provide effective statistical inference in this challenging context of microarray data analysis. Other methods used in this field included linear model analysis [9–14], Bayesian method [15–17] and network component analysis (NCA) [18–20]. Thus, statistics is playing a critical role through the development of methodologies that give high statistical power (SP) (i.e., accurate identification), and low false discovery rate (FDR)[21] (i.e. low misidentification). To this end, this article introduces two new PCA based statistics for determining gene rank for differential expression between two PCA identified assay groups. This work extends the technique introduced by Rollins et al.[2] that determines gene rank for a single PCA identified assay group. Thus, the proposed method (PM) in this work is aimed at finding the genes with high expression levels in one group and low expression levels in the other group.
The PM uses PCA to first establish the existence of the assay groupings of interest. Then using the results that established the grouping, the differential contribution for each gene is determined using a statistic based on eigenvalues. This article proposes and evaluates two statistics. The first one is the group averaged difference of eigenvalue linear combinations that we call T _{ diff }. The second one divides T _{ diff } by its estimated pooled standard deviation that we call T _{ scaled }. The genes are ranked based on the largest absolute value of these statistics. The PM is evaluated against the ranking determined by the well known Student's tstatistic [14] that we call T _{ pooled } in this work. We will refer to T _{ pooled } as the current method (CM) which is actually a subclass of the PM that weighs each assay equally in each group. Note that for the CM the assay members in each group is not established based on the data but by á priori considerations. In contrast, for the PM the data drives the assay weight as well as group assignment of the assays.
The CM and PM are applied in the following three case studies to compare their effectiveness (i.e., power) in identifying assayspecific signature: (i) exposure of E. coli cells to two different levels of ethanol concentration [22]; (ii) the use of myostatin as inhibitor of skeletal muscle growth for five 5weeksold myostatin and nontreated mice [10]; and (iii) a simulation study based on statistical properties of the second case study.
This work is organized into the following sections. The Background Section gives a brief review of PCA and connects it to our application in FG's data analysis. This section is followed by the Methods Section that derives and presents the test statistics of the CM and PM. These test statistics are evaluated and compared in three studies in the Results and Discussion Section. The final section summarizes the results and gives concluding remarks on the contribution of this work.
Background
The microarray data set is given as an m by n matrix X where n is the number of assays expressed along columns (i.e. variables) and m represents the number of genes expressed along rows. The cells in this matrix are given as x _{ ij } which is the expression level of the i ^{th} gene for the j ^{th} assay (i.e. condition). Principal component analysis (PCA) is a multivariate technique that mathematically transforms (rotates) the original coordinate system to a new orthogonal coordinate system based on correlations among the variables [23]. The principal components (PCs) are eigenvectors generated from either the covariance matrix (scaled sum of squares and cross products) or the correlation matrix (sums of squares and cross products from standardized data) of the variables involved. They are used to construct n linear combinations of the n variables that can be thought of as n pseudo variables [23]. A PC is rank ordered by the amount of variation in the original data set that it captures.
Methods
Eigengene Contribution Approach
The basic difference between the method in Rollins et al. [2] and this extension is that work developed gene signatures for individual groups using equations of the form given by (3) and (4) and this work uses equation of the form given by Eq. 5.
Eigenassay Contribution Approach
Test Statistics
From Eq. 23 it is clear that scaling the EA differential contribution is not providing any new technique in PCA and therefore is not a useful result under the PM. Thus, we do not propose scaling for the EA approach.
 1.
Standardize X to obtain Z.
 2.
Obtain the loading and scores matrices for X (EG) based on correlation.
 3.
Obtain the loading and scores matrices for X ^{ T } (EA) based on covariance.
 4.
For each of the n EG loading vectors, plot its loadings against the assay number. Select the plot(s) that separate the assays into desired or interesting groups for further analysis.
 5.
For each n EA score vectors, plot its scores against the assay number. Select the plot(s) that separate the assays into desired or interesting groups for further analysis.
 6.
For each selected EG loading vector in Step 4, using Z and Eq. 5 determine the differential EG contribution for each gene.
 7.
For each selected EA loading vector in Step 5, using X and Eq. 9 determine the differential EA contribution for each gene.
 8.
For each case in Steps 6 and 7, rank order the differential contribution and then table (with the corresponding gene) and plot these values against the rank. These signature plots can be used to determine where to make cutoffs as described in Rollins et al. [2].
In the next section we evaluate the proposed test statistics that we have derived in this section against a current method that uses the Student's ttest statistic. This work also includes an evaluation to determine when it is better to choose T _{ diff } or T _{ scale }.
Results and Discussion
The best choice for a test statistics is the one that has the highest statistical power (SP) and the lowest false discovery rate (FDR) [21]. This section presents three case studies to evaluate the proposed test statistics against one another and against a current method (CM) that uses T _{ pooled }. The first study revisits the single group analysis in Rollins et al. [2] involving exposure of E. coli cells to two different levels of ethanol concentration [2, 22]. The second study applies the proposed method (PM) to data from Steelman et al. [10]. This data set involves the use of myostatin as an inhibitor of skeletal muscle growth for five 5weeksold myostatin (called "mutant") and nontreated (called "wildtype") mice in each group. The third study is a mathematically simulated data study using characteristics of the data from Study 2.
Exposure of E. coli cells Study
The data set for the first case study contains E. coli cells that were exposed to two different ethanol concentrations. In Rollins et al.[2] ranked signatures were obtained for nonethanol (i.e., nontreated) (Group A) and ethanol (Group B) separately. Thus, these signatures ranked the genes based on their contribution to the score of their group. However, the goal of this work is to obtain a ranked signature of the genes that is based on the difference of gene contribution between the two groups. Therefore, under this objective, genes with high contribution in both groups would not be ranked high; whereas, genes with low contribution in one group and high contribution in the other group could be ranked high based on the greatest negative, positive, or absolute difference, depending on the interests of the experimenter. For this study, we ranked the genes based on absolute difference for evaluative purposes.
Top 20 genes that showed distinct difference between ethanol and nonethanol along with their ranking
Rank  Gene Name  EtOH Rank*  NonEtOH Rank*  Rank  Gene Name  EtOH Rank*  NonEtOH Rank* 

1  b2387  729  2001  11  argT  925  2330 
2  ybdO  558  2182  12  argH  2  4286 
3  b1455  959  2120  13  ycbE  317  3626 
4  gltD  2884  151  14  b0538  328  2658 
5  appY  360  2457  15  citB  372  2642 
6  caiA  5  3810  16  wbbH  2952  408 
7  b0960  2664  787  17  ccmD  2605  1083 
8  yaiD  1  4284  18  agaA  885  2483 
9  b1815  3178  43  19  ymcC  568  2587 
10  ydaK  375  2548  20  abc  389  2705 
Skeletal Muscle Growth in Mice Study
The second study is a data set that involved the use of myostatin as inhibitor of skeletal muscle growth for five 5weekold myostatin (called "mutant") and nontreated (called "wildtype") mice in each group. A powerful method for ranking genes and determining the size of signatures is the Qmethod developed by Storey and Tibshirani [12]. The Qmethod uses T _{ pooled } and a novel method for achieving high SP and low FDR. The Qmethod first uses T _{ pooled } to obtain pvalues then convert to qvalues to determine where to cutoff signatures based on a maximum qvalue. Given that the qvalue is related to the pvalue, one could also rank genes based on pvalues or their T _{ pooled } values which are inversely related. Since we are primarily interested in ranking genes in this work, we will compare the techniques based on the abilities of T _{ pooled } and the PM to find top ranked genes.
Top ranked genes of one method in the top 200 genes of the other method in study (ii)
x  x in top 200 CM genes  %x in Top 200 CM genes  y  y in top 200 PM genes  %y in Top 200 PM genes 

10  7  70  10  5  50 
20  9  45  20  7  35 
30  10  33  30  7  23 
40  10  25  40  11  28 
50  11  22  50  15  30 
60  17  28  60  17  28 
70  19  27  70  17  24 
80  20  25  80  19  24 
90  20  22  90  20  22 
100  21  21  100  21  21 
Simulation Study
Thus, 200 of the genes for each of the assays in Group A had the largest mean and were significantly different than all the other genes that had a mean of 5.3. The study will evaluate the ability of the CM and PM to identify these 200 genes when the variance for all the data in the data matrix is the same (Part 1) and when the variance differs from gene to gene (Part 2). Each result in the simulation study is an average of five trials. For simplicity, all the results in this study will be based on eigengene (EG) principal components (PCs) as it gave strong separation of the groups.
Simulation Study  Part 1
Simulation Study  Part 2
In the second simulation study we evaluated the techniques by varying levels of ${\sigma}_{{x}_{i}}^{2}$ for each gene and two levels of δ: 1 and 3. More specifically, the distribution for ${\sigma}_{{x}_{i}}$ was log normal with mean 0.37 and variance 0.37^{2}. Thus, for each data table a ${\sigma}_{{x}_{i}}$ was randomly generated for each gene i, i = 1, ..., m, and then ten simulated expression values, one for each assay, were generated according to Eqs. 24 and 25 for the given level of δ.
Our final analysis in this study evaluated performance in signature size determination. The CM is the Qmethod developed by Storey and Tibshirani [12] that uses the pvalues of the ttest (i.e., T _{ pooled }) and cuts the list off at a maximum Qvalue, commonly 0.05, the value used in this analysis. The PM is the Inflection Method (IM) that is described in Rollins et al.[2] that cuts the list off at the greatest change in the signature plot of the ranked genes. The results are from Part 1 of the simulation study with a constant ${\sigma}_{{x}_{i}}$ for all the genes in a data table.
Conclusion
This work proposed a new principal component analysis (PCA) method for analyzing large dimensional data set such as gene expressions data set. The strength of the proposed method (PM) comes from its data driven nature. It is data driven because the relationships obtained by PCA are only determined by those that exist in the data. Thus, no predetermined grouping or any á priori knowledge has influence on the principal components (PCs) obtained. After obtaining the PCs, they are used to match and verify the existence of the assay groups of interests. From the PCs that have the strongest match, the contribution of each gene providing the greatest differential expressions are identified and ranked. Thus, a PM signature is not just a difference of expression levels for genes but differences in a direction verified to have the characteristics of interests. This approach distinguishes PM from methods that do not form groups on the basis of data analysis and develop signatures from the differences between two groups in the original data space. One should be cautioned that as the number of members in the groups becomes smaller, the probability a particular order of the assays increases. Thus, for a small number of assays, one should require greater separation of groups for high confidence in the true existence of the groups.
Following Rollins et al. [2], the PM develops test statistics treating the assays as variables (eigengenes, EG) and the genes as variables (eigenassays, EA). These test statistics are linear combinations of these variables (i.e., pseudo variables) as determined from the elements of the eigenvectors. One test statistic, called T _{ diff } is the difference of the average expression levels between two groups of pseudo variables. The other test statistics, called T _{ scaled }, is T _{ diff } divided by the estimated pooled standard deviation. We compared the performance of these two test statistics with the common and popular Student's tstatistic, T _{ pooled } that we called the current method (CM). Two real data studies provided evidence in support of the PM as a viable technique. A simulation study provided the strongest supportive evidence for the use of T _{ diff } when the gene variability is fairly uniform throughout a data table and for T _{ scaled } when the variability is not fairly uniform. However, one should note that this study was done under a particular set of model assumptions. The most critical one is independence. If the data have a particular correlation structure, which is not uncommon given that all the genes in an assay experience the same set of conditions, the results in this article may not be supported. Future work will include evaluating the PM under the kinds of correlation structures found in real expression data. Finally, with the PM test statistics, the inflection method (IM) introduced by Rollins et al. [2], indicated strong promise in determining signature cutoffs in terms of statistical power and false discovery rate (FDR) as compared to CM.
We are applying the PM in a variety of applications involving biological as well as physical phenomenon, with promising results. These applications include: 1. Nitric Oxide and Snitrosoglutathioneresponsive genes in Ecoli; 2. analysis of DNA microarray data for juvenile small round blue cell tumors; 3. analysis of metabolite data from corn tissues (silk, pollen, coleoptile, and seedlings) for differential expression levels between the wild type and genetic mutations; 4. analysis of spectroscopy data for super alloys; and 5. the enhancement of nondestructive tests for ceramic armor in the resistance of ballistic penetration. Thus, the PM has potential application in a variety of situations where differential analysis is needed on large data sets with a relatively small number of different conditions or assays. It appears to have promise for these applications for high SP and low FDR as compared to other currently available methods.
Abbreviations
 CM:

current method
 EA:

Eigenassay
 EG:

Eigengene
 FDR:

false discovery rate
 FG:

functional genomics
 IM:

inflection method
 g _{ i } gene contribution for i ^{ th } :

gene
 l :

loading
 PCA:

Principal Component Analysis
 PC:

Principal Component
 PM:

proposed method
 S :

score
 SDG:

statistically different genes
 SP:

statistical power
 SS:

signature size
 T _{ diff } :

difference of the average expression levels between two groups of pseudo variables
 T_{pooled} :

Student's t statistics
 T _{ scaled } scaled statistics by dividing T _{ diff } :

by its estimated pooled standard deviation
 δ:

differential effect
Declarations
Acknowledgements
We would like to thank Professor Daniel S. Nettleton for his valuable advice and suggestions in conducting the simulation study. We are also grateful for his assistance in statistical software package used in this work. Funding for this work was provided through a grant from the Army Research Laboratory under Cooperative Agreement W911NF080036.
Authors’ Affiliations
References
 Dharmadi Y, Gonzalez R: DNA Microarrays: Experimental Issues, Data Analysis, and Application to Bacterial Systems. Biotechnol Progress. 2004, 5: 13091324. 10.1021/bp0400240.View ArticleGoogle Scholar
 Rollins DK, Zhai D, Joe AL, Guidarelli JW, Murarka A, Gonzalez R: A novel data mining method to identify assayspecific signatures in functional genomic studies. BMC Bioinformatics. 2006, 7: 37710.1186/147121057377.View ArticlePubMedPubMed CentralGoogle Scholar
 Zhang W, Carriquiry A, Nettleton D, Dekkers JCM: Pooling mRNA in Microarray Experiments and Its Effect on Power. Gene Expression. 2007, 23: 12171224.Google Scholar
 Shendure J, Ji H: Nextgeneration DNA sequencing. Nature Biotechnology. 2008, 26: 1010.1038/nbt1486.View ArticleGoogle Scholar
 Morozova O, Marra MA: Applications of nextgeneration sequencing in functional genomics. Genomics. 2008, 92: 255264. 10.1016/j.ygeno.2008.07.001.View ArticlePubMedGoogle Scholar
 Raychaudhuri S, Stuart JM, Altman RB: Principal Component Analysis to Summarize Microarray Experiments: Application to Sporulation Time Series. Pacific Symposium on Biocomputing. 2000, 5: 452463.Google Scholar
 Misra JW, Hwang D, Hsiao LL, Gullans S, Stephhanopoulos G: Interactive exploration of microarray gene expression patterns in a reduced dimensional space. Genome Res. 2002, 12: 11121120. 10.1101/gr.225302.View ArticlePubMedPubMed CentralGoogle Scholar
 Sharov AA, Dudekula DB, Ko MSH: A webbased tool for principal component and significance analysis of microarray data. Bioinformatics. 2005, 21: 25482549. 10.1093/bioinformatics/bti343.View ArticlePubMedGoogle Scholar
 Nettleton D: A discussion of Statistical Methods for Design and Analysis of Microarray Experiments for Plant Scientists. The Plant Cell. 2006, 18: 21122121. 10.1105/tpc.106.041616.View ArticlePubMedPubMed CentralGoogle Scholar
 Steelman CA, Recknor JC, Nettleton D, Reecy JM: Transcriptional profiling of myostatinknockout mice implicates Wnt signaling in postnatal skeletal muscle growth and hypertrophy. The FASEB Journal. 2006, 20: 580582.PubMedGoogle Scholar
 Dudoit S, Yang YH, Callow MJ, Speed TP: Statistical Methods for Identifying Differentially Expressed Genes in Replicated cDNA Microarray Experiments. Stat Sinica. 2002, 12: 111139.Google Scholar
 Storey JD, Tibshirani R: Statistical significance for genomewide studies. Proc Natl Acad Sci USA. 2003, 100: 94409445. 10.1073/pnas.1530509100.View ArticlePubMedPubMed CentralGoogle Scholar
 Ge Y, Sealfon SC, Speed TP: Statistical Methods in Medical Research. Statistical Method in MedicalResearch. 2009, 18: 54310.1177/0962280209351899.Google Scholar
 Devore JL: Probability and Statistics for Engineering and the Sciences. 2007, Duxbury Press, 7Google Scholar
 Sartor MA, Tomlinson CR, Wesselkamper SC, Sivaganesan S, Leikauf GD, Medvedovic M: Intensitybased Hierarchical Bayes Method Improves Testing for Differentially Expressed Genes in Microarray Experiments. BMC Bioinformatics. 2006, 7: 53810.1186/147121057538.View ArticlePubMedPubMed CentralGoogle Scholar
 Kendziorski CM, Newton MA, Lan H, Gould MN: On Parametric Empirical Bayes Methods for Comparing Multiple Groups Using Replicated Gene Expression Profiles. Statistics in Medicine. 2003, 22: 38993914. 10.1002/sim.1548.View ArticlePubMedGoogle Scholar
 Efron B, Tibshirani R, Storey JD, Tusher V: Empirical Bayes Analysis of a Microarray Experiment. Journal of American Statistical Association. 2001, 96: 45610.1198/016214501753382129.View ArticleGoogle Scholar
 Liao JC, Boscolo R, Yang YL, Sabatti C, Rpoychowdhury VP: Network component analysis: Reconstruction of regulatory signals in biological systems. The National Academy of Sciences. 2003, 100: 1552215527. 10.1073/pnas.2136632100.View ArticleGoogle Scholar
 Tran LM, Brynildsen MP, Kao KC, Suen JK, Liao JC: gNCA: A framework for determining transcription factor activity baed on transcriptome: identifiability and numerical implementation. Metabolic Engineering. 2005, 7: 128141. 10.1016/j.ymben.2004.12.001.View ArticlePubMedGoogle Scholar
 Hyduke DR, Jarboe LR, Tran LM, Chou KJY, Liao JC: Integrated network analysis identifies nitric oxide response networks and didihydratase as crucial target in. Proceedings of the National Academy of Sciences of the United States of America. 2007, 104: 84888489. 10.1073/pnas.0610888104.View ArticleGoogle Scholar
 Benjamini Y, Hochberg Y: Controlling the false discovery rate: apractical and powerful approach to multiple testing. Journal of Royal Statistical Society Series B. 1995, 57: 289300.Google Scholar
 Gonzalez R, Tao H, Purvis JE, Shanmugam KT, York SW, Ingram LO: Gene ArrayBased Identification of Changes That Contribute to Ethanol Tolerance in Ethanologenic Escherichia coli: Comparison of KO11 (Parent) to LY01 (Resistant Mutant). Biotechnol Prog. 2003, 19: 612623. 10.1021/bp025658q.View ArticlePubMedGoogle Scholar
 Johnson RA, Wichern DW: Applied Multivariate Statistical Analysis. 2008, Pearson Prentice Hall, 430470. 6Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.