Biological Networks for Predicting Chemical Hepatocarcinogenicity Using Gene Expression Data from Treated Mice and Relevance across Human and Rat Species

Background Several groups have employed genomic data from subchronic chemical toxicity studies in rodents (90 days) to derive gene-centric predictors of chronic toxicity and carcinogenicity. Genes are annotated to belong to biological processes or molecular pathways that are mechanistically well understood and are described in public databases. Objectives To develop a molecular pathway-based prediction model of long term hepatocarcinogenicity using 90-day gene expression data and to evaluate the performance of this model with respect to both intra-species, dose-dependent and cross-species predictions. Methods Genome-wide hepatic mRNA expression was retrospectively measured in B6C3F1 mice following subchronic exposure to twenty-six (26) chemicals (10 were positive, 2 equivocal and 14 negative for liver tumors) previously studied by the US National Toxicology Program. Using these data, a pathway-based predictor model for long-term liver cancer risk was derived using random forests. The prediction model was independently validated on test sets associated with liver cancer risk obtained from mice, rats and humans. Results Using 5-fold cross validation, the developed prediction model had reasonable predictive performance with the area under receiver-operator curve (AUC) equal to 0.66. The developed prediction model was then used to extrapolate the results to data associated with rat and human liver cancer. The extrapolated model worked well for both extrapolated species (AUC value of 0.74 for rats and 0.91 for humans). The prediction models implied a balanced interplay between all pathway responses leading to carcinogenicity predictions. Conclusions Pathway-based prediction models estimated from sub-chronic data hold promise for predicting long-term carcinogenicity and also for its ability to extrapolate results across multiple species.


Introduction
There is push to use a broad array of biological data in toxicity testing to improve reliability and provide reasonably quick indications of animal and human toxicity of chemical compounds. The National Toxicology Program (NTP) 's High Throughput Screening (HTS) [1] and the Environmental Protection Agency's ToxCast programs [2] are efforts in this direction. These efforts have been initiated in light of the fact that for a lot of the 100,000 chemicals in commerce in the U.S. and Europe, information on toxicity is extremely limited [3]. Data provided from the long term chronic testing of chemicals on animals, like those from NTP's two-year cancer bioassay though invaluable in providing carcinogenicity information in vivo over a wide range of animal tissues and across multiple doses, is not suitable when one needs to rapidly identify potential harm from chemical exposure. Over the past approximately 40 years, only around 600 chemicals have been tested by the NTPs in their two-year cancer bioassay and only about 5-6 times this number have been tested worldwide.
The goal of toxicity testing is ultimately to protect human health. Even though the two-year cancer bioassay is performed on rats and mice, there is considerable data supporting its use to screen for carcinogens in humans [4]. The HTS and ToxCast programs are performed on cell lines; these can be both human and animal and thus provide a logical path for extrapolation between species. However, cell lines are distantly removed from a functional whole organism so it becomes necessary to extrapolate the in vitro results from the various cell lines to the human in vivo context.
An alternative to the cell-based assays in the HTS and ToxCast programs is the use of short term molecular data (typically gene expression) from exposed animals or humans to predict future toxicity or carcinogenicity. Others have used gene expression data based classifiers that would distinguish a toxin or carcinogen from a non-toxin or non-carcinogen [5,6,7,8,9,10,11,12,13,14]. These representative examples demonstrate different approaches to predictive toxicology using gene expression data obtained from rodents dosed over a period of time ranging from a day to 90 days. In all of the published studies, chemicals were classified in a binary fashion (for example, attempts were made to classify chemicals as being carcinogenic or not). Such an approach underutilizes the results from the chronic studies in light of the fact that a continuous statistic was used to evaluate the data originally (e.g. pvalue from the Fisher exact test or the trend statistic from the Poly-3 test [15]) and a considerable amount of other biological information was brought to bear on a final evaluation. This leads to obvious loss of information regarding the differential toxicity of different compounds at their relevant doses.
A characteristic of gene expression data is the very large features (genes) to samples (animals) ratio. This characteristic contributes to the lack of identification of robust classifiers as demonstrated by [16] -multiple classifiers give the same prediction accuracies which quite often are 100 percent. This could mean that the data was over-fit or there are correlations between genes that allow for alternate sets of genes to server as predictors. In the second scenario of correlated gene expressions, this would suggest the use more of mechanistically-relevant prediction models where the data is examined at the biological pathway level rather than at the gene level.
In order to gauge quantitative predictive accuracy of a pathwaybased prediction model we obtained gene expression data from the livers of mice dosed with 26 chemicals over a period of 90 days. All the 26 chemicals were tested in a two-year bioassay in the same strain of mice dosed via the same respective routes. Thus, the differential liver carcinoma and adenoma rates between treated and control animals from the results of the two-year bioassays are known. Using the data, we derived pathway-based models to predict the differential tumor rates at the end of a two-year bioassay. The use of a pathway-based model would result in a reduction of a model with around 6000 gene features to one with around 200 pathways in the KEGG database [17,18,19]. The model derived using the mouse data was then used to predict carcinogenicity in appropriately chosen scenarios in rats and humans.

