Skip Navigation


ToxSci Advance Access originally published online on April 6, 2006
Toxicological Sciences 2006 92(2):560-577; doi:10.1093/toxsci/kfj184
This Article
Right arrow Abstract Freely available
Right arrow FREE Full Text (PDF) Freely available
Right arrow All Versions of this Article:
92/2/560    most recent
kfj184v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Add to My Personal Archive
Right arrow Download to citation manager
Right arrowRequest Permissions
Right arrow Disclaimer
Google Scholar
Right arrow Articles by Yu, X.
Right arrow Articles by Faustman, E. M.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Yu, X.
Right arrow Articles by Faustman, E. M.
Social Bookmarking
 Add to CiteULike   Add to Connotea   Add to Del.icio.us  
What's this?

© The Author 2006. Published by Oxford University Press on behalf of the Society of Toxicology. All rights reserved. For Permissions, please email: journals.permissions@oxfordjournals.org

HIGHLIGHTED ARTICLE

A System-Based Approach to Interpret Dose- and Time-Dependent Microarray Data: Quantitative Integration of Gene Ontology Analysis for Risk Assessment

Xiaozhong Yu*, William C. Griffith*, Kristina Hanspers{dagger}, James F. Dillman, III{ddagger}, Hansel Ong*, Melinda A. Vredevoogd§ and Elaine M. Faustman*,1

* Department of Environmental and Occupational Health Sciences, University of Washington, Seattle, Washington 98105; {dagger} GenMAPP Development Team, Bioinformatics Research Associate/Conklin Lab Gladstone Institute of Cardiovascular Disease/UCSF, San Francisco, California 94158; {ddagger} Cell and Molecular Biology Branch, U.S. Army Medical Research Institute of Chemical Defense, Aberdeen Proving Ground, Maryland 21010; and § Institute for Risk Analysis and Risk Communication, Environmental and Occupational Health Sciences, Seattle, Wasington

1 To whom correspondence should be addressed at Institute for Risk Analysis and Risk Communication, Department of Environmental and Occupational Health Sciences, University of Washington, 4225 Roosevelt Way NE, Suite #100, Seattle, WA 98105. Fax: (206) 616-4875. E-mail: faustman{at}u.washington.edu.

Received March 30, 2006; accepted April 1, 2006


    ABSTRACT
 TOP
 ABSTRACT
 INTRODUCTION
 SYSTEMS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSION
 REFERENCES
 
Although microarray technology has emerged as a powerful tool to explore expression levels of thousands of genes or even complete genomes after exposure to toxicants, the functional interpretation of microarray data sets still represents a time-consuming and challenging task. Gene ontology (GO) and pathway mapping have both been shown to be powerful approaches to generate a global view of biological processes and cellular components impacted by toxicants. However, current methods only allow for comparisons across two experimental settings at one particular time point. In addition, the resulting annotations are presented in extensive gene lists with minimal or limited quantitative information, data that are crucial in the application of toxicogenomic data for risk assessment. To facilitate quantitative interpretation of dose- or time-dependent genomic data, we propose to use combined average raw gene expression values (e.g., intensity or ratio) of genes associated with specific functional categories derived from the GO database. We developed an extended program (GO-Quant) to extract quantitative gene expression values and to calculate the average intensity or ratio for those significantly altered by functional gene category based on MAPPFinder results. To demonstrate its application, we applied this approach to a previously published dose- and time-dependent toxicogenomic data set (J. F. Dillman et al., 2005, Chem. Res. Toxicol. 18, 28–34). Our results indicate that the above systems approach can describe quantitatively the degree to which functional gene systems change across dose or time. Additionally, this approach provides a robust measurement to illustrate results compared to single-gene assessments and enables the user to calculate the corresponding ED50 for each specific functional GO term, important for risk assessment.

Key Words: microarray; GO analysis; dose and time dependent; quantitative.


    INTRODUCTION
 TOP
 ABSTRACT
 INTRODUCTION
 SYSTEMS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSION
 REFERENCES
 
