Carcinogens are a large group of substances, organic and inorganic, that are directly involved in causing cancer. According to their mode of action they can be categorized as either genotoxic (GTX) or not genotoxic (NGTX). The former induces specific mutations or chromosome aberrations through direct interaction with DNA, usually by formation of covalent bonds (1). Such alterations are detected by a battery of tests (Ames test, chromosomal aberration, micronucleus assays) that measures the integrity and the structure of the DNA. Non-genotoxic drugs, on the other hand, represents chemicals capable of producing tumorigenesis by some secondary mechanism not directly related to DNA damage (2). Their activities are so diverse, that it is easier to define the properties they lack rather than the properties they possess. In general they are chemicals that do not induce DNA repair, and are negative in in vivo and in vitro tests for mutagenicity.

Non-genotoxic carcinogens have a wide variety of mechanisms of cancer induction including receptor mediated endocrine modulation, non-receptor mediated endocrine modulation, regenerative proliferation, oxidative stress, xenobiotic receptor activation, peroxisome proliferation, induction of inflammatory response and/or gap junction intercellular inhibition (3). Free radical production (particularly ROS) is a common sub-mechanism enhanced by several non-genotoxic carcinogens. Basically, cellular damage is promoted when the balance between pro and anti-oxidants is disturbed and the oxidants are not properly neutralized.

The diverse mechanisms of action, the tissue specificity and the lack of genotoxicity make non-genotoxic identification a challenging task. Rodent bioassays are considered the best available method for detecting such carcinogens. Risk assessment is done combining data from bioassays, epidemiological data, toxico-kinetic and disposition studies (3). The rationale behind this approach is that many of the drugs known to be carcinogens to humans are also carcinogens to animals. Classical studies in rats involve exposures for periods that range from 13 to 14 weeks. However, a proportion of chemicals are detected at the end of a 2 year period, making the animal chronic exposure assay elaborate and costly intensive.

A rapid and sensitive method for detecting hepatocarcinogenicity in drug screening is a long sought target. Control of gene transcription is the main regulatory mechanism of biological systems. Gene expression precedes protein synthesis, cell proliferation and ultimately pathological modifications. Therefore, it should be the most sensitive point to detect early changes (4). The aim of this analysis was to build a model that distinguishes non-genotoxic liver carcinogens by using expression profiles from short term exposure chemical treatments in rodents. Experimental data were obtained from the toxicogenomic database DrugMatrix™, The National Toxicology Program (U.S. Department of Health and Human Services) and The Toxicogenomics Project Genomics Assisted Toxicity Evaluation system (TG-GATEs) (5,6).

Experimental design and compounds

To evaluate molecular profiles, public available data from the National Toxicological Program (NTP) was selected (GEO Accesion number GSE57822). This entity performs pre-chronic and two year studies in laboratory animals in order to assess specific needs in toxicology, yielding the largest molecular toxicology reference. Briefly, arrays corresponding to 77 chemicals and their respective controls were downloaded from DrugMatrix (Table 1). Total data points were 363: three repeats per treatment involved 231 arrays and every treatment had 4 controls, in total 132 control arrays in 22 control groups. The Carcinogenic Potency Database was used as a first option to label the chemicals (7). Each array was obtained from test-compound treated and vehicle control-treated male rats after 72 hr of exposure with daily dosing (Sprangle-Drawley, 6~8 weeks old). Liver tissues (medial lobe) from three rats per chemical was collected and submitted to array processing. More data on the original experiments can be found in (8). Concentrations selected for each compound are summarized in Table 2.

The validation set was extracted from The Toxicogenomics Project Genomics Assisted Toxicity Evaluation system (TG-GATEs) (ExpressArray E-MTAB-800), a large-scale database of transcriptomics and pathology data useful for predicting the toxicity of new chemical entities (6). We downloaded data from rats exposed daily for 4 days at three doses (low, middle and high). Four hundred and seven arrays corresponding to 34 chemicals and their corresponding controls were obtained. Table 1 shows selected drugs and Table 3 their dose levels.