Mice chemicals, animals and treatments
The chemical treatments used in the experiments on female B6C3F 1 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.

Rat liver gene expression data
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).

Human gene data
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 2 -transformed. The entire set of data sets associated with risk factors for liver cancer is summarized in Table 3.

Biochemical pathways used in the analysis
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.
Structurally enhanced pathway enrichment analysis (SEPEA) 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 5610 4 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, where W refers to the standard normal distribution and P ij denotes the p-value obtained using SEPEA for pathway j as a result of the The results for liver tumors in this study were considered equivocal or borderline significant. Combined increase in hepatocellular adenomas or carcinomas resulted in p = 0.075 and p = 0.084 for benzene and coumarin respectively. d Due to signs of toxicity, the 16,000 ppm dose was reduced to 0 ppm on day 9 for a period of 2 days. The dose was raised to 8,000 ppm for a period of 9 days and returned to 16,000 ppm for the remainder of the study. The time weighted average dose was 14,800 ppm. e The initial dose of 3,000 ppm was reduced to 2,000 ppm in week 2 of the study due to taste aversion and weight loss. The 2,000 ppm dose is the same as the low dose in the original bioassay. f Chemical not evaluated by the NTP. Bioassay performed by Alexander et al. [20]. doi:10.1371/journal.pone.0063308.t001 mice being dosed with chemical i at a given dose or human data set i.

Evaluation of carcinogenicity for mice chemical data
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 twoyear cancer bioassay developed by the NTP(see Text S1 for reference to the NTP technical reports). The poly-3 survivaladjusted 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, 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 2e25 was used when z i $z 2e-5 where z 2e25 is the critical value of the standard normal distribution corresponding to 2610 25 significance.

Prediction model relating carcinogenicity to either pathway features
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): 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 neuralnetworks -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 Super-Learner 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.

Evaluation of the results of the predictions
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-underthe-curve (AUC) measures [40] of the ROC curves.

Choice of prediction algorithm
Different prediction models were evaluated in terms of their 5fold 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.

Receiver-operator curves for the predictions in the three species
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.

Features of the prediction models
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 fibrogen-esis, 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 ztransformed 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 Humans) The legend in the sub-plots provides the area-under-thecurve (AUC) for the corresponding ROC curve. The curves for mice, rats and humans are based on datasets corresponding to those in Tables 1,  Table S1 and Table 3 respectively. doi:10.1371/journal.pone.0063308.g001 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.

Dose-response predictions
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 pvalues) 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.

Conclusions
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 Table 4. False positives and false negatives predictions of predictors across the three species at appropriately chosen points on the receiver-operator curves in Figure 1. 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.

ROC curve
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 Figure 2. Dose-response predictions using the pathway-based prediction models for data of the chemicals tabulated in Table 2. Dose response predictions for five chemicals (a. dcbz, b. mecl, c. npth, d. pgbe and e. tcpn) treated at four different doses in mice. doi:10.1371/journal.pone.0063308.g002 Table 6. Slopes of the dose-response curves in Figure 2. 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. Figure S1 Clustergram of transformed p-values (Equation (1)) representing the enrichment of the 216 pathways across the 26 chemicals treatments in mice, was generated using hierarchical clustering with the euclidean distance metric and average linkage to generate the hierarchical trees of pathways and chemicals using the Cluster and Tree view programs [58].

Supporting Information
(PNG) Text S1 Provides additional description of the microarray experiments, involved in mice and rats, the data for which are used in the analysis. Also provided are details of generation of Figure S1, a discussion of the identified false positives and false negatives by the prediction models, primers on the technique of cross-validation , receiver-operator curves and references to the NTP technical reports of the two-year cancer bioassays for each of the 26 mice treated chemicals.
(DOC)  Table 1 that were chosen to be in each of the 5 folds in the cross-validation analysis performed.