The chemical treatments used in the experiments on female B6C3F₁ mice are summarized in Table 1. Note the treatments are the same as the maximum tolerable dose (MTD) used for the corresponding chemical in the two-year bioassay. Treatments involving lower doses of a subset of 5 chemicals from Table 1 are provided in Table 2. More detailed descriptions of the chemicals, treatments and gene expression analysis are provided in the Text S1.

To address the question of species extrapolation, gene expression and tumor data from chemical exposures in Fischer 344 rats were obtained from the literature [5]. The chemical treatments are summarized in Table S1. Additional information is provided in the Text S1. Normalized probe intensity data for each treatment were grouped separately with the control and normalized using the quantile normalization function in MATLAB (2008a, The MathWorks, Natick, MA).

The pathology review of the treated mice after 90 days of exposure to hepatocarcinogens did not identify liver tumors (see Text S1). Hence, the 90 day gene expression patterns are potentially reflective of a pre-neoplastic state in the case of the chemicals that are known hepatocarcinogens (as characterized by the 2 year bioassay). We considered known risk factors for various human diseases including liver cancer. Both genetic (i.e., single nucleotide polymorphisms) and non-genetic (i.e., disease) risk factors were considered.

Genetic risk factors in terms of single nucleotide polymorphisms associated with various human diseases including liver cancer are tabulated on databases like the Genetic Association Database [20]. Genes having polymorphisms associated with various human cancers,Alzheimer's disease and Schizophrenia were downloaded from the database. Identifying non-genetic risk factors for liver cancer has been an active research topic [21], [22], [23], [24], [25], [26], [27], [28]. Among the main factors contributing to the risk are Hepatitis C and B virus [22], [24], [25], cirrhosis [21] and diabetes [27], [28]. The risks of liver cancer associated with cirrhosis induced by Hepatitus C virus was estimated to be comparable with cirrhosis induced by non-alcoholic steatohepatitis (NASH) [29]. In fact NASH is correlated with characteristics of metabolic syndrome like obesity and diabetes. Various human gene expression datasets associated with these risk factors are identified on the Gene Expression Omnibus (GEO) database [30]. The processed and normalized gene expression data were used as is from the database except for cases where the data were presented as intensity values. In these cases, the intensity values were log₂-transformed. The entire set of data sets associated with risk factors for liver cancer is summarized in Table 3.

The biochemical pathways used in the analysis in this paper were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway database [17], [18], [19]. The data for the set of genes involved in each of the pathways and their associated interactions were downloaded via the KEGG API. All pathways in the database were used. The list of 216 pathways that were homologous across human, mice and rats species, along with their mouse KEGG ids are reported in Table S2.

SEPEA, a network based pathway enrichment method described in detail in [31], was used to evaluate the linkage between the gene expression data and the KEGG pathways. Unlike traditional pathway enrichment methods that treat pathways as sets of genes, SEPEA treats pathways as networks of interacting proteins and/or enzymes. The genes corresponding to the proteins in the signaling network are given more weight according to whether they are at the receptor or the terminating end of the pathway that typically signals for transcription in a number of genes. Further, pathways where the perturbed genes are close relative to each other on the associated network are modeled as being more likely to be affected than pathways where the perturbed genes occur further apart over the network. There were three analytic methods described in [31]; in the work described here, SEPEA_NT2 is used [31] for the gene expression analysis and SEPEA_NT3 for the gene polymorphism data. The goal of the enrichment analysis is to assign significance (in terms of p-values) to all chosen KEGG [17], [18], [19] human, rat and mouse pathways for the increased likelihood of being affected in the livers of the treated animals over the matched controls. The significance obtained by SEPEA_NT2 were based on 5×10⁴ randomizations (see [31] for details). For the analysis performed here, the SEPEA p-value was converted to a z-score, Z_ij using the equation,(1)where Φ refers to the standard normal distribution and P_ij denotes the p-value obtained using SEPEA for pathway j as a result of the mice being dosed with chemical i at a given dose or human data set i.

