Scholarly article on topic 'BMDExpress Data Viewer - a visualization tool to analyze BMDExpress datasets'

BMDExpress Data Viewer - a visualization tool to analyze BMDExpress datasets Academic research paper on "Biological sciences"

Share paper
Academic journal
Journal of Applied Toxicology
OECD Field of science

Academic research paper on topic "BMDExpress Data Viewer - a visualization tool to analyze BMDExpress datasets"

Research article

Received: 23 September 2015, Revised: 15 October 2015, Accepted: 16 October 2015 Published online in Wiley Online Library: 15 December 2015

( DOI 10.1002/jat.3265

BMDExpress Data Viewer - a visualization tool to analyze BMDExpress datasets

Byron Kuoa, A. Francina Webstera,b, Russell S. Thomasc and Carole L. Yauka*

ABSTRACT: Regulatory agencies increasingly apply benchmark dose (BMD) modeling to determine points of departure for risk assessment. BMDExpress applies BMD modeling to transcriptomic datasets to identify transcriptional BMDs. However, graphing and analytical capabilities within BMDExpress are limited, and the analysis of output files is challenging. We developed a web-based application, BMDExpress Data Viewer (, for visualizing and graphing BMDExpress output files. The application consists of "Summary Visualization" and "Dataset Exploratory" tools. Through analysis of transcriptomic datasets of the toxicants furan and 4,4'-methylenebis(N,N-dimethyl)benzenamine, we demonstrate that the "Summary Visualization Tools" can be used to examine distributions of gene and pathway BMD values, and to derive a potential point of departure value based on summary statistics. By applying filters on enrichment P-values and minimum number of significant genes, the "Functional Enrichment Analysis" tool enables the user to select biological processes or pathways that are selectively perturbed by chemical exposure and identify the related BMD. The "Multiple Dataset Comparison" tool enables comparison of gene and pathway BMD values across multiple experiments (e.g., across timepoints or tissues). The "BMDL-BMD Range Plotter" tool facilitates the observation of BMD trends across biological processes or pathways. Through our case studies, we demonstrate that BMDExpress Data Viewer is a useful tool to visualize, explore and analyze BMDExpress output files. Visualizing the data in this manner enables rapid assessment of data quality, model fit, doses of peak activity, most sensitive pathway perturbations and' other metrics that will be useful in applying toxicogenomics in risk assessment. © 2015 Her Majesty the Queen in Right of Canada. Journal of Applied Toxicology published by John Wiley & Sons, Ltd.

H Additional supporting information may be found in the online version of this article at the publisher's web-site

Keywords: microarray; gene expression; genomics; transcriptomics; benchmark dose (BMD); BMDExpress; bioinformatics; data visualization; human health risk assessment; furan; 4,4'-methylenebis(N,N-dimethyl)benzenamine (MDMB); dose-response


Chemical risk assessment aims to establish acceptable levels of exposures based on toxicological dose-response studies. Traditional methods that apply the lowest-observed-adverse-effect-level or no-observed-adverse-effects-level, may be limited by the selection of doses, sample sizes required to detect subtle effects and by technical and biological variability that limits ability to detect significant changes (Crump, 1984). In contrast, benchmark dose (BMD) modeling fits experimental dose-response data with a statistical model to identify a defined level of response relative to a control group. BMD was developed to overcome the limitations of the lowest-observed-adverse-effect-level/no-observed-ad-verse-effects-level approach (Crump, 1984). Regulatory agencies have increasingly adopted BMD modeling for human health risk assessment (Budtz-Jorgensen et ai, 2013; Health Canada, 2013).

Toxicogenomic studies use genomics technologies, such as DNA microarrays and RNA-sequencing, to investigate global transcriptional responses following chemical exposures. It has been proposed that BMD modeling be applied to model the dose-response of toxicogenomics data to derive BMD values for transcriptional endpoints, including genes, pathways and gene ontologies (Thomas et ai. 2007). This approach has already been applied to a variety of chemicals (Bhat et ai., 2013; Moffat et ai., 2015; Thomas et ai., 2013b). Importantly, these studies have shown concordance between the BMD values modeled for transcriptional and apical endpoints. However, this type of work is very much in its infancy, and a significant amount of research is required to

determine the best approaches for use in practical regulatory science or product development settings.

The standard software BMDS (US EPA, 2015), which was developed to model the dose-response of apical data, is inefficient for modeling the response of transcriptomics datasets because these contain gene expression profiles of tens of thousands of genes. To address this issue, Yang et ai. (2007) developed a tool called "BMDExpress" that automates the statistical analysis and BMD modeling computational steps to facilitate timely analysis of very

Correspondence to: Caroie L. Yauk, Environmentai Heaith Science and Research Bureau, Heaith Canada, Tunne/s Pasture, Bidg. 8 (P/L 0803A), 50 Coiumbine Driveway, Ottawa, Ontario K1A 0K9, Canada E-maii:

aEnvironmentai Heaith Science and Research Bureau, Heaith Canada, Ottawa, ON, Canada, K1A 0K9

bDepartment of Bioiogy, Carieton University, 1125 Coionei By Drive, Ottawa, K1S 5B6, Canada

c United States Environmentai Protection Agency, Office of Research and Deveiop-ment, Nationai Center for Computationai Toxicoiogy, Research Triangie Park, NC 27711, USA

Reproduced with the permission of the Minister of Heaith Canada.

This is an open access articie under the terms of the Creative Commons Attribution-NonCommerciai-NoDerivs License, which permits use and distribution in any medium, provided the originai work is properiy cited, the use is non-commerciai and no modifications or adaptations are made.

J. Appl. Toxicol. 2016; 36:1048-1059

BMDExpress Data Viewer


large transcriptomic datasets. In addition to modeling the dose-response of individual genes, BMDExpress also computes BMD values for biological categories (e.g., pathway BMD means, medians, 5th percentiles, etc.) using a feature called "Functional Classifications." Global transcriptional BMD modeling using BMDExpress has been used to demonstrate that toxicogenomic BMDs fall within the range of BMDs derived from apical data, and has been useful for investigating the mode of action of several toxicants (Bhat et al, 2013; Bourdon et al, 2013; Jackson et al., 2014; Moffat et al., 2015). These studies support the notion that transcriptional BMDs can be used to derive points of departure for human health risk assessment, particularly in instances when apical data are not available (Chepelev et al., 2015; Moffat et al., 2015; Thomas et al, 2012, 2013a, b).

BMDExpress computations are presented in tabular format viewable in the BMDExpress software. However, because of the limited capability to perform additional analyses and data visualization in the BMDExpress application, the results are typically exported to separate software (e.g., a spreadsheet) for further exploration. There are two major types of outputs that can be exported, i.e., (1) "BMD Analysis," and (2) "Functional Classifications." A "BMD Analysis" output file contains gene (or microarray probe), BMD and BMD lower confidence (BMDL) values for each statistical model, as well as the information required for model selection. A "Functional Classifications" output file can be exported as "Gene Ontology Analyses," "Signaling Pathway Analyses" or "Defined Category Analyses." This file includes statistical BMD values for each pathway, and lists the transcripts that were used to obtain the pathway BMD. These output files typically contain 30-60 columns and thousands of rows of numbers and text. Processing and analysis of these complex files can be challenging and time-consuming. Therefore, we have developed a tool that summarizes the BMDExpress output files, and allows for quick exploration, analysis and visualization of these complex datasets. Our tool allows users to identify important genes or biological processes based on their BMD values, which we anticipate will be very useful for human health risk assessment. In this paper, we describe the development of BMDExpress Data Viewer, which is a user-friendly web-based tool that reads BMDExpress output files, provides visual (graphical and tabular) data summaries and provides an interactive interface for the analyses. We anticipate that BMDExpress Data Viewer will streamline the analyses of BMDExpress output files and facilitate the use of transcriptional BMD data in human health risk assessment.

Materials and methods

Implementation, compatibility and availability

The core component of BMDExpress Data Viewer was written in the Java programming language, while the browser-based user interface was developed with a combination of hypertext markup language (HTML), Java Servlet and Javascript. Data visualization was created using the Google Charts application programming interface (API) (Google Inc., 2015). All visualization APIs, including histogram, scatter plot, pie chart, column chart, table and bubble chart, were provided by Google Developers, while the heatmap plot was provided by the Institute of Systems Biology (Institute of Systems Biology 2008). Javascripts to generate charts and tables are downloaded at the time the user requests BMDExpress Data Viewer to generate applicable visualizations, and the user's data are processed directly within the web browser.

The application is hosted on a Tomcat server, version 7.0.6, with Java version 1.6.0_37-b06. The Google Charts API is based on

HTML5 and Scalable Vector Graphics technology (Google Inc., 2015). The application is compatible with all modern browsers, including Internet Explorer 9+ (Microsoft Windows), Safari (Mac OSX), Firefox, Chrome and Opera.

BMDExpress Data Viewer is available as a browser-based application, and is freely available at BMDX_Viewer/ (MIT license).

BMD Analysis Summary: probe-to-gene mapping

To convert microarray probes to gene symbols, we first compiled a probe-to-gene symbol file by matching the probe2gene file to the genes2symbols file from each supported platform in the BMDExpress installation directory. A probe-to-accession number file was also compiled by obtaining the respective platform annotation files from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO). This allows the gene symbols that correspond to each probe to be displayed (if available) in the BMD Analysis Summary Dynamic Viewer. In the case where an official gene symbol is not available, corresponding GenBank accession numbers are displayed. In instances when neither gene symbol nor GenBank accession number are available, or if the imported file is not of a supported platform, the original probe identifiers are displayed. Currently, supported platforms include: Agilent human whole genome (4 x 44 K and 8 x 60 K), Agilent mouse whole genome (4 x 44K and 8 x 60 K), Agilent rat whole genome (4 x 44 K and 8 x 60 K), Affymetrix human (HG_Focus, HG_U133A, HG-U133A_2 and HG-U133_Plus_1), Affymetrix mouse (MG_U74A, MG_U74Avs, MOE430A, MOE320B, Mouse430A_2 and Mouse430_2) and Affymetrix Rat (RAE230A, RAE230B, Rat230_2 and RGU34A).

BMDExpress Data Viewer also supports human, mouse and rat RNAseq-based datasets that use Ensembl Gene ID to represent genes. The gene2ensembl file was downloaded from the NCBI's FTP site and was parsed for human (Taxonomy ID: 9606), mouse (Taxonomy ID: 10090) and rat genes (Taxonomy ID: 10116). The parsed file was matched by Entrez Gene ID to the Homo_sapiens. gene_info, Mus_musculus.gene_info and Rattus_norvegicus. gene_info files to compile an ensembl2symbols file. The BMD Analysis Summary Dynamic Viewer uses the ensembl2symbols file to display corresponding gene symbols. In the case where a gene symbol is not available, the Ensembl Gene ID is displayed.

Functional Enrichment Analysis: implementation of the Fisher's exact test

The Fisher's exact test in the Functional Enrichment Analysis tool uses Java Statistical Classes version 1.0 (Bertie, 2002). The Fisher's exact test examines the association between a set of genes to a particular pathway. For each pathway, a 2 x 2 contingency table is constructed. The Functional Enrichment Analysis tool defines the four elements of the contingency table as: (1) number of genes from the input list that are found in the biological process (or pathway, network, etc.) and show a dose response (i.e., genes with goodness-of-fit P > 0.1); (2) total number of genes in the pathway minus the number of significant genes in the pathways; (3) number of input genes minus the number of significant genes in the pathway; and (4) number of background genes minus the total number of genes in the pathway. The Fisher's exact test performed in the Functional Enrichment Analysis tool requires that the genes show a dose-response that can be modeled to be included in the analysis; thus, the approach is slightly different from other pathway enrichment analyses (e.g., the Database for Annotation,

Visualization and Integrated Discovery (DAVID); Huang et al., 2007b). The number of input genes is equal to the total number of genes surveyed by the technological platform, while the number of background genes is set to 30 000, which is the approximate total number of genes in human, mouse and rat.