Microarray technology has become a powerful tool to explore the expression levels of thousands of genes or even complete genomes after exposure to toxicants. A major challenge in the interpretation of toxicogenomic data is the functional interpretation, linking potentially interrelated alterations in gene expression to conventional toxicological end points (Hamadeh et al., 2002Go; Moggs, 2005Go). The functional interpretation of microarray data sets represents a laborious and challenging task. Researchers gather annotations from databases and identify characteristics of extracted sets of genes. The Gene Ontology (GO) Consortium initiated the standardization of annotation terms making them applicable for different organisms, facilitating data exchange among laboratories and databases (Ashburner et al., 2000Go; Camon et al., 2004Go; Harris et al., 2004Go). These characteristics render GO annotations a powerful tool for the interpretation of microarray data. The ability to determine which GO terms apply or which biological pathways are associated with each differentially expressed gene from a microarray experiment provides an ideal model to gain an understanding of what molecular processes are affected by the observed changes in gene expression. Recently, many groups have developed methods and tools for pathway analysis that are compatible with GO mapping and reveal statistically significant annotations associated with microarray data (Al-Shahrour et al., 2004Go; Beissbarth and Speed, 2004Go; Dennis et al., 2003Go; Doniger et al., 2003Go). The combination of using GO and pathway mapping has been proposed as a powerful approach to generate an unbiased view of biological processes and cellular components that are regulated by toxicant exposure at the transcriptional level (Currie et al., 2005Go). In general, these methods are based on a two-step process. First, a list of significantly altered genes is compiled based on statistical analysis (e.g., ANOVA). Second, annotation terms that are significantly over- or underrepresented within this list based on criteria such as fold change or statistical p value are compared to a reference list that usually consists of all genes on the array. An excellent example of the application of this approach for the analysis of di(2-ethylhexyl)phthalate exposure is given (Currie et al., 2005Go). However, these current methods only allow for a comparison of two experimental settings at a time and/or are restricted to one of the main GO categories or to a specific GO level. Furthermore, the resulting annotations are presented in extensive gene lists lacking quantitative information on each regulated gene. This one-dimensional analysis does not allow for quantitative or qualitative evaluation of possible dose- or time-dependent genomic relationships, which is crucial for the application of toxicogenomic data in the field of risk assessment.

An approach using the data average from many different genes categorized into broad "omic categories" has been proposed (Jansen and Gerstein, 2000Go). However, there is no consensus about the omic categories, and no systemic approach has been proposed to identify and calculate the average data. Since then, there has been a rapid development of the GO project, a collaborative effort to address the need for consistent descriptions of gene products in different databases (Ashburner et al., 2000Go; Camon et al., 2004Go; Harris et al., 2004Go). The GO collaborators are developing structured, controlled vocabularies (ontologies) that describe gene products in terms of their associated biological processes, cellular components, and molecular functions in a species-independent manner. The use of GO terms by several collaborating databases facilitates uniform queries across them. The Gene Ontology identifier (GOID) for the term seems to be a perfect candidate for the proposed omic categories. Computer programs such as MAPPFinder are available to link gene expression data to the GO hierarchy and to statistically identify the GO terms altered by toxicant exposure (Doniger et al., 2003Go). To facilitate quantitative interpretation of dose- or time-dependent genomic data, we propose a system-based approach to integrate the raw gene expression data with the GO Ontological analysis. We developed a computer program called GO-Quant that automatically links functional gene category analysis result from MAPPFinder with the original gene expression data and calculates the average intensity or ratio for those significant genes within each GO term. We applied this approach to demonstrate its application in a published dose- and time-dependent toxicogenomic data set (Dillman et al., 2005Go). Our results indicate that the above systematic approach can describe quantitatively the degree to which functional gene systems change across dose or time course. Additionally, this approach provides a global measurement to illustrate our results compared to single-gene assessments and enables the user to calculate the corresponding ED50 for each specific functional GO term.


    SYSTEMS AND METHODS
 TOP
 ABSTRACT
 INTRODUCTION
 SYSTEMS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSION
 REFERENCES
 
Gene expression data set.
In this paper, dose and temporal genomic data were obtained from Dillman et al. (2005)Go. The complete raw data set was retrieved from NCBI's gene expression omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) through GEO Series accession number GDS 1027. As described in the study of Dillman et al. (2005)Go, rats were injected in the femoral vein with liquid bis(2-chloroethyl)sulfide (sulfur mustard, SM), which circulates directly to the pulmonary vein and then to the lung. Rats were exposed to saline, isopropyl alcohol (vehicle control), 1, 3, or 6 mg/kg of SM in vehicle (isopropyl alcohol) and lungs harvested at 0.5, 1, 3, 6, and 24 h postinjection. Three biological replicates were used for each time point and dose tested. Physiological saline administrated in a similar manner was also included as an additional control. RNA was extracted from the lungs and used as the starting material for the probing of replicate oligonucleotide microarrays. All experiments were performed using Affymetrix Rat RAE230A oligonucleotide arrays (Affymetrix, Santa Clara, CA).

