Introduction

Top Abstract Introduction Materials and Methods Results and Discussion References

Toxicologists employ a battery of tests to identify chemicals with potential for human toxicity or that might cause environmental harm. According to the United States National Toxicology Program (NTP), a thorough analysis of each chemical requires $2 to 4 million and several years to complete (National Toxicology Program, 1996). Because of the cost- and labor-intensive nature of these studies, the number of chemicals currently tested by the NTP stands at 505 in long-term studies, 66 in short-term tests, and one subchronic study (http://ntp-server.niehs.nih.gov/). Given that there are approximately 70,000 chemicals in commerce today (National Toxicology Program, 1996), it is clearly impossible to apply current testing methods to all chemicals of concern. It is apparent that alternative testing approaches must be developed if science is going to maintain a significant role in environmental and public health policy.

The development of a screen that would allow prioritization of untested chemicals based upon their toxic potential would have a significant impact on how efficiently we evaluate both synthetic and naturally occurring compounds. One approach for predicting toxic potential is to classify chemicals based upon their capacity to alter transcriptional programs in a manner that is similar to known toxicants (Nuwaysir et al., 1999). Test chemicals that induce transcriptional responses in a manner similar to those induced by a known poison could then be classified as harboring toxic potential and examined carefully by more thorough toxicological means. This approach has two underlying assumptions: 1) that we have enough scientific information to allow proper classification of prototype toxicants and 2) that most if not all toxic chemical exposures will alter gene expression at some level. In support of this second assumption, signal transduction pathways that culminate in a transcriptional response mediate the toxicity of many chemicals. In addition, toxicity is commonly manifested as inflammation, proliferation, apoptosis, necrosis, and/or cellular differentiation. All of these toxic endpoints are intimately linked to specific alterations in gene expression.

To test the hypothesis that toxicants can be classified according to their influence on global gene expression profiles, we employed cDNA microarray technology (Schena et al., 1995; Golub et al., 1999) and attempted to classify 24 prototype chemical treatments that fall into five well-characterized toxicological classes. Three of the classes include environmental pollutants that are targets of regulatory concern [i.e., peroxisome proliferators, AHR agonists, and noncoplanar PCBs (Schmidt and Bradfield, 1996; Carpenter, 1998; Vanden Heuvel, 1999)]. The remaining two classes consist of agents that stimulate other common toxic endpoints, such as treatments that induce the inflammatory response and treatments that stimulate the hypoxia signal transduction pathway.

	Materials and Methods

Top Abstract Introduction Materials and Methods Results and Discussion References

Animals and Treatment. All animals treated were male C57BL/6J mice. The treatment, dose, vehicle, and time at sacrifice were selected from the literature and are shown in Table 1. All treatments were compared with their respective vehicle controls at the corresponding time points. A summary of these agents and their general toxicological category is provided in Table 2.

RNA Isolation. Total RNA from cells and frozen liver tissue was isolated using the RNeasy system (QIAGEN, Valencia, CA). Poly-A selected RNA was purified using the Oligotex kit (QIAGEN) and checked using gel electrophoresis and/or spectrophotometric measurements. mRNA from the livers of three or more mice were pooled for analysis on the microarrays.