The carcinogenicity rate of a chemical is defined using the survival-adjusted proportions of animals treated with the chemical that developed liver adenomas or carcinomas and the corresponding proportion for the control animals at the end of the two-year bioassay. The data for all chemicals except tetrachloroethane [32] were obtained from their respective technical reports for the two-year cancer bioassay developed by the NTP(see Text S1 for reference to the NTP technical reports). The poly-3 survival-adjusted numbers [15], [33] were used when available, else the survival adjustment provided in the technical reports were used. The poly-3 statistic, z_i, for chemical i used in this analysis is defined as,(2)where n_i^treat and p_i^treat are the survival adjusted number and proportion, respectively of animals treated with chemical i that developed adenomas or carcinomas, n_i^cont and p_i^cont are the corresponding terms for the control animals. To avoid outlier effects in the analysis, z_i = z_2e−5 was used when z_i≥z_2e-5 where z_2e−5 is the critical value of the standard normal distribution corresponding to 2×10⁻⁵ significance.

The critical hypothesis of this manuscript is that the observed changes in gene expression manifest at the pathway level have adequate information to predict future carcinogenicity. However, it is not obvious what mathematical functional form relates the observed changes in pathways to the carcinogenicity predictions.

Specifically, the poly-3 statistic for a chemical obtained from the two-year cancer bioassays is assumed to be functionally related to the 90 day gene expression-based perturbation of the pathways (obtained from the SEPEA analysis):(3)where N_p represents the number of pathways (216, in the analysis in this paper), N_c the number of chemicals (26 chemicals were tested in mice). z_i is given by Equation (2) and Z_ij is given by Equation (1) for the pathway-specific predictor. Note, the model is fit to data sets corresponding to chemical treatments in Table 1.

The framework of the so-called Super Learner [34] coded in the SuperLearner package [35] in R [36] provides a reasonable way of evaluating alternate functional forms in a cross-validation framework. The functional forms tried include tree-based methods (random forests [37], bagging, conditional tree forests), support vector machines, loess polynomial regression, bayes generalized linear models, sparse partial least squares regression and neural-networks – the book [38] describes details of all of the these algorithms except random forests. Unless otherwise specified the default settings of these learning algorithms in the SuperLearner package are used. The SuperLearner algorithm that is based on the optimal continuous combination of the predictions of the other algorithms and the Discrete SuperLearner that picks the best predictor at each fold of cross-validation were used. The performances of these algorithms were evaluated in terms of their cross-validated risk [34].

Five-fold cross-validation was chosen where the chemicals in each of the five folds were fixed (Table S3) and chosen so that distribution of high and low carcinogenicity values in each fold was more or less the same. Predictions for all the data sets were based on this 5-fold cross validation framework. In order to get ‘honest’ predictions for the chemicals with data from mice, because these chemical data are used to make predictions for the very same set of chemicals, two levels of 5-folds cross-validation was implemented using the CV.SuperLearner function implemented in the SuperLearner package. Specifically, the predictions for the chemicals in each of the 5 test sets (corresponding to the 5-folds) are derived from an additional 5-fold cross-validation of the remaining 20 or 21 chemicals in the 5 training tests.

The two year cancer bioassay carcinogenicity calls of the chemicals used in the mice and rats experiments along with the continuous predictions from the model (Equation (3)) derived using the mouse data are used to generate receiver-operator (ROC) curves using the pROC package [39] in R [36]. The carcinogenicity calls (carcinogenic or non-carcinogenic) are based on a significance p-value threshold of 0.01 using the statistic given in Equation (2)) Myristicin and isosafrole were untested chemicals among the rat data and predictions for these were not used in generating the ROC curves. For the case of the human data, all gene expression data sets associated with risks for liver cancer were considered positive for carcinogenicity. The data set of gene polymorphisms associated with liver cancer was also considered positive while all other gene polymorphism data sets were considered negative. The ROC curve naturally defines notions of false positives and false negatives at chosen levels of specificity and sensitivity on the ROC curves. Additionally the area-under-the-curve (AUC) measures [40] of the ROC curves.

Different prediction models were evaluated in terms of their 5-fold cross-validated risk for fitting the model given by Equation (3) (see Table S4). Random forest and support vector machines using all the pathway features displayed the lowest cross-validated risk. Random forests [37] (denoted by SL.randomForest.1_All in Table S4) were because they provide a rigorous metric called importance for each of the biochemical pathways used in the prediction model. The importance of a pathway is a measure of how much it contributes to increases the accuracy of the predictions.