Statistical analysis and comparisons.
Raw gene expression data were imported to BRB Array Tools v3.3 for statistical analysis (Wright and Simon, 2003Go). A log2 transformation was applied to the data set. Each array was normalized by using the median intensity over the entire set of arrays. Since in this paper we are demonstrating our system approach by applying it to dose or time course microarray data, we conducted one-way ANOVAs across doses or time points. Statistical comparisons were conducted across doses (vehicle control, 1, 3, or 6 mg/kg) for each time point (0.5, 1, 3, 6, and 24 h) or across time points for each dose by using the "randomized variance" F test (Long et al., 2001Go; Wright and Simon, 2003Go). Significantly changed genes across doses (dose response) at specific time points (6 and 24 h) were selected based on one-way ANOVA test at p < 0.005 while significantly changed genes across time points (time course) at a dose of 6 mg/kg were selected based on one-way ANOVA test at p < 0.001 since there was a large number of genes changed with time. Hierarchical clustering analysis of functional gene category based on MAPPFinder results was performed across doses at two time points (6 and 24 h) or across time points at a dose of 6 mg/kg for the significantly changed genes. For the selected significant genes, ratios were derived by dividing each normalized gene value by the average value of the controls (vehicle) and then transforming the ratios to their respective log2 ratios. These log2-transformed ratios were used as input for the TIGR Multiexperiment Viewer (MEV) (Saeed et al., 2003Go). Hierarchical clustering analysis was conducted on the output genes using average linkage and Euclidean dissimilarity methods (Eisen et al., 1998Go). The complete output of the cluster analysis is shown with the fold change indicated colorimetrically.