Data pre-processing

Complete “.CEL” files were downloaded from the National Toxicological Program, Department of Health and Human Services (USA) and the National Bioscience Database Centre (National Bioscience Database Center, Tokyo, Japan). Files belonged to the Affymetrix Rat Genome 230 2.0 GeneChip Array (Affymetrix, Santa Clara, CA, USA). Preprocessing adjustments were performed with Expression Console (Affymetrix). Additional information and raw data from the public repositories can be found at https://ntp.niehs.nih.gov/drugmatrix/index.html and http://dbarchive.biosciencedbc.jp/en/open-tggates/download.html/.

Each array (CEL file) was preprocessed and background corrected, normalized and summarized using RMA (Robust Multiarray Average) using Expression Console (Affymetrix) and Bioconductor packages of the R software (Fred Hutchinson Cancer Research Center, Seattle, USA) (9). For probe filtering, unspecific selection was carried out according to the interquartile range (IQR) (cut off value according to the IQR density plot).

Genes that had more than 1.5 fold increase/decrease relative to controls were chosen for further analysis. Additionally, each drug treatment was compared with its respective control using t test statistic. Results were further corrected for multi-testing using the Benjamini & Yekutieli (2001) procedure for (conservative) control of the false discovery rate (FDR), with 0.05 as the significance level (10). Features that met both criteria (t-test and 1.5 fold change) were combined in a single list of differentially expressed genes that resulted in 3778 probes that underwent classification analysis.

Class discrimination

Pre-processing

After normalization and log transforming the data, unspecific filtering was applied to the 31099 probes, leaving 10091 features. Differential expressed genes were identified for each treatment when compared with the set of corresponding controls by t-test and fold change. A unique DEG list from all treatments consisting of 3778 probes was built. The filtered DrugMatrix set was used for Principal Component Analysis (PCA). Fig. 1 shows a two-dimensional plot of the data. Each color represents a single compound defined by the rainbow pallet. Variance was relatively low for the first four components: PCA1 (12%), PCA2 (9%), PCA3 (8%), PCA4 (4%). Overall, there was no significant clustering among NGTX and other compounds, suggesting that further filtering steps were required to make a successful discrimination.

Feature selection by random forest

Random Forests analysis combined with a feature selecting algorithm was used to build the predictor. In the training set, the out of the bag error (OOB) for the initial random forest was 0.09 (81.5% sensitivity, 100% specificity). Fifty four variables were selected by variable selection without dropping the OOB error substantially (Table 4). Fig. 2 shows a Principal Component Analysis (PCA) constructed with the selected genes as variables. Each compound was averaged among replicates. The proportions of variance for the first three components were: PC1 39%, PC2 12% and PC3 9%. As can be seen on the plot, NGTX carcinogens were clustered at the right of the PC1 axis, with the exception of Carbon Tetrachloride. Fig. 3 represents the mean expression difference for each one of the predictor genes.

In order to assess an honest prediction of the error rate from the training data we performed bootstrapping (0.632) in which the random forest constructed for a certain number of variables was subsampled and compared. The prediction error rate among the bootstrap samples was 0.15 (Supplemental Fig. 1). Thirty genes were consistently selected (stability) in all the sub-samples (above 20%), with CYP1A2 being selected in 65%, Prodh in 32% and Itgb3 in 30%.

Prediction in the test set

The expression profiles of the test sample group were run through the obtained random forest. Overall, prediction error was lower than expected, 0.13. Sensitivity was 97% and specificity 81%. Three agents were misclassified as NGTX carcinogens: Clotrimazole, Stavudine and Dexamethasone. We found that six drugs had similar profiles to our predictor in the unclassified group: Aminoglutethimide, Crotamiton, Finasteride, Modafinil, Rosiglitazone and Simvastatin. Table 2 shows the results of the random forest in the DrugMatrix set.

Validation set

