Shiny GATOM: omics-based identification of regulated metabolic modules in atom transition networks.

Multiple high-throughput omics techniques provide different angles on systematically quantifying and studying metabolic regulation of cellular processes. However, an unbiased analysis of such data and, in particular, integration of multiple types of data remains a challenge. Previously, for this purpose we developed GAM web-service for integrative metabolic network analysis. Here we describe an updated pipeline GATOM and the corresponding web-service Shiny GATOM, which takes as input transcriptional and/or metabolomic data and finds a metabolic subnetwork most regulated between the two conditions of interest. GATOM features a new metabolic network topology based on atom transition, which significantly improves interpretability of the analysis results. To address computational challenges arising with the new network topology, we introduce a new variant of the maximum weight connected subgraph problem and provide a corresponding exact solver. To make the used networks up-to-date we upgraded the KEGG-based network construction pipeline and developed one based on the Rhea database, which allows analysis of lipidomics data. Finally, we simplified local installation, providing R package mwcsr for solving relevant graph optimization problems and R package gatom, which implements the GATOM pipeline. The web-service is available at https://ctlab.itmo.ru/shiny/gatom and https://artyomovlab.wustl.edu/shiny/gatom.

INTRODUCTION

The study of metabolism has emerged as one of the important and promising areas of research in biology, with applications in the fields like immunology, cancer and ageing (1–5). As many of the cell functions are tied to metabolic processes (6–8), understanding metabolic regulation is indispensable for understanding cell biology.

Systematically, metabolism can be studied with a number of experimental omics: transcriptomics, proteomics, metabolomics and lipidomics. These techniques provide views on different aspects of metabolic regulation and can be useful in separate, as well as in integration. For the analysis of these data in the context of metabolism various tools were developed, such as IMPaLa (9), 3Omics (10) and MetaboAnalyst (11), or more focused for the analysis of lipidomics data like LipidSuite (12) or LipidSig (13). However, these tools heavily rely on the presence of metabolic pathway annotations which limits the analysis, in particular when only parts of the pathways are regulated or when complex pathway interactions are involved.

Network analysis methods, on the other hand, can be utilized for an unbiased analysis of metabolic interaction, without requirement for predefined pathways. In particular, so called active module approaches allow to find the most regulated subnetwork based on the provided omics data. These methods were successfully used to interpret transcriptomics data (jActivemodules (14), BioNet (15), KeyPathwayMineR (16)) and genomics data (HotNet2 (17), NetSig (18)). However, straightforward application of these approaches for studying metabolism is difficult due to unique features of metabolic networks.

