A predictive toxicogenomics signature to classify genotoxic versus non-genotoxic chemicals in human TK6 cells

Genotoxicity testing is a critical component of chemical assessment. The use of integrated approaches in genetic toxicology, including the incorporation of gene expression data to determine the DNA damage response pathways involved in response, is becoming more common. In companion papers previously published in Environmental and Molecular Mutagenesis, Li et al. (2015) [6] developed a dose optimization protocol that was based on evaluating expression changes in several well-characterized stress-response genes using quantitative real-time PCR in human lymphoblastoid TK6 cells in culture. This optimization approach was applied to the analysis of TK6 cells exposed to one of 14 genotoxic or 14 non-genotoxic agents, with sampling 4 h post-exposure. Microarray-based transcriptomic analyses were then used to develop a classifier for genotoxicity using the nearest shrunken centroids method. A panel of 65 genes was identified that could accurately classify toxicants as genotoxic or non-genotoxic. In Buick et al. (2015) [1], the utility of the biomarker for chemicals that require metabolic activation was evaluated. In this study, TK6 cells were exposed to increasing doses of four chemicals (two genotoxic that require metabolic activation and two non-genotoxic chemicals) in the presence of rat liver S9 to demonstrate that S9 does not impair the ability to classify genotoxicity using this genomic biomarker in TK6cells.


a b s t r a c t
Genotoxicity testing is a critical component of chemical assessment. The use of integrated approaches in genetic toxicology, including the incorporation of gene expression data to determine the DNA damage response pathways involved in response, is becoming more common. In companion papers previously published in Environmental and Molecular Mutagenesis,   [6] developed a dose optimization protocol that was based on evaluating expression changes in several well-characterized stress-response genes using quantitative real-time PCR in human lymphoblastoid TK6 cells in culture. This optimization approach was applied to the analysis of TK6 cells exposed to one of 14 genotoxic or 14 non-genotoxic agents, with sampling 4 h post-exposure. Microarray-based transcriptomic analyses were then used to develop a classifier for genotoxicity using the nearest shrunken centroids method. A panel of 65 genes was identified that could accurately classify toxicants as genotoxic or non-genotoxic.  [1], the utility of the biomarker for chemicals that require metabolic activation was evaluated. In this study, TK6 cells were exposed to increasing doses of four chemicals (two genotoxic that require metabolic activation and two non-genotoxic chemicals) in the presence of rat liver S9 to demonstrate that S9 does not impair the ability to classify genotoxicity using this genomic biomarker in Value of the data The data were integral in developing a genomic biomarker-based approach to predict genotoxicity in human TK6 cells based on expression profiles induced by 28 chemicals that span a variety of well-defined genotoxic and non-genotoxic modes of action.

Contents lists available at
The data also demonstrate that the use of S9 does not alter gene expression changes used to classify genotoxicity in TK6 cells, expanding on the test agents applied using the biomarker.
The biomarker and these original training data sets can serve as a basis for testing new chemicals for genotoxicity using DNA microarray and other genomics platforms, in other cell types, and potentially in alternative organisms.
The database can be expanded to develop signatures that delve into more detailed genotoxic modes of action (e.g., signatures for cross-linking agents or other types of DNA lesions represented in the training set).
We anticipate that the importance of toxicogenomics studies in chemical risk assessment will continue to increase in the coming years and believe that the rate at which this occurs will be highly dependent upon ensuring public availability of these very powerful datasets sets and tools such as those described.

Data
The training set (GSE58431) consists of transcriptional changes in TK6 cells following exposure to a diverse set of model agents that include: DNA alkylators, DNA strand breaking agents, topoisomerase inhibitors, nucleotide antimetabolites, endoplasmic reticulum (ER) stressors, energy metabolism inhibitors, histone deacetylase (HDAC) inhibitors, microtubule inhibitors, and heavy metals. From these data, a panel of 65 genes was identified that could accurately classify toxicants as genotoxic or non-genotoxic.
As TK6 cells are not metabolically competent, the test data (GSE51175) examine the utility of the biomarker for use with chemicals requiring metabolic activation in order to broaden the biomarker's application [1]. Here, chemical exposures were conducted in the presence of rat liver S9.

Experimental design
Transcriptome measurements were performed using a two-color dye swap design [4]. For each agent, cyanine-3 and cyanine-5 labeled cRNA for treated and untreated samples were co-hybridized to two Agilent oligonucleotide microarrays. One microarray was hybridized with the treated cyanine-3 labeled cRNA co-hybridized with the untreated cyanine-5 labeled cRNA. The second slide was hybridized with the treated and untreated samples having the reversed cyanine-3 and cyanine-5 labeling. A summary of the workflow for the subsequent steps are presented as Fig. 1.