cDNA Microarray Construction and Analysis. Approximately 1200 minimally redundant cDNAs from internal expressed sequence tag projects (http://edge.oncology.wisc.edu/) and the public expressed sequence tag effort were identified and the glycerol stocks were rearrayed into a separate set of clones. An aliquot from these clone sets was amplified by polymerase chain reaction and used to construct custom cDNA microarrays using methods described previously (Worley et al., 2000). Each cDNA clone was spotted four times on each slide for replicate analysis. Labeled cDNA probe was produced from 1 µg of poly-A RNA by incorporation of Cy3- or Cy5-dCTP (Amersham Pharmacia Biotech, Piscataway, NJ) during a standard reverse transcriptase reaction as described previously (Penn et al., 2000). The slides were scanned using a microarray scanner (Molecular Dynamics, Sunnyvale, CA) and the fluorescence data were analyzed using ArrayVision software package (Imaging Reseach, St. Catharines, ON, Canada). For each hybridization, the four spots corresponding to each cDNA were averaged and normalized to an internal suite of 15 housekeeping transcripts coding for ribosomal proteins. Normalization was carried out using methods described previously (Chen et al., 1997).

Data Reduction. Before statistical analysis, the data set was screened for transcripts that did not respond to any of the treatments used in the study and would not contribute significantly to the classification. A threshold of 2-fold change in gene expression was used as the cut-off value. The 2-fold cut-off value was slightly more conservative than the confidence limit calculated using techniques published previously (i.e., 1.4-fold) (Chen et al., 1997) and is similar to a standard threshold level used in other studies (e.g., Schena et al., 1996). Using this threshold value, we determined that approximately 500 of the 1200 transcripts changed significantly in response to at least one treatment.

Statistical Classification Analysis. The Naïve Bayesian classification used in this article is based on previous work by Kontkanen et al. (1998). Briefly, the predictor variables X₁,... ,X_k are assumed to be independent of each other when conditioned on the class variable C. Our model M is constructed by the joint probability distribution for a data vector (x, c) = (X₁ = x₁, ... , X_k = x_k, C = c) and can be written as follows:

(1)

Before incorporating any data, we also assume that all classes are equally probable (i.e., the probability that a chemical will belong to a certain class is the same for each class) and that within each class, the gene expression values of each transcript are also equally probable. Given these assumptions, we can use Bayesian probability theory to calculate the conditional predictive distribution for the class c given x and the data set D by:

(2)

where the numerator is calculated as:

(3)

where t_c is the number of treatments in class c, N_t is the total number of treatments, NC is the total number of classes, f_cxi is the number of cases in class c having a value equal to x_i, F_c is the number of treatments in class c, and V_xi is the number of possible values of x_i. The denominator is the same for all c and is calculated as:

(4)

The result of this conditional predictive distribution is then used to classify the data vector. For the parameter selection, a modification of forward parameter selection process outlined in Huberty (1994) was employed. Specifically, an iterative process was used in which transcripts were run individually using the Naïve Bayes model and the transcript with the best internal classification rate and highest confidence (represented by the sum of all probabilities for correctly classified treatments) was selected. In the subsequent round, the selected transcript was fixed and the remaining parameters were added individually to find which transcript, along with the first selected transcript, produced the highest internal classification rate and confidence. This process was repeated until all 74 transcripts were added to the model in the order of their internal classification rate. This type of selection ranks the transcripts in the order of their estimated predictive value and adds them sequentially to the model. It should be noted that the forward selection approach does not necessarily yield the best set or even the smallest set. In addition, genes with similar expression profiles may also be added to the set if they significantly increase the accuracy or confidence. The approach is simply a heuristic to look for a diagnostic set of predictors with a high accuracy. A flow chart describing the data reduction and classification analysis is provided in Fig. 1.

View larger version (22K): [in this window] [in a new window]   Fig. 1.   A flow-chart illustrating the method for the data reduction and classification analysis leading to the diagnostic set of transcripts.

	Results and Discussion

Top Abstract Introduction Materials and Methods Results and Discussion References

To allow the accurate classification of large numbers of toxicants based on gene expression, several important factors were considered in the experimental design. First, exposures were performed in inbred mice in an effort to minimize the influence of genetic polymorphism on transcriptional responses to toxicants. Second, gene expression monitoring was focused on the liver, because this organ is often the first significant site of chemical exposure and exhibits a wide array of pathological and adaptive responses to a broad spectrum of toxicants. Third, gene expression was characterized using mRNA obtained from whole liver after an acute exposure to a test chemical because the whole-organ response better represents the full range of transcriptional changes that can result from toxicity. In contrast to cell culture studies, this in vivo approach is particularly important when multiple cell-types and paracrine signaling are required for the complete toxic response (e.g., inflammation). Finally, in an effort to understand and adjust for the influence of temporal and sample acquisition variables, multiple time-points were analyzed for several of the treatments.

To represent the transcriptional response as a whole, a two-dimensional hierarchical clustering method was employed (Eisen et al., 1998) and the relationships between treatments based on gene expression profiles are highlighted in the adjacent dendrogram (Fig. 2). A visual inspection of the treatment-dendrogram indicates that individual chemicals, with a few exceptions, generally fall into their anticipated toxicological classes. Specifically, the hypoxia treatments (i.e., phenylhydrazine and cobalt) are intermixed with the noncoplanar PCBs showing some similarity among the gene expression profiles for these treatments (Fig. 2). Application of a more formal nearest-neighbor classification analysis using a similar correlation based metric indicated that 7 of 24 treatments were closer to treatments outside of their class as opposed to within (data not shown).

View larger version (59K): [in this window] [in a new window]   Fig. 2.   Variation in expression of approximately 500 transcripts in 24 experimental treatments. The data are presented after two-dimensional hierarchical clustering, which organizes transcripts and treatments on the basis of similarity. Each row represents a single transcript and each column an experimental treatment. For each treatment, the ratio of the expression of each transcript to the expression of the transcript in control samples is represented by the color of the corresponding cell in the graph. Green represents down-regulation of the transcript, black means no change, and red represents up-regulation of the transcript. Color saturation reflects the magnitude of the change. Clustering was performed using complete linkage clustering with a centered correlation coefficient as the similarity metric (Eisen et al., 1998). A dendrogram representing similarities between experimental treatments is shown at the top and the approximate locations of the genes in the final diagnostic set are noted with bars to the right of the figure. Wyeth, Wy-16,463.

The imperfect classification scheme attained by the initial clustering analysis on the transcriptional response across the whole microarray was not surprising. Although our classification scheme assumes that chemicals in each toxicant class act through a common mechanism, individual members of each class may stimulate unique pathways that may be secondary or unrelated to toxicity. For example, although considerable pharmacological and genetic evidence indicates that halogenated-dioxins lead to toxicity through their binding to the Ah-receptor, the agonist BNF is also known to stimulate the antioxidant response pathway (Poland and Glover, 1980; Radjendirane and Jaiswal, 1999). Similarly, although agents that stimulate the transcriptional response to hypoxia act through up-regulation of HIF1, some of these agents have also been shown to induce the acute phase response (Wenger et al., 1995). The resulting pattern of expression across the whole microarray reflects small, chemical-specific differences at the molecular level that result in a virtual `fingerprint' of expression for individual compounds. Therefore, that subset of transcripts that allows classification of these treatments must be defined according to our understanding of the primary toxic mechanism. In other words, predictor variables (in our case transcripts) must be screened for their ability to discriminate between groups, and a subset of these variables must be derived that allows accurate predictions.

As a method of identifying a diagnostic set of predictor transcripts, we applied a probabilistic approach based upon Bayesian statistics (Duda and Hart, 1973; Kontkanen et al., 1998). Here, gene expression values were discretized and a standard forward parameter selection algorithm was employed to select predictor variables to be added to the model (Huberty, 1994). This type of selection ranks the transcripts in the order of their estimated predictive value and adds them sequentially to the model. For example, in the first round of selection, all transcripts were run individually using the Naïve Bayes model, and the transcript with the best internal classification rate and highest confidence (represented by the probabilities for all correctly classified treatments) was selected. The internal classification rate for a given set of predictor variables is measured based on the classification rate within the training set. In the next round, the selected transcript was fixed and the remaining parameters were added individually to find which transcript, along with the first selected transcript, produced the highest internal classification rate and confidence. This process was repeated until all transcripts were added to the model in the order of their internal classification rate. To estimate the predictive accuracy of this approach, the process was integrated with leave-one-out cross-validation, in which one of the treatments is removed from the analysis, then the model is constructed and used to predict the left-out treatment.

The results of this analysis show the predictive accuracy as a function of the number of transcripts added to the model during the forward parameter selection process (Fig. 3). Based on this analysis, we found that the predictive accuracy and confidence of the model began to level out after the addition of approximately a dozen transcripts and even began to drop soon thereafter. Consequently, a set of 12 transcripts was considered `diagnostic' for the classification of treatments into the toxicological classes investigated in this study. The final `diagnostic set' was derived by following the same procedure on the complete data set (i.e., no treatment left out). The transcriptional profile of the 12 diagnostic transcripts, and their order of their addition to the model, is shown in Fig. 4.

View larger version (14K): [in this window] [in a new window]   Fig. 3.   Estimated predictive accuracy and relative confidence of the classification model with the addition of selected transcripts. The order of transcripts added to the model and the predictive accuracy were performed using a forward selection scheme and leave-one-out cross-validation, respectively (see text for details).

View larger version (48K): [in this window] [in a new window]   Fig. 4.   Variation in expression of the diagnostic set of transcripts used to classify the 24 treatments into the five toxicological classes. Each row represents a single transcript and each column an experimental treatment. The order of the transcripts, from top to bottom, is in the order of addition to the probability model. For each treatment, the expression ratio of each transcript in a treated sample to the expression in the control sample is represented by the color of the corresponding cell in the graph. Green represents down-regulation of the transcript, black means no change, and red represents up-regulation of the transcript. Color saturation reflects the magnitude of the change and can be compared with the ratio key. Wyeth, Wy-16,463.

The forward selection and cross-validation process identified several properties in our toxicological profiles. First, a `diagnostic set' of 12 transcripts was identified that allows an estimated 100% predictive accuracy for the toxicological classes chosen for this study. Interestingly, the transcripts in this diagnostic set are a mix of transcripts previously known to be altered by these treatments and relatively uncharacterized transcripts with respect to toxicant regulation. For example, CYP1A2 and CYP4A10 are known to be up-regulated by TCDD and peroxisome proliferators, respectively (Bell et al., 1993; Schmidt and Bradfield, 1996). In contrast, for IL-18 and betaine homocysteine methyltransferase, we have little information regarding their response after toxicant exposure. Overall, about half of the changes in the optimal set were described previously at some level in the literature. Second, the forward selection approach also explains the previously described 50 to 70% predictive accuracy attained when using the whole data set (500 transcripts), in that soon after the addition of the diagnostic set transcripts, the estimated predictive accuracy begins to decline significantly. In other words, the further addition of transcripts beyond the diagnostic set begins to split out treatments and the `individuality' of the treatment profiles begins to take over.

Despite the apparent success of our study in establishing the potential for classification of toxic chemicals according to diagnostic sets of genes, the success should be framed in a couple ways. It should be noted that the exposures chosen in this study represent single acute doses and may not reflect the gene expression changes at lower, environmentally relevant doses. In an organism, multiple factors converge to ultimately influence the manifestation of toxicity and the associated gene expression patterns. Among these factors are time, dose, route of administration, age of the animal, and sex. Characterizing the influence of all of these variables on transcript profiles with even a small number of treatments would require considerable resources (e.g., 12 treatments × 3 time points per treatment × 3 doses per time point × 3 routes per dose = 324 microarray studies). Although in this study we chose to primarily address the factor of time, additional experimentation was also performed to look at how different doses of TCDD affected classification. In these preliminary dose-response studies, our statistical model proved to be resistant to the variation in TCDD doses with correct classification at doses as low as 0.05 µg/kg and as high as 100 µg/kg (data not shown). Arguably, other chemicals may not be as easy to classify over such a large dose range and additional studies will be needed to address this issue and other factors that may affect the predictive accuracy of the model. To understand how our statistical model performed with a treatment that was not in our five toxicological categories, results from arsenic-treated mice were also analyzed. Interestingly, the model did not classify the arsenic treatment with any degree of confidence, with inflammatory as the closest category and AHR agonists as the least similar category (data not shown). This adds additional confidence that accurately predicting toxicological endpoints based on gene expression is achievable.

Several interesting conclusions can be drawn from this work. Most importantly, we have presented evidence that the accurate classification of toxic chemicals according to their transcript expression profiles is possible. This opens the door to a new era of toxicological testing where relatively short and inexpensive studies using transcript expression as an endpoint allow the prioritization of untested chemicals based upon their classification. This would mean significant savings in both animal usage and financial resources and would reduce the disparity between the number of tested and untested chemicals in commerce today. However, this study is just the first step toward this goal. The toxicological categories selected in our study primarily reflect the model compounds that toxicologists have studied extensively over the last decade and represent only a small percentage of the 70,000 chemicals in commerce today.

Obviously, as the public gene expression database grows, more toxicological categories can be added to the model and the more predictive our model and those like it will become. Another conclusion from this work is that large arrays with thousands of transcripts are unnecessary to make these classifications. Although the large arrays are necessary to initially identify the diagnostic gene set, once these `diagnostic' sets of indicator transcripts are identified, simple measurements of only one or two dozen transcripts may allow the average investigator the ability to make judgments as to the relative toxicity of a particular chemical. Finally, we have purposely chosen to use a robust set of statistical methods and conservative assumptions to develop this predictive set. The discrete nature of the approach does not require an accurate measurement of transcriptional changes beyond the assessment of only significant up- or down-regulation. This should make subsequent evaluations more resistant to inter- and intralaboratory variability such as that observed when switching arrays or performing hybridizations in multiple laboratories.