Here, we present Shiny GATOM: an update of a previously published web-service GAM (https://academic.oup.com/nar/article/44/W1/W194/2499310) for integrated network analysis of metabolic and transcriptional data. While the overall approach remains the same: the method takes as input transcriptional and/or metabolic data and finds a metabolic subnetwork most regulated between the two conditions of interest, we introduced several major improvements both conceptual and technical. First, in GATOM we introduce the atom transition network structure which better captures pathway-like structure of the metabolic interactions. Second, to address the changes in the network structure we introduced a new computational approach for finding an active metabolic module. Finally, GATOM covers a wider range of reference networks, thus expanding the scope of its applicability, in particular for analysis of lipidomics data.

MATERIALS AND METHODS

Analysis pipeline overview

The GATOM network analysis pipeline is similar to the previous version implemented in GAM web-service (19). The key steps are:

The input table with differential expression results for genes and/or metabolites are overlaid on a global metabolic network to construct an experiment-specific network.

The obtained network is scored, so that nodes (metabolites) and edges (genes) with highly significant differential expression have positive scores and ones with not significant results have negative scores.

A connected subnetwork with the best total score is found by solving an optimization problem.

The found subnetwork represents an active metabolic module, capturing the most regulated reactions. It can be explored, annotated and downloaded.

Global network construction

The GATOM pipeline relies on three types of database-like objects:

A reaction network – contains a list of reactions, their reactants and mappings to enzyme classes.

A metabolite annotation – contains mappings between different metabolite identifiers.

An organism annotation – contains an organism-specific mapping between enzyme classes and genes, as well as mapping between different types of gene identifiers.

There are two major network sources used for global network construction: KEGG (20) and Rhea (21). The KEGG database is more mature, but the Rhea database provides significantly better coverage of the lipid species. For KEGG database, KEGG API is used for downloading the necessary information. For Rhea, the information is extracted from the files available for download.

An important step for both network databases is the atom mapping in reactions. For this purpose we use the Reaction Decoder Tool (RDT) (22). RDT takes as an input a reaction in RXN format and returns one-to-one correspondence between atoms in a substrate and atoms in a product. For Rhea database RXN files are immediately available for download, for KEGG database RXN files are generated by our pipeline from individual compound MOL files downloaded from KEGG database. The output of RDT is processed using the ChemmineR package (23).

Metabolite annotations are constructed from KEGG and ChEBI databases. Mappings between ChEBI and HMDB were obtained with metaboliteIDmapping package. Mappings between ChEBI and KEGG databases were extracted via KEGG API.

GATOM includes four organism annotations built-in: for human, mouse, arabidopsis and yeast. The mappings are generated based on the corresponding organism annotations packages from Bioconductor (org.Hs.eg.db, org.Mm.eg.db, org.Sc.sgd.db, org.At.tair.db) with additional enzyme mapping obtained from KEGG via its REST API. Additionally, selected metabolic pathways from KEGG and Reactome (24) databases are stored for module annotation purposes.

Finally, a lipid specific network is generated based on the Rhea network. SwissLipids database (25) is used to identify lipids among the reactants, and only the reactions with at least one lipid species are kept. SwissLipids and ChEBI databases are used to generate identifier mapping for lipid species. Mappings between lipid reactions and genes are enriched using direct UniProt references from the Rhea database.

Input data

As the input GATOM takes tables with differential expression results for genes and/or metabolites in text or XLSX formats. For genes, RNA-seq or microarray-based data tables are supported, in particular coming from tools like DESeq2 (26) or limma (27). For metabolites, metabolomics and lipidomics data tables are supported.

The following columns are expected:

‘ID’: RefSeq mRNA transcript ID, Entrez ID or symbols for genes; HMDB ID, KEGG ID or ChEBI ID for metabolites; SwissLipids ID, LipidMaps ID, ChEBI ID or Species name for lipids (supports LipidMaps and SwissLipids nomenclature). Multiple annotations for features are supported (e.g. multiple genes for a microarray probe or multiple metabolites for a mass-spectrometry peak), in that case IDs should be separated with ‘///’.

‘pval’: Differential expression P-value (non-adjusted).

‘log2FC’: Base 2 logarithm of the fold-change.

‘baseMean’ (for genes): average expression level.

‘signal’ (optional): ID of the measured entity such as probe ID for gene expression microarrays and peak ID for mass-spectrometry data. If absent, the ‘signal’ column will be generated based on unique ‘pval’ and ‘log2FC’ combinations.

Three example datasets are provided. Two of them with gene and metabolite differential expression between unstimulated and LPS + IFNg-stimulated macrophages are the same as described previously in GAM. Additionally we provide a lipidomics dataset with comparison of high fat and normal fat diet in peritoneal macrophages. For this we obtained the dataset ST001289 (28) from Metabolomic Workbench. Only wildtype samples were selected (without LDLR knock-out), and the differential expression analysis was performed between normal cholesterol, normal fat (NCNF) diet and high cholesterol, high fat (HCHF) diet samples using LipidSig web-service (13). The corresponding file is available as Supplementary Table S1.

Finding active metabolic module

Following the GAM pipeline, GATOM is built on the maximum-likelihood formulation of the active module problem introduced by Dittrich et al. (29). However, instead of directly scoring genes and metabolites, as was done in GAM, in GATOM we introduce a notion of signals: the entities that were measured and analysed (for example, microarray probes or mass-spectrometry ions). Consequently, an individual signal can be associated with multiple genes or metabolites and assigned to multiple reaction network elements. Signals for each data type are scored independently.

Accordingly, we define a new variant of the maximum weight connected subgraph (MWCS) problem - the signal generalized maximum weight connected subgraph (SGMWCS) problem. Let consider a graph G with vertices V, and edges E, and a set of signals S. Each signal has a score w(S) and can be assigned to multiple graph elements (both vertices and edges). The score of a subgraph is defined as a sum of scores of unique signals assigned to its vertices and edges. The SGMWCS problem consists in finding a subgraph with a maximal score.

The SGMWCS problem is NP-hard as well as MWCS, as the former is a more general case. To be able to solve instances of the problem, we developed an exact solver Virgo-solver (https://github.com/ctlab/virgo-solver) based on reduction to mixed integer linear programming (MILP) problem and solving it using IBM CPLEX library, extending the solver previously developed by our group (30).

Postprocessing

After the construction of the module post-processing steps can be done within the web-service.

First, reactions not controlled by products of highly changing genes but controlled by products of genes with high average expression can be added, using the baseMean input column.

Second, as GATOM works with atom-transition graphs, when a module is found the same metabolite can appear multiple times in the module via different atoms. In such cases it may be useful to connect atoms belonging to the same metabolite with edges. The other way is to collapse them into one vertex.

Finally, canonical metabolic pathways can be highlighted in the module. Hypergeometric test for enrichment is applied with genes from the scored network used as a universe and genes from the module as the query gene set. Pathways reaching an adjusted P-value of 0.05 are shown, and reactions with the corresponding enzymes can be selected.

The resulting module can be downloaded in SVG format and can easily be later edited in any vector image editor. Alternatively, the module can be downloaded in XGMML format and imported into Cytoscape desktop application.

Implementation details

Shiny GATOM is implemented in R and uses shiny library for interaction with the user. The source code is available at https://github.com/ctlab/shinyGatom. shinyCyJS and cytoscape.js (31) libraries are used for interactive module display. The analysis pipeline itself is implemented as an R package gatom (https://github.com/ctlab/gatom). Gatom package depends on the ability to solve the SGMWCS problem. For this we developed an R package mwcsr (https://cran.r-project.org/package=mwcsr), which implements an interface to exact SGMWCS solver Virgo-solver, along with two heuristic solvers: one based on simulated annealing and one based on relax-and-cut procedure.

RESULTS

Atom-based network topology

In GATOM we introduce an atom-based metabolic network topology. There, each vertex corresponds to an individual carbon atom of a metabolite. Two atoms A1 and A2 of metabolites M1 and M2 are connected if there is a reaction R with M1 and M2 on different sides of the equation and atom A1 of M1 transforms into atom A2 of M2 in that reaction. Useful notion is that if atom A1 of metabolite M1 is heavy-labeled (for example replaced with C13) then atom A2 will be also labeled.

Atom-based network topology has a useful property: any path from one vertex to another corresponds to a valid sequence of reactions. This property has been previously used for identification of pathways (32,33) as it forces a pathway-like structure of the subgraph. In the context of identifying active metabolic modules it simplifies the interpretation. In particular, linear pathways such as glycolysis are naturally recovered.

Figure 1 shows the active module found for macrophages activated with IFNg + LPS based on the atom transition network. The module displays major metabolic regulation of glycolysis, TCA cycle, pentose phosphate pathway and fatty acid metabolism. For comparison, the module obtained for the metabolite-level network is shown on Supplementary Figure S1. While the overall present pathways are similar, the connections between pathways in the atom based module are more canonical and consistent with the findings from the original publication (34): glycolysis pathway directly feeds into TCA cycle, and arginine metabolism fragment is connected to TCA cycle via glutamate.

Figure 1.

Metabolic module obtained for the example data of macrophages activated by IFNg + LPS. The KEGG network with atom topology was used. The module layout was done in Shiny GATOM. The shown pathway annotations were added separately.

Signal maximum weight connected subgraph problem

As in the previous version, the key step of the GATOM analysis is finding a connected subnetwork of a metabolic network that captures the most changing genes and metabolites. Mathematically, this corresponds to solving a Maximum Weight Connected Subgraph (MWCS) problem. However, modifications from a classic problem where only nodes are weighted are required to better capture features of metabolic networks.

Here, we introduce a new variant of the MWCS problem, the signal generalized maximum weight connected subgraph (SGMWCS) problem. There, as in the previously defined generalized MWCS (GMWCS) problem, both nodes and edges can be scored (e.g. from metabolomic and transcriptomic data respectively). However, the scores are assigned indirectly via a notion of signals. Intuitively, a signal in SGMWCS corresponds to a measured entity (e.g. mass-spectrometry ion or a gene) with its score representing a significance of the change between two conditions (as in GAM). These signals can be repetitive in the metabolic network, such an ion corresponding to multiple metabolites or a gene participating in multiple reactions. Consequently, the weight of the subgraph in SGMWCS is defined as a sum of scores of unique signals.

To solve SGMWCS instances we extended the previously developed GMWCS solver (30). Briefly, the solver uses a reduction to the mixed integer linear programming (MILP), which then can be solvedto a provable optimality with IBM ILOG CPLEX library. The core of the MILP formulation and used heuristics remain the same as in GMWCS, with the necessary adaptations to account for the signal-based scoring. The solver supports setting a time limit, after which the current best solution is returned. Similar to the GMWCS solver performance benchmark (30) we have run SGMWCS solver on instances generated from the datasets submitted to Shiny GAM during the testing phase. For gene only inputs (66 instances) median time to solve the SGMWCS problem to optimality was 6 s for metabolite topology and 25 s for atom topology. For simultaneous gene and metabolite inputs (27 instances) median time to solve the SGMWCS problem to optimality was 18 s for metabolite topology and 48 seconds for atom topology. The benchmarks were run on a laptop with Intel Core i5-7200U CPU @ 2.50GHz.

Formulating the active module problem as the SGMWCS problem prevents artifacts, which happens when multiple instances of the same significantly changed signal is repeated multiple times within a small region of the network, effectively multiplying its score. To illustrate the problem we have applied the GMWCS solver for finding an active module for the example macrophage activation dataset, with the only difference from the GATOM pipeline in counting each gene and metabolite as many times as it appears in the module. The resulting module for metabolite topology (Supplementary Figure S2) contains two artifacts that add major bias for the module score calculation: repeated appearance of Mgst3 gene and repeated appearances of phospholipid metabolism genes (Pla2g4a, Pla2g7, Plaat3, Lipg, Cept1). In the module for the atom topology (Supplementary Figure S3) these artifacts are multiplied further, as each metabolite appears in the atom transition network several times, creating even further bias. Previously, in GAM, we addressed this problem by using different network structures (interpreting reactions as nodes or as edges) depending on the input data types. Introduction of SGMWCS allows the user to use the desired network structure as it is, independent of the provided types of data.

Updated networks

Compared to the previous release, we have updated the network construction pipeline so that up-to-date reaction databases are used. Importantly, on top of the upgrade to the recent versions of KEGG database, we have added support for Rhea database. Overall the KEGG-based networks and Rhea-based networks are similar, producing the similar analysis results. However, the Rhea-based networks provide much better coverage of lipid species and thus allows focused analysis of lipidomics data. For that purpose we created a lipid-specific Rhea subnetwork. The sizes of the networks (in terms of associated genes) are shown in Table 1.

Table 1.

Sizes of Shiny GATOM built-in metabolic networks. Number of genes corresponding to reactions in the networks are shown

	Organism
Network type	Human	Mouse	Arabidopsis	Yeast
KEGG	2641	2713	3212	1072
Rhea	2719	2799	3270	1153
Rhea-lipidomic	626	624	582	260

An important aspect of working with lipids is the existence of complex lipid hierarchy spanning from high-level lipid classes to isomeric subspecies. For conciseness, reactions in Rhea databases use lipid identifiers from all of these levels. To account for that we use the following scheme of conversion between input lipids and ChEBI IDs used in the Rhea network. First, using SwissLipids hierarchy we group the lipids subspecies based on their species. We denote the corresponding species as a representative of the group. Next, data entries for any of the group elements are assigned to all of the group members. After that the data entries are propagated up by the hierarchy to lipid classes. Finally, thus expanded entries are mapped from SwissLipids to ChEBI using the direct references available at these databases.

An R package

The core functions of the web-service were implemented as an R package gatom available at GitHub repository https://github.com/ctlab/gatom. The package can be easily installed locally if an R environment is available. Using the package one can carry the analysis including input-data handling, finding the active metabolic modules and visualizing the result. Additionally, the package can be used to create new organism annotations based on corresponding Bioconductor packages.

GATOM analysis depends on the ability to solve the maximum weight connected subgraph problem. To simplify the set up we additionally developed an R package mwcsr (https://cran.r-project.org/package=mwcsr) containing an R interface for the exact Virgo solver (https://github.com/ctlab/virgo-solver) which requires IBM ILOG CPLEX libraries, as well as a number of heuristic solvers which do not require any dependencies.

Case study: Trem2 deficiency in bone-marrow derived macrophages

As a case study we applied GATOM to transcriptomics and metabolomics profiling data of bone-marrow derived macrophages from WT and Trem2-deficient mice (35). The input data is provided as Supplementary Tables S2 and S3.

The differential expression tables were loaded into Shiny GATOM web-service and analysed using atom and metabolite topologies with the default parameters, as well as loaded into previously developed Shiny GAM web-service. The module obtained for the atom transition network is shown on Supplementary Figure S4, for the metabolite topology the module is shown on Supplementary Figure S5, and the module obtained with Shiny GAM is shown on Supplementary Figure S6. Importantly, only the atom based module shows deregulation of the glycolysis pathway, with the other modules containing only individual glycolysis reactions. In the study (35) deficiency of the glycolysis pathway in Trem2 KO macrophages was validated by measurement of extracellular acidification rate (ECAR).

Case study: lipidomics

As a case study of working with lipidomics data we compared peritoneal macrophages from mice fed with normal fat and high fat diet (28). The dataset ST001289 was obtained from Metabolomic Workbench. The differential abundance analysis was performed In LipidSig and the results were uploaded into Shiny GATOM. The input data are provided as Supplementary Table S1 and as the example dataset in the web-service.

The resulting module is shown on Figure 2. The module shows up-regulation of cholesterol derivatives in macrophages from high fat diet samples, as well as major deregulation of glycerphospholipid metabolism. Upregulated cholesterol is naturally expected under the high-cholesterol diet and might be considered as a positive control. As many other macrophages, peritoneal macrophages have high levels of base metabolism of phospholipids, comparing to other myeloid cells (monocytes and dendritic cells) (36), thus glycerophospholipid deregulation on high-fat diet might reflect the response to maintain phospholipid homeostasis. Interestingly, arachidonoyl-CoA (labeled as ‘(5Z,8Z,11Z,14Z)-eicosatetraenoyl-CoA’ in the module) appears in the module surrounded by changing phospholipids, without being measured. This is consistent with arachidonic metabolism being important in atherosclerosis and being a viable therapeutic target (37).

Figure 2.

Metabolic module obtained for the example lipidomics data, comparing peritoneal macrophages from mice fed with normal fat and high fat diet.

DISCUSSION

In this study, we introduce a novel GATOM approach representing an updated version of previously published GAM method for analysis of ‘omics datasets (19). First, being free from predefined annotation of reference metabolic network, GATOM due to its atom based approach still produced final modules having pathway-like structure. Second, a signal variant of maximum weight connected subgraph problem (SGMWCS) was formulated to account for the updated network structure, and a corresponding exact solver was developed. Taken together, the developed approach simplifies the interpretation of the obtained results as any path from one vertex to another corresponds to a valid sequence of reactions.

Additionally, the Rhea network was specifically preprocessed and integrated into the GATOM analysis pipeline. Rhea network has much better coverage of lipid species, compared to KEGG, and allows analysis of lipid-rich samples like brain (38) or adipose tissues (39). Lipid metabolism also is recognized as an important regulatory circuit in immunity (40).

The GATOM method has been tested in multiple studies. For example, in (35) it was used to show deregulation of energy metabolism in Trem2-deficient macrophages. In (41) the GATOM analysis highlighted activation of inositol-triphosphate metabolism in fasting mice, indicating a potential signalling mechanism.

Finally, the Shiny GATOM web-service provides easy access to the GATOM method. Still, GATOM analysis can be carried locally in an R environment. Since the previous version, the underlying implementation was greatly reworked, simplifying the installation. GATOM is available as an R package on GitHub and depends on the mwcsr R package developed by our group with interfaces to a number of solvers, both exact and heuristic.

DATA AVAILABILITY

Source code of the Shiny GATOM is available in several GitHub repositories:

Web-service implementation: https://github.com/ctlab/shinyGatom

R package with a GATOM pipeline: https://github.com/ctlab/gatom

Global network construction pipelines for Rhea and KEGG: https://github.com/ctlab/Rhea-network-pipeline, https://github.com/ctlab/KEGG-network-pipeline

Click here for additional data file.