Depending on the analysis approach, the user may apply pre-filtering (e.g., statistical presence or fold change) to the microarray probes before importing datasets to BMDExpress. If no pre-filtering or only the built-in ANOVA in BMDExpress is applied, all of the four Fisher's exact test parameters can be obtained from the two files uploaded to generate the Bubble Chart view (described below). In the case where pre-filtering other than the built-in ANOVA in BMDExpress was applied, the user needs to specify the platform and the pathway database, to determine the correct number of input genes and sizes of pathways.

Case study datasets

To demonstrate the functionality and features of BMDExpress Data Viewer, two microarray datasets were downloaded from GEO (, derived from two microarray platforms: Agilent mouse DNA microarrays and Affymetrix rat microarrays. In the first dataset (GSE48644), mice were exposed to increasing doses (0, 1, 2, 4, 8 mg kg-1 day-1 body weight) of the hepatocarcinogen furan (the liver is the target organ for cancer in rodents following furan exposure) by oral gavage for 21 days. Dose and tissue (liver) selection was based on a previous rodent cancer bioassay (Haseman et al., 1998), and included two doses that did not induce liver cancer in the same strain of mice (1 and 2 mg kg-1 day-1) and two doses that caused liver cancer (4 and 8 mg kg-1 day-1). Agilent DNA microarrays were used to examine hepatic transcriptional profiles for these mice. Detailed methods for the preparation of samples, microarray gene expression and normalization of the data are described in the publications (Jackson et al., 2014; Webster et al., 2014). In the second dataset (GSE45892), rats were exposed to increasing doses of one of six chemicals for 5,14, 28 and 90 days. Affymetrix HT RG-230 PM microarrays were used to examine transcriptional profiles for these rats (Thomas et al., 2013b). We obtained the datasets for the chemical 4,4'-methylenebis(N,N-dimethyl)benzenamine

(MDMB), which examined transcriptional changes in the thyroid following exposure to 50, 200, 375, 500 and 750 ppm in drinking water. All timepoints were considered. The chosen target tissue, strain, sex and route of exposure were based on the critical effect in the IRIS toxicological review for MDMB (US EPA, 1989). The MDMB datasets were normalized using the Robust Multi-array Average (Irizarry et al., 2003).

Benchmark dose analyses for the case study datasets

BMDExpress version 1.4.1 (Yang et al., 2007) was used to perform BMD analyses. For the Agilent dataset (furan study), only the probes that had signal intensities that were significantly above background probe (non-murine control oligonucleotides) intensities (at least three standard deviations above the mean of the background signal intensities) in at least one dose group were imported into BMDExpress for analysis. These datasets were pre-filtered using the built-in one-way ANOVA function for differential gene expression in three ways: (1) unfiltered; (2) ANOVA P < 0.05 (in at least one dose); and (3) false discovery rate (FDR) P < 0.05 (in at least one dose). Detailed analyses of these results were published previously (Webster et al., 2015). For the Affymetrix dataset

(MDMB study), all probes were imported into BMDExpress and the ANOVA P < 0.05 was applied. In the BMDExpress analysis for both furan and MDMB datasets, the best fitting model for each probe was selected based on: (1) a nested chi-squared cutoff value of 0.05 to select between linear and polynomial models; (2) lowest Akaike information criterion value for the Hill and Power models; and (3) likelihood ratio test goodness-of-fit P > 0.1. Other modeling parameters included maximum iterations of 250, confidence interval of 0.95, benchmark response value of 1.349 and restricting power parameter to > 1. A Hill model was flagged if the k parameter was greater than one-third of the lowest dose, and the next best model was selected if it had a goodness of fit P > 0.05. The Hill model was only used and the BMD was modified to 0.5 of the lowest BMD value if no other model had a P > 0.1. Using the "Defined Category Analyses" feature in BMDExpress, the BMD analyzed datasets were mapped to the Ingenuity Pathway Analysis (IPA) mouse (Furan) and rat (MDMB) canonical pathways, which were downloaded on April 24, 2014. In mapping to the IPA pathways, probes were removed if the BMD value was greater than the highest dose used in the experiment. Probes were also removed if the goodness-of-fit P<value was less than 0.1. Conflicting probe sets were identified and a correlation cutoff of 0.5 was applied to the conflicting probe sets. Both the BMD analyzed dataset and the IPA mapped dataset were exported using the "Export" function.


We developed BMDExpress Data Viewer, a web-based tool designed for visualization and exploration of BMDExpress output files. Our application provides graphical summaries and statistics of BMDExpress-generated output files. We demonstrate through two case studies the features and potential use cases of the application.

Overview of BMDExpress Data Viewer

BMDExpress Data Viewer consists of two collections of tools: "Summary Visualization Tools," and "Dataset Exploratory Tools" (Fig. 1). Each collection contains tools that summarize datasets in addition to allowing users to explore datasets and perform specific analyses. The workflow for BMDExpress Data Viewer begins with uploading unmodified BMDExpress output files, followed by interfaces to guide the user through to display results or perform analyses (Fig. 1).

Example analysis with furan-treated mouse liver and 4,4-methylenebis(N,W-dimethyl)benzenamine-treated rat thyroid datasets

To illustrate the features and functionality of BMDExpress Data Viewer, BMDExpress output datasets were obtained from two case studies, i.e., (1) hepatic gene expression (Agilent) profiles of mice orally exposed to increasing doses of furan by oral gavage (Jackson et al., 2014), and (2) thyroid gene expression (Affymetrix) profiles from rats that were orally exposed to increasing doses of the chemical MDMB (Thomas et al, 2013b).

Three statistical pre-filtering stringencies were applied to the furan data before BMD analysis: (1) no filtering (removing only the probes with signal intensities within the background range across all conditions, i.e., genes with low or no expression); (2) genes that reached an ANOVA P < 0.05 compared to controls in at least one dose group; and (3) genes that reached an FDR P < 0.05 in at least one dose group compared to controls. Overall, the © 2015 Her Majesty the Queen in Right of Canada. Journal of Applied Toxicology J. Appl. Toxicol. 2016; 36:1048-1059

published by John Wiley & Sons, Ltd.

Dataset Summary

• Inspect probe & pathway BMD(L) distributions

• Inspect probe fit p-value distributions

• Inspect model distributions

• Inspect ratios of BMD vs. BMDL

• Identify noise and outliers

• Identify modes & predict POD BMD(L)s

Dynamic Viewer

• Inspect probe or pathways within specific regions (e.g. peaks or noise) with zoomable slider bars

Functional Enrichment Analysis

• Visualize pathways (or GO) as bubbles

• Filter with Fisher's Exact Test p-value, pathway size, # of significant genes & BMD(L).

• Select sensitive pathways for further investigations

Multiple Dataset Comparison

• Compare pathway BMD(L)s across multiple experimental conditions (e.g. time points, tissues & platforms)

• Visualize trends across multiple datasets

BMDL-BMD Ranee Plotter

Visualize trends of BMDL/BMD ranges in selected pathways Identify pathways with low and/ or high BMD(L) values Predict POD BMD(L)s

Figure 1. Workflow demonstrating the features and functionality of BMDExpress Data Viewer. BMD, benchmark dose; BMDL, benchmark dose lower confidence (values); GO, gene ontology; POD, point of departure.

filtered datasets contained 29847, 2597 and 362 probes, for the un-filtered, ANOVA-filtered and FDR-filtered datasets, respectively. For the MDMB study, we applied an ANOVA P < 0.05 filter on all four timepoints before BMD analysis. The resulting dataset contained 6006, 7130, 8981 and 2697 probes for the 5-, 14-, 28- and 90-day timepoints, respectively.

To demonstrate the features of each tool, we compared the effects of statistical pre-filtering in the furan study, as well as changes in the BMD over the four timepoints in the MDMB study. We analyzed the summaries of probes and pathways for the three pre-filtered furan datasets and the four timepoints in the MDMB study. We applied both studies to the "Functional Enrichment Analysis" tool to identify potentially sensitive pathways. The pathways that were common across timepoints in the MDMB study were further examined in the "Multiple Dataset Comparison" and "BMDL-BMD Range Plotter" tools.

BMD Analysis Summary

The BMD Analysis Summary tool was applied to provide visual summaries of the BMD, BMDL and fit P-values, and the best-fit models, which are the commonly examined parts of the probe BMD datasets (Fig. 2). Two components, Dataset Summary and Dynamic Viewer, are available in the BMD Analysis Summary tool. The dataset summary component presents overall summaries of the BMD and BMDL values as histograms and scatterplots to illustrate distributions, enabling inspection of main trends and outliers. A histogram of fit P-values is also presented for inspection of the overall quality of modeling. Summary statistics are also computed for maximum, minimum, mean and median BMD and BMDL values, as well as fit P-values. Pie charts are plotted to illustrate the distribution of best-fit models selected by BMDExpress. In the Dynamic Viewer component, the complete

Figure 2. Summaries of the analyses performed using the "BMD Analysis Summary" tool. (A) Effects of pre-filtering methods on the distributions of probe BMD values (furan dataset on the Agilent platform). (B) Distributions of probe BMD values over timepoints (MDMB dataset on the Affymetrix platform). (C) Distributions of fit P-values (MDMB dataset on the Affymetrix platform). (D) BMD vs. BMDL over time (MDMB dataset on the Affymetrix platform). (E) Distribution of selected models over time (MDMB dataset on the Affymetrix platform). BMD, benchmark dose; BMDL, benchmark dose lower confidence (values); FDR, false discovery rate; MDMB, 4,4'-methylenebis(N,N-dimethyl)benzenamine.

list of probes uploaded is presented in both a stepped chart and a table, illustrating BMD and BMDL values for each probe. Slider bars for BMD and BMDL values allow the user to set ranges to focus on specific parts of the dataset.

We used the BMD Analysis Summary tool to investigate the distributions of BMD and BMDL values in both datasets. As illustrated in Fig. 2(A), the BMD histogram of the unfiltered furan dataset showed an even distribution of BMDs across dose. In contrast, peaks of BMD values were clearly visible in both the ANOVA and FDR-filtered datasets. The even distribution in the unfiltered dataset and the bimodal distribution in the ANOVA-filtered dataset were reduced to approach a unimodal distribution in the FDR-filtered dataset. The main mode (peak) in the

ANOVA-filtered dataset occurred at approximately 4.2 mg kg-1 day-1, while for the FDR-filtered dataset, it was approximately 3 mg kg-1 day-1. In the MDMB study (Fig. 2B), ANOVA pre-filtering was applied and all four timepoints displayed clearly visible peaks, with unimodal distributions at 14 and 28 days, and bimodal distributions at 5 and 90 days. The main modes occurred at approximately 440, 190, 450 and 460 ppm for 5, 14, 28 and 90 days, respectively. Fit P-value distributions of the MDMB datasets were examined (Fig. 2C). No clear peaks were observed for any of the four timepoints. However, the distribution showed a slightly increasing trend towards higher fit P-values at 5 days, while decreasing trends were observed at 14, 28 and 90 days. © 2015 Her Majesty the Queen in Right of Canada. Journal of Applied Toxicology J. Appl. Toxicol. 2016; 36:1048-1059

published by John Wiley & Sons, Ltd.

We also examined the BMD versus BMDL scatterplots for the MDMB dataset (Fig. 2D). The BMD/BMDL ratios were approximately 1.2 with few outliers (two at 14 days and one at 90 days). The pie chart of models suggested that the Power model was the most frequently applied (36.7-44.7%), followed by the Hill model (27.2-31.7%) in 5-, 28- and 90-day datasets (Fig. 2E). For 14 days, the 2° Polynomial model was the predominant best-fit model at 55.9%, followed by the Hill model at 17.8%.

Functional Classification Summary

The Functional Classification Summary tool visually summarizes the BMD means, BMD medians, BMDL means, BMDL medians and BMD 5th and 10th percentiles from the export files of Functional Classification Analyses, which may be one of: Signaling Pathway Analyses, Gene Ontology Analyses or Defined Category Analyses. Similar to the BMD Analysis Summary tool, two components are available: Dataset Summary and Dynamic Viewer. In the dataset summary component, histograms are plotted to illustrate the distributions. Scatterplots of BMD mean versus BMDL mean, BMD median versus BMDL median and BMD 5th

percentile versus BMD 10th percentile illustrate the relationships of these values for each pathway. In the Dynamic Viewer component, a complete list of pathways and their respective BMD and BMDL means and medians, BMD 5th and 10th percentiles, as well as the number of significant genes and the percentage of these genes relative to the total number of genes in the pathway are also plotted in a stacked chart and a data table. Slider bars for these values are available for setting ranges to focus on the desired part of the dataset. By default, the number of significant genes is "Items with BMD P > 0.1" column, where the value 0.1 can be any cutoff the user used in BMDExpress analysis. This ensures only genes with the BMD less or equal to the highest dose and with a P-value greater or equal to the cutoff for the model fit is examined.

Using the Functional Classification Summary tool, we investigated the distributions of BMD/BMDL mean and median values, as well as BMD 5th and 10th percentile values in both the furan and the MDMB studies. We highlight the results of BMD mean distributions in both studies (Fig. 3A,B). In both studies, BMD means for the pathways displayed unimodal distributions, with minimal noise on either side of the histograms. For the furan study, the

Figure 3. Summaries of the analyses performed using the "Functional Classification Summary" tool. (A) Effects of pre-filtering methods on the distributions of pathway BMD mean values (furan dataset on the Agilent platform). (B) Distributions of pathway BMD mean values over time (MDMB dataset on the Affymetrix platform). (C) Twenty-four pathways identified by the Dynamic Viewer in the peak regions of the MDMB 90-day timepoint. Slider bar ranges for BMDL mean, median, BMD fifth and 10th percentile values were selected based on the main modes of their respective histograms. A minimum of five genes was applied as suggested in previous studies (Moffat et al., 2015; Thomas et al., 2011,2012). BMD, benchmark dose; FDR, false discovery rate; MDMB, 4,4'-methylenebis(N,N-dimethyl)benzenamine.

B. Kuo et ai.

modes were approximately 4.5, 4 and 2.5 mg kg-1 day-1 for the unfiltered, ANOVA-filtered and FDR-filtered datasets, respectively. For the MDMB study, these numbers were approximately 375, 315, 435 and 440 ppm for the 5-, 14-, 28- and 90-day timepoints, respectively. Compared to the bimodal distributions in the ANOVA-filtered furan dataset observed at the gene level, as well as in the 5- and 90-day MDMB timepoints, the secondary modes observed at the lower doses had largely been reduced.

Using the Dynamic Viewer, we investigated the pathways that were at the modes of the MDMB 90-day histograms (see Supplementary Fig. S1 for all histograms). By setting the sliders on top of the Dynamic Viewer to where the modes of each histogram occurred, with a minimum of five genes, 24 pathways were identified (Fig. 3C). The majority of these pathways had more than 10 significant genes at approximately 15% mapping percentage.

Functional Enrichment Analysis

Although some studies have suggested that application of BMD gene expression data focus on selecting biological processes, such as gene ontology terms and pathways, which are associated with disease processes or known mechanisms of action in toxicity, it is valuable to consider more unbiased approaches to filter through and explore pathways of potential toxicological significance. The Functional Enrichment Analysis tool displays BMDs for biological processes in a bubble chart view to enable exploration of the dataset and selection of perturbations that may have significant biological impacts.

In the Functional Enrichment Analysis tool, each biological process is represented as a bubble in a two-dimensional space, with enrichment probability on the Y-axis and BMD value on the X-axis (Fig. 4). Enrichment probability is computed using a Fisher's exact test. The BMDs that are shown can be BMD mean, BMD median, BMDL mean, BMDL median or BMD 5th or 10th percentiles. In addition to position, a bubble is also defined by its size and color.

The size of a bubble represents the total number of genes in the biological process, whereas the depth of color indicates the number of genes from the user's input list that mapped to the pathway (i.e., more affected genes in a pathway will mean a darker color). Slider bars are available to set ranges to limit the number of significant genes, size of pathways and Fisher's exact test P-value in percentages. Altering the ranges of these values dynamically refreshes the bubble chart and the data table by displaying only the pathways fitting the criterion. The user can mouse-over to the pathway (bubble) of interest and a dialogue box will automatically appear to show the respective information on the number of genes, pathway size, Fisher's exact test P-value and BMD.

To demonstrate the ability to filter for sensitive pathways, we applied the FDR-filtered dataset of the furan study to the Functional Enrichment Analysis tool and compared the resulting pathways to the pathways identified previously using IPA (Jackson et ai., 2014). By applying the Functional Enrichment Analysis tool with Fisher's exact test P < 0.05 (5%), and the minimum number of significant genes equal to three (in accordance with Jackson etai. [2014]), we obtained 24 pathways (Table 1). Fourteen (58%) of the identified pathways were also found in our previous IPA analysis (Jackson etai., 2014). The mean and median values for the BMDL means for these pathways were 1.77 and 1.63 mg kg-1 day-1, respectively. The results show that by applying a cut-off with the Fisher's exact test P-value, the Functional Enrichment Analysis tool is capable of identifying known pathways that are potentially important for the mode of action of furan, and that exhibits a dose-response.

We also applied the Functional Enrichment Analysis tool to the four timepoints of the MDMB study. We applied a Fisher's exact test P < 0.01 (1%) and the minimum number of significant genes equal to five (in accordance with Thomas et ai., 2013b) to all four timepoints. For each of the 5-, 14-, 28- and 90-day timepoints, we obtained 91, 67, 242 and 44 pathways, respectively. Comparison

Figure 4. Screenshot of the analysis on the 90-day timepoint of the MDMB study using the "Functional Enrichment Analysis" tool. Each pathway is represented as a bubble, with size of bubble representing the total number of genes in the biological process, and the depth of color indicating the number of genes from the user's input list that mapped to the pathway. A minimum of five genes and a Fisher's exact test P-value cutoff of 5% (0.05) were applied in our example analyses. BMD, benchmark dose; MDMB, 4,4'-methylenebis(N,N-dimethyl)benzenamine.

J. Appi. Toxicoi. 2016; 36:1048-1059

Table 1. Pathways identified by the Functional Enrichment Analysis tool for the ANOVA-filtered furan-treated mouse liver dataset. A minimum of three significant genes and a Fisher's exact test P < 0.05 (5%) were applied for the analysis

Pathway BMDL mean Fisher's exact test Number of

(mg kg-1 day-1) P-value (%) significant genes Pathway size

14-3-3-mediated signalinga 1.63 1.04 5 110

ATM signaling 1.67 0.34 4 52

Aldosterone signaling in epithelial cellsa 1.40 0.78 6 144

Aryl hydrocarbon receptor signalinga 2.15 0.07 7 124

BMP signaling pathway 2.01 4.95 3 69

Bladder cancer signaling 2.43 1.42 4 78

Breast cancer regulation by stathmin1a 1.67 1.98 6 177

Germ cell-Sertoli cell junction signalinga 1.42 3.10 5 146

Glioblastoma multiforme signaling 1.31 2.80 5 142

Glucocorticoid receptor signalinga 1.41 0.26 9 242

Hypoxia signaling in the cardiovascular systema 1.67 0 7 57

NGF signaling 1.98 3.50 4 103

NRF2-mediated oxidative stress responsea 2.58 0.02 9 165

PI3K signaling in B lymphocytes 2.18 0.06 7 118

PI3K/AKT signalinga 1.52 0.27 6 116

PPAR signalinga 1.30 1.68 4 82

PTEN signaling 1.88 4.91 4 115

PXR/RXR activationa 3.06 1.55 3 44

Prostate cancer signaling 1.36 1.30 4 76

Protein ubiquitination pathwaya 1.61 0.06 10 239

Remodeling of epithelial adherens junctions 1.54 3.49 3 60

Retinol biosynthesis 1.50 0.36 3 26

Sertoli cell-Sertoli cell junction signalinga 1.74 1.16 6 157

p53 Signalinga 1.64 0.34 5 84

BMDL, benchmark dose lower confidence (values).

aPathways that matched to the Ingenuity Pathway Analysis (IPA) result performed in a previous study (Jackson etal, 2014).

of these four lists of pathways showed that nine were in common, with means (medians) of 351.50 (352.16), 350.27 (340.27), 404.50 (404.67) and 415.48 (417.12) ppm for 5,14,28 and 90 days, respectively (Table 2).

Multiple Dataset Comparison

The Multiple Dataset Comparison tool lets the user compare how BMD values vary over different experimental conditions, such as

across timepoints, tissues or treatments. This tool lets the user select the number of datasets to compare, followed by uploading the files accordingly. It is important that the files uploaded in a comparison analysis are analyzed using the same pathway or defined category source, as the tool looks for exact matches in names between files and generates a list of biological processes to be selected for comparison. The tool displays only the biological processes common to the uploaded files for selection. When only one file is selected and uploaded, the user can compare pathway

Table 2. Common pathways identified using the Functional Enrichment Analysis tool on the four timepoints of the MDMB study (values for each timepoint are BMD means) Pathway 5 days 14 days 28 days 90 days

Colorectal cancer metastasis signaling 332.01 342.84 394.15 396.45

Endothelin-1 signaling 358.24 325.90 410.97 419.76

Glioblastoma multiforme signaling 344.00 321.17 404.83 396.32

HMGB1 signaling 354.35 360.03 404.39 423.17

Hepatic fibrosis/hepatic stellate cell activation 352.81 432.79 421.11 457.97

IL-8 signaling 343.88 337.70 387.05 418.76

ILK signaling 358.97 379.46 408.96 387.71

Nitric oxide signaling in the cardiovascular system 348.94 337.21 410.81 379.97

Regulation of the epithelial-mesenchymal transition pathway 370.30 315.33 398.27 459.18

Mean 351.50 350.27 404.50 415.48

Median 352.16 340.27 404.67 417.12

BMD, benchmark dose; IL, interleukin; MDMB, 4,4'-methylenebis(N,N-dimethyl)benzenamine.

O Ui U1

B. Kuo et al.

Figure 5. Comparison of BMD mean values for the nine common pathways identified across four timepoints in the MDMB study using the "Multiple Dataset Comparison" tool. BMD, benchmark dose; IL, interleukin; MDMB, 4,4'-methylenebis(N,N-dimethyl)benzenamine.

BMDs within the same dataset. BMD comparisons can be performed on BMD mean, BMD median, BMDL mean, BMDL median and BMD fifth and 10th percentiles. Four visualization options are available for comparisons: column chart, stepped chart, heat map and data table.

To demonstrate the Multiple Dataset Comparison tool, we compared nine pathways that we found were enriched using the Functional Enrichment Analysis tool in each of the four timepoints for the MDMB study (Table 2 and Fig. 5). In general, BMD mean values were considerably lower at the 5- (average of 352 ppm)and 14-day (average of 350 ppm) timepoints than at the 28- (average of 405 ppm) and 90-day (average of 415 ppm) timepoints. All but the giio-biastoma muitiforme signaiing, ILKsignaiing and nitric oxide signai-ing in the cardiovascular system pathways showed increasing BMD mean values from 28 to 90 days. For the 14-day timepoint, the hepatic fibrosis/hepatic steiiate ceii activation pathway had a BMD mean value of 433 ppm, while the other eight pathways were within 310-380 ppm. Overall, the BMD mean values for these pathways displayed an increasing trend over time.

BMDL-BMD Range Plotter

The BMDL-BMD Range Plotter tool displays both BMD and BMDL mean or median values as upper and lower bound ranges for user-selected pathways in a single figure, similar to approaches

described previously (Slob and Setzer, 2014). As BMDExpress does not provide BMD upper confidence limit (BMDU) values at this time, we applied BMD as the upper bound in the figure at this time. We define the value bound between the BMDL and the BMD as a "range." The tool displays the range of BMDs for the selected pathways, enabling assessment of these ranges across pathways, and providing information on trends that are useful in determining what BMD may be relevant for regulatory purposes.

To demonstrate the BMDL-BMD Range Plotter tool, we selected the nine pathways identified by the Functional Enrichment Analysis tool that were common between the four timepoints of the MDMB study. We applied the tool to the datasets for each of these timepoints. The analysis revealed consistent ranges between the BMDL and BMD (approximately 110 ppm) for all timepoints (Fig. 6). In each of the 5-, 28- and 90-day timepoints, the nine pathways overlapped with each other, and the lengths of overlaps were approximately 60-70 ppm. In contrast, the ranges of all of the pathways did not overlap with each other in the 14-day timepoint. In general, the ranges extended to lower doses for the earlier timepoints (BMDLs approximating 250-350 ppm, and 220-350 ppm, for 5 and 14 days, respectively), while they were bound at higher doses for the later timepoints (approximately 300-400 ppm, and 250-450 ppm, for 28 and 90 days, respectively). For the 14-day timepoint, the range of the hepatic fibrosis/hepatic steiiate ceii activation pathway had a higher BMDL and BMD

Figure 6. BMDL-BMD range plots for the nine common pathways identified across four timepoints in the MDMB study using the "BMDL-BMD Range Plotter" tool. BMD, benchmark dose; BMDL, benchmark dose lower confidence (values); MDMB, 4,4'-methylenebis(N,N-dimethyl)benzenamine.

J. Appl. Toxicol. 2016; 36:1048-1059

bound, while the other eight pathways extended to lower doses. In general, the centers of overlaps were lower (290 and 300 ppm) for 5 and 14 days, respectively, and higher (340 and 360 ppm) for 28 and 90 days, respectively.


We developed BMDExpress Data Viewer to assist in the analysis of BMDExpress output files by providing visualization tools to summarize datasets, and data exploratory tools to identify sensitive biological processes, to compare BMDs across experiments and to determine BMDs of potential regulatory use. In this paper, we analyzed DNA microarray datasets from rodent tissues generated from two separate dose-response experiments to demonstrate the features and functionality of BMDExpress Data Viewer. We show that BMDExpress Data Viewer can be used to visualize BMDExpress output files and provide novel ways to explore and analyze these datasets.

By applying the furan and MDMB datasets to the BMD Analysis Summary and the Functional Classification Summary tools, we demonstrate some variables that may influence the BMDs derived from whole genome transcriptional profiles. For example, statistical pre-filtering of the data influences the distribution of gene BMDs leading to more normal distributions that may yield more robust BMD estimates. In addition, a general trend of declining BMDs with increasing statistical stringency (i.e., unfiltered versus ANOVA-filtered versus FDR-filtered) was found. However, secondary modes and the "noise" at extreme ends of the probe BMD histograms (Fig. 2A,B) were largely minimized in the histograms of pathway BMDs (Fig. 3A,B). Our results also showed that the peaks (main modes) of BMD were very similar to the BMD values reported in previous studies. It has previously been proposed that the lowest BMD be used as a point of departure for risk assessment (Moffat et ai., 2015; Thomas et ai., 2011, 2012). Our results suggest that the main modes of gene expression BMDs may be useful surrogates to estimate the BMDs at which adverse apical effects (such as cancer in the 2-year cancer bioassay) may occur (Webster et ai., 2015). More work is needed across a diverse array of toxicological endpoints to determine the most effective application of gene expression BMDs as points of departure.

In the furan ANOVA dataset, as well as in the MDMB 5- and 90-day datasets, the BMDs were bimodally distributed. As an example of how this distribution might be used to inform risk assessment, we extracted the genes underlying each mode from the three BMDExpress output files and performed pathway enrichment analysis using IPA. There was no overlap in the top 10 pathways from these analyses, with a single exception (mitotic roies of poio-iike kinase pathway in the MDMB 5-day dataset) (Supplementary file). Given that the pathways in the lower mode are activated by a lower dose than the pathways in the higher mode, these processes may represent the most sensitive drivers of the underlying toxicological effects. We speculate that these processes will be important considerations for establishing the mode of action and point of departure for this chemical. However, detailed analysis and interpretation of the biological effects observed here are beyond the scope of this paper. More work is needed to investigate further the mechanisms of interactions between the groups of pathways.

The histogram of fit P-value distributions provides an opportunity to investigate the overall quality of the models. Our results from the four MDMB timepoints suggest an overall lower confidence of fit for the 14-, 28- and 90-day datasets (fit P-values closer

to cutoff of 0.1) relative to the 5-day dataset where more probes were modeled with higher confidence (fit P-values closer to 1). To demonstrate an example of further analysis, we obtained the number of probes before and after the fit P-value filter (by checking and unchecking the option in the BMD Analysis Summary tool) and computed the percentage of probes that could be modeled (Supplementary File 2). For all timepoints, except for 14 days, over 80% of the probes that passed the initial filter of one dose with ANOVA P < 0.05 relative to control could be modeled. In contrast, only ~50% of probes could be modeled at 14 days. Interestingly, the 14-day timepoint also exhibited a different model distribution from the other three timepoints (Fig. 2E). More research is needed to understand the biological basis for these differences, but this is beyond the scope of the present study.

The scatterplot of BMD versus BMDL provides the capability to identify quickly the outliers and allows the user to investigate ratios of BMD/BMDL. We note that it has been suggested that since BMDL calculation takes into account measures of data uncertainty, a higher BMD/BMDL ratio may indicate higher uncertainty; thus, it has been proposed that a BMD/BMDL ratio < 2 indicates low uncertainty (Muri etai., 2009). BMDExpress Data Viewer enables users to assess rapidly which genes and pathways are not in agreement with these principles, facilitating the filtering of poor fitting data before subsequent analyses.

The pie chart of models in the BMD Analysis Summary tool displays an overall representation of the models used. Where a specific model is preferred, such as the Hill model for receptor-mediated responses (Zhao et ai., 2010), the pie chart can be used for a quick assessment of the suitability of the models used. If an unexpected skew is observed in the types of models used, the model selection parameters in BMDExpress may then be adjusted. In our analysis with the four MDMB timepoint datasets, we noted a considerably different distribution of models in the 14-day timepoint in comparison to other timepoints (Fig. 2E). The dissimilar model distribution in the 14-day timepoint may be the cause for some divergent results observed in the main modes of probe and pathway BMD histograms at this timepoint, as well as in the reduced degree of overlap in the BMDL-BMD range plots. More work is needed to assess the source of this observed dissimilarity.

The Dynamic Viewers in the BMD Analysis Summary and Functional Classification Summary tools provide an opportunity for the user to explore graphically the input datasets with the flexibility to zoom in and out for more or less detail. If the user has an expected range of BMDs in mind, or if pathways within a particular range are of interest, the Dynamic Viewers will reveal the genes or pathways in the set ranges. In addition, the Dynamic Viewers can be used to examine the specific ranges of BMDs within a region (e.g., within a mode, or within an outlier regions) identified in the histograms or scatterplots. Overall, this viewer enables a dynamic filtering capability for comparison of BMD values across pathways. In our investigation of the modes of the 90-day timepoint of the MDMB study, we identified 24 pathways (Fig. 3C); of these, only the endotheiin-1 signaiing, HMGB1 signaiing, IL-8 sig-naiing and RAR activation pathways were among the nine common pathways later identified using the Functional Enrichment Analysis tool. This low number of matches is not unexpected, as the Dynamic Viewer simply presents the pathways within the specified BMD ranges without utilizing probability values of the statistical enrichment test.

The Functional Enrichment Analysis tool represents biological processes as bubbles in a two-dimensional space. The main

difference between this analysis and other pathway or functional enrichment approaches is that it only includes genes that were fit to a mathematical model and therefore have a BMD. Thus, the test requires that the gene was significant in at least one dose, in addition to showing a dose-response that could be modeled. Overall, it makes use of dose-response information that is not considered in the other pathway enrichment approaches. With this tool, users can visualize and evaluate the relevance of each pathway. Larger bubbles represent larger pathways (pathways with more genes) and the darker the color, the more genes perturbed in the pathway. Thus, large and dark bubbles that are near the bottom left (i.e., lower BMDs and P-values) may be important biological perturbations. Applications such as DAVID have integrated a Fisher's exact test to determine the probability of a pathway matching from a short list of genes as opposed to randomly picking it up from a full list of genes in the organism (Huang et al, 2007a). We found that by applying a Fisher's exact test P < 0.05 and a minimum of three significant genes, similar pathways were identified to our previous IPA approach (Jackson et al. 2014) (Table 1). However, our approach here has the added benefit of including the BMDs within the pathway analysis. As we decreased the Fisher's exact test P-value and increased the minimum number of significant genes in the Functional Enrichment Analysis tool, more pathways in the upper range of the BMDL mean axis were filtered out. The remaining pathways that matched the previous IPA analysis were generally larger. This shows how the Functional Enrichment Analysis tool presented here is capable of identifying potentially important pathways for further investigations.

By using the Multiple Dataset Comparison tool, the nine common pathways in the MDMB study showed an increasing BMD trend over time (Table 2 and Fig. 5). A potential hypothesis is that over time the organism may be adapting to the chemical exposure, hence requiring a higher dose to activate these pathways. We also have applied this tool to compare datasets across three gene expression technology platforms: reverse transcription-polymerase chain reaction, microarray and RNAseq (Webster et al., 2015). Overall, the "Multiple Dataset Comparison" tool enables an analysis of important trends across different datasets to provide support in identifying the most relevant pathways.

Application of the BMDL-BMD Range Plotter tool to nine common pathways between the four MDMB timepoints revealed lower BMD values at two earlier timepoints than the later timepoints (Fig. 6). This observation was consistent with results of the Multiple Dataset Comparison analysis. The center values of overlaps for the timepoints were very similar to the 2-year rodent cancer bioassay BMD of 381 ppm (Thomas et al., 2013b). The higher BMDL-BMD range of the hepatic fibrosis/hepatic stellate cell activation pathway in 14, 28 and 90 days may indicate that the pathway requires a higher dose to be activated than the other pathways. Slob and Setzer (2014) demonstrated that within a BMDL-BMDU range, a true BMD can be defined with 90% confidence and thus true differences across pathway BMDs can be assessed. However, in the BMDL-BMD Range Plotter, similar confidence cannot be assumed. When two BMDL-BMD ranges overlap, their true BMD values are more likely similar if the overlapping interval is larger. Such confidence will increase if the ranges include BMDU values. Non-overlapping BMDL-BMD ranges, in addition, may or may not overlap with BMDU values. BMDExpress does not compute BMDU. Enabling this feature would improve the utility of this tool and will be addressed in future work. At present, the BMDL-BMD Range

Plotter provides an effective means to visualize trends and assess potential pathways for use in determining a transcriptional point of departure.

In conclusion, BMDExpress is an important software application used to compute BMDs from high-throughput transcriptional datasets. We anticipate that the application of such approaches will grow as higher-throughput and less expensive methodologies for global gene expression analysis are produced. Indeed, it has been proposed that BMDs from toxicogenomic experiments be used as points of departure in human health risk assessment in the absence of data from traditional methods (Chepelev et al. 2015; Moffat et al. 2015; Thomas et al. 2012, 2013a, b). The application of toxicogenomics in risk assessment is still in an early phase of development. We hope BMDExpress Data Viewer can help advance the utility of toxicogenomics in developing effective approaches to determine the most suitable BMDs in a whole genome analysis. Previous work has shown that toxicogenomic experiments with sufficient biological replicates and doses (preferably encompassing doses both above and below the BMD) can be modeled to derive the BMDs that are comparable to apical BMDs (Black et al., 2014; Webster et al., 2015). We highly recommend that the designs of experiments in this field follow best practices that are already established for BMD modeling of apical data, for example, having at least one dose near the anticipated BMD as well as including the lower doses that should be closer to point of departure for sensitive endpoints (US EPA, 2012). We demonstrate that BMDExpress Data Viewer can be used to assist in the interpretation of toxicogenomic BMDs derived using BMDExpress, and hope that it will facilitate further work to establish the best approaches for BMD applications to toxicogenomic data. The application provides a means to evaluate the quality of study data, explore BMD ranges across pathways, experiments and experimental conditions, and identify the most relevant perturbations. In the example analyses, we showed how this tool could be used to explore the most sensitive biological perturbations occurring in a transcriptomics study through analysis of gene and pathway BMDs. These features and results illustrate that BMDExpress Data Viewer is a useful tool to visualize, explore and analyze BMDExpress output files.


We thank Jill Franzosa, Reza Farmahin, Nikolai Chepelev, John Wills and Sarah Labib for providing useful comments and suggestions. We thank Reza Farmahin for identifying the MDMB datasets to be included in the analysis. We thank Andrew Williams for providing useful advice on the statistical methods. We acknowledge Scott Auerbach of the National Toxicology Program along with Ruchir Shah, Dan Svoboda and Jason Phillips of Sciome LLC for technical support. This research was funded by the Health Canada Genomics Research and Development Initiative and the Chemicals Management Plan. AFW was supported by the Natural Sciences and Engineering Research Council and the Ontario Graduate Scholarship.

Conflict of interest

The authors did not report any conflict of interest.


Bertie AJ. 2002. Java applications for teaching statistics. MSOR Connect 2: ui-81. © 2015 Her Majesty the Queen in Right of Canada. Journal of Applied Toxicology J. Appl. Toxicol. 2016; 36:1048-1059

published by John Wiley & Sons, Ltd.

Bhat VS, Hester SD, Nesnow S, Eastmond DA. 2013. Concordance of tran-scriptional and apical benchmark dose levels for conazole-induced liver effects in mice. Toxicol. Sci. 136: 205-215.

Black MB, Parks BB, Pluta L, Chu TM, Allen BC, Wolfinger RD, Thomas RS. 2014. Comparison of microarrays and RNA-seq for gene expression analyses of dose-response experiments. Toxicol. Sci. 137:385-403.

Bourdon JA, Williams A, Kuo B, Moffat I, White PA, Halappanavar S, Vogel U, Wallin H, YaukCL. 2013. Gene expression profiling to identify potentially relevant disease outcomes and support human health risk assessment for carbon black nanoparticle exposure. Toxicology 303: 83-93.

Budtz-Jorgensen E, Bellinger D, Lanphear B, Grandjean P, International Pooled Lead Study Investigators. 2013. An international pooled analysis for obtaining a benchmark dose for environmental lead exposure in children. Risk Anal. 33: 450-461.

Chepelev NL, Moffat ID, Labib S, Bourdon-Lacombe J, Kuo B, Buick JK, Lemieux F, Malik AI, Halappanavar S, Williams A, YaukCL. 2015. Integrating toxicogenomics into human health risk assessment: lessons learned from the benzo[a]pyrene case study. Crit. Rev. Toxicol. 45:44-52.

Crump KS. 1984. A new method for determining allowable daily intakes. Fundam. Appl. Toxicol. 4: 854-871.

Google Inc. 2015. Google Charts - Google Developers. https://developers. (25 March 2015).

Haseman JK, Hailey JR, Morris RW. 1998. Spontaneous neoplasm incidences in Fischer 344 rats and B6C3F1 mice in two-year carcinogenicity studies: a National Toxicology Program update. Toxicol. Pathol. 26:428-441.

Health Canada. 2013. Final Human Health State of the Science Report on Lead. Government of Canada.

Huang DW, Sherman BT, Tan Q, Collins JR, Alvord WG, Roayaei J, Stephens R, Baseler MW, Lane HC, Lempicki RA. 2007a. The DAVID Gene Functional Classification Tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 8: R183.

Huang DW, Sherman BT, Tan Q, Kir J, Liu D, Bryant D, Guo Y, Stephens R, Baseler MW, Lane HC, Lempicki RA. 2007b. DAVID Bioinformatics Resources: expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res. 35: W169-175.

Institute of Systems Biology.2008. Visualization: BioHeatMap - A Heatmap for Gene Expression and Other Data http://informatics.systemsbiology. net/visualizations/heatmap/bioheatmap.html (25 March 2015).

Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. 2003. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics 4: 249-264.

Jackson AF, Williams A, Recio L, Waters MD, Lambert IB, Yauk CL. 2014. Case study on the utility of hepatic global gene expression profiling in the risk assessment of the carcinogen furan. Toxicol. Appl. Pharmacol. 274:63-77.

Moffat I, Chepelev NL, Labib S, Bourdon-Lacombe J, Kuo B, Buick JK, Lemieux F, Williams A, Halappanavar S, Malik AI, Luijten M, Aubrecht J, Hyduke DR, FornaceAJJ,SwartzCD, Recio L,YaukCL. 2015.Comparison of toxicogenomics and traditional approaches to inform mode of action and points of departure in human health risk assessment of benzo[a]pyrene in drinking water. Crit. Rev. Toxicol. 45:1-43.

Muri SD, Schlatter JR, Bruschweiler BJ. 2009. The benchmark dose approach in food risk assessment: is it applicable and worthwhile? Food Chem. Toxicol. 47: 2906-2925.

Slob W, Setzer W. 2014. Shape and steepness of toxicological dose-response relationships of continuous endpoints. Crit. Rev. Toxicol. 44: 270-297.

Thomas RS, O'Connell TM, Pluta L, Wolfinger RD, Yang L, Page TJ. 2007. A comparison of transcriptomic and metabonomic technologies for identifying biomarkers predictive of two-year rodent cancer bioassays. Toxicol. Sci. 96: 40-46.

Thomas RS, Clewell HJ3, Allen BC, Wesselkamper SC, Wang NC, Lambert JC, Hess-Wilson JK, Zhao QJ, Andersen ME. 2011. Application of transcriptional benchmark dose values in quantitative cancer and noncancer risk assessment. Toxicol. Sci. 120:194-205.

Thomas RS, Clewell HJ3, Allen BC, Yang L, Healy E, Andersen ME. 2012. Integrating pathway-based transcriptomic data into quantitative chemical risk assessment: a five chemical case study. Mutat. Res. 746:135-143.

Thomas RS, Philbert MA, Auerbach SS, Wetmore BA, Devito MJ, Cote I, Rowlands JC, Whelan MP, Hays SM, Andersen ME, Meek ME, Reiter LW, Lambert JC, Clewell HJ3, Stephens ML, Zhao QJ, Flowers L, Carney EW, Pastoor TP, Petersen DD, Yauk CL, Nong A. 2013a. Incorporating new technologies into toxicity testing and risk assessment: moving from 21st century vision to a data-driven framework. Toxicol. Sci. 136: 4-18.

Thomas RS, Wesselkamper SC, Wang NC, Zhao QJ, Petersen DD, Lambert JC, Cote I, Yang L, Healy E, Black MB, Clewell HJ3, Allen BC, Andersen ME. 2013b. Temporal concordance between apical and transcriptional points of departure for chemical risk assessment. Toxicol. Sci. 134: 180-194.

US EPA. 1989. Integrated Risk Information System (IRIS) Chemical Assessment Summary: 4,4'-Methylene bis(N,N'-dimethyl)aniline (CASRN 10161-1). 0386_summary.pdf (5 October 2015).