The performance of the classifier was tested in an independent dataset. We applied the random forest model to assess treated rat livers for 4 days and 34 drugs from the TG-GATE database. In contrast to DrugMatrix, three different doses were available at this time point (low, middle and high). Prediction error varied significantly among dose levels. While all concentration levels yielded specificities of 100% (no false positives), sensitivity increased from 10% at low doses and 38% at medium doses, to 61% at high doses. Sensitivity could be increased at expense of specificity. Areas under the curve (AUC) for the ROC curve (Receiver Operating Characteristic) were 0.73 for low doses, 0.83 for middle doses and 0.87 for high doses. The results of the classification model applied to the validation set are described in Table 3.

Pathway analysis

In order to evaluate the biological pathways involved in the expression of the 54 highest scoring genes, we performed functional analysis of the GO terms that were significantly represented. Differential expressed genes were used as background list. We have found that significant genes were scattered among diverse GO Biological Process, most of them related to mitochondrial respiration, energy production and lipoprotein metabolism. Detoxification was located later in the list. Molecular functions included several terms of signal transducing and receptor activities. Table 5 shows the top eight GO annotations according to Fisher’s exact test for BP (Biological Processes) and MF (Molecular Functions). There is a representation of how these terms are distributed over the GO graph in Supplemental Fig. 2.

Toxicogenomics is defined as the application of genomics technologies to study the adverse effects of pharmaceutical and environmental chemicals, with the hope of improving risk assessment and hazard screening (5). In this paradigm, public databases are essential tools for multiple actions, such as comparing profiles, discovering patterns or integrating networks. The largest databases available to date are the Japan Project Database (TG-GATE) and DrugMatrix (National Toxicological Program, NTP, USA) (5,6). Their diversity and standardized protocol make them the current reference. In recent years other databases were made available (CBES, CTD, etc.) and the challenge has been extended to data mining (15).

The mechanisms of non-genotoxic carcinogenicity are well described, although not well understood. It is assumed that most non-genotoxic carcinogens induce neoplasm and exhibit threshold tumor dose-responses. Classically, NGTX are regarded as tumor promoters but mechanisms such as regenerative hyperplasia, cytotoxicity or induction of oxidative damage are also key events for tumor initiation (16). There is still uncertainty whether NGTX are capable of initiate tumor events by themselves or whether they need coincidental factors. In this study, we focused our efforts to rodent hepatocarcinogens, a common target of safety assessment in drug development.

One of the most, if not the most, important step in class prediction is the correct assignment of the training set. New experimental data may change the status of a substance to another category limiting applicability of the predictor. Generalization of data is also critical, and in agreement with previous studies, false positives are often difficult to avoid (or “unavoidable”). In order to evaluate a complex process like carcinogenesis, we have selected compounds with three modes of action to train our model: oxidative stressors, peroxisome proliferators and hormone modulators. Our hypothesis was that in some point, similarities among gene profiles would allow discrimination between non-genotoxic carcinogens and other class of chemicals.

We chose Random Forest as our machine learning method because is a robust classification algorithm. Feature reduction was performed by recursive variable elimination while maintaining class error (11). We report a predictor of 54 variables, a number that may be reduced to 24 without losing much prediction accuracy. After training the model, we applied our classifier to the DrugMatrix test set, which resulted in a relatively low total error. Only three false positives were detected and no false negatives among the OOB samples. The small sample size could explain why the prediction error of the test set was lower than the estimation of the training set (0.13). It is interesting that prediction accuracy increased with increasing doses in the validation set. This behavior could be explained by the fact that DrugMatrix doses for 3 day treatments were higher than TG-GATE in almost all compounds, frequently several times higher than the therapeutic dose.

