Application of a metabolic network-based graph neural network for the identification of toxicant-induced perturbations

Abstract

Transcriptomic analyses have been an effective approach to investigate the biological responses and metabolic perturbations by environmental contaminants in rodent models. However, it is well recognized that metabolic networks are highly connected and complex, and that traditional gene expression analysis methods, including pathway analyses, have a limited ability to capture these complexities. Given that metabolism can be effectively represented as a graph, this study aims to apply a network-based graph neural network (GNN) to uncover novel or hidden metabolic perturbations in response to a toxicant. A GNN model based on the mouse Reactome pathways was trained and validated on 7,689 transcriptomic samples from 26 mouse tissues curated from Recount3. This model was then used to identify important reactions in publicly available data from livers of mice treated with the environmental contaminant 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) achieving a performance of 100% when comparing a single dose to a control group. Integrated gradients and centrality analyses identified perturbation of the SUMOylation, cell cycle, P53 signaling, and collagen biosynthesis pathways by TCDD which were not identified using a pathway analysis approach. Collectively, our results demonstrate that GNNs can reveal novel mechanistic insights into toxicant-mediated metabolic disruption, presenting a putative strategy to characterize biological responses to toxicant exposures. Our studies illustrate how the use of a reaction-based graph neural network can support the discovery of toxicant-induced metabolic perturbations, and highlight strengths and challenges in the application of artificial intelligence methods for environmental health research.

Keywords: artificial intelligence; graph neural networks; metabolism; mice; toxicogenetics.

Fig. 1.

Training and validation of a reaction-based GNN on bulk RNA sequencing datasets from disparate mouse tissues. (A) Distribution of reaction-wise principal component score distribution for the first 10 components. A total of 9,844 reactions are shown and the box represents the median, 25th, and 75th percentiles. (B) Reaction-wise ARIs for each tissue were calculated following clustering using a KNN algorithm. Bars represent the average of the top 10 ARIs for each tissue. (C) Training performance at each epoch for GNN and ResNet18 models, as well as control datasets, using reaction-wise PC1 values. The final accuracy is shown in the figure legend. ROC curves are shown for validation datasets using the (D) GNN and (E) ResNet18 models. AUROCs are indicated next to sample labels for each tissue type. Additional performance metrics are in Table S4.

Fig. 2.

Testing performance for classification of TCDD test datasets by GNN model. Sample dose classification compared with ground truth from the GNN model (A) including or (B) excluding the time-course data (SRP131784). The point size represents the total number (log₁₀ scale) of samples in each classification and true dose category. (C) Final accuracy at 500 epochs for the classification of the dose–response TCDD test dataset.

Fig. 3.

Visualization of key reaction nodes identified using IG values and centrality measurements. (A) Reactions identified as the highest IG values without centrality measurements found R-MMU-8956040 (COP9 signalosome deneddylates cytosolic CRL E3 ubiquitin ligase complexes”) to R-MMU-8956140 I (NEDD8 and UBD bind NUB1 and the 26S proteasome) and R-MMU-2002460 (P4HB binds Collagen chains) to R-MMU-1650808 (Prolyl 4-hydroxylase converts collagen prolines to 4-hydroxyprolines) with IG values larger than 0.01. (B) Using closeness centrality analysis, key reactions found were R-MMU-6799246 (CHEK1 phosphorylates TP53) and R-MMU-5693609 (ATM phosphorylates TP53 at S15). (C) Eigenvector centrality found R-MMU-2993780 (Transfer of SUMO1 from E1 to UBE2I (UBC9)).

Fig. 4.

The number of genes identified uniquely or commonly by Differential Expression Analysis and the GNN. (A) Number of different and common genes between Differential Expression Analysis and GNN. (B) Number of different and common genes between Enrichment Analysis and GNN.