US EPA. 2012. Benchmark Dose Technical Guidance (EPA/100/R-12/001). /documents/ benchmark_dose_guidance.pdf (7 October 2015).

US EPA. 2015. Benchmark Dose Software (BMDS). Environmental Assessment. US EPA (16 April 2015).

Webster AF, Williams A, Recio L, Yauk CL. 2014. Gene expression analysis of livers from female B6C3F1 mice exposed to carcinogenic and non-carcinogenic doses of furan, with or without bromodeoxyuridine (BrdU) treatment. Genom. Data 2:117-122.

Webster AF, Chepelev N, Gagne R, Kuo B, Recio L, Williams A, Yauk CL. 2015. Impact of Genomics Platform and Statistical Filtering on Transcriptional Benchmark Doses (BMD) and Multiple Approaches for Selection of Chemical Point of Departure (PoD). PLoS One 10: e0136764.

Yang L, Allen BC, Thomas RS. 2007. BMDExpress: a software tool for the benchmark dose analyses of genomic data. BMC Genomics 8:387.

Zhao QJ, Haber L, Kohrman-Vincent M, Nance P, Dourson M. 2010. Quantitative modeling in noncancer risk assessment. In Quantitative Modeling in Toxicology, Krishnan K, Andersen ME (eds). John Wiley & Sons Ltd: West Sussex; 371-398.

Supporting information

Additional supporting informationmay be found in the online version of this article at the publisher's web-site