Measurement(s)gene expression, metabolic relationshipsTechnology Type(s)Genome Scale Metabolic Models; Computational network biologySample Characteristic - OrganismHomo sapiens
Studies about the metabolic alterations during tumorigenesis have increased our knowledge of the underlying mechanisms and consequences, which are important for diagnostic and therapeutic investigations. In this scenario and in the era of systems biology, metabolic networks have become a powerful tool to unravel the complexity of the cancer metabolic machinery and the heterogeneity of this disease. Here, we present TumorMet, a repository of tumor metabolic networks extracted from context-specific Genome-Scale Metabolic Models, as a benchmark for graph machine learning algorithms and network analyses. This repository has an extended scope for use in graph classification, clustering, community detection, and graph embedding studies. Along with the data, we developed and provided Met2Graph, an R package for creating three different types of metabolic graphs, depending on the desired nodes and edges: Metabolites-, Enzymes-, and Reactions-based graphs. This package allows the easy generation of datasets for downstream analysis.
Background & SummaryCancer is a complex disease caused by a myriad of factors and characterized by an astonishing complexity of phenotypes and traits, which determine its wide heterogeneity, even among cells of a single tissue. Nonetheless, three key processes are shared by all cancer cells: proliferation, invasion, and metastasis. To fulfill these tasks, cancer cells need to reprogram their metabolic activities and cross-talk with their neighborhood1,2. This evidence gives the metabolism and its players a crucial role in cancer progression and, consequently, cancer research.Among all the biological networks, the metabolic ones are particularly complex and highly interconnected. Still, they probably are the best characterized in terms of connections and those that better represent the genotype-phenotype associations3. According to this, the reconstruction of comprehensive networks through the integration of omics data into metabolic scaffolds is one of the tools preferred by the systems biology approach for investigating biological phenomena from a holistic point of view. The metabolic scaffolds are given by the Genome-Scale Metabolic Models (GSMs), built from multi-omics data integration, and carrying information concerning the genes/proteins with enzymatic activity, how they interact with bioactive compounds in the context of biochemical reactions, and how the metabolic interconnections change in different cells, tissues or specific conditions4. There is a great interest in exploiting these models to generate condition-specific graphs at the service of machine learning approaches. In the era of precision medicine, the main goal is to develop approaches and tools to face the well-known heterogeneity of physiological and pathological manifestations and provide focused solutions for specific conditions. Considering the disease cohort as a single group, including all the diagnosed patients, is a simplistic approach that does not contemplate any inter-samples heterogeneity due to genetic and environmental factors. While modern biology has accepted the intra-sample heterogeneity of single cells, it seems anachronistic to talk about disease- instead of patient-specific conditions. There are several studies that address the problem of heterogeneity by exploiting network-structured approaches5–7.Metabolic networks are complex and can involve different metabolic players (i.e., metabolites, enzymes, reactions). Machine and deep learning frameworks allow extracting knowledge from the metabolic networks while dealing with their structural and relational complexity4. In the context of findability, accessibility, interoperability, and reusability (FAIR) principles8, providing benchmark datasets for comparing novel approaches and for the general advancement of a specific research domain is extremely important. Graph-structured data coupled with machine learning approaches are receiving growing interest9–13, and many benchmark datasets have been proposed in the context of biomedical graphs, especially derived from protein-protein interaction, chemical, imaging data14–18. To the best of our knowledge, metabolic networks based on context- and patient-specific metabolic models have not been provided so far. To fill this gap, here, we provide the TumorMet repository. TumorMet contains two main sets of networks depending on the models from which they derive: Tissue-derived networks generated starting from tissue-specific models and PDGSMMs-derived networks obtained using Patient-Derived Genome-Scale Metabolic Models (PDGSMMs). The interesting implications of using the metabolic networks are twofold, from both a computational and biological perspective. Their complexity in terms of nodes and connections, and the plasticity given by the multiple ways in which they can be generated, make them appealing for the proposal and validation of novel approaches in the context of computational graph-based research. In this work, we presented three alternatives, each focused on a specific set of metabolic players (i.e., metabolites, enzymes, and reactions). As demonstrated by19, reconstruction algorithms used to generate context-specific models present a bug which determines an underestimation of the molecular context. The model’s conversion into a network allows further contextualization by integrating context-specific data. Being aware that the networks we generated for TumorMet are just a portion of the possibilities, we provided the Met2Graph package to give the user the freedom to build the networks depending on specific needs. Met2Graph indeed implements a flexible process flow to build the metabolic graphs, can be easily integrated with user-customized functions, and provides several arguments to personalize the networks. Some of the networks in this dataset were used for assessing graphs classification, clustering, and embedding20–23, as well as for multimodal data analysis24,25, demonstrating their benefits. An exciting field of biological network usage is also represented by the application of node classification approaches aimed at predicting the essential genes, namely those genes crucial for an organism’s viability. Usually, the Protein-Protein Interaction (PPI) networks are exploited to this extent, based on the assumption that the topological centrality is correlated to a functional centrality. As hypothesized in26, one of the reasons why the PPI are the most used networks for this purpose could be their abundance compared to the other types, such as Metabolic networks, highlighting the importance of providing network datasets. Still, only physical interactions, additionally not contextualized, are insufficient to represent the genetic connections’ complexity27. Modern biology extensively uses networks to integrate and analyze data in a way in which organisms, tissues, or cells are considered systems. This perspective gives a crucial role to the connections among biological components, and the network-based analyses are exploited for making relevant biological inferences. The central role of metabolism in different aspects of pathophysiological mechanisms and their tune regulation make these networks particularly interesting for extracting knowledge and making predictions. For example, the analysis of hub nodes28 and the comparison of topological properties between different context-specific networks29 are valuable resources in diagnostic and prognostic markers investigation for precision medicine. Along with the data, we also provide an R package, Met2Graph, to create metabolic graphs starting from GSMs and gene expression data. The package can generate three types of graphs, depending on the desired nodes and edges: Metabolites-based graphs, where metabolites are nodes connected by reactant-product relationships and the edges can be weighted by expression values of the enzymes catalyzing the corresponding reactions; Enzymes-based graphs, where enzymes are nodes that are connected if they catalyze two reactions, each producing and consuming a specific metabolite; and Reactions-based graphs, with reactions as nodes connected if the metabolite produced by one is consumed by the other. TumorMet is deposited at figshare repository30 and the Met2Graph package used to generate it is available at the Met2Graph Github repository (https://github.com/cds-group/Met2Graph).MethodsThe metabolism involves several players, and focusing on one or another influences the type of analysis and the knowledge that can be extracted. The metabolites and the enzymes represent the main molecular components. A biochemical reaction is a transformation process that uses/consumes some metabolites (reactants) to produce new ones (products). The enzymes can facilitate these transformations as they are particular proteins having catalytic activity and the ability to speed up the rate of a reaction binding the substrate by a lock-key or induced-fit model. Not all the reactions are catalyzed by enzymes, as some of them can occur spontaneously. The enzymes are selective; this means that one binds specifically one or few substrates and, consequently, can catalyze one or more reactions, while the same reaction can be catalyzed by more enzymes acting as complex or as mutually exclusive catalyzers. This information is crucial in defining the rules to design a metabolic network since the connections between the metabolic players can be multiple and of different nature when involving the enzymes. In order to manage this issue, we defined some simplification strategies when enzymes represent edges and give rise to multiple connections (as in the case of Metabolites-based networks) and a different consideration of complex and mutually exclusive relationships when enzymes represent the nodes (as in the case of Enzymes-based networks). Further details are provided below in the network construction sections. The repository we provide contains different types of metabolic networks, depending on the nodes and the rules behind the connections: Metabolites-, Enzymes- and Reactions-based networks. A graphical overview of the metabolic networks construction is provided in Fig. 1.Fig. 1Overview of the Metabolic networks construction. The context-specific GSMs used in this study derive from the human generic GSM through the integration of tissue-specific multi-omics data (tissue-specific GSMs from Human Metabolic Atlas) or by integration of TCGA transcriptomics data (PDGSMMs from Biomodels). The context-specific GSMs carrying information about biochemical reactions are the input to create the context-specific metabolic networks of the TumorMet repository. Metabolites-based_tissue networks are generated by integrating TCGA gene/enzyme-expression data into the tissue-specific GSMs to weight the edges represented by enzymes connecting two metabolites. Networks of different patients have the same structure with different edge weights depending on patient expression profile. Enzymes-, Reactions and Metabolites-based_PDGSMMs networks are created from PDGSMMs and have enzymes/reactions as nodes connected by metabolites or metabolites as nodes connected by enzymes. Networks of different patients have different structures and no weights.Metabolic modelsTissue-specific GSMs for 5 of the different origin sites of cancer (lung, kidney, brain, ovary, prostate)31 and breast cancer INIT model32 were downloaded from the Metabolic Atlas repository (http://www.metabolicatlas.org) in the compressed Systems Biology Markup Language (SBML) format33 to create the Metabolites-based graphs. PDGSMMs from the Biomodels repository (https://www.ebi.ac.uk/biomodels/pdgsmm/) have been downloaded to generate Metabolites-, Enzymes- and Reactions-based_PDGSMMs graphs for each patient. The Gene-Protein-Reaction (GPR) relationships were extracted from version 1.4.1 of the human generic GSM (https://github.com/SysBioChalmers/Human-GEM/tree/master/model).Gene expression dataGene expression data from 6 different tumor primary sites were used to create context-specific Metabolites-based metabolic networks. FPKM (fragments per kilobase per million reads mapped) normalized and log-transformed read counts from RNA sequencing experiments of the breast (TCGA-BRCA), lung (TCGA-LUAD and TCGA-LUSC), kidney (TCGA-KIRC and TCGA-KIRP), brain (TCGA-GBM and TCGA-LGG), ovary (TCGA-OV), and prostate (TCGA-PRAD) cancers were obtained from the Genomic Data Commons (GDC) data portal (https://portal.gdc.cancer.gov). GDC includes several cancer projects, among which The Cancer Genome Atlas (TCGA), which we selected to download the data. Each of them represents a dataset of the repository. Clinical annotations of the samples were also extracted from the database and included in each dataset as sample-sheets.Metabolites-based_tissue networks constructionThe metabolites are the nodes of the network, labeled by the corresponding ID, connected if they are involved in the same reaction, one as a reactant and one as a product. The connections have been created using the information from the relative context-specific metabolic model. Recurrent metabolites (e.g, ATP, CO2, H2O) have been removed to avoid redundant connections and unrealistic definition of paths34. The small molecules such as H2O, NH3, O2, CO2, phosphate, and cofactors are generally considered recurrent metabolites. The recurrent metabolites list we used is provided as external data of the package Met2Graph; the argument rmMets can be set to FALSE to avoid removal, or the list can obviously be personalized by the user. The GPR associations have been derived from the generic human GSM. Each edge is labeled by the Ensembl stable ID (in the form of ENS[species prefix][feature type prefix][a unique eleven-digit number]) of the enzyme/s catalyzing the reaction, when present, and weighted by the expression value/s of the corresponding gene/s obtained by the GDC Portal. Each resulting graph corresponds to a specific sample of the GDC tumor dataset considered. These rules create graphs where a couple of nodes can have multiple edges since multiple enzymes are involved in the same reaction and/or because the same nodes pair can be present in different reactions. Multiple edges have been simplified by averaging the expression values of enzymes acting in the same reaction and then summing up these averages corresponding to different reactions with the same nodes pair. Thus, all the graphs resulting from the same metabolic model have the same number of nodes and edges but different edge weights. The networks are then personalized for each patient by using the expression values and as a consequence, the gene context mentioned by19 is met. Based on the rules defining the edges, these networks are directed. The properties of these networks are summarized in Table 1.Table 1Properties of the Metabolites-based networks derived from tissue models.KidneyLungBrainBreastOvaryProstate# Graphs92811357021217379551# Vertices403439903922339438273939# Edges921090588914654885338747Edge density0.0005660.000560.000580.000570.000580.00056Avg. network degree4.564.544.543.864.464.44Edge weights√√√√√√Assortativity degree−0.038−0.035−0.034−0.049−0.027−0.03Global transitivity0.120.130.130.0530.1350.132Avg. local transitivity0.140.140.150.130.140.15Minimum diameter134.75134.15143.41141.43131.15146.73Maximum diameter243.08206.47200.47236.41188.13225.4For each tissue dataset (along the columns), we report the number of graphs (first row) and the corresponding networks topological properties, such as the number of vertices and edges, edge density, average network degree, eventual presence of edge weights, assortativity degree, global transitivity, average local transitivity, minimum and maximum diameter (second through and eleventh rows). Observe that, for each tissue dataset, Metabolites-based networks share the same network structure, and thus topological properties, for all the samples since they derive from the same tissue metabolic model personalized with gene expression values.Metabolites-based_PDGSMMs networks constructionThe logic behind the generation of Metabolites-based_PDGSMMs networks is the same as that of the networks derived from tissue models described in the previous paragraph, with the only difference that here each patient-specific network is derived from the corresponding PDGSMM downloaded from the BioModels repository. The edges are weighted using the patient’s gene expression data from the GDC repository. Therefore, each patient-specific network has a different structure and different edge weights. These graphs are directed and weighted. The properties of these networks are summarized in Table 2a.Table 2For each tissue dataset of the Metabolites- (a), Enzymes- (b), and Reactions-based_PDGSMMs (c) networks (along the columns), we report the number of graphs (first row) and the corresponding networks topological properties, such as the number of vertices and edges, edge density, average network degree, eventual presence of edge weights, assortativity degree, global transitivity, average local transitivity, minimum and maximum diameter (second through and eleventh rows).KidneyLungBrainBreastOvaryProstate(a) Properties of the Metabolites-based networks from PDGSMMs# Graphs737829138920295470# Vertices2679.05 ± 316.112619.5 ± 310.492634.49 ± 277.22576 ± 303.922576.93 ± 307.472676.14 ± 300.88# Edges6121.64 ± 839.576008.53 ± 908.156074.34 ± 783.775870.26 ± 841.165729.2 ± 837.645799.38 ± 769.08Edge density0.00086 ± 0.0000090.0009 ± 0.00000090.0009 ± 0.00000090.0009 ± 0.00010.00087 ± 0.00010.0008 ± 0.000009Avg. network degree4.56 ± 0.234.57 ± 0.34.6 ± 0.284.54 ± 0.284.44 ± 0.294.33 ± 0.27Edge weights✓✓✓✓✓✓Assortativity degree−0.01 ± 0.020.006 ± 0.029−0.004 ± 0.0330.012 ± 0.031−0.008 ± 0.034−0.018 ± 0.03Global transitivity0.16 ± 0.020.17 ± 0.0260.16 ± 0.030.17 ± 0.0290.14 ± 0.0350.12 ± 0.035Avg. local transitivity0.12 ± 0.020.12 ± 0.020.12 ± 0.020.12 ± 0.020.11 ± 0.0210.11 ± 0.02Minimum diameter134.8134.01140.76118.44120.84145.97Maximum diameter302.7241.72217.06255.8211.24282.75(b) Properties of the Enzymes-based networks from PDGSMMs# Graphs737829138920295470# Vertices1941.256 ± 300.921859.76 ± 317.841911.35 ± 274.351846.58 ± 305.481859.98 ± 309.681934.25 ± 266.7# Edges63906.79 ± 18916.4959341.79 ± 20947.8863485.08 ± 17933.1959530.08 ± 19744.6759316.15 ± 202888.0663922 ± 16898.25Edge density0.016 ± 0.0020.016 ± 0.0020.07 ± 0.0020.016 ± 0.0020.016 ± 0.0020.016 ± 0.002Avg. network degree63.8 ± 14.1661.23 ± 16.364.63 ± 14.0762.15 ± 15.6761.39 ± 15.8664.49 ± 13.05Edge weightsxxxxxxAssortativity degree0.25 ± 0.040.25 ± 0.040.25 ± 0.040.26 ± 0.040.24 ± 0.0460.25 ± 0.038Global transitivity0.18 ± 0.040.19 ± 0.0460.18 ± 0.0390.19 ± 0.0460.19 ± 0.0470.18 ± 0.035Avg. local transitivity0.29 ± 0.0180.3 ± 0.020.03 ± 0.0180.3 ± 0.020.298 ± 0.020.29 ± 0.015Minimum diameter141313141314Maximum diameter343628333530(c) Properties of the Reactions-based networks from PDGSMMs# Graphs737829138920295470# Vertices3578.24 ± 595.0373511.46 ± 637.323560.4 ± 543.413431.28 ± 591.123327.51 ± 582.493398 ± 527.44# Edges54823.89 ± 16130.960808.68 ± 19146.2260137 ± 1774959467 ± 1833049776.08 ± 17158.546345.11 ± 14345.88Edge density0.0043 ± 0.00080.0048 ± 0.00070.004 ± 0.000850.0049 ± 0.000860.004 ± 0.000920.004 ± 0.0008Avg. network degree30.2 ± 6.1333.74 ± 7.0533.17 ± 7.133.9 ± 7.1729.3 ± 6.9626.91 ± 5.79Edge weightsxxxxxxAssortativity degree0.027 ± 0.0160.052 ± 0.180.023 ± 0.150.048 ± 0.170.065 ± 0.20.06 ± 0.17Global transitivity0.028 ± 0.0150.038 ± 0.0170.037 ± 0.0160.038 ± 0.0170.028 ± 0.0160.024 ± 0.013Avg. local transitivity0.038 ± 0.0060.04 ± 0.0060.043 ± 0.00590.04 ± 0.0060.04 ± 0.0060.04 ± 0.006Minimum diameter484848484849Maximum diameter10311397104102101Observe that each network derived from PDGSMMs and corresponding to each patient sample has a different structure since the starting models are patient-specific (see Paragraphs on Metabolites-, Enzymes-, and Reactions-based PDGSMM networks). Therefore, values for network properties are reported as average ± standard deviation across all the networks of each dataset.Enzymes-based_PDGSMMs networks constructionThese networks have enzymes as nodes connected if one catalyzes a reaction producing a metabolite consumed in a reaction catalyzed by the other. The recurring metabolites have also here been removed. According to the GPR, the enzymes involved in each reaction are associated by AND or OR logical relationship, indicating an enzymatic complex or an alternative activity, respectively. Based on this, enzymes related by AND have been considered as a single node, while OR relationships have been split into different nodes. To create patient-specific networks, PDGSMMs have been used as starting models for Metabolites-, Enzymes-, and Reactions-based_PDGSMMs datasets and downloaded from the BioModels repository. Each sample graph has then a different structure deriving from a different model. These graphs are directed and not weighted. The properties of these networks are summarized in Table 2b.Reactions-based_PDGSMMs networks constructionThe rules behind these networks are similar to those of Enzymes-based networks, with the difference of having reactions as nodes, connected if one produces a metabolite consumed by the other. Recurring metabolites have been removed as well. To have sample-specific graphs also in this case we used the PDGSMMs from Biomodels. The resulting graphs are unweighted and directed, and each sample has a different structure determined by the different starting models. The properties of these networks are summarized in Table 2c.Simplified networks constructionGiven the complexity and the size of these networks, we also provided a set of Metabolites-based sub-networks of a subset of kidney and lung samples, simplified according to the approach described in21. Briefly, central nodes have been selected by the Eigen centrality score, a measure describing the importance of a node in a graph that depends on that of its neighbors. The classification tests performed to demonstrate the reliability of these sub-networks compared to the whole networks gave comparable accuracy results (see Tables 3 and 4 in21). For each tissue, two sets of networks with a different number (#) of resulting nodes are provided. The properties of these networks, forming the Simpl-Kidney-# and Simpl-Lung-# datasets, are summarized in Tables 3 and 4.Table 3Properties of the Simplified Networks. See the caption of Table 1 for details.Simpl-Kidney-441Simpl-Kidney-1034Simpl-Lung-312Simpl-Lung-1017# Graphs299299337337# Vertices44110343121017# Edges1585322610903102Edge density0.01630.0060.0220.006Avg. network degree7.186.246.986.1Edge weights✓✓✓✓Assortativity degree−0.22−0.13−0.11−0.12Global transitivity0.30.210.450.23Avg. local transitivity0.230.220.290.22Minimum diameter15.52125.9916.8879.7Maximum diameter39.37455.3632.14267.6Table 4Classes per dataset for usage validation of Metabolites-based networks through classification. Only primary tumors have been selected.KidneyLungBrainCases822Cases1025Cases666Kidney Renal Papillary cell carcinoma (KIRP)288Adenocarcinoma (LUAD)524Glioblastoma multiforme (GBM)155Kidney Renal Clear cell carcinoma (KIRC)534Squamous cell carcinoma (LUSC)501Lower grade glioma (LGG)511BreastOvaryProstateCases1085Cases290Cases497High-grade serous ovarian cancers subtypes40Gleason scoreBasal-like192Differentiated75Pattern 3199HER2-enriched82Mesenchymal75Pattern 4249Luminal A564Proliferative75Pattern 549Luminal B207Immunoreactive65Normal-like40ClassificationMetabolites-based_tissue datasetsIn previous works, we have demonstrated the utility of the network datasets in classification and clustering tasks using subsets of some of the Metabolites-based graph datasets now included in the TumorMet repository20,21,35–37. Here, we extend to the entire repository the usage validation introduced in20, wherein we classify whole graphs sharing the same set of nodes. The basic idea is to 1) represent each graph of a dataset using probability distributions describing the topological properties of each node; 2) extract the distance matrix (Gram matrix), i.e., the symmetric square matrix containing the distances, taken pairwise, between the networks of the dataset; and 3) classify the networks based on the obtained distance vectors.Based on the performance results achieved in20,21,35–37, here we selected the Transition Matrix of order one \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathscr{T}}}^{r}$$\end{document}Tr for representing each graph \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathscr{G}}}^{r}$$\end{document}Gr, whose generic element \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathscr{T}}}_{i,j}^{r}$$\end{document}Ti,jr is the probability of a node i to be reached in one step by a random walker located in node j. Each row \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathscr{T}}}_{i}^{r}$$\end{document}Tir of this matrix includes local information on the connectivity of node i.For computing the distance between two networks \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathscr{G}}}^{p}$$\end{document}Gp and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathscr{G}}}^{q}$$\end{document}Gq, we selected the network distance:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr{M}}({{\mathscr{G}}}^{p},{{\mathscr{G}}}^{q})=\frac{1}{l}{\sum }_{i=1}^{l}{d}_{JS}({{\mathscr{T}}}_{i}^{p},{{\mathscr{T}}}_{i}^{q}),$$\end{document}M(Gp,Gq)=1l∑i=1ldJS(Tip,Tiq),obtained by averaging over all the l graph nodes the Jensen-Shannon distances dJS of the probability distributions of their nodes38.For classification, we considered the primary tumor classes described in Table 6. In particular, for Kidney, Lung, and Brain, the Primary-Tumor diagnoses indicated in the GDC sample metadata file, downloaded along with the gene expression files, have been used to label the samples and fulfill the classification task. For Breast, the 5 subtypes have been derived from the PAM50 classification39. As the Normal-like subtype has only 40 samples and is very similar to the Luminal A subtype, we performed the tests both including (Breast_5cl) and excluding (Breast_4cl) this class. For Prostate, as having only one class of diagnosis, the Gleason pattern score, an indicator of different grades of malignancy, has been used. Among the possible four classes (Pattern from 2 to 5), we excluded the Pattern 2 class (not shown in Table 6), as it is made of only one sample. Moreover, we considered two different classification problems: the Prostate1 case, that aims at discriminating the Pattern 3 samples (199) from the Pattern 4 ones (249); and the Prostate2 case, that consists in discriminating the Pattern 3 samples from the samples being assigned to Pattern either 4 or 5 (289). For Ovary, the subtype assignment of High-Grade Serous Ovarian Cancer (HGSOC) has been taken from40.Metabolites-, Enzymes-, and Reactions-based_PDGSMMs datasetsThe graph2vec framework41 is a neural method for learning graph-level embeddings in an unsupervised manner. It describes nodes through a recursive node relabeling algorithm assigning to each node a label uniquely representing its rooted subgraph (neighborhood). These labels form a vocabulary of words, and graphs are represented in the form of documents. Then, the Distributed Bag of Words doc2vec approach42 is used to learn the graph (document) embeddings. The performance has been evaluated by means of a stratified 10-fold Cross-Validation (CV) in which a SVM classifier, with a linear kernel, was applied to train and make predictions on 64-sized vectorizations of graphs (embeddings) produced by graph2vec with a recursive depth of 3 and a training duration of 200 epochs. The class labels used for the classification task are specified in Table 5.Table 5Classes of PDGSMMs used to accomplish the classification task of Kidney and Lung PDGSMMs derived networks.KidneyLungCases737829ClassesKIRC 484LUAD 429KIRP 253LUSC 400Data RecordsThe network files and associated metadata composing the repository TumorMet are available at figshare repository30. The file TumorMet-repository.pdf summarizes the content of the repository. For easy access to the files, the repository is organized into seven datasets, each in a separate folder, representing the six tumor tissues and the simplified networks (i.e., Prostate, Lung, Kidney, Breast, Ovary, Brain, and Simplified networks). In each main tissue dataset folder, the sample-sheet file reporting the sample metadata as downloaded from GDC (i.e. Sample sheet.tsv) and an excel file reporting the correspondences between PDGSMM ids and TCGA ids (Dictionary_ids.xlsx) are provided. Each tissue dataset folder contains subfolders for the different types of networks, namely Metabolites-, Enzymes-, and Reactions-based, compressed in.zip format. The Metabolites-based folder is further subdivided into folders containing the Metabolites-based networks deriving from tissue models (Metabolites-based_tissue) and BioModels PDGSMMs (Metabolites-based_PDGSMMs). Enzymes- and Reactions-based networks are only derived from PDGSMMs. Simplified networks are provided for Kidney and Lung tissues. Each tissue folder contains the sample-sheet file reporting the sample metadata as downloaded from GDC (i.e., Sample sheet.tsv) and two subfolders for the networks files based on the number of nodes retained after the simplification process (for Kidney eigen_simplified_441_nodes and eigen_simplified_1034_nodes; for Lung eigen_simplified_312_nodes and eigen_simplified_1017_nodes). All the network files are provided in GraphML format. GraphML is a flexible and convenient XML format for storing network information. It supports unweighted, weighted, undirected, and directed networks and allows for the definition of node and edge attributes (http://graphml.graphdrawing.org/). A scheme of the repository content is illustrated in Fig. 2, while a summary of the networks features in terms of starting material and number of networks is provided in Table 6.Fig. 2Scheme of the content of the TumorMet repository.Table 6Networks provided in the TumorMet repository.Type of networkData used to build the networksNumber of networksKidneyMetabolites-based_tissue• Tissue-Specific Model - Kidney• TCGA-KIRC & TCGA-KIRP GE928:607 TCGA-KIRC321 TCGA-KIRPMetabolites-, Enzymes-, Reactions-based_PDGSMMs• PDGSMMs from TCGA-KIRC & TCGA-KIRP• TCGA-KIRC & TCGA-KIRP GE (only for Metabolites-based)737:484 TCGA-KIRC253 TCGA-KIRPSimplified• Tissue-Specific Model - Kidney• TCGA-KIRC & TCGA-KIRP GE299 for each simplification:193 TCGA-KIRC106 TCGA-KIRPLungMetabolites-based_tissue• Tissue-Specific Model - Lung• TCGA-LUAD & TCGA-LUSC GE1135:585 TCGA-KIRC550 TCGA-KIRPMetabolites-, Enzymes-, Reactions-based_PDGSMMs• PDGSMMs from TCGA-LUAD & TCGA-LUSC• TCGA-LUAD & TCGA-LUSC GE (only for Metabolites-based)829:429 TCGA-LUAD400 TCGA-LUSCSimplified• Tissue-Specific Model - Lung• TCGA-LUAD & TCGA-LUSC GE337 for each simplification:174 TCGA-LUAD163-TCGA-LUSCBrainMetabolites-based_tissue• Tissue-Specific Model - Brain• TCGA-GBM & TCGA-LGG GE702:173 TCGA-GBM529 TCGA-LGGMetabolites-, Enzymes-, Reactions-based_PDGSMMs• PDGSMMs from TCGA-GBM• TCGA-GBM GE (only for Metabolites-based)138 TCGA-GBMBreastMetabolites-based_tissue• INIT Cancer Model - Breast TCGA-BRCA GE1217 TCGA-BRCAMetabolites-, Enzymes-, Reactions-based_PDGSMMs• PDGSMMs from TCGA-BRCA• TCGA-BRCA GE (only for Metabolites-based)920 TCGA-BRCAOvaryMetabolites-based_tissue• Tissue-Specific Model - Ovary379 TCGA-OV• TCGA-OV GEMetabolites-, Enzymes-, Reactions-based_PDGSMMs• PDGSMMs from TCGA-OV• TCGA-OV GE (only for Metabolites-based)295 TCGA-OVProstateMetabolites-based_tissue• Tissue-Specific Model - Prostate• TCGA-PRAD GE551 TCGA-PRADMetabolites-, Enzymes-, Reactions-based_PDGSMMs• PDGSMMs from TCGA-PRAD• TCGA-PRAD GE (only for Metabolites-based)470 TCGA-PRADFor each tumor tissue: the type of networks, the data used to generate the networks in terms of metabolic models and Gene Expression (GE) data from TCGA projects, and the number of networks, eventually subdivided by TCGA project ID. Observe that in the case of PDGSMMs derived networks, only for Metabolites-based_PDGSMM networks the GE data have been used to weight the edges.Technical ValidationOur validation process consisted of data-type and structural validation, as well as usage validation through downstream applications.Data-type and structural validationThe quality of the original data used to generate the networks is given by the reliability of the data sources repositories, i.e., GDC, Human Metabolic Atlas, and BioModels. Node IDs were verified to be of the same type. All edges were verified to be between nodes in the node list. All attribute data were verified to correspond to an existing node or edge. The structural integrity of the networks has been assessed by removing self-loops. Any duplicate edges were also removed. We further checked that nodes with no edges were not present in the networks.Usage validationThe tumor metabolic networks can be exploited in several downstream applications, ranging from pure network analysis to multi-level integration with other biological networks or data, to machine and deep learning approaches for unraveling the complex metabolic machinery and its role in precision medicine. In this section, we show the usage of TumorMet networks in classification of tumor samples, thus giving an idea of one of their potential applications. To furnish a baseline for comparing methods and approaches, we give several details of the two different workflows used for Metabolites-based networks derived from tissue models and Metabolites-, Enzymes-, Reactions-based networks derived from PDGSMMs.Metabolites-based_tissue datasetsFor the evaluation of classification performance, i) each of the Metabolites-based datasets was subdivided into a training and a test set; ii) a statistical validation was obtained on the training sets using a 10-fold CV, to ensure that the results were not biased to a specific training subset; iii) finally, the classification performance on the test datasets was evaluated using the models built on the training datasets.i).In the case of Kidney, Lung, Breast, and Brain tissue datasets, the choice of the training sets was driven by our previous work36, where subsets of these datasets were already adopted for classification. Therefore, those subsets have been adopted here as training sets, while the newly added samples were assigned to the test sets. For the tissues not used previously (Ovary and Prostate), we obtained the training and test sets by using a 70:30 split ratio. The sample partitioning for each tissue is reported in Supplementary Table 1, while Figs. 3–4 provide the t-distributed Stochastic Neighbor Embedding (t-SNE) plots for the test sets.Fig. 3t-SNE representations of the Gram matrices of the test sets of the Kidney (a), Lung (b), Brain (c), and Ovary (d) Metabolites-based_tissue datasets. The TSNE function of the sklearn.manifold library has been used to generate the plots.Fig. 4t-SNE representations of the Gram matrices of the test sets of the Breast_4cl (a), Breast_5cl (b), Prostate1 (c), and Prostate2 (d) Metabolites-based_tissue datasets. The TSNE function of the sklearn.manifold library has been used to generate the plots.ii).For the statistical validation on the training sets, the data were min-max normalized and a Support Vector Machine (SVM) classifier with linear kernel was adopted using the libsvm implementation43 available in scikit-learn44. The one-vs.-rest strategy was used to classify the multi-class datasets. To account for unbalanced datasets, the “balanced” mode in sklearn was used to set the class weights; this parameter penalizes the wrong prediction of the classes having a number of instances lower than the others. The 10-fold CV on the training datasets was repeated 10 times, and the average of the CV scores are reported in Table 9 (top); these scores are also shown in the form of box plots in Fig. 5.Fig. 5Classification scores on the Metabolites-based_tissue datasets. The box-plots show the classification scores obtained from the 10 iterations of the evaluation procedure on the training sets of the six Metabolites-based_tissue datasets. (a–c) report Accuracy, Precision, Recall, and F1 as percentages; (d) reports MCC values.iii).The classification performance on the test sets was computed using the same SVM classifier learned on the training sets. The obtained results are reported in Table 9 (bottom). Kidney, Lung and Brain graphs are well classified, as shown by accuracy scores both in CV on training sets and using new samples as testing data (Table 9 and Figs. 3, 5). More challenging tasks are instead given by the classification of Breast, Ovary and Prostate samples.Regarding Breast, the inclusion of the Normal-like subtype into the classification does not dramatically change the results; however, compared to the tissues mentioned above, the results are worse, having an accuracy of around 80%. Looking at the t-SNE plots (Fig. 4a,b), it is evident how the Basal is the best discriminated and most homogeneous subtype, while some samples of Luminal A, Luminal B, and Her2 are overlapped, especially the latter two. Normal-like samples, as expected, are difficult to separate from Luminal A ones. Ovary samples are completely overlapping (Fig. 3d) and lead to poor accuracy percentage (around 70%, as reported in Table 9). Finally, the CV scores reported in Table 9 (top) and plotted in Fig. 5c, as well as the test samples validation results reported in Table 9 (bottom), indicate that Prostate samples are generally poorly discriminated and the results are slightly better for the Prostate2 classification task (when the Gleason Pattern 5 is assimilated to Pattern 4). Prostate cancer is characterized by a high molecular heterogeneity45 which is evidently not caught considering only the Gleason score, as also highlighted by the t-SNE plots reported in Fig. 4c,d.Metabolites-, Enzymes-, Reactions-based_PDGSMMs datasetsAs detailed in the Section on Metabolic networks construction, these PDGSMMs derived graphs differ from the Metabolites-based graphs in that they do not share a common set of nodes across all patients. Therefore, we decided to accomplish the classification task on these datasets through a whole-graph embedding framework. Classification results based on these embeddings using the class labels specified in Table 5 for the Kidney and Lung PDGSMMs derived network datasets are reported in Table 8.It is evident that the performance for these types of networks is not as good as the one obtained with Metabolites-based graphs, but it is worth pointing out that the two approaches to the classification task are completely different due to the different nature of the networks. Enzymes- and Reactions-based networks are indeed not weighted and have different structures being generated from different models. The complexity and density of these networks surely require a deeper investigation of the best suitable approach and parameters tuning to discriminate the differences among the samples, which is not the aim of this paper. As mentioned previously, one of the interesting aspects of the metabolic networks is their plasticity since different types of graphs can be generated depending on the desired nodes and connections. In future work, we will consider generating unique tri-partite graph for each patient to investigate the possibility to reduce classification performance differences. As for the networks extracted from tissue-specific models, the Metabolites-based_PDGSMMs networks are weighted by gene expression values. Comparing weighted vs. non-weighted networks in terms of classification performance, it is evident that the weights do not add any crucial information for discriminating the classes (Table 9). These networks derive from PDGSMMs reconstructed through the tINIT algorithm integrating TCGA gene expression data. Adding expression values to edges is therefore redundant and likely the models are already well contextualized. Instead, the weights have a different role in Metabolites-based_tissue networks, where are crucial for personalizing the networks in terms of patients. Furthermore, even if tested with different methods, the patients-specific Metabolites-based networks derived from tissue models seem to well contextualize the tissue models in terms of patients resulting as more representative of the tumor classes and with a higher discriminative power, as highlighted by classification performances (Table 7).Table 7Classification scores on Metabolites-based_tissue datasets.# ClassesKidneyLungBrainBreast_4clBreast_5clOvaryProstate1Prostate222245422Cross-validation on training sets# Samples per class159/90158/150109/358135/58/395/145135/58/395/145/2853/46/53/53140/172140/209Accuracy avg %92.80 ± 4.8794.87 ± 3.6895.83 ± 2.6584.91 ± 4.1581.02 ± 4.2979.78 ± 7.7971.83 ± 8.1775.086.17Precision avg %91.97 ± 5.594.94 ± 3.8593.63 ± 4.6981.60 ± 4.9972.30 ± 6.679.83 ± 8.5772.23 ± 8.3274.84 ± 6.11Recall avg %92.99 ± 5.194.95 ± 3.5495.23 ± 3.5785.93 ± 5.3678.55 ± 7.379.93 ± 8.8672.25 ± 8.2275.99 ± 6.39F1 avg %92.12 ± 5.394.74 ± 3.8394.15 ± 3.7682.66 ± 4.7673.36 ± 6.1178.09 ± 8.8371.14 ± 8.4174.31 ± 6.35MCC avg0.85 ± 0.10.90 ± 0.070.89 ± 0.070.77 ± 0.060.73 ± 0.060.73 ± 0.10.44 ± 0.160.51 ± 0.12Test samples validation# Samples per class375/198366/35146/51157/24/169/6257/24/169/62/1222/19/22/2259/7759/89Accuracy %97.0393.7291.0085.2683.6470.5973.5373.00Precision %96.4093.7285.9280.6274.4673.4973.6372.60Recall %97.1393.7293.3683.9981.8370.3374.0573.53F1%96.7593.7288.5682.0277.4871.0573.4472.57MCC0.940.870.790.770.760.610.480.46Top: CV on training sets; Bottom: Validation on test sets.Table 8Classification scores on Enzymes- and Reactions-based_PDGSMMs Kidney and Lung datasets.KidneyLungEnzymes-based_PDGSMMsReactions-based_PDGSMMsEnzymes-based_PDGSMMsReactions-based_PDGSMMs# Classes2222# Samples per class484/253484/253429/400429/400Accuracy avg %78.97 ± 5.1583.44 ± 4.3278.17 ± 2.8977.93 ± 2.62Precision avg %77.35 ± 6.0081.99 ± 4.9178.57 ± 2.6978.36 ± 2.44Recall avg %75.72 ± 5.3681.16 ± 5.0078.00 ± 3.0477.83 ± 2.75F1 avg %76.18 ± 5.4281.39 ± 4.9478.05 ± 3.0477.76 ± 2.75MCC avg0.53 ± 0.110.63 ± 0.100.57 ± 0.060.56 ± 0.05Table 9Classification scores on weighted and unweighted Metabolites-based_PDGSMMs networks of Kidney samples.Metabolites-based_PDGSMMs Kidney# Classes2# Samples per class484/253weightedunweightedAccuracy avg %83.45 ± 4.5885.48 ± 3.12Precision avg %82.28 ± 5.1384.43 ± 3.71Recall avg %80.87 ± 4.8282.99 ± 3.36F1 avg %81.32 ± 5.0883.59 ± 3.47MCC avg0.63 ± 0.100.67 ± 0.07Usage NotesThe networks presented here have been generated using the Met2Graph R package we developed (see the paragraph on “Code availability”). The model in SBML format is imported and read by the Met2Graph package through the function readSBMLmod from the sybilSBML46 package. Several checkpoints are included in the function to validate the model object before importing it, such as check of upper and lower bounds, GPR mapping, reactions’ ids, and presence of list of reactants and products. The code snippets of Listings 1–4 show Met2Graph functions and arguments used to obtain the different networks:Listing 1 Metabolites-based_tissue networks.Listing 2 Metabolites-based_PDGSMMs networks.Listing 3 Enzymes-based_PDGSMMs networks.Listing 4 Reactions-based_PDGSMMs networks.There are several open-source network libraries that can be used to analyze and visualize the networks provided in GraphML format. Examples of network analysis and visualization software include NetworkX, igraph, Cytoscape, yEd and Gephi.Supplementary information
Supplementary Table 1