Development of a prioritization method for chemical-mediated effects on steroidogenesis using an integrated statistical analysis of high-throughput H295R data

Abstract

Synthesis of 11 steroid hormones in human adrenocortical carcinoma cells (H295R) was measured in a high-throughput steroidogenesis assay (HT-H295R) for 656 chemicals in concentration-response as part of the US Environmental Protection Agency's ToxCast program. This work extends previous analysis of the HT-H295R dataset and model by examining the utility of a novel prioritization metric based on the Mahalanobis distance that reduced these 11-dimensional data to 1-dimension via calculation of a mean Mahalanobis distance (mMd) at each chemical concentration screened for all hormone measures available. Herein, we evaluated the robustness of mMd values, and demonstrate that covariance and variance of the hormones measured appear independent of the chemicals screened and are inherent to the assay; the Type I error rate of the mMd method is less than 1%; and, absolute fold changes (up or down) of 1.5 to 2-fold have sufficient power for statistical significance. As a case study, we examined hormone responses for aromatase inhibitors in the HT-H295R assay and found high concordance with other ToxCast assays for known aromatase inhibitors. Finally, we used mMd and other ToxCast cytotoxicity data to demonstrate prioritization of the most selective and active chemicals as candidates for further in vitro or in silico screening.

Keywords: Endocrine disruption; Prioritization; Steroidogenesis; ToxCast.

Figure 1.. Workflow of data simulation study.

(A) Original level 0 data from the ToxCast database, invitrodb, were downloaded and converted into micromolar concentrations and were subset into four datasets to represent vehicle controls, low, medium, and high responses in the HT-H295R assay based on the maxmMd values from Haggard et al. (2018). These subsets were used to calculate 4 covariance matrices for each block of assay data (for a total of 32 covariance matrices) for MVN sampling. Using mean hormone response values by chemical sample and the appropriate covariance matrices, 2000 unique simulated datasets were generated. (B) For each simulated dataset, three sets of mMds were calculated for all simulated chemical samples using three covariance matrices estimated from either low responding chemicals, high responding chemicals, or the global dataset, for a total of 6000 sets of mMd values.

Figure 2.. Example of Mahalanobis distance.

The difference between Euclidean distance and Mahalanobis distance using an example to two theoretical hormone measures with positive covariance. (A) Scatterplot for two theoretical hormone measures X and Y in Euclidean space. Grey points denote mean hormone measures for theoretical chemical samples. The error distribution, denoted by the dashed ellipse, shows the variance among measurements of hormone Y (standard deviation = 0.16) is greater than measurements of hormone X (standard deviation = 0.08), and these measures show positive covariance (correlation = 0.8). The points labeled conc 1, conc 2, and conc 3 represent the mean (natural logarithm) concentrations of three concentrations of an example test chemical. In terms of hormone concentrations, the response at conc 3 for hormone X is roughly twice as far from conc 1 as is the response at conc 2 (and is a larger number of standard deviations away from the mean of the distribution for hormone X); however, the Euclidean distances of conc 2 and conc 3 to conc 1 are the same. (B) Scatterplot of the same hormone measures as in A after Mahalanobis distance transformation. Mahalanobis distance adjusts for the covariance between measurements, which rotates and scales the data so that the data are uncorrelated and have the same variance (as denoted by the error ellipse now being a circle). After transformation, we see that conc 3 is now around four times as distant from conc 1 as is conc 2, better reflecting the difference of conc 3 from the overall distribution of measures in hormones X and Y compared to conc 2.

Figure 3.. Hormone variance and pairwise correlation of covariances.

The pairwise correlation of the mean pooled covariance matrix of the hormones for the low (italics), global (boldface), and high (normal) covariance types are shown. The diagonal of the matrix indicates the mean variance for each measured hormone. The correlation matrix is colored (red to blue; −1.0 to 1.0) based on the correlation coefficients of the global covariance matrices (boldface) and are clustered using Ward’s method.