The model derived using the mouse data was used to predict the carcinogenicity in mice, rats and humans for the corresponding data sets. In this manuscript, the predicted values were treated as proxies for the continuous NTP carcinogenicity values in Equation (2). The receiver-operator curves for the three species are shown in Figure 1. The Area-Under-the-Curve (AUC) metric for both mice and rats (0.66 and 0.74) are reasonably good while it was 0.91 for humans. This is a very encouraging result suggesting the utility of extrapolation of stated mice results to the case of humans. The fact that the rat and human prediction performances were better than that for the mice may seem surprising and counter-intuitive. However, it should not be. For the data sets under consideration, the performance of predictions from test sets (in this case those from rats and humans) using the inferred prediction model depends on how closely the majority of positive (negative) test sets were to a majority of the elements of the positive (negative) training set (in this case those from mice). The ‘closeness’ here is measured in the pathway-response feature space. Therefore, the mechanisms of hepatocarcinogenicity implied by the pathway responses of the 26 chemicals in mice were enough to adequately cover the corresponding mechanisms seen in the data sets from rats and humans.

The false positives and false negatives in each of the three species at a level of specificity between 56–75% (i.e., comparable Type I error rate) is summarized are Table 4. A further discussion of the identified false positives and false negatives is provided in Text S1.

This paper presents results only for pathway-based prediction models. Gene-based prediction models were also evaluated and they showed a similar performance in terms of the AUC metric (data not shown) [12]. Evaluated gene-based models for predicting lung carcinogenicity using the same set of 26 chemicals as used in this paper. The average AUC metric they observed using a range of learning models was around 0.7.

The importance measures for each of the pathway-based features are reported in Table S2. The top fifteen pathways ranked by their importance measures output from the random forests learning algorithm are provided in Table 5. From the KEGG pathway database [17], [18], [19], the broad categories in which each of the pathways lie are also provided. Among these broad categories, altered fatty acid metabolism is associated with liver steatosis, leading to steatohepatitis and subsequently to an inflamed liver, liver cell death including apoptosis, inflammation, hepatocellular regeneration, stellate cell activation, and fibrogenesis, events that culminate in cirrhosis and liver cancer [41], [42]. This could also explain the reason why the gene expression data based on non-genetic risk factors like obesity, cirrhosis and diabetes were predictive of chemical-induced hepatocarcinogenesis. Calcium signaling is associated proliferating cells [43] and also along with cytochrome c in programmed cell death, apoptosis [44], [45]. Perturbed glycan synthesis has been found in ovarian cancer [46] and human mammary, colon [47], hepatic [48], [49], and glial tumors [50]. Other pathways included those associated with altered gene transcription and translation, xenobiotic, vitamin and carbohydrate metabolism. A clustergram of the z-transformed SEPEA p-values for all the pathways across the 26 chemicals used for the mouse data is shown in figure S1 (see Text S1 for details of generation of the clustergram). The lack of predictability especially for the case of mice (because prediction was based on mice data) suggests plausible new features that are not captured by the the pathway-based ones or that the set of 26 chemicals was not diverse enough in terms of having sufficient number of chemicals with alternate mode-of-actions. Some pathways such as Amyotrophic lateral sclerosis (ALS) may have relatively large importance measures but may have no direct significance to the biological processes in the liver. The interpretation one should have in such cases is that the responses of the set of genes associated with ALS disease process in the liver are relevant to the prediction of hepatocarcinogenicity.

The predicted dose response for the chemicals in Table 2 is shown in Figure 2. The slopes of the curves and the p-values of the alternate hypothesis that the slopes are positive are given in Table 6. NPTH is the only non-liver carcinogen among the five chemicals and had a negative slope. The responses derived from the model produce larger positive and more significant (lower p-values) slopes for the carcinogens. The fall in the carcinogenicity responses across the two lowest doses for PGBE and TCPN could be suggestive of alternate mode-of-actions at lower doses.

