Latent Variables Capture Pathway-Level Points of Departure in High-Throughput Toxicogenomic Data

Abstract

Estimation of points of departure (PoDs) from high-throughput transcriptomic data (HTTr) represents a key step in the development of next-generation risk assessment (NGRA). Current approaches mainly rely on single key gene targets, which are constrained by the information currently available in the knowledge base and make interpretation challenging as scientists need to interpret PoDs for thousands of genes or hundreds of pathways. In this work, we aimed to address these issues by developing a computational workflow to investigate the pathway concentration-response relationships in a way that is not fully constrained by known biology and also facilitates interpretation. We employed the Pathway-Level Information ExtractoR (PLIER) to identify latent variables (LVs) describing biological activity and then investigated in vitro LVs' concentration-response relationships using the ToxCast pipeline. We applied this methodology to a published transcriptomic concentration-response data set for 44 chemicals in MCF-7 cells and showed that our workflow can capture known biological activity and discriminate between estrogenic and antiestrogenic compounds as well as activity not aligning with the existing knowledge base, which may be relevant in a risk assessment scenario. Moreover, we were able to identify the known estrogen activity in compounds that are not well-established ER agonists/antagonists supporting the use of the workflow in read-across. Next, we transferred its application to chemical compounds tested in HepG2, HepaRG, and MCF-7 cells and showed that PoD estimates are in strong agreement with those estimated using a recently developed Bayesian approach (cor = 0.89) and in weak agreement with those estimated using a well-established approach such as BMDExpress2 (cor = 0.57). These results demonstrate the effectiveness of using PLIER in a concentration-response scenario to investigate pathway activity in a way that is not fully constrained by the knowledge base and to ease the biological interpretation and support the development of an NGRA framework with the ability to improve current risk assessment strategies for chemicals using new approach methodologies.

Figure 1

Computational workflow overview. First, a concentration–response HTTr data set is used as the input into PLIER along with a collection of predetermined gene sets (step 1). PLIER returns a collection of LVs, which represent patterns of co-regulated genes that may or may not align with any of the gene sets supplied. Next, the concentration–response activity of each LV in each compound is investigated using a new updated version of the ToxCast pipeline (tcplfit2).

Figure 2

Transcriptional activity after compound application. (A) Magnitude of activity in terms of active LVs found by the workflow developed in this work. In cyan, we report those LVs that are found to align with prior knowledge and in red the ones that are found to not align with prior knowledge. (B) Magnitude of perturbed genes across concentrations as identified by Harrill et al.

Figure 3

Concentration-dependent activity of LVs for estrogenic (bisphenol A, bisphenol B, 4-nonylphenol, and 4-cumylphenol) and antiestrogenic (4-hydroxytamoxifen, clomiphene citrate, and fulvestrant) compounds. (A) The plots show LVs found to be active in each compound along with their PoD across the dose-range tested. The presence or absence of association with biology is reported in cyan and red, respectively. (B) Gene sets aligning with LV 30, which was found to be active in all the estrogenic and antiestrogenic compounds along with the AUC and the FDR set at 0.7 and 0.05, respectively. (C) Concentration–response models for LV30 in response to ER agonists and ER antagonists. For each panel, the noise band is represented by the gray band spanning zero, while the estimated PoD (green line) is reported along with its confidence intervals (green box). mthd = best-fit concentration–response model; Hitcall = confidence for the fitted model; top = top parameter for the fitted model.

Figure 4

PoD estimates of LVs not aligning with prior knowledge. (A) Distribution of PoDs for compounds having the lowest LV not associated with prior knowledge. The plots display the distribution of the different PoDs (each dot, ordered from lowest to highest) along with their confidence intervals, estimated using the approach developed in this study. PoDs are colored in cyan or red depending on whether they have an association with existing gene sets (aligning) or not (not aligning), respectively. (B) The scatterplot shows the agreement between PoD estimates across the two methods: the gene signature (GS) and PLIER. Chemicals are colored depending on whether the LV with the lowest PoD aligns (cyan) or not (red) with prior knowledge.

Figure 5

LV activity. The plot shows LVs found to be active in each compound across the different cell lines. The y-axis reports the different LVs found to be active, while the x-axis shows the concentration range tested in the log scale. The presence or absence of association with biology is reported in cyan and red, respectively.

Figure 6

Summary of the biological activity for andrographolide, flutamide, and triclosan. (A) The table displays LVs found to be active in each treatment/cell line along with the AUC and the FDR. (B) The concentration–response plots show example activities for some of the LVs listed in (A), modeled using tcplfit2. For each plot, the noise band is represented by the gray band spanning zero, while the estimated PoD (green line) is reported along with its confidence intervals (green box). mthd = best-fit concentration–response model; Hitcall = confidence for the fitted model; top = top parameter for the fitted model.

Figure 7

Agreement between PoD estimates across methods. The CCC is reported along with the CI.