One of the false positives detected was clotrimazole (CTL), an imidazole antimitotic. It is an inhibitor of p450, blocker of intracellular Ca⁺⁺ stores and activator of the xenobiotic response (17). It has important correlation with our predictor, particularly in the induction of ZDHHC2 (Palmitoyltransferase ZDHHC2), a protein responsible for membrane binding and protein localization (receptors). Imidazole agents modulate p450 response by different gene subfamilies and display different actions (18). The other two miss-classified drugs were Dexamethasone (DXM) and Stavudine (d4T). When compared to our predictor, we found proximity in genes related to stress and inflammation (Cyp1a2, IL1A, COL4A3BP, SDC2, PTTG1, C5ar1, CFH, Cfi).

In line with previous studies, a mechanistic approach is a promising strategy for prediction as well as for pathway or functional category analysis. In this study, the most significant biological processes comprised those related to cellular respiration, energy derivation by oxidation of organic compounds and generation of precursor metabolites. This observation is in agreement with the known action of nongenotoxins and the role oxidative stress in carcinogenesis (3,19). Additionally, the presence of genes related to lipid metabolism is typical of profiles delivered by peroxisome proliferators. Others genes involved complement cascade, inflammation response, and some of them were related to cellular attachment.

Our signature included synthetic sex steroids (cyperone acetate, 17-methyltestosterone, diethylstilbestrol and ethinylestradiol). Several of them were shown to produce liver tumors in rats when given at therapeutic doses (20). Initially, the non-genotoxic effects of steroids were thought to be triggered by downstream genes of specific receptors, although increasing data suggest that they possess genotoxic action as well. We decided to keep them as NGTX carcinogens because the effect was complementary. Increased CYP activities and alterations in sterol metabolism are frequently associated with hepatomegaly in a non-genotoxic manner (3).

A steroid responsive gene, Interleukin-1a (IL-1a), is a high scoring gene in our signature. It is a critical cytokine whose expression is related with various aspects of human reproduction and expressed in a number of solid tumors. IL-1a serves as attractant by lymphocytes that keep the inflammation state that precedes malignancy. Sex steroid receptors downregulate and confine IL-1a expression, but its deregulation is a key inducer of proteolytic enzymes that degrade the extracellular matrix and remodel tissue (21).

Several studies have provided gene predictors for nongentoxic hepato-carcinogens (22?28). Most of the methodologies used were based on support vector machines (SVM) coupled with a feature reduction algorithm. Published gene classifiers ranged from 9 to more than 100 probes and prediction accuracies were roughly equivalent. In general, data from rats under longer exposures provided higher accuracies. The composition of predictors differed among studies in length and diversity and only few genes overlapped. The reason for such heterogeneity is given by the differences in experimental designs. However, all these factors make biological interpretability more complex.

A recent study developed by Gusenleitner et al. (2014) scanned most of the NTP database and built a classifier that englobed liver carcinogenesis (combining both genotoxic and non-genotoxic carcinogens) (29). They validated their model with TG-GATE liver data and estimated liver carcinogenesis with an AUC of 0.78 (56.8% sensitivity and 82.91% specificity). It is interesting that our signature, although confined to a non-genotoxic subgroup of short term exposure, had in common CYP1a2 and ZDHHC2 as high scoring genes. Also, the presence of members of the superfamily of solute carrier proteins (SCL genes) in both signatures indicates the importance that membrane transporters have for drugs that “hitch-hike” one or another to enter the cells (30). On the other hand, the inclusion of genotoxic compounds to predict global carcinogenesis resulted in a signature with several cell cycle control genes, essential to the regulation of growth and apoptosis during mutagenic stress.

Toxicogenomics is a promising field with the hope of aid both the earlier elimination of toxic compounds in the drug pipeline and the discovery of new toxicity mechanisms. The biggest target is accurately predict the toxicity of unknown drugs. In this analysis, we presented a classifier that can predict non-genotoxic carcinogenicity by using short term exposure treatments from the NTP database. Although we are aware that the classifier does not have the prediction accuracy of signatures of long term exposure, early screening is an advantage that would allow prioritizing compounds for further testing.

We would like to thank to M.Sc. Patricia Dell’Arciprete (National Patagonic Center-CONICET) for help in data processing.