GO analysis.
To establish associations between treatment and the affected GO terms and pathways, we used MAPPFinder (http://www.GenMapp.org) (Doniger et al., 2003Go), which links gene expression data to the GO hierarchy (Ashburner et al., 2000Go; Camon et al., 2004Go; Harris et al., 2004Go). MAPPFinder calculates a Z score as well as a permutation test p value for the genes that are significantly altered across doses (p ≤ 0.005) at 6 or 24 h or across time points (p ≤ 0.001) at a dose of 6 mg/kg doses. Since we proposed to use combined average raw gene expression values (e.g., intensity or ratio) of genes associated with specific functional gene categories derived from the MAPPFinder results and there are two directions of gene expression changes at a certain dose (or time), either upregulation or downregulation within a specific gene category, the average of these two different directions of gene expression alteration will mask the degree of changes. Therefore, the upregulated or downregulated genes will be calculated separately in the MAPPFinder. In this study we found simple patterns of expression of a consistent upregulation or consistent downregulation across dose and time based on K-means cluster analysis (Soukas et al., 2000Go). K-means cluster analysis was used to group individual genes across doses (or times) into two groups classified as up- or downexpression based on the trend with dose (or time) of the mean expression pattern. The number of clusters used in our analysis was chosen based upon having the mean expression patterns against dose (or time), and this assumption of using only two clusters was tested by examining using more than two clusters. To test the robustness of this method, we compared the K-means cluster analysis with a manual analysis of the trend. In the manual analysis of trend, we assigned an upregulation trend of a gene when the gene expression ratio was significantly elevated in at least three of the four groups or elevated at both SM 3 and 6 mg/kg groups. The trend of downregulation of a gene across groups was assigned when the gene expression level was significantly downregulated in at least three of the four groups or downregulated at both SM 3.0 and 6 mg/kg groups. The manual analysis of trend tended to create isolated clusters with few genes and larger clusters poorly represented by their mean expression patterns.

Each GO node and the cumulative total of the number of genes either upregulated or downregulated in a parent GO term combined with all its children were calculated. Z scores, a statistical measure of significance for gene expression in a given group, were calculated by subtracting the number of genes expected to be randomly changed in a GO term from the observed number of changed genes in that GO term. This value was then divided by the SD of the observed number of genes under a hypergeometric distribution. The equation can be written out as

Formula 1(1)
where N is the total number of genes measured, R is the total number of genes meeting the criterion that the gene be significantly changed based on an F test at significant p value, n is the total number of genes in each specific GO term, and r is the number of genes meeting the criterion in this specific GO term. A positive Z score indicates that there are more genes meeting the criterion in a GO term than would be expected by random chance. If the MAPPFinder data truly obeyed the assumptions of the hypergeometric distribution, then a Z score of 1.96 or –1.96 would correlate with a p value of 0.05. In order to address the multiple testing that occurs in GO database in using Z score, a nonparametric permutation test of the data was calculated (see Doniger et al., 2003Go, for more detail). This Z score allows us to rank and describe alterations in gene expression by each GO term based on the relative amount of altered gene expression within each category. For a GO term to be included in Tables 24, we set that at least three genes changed significantly (nested results) and that the permutation p value ≤ 0.05.


View this table:
[in this window]
[in a new window]
 
TABLE 2 Statistically Significantly Changed Gene Ontology Terms across Dose at 6 h

 

View this table:
[in this window]
[in a new window]
 
TABLE 4 Statistically Significantly Changed GO Terms across Time at a Dose of 6 mg/kg of SM

 
Quantitative integration of GO analysis to evaluate dose- or time-dependent genomic data.
MAPPFinder links gene expression data to the GO hierarchy and statistically identifies the GO term altered by the toxicant exposure (Doniger et al., 2003Go). Selecting genes for each category based on the arbitrary cutoff of fold change or p value change does not take into account the magnitude of the fold change of each gene, which is critical in the evaluation of dose- or time-dependent genomic data. To facilitate quantitative interpretation of dose- or time-dependent genomic data, we propose a system-based approach to integrate the raw gene expression data with the results from GO analysis. We developed a program called GO-Quant which automatically links results from functional GO analysis (MAPPFinder) with the original gene expression data and calculates the average intensity or ratio for those significantly genes within each GO term. For example, the average intensity or ratio of all significantly changed genes under statistically identified GOID was automatically calculated for each dose. Within a GOID, both upregulated and downregulated genes exist, and the calculation of this kind of average can be complicated. Thus, in our MAPPFinder analysis, the upregulated and downregulated genes were analyzed separately, and the average of intensity or ratio within a GOID in the upregulation or downregulation was calculated separately. The output file from the GO-Quant includes GO term, Z score, or average intensity or ratio, as well as gene ID. In the further analysis, the above output file was imported to MEV (Saeed et al., 2003Go) and a hierarchical clustering analysis on each GO term conducted using average linkage and Euclidean dissimilarity methods (Eisen et al., 1998Go). This hierarchical clustering analysis allows one to visualize the quantitative GO results. A stand-alone application for the above calculation (GO-Quant) is available now from the authors (http://depts.washington.edu/irarc/Go-Quant/).


    RESULTS AND DISCUSSION
 TOP
 ABSTRACT
 INTRODUCTION
 SYSTEMS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSION
 REFERENCES
 
Framework for the Quantitative Functional Interpretation of Dose- and Time-Dependent Genomic Data
As shown in Figure 1, typically three main steps are included in the framework for the quantitative functional interpretation of dose- and time-dependent genomic data. The first step is a traditional microarray analysis, aimed at identifying the statistically significant genes across dose or time course by using the ANOVA method. The next step is to identify significantly altered functional GO gene categories or pathways by using unsupervised GO analysis or a pathway-mapping tool such as MAPPFinder. The last step is to extract the quantitative gene expression data and calculate the average intensity or ratio for those genes significantly altered within each functional gene category based on MAPPFinder by using our GO-Quant program. In order to demonstrate its power in quantitative evaluation of dose- and time-dependent microarray data, we applied this approach to a published dose- and time-dependent toxicogenomic data set in rats treated with SM.


Figure 1
View larger version (31K):
[in this window]
[in a new window]
 
FIG. 1. System-based schematic framework for the interpretation of dose- and time-dependent genomic data.

 
Table 1 shows the outline of gene alteration across dose and time. Time course of the total number of significantly changed genes across doses (0, 1, 3, or 6 mg/kg) for each time point (0.5, 1, 3, 6, and 24 h) by one-way ANOVA comparisons (p ≤ 0.005) is listed at the right side of the table. The number of significantly changed genes increased in a time-dependent manner, peaking at 6 h. Dose-dependent increase in the number of significantly changed genes (p ≤ 0.005) across time points (0.5, 1, 3, 6, and 24 h) was listed in the left side of the table. Less than 137 genes changed with isopropyl alcohol across the time period observed, and more than 2600 genes changed with a dose of 6 mg/kg of SM.


View this table:
[in this window]
[in a new window]
 
TABLE 1 Dose- and Time-Dependent Gene Alterations after Injection of SM

 
Dose-Dependent Alteration of Gene Expression after Injection of SM
Hierarchical clustering analysis.
Figure 2 shows the hierarchical cluster analysis of gene expression alterations across dose. Significant gene expression changes were seen in a dose-dependent manner after injections of SM at 6 (A) and 24 h (B). Genes significantly changed across doses (p ≤ 0.005) at 6 (737 genes) or 24 h (441 genes) after injection were chosen to conduct the hierarchical clustering analysis. Examples of dose-dependent changes included downregulation of protein-tyrosine phosphatase receptor type G and upregulation of cyclin-dependent kinase inhibitor 1A and cyclin G1 transcription after 6 h postinjection (A and B clusters of Fig. 2A). Additionally, dose-dependent downregulation of T-cell receptor beta-chain and granzyme M and upregulation of cyclin-dependent kinase inhibitor 1A and cyclin G1 transcription were observed 24 h postinjection (A and B clusters in Fig. 2B). However, based solely on the above statistical analysis, it is difficult to interpret these example changes in gene expression within a biological context.


Figure 2
Figure 2
View larger version (142K):
[in this window]
[in a new window]
 
FIG. 2. Hierarchical cluster analysis of dose response gene expression alteration after injection of SM at 6 (A) and 24 h (B). Genes significantly changed across doses at 6 (717 genes) or 24 h (441 genes) after injection were included in hierarchical clustering analysis using average linkage and Euclidean dissimilarity methods. Clusters A and B show the details of genes including Affymetrix probe ID and gene name.

 
GO analysis of SM-induced genes alterations.
We used MAPPFinder to find significant changes within each GO functional category. MAPPFinder linked the 11,212 probe sets measured in this experiment to the GO database. A total of 4485 probes were linked to an Ensembl ID and 3144 genes linked to a GO term in the hierarchy. MAPPFinder calculates a Z score as well as a permutation test p value for the genes that are significantly altered across doses (p ≤ 0.005) at 6 or 24 h, or across time points (p ≤ 0.001) at a dose of 6 mg/kg. MAPPFinder calculated the percentage of genes meeting the criteria of either a statistically significant increase or decrease in gene expression across dose (or time). In addition, a Z score for each GO term was calculated based on Formula 1 to identify which GO terms had a significant number of altered genes. Since there are two directions of gene expression alteration after treatment at certain a gene category, the separation of the genes based on their trend (either upregulation or downregulation) of average gene expression pattern across dose (or time) allows to quantitatively evaluate the degree of changes across dose or time. In this study we identified two simple patterns of either consistent upregulation or consistent downregulation for gene expression across doses (or time) within GO terms based on K-means cluster analysis. MAPPFinder results represent a global picture of biological processes, cellular components, and molecular functions that are significantly altered (upregulated or downregulated) at the transcriptional level after treatment (Tables 2 and 3). Significant downregulations of genes involved in protein modification (GOID 6464), such as phosphorylation (GOID 6468) and signaling transduction (GOID 7165, 8277, 7242, 9966, and 9968), were observed based on biological process gene category 6 h after SM injection (Table 2). In contrast, significant upregulation of genes involved in the cell cycle (GOID 7049) was observed based on biological process gene category 6 h after SM injection (Table 2). As shown in Table 3, in the biological process, genes in cell cycle, especially in M phase, such as M phase of mitotic cell cycle (GOID 87), mitosis (GOID 7067), cell division (GOID 51301), M phase (GOID 279), and mitotic cell cycle (GOID 278), were upregulated. Transferase activity (GOID 16746) was also upregulated. Significant downregulation in immune responses (GOID 6955) such as cell activation (GOID 1775) and cell adhesion (GOID 7155) was observed.


View this table:
[in this window]
[in a new window]
 
TABLE 3 Statistically Significantly Changed GO Terms across Dose at 24 h

 
Consistent with the previously published analysis (Dillman et al., 2005Go), we identified alterations in the cell cycle regulation functional gene category as the most striking changes resulting from SM exposure. Furthermore, GO analysis suggests that multiple biological pathways are involved in SM-induced toxic effects. A recent report comparing global gene expression analyses between laboratories and across platforms suggested that GO analysis is a more biologically meaningful and statistically robust approach to data analysis than relying on single-gene analyses (Bammler et al., 2005Go). Our GO analysis demonstrates that unsupervised GO mapping and subsequent statistical significance analysis constitute a powerful unbiased approach to defining the biological pathways of known and novel SM-induced molecular changes in the rat lung. We noticed in MAPPFinder analysis that among the 11,212 probe sets measured in this experiment, only 4485 probe sets linked to an Ensembl ID and 3144 genes linked to a GO term. Thus, the biological significance of a large number of rat genes still needs to be defined. Once the annotation of these rat genes is completed, GO mapping and subsequent statistical significance analysis of microarray data will prove to be a powerful unbiased approach in defining biological pathways.

Quantitative Integration of GO Analysis to Evaluate Dose Response
Toxicogenomic data can provide profoundly important and novel information to the field of risk assessment. However, we are currently facing the difficult challenge of how to quantitatively interpret dose- or time-dependent relationships resulting from genomic approaches. As presented, the conventional GO analysis does not provide any quantitative information about microarray data in a dose-dependent manner, and no methodology is available to achieve quantitative analysis of pathway-specific information. Thus, we developed a framework to evaluate functional genes from microarray data, as illustrated in Figure 1, and propose the use of such data within a statistically significant functional gene category, such as a specific GO term in the GO database. The above calculation includes a large number of genes within each functional gene category, and there is currently no tool available to conduct such a task. We developed a computer program called GO-Quant, which automatically links functional GO analysis results from MAPPFinder with the original microarray data set, to calculate the average intensity or ratio for those significant genes within each GO term. Figure 3 shows the hierarchical cluster analysis of quantitative GO results across SM treatments after 6 (Fig. 3A) or 24 h (Fig. 3B). Compared to the conventional qualitative GO analysis shown in Tables 2 and 3, Figure 3 shows the dose-dependent alterations of each functional gene category (GOID) following SM treatment. Cell cycle genes (GOID 7049) were dose dependently upregulated 6 h after SM injection (Fig. 3A), and cell cycle alterations, especially in M phase pathway genes (GOID 87, 279, 278, 7067, and 51301), were dose dependently upregulated 24 h after SM injection (Fig. 3B). Dose-dependent downregulation of protein modification gene pathways, such as phosphorylation, were observed at the 6 h time point, while dose-dependent downregulation of immune response gene pathway (GOID 6952, 1775, and 45321) was observed at 24 h after SM injection (Fig. 3B). Figures 4 and 5 further demonstrate the dose response of quantitative gene expression pathway alterations within each functional gene category in biological process (Fig. 4) and molecular function genes (Fig. 5). Significant upregulation in cell cycle gene pathways becomes obvious at 3 mg/kg both at 6 (Fig. 4A) and at 24 h (Fig. 5A) after injection. Our results confirm that the above system approach is a more robust and clear measurement to determine the degree to which a functional gene system changes across dose. As an alternative to representing annotations by the average of the annotated genes, the sum or the median expression profile could also be calculated using the GO-Quant. In most cases, the resulting plots are virtually identical (data not shown). Such an approach gives one the ability to calculate the ED50 for a specific functional GO term, which is important for risk assessment.


Figure 3
Figure 3
View larger version (107K):
[in this window]
[in a new window]
 
FIG. 3. GO-based hierarchical cluster analyses of gene expression data across dose after injection of SM at 6 (A) and 24 h (B). GO analysis was performed using MAPPFinder. The average ratio of the genes in a statistically significant functional gene category (GOID in GO database) was calculated with GO-Quant. The output file, including GO term, Z score, average intensity or ratio, and gene ID, was imported to TIGR MEV (Saeed et al., 2003Go) and a hierarchical cluster analysis conducted using average linkage and Euclidean dissimilarity methods (Eisen et al., 1998Go).

 

Figure 4
View larger version (30K):
[in this window]
[in a new window]
 
FIG. 4. Dose response of GOID in biological process (A) and molecular function gene categories (B) 6 h after SM injection. GO analysis was performed using MAPPFinder. The average ratio of the genes in a statistically significant functional gene category (GOID in GO database) was calculated with GO-Quant.

 

Figure 5
View larger version (24K):
[in this window]
[in a new window]
 
FIG. 5. Dose response of GOID in biological process (A) and molecular function gene categories (B) 24 h after SM injection. GO analysis was performed using MAPPFinder. The average ratio of the genes in a statistically significant functional gene category (GOID in GO database) was calculated with GO-Quant.

 
Time Course Gene Alteration after a Single Injection of SM at a Dose of 6 mg/kg
In order to demonstrate the above application to the time course within this microarray data set, we further applied the above approach as illustrated in Figure 1 to analyze the time course data after the single injection of SM at a dose of 6 mg/kg. Figure 6 shows the hierarchical clustering analysis of 1111 genes whose transcription was significantly (p ≤ 0.001) altered across the selected time points (6 mg/kg SM). Using cluster analysis, we identified changes in transcription in specific genes, e.g., upregulation of heme oxygenase 1, cyclin G, and cyclin-dependent kinase inhibitor 1A and downregulation of protein kinase C and chemokine ligand 12. As discussed, it is challenging to interpret these previous results in a biological context. Table 4 shows the list of GOID from the MAPPFinder analysis representing a global picture of biological processes, cellular components and molecular function genes that are significantly altered (upregulated or downregulated) at the transcriptional level following SM exposure at 6 mg/kg (Table 4). As shown in Figures 7 and 8, GO-Quant analysis revealed time-dependent alterations in various functional gene categories within biological process, cellular component, and molecular function genes. In the biological process gene category, time-dependent downregulation in immune cell activation (GOID 45321, 46649, and 1775), cell division (GOID 51310), and response to drug (GOID 42493) were observed, especially at 6 h after SM exposure. In addition, time-dependent upregulations in the regulation of cell cycle (GOID 7050 and 74), protein ubiquitination (GOID 16567) gene pathways were also observed, particularly at 6 h postexposure. Our results confirm that the above systems approach allows us to demonstrate time-dependent alteration in cellular and molecular gene function based on GO categorization. It is not our focus in this paper to describe why these functional transcriptional alterations after SM injection are observed. Our current analysis is intended to demonstrate that our systematic approach to applying the average data within the statistically significant GOID is a powerful tool to quantitatively describe and interpret dose- or time-dependent relationships discovered using microarray data.


Figure 6
View larger version (92K):
[in this window]
[in a new window]
 
FIG. 6. Hierarchical cluster analysis of altered gene expression in a time-dependent manner after injection of SM at a dose of 6 mg/kg. Genes significantly changed across time points after exposure were included in hierarchical cluster analysis using average linkage and Euclidean dissimilarity methods. Clusters A and B show the details of genes including Affymetrix probe ID and gene name.

 

Figure 7
View larger version (64K):
[in this window]
[in a new window]
 
FIG. 7. GO-based hierarchical cluster analyses of gene expression data across time after injection of SM at a dose of 6 mg/kg. GO analysis was performed using MAPPFinder. The average ratio of the genes in a statistically significant functional gene category (GOID in GO database) was calculated with GO-Quant. The output file, including GO term, Z score, average intensity or ratio, and gene ID, was imported to TIGR MEV (Saeed et al., 2003Go) and a hierarchical cluster analysis conducted using average linkage and Euclidean dissimilarity methods (Eisen et al., 1998Go).

 

Figure 8
View larger version (34K):
[in this window]
[in a new window]
 
FIG. 8. Time-dependent alteration of GOID in biological process (A) and molecular function genes (B) after SM injection at dose of 6 mg/kg. GO analysis was performed using MAPPFinder. The average ratio of the genes in a statistically significant functional gene category (GOID in GO database) was calculated with GO-Quant.

 

    CONCLUSION
 TOP
 ABSTRACT
 INTRODUCTION
 SYSTEMS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSION
 REFERENCES
 
Using gene expression profiling following SM exposure, we have conducted a conventional GO analysis and compared it with our proposed quantitative approach for functional GO analysis. Our results confirmed that GO analysis is a powerful approach to generate an unbiased view of the functional gene alterations by SM. However, the quantitative approach we propose using combined average raw gene expression values (e.g., intensity or ratio) of genes associated with specific functional gene categories derived from the MAPPFinder results in conjunction with the GO-Quant algorithm provides a powerful tool to interpret microarray data in a dose- or time-dependent manner. Additionally, this approach provides a global measurement based on pathway response as compared to single-gene assessments and enables one to calculate the corresponding ED50 for each specific functional GO term, which is important for risk assessment. The application of such quantitative interpretation of toxicogenomic data is likely to become increasingly useful for the interpretation of toxicogenomic data generated for evaluating mechanistic similarity of novel chemicals.


    ACKNOWLEDGMENTS
 
This work was supported in part by a toxicogenomic grant from National Institute of Environmental Health Science (NIEHS) (U10 ES 11387), the US Environmental Protection Agency-NIEHS UW Center for Child Environmental Health Risks Research (EPA R826886 and NIEHS 1PO1ES09601), the NIEHS grant R01-ES10613, and Center for Oceans and Human Health Research (NIEHS: P50 ES012762 and NSF (National Science Foundation): OCE-0434087), and UW NIEHS Center for Ecogenetics and Environmental Health (5 P30 ES07033).


    REFERENCES
 TOP
 ABSTRACT
 INTRODUCTION
 SYSTEMS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSION
 REFERENCES
 
Al-Shahrour, F., Diaz-Uriarte, R., and Dopazo, J. (2004). FatiGO: a web tool for finding significant associations of Gene Ontology terms with groups of genes. Bioinformatics 20, 578–580.[Abstract/Free Full Text]

Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., Davis, A. P., Dolinski, K., Dwight, S. S., Eppig, J. T., et al. (2000). Gene ontology: Tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet. 25, 25–29.[CrossRef][ISI][Medline]

Bammler, T., Beyer, R. P., Bhattacharya, S., Boorman, G. A., Boyles, A., Bradford, B. U., Bumgarner, R. E., Bushel, P. R., Chaturvedi, K., Choi, D., et al. (2005). Standardizing global gene expression analysis between laboratories and across platforms. Nat. Methods 2, 351–356.[CrossRef][ISI][Medline]

Beissbarth, T., and Speed, T. P. (2004). GOstat: Find statistically overrepresented gene ontologies within a group of genes. Bioinformatics 20, 1464–1465.[Abstract/Free Full Text]

Camon, E., Magrane, M., Barrell, D., Lee, V., Dimmer, E., Maslen, J., Binns, D., Harte, N., Lopez, R., and Apweiler, R. (2004). The gene ontology annotation (GOA) database: Sharing knowledge in Uniprot with gene ontology. Nucleic Acids Res. 32, D262–D266.[Abstract/Free Full Text]

Currie, R. A., Bombail, V., Oliver, J. D., Moore, D. J., Lim, F. L., Gwilliam, V., Kimber, I., Chipman, K., Moggs, J. G., and Orphanides, G. (2005). Gene ontology mapping as an unbiased method for identifying molecular pathways and processes affected by toxicant exposure: Application to acute effects caused by the rodent non-genotoxic carcinogen diethylhexylphthalate. Toxicol. Sci. 86, 453–469.[Abstract/Free Full Text]

Dennis, G., Jr, Sherman, B. T., Hosack, D. A., Yang, J., Gao, W., Lane, H. C., and Lempicki, R. A. (2003). DAVID: Database for annotation, visualization, and integrated discovery. Genome Biol. 4, 3.

Dillman, J. F., III, Phillips, C. S., Dorsch, L. M., Croxton, M. D., Hege, A. I., Sylvester, A. J., Moran, T. S., and Sciuto, A. M. (2005). Genomic analysis of rodent pulmonary tissue following bis-(2-chloroethyl) sulfide exposure. Chem. Res. Toxicol. 18, 28–34.[CrossRef][ISI][Medline]

Doniger, S. W., Salomonis, N., Dahlquist, K. D., Vranizan, K., Lawlor, S. C., and Conklin, B. R. (2003). MAPPFinder: Using gene ontology and GenMAPP to create a global gene-expression profile from microarray data. Genome Biol. 4, R7.[CrossRef][Medline]

Eisen, M. B., Spellman, P. T., Brown, P. O., and Botstein, D. (1998). Cluster analysis and display of genome-wide expression patterns. Proc. Natl. Acad. Sci. U.S.A. 95, 14863–14868.[Abstract/Free Full Text]

Hamadeh, H. K., Bushel, P. R., Jayadev, S., Martin, K., DiSorbo, O., Sieber, S., Bennett, L., Tennant, R., Stoll, R., Barrett, J. C., et al. (2002). Gene expression analysis reveals chemical-specific profiles. Toxicol. Sci. 67, 219–231.[Abstract/Free Full Text]

Harris, M. A., Clark, J., Ireland, A., Lomax, J., Ashburner, M., Foulger, R., Eilbeck, K., Lewis, S., Marshall, B., Mungall, C., et al. (2004). The gene ontology (GO) database and informatics resource. Nucleic Acids Res. 32, D258–D261.[Abstract/Free Full Text]

Jansen, R., and Gerstein, M. (2000). Analysis of the yeast transcriptome with structural and functional categories: Characterizing highly expressed proteins. Nucleic Acids Res. 28, 1481–1488.[Abstract/Free Full Text]

Long, A. D., Mangalam, H. J., Chan, B. Y., Tolleri, L., Hatfield, G. W., and Baldi, P. (2001). Improved statistical inference from DNA microarray data using analysis of variance and a Bayesian statistical framework. Analysis of global gene expression in Escherichia coli K12. J. Biol. Chem. 276, 19937–19944.[Abstract/Free Full Text]

Moggs, J. G. (2005). Molecular responses to xenoestrogens: Mechanistic insights from toxicogenomics. Toxicology 213, 177–193.[CrossRef][ISI][Medline]

Saeed, A. I., Sharov, V., White, J., Li, J., Liang, W., Bhagabati, N., Braisted, J., Klapa, M., Currier, T., Thiagarajan, M., et al. (2003). TM4: A free, open-source system for microarray data management and analysis. Biotechniques 34, 374–378.[ISI][Medline]

Soukas, A., Cohen, P., Socci, N. D., and Friedman, J. M. (2000). Leptin-specific patterns of gene expression in white adipose tissue. Genes Dev. 14, 963–980.[Abstract/Free Full Text]

Wright, G. W., and Simon, R. M. (2003). A random variance model for detection of differential gene expression in small microarray experiments. Bioinformatics 19, 2448–2455.[Abstract/Free Full Text]


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


This article has been cited by other articles:


Home page
Toxicol SciHome page
R. S. Thomas, B. C. Allen, A. Nong, L. Yang, E. Bermudez, H. J. Clewell III, and M. E. Andersen
A Method to Integrate Benchmark Dose Estimates with Genomic Data to Assess the Functional Effects of Chemical Exposure
Toxicol. Sci., July 1, 2007; 98(1): 240 - 248.
[Abstract] [Full Text] [PDF]


This Article
Right arrow Abstract Freely available
Right arrow FREE Full Text (PDF) Freely available
Right arrow All Versions of this Article:
92/2/560    most recent
kfj184v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Add to My Personal Archive
Right arrow Download to citation manager
Right arrowRequest Permissions
Right arrow Disclaimer
Google Scholar
Right arrow Articles by Yu, X.
Right arrow Articles by Faustman, E. M.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Yu, X.
Right arrow Articles by Faustman, E. M.
Social Bookmarking
 Add to CiteULike   Add to Connotea   Add to Del.icio.us  
What's this?