Hybridization and data quantification
Hybridization and washing was performed according to the manufacturer's protocol. Arrays were scanned with an Agilent DNA microarray scanner. Feature Extraction (Version 9.1; Agilent) was used to quantify the scanned images for the 4x44k platform, to calculate the processed signal intensity and to estimate the log ratios (log 10). Similarly for the 8x60k platform, Agilent Feature Extraction software (version 11.0.1.1) was used for quantification of the generated image files, to generate the QC reports as well as the processed signal intensity and to estimate the log ratios (log 10).

Data processing and normalization
Microarray annotation files for the Agilent human 1x44k (eArray 1x44K human_20120130.txt), 4x44k (eArray 4x44K human_20120130.txt), 4x44k version 2 (eArray 4x44K human_v2_20120130. txt), 8x60k (eArray 8x60K human_20120411.txt), and the 8x60k version 2 (eArray 8x60K human_v2_20120628.txt) platform were obtained from eArray (https://earray.chem.agilent.com/ earray/). These annotation files were read into the R statistical environment [6]. All unique Agilent probe ids from these platforms were used to create an annotation file using the Probes.R script named Agilent annotation.txt. This annotation file was then merged with the Feature Extraction data files to collapse the probe ID's to the gene symbol.
Agilent Feature Extraction.TXT files were read into the R environment using the Read.R script. This R script extracts the normalized relative intensities, merges the Agilent annotation.txt file using the probe ids and uses the weighted mean to average the probes with multiple gene symbols. Once all the data files have been processed, the individual data sets were then merged together using the gene symbol. The data were merged together such that the dye swaps were in adjacent columns. The dye swaps were then averaged and the processed data file named LogRatio.txt was written out to disk.

Development of the TGx-28.65 genotoxicity classifier
The original 65-gene signature, as published in Li et al. [5], referred to as TGx-28.65 (28 refers to the use of 28 chemicals in the training set) was developed using the nearest shrunken centroids approach [7]. Gene expression data were exported from Rosetta Resolver based on Entrez Gene Identifiers. The data were then read into the R statistical environment employing the pamr package [3]. The standardized centroids to predict whether an agent is genotoxic or not (either directly or indirectly) was computed for each class using a training set of agents, as described by Li et al. Briefly, the standardized centroid is the mean expression level for each gene in a class divided by its within-class standard deviation. The standard centroid for each class is then shrunken toward the overall centroid to produce the nearest shrunken centroid. The method employs a shrinkage parameter that is used to control the number of features used to construct the classifier.
The shrinkage parameter was identified by employing 10-fold cross validation [2]. This is done by randomly assigning samples to one of ten approximately equal-sized sets, which is also roughly balanced for the two classes. Prediction accuracy is assessed for each set with the other nine sets used to construct a classifier. From this analysis a shrinkage threshold of 2.2 produced a 65-gene panel with 100% accuracy based on 10-fold cross-validation.
2.5. Application of the TGx-28.65 signature to TK6 cells co-exposed with rat liver S9 Due to revisions in the gene annotation and differences between the 4x44k and 8x60k platforms, two of the 65 genes in the TGx-28.65 biomarker were unavailable when the þS9 data were incorporated into the study. As a result, MGC5370 (due to the annotation update) and USP41 (not present on the 8x60k platform) were removed from the TGx-28.65 biomarker.
Using the remaining 63 genes, the centroids were re-estimated using the PAM.R script. This R script reads in the LogRatio.txt data into the R statistical environment. Using the training set, the pamr package was used to update the classifier. These were written to disk using the filename Classifier.txt. The summarized data are presented in Table 1. Removing MGC5370 and USP41 did not impact the performance of the classifier as the accuracy remained 100% (estimated using 10-fold cross-validation).

PCA analysis
A principle component analysis (PCA) was performed using the prcomp function [8] in R on the training set data (Li et al. data) from the LogRatio.txt data file. Using the 63 genes in the TGx-28.65 classifier the PCA Table 1 TGx-28.65 classifier.
Presented are the estimated standard deviation and the shrunken centroids for the genotoxic and non-gentoxic classes.   was conducted using the PCA.R script. A summary of the PCA results are presented in Table 2 and the scatter plot of the first two principle components for the training set are presented as Fig. 2. The red line in Fig. 2 was added as the first principle component and could be interpreted as a contrast between the genotoxic (red) and non-genotoxic (blue) classes as there is a clear separation between the two classes. The PCA loadings obtained from this analysis was applied to the þS9 data. A scatter plot of these data for the first two principle components is presented as Fig. 3. The data from the training set are represented by the red and blue circles, whereas the labels for the þS9 data are displayed with red font for the genotoxic agents and blue font for the non-genotoxic agents. For the þS9 data, the low dose samples for the genotoxic agents were displayed in gray as their genotoxicity was uncertain. Similar to Fig. 2, the genotoxic agents and non-genotoxic agents separate into two groups with the low dose genotoxic agents falling in the middle.

Conflicts of interest
None.