BMDExpress Data Viewer - a visualization tool to analyze BMDExpress datasets.

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 (http://apps.sciome.com:8082/BMDX_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.

Introduction

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‐adverse‐effects‐level approach (Crump, 1984). Regulatory agencies have increasingly adopted BMD modeling for human health risk assessment (Budtz‐Jorgensen et al., 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 al. 2007). This approach has already been applied to a variety of chemicals (Bhat et al., 2013; Moffat et al., 2015; Thomas et al., 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 al. (2007) developed a tool called “BMDExpress” that automates the statistical analysis and BMD modeling computational steps to facilitate timely analysis of very 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 http://apps.sciome.com:8082/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 × 44 K and 8 × 60 K), Agilent mouse whole genome (4 × 44K and 8 × 60 K), Agilent rat whole genome (4 × 44 K and 8 × 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 × 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 (http://www.ncbi.nlm.nih.gov/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.

Results

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).

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.

Example analysis with furan‐treated mouse liver and 4,4′‐methylenebis(N,N‐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 filtered datasets contained 29847, 2597 and 362 probes, for the unfiltered, 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 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.

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.

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.

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 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.

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.

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.

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.

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 al., 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 et al. [2014]), we obtained 24 pathways (Table 1). Fourteen (58%) of the identified pathways were also found in our previous IPA analysis (Jackson et al., 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.

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 (mg kg–1 day–1)	Fisher's exact test P‐value (%)	Number of 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).

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

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 al., 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 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).

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.

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 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 glioblastoma multiforme signaling, ILK signaling and nitric oxide signaling in the cardiovascular system pathways showed increasing BMD mean values from 28 to 90 days. For the 14‐day timepoint, the hepatic fibrosis/hepatic stellate cell 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.

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.

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 stellate cell activation pathway had a higher BMDL and BMD 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.

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.

Discussion

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 al., 2015; Thomas et al., 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 al., 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 roles of polo‐like 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 et al., 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 al., 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 endothelin‐1 signaling, HMGB1 signaling, IL‐8 signaling 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.

Conflict of interest

The authors did not report any conflict of interest.

Supporting info item

Click here for additional data file.

Supporting info item

Click here for additional data file.

Supporting Info item

Click here for additional data file.