A novel data mining method to identify assay-specific signatures in functional genomic studies
BMC Bioinformatics volume 7, Article number: 377 (2006)
The highly dimensional data produced by functional genomic (FG) studies makes it difficult to visualize relationships between gene products and experimental conditions (i.e., assays). Although dimensionality reduction methods such as principal component analysis (PCA) have been very useful, their application to identify assay-specific signatures has been limited by the lack of appropriate methodologies. This article proposes a new and powerful PCA-based method for the identification of assay-specific gene signatures in FG studies.
The proposed method (PM) is unique for several reasons. First, it is the only one, to our knowledge, that uses gene contribution, a product of the loading and expression level, to obtain assay signatures. The PM develops and exploits two types of assay-specific contribution plots, which are new to the application of PCA in the FG area. The first type plots the assay-specific gene contribution against the given order of the genes and reveals variations in distribution between assay-specific gene signatures as well as outliers within assay groups indicating the degree of importance of the most dominant genes. The second type plots the contribution of each gene in ascending or descending order against a constantly increasing index. This type of plots reveals assay-specific gene signatures defined by the inflection points in the curve. In addition, sharp regions within the signature define the genes that contribute the most to the signature. We proposed and used the curvature as an appropriate metric to characterize these sharp regions, thus identifying the subset of genes contributing the most to the signature. Finally, the PM uses the full dataset to determine the final gene signature, thus eliminating the chance of gene exclusion by poor screening in earlier steps. The strengths of the PM are demonstrated using a simulation study, and two studies of real DNA microarray data – a study of classification of human tissue samples and a study of E. coli cultures with different medium formulations.
We have developed a PCA-based method that effectively identifies assay-specific signatures in ranked groups of genes from the full data set in a more efficient and simplistic procedure than current approaches. Although this work demonstrates the ability of the PM to identify assay-specific signatures in DNA microarray experiments, this approach could be useful in areas such as proteomics and metabolomics.
The availability of gene structure data for many organisms  has paved the way for the challenging task of assigning biological functions to each individual gene, and more challenging still, explaining the highly complex metabolic and regulatory networks inside the living cell, where genes, proteins, and metabolites all interrelate. The wealth of information and technologies created by the availability of these genome sequences ushered in what frequently is called the post-genomic era, along with the appearance of a new field called functional genomics (FG). FG refers to "the development and application of system-wide experimental approaches to assess gene function by making use of the information and reagents provided by structural genomics" . There are at least three areas of FG for which experimental techniques are currently well developed: transcriptomics, proteomics, and metabolomics. Using a combination of these techniques with mathematical and computational tools for data analysis, the cell transcriptome, proteome, and metabolome can be identified (which refers to the inventory of all transcripts, proteins, and metabolites, respectively). Typical studies in these areas consist of surveying the levels of these species under a variety of environmental conditions and/or genetic backgrounds (referred to as assays in this article). The ultimate goal of these studies is therefore the identification of assay-specific signatures in the surveyed domain; e.g. which transcripts, proteins, or metabolites are associated with a given physiological condition such as health or disease or the response of an organism to an environmental challenge.
By far the most valuable experimental technique in FG is DNA microarrays, which have been used to study the transcriptional response of many organisms to genetic and environmental perturbations . As with any other experimental tool in FG, the use of DNA microarrays to study gene expression results in large data sets, which can consist of measurements for thousands of genes. In contrast, the number of assays is typically less than hundred. Therefore, since the number of genes is much greater than the number of assays, efficient information extraction and dimensionality reduction methods are needed to obtain assay-specific gene signatures of ranked order, which is the objective of the method proposed in this work. Several statistical methodologies have been proposed to achieve this goal including principal component analysis (PCA) – a specific case of singular value decomposition (SVD), multidimensional scaling, cluster analysis, self-organizing maps, and Fisher discriminant analysis ( and references therein). Among them, PCA has been widely used not only in the areas of transcriptomics [4–6], but proteomics [7–9] and metabolomics [10–12] as well, and it has shown much promise in the fulfillment of our classification objective .
This article introduces a new and powerful PCA-based approach that differs from current PCA approaches in several critical ways. First, our proposed method (PM) determines and exploits assay-specific gene contribution from the complete set of PCs, that is, both eigenassays (EG) and eigengenes (EA) PCs . To our knowledge this is the first application of an assay-specific gene contribution approach in this context. The score for gene i obtained from an EG is equivalent to contribution over all the assays for this gene. However, we are not aware of work in this context that uses a predefined subset of the score to obtain the signature as we do in this work. In contrast, for EA, gene contribution is not related to an EA score, and thus, a totally new application. Secondly, in contrast to current PCA methods, the PM does not rely on the existence of structure (i.e., gene clustering) in two dimensional score plots of dominant principal components (PCs) because they can be weak or even absent. Thirdly, this is the only PCA method, to our knowledge, that determines gene-ranked assay-specific signatures in a final step from the set of all genes, in contrast to methods that use reduced sets that have the possibility of removing critical genes at some intermediate step. This new PCA-based approach is presented using the following outline. The next section discusses concepts of PCA important to the understanding of the PM. The PM is then described in a separate section. Finally, three case studies are presented aiming at: (1) demonstrating the application and strengths of the PM along with limitations of current methods; (2) comparing the PM to an approach with similar features on an actual data set; and (3) applying the PM to study a dataset not explored before with this type of tool. Concluding remarks are given in the final section.
A DNA microarray context for PCA
The purpose of this section is to discuss fundamental concepts of PCA in the context of microarray gene expression data that are important to the understanding of the PM. In a typical DNA microarray experiment the expression level of m genes is surveyed for n environmental conditions and/or genetic backgrounds (called assays in this article following the work of Wall et al ). We define X as an m by n matrix of expression data with the n assays expressed along the columns and the m genes expressed along the rows with m > > n Thus, x ij is the expression level of the ith gene in the jth assay. For convenience of discussion we let the rank of X = n. In this article, the variables of any matrix are expressed along the columns and the measurements are expressed along the rows. Thus, for X, the variables are the assays and the genes are the measurements. In contrast, for XT, the variables are genes and the measurements are the assays.
PCA on X has the effect of obtaining an n dimensional orthogonal space such that the composite variability (i.e., change) of the genes is maximized in the first principal direction, and maximized in the second principal direction, with the variability captured in the first principal direction removed, and so on, until the last principal direction has the smallest variability. The n principal components (PCs) can be determined from either the covariance matrix or the correlation matrix of the variables of X. (This work uses only the correlation matrix.) The PCs are the eigenvectors of the covariance (scaled sums of squares and cross products) or the correlation (sums of squares and cross products from standardized data) matrix, and are ranked by their corresponding eigenvalues with the largest one corresponding to the first principal component (PC1). The ratio of eigenvalue i to the sum of all the eigenvalues is the proportion of the total variability captured by the i th PC. Note that the n PCs based on X have dimensions n by 1 (i.e., n elements) and following Alter et al.  they will be called eigengenes (EGs) in this article to distinguish them from the PCs based on XT, which will be called eigenassays (EAs) following Wall et al.  who departed from the "eigenarray" terminology introduced by Alter et al. . Note that the EAs have dimensions m by 1. The elements (i.e., coefficients) of a PC are called loadings. A PC is a linear combination of the variables and the projection of the i th row of measurements onto a PC (i.e., the scalar product) gives its i th score. Let V (n by n) and U (m by n) be the loading matrices for X and XT, respectively. Also, let Sv (m by n) and Su (n by n) be the score matrices for X and XT, respectively. Figure 1 illustrates these concepts and matrices visually. As shown, the EG and EA are given by vj and uj, respectively. Also note that, following Wall et al. , gi is the transcriptional response of the i th gene and the aj's are the assay expression profiles, j = 1, . . ., n. Next we present the PM using the notation and description shown in Figure 1.
Results and discussion
Proposed method (PM)
In microarray data studies, the assay profiles are part of the experimental design. More specifically, most DNA microarray studies seek to identify assay-specific signatures (i.e., which genes determine a given phenotype as expressed by the assay). The experimenter chooses the assays because of natural relationships, a hypothesis, or both. Consequently, a considerable amount of á priori knowledge is known about the assay profiles before any analysis on the data. However, often times hidden relationships that advance knowledge are brought to light by PCA. The PM exploits the expected assay profile behavior and new knowledge to obtain assay-specific gene expression profiles or signatures with a ranking of all the genes in a profile. Our PM seeks to accomplish this objective using the following basic steps after obtaining the PCs (i.e., EGs and EAs):
Determine the EGs for X (V) and EAs for XT (U) based on their correlation matrices.
Determine the pseudo score matrix for EA (Su) using U in Step 1 and X, i.e., XTU. X is not standarized.
Plot the EG loadings (the row coefficients of V) and pseudo EA scores (the row coefficients of XTU) against the assay index starting from the most dominant to the lesser dominant PCs until the assay signature of interest is revealed as a group of large outliers.
Determine the contribution for each gene using the selected PC in Step 1 for the assay signature group.
Rank the genes by contribution and plot the contributions in ranked order against an index that increases by one unit for each gene. Gene identification with rank should be maintained.
Obtain the assay-specific gene signatures from the ranked assay profiles using the plot in Step 3 to determine cut-offs or limits on the size of profile groups.
Most PC software packages give an option for selecting between the covariance or correlation matrix of the variables. For standardization of scale, we recommend using the correlation option in Step 1. However, for obtaining the score matrices in Step 2, our procedure is to standardize X for obtaining EG scores but not for EA scores since we do not standardized X when determining EA gene contribution. Even though the EA are determined from the correlation matrix, we do not standardize X because under EA standardization X is row centered and since the contribution of each gene i is proportional to the sum across row i, these contribution values tend to disappear under standardization. The EG standardization of X is given below:
Note that the scores matrix for EG (Sv) is determined by Sv = ZV.
The objective of Step 3 is to determine and select the PC with the strongest character of the assay group of interest using loading plots for the EGs and score plots for the EAs. These plots are compared against the á priori knowledge of the assays. In addition, these plots can also bring to light relationships among assays not known or hypothesized á priori. The strength of a specific profile is determined by its PC association, the amount of variation explained by the associated PC, and its separation from the other assay profiles. The PC associated with the plot with the strongest profile is selected to determine the signature. This can be either an EG or an EA. Our experience have found that a number of EA and EG PCs provide similar information but that either one can provide unique information. This is why we give the recommendation to examine both types of eigenvectors.
Steps 4–6 involve the development, analysis, and use of gene contribution plots. To our knowledge this is the first use of contribution plots for determining signatures. This approach has the advantage over a loading approach because the loadings only give the relative weight for a specific PC but the contribution gives the product of the relative weight and the expression level for a particular gene. In addition, our approach has the advantage of determining this contribution for only the members in the assay group of interest whereas the loading gives an effect across all the assays. The contribution for a specified assay group, GP, for gene i, using EG vk is
Similarly, using EA uk, is
If the absolute value of
is high, then we are assuming that it is critical to the assay-specific signature and the greater , the more critical the gene. Thus, by ranking the genes from the highest to the lowest , we are able to obtain a ranked list of genes for each assay-specific signature. The cut-off point is guided by the level of change across the gene index revealed in the ranked contribution plots. Nonetheless, since the genes are ranked, no matter where one makes the cut-off, the list will be the strongest one for the selected number of genes. We now apply the PM to three cases; a simulation study of artificial data and two cases of real data and demonstrate its ability to obtain assay-specific gene signatures of ranked order.
The case studies
In this section we present three case studies to evaluate our PM. The first study uses simulated data based on the study of Wall et al. . This study will validate the PM by demonstrating its ability to obtain the correct assay-specific gene signatures. The second study involves applying the PM to a real DNA microarray data set (the case studied by Misra et al. ). This study consisted of 6,972 genes and 40 normal human tissue samples, which they used to develop and apply a novel method for obtaining tissue-specific gene expression signatures. In this case study we discuss the limitations of their method as compared to the strengths of the PM and compare results. The last study uses real expression data of 4,290 identified genes and twelve assays, representing a combination of two recombinant E. coli strains and different cultivation conditions including the alternative use of two sugars and the exposure of the cells to different ethanol concentrations (the case studied by Gonzalez et al. ). We analyze this data set using our PM and present gene expression profiles for the most relevant assay-specific signatures: i.e., the response of the cells to a 2% ethanol challenge.
The simulation data study
This study generates artificial gene expression data mimicking the simulation study in Wall et al.  to introduce and illustrate the PM. This data have three kinds of transcriptional responses: 1. noisy (genes 1 to 1,600); 2. noisy sinusoidal behavior (genes 1601 to 1,800); and 3. noisy exponential pattern (1,801 to 2,000). The added noise was distributed normally with mean 0 and standard deviation 0.5. The sine pattern was
where a is distributed uniformly on the interval (1.5, 3). The exponential pattern was
where b is distributed uniformly on the interval (4, 8) and t is the time (in minutes). The number of assays is fourteen with assay i corresponding to sampling at time 10i min, i = 0, 1, . . ., 13. To verify that our data were in agreement with the Wall et al.  data, we plotted the graphs in Figure 2. These plots are in agreement with their Figure 5.3. As we discussed in previous section, a two-dimensional score scatter plot based on EG1 and EG2 is common in PCA analysis. Although Figure 3 separates three clusters of genes in agreement with the nature of this data set, it does not provide knowledge of the gene or assay characteristics for the clusters and is therefore not used by the PM. In this study we are assuming that the experimenter has designed the sinusoidal and exponential patterns into the assays. Thus, the experimenter has a complete knowledge of this behavior. For space considerations we will not show these plots but they would look like the ones in Figure 2, except that there would be one point plotted for each assay number. The PM will be able to effectively accomplish its objective if it can capture and reveal these patterns in a few dominant PCs, preferably PC1 and PC2.
In applying the PM, we first determine the PC matrices U and V as required by Step 1. We then determine the pseudo scores for the EAs per Step 2 and produce Figure 4 under Step 3. This figure plots the loading for EG and the pseudo scores for EA against the assay time index. Sinusoidal patterns are clearly revealed by all the curves although the one for EG2 appears to be distorted. The EA2 plot appears to match the sinusoidal behavior in Figure 2 the best. The EG1 loading plot and the EA1 score plot appear to have the right sinusoidal shape but with a larger period. Closer examination of these plots, especially the EA1 score plot, indicates evidence of a composite behavior of both patterns. EG1 and EG2 explained about 51% and 13% of the total variation, respectively, and EA1 and EA2 explained about 21% and 9%, respectively. In our research, we have found much lower amounts for EAs to be common. However, one should not use this as a criterion for choosing an EG analysis over an EA analysis. As discussed previously, the PC selection criteria that we recommend are the separation and grouping of assay groups as well as the order of the PC.
In applying Step 4, at least two types of assay profiles or groups need to be defined to identify the sinusoidal and exponential genes. For this particular case, we chose to do this by defining these groups based on positive and negative loadings and scores in Figure 4. The results of this grouping are given by the score contribution plots in Figure 5. In these plots, the contribution of gene i, for a given PC, is the sum of the scores for all the genes for the particular grouping (either positive or negative). EA1 plot is very effective separating both induced gene types (1,601–2,000) from the noisy genes (Figure 5C), while the EA2 plot confirms that this PC represents only the sinusoidal genes (Figure 5D), although the identification is not so distinct for a significant number of genes. Although EG1 appears to represent both gene types and does an excellent job of identifying as one group all 400 genes (Figure 5A), the EG2 score contribution plot appears to represent the sinusoidal genes for the positive (+) assay group and the exponential genes for the negative (-) assay group (Figure 5B). In summary, the EA analysis is a better choice for this data set than an EG analysis since it better identifies both types of genes using EA1 and sinusoidal genes using EA2.
After separating the genes by their PC assays groups, Step 5 is applied. The objective of this step is to develop a ranked-order list for all the groups identified in Step 4. To do this, we again use the contributions obtained in Step 4 and appearing in Figure 5 but now we plot them in ascending or descending order against an index that increases by one unit for each gene. The ranked gene plots for both EG and EA results are given in Figure 6. Both plots look very similar except for low index numbers. For the EG plot, the lines for the two groups cross and separate, thus creating a false gene signature for the lower ranking genes (i.e., this is not true since they are all noisy genes). The EA plot however does not create such a false signature since both lines are close to zero as the rank decreases and is thus the better choice here also. As the rank increases, the shape of the curves permit the identification of two groups of genes clearly separated by two inflection points around gene numbers 1,600 and 1,800. These inflection points can be obtained by identifying the optimum (minimum or maximum) of the first derivative. We show this calculation for the EAs plot in Figure 6 and, as expected, two inflection points are observed at gene numbers 1,600 and 1,800. This feature of the gene contribution plots can be used to precisely define complete sets of genes that represent an assay signature. It is noteworthy to mention that the only sharp changes in this curve take place in the transition from the noisy to induced genes thus clearly indicating the start of the signature (i.e., induced genes). No sharp breaks occurred within the region of induced genes, indicating the truth that no gene in this expression signature is dominant.
Human tissue expression data case study
Our second study compares the results and methodology of the PM with another one (Misra et al. , referred to as the "Misra et al." method in what follows) of similar ability (i.e., a PCA-based method) in a study of real microarray data. Their approach represents the most refined use of PCA for the analysis of FG data that we found in the literature. This data set consists of assays from several human tissues from the brain, kidney, lung, esophagus, skeletal muscle, breast, stomach, colon, blood, spleen, prostate, testes, vulva, proliferative endometrium, myometrium, placenta, cervix and ovary [see Additional file 1]. We will first briefly describe their approach and compare and contrast it to the PM. Then we apply the PM and determine tissue-specific gene expression signatures for brain, liver, and muscle assays. We close this section by comparing the signature results of both methods.
Although it is not clear from their loading and score plots, we determined that the PCA method of Misra et al. is based on eigengenes (EGs). Their method is basically a gene screening method where genes are removed at various reduction steps while seeking to maintain the structure of the loading plot. Their first reduction step is a course screening and reduces the total set of genes to a workable subset (from about 7,000 genes to 425) using a filtering method that eliminates genes with loadings below a threshold. One limitation is that this threshold is not based on statistical significance but is somewhat subjective, although quantitative. This set of genes is then reduced to unrefined signatures without identities at this point using a histogram and visually determining the classes (i.e., the gene groups). These unrefined gene groups are refined (i.e., further reduced) using cluster analysis and visually eliminating genes outside subjectively chosen cluster groups. The final step uses score plots to "reveal" the "nature of each gene group" and to rank the genes within each group.
The Misra et al. method is similar to our PM in that the objective is to identify signature groups and prioritize the genes in these groups. However, there are some very critical differences. First, they only make use of EGs in contrast to the PM that compare EGs to EAs and selects the one that best represents the assay group of interest. Another critical difference is that their method relies on structure in two dimensional loading plots but the PM does not since, as we determined by simulation (previous section) and will show in the ethanol study (next section), this type of plot will not always separate genes into identifiable groups (i.e., have useful structure). Furthermore, for the Misra et al. method, gene contribution is based on the loading of a gene whereas in the PM gene contribution is the sum of the product of the loading and the expression level for only the assays in the signature group (see Eqs. 2–3). Another critical difference is that when the Misra et al. method ranks the genes within a gene group only the "refined" genes are candidates; and thus, not giving a chance to those not in this group that could have been incorrectly eliminated. In contrast, when the genes are ranked by the PM all the genes are considered. Finally, the PM seems to be much simpler to apply as we now demonstrate on this study.
Figure 7 gives the loading plots for EA1–EA3 and Figure 8 gives the score plots for EA1 to EA3. EG1–EG3 explained about 60%, 7%, and 6% of the total variation, respectively, and EA1–EA3 explained about 12%, 8% and 7%, respectively. EG1 separates the muscle and liver assays well, the only drawback is that they are on the same side of the cluster of points. This drawback is also true of EG2 in the separation of the muscle and brain assays. However, this is not true for EG3 in the separation of the muscle and liver assays but the drawback of this plot is that the separations are not very prevalent. We did not choose any of these PCs for signature definition because we found some EAs with better properties.
Figure 8 contains the score plots for the first three EAs. For the EA1 plot, from the positive values, the first PC appears to reflect a combination of Placenta 2 (gray dot), the three muscles (squares), the six brain (diamonds), Breast 9 (the large black dot), and the two liver (the triangles) tissues. All these groups are on the same side of the cluster and, thus we did not use this PC for assay development. In contrast, EA2 separates nicely into muscle and brain assays on opposite sides of the cluster. Similarly, EA3 separates nicely into only a liver assay group. Hence, we selected EA2 to develop the muscle and brain signatures and EA3 to develop the liver signatures.
The contribution plots ranking the genes within the signatures using EA2 and EA3 are shown in Figure 9A–C. As one can see, the EAs produced well defined signatures with sharp, steep rises that significantly distinguish contributions for the most dominant genes. As in the simulation case study, an inflection point can be identified which marks the start of a tissue-specific signature (see insets in Figure 9A–C). Using this approach, 3,014, 3,492, and 2,647 genes are identified as the complete brain, muscles, and liver signatures, respectively. This list is clearly too comprehensive and of limited practical use due to its length. However, because of the sharp changes over a small number of genes as the highest ranking genes are approached, a subset of genes contributing the most to each tissue-specific signature can be clearly observed.
We identified the cut-off point that marks the limit of this subset of genes by using the curvature (κ) as a metric to assess the steepness of the curve. For a function y = f (x), the curvature is defined as
When the curve is changing sharply κ is large and when the curve changes slowly κ is small. Thus, by calculating the gene index that produces a maximum value of κ, a clear cut-off in the signature can be identified. We calculated κ for the sharp regions in the ranked gene contribution plots, and in all cases there is a clear maximum defining the cut-off points (see Figure 9D). Using this approach the size of the subset of genes defining the brain, liver, and muscles signatures was 67, 107, and 99, respectively. Individual genes in the brain-specific signature are shown in Table 1, and Tables 2 and 3 show the results for the liver and muscle signatures [for complete tables, including gene function, see Additional file 2].
The PM clearly identifies the genes reported by Misra et al.  as part of each tissue-specific signature (genes reported in their study are shown in bold in these tables with their original ranking in parenthesis). We also investigated whether the genes newly identified using the PM are indeed representative of each tissue by reviewing current information available at "Entrez Gene" (http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=gene), a database from the National Center for Biotechnology Information  and its multiple links to the latest reports available in the literature. Genes marked with an asterisk in Tables 1, 2, 3 are deemed to be good markers of the specific tissue because previous studies have shown them to be either functionally associated with the health/disease state and/or expressed at high levels in the tissue. Clearly, a very large proportion of the newly identified genes are in fact good markers of the individual tissues. In addition, the identification of new genes constitute the basis to formulate hypothesis regarding their involvement in the functioning of the specific tissue, thus evidencing the potential use of the PM to identify new functions. For example, a very large proportion of newly identified genes in the muscle signature (see Table 3) encode ribosomal proteins, which would indicate their differential expression in this tissue. This is in fact in agreement with existing evidence that supports the tissue-specific expression of ribosomal proteins . Using additional PCs (actually only EGs), we were able to develop signatures for a number of assay tissue groups. For space considerations, we do not present them in this article but for the interested reader they can be found at the website: http://www.public.iastate.edu/%7edrollins/data.html .
The ethanol response in E. coli case study
After demonstrating the superior performance of the PM in both a simulation study and a typical analysis of DNA microarray data using existing PCA methodologies (two previous sections), we now present the application of the PM to investigate a data set that had not been analyzed using any PCA-based method; i.e., the study of ethanol tolerance in ethanologenic E. coli strains KO11 and LY01 by Gonzalez et al.  [see Additional file 3]. The aim of their study was to identify the basis for ethanol tolerance in ethanol-resistant strain LY01 by comparison to its parent strain KO11 using genome-wide transcriptional profiling. This data set consists of 12 assays, which are a combination of two genetic backgrounds (strains KO11 and LY01) and different variations of cultivation conditions such as presence/absence of ethanol (as 1% in the initial medium and as a 2% ethanol challenge), and two different types of sugars (glucose or xylose). Table 4 summarizes the most important features of each assay in this data set. The expression levels for 4,290 genes were surveyed for each assay in this study. It is noteworthy to mention that their analysis focused on comparing the two genetic backgrounds, thus trying to establish strain-specific signatures.
The two-dimensional score plot with EG2 vs. EG1 is shown in Figure 10. From a visual examination of this plot, there does not appear to be any grouping of genes. Therefore, this appears to be another case where this plot is not useful due to a lack of gene clusters. We apply the PM beginning with Step 1. The EG loading and EA score plots are given in Figure 11 for the first three PCs. EG1, EG2, and EG3 explained about 80%, 13%, and 2% of the total variation, respectively, and EA1, EA2, and EA3 explained about 55%, 13%, and 10%, respectively. Comparison of the two plots suggests that the EA does the best job of grouping the assays. Both EG1 and EG2 separate assays 1–8 from assays 9–12: i.e., distinguish the 2% ethanol-challenged cells (9–12) from the non-ethanol challenged cells (1–8) (see also Table 4). Although EG1 explains about 80% of the total variation it just marginally separated these two groups of cultures. EG3, in a very undistinguishable manner, appears to account for a combination of the effects of strains and type of sugar.
The EA plot (Figure 11B) however distinguishes three different characteristics of the assays in a very clear way. EA1, which explains the largest fraction of the total variation, made a clear distinction between ethanol and non-ethanol-challenged cultures. EA2 on the other hand, further separates non-ethanol challenged cultures into those with no ethanol present in the initial medium (1–4) from those containing 1% ethanol in the initial medium (5–8) (see also Table 4). The latter group is further separated into two cultures containing the sugar glucose (5–6) and two others containing xylose (7–8). Finally, EA3 represents a strain-specific signature as it clearly separates strain KO11 (assays 1, 3, 5, 7, 9, 11) from strain LY01 (2, 4, 6, 8, 10, 12). In summary, the EA-based analysis is clearly superior here, and we therefore chose to do only an EA analysis. In addition, we decided to focus on the results for the first PC (i.e., EA1), which represents the assay signature corresponding to the response of the cultures to a 2% ethanol challenge. Therefore, the assay-specific signature we are exploring is clearly different from the one explored by the study of Gonzalez et al. ; i.e., while we are studying the ethanol signature their study focused on the strain signature. From the EA score plot in Figure 11B, the non-ethanol scores are all negative and closer together. Similarly, the ethanol scores are all positive and closer together. Therefore, only one PC (EA1) and two assay profiles from this PC are sufficient for determining the gene signatures for ethanol and non-ethanol challenged cultures. The negative scores for EA1 make up the non-ethanol profile and the positive scores make up the ethanol profile.
In Figure 12 we plot the gene contributions of the non-ethanol assay followed by the ethanol assay against the order of the genes. This plot shows that the non-ethanol profile has about twice as much spread as the ethanol profile. In addition, the non-ethanol spread is skewed towards the negative numbers and the ethanol spread is skewed towards the positive numbers, which explains the negative scores for non-ethanol assays and the positive scores for ethanol assays. The negative values that separate from the non-ethanol cluster make up the gene signature for the non-ethanol. Similarly, the positive values that separate from the ethanol cluster correspond to the gene signature for ethanol. The ranked contribution plots in Figure 13 allow us to find the gene signatures for these two groups of assays under EA1.
As for the two previous case studies, inflection points and sharp regions in the ranked contributions plots in Figure 13 clearly define signature groups. Therefore, the identification of the group of genes composing the ethanol and non-ethanol signatures is straightforward as described in previous sections: i.e., cut-offs are identified by calculating points of inflection and maximum curvature. Using this approach, 702 and 522 genes were identified as part of the complete ethanol and non-ethanol signatures, respectively. In addition, two subsets of genes were identified that contributed the most to each signature (see identification of the point of maximum curvature in insets of Figure 13), which include 126 and 51 genes in the ethanol and non-ethanol signatures, respectively. The top ranked genes are labeled in Figure 13, and Table 5 shows the top ten genes in each signature along with their rank and expression ratios. Note that genes defining the ethanol signature are in fact up-regulated (expression ratios > 1) in response to the ethanol challenge while those composing the non-ethanol signature are down-regulated (expression ratios < -1). Although highly ranked genes frequently exhibit a larger expression ratio, this is not always true because in the PM a gene's contribution is based on product of its loading and expression level and not just its loading. Tables 6 and 7 give the signatures in terms of the genes contributing the most as identified by the points of maximum curvature in Figure 13 [for complete tables, including gene function, see Additional file 2]. We examined the functions encoded by this group of genes (Table 8) and found them very revealing of the metabolic rearrangements associated with the response of the cells to an ethanol challenge.
As expected for the immediate response to an environmental challenge, the cells significantly modified functions involved in regulatory, transport, and general processes (i.e., the largest changes in the "known functions" category in Table 8 are observed for these three functional groups), and important changes appear to take place at the translation and post-translational modification levels. The same trend was observed for the "putative functions" categories, where most populated functions were those involved in regulation and transport. The group of genes composing these signatures also provides very important information to formulate hypothesis about which specific gene(s)/function(s) are involved in the cellular response to an ethanol challenge. For example, many of the top-ranked genes (Table 5) encode functions that one would expect to be involved in the cellular response to an ethanol challenge such as the metabolism and transport of osmolytes (mngA, cvrA, and caiA), the biosynthesis of phospholipids (plsC) which are major constituents of the cell membrane, and the repairing of misfolded proteins (dsbC). In fact, increased tolerance to ethanol in certain E.coli strains is related to the increased availability of osmolytes like betaine and trehalose [16, 19]. In summary, the application of the PM to this data set allowed the identification of a group of genes representing the ethanol signature, which can be hypothesized to be involved in the cellular response to an ethanol challenge. Such hypothesis form the basis of current studies in our groups aimed at elucidating the mechanisms/processes involved in the response to ethanol in E. coli.
This article proposes a new and powerful PCA-based method for the identification of assay-specific gene signatures of ranked order in the analysis of FG data. This method is unique for several reasons. First, it is the only one, to our knowledge, that uses gene contribution, a product of the loading and expression level, to obtain assay signatures. Our proposed method (PM) develops and exploits two types of assay-specific contribution plots. To our knowledge the development and use of these plots is new to the application of PCA in the FG area. The first type plots the assay-specific gene contribution against the given order of the genes and reveals variations in distribution between assay-specific gene signatures as well as outliers within assay groups indicating the degree of importance of the most dominant genes. The second type plots the contribution of each gene in ascending or descending order against a constantly increasing index. This type of plots reveals assay-specific signatures defined by the inflection points in the curve. In addition, sharp regions within the gene signature in the ranked contribution plots define the genes that contribute the most to the signature. We proposed and used the curvature as an appropriate metric to characterize sharp regions of these plots, thus identifying the subset of genes contributing the most to the signature. Secondly, we know of no other method that selects an assay group by comparing EG loadings against the assays and EA scores against the assays. It is worth noting that no set of PCs (EGs or EAs) will necessarily contain all the information for grouping assays. For example, for the human tissue data in our second study, although EA2 and EA3 were best for obtaining the muscle, brain and liver signatures, the EGs provided better assay groupings for Colon 5; Lungs 1 and 2; Lungs 4 and 5; Stomach 1; Breast 7 and 9; Kidneys 7–10; Placenta 2 and 3; Vulva "*", 1 and 2; Blood; and Cervix 2. As mentioned earlier, these signatures are available at the website: http://www.public.iastate.edu/%7edrollins/data.html . Finally, the PM uses the full dataset to determine the final gene signature, thus eliminating the chance of gene exclusion by poor screening in earlier steps.
This article presented the PM in six basic steps and applied it to three different case studies: one artificial and two real data sets. In the artificial data study, which consisted of 1,600 noisy genes, 200 noisy sinusoidal genes, and 200 noisy exponential genes, the PM identified all 400 genes with patterns in a sinusoidal/exponential-specific gene signature and most of the genes with sinusoidal patterns in a sinusoidal-specific gene signature. In the first study involving real data, we compared the results and methodology of the PM with another one (the method of Misra et al. ) of similar ability (i.e., a PCA-based method) in a study of real DNA microarray data. Their approach represents the most refined used of PCA for the analysis of FG data that we have found in the literature. The tissue-specific signatures identified by the PM not only included the genes previously identified by Misra et al.  but also added genes that are known to be linked to the specific tissue and other of unknown connection, thus establishing hypothesis regarding their potential involvement in the functions of those tissues. In the second study involving real data, the ethanol case, we were able to explore a data set that had not been previously analyzed using any PCA-based approach. The PM was able to identify different assay-specific signatures including the ethanol- and strain-signature. A detailed study of the ethanol signature using the PM resulted in the clear identification of a group of genes that are relevant to the response of E. coli to an ethanol challenge. These findings confirm the capability of the PM to generate testable hypothesis that will contribute to elucidating different biological processes (in this example the basis of the response of the cells to an ethanol challenge).
principal component analysis
Bernal A, Ear U, Kyrpides N: Genomes OnLine Database (GOLD): a monitor of genome projects world-wide. Nucleic Acids Res 2001, 29: 126–127.
Hieter P, Boguski M: Functional genomics: it's all how you read it. Science 1997, 278: 601–602.
Dharmadi Y, Gonzalez R: DNA Microarrays: Experimental Issues, Data Analysis, and Application to Bacterial Systems. Biotechnol Prog 2004, 5: 1309–1324.
Stoyanova R, Querec TD, Brown TR, Patriotis C: Normalization of single-channel DNA array data by principal component analysis. Bioinformatics 2004, 20: 1772–1784.
Stevenson MD, Chan V, Gustafson S, Kelley-Loughnane N, Harker B, Wang C, Rudnicki D, Hussain S, Frazier J: Comparative study of DNA microarray data analysis: Principal component analysis verses fisher linear discriminant analysis. Toxicol Sci 2003, 72: 92.
Crescenzi M, Giuliani A: The main biological determinants of tumor line taxonomy elucidated by a principal component analysis of microarray data. FEBS Lett 2001, 507: 114–118.
Verhoeckx KC, Bijlsma S, de Groene EM, Witkamp RF, van der Greef J, Rodenburg RJ: A combination of proteomics, principal component analysis and transcriptomics is a powerful tool for the identification of biomarkers for macrophage maturation in the U937 cell line. Proteomics 2004, 4: 1014–1028.
Marengo E, Leardi R, Robotti E, Righetti PG, Antonucci F, Cecconi D: Application of three-way principal component analysis to the evaluation of two-dimensional maps in proteomics. J Proteome Res 2003, 2: 351–360.
Bryant DK, Monte S, Man WJ, Kramer K, Bugelski P, Neville W, White IR, Camilleri P: Principal component analysis of mass spectra of peptides generated from the tryptic digestion of protein mixtures. Rapid Commun Mass Spectrom 2001, 15: 418–427.
Desbrosses GG, Kopka J, Udvardi MK: Lotus japonicus metabolic profiling. Development of gas chromatography-mass spectrometry resources for the study of plant-microbe interactions. Plant Physiol 2005, 137: 1302–1318.
Cheng LL, Burns MA, Taylor JL, He W, Halpern EF, McDougal WS, Wu CL: Metabolic characterization of human prostate cancer with tissue magnetic resonance spectroscopy. Cancer Res 2005, 65: 3030–3034.
Wu H, Zhang X, Li X, Li Z, Wu Y, Pei F: Comparison of metabolic profiles from serum from hepatotoxin-treated rats by nuclear-magnetic-resonance-spectroscopy-based metabonomic analysis. Anal Biochem 2005, 340: 99–105.
Misra JW, Hwang D, Hsiao LL, Gullans S, Stephanopoulos G, Stephanopoulos G: Interactive exploration of microarray gene expression patterns in a reduced dimensional space. Genome Res 2002, 12: 1112–1120.
Wall ME, Rechtsteiner A, Rocha LM: Singular Value Decomposition and Principal Component Analysis. In A Practical Approach to Microarray Data Analysis. Edited by: Berrar DP, Dubitzky W, Granzow M. Massachusetts: Kluwer Academic Publishers; 2003:91–109.
Alter O, Brown PO, Botstein D: Singular Value Decomposition for Genome-Wide Expression Data Processing and Modeling. Proc Natl Acad Sci USA 2000, 97: 10101–10106.
Gonzalez R, Tao H, Purvis JE, Shanmugam KT, York SW, Ingram LO: Gene Array-Based Identification of Changes That Contribute to Ethanol Tolerance in Ethanologenic Escherichia coli : Comparison of KO11 (Parent) to LY01 (Resistant Mutant). Biotechnol Prog 2003, 19: 612–623.
Maglott D, Ostell J, Pruitt KD, Tatusova T: Entrez Gene: gene-oriented information at NCBI. Nucleic Acids Res 2005, 33: D54-D58.
Bortoluzzi S, d'Alessi F, Romualdi C, Danieli GA: Differential expression of genes coding for ribosomal proteins in different human tissues. Bioinformatics 2001, 17: 1152–1157.
Purvis JE, Yomano LP, Ingram LO: Enhanced Trehalose Production Improves Growth of Escherichia coli under Osmotic Stress. App Environm Microbiol 2005, 71: 3761–3769.
Riley M, Serres MH: Interim Report on Genomics of Escherichia coli . Annu Rev Microbiol 2000, 54: 341–411.
Bebien M, Kirsch J, Mejean V, Vermeglio A: Involvement of a putative molybdenum enzyme in the reduction of selenate by Escherichia coli . Microbiology 2002, 148: 3865–3872.
Brokx SJ, Ellison M, Locke T, Bottorff D, Frost L, Weiner JH: Genome-wide analysis of lipoprotein expression in Escherichia coli MG1655. J Bacteriol 2004, 186: 3254–3258.
Cusa E, Obradors N, Baldoma L, Badia J, Aguilar J: Genetic analysis of a chromosomal region containing genes required for assimilation of allantoin nitrogen and linked glyoxylate metabolism in Escherichia coli . J Bacteriol 1999, 181: 7479–7484.
Dolzan MK, Roig-Zamboni V, Campanacci V, Tegoni M, Schneider G, Cambillau C: Crystal structure and reactivity of YbdL from Escherichia coli identify a methionine aminotransferase function. FEBS Lett 2004, 571: 141–146.
Drew D, Sjostrand D, Nilsson J, Urbig T, Chin CN, de Gier JW, von Heijne G: Rapid topology mapping of Escherichia coli inner-membrane proteins by prediction and PhoA/GFP fusion analysis. Proc Natl Acad Sci USA 2002, 99: 2690–2695.
Gardner AM, Gessner CR, Gardner PR: Regulation of the nitric oxide reduction operon (norRVW) in Escherichia coli . Role of NorR and sigma54 in the nitric oxide stress response. J Biol Chem 2003, 278: 10081–10086.
Gomes CM, Giuffre A, Forte E, Vicente JB, Saraiva LM, Brunori M, Teixeira M: A novel type of nitric-oxide reductase. Escherichia coli flavorubredoxin. J Biol Chem 2002, 277: 25273–25276.
Heidrich C, Templin MF, Ursinus A, Merdanovic M, Berger J, Schwarz H, de Pedro MA, Holtje JV: Involvement of N-acetylmuramyl-L-alanine amidases in cell separation and antibiotic-induced autolysis of Escherichia coli . Mol Microbiol 2001, 41: 167–178.
Kim C, Song S, Park C: The D-allose operon of Escherichia coli K-12. Escherichia coli 1997, 179: 7631–7637.
Kim GJ, Lee DE, Kim HS: Functional expression and characterization of the two cyclic amidohydrolase enzymes, allantoinase and a novel phenylhydantoinase, from Escherichia coli . J Bacteriol 2000, 182: 7021–7028.
Ladner JE, Obmolova G, Teplyakov A, Howard AJ, Khil PP, Camerini-Otero RD, Gilliland GL: Crystal structure of Escherichia coli protein YbgI, a toroidal structure with a dinuclear metal site. BMC Struct Biol 2003, 3: 7–14.
Lubitz SP, Weiner JH: The Escherichia coli ynfEFGHI operon encodes polypeptides which are paralogues of dimethyl sulfoxide reductase (DmsABC). Arch Biochem Biophys 2003, 418: 205–216.
Miller BG, Raines RT: Identifying latent enzyme activities: substrate ambiguity within modern bacterial sugar kinases. Biochemistry 2004, 43: 6387–6392.
Oussenko IA, Sanchez R, Bechhofer DH: Bacillus subtilis YhaM, a member of a new family of 3 '-to-5 ' exonucleases in gram-positive bacteria. J Bacteriol 2002, 184: 6250–6259.
Pernestig AK, Melefors O, Georgellis D: Identification of UvrY as the cognate response regulator for the BarA sensor kinase in Escherichia coli . J Biol Chem 2001, 276: 225–231.
Poulsen TS, Chang YY, Hove-Jensen B: D-Allose catabolism of Escherichia coli : involvement of alsI and regulation of als regulon expression by allose and ribose. J Bacteriol 1999, 181: 7126–7130.
Reed JL, Vo TD, Schilling CH, Palsson BO: An expanded genome-scale model of Escherichia coli K-12 (iJR904 GSM/GPR). Genome Biol 2003, 4: R54.
Saurin W, Hofnung M, Dassa E: Getting in or out: early segregation between importers and exporters in the evolution of ATP-binding cassette (ABC) transporters. J Mol Evol 1999, 48: 22–41.
Utsumi R, Kawamoto K, Yamazaki K, Taniguchi M, Yoshioka S, Tanabe H: Characterization of the signal transduction via evgS and evgA in Escherichia coli . J Gen Appl Microbiol 1996, 42: 155–162.
Verheul A, Wouters JA, Rombouts FM, Abee T: A possible role of ProP, ProU, and CaiT in osmoprotection of Escherichia coli by carntine. J Appl Microbiol 1998, 8: 1036–1046.
Verkhovskaya ML, Barquera B, Wikstrom M: Deletion of one of two Escherichia coli genes encoding putative Na+/H+ exchangers ( ycgO ) perturbs cytoplasmic alkali cation balance at low osmolarity. Microbiology 2001, 147: 3005–3013.
Yew WS, Gerlt JA: Utilization of L-ascorbate by Escherichia coli K-12: assignments of functions to products of the yjf-sga and yia-sgb operons. J Bacteriol 2002, 184: 302–306.
Yim L, Martinez-Vicente M, Villarroya M, Aguado C, Knecht E, Armengod ME: The GTPase activity and C-terminal cysteine of the Escherichia coli MnmE protein are essential for its tRNA modifying function. J Biol Chem 2003, 278: 28378–28387.
We would like to thank Profs. G. Stephanopoulos and L. O. Ingram for making available their data sets. The assistance of Rhonda DeCook and Tanzy Love in the Iowa State University Department of Statistics in using statistical software packages is acknowledged. This project was supported by grants from the U.S. National Science Foundation (BES-0331388 and BES-0601549).
D.K.R. and R.G. designed research; D.K.R., R.G., D.Z., A.L.J., A. M., and J.W.G. performed research; D.K.R., R.G., J.W.G., and A. M. wrote the paper.
Electronic supplementary material
Additional File 1: Human tissue expression data case study. Data sets used in the "human tissue expression data case study". (XLS 2 MB)
Additional File 2: Gene signatures for "case studies". Complete gene signatures, including gene names and a brief description of their function, are provided for the "human tissue expression data case study" and "the ethanol response in E. coli case study". (XLS 204 KB)
Additional File 3: Ethanol response in E. coli case study. Data sets used in the "ethanol response in E. coli case study". (XLS 1 MB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
About this article
Cite this article
Rollins, D.K., Zhai, D., Joe, A.L. et al. A novel data mining method to identify assay-specific signatures in functional genomic studies. BMC Bioinformatics 7, 377 (2006). https://doi.org/10.1186/1471-2105-7-377