Figure 4.. Stability of mMd values calculated with different covariance matrices.

For each data simulation, three sets of mMds were calculated, one set for each of the three covariance types. (A) Scatterplot to compare the original mMd values calculated in Haggard et al. (2018) to the simulated mMd values calculated for each covariance type. Closed circles represent the median and the bars represent the 95% confidence interval of the mMd values for each covariance type. Dashed line indicates the identity line. (B) Boxplots of mMd concentration response values after prochloraz exposure in HT-H295R across the 2000 simulation and each covariance type. (C) Boxplots of mMd responses after butylparaben exposure in HT-H295R across the 2000 simulations and each covariance type. (D) Boxplots of mMd concentration response values after ethylene dimethanesulfonate exposure in HT-H295R across the 2000 simulation and each covariance type. For all plots, red, green, and blue correspond to the low, global, and high covariance type, respectively. Filled triangles represent mMd values estimated using the original dataset using the three covariance types.

Figure 5.. Power analysis.

As part of the MVN simulation, we generated ten different theoretical hormone response profiles at varying effects sizes (1.1, 1.5, 2, and 2.5-fold increases compared to a DMSO control sample), and calculated mMds for these profiles using the three covariance types. The percentage of maxmMd values above the pre-defined critical value at α = 0.01 (i.e. maxmMd ≥ 1.64) was considered the power to detect a maxmMd with experimental significance. Red, green, and blue correspond to the low, global, and high covariance type, respectively. Circles and triangles denote effect sizes of 1.1- and 1.5-fold, respectively. Note that effect sizes of 2.0- and 2.5-fold are not included in the figure due to these effect sizes having a power of ≥0.999 for all response profiles. Dashed line indicates a power of 0.8. Definitions for the response profile types are as follows: All: all measured hormones increased; Some: CORT, 11-DCORT, ANDR, TESTO, E1, and E2 increased; Glucocorticoids: 11-DCORT and CORT increased; Mineralocorticoids: DOC and CORTIC increased; Androstenedione: ANDR increased; Testosterone: TESTO increased; Androgens: ANDR and TESTO increased; Estrone: E1 increased; Estradiol: E2 increased; Estrogens: E1 and E2 increased.

Figure 6.. Heatmap of hormone response profiles.

Putative aromatase inhibitors were identified after filtering HT-H295R response profiles for chemical samples with statistically significant decreases in E1 and E2 concentrations compared to plate-level DMSO controls (as measured by one-way ANOVA) with a fold decrease of ≥1.5. Rows denote chemical samples and columns are hormones. Rows were clustered using Ward’s method only considering ANDR, TESTO, E1, and E2, and four prominent clusters were identified and are labeled. The aromatase relevant hormones, specifically ANDR, TESTO, E1, and E2, are presented to the left whereas the other hormones are presented to the right. The log₁₀-maxmMd value is also denoted with an annotation bar to the right to indicate effect size for the 11 hormones. Reference chemicals (chrysin, clotrimazole, metyrapone, letrozole, and fadrozole hydrochloride) known to inhibit aromatase are specifically labeled to the right (Judson et al., 2018b). Grey denotes missing data due to sample loss during screening. Cluster numbers 1-4 referenced in the text are labeled.

Figure 7.. Selectivity, area under the curve, and maxmMd as context for prioritization.

The most (top 25) and least (bottom 25) efficacious selective chemicals are illustrated in Figure 7 (for all, see Supplemental Figure 5 and Supplemental Table 1). The BMD is the potency at the critical value for the mMd curve (red); the MTT_acc is the potency at the threshold for a positive response in that assay (green); and the cyto_burst is the potency value for the lower bound of a cytotoxicity distribution (blue). Selectivity (minimum of cyto_burst or MTT_acc minus the BMD), AUC (area under the curve for mMd versus concentration), and maxmMd are shown to the right.