The analyses involved in this paper worked with hepatocarcinogenicity prediction models trained using data from mice and extrapolated to data sets from rats and humans. Alternate models derived from rats to predict carcinogenicity in other species could be attempted. This was not attempted here because of lack of sufficient analogous two-year cancer assay data for the rat and human data sets. The model derived here used data for chemicals at their maximum tolerable dose (MTD) levels. It could be argued that responses at lower doses are the ones with the most human relevance. Again this was not attempted here because of lack of sufficient number of chemicals dosed at lower levels. One possibility of understanding chemical-specific low-dose mechanisms would be the availability to a relatively large epidemiological data set of exposed people (e.g., [51] demonstrates an example of such a study for the case of workers exposed to benzene).

The mice used in the study lacked the genetic diversity of humans. So, the extrapolations of carcinogenicity predictions to human data has also to be viewed in the light of limitations in addressing concerns of identifying and characterizing the risk to susceptible populations. Populations could be considered susceptible either based on their genetics, age, nutritional and physical activity status and possibly other risk factors.

The B6C3F1 mice differ significantly from humans in the etiology of hepatocellular carcinoma (HCC). The mice have a relatively high background rate of developing HCC. In humans, HCC is thought to arise in a background of chronic inflammation, necrosis and regeneration, fibrosis and extracellular matrix deposition [52]. In comparison, HCC in mice are not known to arise under this background. Rather, a genetic event is assumed to precede a stepwise progression to HCC in mice. There are also known differences and similarities for the initiating genetic events for HCC between mice and humans. Therefore there is certainly a point in questioning the use of B6C3F1 mice when attempting to understand the etiology of human HCC. The utility of using the carcinogenicity results from existing rodent bioassays to the context of humans have been questioned by Bruce Ames and colleagues [53], [54], [55], [56]. This concern was in light of the relatively large proportion of tested chemicals being declared as rodent carcinogens. Some of these chemicals occur naturally in human dietary sources and others are being prescribed as pharmaceutical drugs. The exposures to humans from these sources were at doses much lower than what were used in the rodent bioassays. In fact the high positive rates from the rodent bioassays were hypothesized to have been caused by increased cell proliferation induced by the relatively high doses of the tested chemicals. Specifically for the data in this manuscript, none of the rodent liver carcinogens are currently known to be associated with human liver cancer.

However, in the predictive toxicology context of the manuscript, in spite of the known differences in etiologies, our hypothesis is that gene expression levels of the precancerous lesions are the similar in mice and humans. This hypothesis is validated by the relatively high predictability of human hepatocarcinogenicity using an independently derived predictive model using 90 day gene expression data from B6C3F1 mice. Further, the work in Hoeneroff et al suggests that this hypothesis is not without justification. In their work the authors found similar gene expression profiles in cells from HCCs obtained from mice and human samples. The situation is analogous to the case of azoxymethane induced colon cancer [57] in rodents being used as a model of human colon cancer. So the mode of action (chronic inflammation, increased proliferation or specific genetic events) leading to the pre-cancerous lesions in humans and rodents may be similar or different but once the precancerous state is reached then one could hypothesize that (at least evolutionarily) both species may follow similar paths to tumor formation and progression.

In summary, we rigorously derive and evaluate a biochemical pathway based hepatocarcinogenicity prediction model. Among the set of alternate prediction models , random forests were found to perform the best in terms of cross-validated risk. The model used gene expression data from a given tissue at the end of a short term study to predict the risk for development of tumors at the end of a longer period of time. Specifically, the model is evaluated using gene expression data obtained from the liver of mice treated with a range of 26 chemicals over a period of 90 days. The model with the information on affected pathways derived using these gene expression data had sufficient signal to adequately predict the two-year liver carcinogenicity risk of the same chemicals as evaluated in the National Toxicology Program's two-year cancer bioassay. The fact that the model was developed at the biochemical pathway level allows one to reasonably expect conserved behaviors of the chemicals at the pathway level across multiple species. This belief was validated using the model developed for mice to predict results in humans and rats. This fact was observed in Figure 1.The use of pathways in the carcinogenicity prediction model allowed a biologically based reduction of the feature space of the classifier. The responses of the various pathways suggested a complex interplay between them leading to the carcinogenicity prediction. What was encouraging was the excellent extrapolation for the case of human data.