Leveraging high-throughput data and unsupervised learning to characterize cancer metabolic heterogeneity

University of Milano - Bicocca Department of Informatics, Systems and Communication Master’s degree program in Data Science Leveraging high-throughput data and unsupervised learning to characterize cancer metabolic heterogeneity Supervisor: Prof. Dario Pescini Co-supervisors: Prof.ssa Chiara Damiani, Dott. Bruno Giovanni Galuzzi Master’s degree thesis by: Luca Milazzo ID number 856702 Academic Year 2023-2024 Contents Abstract 1 Introduction 3 1 Scientific background 7 1.1 Metabolic Networks . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Linear programming . . . . . . . . . . . . . . . . . . . . . . . 9 1.3 Constraint-based modelling . . . . . . . . . . . . . . . . . . . 11 1.3.1 FBA............................ 11 1.3.2 pFBA ........................... 11 1.3.3 FVA............................ 12 1.3.4 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . 13 1.4 FluxSampling .......................... 14 1.4.1 Markov Chains . . . . . . . . . . . . . . . . . . . . . . 16 1.4.2 Internal Sampling of the Polytope . . . . . . . . . . . 17 1.4.3 Corner-Based Sampling (CBS) . . . . . . . . . . . . . 19 1.4.4 False discoveries and sampling best practises . . . . . 20 1.4.5 Samples statistics . . . . . . . . . . . . . . . . . . . . . 20 1.5 Omics data integration . . . . . . . . . . . . . . . . . . . . . . 23 2 Methods 27 2.1 Cancer Cell Line Encyclopedia . . . . . . . . . . . . . . . . . 27 2.2 ENGRO2............................. 30 2.3 Dataintegration ......................... 31 2.3.1 Cell-specific metabolic models . . . . . . . . . . . . . . 31 2.3.2 scFEA........................... 35 2.4 Preprocessing........................... 41 2.4.1 Transcriptomics denoising: MAGIC . . . . . . . . . . 41 2.4.2 Scaling .......................... 47 2.4.3 Principal Component Analysis (PCA) . . . . . . . . . 48 2.5 Clustering Analysis . . . . . . . . . . . . . . . . . . . . . . . . 50 2.5.1 Hyperparameters tuning . . . . . . . . . . . . . . . . . 50 2.5.2 Clusterings comparison . . . . . . . . . . . . . . . . . 54 ii 2.5.3 Internal Evaluation Metrics . . . . . . . . . . . . . . . 56 2.5.4 External Evaluation Metrics . . . . . . . . . . . . . . . 59 2.5.5 Biomarkers ........................ 62 3 Results and discussion 63 3.1 Metabolomics as clustering benchmark . . . . . . . . . . . . . 63 3.2 Transcriptomics and metabolomics share similar metabolic phenotypes ............................ 66 3.3 INTEGRATE fluxomics does not reflect metabolic profiles . . 68 3.4 MAGIC mitigates RNA-seq noise . . . . . . . . . . . . . . . . 71 3.5 MAGIC improves fluxomics prediction . . . . . . . . . . . . . 77 3.6 CHRR outperforms CBS in samples clustering . . . . . . . . 83 3.7 scFEA as promising integration framework . . . . . . . . . . 88 4 Conclusions 93 A Supplementary figures and tables 105 A.1 PCA preprocessing . . . . . . . . . . . . . . . . . . . . . . . . 105 A.2 Metabolomics biomarkers and meta-data . . . . . . . . . . . . 108 A.3 Clustering hyperparameters . . . . . . . . . . . . . . . . . . . 112 iii Abstract Diseases such as obesity, diabetes and cancer are influenced by various factors, including genetics and environmental conditions. Recent research indicates that metabolic alterations play a significant role in the development and progression of these diseases [ 1 , 2 , 3 , 4 ]. However, directly measuring metabolic fluxes remains challenging due to technical and financial constraints. To address this problem, genome-scale metabolic models (GEMs) provide a comprehensive computational framework for predicting reaction fluxes. These models use techniques such as linear programming to numerically simulate metabolism. Consequently, integration methods are employed to incorporate high-throughput omics data within metabolic models to predict fluxomics with higher precision. Unfortunately, their evaluation is limited to just a few fluxes that can be directly measured in the laboratory. This thesis aims to overcome the limitations of benchmarking integration methods by using intracellular metabolomics as a reference. Two integration methods were applied to 513 different cancer cells from the Cancer Cell Line Encyclopedia (CCLE) [ 5 ]: INTEGRATE [ 6 ], a constraint-based steady-state method, and scFEA [ 7 ], a novel algorithm based on artificial neural networks. Using unsupervised learning techniques, clusters derived from metabolomics, transcriptomics, and fluxomics were compared to evaluate their concordance in recognizing metabolic phenotypes. A high level of coherence was found between metabolomics and transcriptomics, identifying similar metabolic subpopulations and suggesting the robustness of metabolomics as a benchmark. The results highlighted the necessity, even for bulk samples, of preprocessing the RNA-seq matrix with a denoising algorithm, specifically MAGIC [ 8 ]. MAGIC significantly decreased matrix sparsity by revealing lost genetic information and increased the concordance level between metabolomics and inferred INTEGRATE fluxomics. Conversely, scFEA inherently managed to reduce noise and achieved similar performance to the INTEGRATE framework. Overall, neither INTEGRATE nor scFEA drastically outperformed the other in comparing their metabolic clusters with the benchmark. Eventually, clusters identified among metabolomics, transcriptomics and predicted fluxomics can be analyzed by life scientists to extract valuable insights about the biomarkers of the subpopulations. 1 Introduction Many physiological and pathological conditions, including cancer, neurodegenerative diseases and aging, exhibit distinct metabolic features. Metabolism is intricately linked to nearly all cellular processes, making it a key indicator of a cell’s or organism’s physiological state [ 9 ]. Genome-scale metabolic models (GEMs), such as Recon [ 10 ], encompass all cellular reactions and their connections, allowing for the study of cellular metabolism and physiology. Although GEMs have provided a solid understanding of the fundamental structure of metabolic pathways, their systemic regulation remains elusive due to technical challenges. Metabolism can be fully characterized by reaction fluxes, which represent the rates at which metabolites are consumed in biochemical reactions. Directly measuring metabolic fluxes using labeled substrates is more complex compared to other omics technologies, such as metabolomics and transcriptomics [ 11 ]. This process is expensive, labor-intensive, and difficult to apply at single-cell or spatial resolution. Consequently, metabolic fluxes can currently only be inferred indirectly through experimental techniques and only for a limited set of reactions. To overcome these limitations, significant efforts have been made to predict metabolic fluxes from high-throughput data, particularly gene expression data from both bulk and single-cell analyses. Transcriptomics technologies, such as RNA-seq [ 12 ], measure gene expression at either the single-cell or bulk level and are typically used to predict metabolic fluxes, as gene expression regulates the production of enzymes that catalyze specific cellular reactions. However, cell metabolism is also regulated by intracellular metabolomics, which involves measuring the concentrations of metabolites, the reactants in these reactions, within the cell. Despite the availability of extensive transcriptomics, proteomics, and metabolomics datasets, a complete understanding of the regulatory mechanisms controlling metabolism requires appropriate integration of these data into predictive models. Constraint-based steady-state models offer a valuable framework for predicting metabolic fluxes using high-throughput omics data [ 13 , 14 ]. These models are analyzed using the constraint-based modeling (CBM) paradigm, which constructs a mathematical model by constraining network flows through a set of equations. The primary constraint is the steady-state constraint, which assumes that the rates of metabolite consump3 tion and production are balanced. The core equation of these models is Sv =  0 , where S is the stoichiometric matrix of size m×r (with m metabolites and r reactions) and v represents the unknown fluxes. To predict model fluxes, methods such as Flux Balance Analysis (FBA) and Flux Sampling (FS) are employed. FBA solves a linear programming problem to optimize the network for a specific objective function, which represents the cell’s metabolic goal, although assuming a single objective can be quite restrictive. Flux Sampling predicts fluxes by finding feasible solutions within the convex polytope defined by the constraints using Monte Carlo techniques, such as Coordinate Hit and Run with Rounding (CHRR), which samples the entire feasible region, or Corner-Based Sampling (CBS), which samples only the polytope’s boundary based on random objective coefficients. Techniques such as INTEGRATE [ 6 ], iMAT [ 15 ], and FASTCORE [ 16 ] have been developed to integrate transcriptomics data into GEMs within the CBM framework through Gene-Protein-Reaction (GPR) associations, thereby enhancing the characterization of cellular metabolism. GPRs are logical expressions (e.g., AND, OR rules) that capture the interplay between genetic information, the encoded proteins, and the biochemical reactions they catalyze. Novel methods like scFEA [ 7 ], based on artificial neural networks, aim to overcome the limitations of CBM models (e.g., linearity and steady-state assumptions) by directly inferring metabolic fluxes from transcriptomics data and GPRs without the need for a constraint-based framework. Given the technical difficulties associated with directly measuring metabolic fluxes, validating integration methods has been limited to a restricted set of fluxes that can be compared to experimental data. Consequently, benchmarking such models remains an open challenge. The aim of this thesis is to establish an evaluation pipeline for integration methods to address the benchmarking problem and to leverage high-throughput data for characterizing metabolic phenotypes for further biological investigations. The analysis is based on the assumption that even if a specific enzyme is overexpressed in a cell, the absence of the associated metabolite, the reactant, reduces the overall reaction fluxes. Therefore, intracellular metabolomics can be treated as a proxy to identify metabolic similarities among cells. Moreover, it is the only omic that is not exploited to predict fluxomics by either INTEGRATE or scFEA. For these reasons, intracellular metabolomics was selected as the benchmark for comparing, through unsupervised learning techniques, the cancer phenotypes identified by INTEGRATE, and its possible technical variants, and scFEA. The former employs a constraint-based steady-state framework, while the latter uses a novel computational approach based on artificial neural networks. The second objective of the thesis is to identify optimal clustering solutions for each biological data source to assign cluster labels to each omic dataset, facilitating further investigations into biomarkers for these subpopulations. Life scientists can analyze biomarkers for each metabolic phenotype, 4 1.3 Constraint-based modelling 1.3.1 FBA Flux Balance Analysis (FBA) is a widely used linear optimization method for analyzing genome-scale metabolic models. These networks encompass all the metabolic reactions of an organism and the genes encoding the corresponding enzymes [ 26 ]. The objective of FBA is to calculate the rates, or fluxes, of each reaction in the network by extracting the information encoded within the model. To achieve this, it is necessary to translate the model into a mathematical structure: the stoichiometric matrix S (see step 2 in Figure 1.3). The matrix has dimensions m×r , where m is the number of metabolites and r is the number of reactions. Each column represents a reaction and contains the stoichiometric coefficients of the participating metabolites [ 17 ]. The fluxes throughout the entire model are represented by a vector v , and the set of v such that Sv =  0 , i.e., the null space of the matrix, expresses all feasible fluxes of the model. This representation imposes mass balance constraints on the reaction fluxes, assuming a steady state, which refers to the particular state in which metabolite concentrations are constant, balancing production and consumption. Given the solution space, i.e. a polytope identified by the imposed constraints, FBA defines an optimization problem on the fluxes to study the model’s behavior under a specific objective function. FBA aims to maximize or minimize the objective function Z = cTv , where c is a weight vector indicating how much a reaction should contribute to the calculation of Z . Often, the objective function corresponds to a single flux to be optimized, vi , resulting in the following linear programming problem: max /min vi(1.2) S·v = 0 vL≤v ≤vU, where vL and vU represent the lower and upper bounds for the model fluxes. FBA also allows for the comparison of different experimental conditions with varying model constraints. FBA provides a single solution to the optimization problem by visiting the vertices of the feasible region (linear programming), but alternative solutions may exist. 1.3.2 pFBA Parsimonious Flux Balance Analysis (pFBA) is an extension of standard FBA that incorporates an additional objective: minimizing the total flux through the network. While FBA determines the optimal flux distribution to maximize or minimize a given objective function, such as biomass production, it often results in multiple equally optimal solutions due to the redundancy 11 and flexibility inherent in metabolic networks. pFBA addresses this by introducing the principle of parsimony, which assumes that biological systems are likely to minimize the use of resources, leading to the selection of a solution that requires the least overall enzyme usage. The pFBA method operates in two main steps: FBA Optimization and Secondary Optimization for Parsimony. In the first step, standard FBA is performed to determine the maximum or minimum value of the objective function, such as the biomass production rate Zmax: max Z=cTv (1.3) S·v = 0 vL≤v ≤vU, where S is the stoichiometric matrix, v is the flux vector, and c is the coefficient vector for the objective function. Once the optimal value Zmax is determined, a secondary optimization is performed to minimize the total flux through the network while still achieving the objective value found in the first step: min X i|vi|(1.4) S·v = 0 vL≤v ≤vU cTv =Zmax, where the objective is to minimize the sum of absolute flux values, ensuring that the solution is parsimonious by reducing the overall flux through the network. By incorporating this secondary optimization, pFBA not only finds a feasible solution but also ensures that the solution is biologically plausible, reflecting the likely preference of organisms to minimize metabolic costs. 1.3.3 FVA The Flux Variability Analysis (FVA) is used to calculate the optimality ranges of each flux, defining their limits. FVA solves the following two optimization problems, one of minimization and one of maximization, for each desired flux vj, with j= 1, . . . , r: max /min vj(1.5) S·v = 0 vL≤v ≤vU vj≥γZ0, where Z0 is an optimal solution for 1.2 and γ is the parameter used for performing FVA referring to a suboptimal (0 ≤γ < 1) or optimal ( γ = 1) state. 12 Figure 1.3: The steps of FBA. 1.3.4 Sensitivity analysis Mathematical models, such as constraint-based metabolic networks, can be highly complex, making the relationships between inputs and outputs difficult to discern. Consequently, sensitivity analysis is employed in the context of metabolic models to assess the influence of alterations in constraints on a predefined output feature. For instance, the deletion, or knockout, of a model reaction from the stoichiometric matrix S is equivalent to setting its upper bound to zero. By doing this, it becomes possible to analyze the effect of such deletions by performing Flux Balance Analysis (FBA) using a key reaction, such as Biomass, as the objective function. The sensitivity coefficient of reaction j indicates the fraction of preserved Biomass flux after the deletion of reaction jand it is defined as follows: Sensitivityj=ZKO j ZBiomass (1.6) where ZKO j is the Biomass optimized flux after the deletion of reaction j and ZBiomass is the Biomass optimized flux of the original model. This procedure is repeated for all model reactions, resulting in an array of sensitivity coefficients of length r, where rrepresents the number of model reactions. 13 1.4 Flux Sampling Flux Balance Analysis (FBA) imposes an objective function to explore metabolic models, which introduces a bias because it assumes a specific cellular goal. For example, studying how a bacterium’s metabolism changes under high temperature versus normal conditions requires identifying an objective to optimize, which is often unknown and can lead to inaccurate conclusions if chosen incorrectly. On the other hand, FVA returns only the range of each flux, while it does nor provide specific flux distributions. Flux sampling (FS) is used to uniformly sample the entire feasible region, unlike FBA which focuses on a single vertix, creating a set of feasible solutions. This sampling provides a finite representation of the metabolic network, allowing the extraction of key insights about the entire model. As compared to FVA and FBA, FS also offers information on the range and probability distributions of fluxes. There are two formulations for each metabolic model in the context of flux sampling: deterministic, used in this thesis, and stochastic. The deterministic approach follows FBA, applying steady-state constraints, while the stochastic formulation relaxes these constraints, allowing for measurement errors. Figure 1.4 shows how these methods differ in sampling fluxes v1 and v2 . In the deterministic formulation (panels a and b), sampling v1 determines v2 ’s values, producing probability distribution functions (pdf) that reflect their mass points. The stochastic approach (panels c and d) introduces measurement uncertainty in v1 , leading to different distributions (panel c). Adding stoichiometric constraint relaxations (panel d) provides a more realistic stochastic metabolism model, avoiding strong steady-state assumptions. Considering measurement uncertainty and stoichiometric relaxations results in a broader pdf and a feasible region (polytope) that is not necessarily convex. 14 Figure 1.4: Solution space and pdfs associated with sampling. (a): deterministic formulation without sampling. (b): deterministic approach with sampling. (c): stochastic approach with uncertainty in v1 . (d): stochastic approach with additional stoichiometric relaxation [30]. 15 The quality of sampling affects conclusions, making the choice of flux sampling algorithm crucial, whether deterministic or stochastic (e.g., Gibbs sampler [ 31 ], not covered here). Figure 1.5 shows two examples of sampling the same feasible region. In the left panel, the sampling is ineffective, limited to one area, while the right panel shows well-representative sampling of the entire model. Figure 1.5: Difference between ineffective sampling of the polytope (left panel) and well-representative sampling (right panel). Sampling algorithms are used in flux sampling to uniformly sample the polytope defined by model constraints, creating a sequence of points called a chain. This thesis used the Coordinate Hit-and-Run with Rounding (CHRR) to sample the entire feasible region and the Corner-Base Sampling (CBS) to sample the vertices of the polytope. CHRR is based on Markov chains, necessitating an introduction to their fundamental properties. 1.4.1 Markov Chains A Markov chain X is a finite-time stochastic process {X0, X1,···} characterized by the Markov property, or memorylessness: P(Xt∈A|X0, X1,··· , Xt−1) = P(Xt∈A|Xt−1),(1.7) where A is any set containing states. This equation defines the dependence of state Xt only on Xt−1 , not on past states. In this thesis, all possible states of the chain belong to Rn , where n is the number of fluxes in the model. Thus, it is a general, non-discrete state space. All presented concepts will pertain to the discrete formulation, with extensions to general spaces requiring only technical adjustments [32]. Given the current state of the chain, there is a probability distribution for the next state. If there are n possible states, the transition matrix P, of dimensions n·n , fully describes the chain by indexing various states in rows 16 and columns. The value P( x, y ) represents the probability of transitioning to y at time t + 1 from x at time t . The matrix Pis stochastic because each row contains positive values summing to one, representing a probability distribution in Rn . The most important concept, especially for the CHRR algorithm, is the stationary distribution π . A distribution π is stationary for a Markov chain Pif π P=P. Specifically, if the stationary distribution of Xt is π , it will be the same for Xt+1 , remaining unchanged over time. A Markov chain converges to the stationary distribution if it is: • Irreducible: for all states x and y , there exists t≥ 0 such that Pt(x, y)>0. • Aperiodic: the chain does not periodically oscillate between state groups. Let τ(x) = t≥1 : Pt(x, x)>0 be the set of times for the chain to return to x . The period of state x is the greatest common divisor of τ(x) . The chain is aperiodic if for every state x , gcd τ(x) = 1. Irreducibility also ensures the uniqueness of the stationary distribution π , which is representative of all reachable states, i.e., the entire polytope in our case, because it remains unchanged over time regardless of new points added. 1.4.2 Internal Sampling of the Polytope Hit and Run (HR) The first approach for sampling the feasible region is based on the Hit and Run (HR) algorithm. It follows the Markov Chain Monte Carlo (MCMC) methodology, characterized by Markov chains where points are generated through a random process. The Markov property is respected, meaning the next sampled point depends only on the last point added to the chain, the current state, and not on any previous states. HR collects samples from an N -dimensional convex space called P , the polytope generated by the stoichiometric matrix, arbitrarily identifying the point v0∈P as the first state of the chain (see Figure 1.6). Additionally, the variable i is initialized as an iteration counter, defined as follows: 1. Choose an arbitrary direction Θ i from a uniform distribution limited by the unit sphere in RN. 2. Find the minimum λmin and the maximum λmax such that vi + λmin Θ i∈ Pand vi+λmaxΘi∈P. 3. Randomly choose the step size λiwithin the interval [λmin, λmax]. 4. Generate the new sample vi+1 =vi+ Θiλisuch that vi+1 ∈P. 17 HR can be efficient even in high-dimensional polytopes, but despite the convergence to the target distribution being ensured by the Markov property, irreducibility, and aperiodicity of the chains, it is generally not used due to the phenomenon of slow mixing. This occurs when, during exploration, one approaches corners or narrow regions of the polytope. Consequently, the calculated step sizes are short, leading to the generation of very similar samples, risking entrapment in the sector for many iterations. There is a parameter, called the seed, used by the algorithm to generate all the random values used in the various iterations. Indeed, the same seed value results in identical chains. This parameter is commonly included in all algorithms with a stochastic component. Figure 1.6: An example of polytope Psampling using the HR algorithm. Coordinate Hit-and-Run with Rounding (CHRR) ACHR [ 33 ] was developed to address the slow mixing issue of the HR algorithm. The goal is to identify optimal directions to utilize larger step sizes. Setting M as the number of points called warm-up points, M≥N iterations are performed as defined by HR. During the first iteration, a point v0∈P is determined as the first element of the chain, and the artificial center ˆc = v0 is initialized. After the warm-up phase, the actual sampling begins: 1. Choose a number a uniformly distributed in the interval of integers [0, i] and calculate the direction Θi=va−ˆc ∥vi−ˆc∥. 2. Determine the step size λias in HR. 3. Generate the new sample vi+1 =vi+ Θiλi. 4. Update the estimated center ˆc=iˆc+vi i+1 . 18 Samples generated by ACHR depend on all the points already present in the chain, which makes the algorithm non-Markovian. Without the Markov property, convergence to the stationary distribution π is not guaranteed. Additionally, ACHR exhibits slow sampling in high-dimensional polytopes. There is an additional parameter K , called the thinning value, which reduces the autocorrelation of samples by adding points to the chain only every K iterations. As previously described, the shape of the polytope P can negatively affect the HR process (ill-conditioning). To address this issue, CHRR [ 34 ] employs a two-phase strategy: rounding and sampling. The rounding phase involves two operations: 1. Construct an ellipsoid of maximum volume inscribed within the polytope P. 2. Round Pby transforming the inscribed ellipsoid into a unit sphere. The sampling phase is executed using the coordinate hit-and-run (CHR) algorithm [ 35 ], a variant of HR capable of sampling the rounded polytope. At each CHR iteration, the direction Θ i is randomly determined within the space, respecting the Markov property. Finally, it is necessary to transform the coordinates of the sampled points in the rounded polytope back to their original positions in P . Unlike ACHR, CHRR guarantees convergence to the target distribution because it is based on the CHR algorithm, which belongs to the MCMC methods family. The thinning parameter is also present in this case. 1.4.3 Corner-Based Sampling (CBS) The algorithms presented so far sample the entire feasible region, but there is a second approach in the realm of flux sampling based on optimizing fluxes through random objective functions [ 36 ]. Each sample is the result of an optimization process, thus representing a vertex of the feasible region rather than a point within it. In this thesis, the CBS algorithm was implemented as follows [1]: 1. Randomly choose which reactions take part in the objective function Z. 2. Random weights wiare uniformly drawn from the interval [−1,1]. 3. The weights are then divided by the maximal flux obtained using Flux Variability Analysis (FVA). 4. Finally, each objective function is either maximized or minimized with equal probability 0.5. 19 At the end of the initialization procedure, all the previously calculated parameters are used to perform S iterations of the following steps, with S being the desired number of samples. In models of relatively low dimensions, such as ENGRO2, CBS requires only a few samples to identify all the vertices of the polytope, leading to chains of flow solutions that are very similar to each other or duplicates, in the case of sufficient rounding. For flux sampling algorithms based on MCMC methods, such as CHRR it can be useful to sample the model multiple times, with different seeds, increasing the quantity of collected points each time. This allows for analyzing the different behaviors of the algorithm in relation to the number of sampled points. However, concerning CBS, the absence of the concept of a chain, due to the purely stochastic process of the algorithm, renders the increase in the cardinality of sampling meaningless. For instance, performing 21 samplings of 1000 points each, for a total of 21000 samples, is not different from performing samplings from 1,000 to 6,000 points with a step of 1000, again totaling 21000 samples. There is no difference because CBS always behaves the same regardless of the collected samples and their number. This is not the case for the CHRR because the chains depend on at least one past state, so the sampling cardinality can influence the behavior of the algorithm. 1.4.4 False discoveries and sampling best practises HR sampling algorithms can encounter convergence issues, such as insufficient sampling to fully characterize the solution space or exploration limited to a subset of the entire solution space. Consequently, when conducting differential metabolic flux analysis, comparing flux distributions from two patient-specific or tissue-specific models may lead to false discoveries if this bias is not properly addressed [ 17 ]. For instance, rejecting the null hypothesis that two reactions have identical metabolic flux (type I error) in the two models could occur due to sampling bias, as shown by Figure 1.7. Recent findings suggested that CHRR is the most promising Monte Carlo sampling algorithm. Additionally, among the potential configurations of MC sampling parameters, employing a high thinning value is crucial for reducing the FDR [1]. 1.4.5 Samples statistics The post-processing of flux distributions obtained through flux sampling can represent a computational challenge because of the high dimensionality of the feature space and the large sample size required to minimize the risk of false discoveries. For these reasons, it is possible to apply statistical analyses to summarize the flux distributions. For a generic metabolic model, the resulting matrix containing samples of size n×m (where m is the number of 20 Chapter 2 Methods 2.1 Cancer Cell Line Encyclopedia The Cancer Cell Line Encyclopedia (CCLE) [ 5 ] Project is a comprehensive and collaborative effort aimed at cataloging and analyzing the genetic and pharmacological characteristics of cancer cell lines. Initiated by the Broad Institute and Novartis, the project aims to provide a resource for the scientific community to facilitate cancer research and drug development. In this thesis, bulk transcriptomics (RNA-seq [ 12 ]) and intracellular metabolomics data from 512 different cancer cell lines were downloaded from CCLE. Bulk samples differ from single-cell RNA-seq (scRNA-seq) in that gene expression is measured over a collection of cells from the same tissue (bulk) rather than from individual cancer cells. This technique overcomes many of the technical challenges associated with scRNA-seq experiments, where noise can significantly impact data quality due to the reduced sample size. All the tumor samples were cultured in the same RPMI 1640 medium [ 42 ], ensuring a consistent nutrient environment (extra-cellular metabolomics). The cancer cell lines originate from various primary sites, as illustrated in Figure 2.2, with over 25% of the tumors affecting hematopoietic and lymphoid tissues (e.g., blood, bone marrow, lymph, and lymphatic system). In terms of pathology, 68.55% of cells are from primary sites, while the remaining are metastatic. The most common histological class is carcinoma (malignant epithelial tumor), comprising more than 60% of the samples. Moreover, 56.92% of the samples are from male patients. The age of patients ranges from 1 to over 80 years, with a median age of around 55 years (see 2.1). 27 Figure 2.1: Age and gender boxplots of the analysed cell lines. 28 Figure 2.2: Explorative analysis of the analysed cancer bulk samples. A) Primary sites relative frequencies. B) Histologies relative frequencies. 29 2.2 ENGRO 2 The human metabolic network ENGRO2 was used to integrate the cell transcriptomics to better characterize metabolic fluxes. It is the evolution of the ENGRO1 model [ 43 ] and it focuses on central carbon metabolism and essential amino acid metabolism. This network was recently used to create the integration framework, INTEGRATE [6]. ENGRO2 is characterized by 462 reactions (394 irreversible and 67 reversible) 391 metabolites, and 495 genes. It is formed by three metabolic compartments which contain a unique set of enzymes and conditions optimized for particular biochemical reactions: extracellular, cytosol and mitochondrion. The medium of a model refers to the set of external conditions and nutrients available to the system being modelled. It specifies which metabolites can be transported into and out of the system, essentially defining the environment in which the metabolic network operates. This is crucial for simulating realistic biological conditions and understanding how changes in the external environment affect cellular metabolism. ENGRO2 medium contains 36 elements: •20 amino acids plus cystine •Glucose and galactose •Cholesterol and palmitic acid •Pyruvate, lactate •Putrescine •Folic acid (only vitamin present) •Thymidine •Oxygen (essential) •Phosphorus •Hydrogen, oxygen and water •Tetrahydrobiopterin (essential for synthesising tyrosine from phen) Of these, the following are essential for cell growth: •the nine essential amino acids plus arginine •the oxygen. •the phosphate The selected 512 cancer bulk sampeles were cultivated with the same medium called RPMI 1640 [ 42 ]. Over a total of 462 reaction, ENGRO2 has the following GPRs: 30 •196 reactions with a single associated gene •16 reactions with AND rules •96 reactions with OR rules •35 reactions with complex rules •114 reactions without associated genes 2.3 Data integration 2.3.1 Cell-specific metabolic models INTEGRATE (Model-based multi-omics data INTEGRAtion to characterize mulTi-level mEtabolic regulation) is a computational framework intended to precisely characterize the landscape of metabolic regulation in various biological samples, starting from metabolomics and transcriptomics data. It integrates transcriptomics data into a constraint-based steady-state model to predict fluxomics and allows for correlation analysis between transcriptomics and metabolomics. It applies three types of constraints: 1. Medium/nutrient availability constraints: The network can only exploit the supplied metabolites present in the growth medium. This can influence the model uptake upper-bounds. 2. Extracellular fluxes: The ratio between consumed and produced metabolites of each single reaction can be estimated and constrained. 3. RNA-seq transcriptomics constraints: For each reaction, the Reaction Activity Score is computed based on RNA-seq and GPRs. In this thesis, the INTEGRATE pipeline was applied to generate new cellspecific metabolic models by integrating transcriptomics data through the computation of the Reaction Activity Score (RAS).This approach permitted the execution of flux optimization techniques, such as Flux Balance Analysis (FBA) and Flux Variability Analysis (FVA), along with flux sampling algorithms, including Coordinate Hit-and-Run with Rounding (CHRR) and Constraint-Based Sampling (CBS), across different cell profiles. The development environment used to create the cell-specific models and predict fluxes derived from the INTEGRATE framework was based on the COnstraintBased Reconstruction and Analysis Toolbox (COBRA toolbox) [ 18 ] and its Python implementation called COBRApy [ 19 ]. Only INTEGRATE constraints based on RNA-seq and medium were used because exchange fluxes ratio were not available. Moreover, only the data integration step of the INTEGRATE pipeline was applied to create cell-specif models characterized by gene expression level (the downstream analysis to identify the level of 31 regulation of fluxes using intracellular metabolomics is beyond the scope of this work). All the cells considered in the analysis shared the very same growth medium, hence nutrient availability constraints are shared amongst all models. Consequentially, type 1 constraints set the theoretically maximum fluxes upper-bound, while RNA-seq data constraints define the effective bounds for a given reaction and cell. RNA-seq data integration defines a new cell-specific metabolic model tailored specifically for each individual cell (see Figure 2.3). The Reaction Activity Score (RAS) is employed to integrate RNA-seq datasets with GPR rules. The RAS assumes additive contributions from enzyme isoforms to reaction activity, with the least expressed enzyme subunit being the limiting factor. For each cell line c in the set C of cell lines, for each sample ξ , and for each reaction r , the RAS RASc,r,ξ is calculated by resolving the corresponding logical expression. Specifically, for genes connected by an AND operator, the minimum transcript level is taken: RASc,r,ξ = min{Tg,c :g∈Gr}(2.1) where Gr denotes the set of genes encoding the subunits of the enzyme catalyzing reaction r . For genes connected by an OR operator, the sum of their values is taken: RASc,r,ξ =XTg,c :g∈Gr(2.2) In scenarios involving both AND and OR operators, their standard precedence is followed. The reaction activity score for each cell line is then averaged over the samples: RASc,r =⟨RASc,r,ξ⟩ξ(2.3) Subsequently, the scores are normalized against the maximum RAS of all cell lines: RASc,r =RASc,r max{RASc,r}c (2.4) If the RAS is zero for all cell lines, the normalized score is also set to zero. For reactions not associated with a GPR, the RAS is set to 1. The RAS and type 1 constraints of the 513 bulk samples were integrated, as done by [ 6 ], into the constraint-based steady-state ENGRO2 model. Given a stoichiometric matrix S of dimensions M×N , where M represents metabolites and N represents reactions, along with a flux vector v , the steady state 32 assumption dictates that S·v = 0. This implies that the null space of the matrix S corresponds to the set of flux vectors v that satisfy the steady state condition. The space of feasible fluxes can be constrained using convex half-planes. These half-planes are defined by the vectors vL and vU , which represent the lower and upper bounds, respectively, for the components of the flux vector v . Thus, the bounds restrict the flux values to ensure they remain within biologically plausible limits. In mathematical terms: - The steady state condition is given by: S·v= 0 - The flux vector vis bounded by: vL≤v≤vU where vL and vU are vectors specifying the lower and upper bounds for each component of v . This approach ensures that the flux vectors v lie within a feasible region defined by the intersection of the null space of S and the convex polytope defined by the bounds. Type 1 constraint were integrated by constraining exchange fluxes. More in depth, for each medium metabolite, an exchange flux was included in the model having upper-bound set to zero and lower bound proportionally tuned to the uptaken metabolite concentration. RAS constraints were integrated by assuming that a reaction having lower RAS has lower capability to run its flux in within the network. For this reason, INTEGRATE makes the boundaries vL and vU of a given cell c dependent on the RASc,r . The complete INTEGRATE procedure to limit the flux of an internal reaction r was modified as follows (type 2 constraints were not available): 1. FVA was performed to compute the effective minimum and maximum bounds of r across all the cell-lines. FVA was run over the generic ENGRO2 model, constrained by type 1 constraints only which are shared by all the bulk samples. 2. Once the vectors vc MIN,r and vc MAX,r were obtained by FVA, each cellline model was differentiated by integrating RAS-dependent constraints: The reaction r in cell c can be mathematically described using the following constraints: RASc,r ·vc MIN,r ≤vc r≤RASc,r ·vc MAX,r (2.5) 3. Otherwise, if ris not associated with a GPR: vc MIN,r ≤vc r≤vc MAX,r (2.6) 33 Figure 2.3: An example of two cell-specific constraint-based polytopes generated from the same generic metabolic model by INTEGRATE. For example, given that the flux v2 is up-regulated in model A and down-regulated in model B , the upper bound of flux v2 is reduced in model B . The two obtained subregions can be internally sampled. Computationally-derived fluxomics All the different cell-specific models were sampled using the CBS and CHRR algorithms to obtain the predicted fluxes. The CBS algorithm generated 10,000 samples per cell, sharing the same random coefficients among the sampled models to avoid any bias related to coefficient randomness. Best practices for internal sampling of the polytopes were applied to effectively gather a representative set of fluxes using the CHRR algorithm. Specifically, 10 different batches of 1.000 samples each were run, with thinning set to 100. For each cell, a total of 10.000 samples were collected. Given the high dimensionality of the sampled solutions (462 reactions) and the large number of sampled fluxes necessary to avoid false discoveries or weak distributions, various statistics were computed to capture the main characteristics of the sampled space: mean, median, mode and fluxes, and Spearman correlation with respect to the Biomass reaction. INTEGRATE, unlike scFEA, has the capacity to create cell-specific models that can be processed using various techniques such as FBA, FVA, pFBA and sensitivity analysis. In this thesis, all these techniques were utilized. Specifically, FBA, pFBA and sensitivity analysis were performed by selecting the Biomass reaction as the objective function. 34 2.3.2 scFEA Single-cell flux estimation analysis (scFEA) is a novel computational method intended to estimate the relative rate of metabolic flux at single-cell resolution from single cell RNA-seq data [ 7 ]. Steady-state constraints, used by frameworks as INTEGRATE, may not be appropriate for diseases with significant metabolite imbalances, such as cancer. Moreover, constraint-based steady-state integration methods do not directly use gene expression to model cellular fluxes, but they are imposed as new constraints on the model to guide the search for solutions. The generation of the scFEA fluxes has to be addressed to a collaboration with Dot. Fabio Marini ([email protected]). scFEA distinguishes itself from all these methods in the following ways: •Linearity between gene expression and reaction rates is not assumed. • Gene expression data is directly used to model cell fluxes without the need for an intermediary constraint-based model, which typically requires optimization techniques like FBA or flux sampling. • It is not limited to relatively small portions of metabolic maps, such as INTEGRATE with ENGRO2, but can leverage complete metabolic maps. In this thesis, scFEA was performed and validated over bulk samples data, thus bulk RNA-seq data, as done with the INTEGRATE pipeline. It is relevant to notice that INTEGRATE generates cell-specific metabolic models that can be successively sampled and analyzed, while scFEA generates a compact flux distribution for each input cell. This method effectively addresses four computational challenges: 1. The metabolic conditions of cells are influenced by several fundamental factors. Chief among these is the availability of external nutrients, which plays a crucial role. This availability can vary widely and result in significant differences in cell phenotypes and metabolic states. 2. The entire human metabolic network is highly complex. Therefore, an appropriate reduction and reconstruction of the network are necessary to balance the characterization of the metabolic state with computational feasibility. 3. The intricate nonlinear relationship between transcriptomic expressions and metabolic reaction rates requires a more sophisticated model to fully leverage these connections. 4. Alternative enzymes with similar functions can lead to common metabolic phenotypes. However, which enzymes specifically share this common effect on altering metabolic flux remains largely unknown. 35 scFEA implements the following solutions, shown by Figure 2.4: 1. An optimization function derived from a probabilistic model to consider the budget constraints of the flow among numerous cells with different metabolic fluxes. 2. An approach to metabolic map reduction based on network topology and gene expression state. 3. A multilayer neural network model to capture the nonlinear dependency of metabolic flux on the expression of enzymatic genes. 4. A new graph neural network architecture and a solution to maximize the overall flow of intermediate substrates in all cells. scFEA is based on Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic map [ 44 ] which was enriched and optimized and finally to have including 862 genes of 390 enzymes, 1880 reactions and 1219 metabolites. The elements of the metabolic map can be depicted as a factor graph, which is a structure that aids in analyzing the relationships between variables. In this context, a factor graph is a bipartite graph where nodes represent variables and factors. The variables correspond to reactions, while the factors are metabolites. These elements are linked together based on the stoichiometric matrix, a representation that outlines the reactions within a cell, detailing the amounts of each metabolite consumed or produced in a reaction. In the factor graph, an edge connects a reaction to a metabolite if that metabolite is either a substrate or a product of the reaction. Consequently, metabolites act as connecting nodes between different reactions, forming a network of interactions that illustrates the flow and clarifies the relationships between reactions and metabolites. The complexity of the factor graph was reduced by leveraging its topological structure, grouping connected reactions into modules, also known as metabolic modules. This was done by ensuring that no intermediate metabolites had more than one incoming or outgoing reaction corresponding to multiple modules and that no intermediate metabolites had input or output fluxes different from the combined reactions or the module’s inputs and outputs. Under these conditions and the flow balance condition, it was shown that changes in reactions within the module do not affect external reactions. Thus, solving the flow of each individual reaction within a unified module is equivalent to solving the module’s flow. The combined reactions form a variable node containing multiple reactions in the factor graph. Specific classes were identified for metabolites, such as different types of fatty acids, pyrimidines, purines, and steroid hormones, forming highly connected metabolic pathways. Instead of solving the flow for each individual metabolite, metabolites of the same class were treated as a single factor. A truncated Gaussian model was used to determine which genes were 36 • Compute Dpca to retain only a predefined number of principal components ( PCs . This reduced matrix is used to compute the affinity matrix, while the imputation step works the full data matrix. 2. Distance matrix • Compute the cell-cell distance matrix Dist between cells over Dpca . The Euclidean distance between two cells p = ( p1, p2, . . . , pn ) and q= (q1, q2, . . . , qn) in an n-dimensional space is given by: d(p, q) = v u u t n X i=1 (pi−qi)2(2.10) 3. Affinity matrix • Define the affinity matrix A using the Gaussian kernel on the distance matrix Dist: Aij = exp −∥xi−xj∥2 2σ2(2.11) where xi and xj are the expression profiles of cells i and j , and σ is a scaling parameter. The Gaussian kernel is used to convert distances into affinities in such a way that similarity between cells decreases exponentially with their initial distances. A fixed kernel would alter and deform the shape of the data as cells in denser regions have more neighbors, while cells in sparse areas have fewer cells to exchange information. The same data with the adaptive kernel as in MAGIC reveals finer structures in the data. Instead of fixing a single value for the kernel width σ , this value is adapted for each cell, based on its local density. In this way, the kernel is wider in sparse areas and smaller in dense areas. Specifically: σi=distance(i, neighbor(i, ka)) (2.12) where σi is the kernel standard deviation for cell i and neighbor ( i, ka ) represents the distance of cell i from its ka - neighbor. Each cell is allowed to have at most k neighbors where k= 3ka. 4. Markov Transition Matrix • Normalize the affinity matrix A to obtain the Markov transition matrix M , where Mij represents the probability to transit from cell ito jin one step: Mij =Aij PkAik (2.13) 43 5. Diffusion Process • Apply the diffusion operator to propagate information across the graph. This is achieved by raising the Markov transition matrix to a power t , where t is the diffusion time parameter. Raising M to the power t results in a matrix where each entry represents the probability that a random walk of length t starting at cell i will reach cell j . Thus, when a Markov matrix is raised to a power, it reduces the importance of noise dimensions with near-zero explanatory power. In this process, the signal is filtered from the noise. Therefore, as t increases, similarity based on high-frequency trends (which often correspond to technological noise) decreases, and the affinity matrix represents similarity along lower-frequency trends that follow data density. As a result, after raising M to a power, phenotypically similar cells should have a strong weighted entry, while spurious neighbors are down-weighted. • The true manifold structure of the data is represented by the top eigenvectors of M, while the remaining eigenvectors are likely indicative of noise. The eigenvalues, which fall within the range [1,0], are gradually diminished through exponentiation. • The potential diffusion times belong to two regimes: an imputation regime and a smoothing regime. The initial steps of diffusion, referred to as the imputation regime, reduce the noise dimensions, bringing these small eigenvalues to zero and eliminating most of the noise in the data. As t increases, cells learn missing values from their neighbors, and it captures the correct relationships between cells that are biologically very similar but were only separated by collection artifacts. Thus, in the imputation regime, the imputed matrix changes rapidly from one iteration to the next. • In the smoothing regime, t is sufficiently large to have reconstructed the manifold with the majority of the noise eliminated. When diffusion has created a shared support for cells, additional diffusion will smooth out the lower frequency trends in the data, which probably represent genuine biological signals. Thus, optimal tuning of t requires determining the point at which noise removal begins to turn into signal removal. • MAGIC algorithm optimizes t by computing the R2 coefficient between the M matrix at t− 1 and t . If the change in R2 is less than 0.05 (the change in correlation between the two matrices is less than 5%), then t−1 has to be used: R2(Mt, Mt−1) = 1−SSE(Mt, Mt−1) 1−SST(Mt, Mt−1)(2.14) 44 where SSE is the sum of squared error and SST is the sum of squared total. 6. Imputation • Impute the gene expression values by diffusing the original logscaled expression matrix Dusing the diffusion operator: Dimputed =Mt∗D(2.15) where Dimputed is the imputed gene expression matrix.In the imputation step, MAGIC enhances the data by leveraging information from cells within each neighborhood. This is achieved by multiplying the transition matrix with the original data matrix, effectively aligning cells with their underlying manifold. During this data diffusion process, cells exchange information with their local neighbors, analogous to the diffusion of heat through data. Raising the diffusion operator to the t -th power corresponds to a t -step random walk through the data. This exponentiation acts as a low-pass filter on the eigenvalues, removing noise dimensions with small eigenvalues while learning the manifold structure. Although PCA is used to enhance robustness when computing the affinity matrix, the imputation is performed on the count matrix before applying PCA. As a result, while data is averaged across cells, each individual cell retains a unique neighborhood, leading to a distinct expression vector for each cell. MAGIC modifies the data by imputing missing values and smoothing the expression of genes across similar cells. Specifically, it: • Denoising: Reduces the technical noise by averaging the expression values over similar cells, thereby enhancing the signal-to-noise ratio. • Imputation: Fills in the dropout events, which are common in singlecell RNA-seq data, where a gene may appear to be unexpressed in a cell due to technical limitations rather than biological absence. • Smoothing: Ensures that the gene expression levels are more consistent among similar cells, thereby recovering continuous gene expression patterns that may be disrupted by technical noise. • Preservation of Biological Variability: Maintains the underlying biological variability while reducing the random noise, allowing for more accurate downstream analyses such as clustering, trajectory inference, and differential expression. 45 It is relevant to discuss MAGIC robustness relating to its parameters and the capacity to avoid introducing numerical artifacts that could bias the biological clusters. MAGIC requires three inputs: ka (to set the adaptive kernel to the distance of the k -th nearest neighbor), t (the number of times the matrix M is powered), and npca (the number of PCA components used to construct the affinity matrix). MAGIC was tested by [ 8 ] to identify optimal parameter ranges to avoid over-smoothing of the dataset. In particular: • MAGIC is robust with ka ≤ 30, while higher values could lead to oversmoothing. MAGIC sets the maximum number of K-nearest neighbors to 3 ·ka. • MAGIC is robust with t≤ 24, while higher values could also lead to over-smoothing. • MAGIC is highly robust with respect to the number of PCs used. Generally, using too few PCs may fail to capture sufficient variance of the dataset. While MAGIC recovers structure by diffusing values between neighboring cells, it is crucial that values do not transfer between different clusters. In fact, [ 8 ] demonstrated that the cluster structure remains intact even after applying MAGIC. It is important to highlight that MAGIC was originally intended for singlecell denoising, which is characterized by a very high level of noise. This is due to the fact that scRNA sequencing technologies suffer from many sources of technical noise, including dropout, which can obscure important gene-gene relationships [ 8 ]. In contrast, this analysis focused on bulk RNA-seq data, where the RNA-seq data derives from multiple cells of the same tissue, thus reducing the technical difficulties seen in single-cell sequencing. As previously mentioned, too high values of ka and t could lead to over-smoothing and thus compromise the original cluster structure. To avoid over-smoothing, a conservative approach to the algorithm’s parameters was adopted, as shown in Table 2.1. These parameter values are recommended by [ 46 ] and are at the lower end of the ranges suggested by the original MAGIC study. Table 2.1: MAGIC parameters Parameter Symbol Value K-nn for σika 5 Maximum K-nn k3·5 = 15 Diffusion operator power t3 Principal components PCs 100 46 2.4.2 Scaling In preprocessing data for machine learning models, it is often essential to scale the features to ensure they are on a similar scale. Two common scaling techniques provided by [ 22 ] are Min-Max Scaling and Standard Scaling. Each omic dataset was pre-processed using both methods before applying PCA, and subsequently, clustering, to ensure comprehensive coverage. Min-Max Scaling Min-Max Scaling, also known as normalization, scales the data to a fixed range, typically [0, 1]. This scaling method transforms each feature individually such that it is bounded by a specified minimum and maximum value. The formula for Min-Max Scaling is: X′=X−Xmin Xmax −Xmin (2.16) where X is the original feature vector, Xmin and Xmax are the minimum and maximum values of the feature, respectively, and X′ is the scaled feature. This transformation ensures that all features are in the same range, which can be particularly useful for algorithms that rely on distance calculations such as k-means clustering and k-nearest neighbors. It is useful when the data does not have outliers and when the features have a bounded range. It is often preferred for algorithms that do not assume normally distributed data. Standard Scaler Standard Scaling, also known as z-score normalization, transforms the data to have a mean of zero and a standard deviation of one. This method standardizes each feature individually by removing the mean and scaling to unit variance. The formula for Standard Scaling is: X′=X−µ σ(2.17) where X is the original feature vector, µ is the mean of the feature, σ is the standard deviation of the feature, and X′ is the standardized feature. This transformation ensures that each feature contributes equally to the analysis, which can be particularly useful for algorithms that assume normally distributed data, such as linear regression and principal component analysis (PCA).It is useful when the data has outliers or when the features have different units or scales. It is often preferred for algorithms that assume normally distributed data. 47 2.4.3 Principal Component Analysis (PCA) Principal Component Analysis (PCA) is a widely-used dimensionality reduction technique in machine learning and data analysis. It transforms the data into a new coordinate system such that the greatest variances by any projection of the data come to lie on the first principal component, the second greatest variances on the second principal component, and so on. This technique is particularly useful for reducing the dimensionality of large datasets, improving computational efficiency, and removing noise. The cumulative explained variance represents the percentage of the original total variance captured by the first n principal components. Each omic dataset was subjected to dimensionality reduction to facilitate clustering analysis. For each omic dataset, n different components were selected such that their cumulative explained variance was at least 80% (see Table A.1 in appendix). The trasncriptomics dataset, given their extremely high dimensionality, did not reach this threshold and 70 principal components were chosen based on the elbow method. The formula for cumulative explained variance is: Cumulative Explained Variance = n X i=1 λi Pp j=1 λj!×100% (2.18) where λi is the eigenvalue of the i -th principal component, and p is the total number of original variables. The goal of PCA is to identify the principal components that capture the maximum variance in the data. The steps involved in PCA are as follows: 1. Standardize the Data: • Given a data matrix X ∈Rn×p , where n is the number of samples and p is the number of features, standardize the data to have a mean of zero and a standard deviation of one: Xstd =X−µ σ(2.19) where µ is the mean and σ is the standard deviation of each feature. 2. Compute the Covariance Matrix: •Calculate the covariance matrix Cof the standardized data: C=1 n−1X⊤ stdXstd (2.20) 48 3. Eigenvalue Decomposition: • Perform eigenvalue decomposition on the covariance matrix to obtain the eigenvalues and eigenvectors: CV =VΛ (2.21) where Vis the matrix of eigenvectors and Λis the diagonal matrix of eigenvalues. 4. Select Principal Components: • Sort the eigenvalues in descending order and select the top k eigenvectors corresponding to the largest eigenvalues. These eigenvectors form the principal components. 5. Transform the Data: • Project the original standardized data onto the selected principal components to obtain the reduced-dimensionality data: XPCA =XstdVk(2.22) where Vkcontains the top keigenvectors. PCA offers several benefits and applications, including: When using PCA, it is important to consider the following points: • Linear Assumption: PCA assumes that the principal components are linear combinations of the original features. It may not capture complex, non-linear relationships in the data. • Variance Retention: The number of principal components to retain is a critical decision. It is often based on the proportion of variance that needs to be preserved, which can be determined using a scree plot or by setting a cumulative explained variance threshold. • Interpretability: The principal components are linear combinations of the original features and may not always be easily interpretable in the context of the original features. 49 2.5 Clustering Analysis The aim of this thesis is to characterize the metabolic heterogeneity of cancer cells by utilizing various biological data, such as RNA-seq and fluxomics, referring to the same set of cells. To achieve this, clustering analysis was employed to uncover hidden associations and patterns within each of the analyzed omics datasets. Subsequently, the subpopulations identified in each omic were compared. Each data source underwent a clustering grid search for hyperparameter optimization based solely on internal evaluation metrics, such as Silhouette and Davies-Bouldin scores. Following this, biomarkers were identified, along with a graphical association analysis between clusters and cell metadata, including cancer histology, pathology, and primary site. Biomarkers were not extensively interpreted, as it is beyond the original scope and would require expertise from life scientists. 2.5.1 Hyperparameters tuning The hyperparameters optimization was not guided by external evaluation metrics, such as v-measure, to maximize the concordance with a ground truth, but was performed independently. In order to run an exhaustive unsupervised analysis, four main clustering algorithm families were tested: • Centroid-based clustering algorithms, such as K-means [ 20 ] and Kmedoids [ 21 ], are most effective when data is likely to form spherical clusters of similar sizes. These algorithms work well with large datasets and are efficient in terms of computational complexity. K-means is particularly useful when the goal is to minimize variance within clusters, while K-medoids is more robust to noise and outliers. • Hierarchical clustering algorithms, like agglomerative clustering [ 22 ], are ideal when the data has a hierarchical structure. These algorithms do not require the number of clusters to be specified in advance and are beneficial when the dataset contains nested clusters. They are also useful for smaller datasets or when the computational cost is not a primary concern. • Density-based clustering algorithms, such as HDBSCAN [ 47 ], are suitable for identifying clusters of arbitrary shapes and handling datasets with noise. They are particularly effective when the clusters are densely packed and separated by low-density regions. HDBSCAN is beneficial when dealing with spatial data or datasets where clusters are expected to have varying sizes and shapes. • Graph-based clustering algorithms, like the Leiden algorithm [ 23 ], are optimal for large-scale networks or graph-structured data. These 50 algorithms are designed to optimize modularity or other community detection measures and are effective in uncovering complex relationships and communities within the data. The Leiden algorithm, in particular, is well-suited for large datasets as it improves the quality and speed of clustering compared to similar algorithms like Louvain [48]. All the algorithm implementations are sourced from scikit-learn [ 22 ], except for the Leiden algorithm, which is implemented by the scanpy library [46]. Clustering hyperparameter optimization is crucial as it significantly enhances the performance of clustering algorithms. By fine-tuning hyperparameters, one can achieve more meaningful and well-separated clusters, leading to better data segmentation and more reliable insights. This process helps minimize errors, improve computational efficiency, and ensure that the clustering model generalizes well to unseen data. To this end, an extensive hyperparameter grid search was conducted to obtain a well-informed set of clustering settings for each pair of omic dataset and clustering algorithm. The tested parameters and values are reported in Table 2.2. Table 2.2: Clustering hyperparameters grid-search. K-means / K-medoids Hyperparameter Values K From 2 to 11 Agglomerative Hyperparameter Values K From 2 to 11 Distance metric Euclidean, Cosine, Manhattan HDBSCAN Hyperparameter Values Minimum cluster size 5, 10, 15, 25, 30 Distance metric Euclidean, Manhattan Leiden Hyperparameter Values Neighbors 15, 30, 50 Distance metric Euclidean, Cosine, Manhattan Resolution 0.1, 0.2, 0.5, 1, 2, 4 51 K-means K-means is a popular clustering algorithm used to partition a dataset into Kdistinct, non-overlapping subsets. The algorithm operates as follows: 1. Initialize K centroids µ1, µ2, . . . , µK randomly selected from the data points. 2. Assign each data point xi to the nearest centroid Cj based on the Euclidean distance: Cj={xi:∥xi−µj∥2≤ ∥xi−µl∥2∀l, 1≤l≤K}(2.23) 3. Update the centroid of each cluster to be the mean of the data points assigned to it: µj=1 |Cj|X xi∈Cj xi(2.24) 4. Repeat steps 2 and 3 until the centroids converge (i.e., the assignments no longer change) or a maximum number of iterations is reached. The Mini-Batch K-means variant was used to handle the most computationally intensive clustering tasks due to the large number of data points being clustered. This approach was particularly applied to the sampled fluxes, which were processed in batches of 1024 points. Mini-Batch K-means is an efficient version of the standard K-means algorithm. It uses small, random subsets of the dataset (mini-batches) to update the cluster centers, rather than using the entire dataset. This reduces the computational load and speeds up the convergence, making it suitable for large datasets. K-medoids K-medoids is a clustering algorithm similar to K-means, but instead of using the mean of the data points as the cluster center, it uses an actual data point (medoid) that minimizes the sum of dissimilarities within a cluster. The algorithm proceeds as follows: 1. Initialize K medoids m1, m2, . . . , mK randomly selected from the data points. 2. Assign each data point xi to the nearest medoid Cj based on a specified distance metric d(xi, mj): Cj={xi:d(xi, mj)≤d(xi, ml)∀l, 1≤l≤K}(2.25) 52 2.5.4 External Evaluation Metrics External evaluation metrics assess the quality of clustering by comparing the clustering results to external ground truth labels. These metrics help to evaluate the degree to which the clustering corresponds to the known class labels. Homogeneity, completeness, normalized mutual information score (NMI), and adjusted Rand index (ARI) were used to define a robust benchmarking pipeline. Figure 2.5: a) Homogeneity = 1 and completeness < 1. b) Homogeneity = 0 and completeness = 1. c) Homogeneity < 1 and completeness = 1. d) Homogeneity = 1 and completeness = 0. Homogeneity Homogeneity measures whether each cluster contains only members of a single class. It is defined as: h= 1 −H(C|K) H(C)(2.34) 59 where H ( C|K ) is the conditional entropy of the class distribution given the cluster assignment, and H ( C ) is the entropy of the class distribution. A homogeneity score of 1 indicates perfect homogeneity. Completeness Completeness measures whether all members of a given class are assigned to the same cluster. It is defined as: c= 1 −H(K|C) H(K)(2.35) where H ( K|C ) is the conditional entropy of the cluster assignment given the class distribution, and H ( K ) is the entropy of the cluster distribution. A completeness score of 1 indicates perfect completeness. Normalized mutual information (NMI) The NMI, also known as v-measure [ 52 ], is the harmonic mean of homogeneity and completeness. It is the normalized version of the Mutual Information (MI) used to scale the results between 0 (no mutual information) and 1 (perfect correlation). It is defined as: NMI = 2 ×homogeneity ×completeness homogeneity + completeness (2.36) This metric balances both the homogeneity and completeness (see Figure 2.5) scores to provide a single measure of the clustering quality. It is preferred over supervised learning evaluation metrics, such as accuracy, precision, recall, and F-measure because of the following properties: • Symmetry: NMI treats homogeneity and completeness equally, ensuring that both the compactness of clusters (homogeneity) and the coverage of classes (completeness) are taken into account. This symmetry is advantageous because it provides a balanced evaluation of clustering performance. • Independence from Cluster Number: Unlike accuracy or Fmeasure, NMI does not assume a fixed number of clusters, making it more suitable for scenarios where the number of clusters may vary or is not known a priori. This flexibility allows for more general applicability across different clustering tasks. • Penalization of Random Assignments: NMI effectively penalizes random cluster assignments, ensuring that the metric reflects meaningful clustering structures rather than arbitrary groupings. This property 60 helps in distinguishing between genuine clustering performance and coincidental alignments with the ground truth. • Handling of Class Imbalances: NMI is robust to class imbalances, which are common in real-world datasets. Traditional supervised metrics can be heavily influenced by dominant classes, whereas NMI provides a more equitable evaluation by considering the overall distribution of classes within clusters. The Adjusted Mutual Information (AMI) is an adjustment of the MI and NMI scores to account for chance [ 53 ]. It accounts for the fact that the NMI, and MI, are generally higher for two clusterings with a larger number of clusters, regardless of whether there is actually more information shared. The NMI was preferred over homogeneity and completeness for its ability to represent both clustering properties. Adjusted Rand index (ARI) The Adjusted Rand Index (ARI) [ 54 ] is a measure of similarity between two data clustering that corrects for chance. It is an adjusted version of the Rand Index (RI), which considers all pairs of samples and counts pairs that are assigned in the same or different clusters in the predicted and true clustering. The adjustment accounts for the expected similarity of all pairwise comparisons between clustering randomly. Given a set of n elements S = {O1, O2, . . . , On} and two clusterings of these elements, U = {U1, U2, . . . , Ur} and V = {V1, V2, . . . , Vs} , let a be the number of pairs of elements that are in the same set in both U and V , and b be the number of pairs of elements that are in different sets in both U and V. The Rand Index (RI) is defined as: RI =a+b n 2(2.37) where n 2is the total number of pairs of elements. The Adjusted Rand Index (ARI) adjusts this index to account for the expected similarity of all pairwise comparisons between clustering randomly. The ARI is defined as: ARI =RI −E[RI] max(RI)−E[RI](2.38) where E [ RI ] is the expected Rand Index of random clusterings and max(RI) is the maximum Rand Index. 61 2.5.5 Biomarkers Biological markers, or biomarkers, are critical indicators used in the identification and characterization of various biological states, including disease conditions. Biomarkers play a crucial role in interpreting the complex data derived from different omics analyses, facilitating the understanding of underlying biological processes. One of the key motivations for identifying biomarkers is their interpretability. Biomarkers provide a tangible link between the vast and often abstract datasets generated through techniques such as RNA-seq and metabolomics, and the concrete biological phenomena they represent. This interpretability is essential for translating complex omics data into actionable insights that can drive further biological research and potential clinical applications. In this thesis, the identification of biomarkers was performed using the Wilcoxon rank-sum test. This non-parametric test is typically employed to compare the distributions of gene expression levels between different groups. More generally, it can be applied to other omics, such as metabolomics and fluxomics, to identify the metabolites with the highest concentrations or the reactions with the highest fluxes. The test statistic is given by: W= n1 X i=1 R(Xi) (2.39) where n1 is the sample size of the first group, Xi is, for example, the expression value of the i -th gene in the first group and R ( Xi ) is the rank of Xi in the combined sample of both groups. By employing the Wilcoxon rank-sum test, it is possible to identify genes, metabolites, or reactions that exhibit statistically significant differences between groups, thereby highlighting potential biomarkers within clusters. The dimensions of each cluster/group were scaled between 0 and 1 as done by the minmax scaler to facilitate the comparison of markers expressions. 62 Chapter 3 Results and discussion 3.1 Metabolomics as clustering benchmark All frameworks used to predict fluxomics face the challenge of benchmarking due to technical difficulties in directly measuring metabolic fluxes at the genome-wide level to be used as ground truth. To overcome this problem, in this work, I used intracellular metabolomics as an evaluation benchmark to quantify the level of agreement of the computationally derived fluxomics with expected metabolic phenotypes. Intracellular metabolomics describes metabolic profiles of cancer cells by measuring the actual concentration of metabolites (reactants) within a sample. Even if a specific enzyme is overexpressed in a cell, the absence of the associated metabolite, the reactant, reduces the overall reaction fluxes. It is better suited as a benchmark for evaluating the goodness of flux predictions because it is not directly used during the computation of fluxomics datasets. On the contrary, predicted fluxes are directly dependent from the transcriptomics and RAs datasts. Specifically, INTEGRATE uses the RNA-seq to compute the RASs, which in turn are used to create the cell-specific models that define the feasible regions that I sampled to obtain all the flux-derived statistics. Instead, scFEA directly uses the reads counts to predict flux distributions. To obtain the metabolic clusters based on intracellular metabolomics to be used as a sort of ground truth, I tested different clustering techniques and hyperparameters. Table 3.1 presents the best clustering settings obtained using five different internal evaluation metrics. All metrics, except for the Davies-Bouldin index, indicated the presence of two clusters. The Silhouette index, with a medium-low value of 0.30, and the Calinski-Harabasz index both identified the same clustering setting, performed by K-means with K = 2. Figure 3.1 displays the clustering selected by the Silhouette index, shown through both PCA and UMAP representations. The assigned classes are quite unbalanced, with 64.71% of cells assigned to Cluster 0 and 35.29% 63 assigned to Cluster 1. The clustering with the highest Silhouette value was selected as the metabolic benchmark. Although the biological interpretation of the two clusters goes beyond the aim of this thesis, to further verify that the two clusters correspond to distinct metabolic phenotypes I performed a marker analysis. Figure A.1 shows a heatmap representing the biomarkers identified using the Wilcoxon rank test between the two assigned labels. The dimensions of each cluster/group were scaled between 0 and 1 using the min-max scaler to facilitate the comparison of marker expressions across different groups. The heatmap colors indicate the metabolite concentrations within each cluster, with the matrix columns representing top-ranked metabolites and the rows corresponding to individual cells. Overall, both clusters are characterized by a set of highly expressed metabolites: Methionine and Leucine for cluster 0, and Uridine monophosphate (UMP) and Adenosine monophosphate (AMP) for cluster 1. Figures A.2, A.3, A.4, illustrate the composition of the two clusters in terms of cell metadata, including respectively histology, pathology, and primary site. The majority of carcinoma samples were assigned to cluster 0, while cluster 1 is characterized by lymphoid neoplasms and hematopoietic neoplasms. These last two histologies refer to blood, or liquid, cancers that are likely to differentiate from solid tumors, as happened in [ 55 ]. Regarding pathology, no significant polarization of the two clusters was identified. Lastly, the primary site panel again shows a well-defined separation between liquid (lymphoid and hematopoietic tissues) and solid tumors. Table 3.1: Best clustering settings for metabolomics. Metric Value Clusters Scaling Algorithm Distortion 216.02 2 Minmax HDBSCAN Silhouette 0.30 2 Minmax K-means Calinski-Harabasz 214.30 2 Minmax K-means Davies-Bouldin 1.00 7 Minmax Agglomerative DBCV 0.21 2 Minmax HDBSCAN 64 Figure 3.1: Metabolomics clustering selected by Silhouette = 0.30: PCA and UMAP representation. 65 3.2 Transcriptomics and metabolomics share similar metabolic phenotypes Transcriptomics regulates metabolism and reaction fluxes through transcriptional control, which relies on the expression levels of genes encoding specific enzymes responsible for reaction catalysis. However, metabolism is influenced by other factors such as post-transcriptional, post-translational and metabolic controls. Metabolomics measures the concentration of metabolites, the reaction reactants, within the cell and even if a reaction enzyme is over-expressed, the lack of the necessary metabolites reduces the reaction flux anyway. It is interesting to analyze the overlap between the metabolic subpopulations identified by transcriptomics and metabolomics. Given the metabolomics benchmark identified in the previous section, a comparison was performed with the clusterings derived from the RNA-seq matrix. The concordance level between these two omics in capturing cells heterogeneity is crucial, as they are the only two omics derived from experimental analysis rather than in silico methods. Table 3.2 presents the best clustering results, obtained through gridsearch, of RNA-seq data restricted to metabolic genes, evaluated trough external metrics with respect to the benchmark. All metrics, except for homogeneity, indicated that the best clustering was achieved using the Leiden algorithm with two clusters and input data scaled by the Minmax scaler. The discrepancy in homogeneity is likely due to its bias towards a higher number of clusters in an attempt to identify purer groups. Generally, the NMI provides a better summary of homogeneity and completeness. The highest NMI is equal to 0.42, indicating a discrete overlap of the metabolic subpopulations identified by transcriptomics and metabolomics (see Figure 3.2). Table 3.2: Best RNA-seq clusterings with respect to metabolomics benchmark. Metric Value Clusters Scaling Algorithm Homogeneity 0.53 36 Standard Leiden Completeness 0.43 2 Minmax Leiden NMI 0.42 2 Minmax Leiden AMI 0.55 2 Minmax Leiden ARI 0.42 2 Minmax Leiden 66 Figure 3.2: Clustering of RNA-seq metabolic genes selected by highest NMI with intracellular metabolomics (NMI=0.42): a) RNA-seq PCA colored by metabolomics labels. b) RNA-seq UMAP colored by metabolomics labels. c) Metabolomics PCA colored by RNA-seq labels. d) Metabolomics UMAP colored by RNA-seq labels. 67 3.3 INTEGRATE fluxomics does not reflect metabolic profiles The clusterings of predicted fluxes generated through the INTEGRATE pipeline, as well as all derived omics data such as mean, median, FBA and sensitivity analysis, were compared with the metabolomics benchmark and RNA-seq to numerically evaluate the degree of coherence between them. This comparison is crucial to determine if this constraint-based integration framework, which involves in silico experiments, aligns with the metabolomic phenotypes observed in laboratory settings. For a given omic, the clustering with the highest NMI against the metabolomic benchmark was selected as the reference solution. This clustering, defined by a combination of clustering algorithm, scaling, and various hyperparameters, was then compared with all other omics, selected in the same manner. A brief summary of the INTEGRATE-derived omics is provided below: • RAS: The Reaction Activity Score is computed for each reaction and cell line by the INTEGRATE framework. It serves as an intermediary step for integrating RNA-seq data into a constraint-based model. • pFBA: Parsimonious Flux Balance Analysis is an extension of the standard FBA that aims to minimize the total flux through the network. The biomass reaction was selected as the objective function. • Sensitivity analysis: The knockout of a model reaction from the stoichiometric matrix S enables the analysis of the effect of such deletion by performing Flux Balance Analysis (FBA). The sensitivity coefficient is the fraction of Biomass flux preserved after a reaction deletion. • Mean, mode, and median of sampled fluxes: These three statistics were used to summarize the distributions of the fluxes sampled by CBS and CHRR. • Fluxes correlation with biomass: The Spearman correlation coefficient was computed to compare the sampled fluxes of each cell line with the sampled Biomass, for both CHRR and CBS. Figure 3.3 displays the NMI scores for the INTEGRATE-related omics derived from the RNA-seq matrix. All the clustering hyperparameters are detailed in Table A.2 in the Appendix. Surprisingly, all predicted fluxomics demonstrated a very poor ability to represent actual metabolic profiles. In fact, low NMIs were observed between fluxomics and metabolomics. For example, the mean fluxes sampled by CBS obtained an NMI of 0.10, while CHRR achieved an NMI of 0.09. As previously discussed, a moderate NMI 68 Figure 3.6: Scaled heatmap of the metabolic genes RNA-seq matrix markers processed by MAGIC (clustering selected by DBCV). Columns represent genes, rows are associated to cells and colors are proportional to genes expression level. To ensure that no numerical artifacts were introduced into the RNA-seq matrix, I compared the concordance level between the two transcriptomics matrices (before and after MAGIC processing) with metabolomics. The best RNA-seq clusterings with respect to external evaluation metrics with metabolomics are shown by Table 3.2, while Table 3.6 presents the same evaluation metrics after applying MAGIC. Again, all metrics indicated that two clusters were the best option with the K-medoids algorithm, except for homogeneity, which preferred 29 clusters with the Leiden algorithm. For both versions of the RNA-seq data, a moderate level of coherence with the metabolomics benchmark was observed: 0.42 before MAGIC and 0.41 after MAGIC. Moreover, the Silhouette score of the clustering selected by highest NMI improved from 0.08 without MAGIC to 0.38 with MAGIC because MAGIC enhances similarities between cells through the diffusion process, thereby reducing intra-cluster distances. However, no significant changes were observed in the external evaluation metrics. The consistency of NMI between RNA-seq and metabolomics, which remained largely unchanged, validates the effectiveness of the MAGIC algorithm. This indicates, as stated by [ 8 ], that MAGIC did not introduce numerical artifacts into the RNA-seq matrix, as the metabolic profiles remained consistent in both cases, maintaining the same level of agreement 75 with the benchmark. Table 3.6: Best RNA-seq clusterings processed by MAGIC with respect to metabolomics benchmark. Metric Value Clusters Scaling Algorithm Homogeneity 0.52 29 Minmax Leiden Completeness 0.41 2 Standard K-medoids NMI 0.41 2 Minmax K-medoids AMI 0.41 2 Minmax K-medoids ARI 0.54 2 Minmax K-medoids 76 3.5 MAGIC improves fluxomics prediction Figures 3.7 display the NMI scores for the INTEGRATE-related omics, where the RNA-seq matrix was processed by MAGIC before being input into the INTEGRATE pipeline. All the clustering hyperparameters are shown in Table A.3in the Appendix. Significant differences are observed in the NMIs between the predicted fluxomics and the metabolomics benchmark. For example, in the case of CBS, the mean fluxes increased from 0.10 to 0.30, the median from 0.08 to 0.39, the mode from 0.08 to 0.32 and the Biomass correlation from 0.07 to 0.39. Similarly, for CHRR, the mean went from 0.09 to 0.33, the median from 0.09 to 0.32, the mode from 0.08 to 0.29 and the Biomass correlation from 0.08 to 0.28. The NMI improved even for omics not derived from sampled fluxes, such as pFBA, from 0.05 to 0.23, and sensitivity analysis, from 0.09 to 0.26. Higher NMIs are also found between different omics, indicating a higher general overlap. For example, the NMI between the CHRR mode and CBS mean increased from 0.53 to 0.96, and between the CHRR median and CBS mean, it increased from 0.56 to 0.94. This suggests that the preprocessing of RNA-seq data with MAGIC significantly enhances the ability of predicted fluxomics to represent actual metabolic profiles and improves the agreement between different omic datasets. INTEGRATE generated constraint-based models and enabled the computation of all fluxes. The denoising activity of the MAGIC algorithm led to higher coherence scores. However, the RNA-seq and RAS matrices maintained relatively consistent NMIs against metabolomics regardless of MAGIC. As previously explained, the RNA-seq matrix exhibited high sparsity without denoising, while after MAGIC, the percentage of non-expressed genes significantly decreased due to the algorithm’s diffusion process. Moreover, MAGIC preserved the original information in the RNA-seq matrix, evidenced by a stable NMI coefficient with metabolomics. This indicates that adjusting the matrix to remove noise does not alter the genetic information it contains, nor does it affect the reaction scores in the RAS matrix. However, significant differences are observed in predicted fluxomics. This happened because new information emerges from intersecting the flux constraints, where the flux of one reaction can influence another when computing a feasible sample. If the RNA-seq matrix contains substantial noise, a low-quality feasible region is created. In summary, MAGIC enhances data coherence by reducing noise and preserving genetic information. Its effects on flux predictions highlight the importance of noise-free inputs for accurately defining the feasible metabolic space and deriving meaningful solutions. 77 Figure 3.7: NMI between omics with MAGIC denoising. Previously, it was shown that INTEGRATE generated feasible regions with a minimal percentage of shared feasible solutions, supporting the hypothesis that low-quality fluxes are due to the RNA-seq matrix. To reinforce this hypothesis, Figures 3.8 and 3.9 show 1000 random samples generated by CHRR for five different cancer cells, before and after MAGIC respectively. It is clear that in both cases, the samples are well-clustered and respect the ground truth of their respective cells. However, well-clustered samples do not necessarily imply that fluxes strictly represent real metabolic phenotypes. In fact, poor NMI scores were found without MAGIC. In summary, the responsibility for low-quality fluxes lies with the noise in the RNA-seq matrix, which is amplified when combining gene expression into the constraint-based model. 78 Figure 3.8: 1000 random samples generated by CHRR of five different cells, coloured by cell partnerhip: a) PCA. b) UMAP. 79 Figure 3.9: 1000 random samples generated by CHRR of five different cells, coloured by cell partnerhip after RNA-seq denoising: a) PCA. b) UMAP. Figures 3.10 and 3.11 show the impact of MAGIC on each omic’s NMI and silhouette scores, respectively. Although the RNA-seq and RAS did not improve their NMI scores, a significant increase in the silhouette score was observed, going from 0.09 to 0.38. This indicates that the quality of the clusters has improved in terms of interand intra-cluster distances because of the diffusion process that enhances similarities among already near cells. For all other omics, it is clear that both the NMI against metabolomics and the quality of the clusters (as indicated by higher silhouette scores) improved after applying MAGIC. This demonstrates that MAGIC not only enhances data coherence but also significantly improves clustering quality across different omics data, while maintaining stable benchmark concordance for RNA-seq and RAS. 80 Figure 3.10: Barplot of omics NMI against the metabolomics benchmark, both with and without MAGIC denoising. 81 Figure 3.11: Barplot of clusterings silhouette, selected by NMI against the metabolomics benchmark, both with and without MAGIC denoising. Among the summary statistics computed for the fluxes sampled by CBS and CHRR, the median and Biomass correlation were found to be optimal for CBS (NMI = 0.39), while the mean should be preferred for CHRR (NMI = 0.33). Beginning with CBS, as discussed by [ 1 ], the final marginal distributions generated by CBS approximate a discrete probability distribution due to sampling a few values multiple times from finite feasible region vertices. The median was identified as the optimal statistic because it is less sensitive to outliers that might adversely affect other statistics, such as the average fluxes. In contrast, CHRR samples the entire feasible region using algorithms, which is less likely to produce bimodal distributions, as indicated by [ 1 ]. The mean was found to be the best choice for CHRR due to the presence of relatively uniform marginal distributions. Both CHRR and CBS statistics showed high NMIs between their mode and median statistics, suggesting that both sampling strategies captured similar metabolic profiles. However, summary statistics alone may not be exhaustive. Therefore, I clustered real fluxes too, as discussed in the next section. 82 3.6 CHRR outperforms CBS in samples clustering It is useful to investigate the relationship between metabolomics and fluxes sampled by CBS and CHRR without relying on any statistics that might excessively simplify the flux distributions. For each cell, I randomly selected 50 samples only, and I used the mini-batch K-means to reduce the clustering computational cost (50 · 513 = 25 , 650 samples). To validate the subsampling, the Wilcoxon–Mann–Whitney test was performed to compare the subsampling to the original flux distribution of each reaction, which contained 10,000 values. This non-parametric test is not influenced by differences in sample size between the two distributions. Table 3.7 shows the average percentage of reactions where the test found a significant difference from the original distribution with α = 0 . 05. On average, CBS registered 1.62% of rejected reactions without MAGIC and 1.96% with MAGIC. CHRR registered 5.17% without MAGIC and 4.95% with MAGIC. The low values suggested a sufficient reliability of the subsamples. At this point, each cell was represented by 50 samples, but the metabolomics benchmark assigned just one label per sample. To compare the two assignments through NMI, each cell was assigned to the cluster containing the majority of its samples. For example, if celli has 26 samples in cluster 0, then the cell is assigned to cluster 0. Figures 3.12 and 3.13 show the NMI for the sampled fluxes without and with MAGIC, respectively. CHRR samples NMI with metabolomics shifted from 0.05 to 0.25 without and with MAGIC, respectively. However, CBS did not improve as well, shifting from 0.01 to just 0.05. Eventually, CHRR outperformed CBS in coherence with the benchmark, despite this phenomenon not being visible when dealing with summary statistics. Table 3.7: The Wilcoxon–Mann–Whitney test was used to compare the subsampled fluxes against the original flux distribution for each cell and sampling algorithm. This table reports the average percentage of rejected reactions for each sampling algorithm. Algorithm Magic Rejected reactions CBS No 1.62% CHRR No 5.17% CBS Yes 1.96% CHRR Yes 4.95% 83 Figure 3.12: CBS, CHRR and metabolomics NMI without MAGIC denoising. Figure 3.13: CBS, CHRR and metabolomics NMI with MAGIC denoising. 84 possible to integrate growth medium details into the loss function or to use a customized metabolic model to shape the metabolic pathways. Figure 3.19: Predicted fluxomics NMIs with metabolomics: scFEA and INTEGRATE. 91 Figure 3.20: Predicted fluxomics NMIs with metabolomics (MAGIC applied over RNA-seq): scFEA and INTEGRATE. 92 Chapter 4 Conclusions Metabolism is intricately linked to nearly all cellular processes, making it a key indicator of a cell’s or organism’s physiological state. A deep understanding of cancer cell metabolism, or pathological conditions in general, can reveal hidden features that could be exploited to combat these cells. Integrating omics data, such as transcriptomics and metabolomics, within metabolic models allows for the prediction of reaction fluxes by characterizing the metabolism. Unfortunately, evaluating the performance of integration methods is challenging because the ground truth of fluxes is typically limited to a small set of reactions due to high costs and technical difficulties. The objective of this thesis is to develop an evaluation pipeline for integration methods to address the benchmarking issue and leverage high-throughput data to characterize metabolic phenotypes for further biological research. This thesis is based on the assumption that cells with similar metabolite concentrations are likely to exhibit similar fluxes and, consequently, similar metabolic profiles. Transcriptomics data from 513 different cancer bulk samples were obtained from The Cancer Cell Line Encyclopedia (CCLE) Project and integrated using two methods: INTEGRATE, a constraint-based steady-state method, and scFEA, a framework based on artificial neural networks designed to overcome the limitations of constraint-based methods. Subsequently, I used intracellular metabolomics of the same set of cells to establish the evaluation benchmark for the predicted fluxes. Specifically, the metabolic cancer phenotypes identified through cluster analysis were compared to those selected by the benchmarks to evaluate their concordance using external metrics such as NMI. INTEGRATE generated cell-specific metabolic models that were sampled by two sampling algorithms: CHRR, which explores the internal feasible region, and CBS, which focuses on the corners of the feasible region. I used summary statistics, such as mean and median, to condense the information gathered through flux sampling. On the other hand, scFEA directly generated a compact flux distribution for each cell. Special attention 93 was paid to the hyperparameter optimization of the clustering algorithms (K-means, K-medoids, Agglomerative, HDBSCAN, and Leiden) and the preprocessing of the RNA-seq matrix with the MAGIC denoiser. The cluster analysis of the intracellular metabolomics robustly revealed the presence of two subpopulations. This result is supported by optimal values in several internal evaluation metrics (Distortion, Silhouette, CalinskiHarabasz, and DBCV). The final benchmark selection was based on the highest Silhouette score of 0.30 (see Figure 3.1), identified by the K-means algorithm. The clusters are imbalanced, with 64.71% of cells in Cluster 0 and 35.29% in Cluster 1. Metabolism is regulated by both genes encoding enzymes that drive reactions and by intracellular metabolite concentrations necessary for metabolic fluxes. This thesis revealed that, within the analyzed dataset, transcriptomics and metabolomics exhibit similar metabolic phenotypes, as evidenced by an NMI of 0.42 between their clusters. Amongst the RNA-seq clusterings, the one with the highest NMI, characterized by two clusters and identified using the Leiden algorithm, was selected (see Figure 3.2). Despite the high level of coherence between metabolomics and transcriptomics, the fluxomics predicted by INTEGRATE did not accurately reflect the metabolic phenotypes identified by the benchmark, as evidenced by low NMIs across all flux statistics (see Figure 3.3). The fluxes showed low NMI values with respect to both intracellular metabolomics and even RNA-seq and RAS matrix where a strong dependence is trivially defined by the INTEGRATE algorithm. However, this discrepancy was not attributable to limitations of the INTEGRATE method itself. In fact, INTEGRATE generated 513 distinct feasible regions, with less than 10% overlap for CBS and 6% overlap for CHRR in shared feasible solutions, demonstrating that the polytopes derived from the integrated RNA-seq data were well-characterized. The real cause of such poor results was the noise contained in the RNA-seq matrix. Originally designed for single-cell RNA-seq, which is characterized by high dropout rates, the MAGIC algorithm reduced the sparsity percentage from 43% to 0.02% for the full gene matrix and from 0.08% to 0% for the matrix containing only metabolic genes (see Figure 3.4). MAGIC did not introduce numerical artifacts, as the transcriptomics NMI with the metabolomics benchmark remained consistent while improving clustering quality by creating more compact clusters. This improvement is attributed to MAGIC’s ability to enhance similarities and differences between cells through the diffusion process in the K-NN graph. Graphical examples are shown in Figures 3.8 and 3.9 where, even without MAGIC, well-clustered fluxes are shown. Eventually, the low fluxomics NMIs were primarily due to noise present in the RNAseq matrix. Although MAGIC was not applied, a sufficient NMI was still observed between transcriptomics and metabolomics benchmarks. However, the fluxes remained significantly influenced by transcriptomic noise. This 94 issue arises because the flux solutions are derived from a combination of multiple constraints based on the original RNA-seq matrix, and the presence of even small errors is amplified, resulting in poor-quality fluxomics. For CBS, the mean fluxes increased from 0.10 to 0.30, the median from 0.08 to 0.39, and the mode from 0.08 to 0.32. For CHRR, the mean went from 0.09 to 0.33, the median from 0.09 to 0.32, and the mode from 0.08 to 0.29. In general, the most promising fluxomics summary statistics are the fluxes median and Biomass correlation for CBS (NMIs = 0.39), and mode for CHRR (NMI = 0.33). Additionally, MAGIC improved cluster quality, leading to higher Silhouette scores as well (see Figures 3.11). The denoiser effects were also clear over CHRR and CBS subsamples (50 flux samples per cell) where increased the NMI between benchmark and CHRR from 0.05 to 0.25 and 0.01 to 0.05 for CBS (see Figures 3.12 and 3.13). The CBS algorithm collected extremely poor results in both cases, while this did not happen in the case of flux summary statistics. The most probable hypothesis for such behavior is based on the fact that CBS distributions are likely to be multi-modal, as discussed by [ 1 ] and shown by 3.14, and thus the clustering algorithms might clusterize the modalities rather than the whole flux distribution. scFEA, which produces a compact flux distribution for each cell, achieved an NMI of 0.38 with the raw RNA-seq matrix, outperforming INTEGRATE fluxomics, whose highest NMI was 0.10, achieved by the CBS mean. The differences between scFEA and INTEGRATE nearly disappeared when applying MAGIC. Specifically, scFEA achieved an NMI of 0.40, while INTEGRATE fluxomics closely followed with an NMI of 0.39, achieved by CBS median and Biomass correlation. scFEA is immune to MAGIC effects due to its intrinsic algorithm structure, which is based on a loss function that maximizes the correlation between gene expression and fluxes. In essence, scFEA automatically performs the role of MAGIC. However, scFEA has some drawbacks. For example, it computes the fluxes of a restricted set of reactions, and it is not yet possible to use a customized metabolic model or growth medium to shape the metabolic pathway. To summarize, the analyses I performed in this thesis demonstrated a moderate level of concordance between intracellular metabolomics and transcriptomics. These findings led to the proposal of a new evaluation pipeline for selecting the optimal integration framework based on metabolomics as benchmark. Initially, the INTEGRATE framework was damaged by noise in the transcriptomics matrix, resulting in low-quality fluxomics and consequently low concordance levels with both benchmarks and transcriptomics. However, the MAGIC algorithm significantly improved the performance of INTEGRATE fluxomics. On the other hand, scFEA proved to be extremely promising, acting effectively as a denoiser with its innovative approach that overcomes traditional constraint-based limitations, such as linear constraints 95 and steady state. Overall, neither INTEGRATE nor scFEA drastically outperformed the other in terms of NMI when using MAGIC. The final benchmarking pipeline was established by utilizing the best intracellular metabolomics clustering (in this case, the highest Silhouette score) to evaluate the subpopulations identified by all fluxomics datasets produced by integration methods like INTEGRATE and scFEA. The identified clusters can be examined by expert life scientists to extract valuable insights about the biomarkers of each subpopulation. An example of such an operation is shown in the RNA-seq clustering, which revealed an interesting pattern based on cancer histology, with one group containing the majority of solid tumors and another group containing liquid ones. This is a known property of cancer cells, as discussed by [55]. Future research extensions may include the analysis of a wider range of integration frameworks, such as iMAT [ 15 ] and FASTCORE [ 16 ], to increase the robustness of the benchmarking pipeline. Additionally, more advanced methods for dimensionality reduction, such as ISOMAP [ 56 ] and Autoencoders [ 57 ], could be explored. Implementing classifiers to adjust clustering hyperparameters to optimize the concordance level with the benchmark is another promising avenue for future work. The INTEGRATE pipeline and clustering comparison implemented in this thesis are going to be deployed in the Google Summer Of Code 2024, as this project was selected by the National Resource for Network Biology (NRNB) under the title: ”COBRAxy: COBRA and MaREA4Galaxy”. COBRAxy is going to be a new Galaxy tool [ 58 ] (Galaxy is an open-source, web-based platform for data-intensive biomedical research) that will extend the already existing MaREA4Galaxy [ 59 ] (a Galaxy tool that implements, amongst other things, the RAS computation and clustering). COBRAxy deploys the complete INTEGRATE pipeline: cell-specific COBRA models, flux sampling, summary flux statistics, clustering grid search and comparison. 96 Bibliography [1] Bruno G. Galuzzi, Luca Milazzo, and Chiara Damiani. Adjusting for false discoveries in constraint-based differential metabolic flux analysis. Journal of Biomedical Informatics, 150:104597, 2024. ISSN 1532-0464. doi: https://doi.org/10.1016/j.jbi.2024.104597. URL https://www. sciencedirect.com/science/article/pii/S1532046424000157. [2] Chiara Damiani, Daniela Gaglio, Elena Sacco, Lilia Alberghina, and Marco Vanoni. Systems metabolomics: From metabolomic snapshots to design principles. Current opinion in biotechnology, 63:190–199, 2020. [3] Yue Gong, Peng Ji, Yun-Song Yang, Shao Xie, Tian-Jian Yu, Yi Xiao, Ming-Liang Jin, Ding Ma, Lin-Wei Guo, Yu-Chen Pei, et al. Metabolicpathway-based subtyping of triple-negative breast cancer reveals potential therapeutic targets. Cell metabolism, 33(1):51–64, 2021. [4] Zachary E Stine, Zachary T Schug, Joseph M Salvino, and Chi V Dang. Targeting cancer metabolism in the era of precision oncology. Nature Reviews Drug Discovery, 21(2):141–162, 2022. [5] Jordi Barretina, Giordano Caponigro, Sungjoon Kim, Nicolas Stransky, Kavitha Venkhatesan, Venkateshwar Reddy, Michael Berger, Michael Morrissey, Paula Morais, Jodi Meltzer, et al. The cancer cell line encyclopedia project: From integrative cancer genomics to personalized cancer therapy. Cancer Research, 70(8 Supplement):2620–2620, 2010. [6] Marzia Di Filippo, Dario Pescini, Bruno Giovanni Galuzzi, Marcella Bonanomi, Daniela Gaglio, Eleonora Mangano, Clarissa Consolandi, Lilia Alberghina, Marco Vanoni, and Chiara Damiani. INTEGRATE: Model-based multi-omics data integration to characterize multi-level metabolic regulation. PLoS Comput Biol, 18(2):e1009337, February 2022. [7] Norah Alghamdi, Wennan Chang, Pengtao Dang, Xiaoyu Lu, Changlin Wan, Silpa Gampala, Zhi Huang, Jiashi Wang, Qin Ma, Yong Zang, et al. A graph neural network model to estimate cell-wise metabolic flux using single-cell rna-seq data. Genome research, 31(10):1867–1884, 2021. 97 [8] David van Dijk, Roshan M. Sharma, David M. Callahan, Karsten A. Balderrama-Gutierrez, Jonathan M. Pachitariu, Yaniv Zaitlen, Yoav Pe’er, Thomas J. Wolf, and Smita Krishnaswamy. Magic: A diffusionbased imputation method reveals gene-gene interactions in single-cell rna-sequencing data. Cell, 174(3):716–729, 2018. doi: 10.1016/j.cell. 2018.05.061. [9] Jens Nielsen. Systems biology of metabolism. Annual review of biochemistry, 86:245–275, 2017. [10] Elizabeth Brunk, Swagatika Sahoo, Daniel C Zielinski, Ali Altunkaya, Andreas Dr¨ager, Nathan Mih, Francesco Gatto, Avlant Nilsson, German Andres Preciat Gonzalez, Maike Kathrin Aurich, et al. Recon3d enables a three-dimensional view of gene variation in human metabolism. Nature biotechnology, 36(3):272–281, 2018. [11] Marta Cascante and Silvia Marin. Metabolomics and fluxomics approaches. Essays Biochem, 45:67–81, 2008. [12] Zhong Wang, Mark Gerstein, and Michael Snyder. Rna-seq: a revolutionary tool for transcriptomics. Nature reviews genetics, 10(1):57–63, 2009. [13] Daniel R Hyduke, Nathan E Lewis, and Bernhard Ø Palsson. Analysis of omics data with genome-scale models of metabolism. Mol Biosyst, 9 (2):167–174, December 2012. [14] Sjoerd Opdam, Anne Richelle, Benjamin Kellman, Shanzhong Li, Daniel C Zielinski, and Nathan E Lewis. A systematic evaluation of methods for tailoring Genome-Scale metabolic models. Cell Syst, 4 (3):318–329.e6, February 2017. [15] Hadas Zur, Eytan Ruppin, and Tomer Shlomi. iMAT: an integrative metabolic analysis tool. Bioinformatics, 26(24):3140–3142, November 2010. [16] Nikos Vlassis, Maria Pires Pacheco, and Thomas Sauter. Fast reconstruction of compact context-specific metabolic network models. PLOS Computational Biology, 10(1):1–10, 01 2014. doi: 10.1371/journal.pcbi. 1003424. URL https://doi.org/10.1371/journal.pcbi.1003424. [17] Bruno G. Galuzzi, Luca Milazzo, and Chiara Damiani. Best practices in flux sampling of constrained-based models. Accepted for The 8th International Online Onsite Conference on Machine Learning, Optimization, and Data Science (LOD22), 2022. 98 [18] Jan Schellenberger, Richard Que, Ronan M. T. Fleming, Ines Thiele, Jeffrey D. Orth, Adam M. Feist, Daniel C. Zielinski, Aarash Bordbar, Nathan E. Lewis, Sorena Rahmanian, Joseph Kang, Daniel R. Hyduke, and Bernhard Ø Palsson. Quantitative prediction of cellular metabolism with constraint-based models: the cobra toolbox v2.0. Nature Protocols, 6(9):1290–1307, Sep 2011. ISSN 1750-2799. doi: 10.1038/nprot.2011.308. URL https://doi.org/10.1038/nprot.2011.308. [19] Ali Ebrahim, Joshua A. Lerman, Bernhard O. Palsson, and Daniel R. Hyduke. Cobrapy: Constraints-based reconstruction and analysis for python. BMC Systems Biology, 7(1):74, Aug 2013. ISSN 17520509. doi: 10.1186/1752-0509-7-74. URL https://doi.org/10.1186/ 1752-0509-7-74. [20] Stuart Lloyd. Least squares quantization in pcm. IEEE transactions on information theory, 28(2):129–137, 1982. [21] Hae-Sang Park and Chi-Hyuck Jun. A simple and fast algorithm for k-medoids clustering. Expert systems with applications, 36(2):3336–3341, 2009. [22] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825–2830, 2011. [23] Vincent A Traag, Ludo Waltman, and Nees Jan Van Eck. From louvain to leiden: guaranteeing well-connected communities. Scientific reports, 9(1):5233, 2019. [24] Davoud Moulavi, Pablo A Jaskowiak, Ricardo JGB Campello, Arthur Zimek, and J¨org Sander. Density-based clustering validation. In Proceedings of the 2014 SIAM international conference on data mining, pages 839–847. SIAM, 2014. [25] Changdai Gu, Gi Bae Kim, Won Jun Kim, Hyun Uk Kim, and Sang Yup Lee. Current status and applications of genome-scale metabolic models. Genome Biology, 20(1):121, Jun 2019. ISSN 1474-760X. doi: 10.1186/s13059-019-1730-3. URL https://doi.org/10.1186/ s13059-019-1730-3. [26] Jeffrey D. Orth, Ines Thiele, and Bernhard Ø Palsson. What is flux balance analysis? Nature Biotechnology, 28(3):245–248, Mar 2010. ISSN 1546-1696. doi: 10.1038/nbt.1614. URL https://doi.org/10.1038/ nbt.1614. 99 [27] Alex Graudenzi, Davide Maspero, Marzia Di Filippo, Marco Gnugnoli, Claudio Isella, Giancarlo Mauri, Enzo Medico, Marco Antoniotti, and Chiara Damiani. Integration of transcriptomic data and metabolic networks in cancer samples reveals highly significant prognostic power. Journal of Biomedical Informatics, 87:37–49, 2018. ISSN 1532-0464. doi: https://doi.org/10.1016/j.jbi.2018.09.010. URL https://www. sciencedirect.com/science/article/pii/S1532046418301849. [28] Sophia Malandraki-Miller, Colleen Lopez, Heba Al-Siddiqi, and Carolyn Carr. Changing metabolism in differentiating cardiac progenitor cells—can stem cells become metabolically flexible cardiomyocytes? Frontiers in Cardiovascular Medicine, 5, 09 2018. doi: 10.3389/fcvm.2018.00119. [29] Gerald J. Lieberman Frederick S. Hillier. Introduction to Operations Research, ninth edition. McGraw-Hill Education, 2009. [30] Shirin Fallahi, Hans J. Skaug, and Guttorm Alendal. A comparison of monte carlo sampling methods for metabolic network models. PLOS ONE, 15(7):1–24, 07 2020. doi: 10.1371/journal.pone.0235393. URL https://doi.org/10.1371/journal.pone.0235393. [31] Yoann Altmann, Steve McLaughlin, and Nicolas Dobigeon. Sampling from a multivariate gaussian distribution truncated on a simplex: a review. In IEEE Workshop on Statistical Signal Processing Proceedings, pages 113–116, United States, 2014. IEEE. ISBN 9781479949755. doi: 10.1109/SSP.2014.6884588. 17th IEEE Workshop on Statistical Signal Processing 2014, SSP 2014 ; Conference date: 29-06-2014 Through 02-07-2014. [32] Walter R. Gilks, Sylvia Richardson, and David J. Spiegelhalter. Markov chain concepts related to sampling algorithms. 1995. [33] David E. Kaufman and Robert L. Smith. Direction choice for accelerated convergence in hit-and-run sampling. Operations Research, 46(1):84–95, 1998. doi: 10.1287/opre.46.1.84. URL https://pubsonline.informs. org/doi/abs/10.1287/opre.46.1.84. [34] Hulda S Haraldsd´ottir, Ben Cousins, Ines Thiele, Ronan M.T Fleming, and Santosh Vempala. CHRR: coordinate hit-and-run with rounding for uniform sampling of constraint-based models. Bioinformatics, 33 (11):1741–1743, 01 2017. ISSN 1367-4803. doi: 10.1093/bioinformatics/ btx052. URL https://doi.org/10.1093/bioinformatics/btx052. [35] J. Telgen and A. Boneh. Telgen j. private communication with a. boneh., 1980. 100 FVA standard 272 22 FVA magic minmax 254 4 FVA magic standard 254 5 pFBA minmax 161 3 pFBA standard 161 7 pFBA magic minmax 109 2 pFBA magic standard 109 2 fluxes CBS stats magic minmax 727 6 fluxes CBS stats magic standard 727 9 fluxes CBS stats minmax 876 27 fluxes CBS stats standard 876 43 fluxes CHRR stats magic minmax 1181 4 fluxes CHRR stats magic standard 1181 7 fluxes CHRR stats minmax 1181 4 fluxes CHRR stats standard 1181 7 sensitivity CHRR minmax 238 48 sensitivity CHRR standard 238 51 sensitivity CBS minmax 248 24 sensitivity CBS standard 248 29 sensitivity CHRR magic minmax 232 45 sensitivity CHRR magic standard 232 48 sensitivity CBS magic minmax 245 10 sensitivity CBS magic standard 245 12 sensitivityKnockOut magic minmax 11 3 sensitivityKnockOut magic standard 11 3 sensitivityKnockOut minmax 38 5 sensitivityKnockOut standard 38 12 107 A.2 Metabolomics biomarkers and meta-data Figure A.1: Scaled heatmap of metabolomics markers (clustering selected by Silhouette) where columns represent metabolites, rows are associated to cells and colors are proportional to metabolites concentrations. 108 Figure A.2: Histology distribution of metabolomics clustering. 109 Figure A.3: Pathology distribution of metabolomics clustering. 110 Figure A.4: Primary site distribution of metabolomics clustering. 111 A.3 Clustering hyperparameters Table A.2: Clusterings selected by highest NMI with the metabolomics benchmark. Omic Algorithm Clusters Silhouette NMI RNA-seq metabolic Leiden 2 0.08 0.42 RAS Leiden 2 0.07 0.40 pFBA HDBSCAN 2 0.17 0.05 Sensitivity analysis K-means 10 0.24 0.09 CBS mean Leiden 8 0.04 0.10 CBS mode K-medoids 5 -0.02 0.08 CBS median Leiden 32 -0.03 0.08 Correlation CBS Leiden 25 0.11 0.07 CHRR mean Leiden 7 0.05 0.09 CHRR mode Leiden 29 0.03 0.08 CHRR median Leiden 12 0.06 0.09 Correlation CHRR K-means 2 0.07 0.08 Table A.3: Clusterings selected by highest NMI with the metabolomics benchmark. In this case, the RNA-seq matrix was processed by MAGIC. Omic Algorithm Clusters Silhouette NMI RNA-seq metabolic kmedoids 2 0.38 0.41 RAS kmedoids 2 0.36 0.42 pFBA agglomerative 9 0.07 0.23 Sensitivity analysis agglomerative 3 0.29 0.26 CBS mean leiden 7 0.42 0.30 CBS mode leiden 5 0.24 0.32 CBS median agglomerative 2 0.34 0.39 Correlation CBS agglomerative 2 0.28 0.39 CHRR mean agglomerative 3 0.35 0.33 CHRR mode leiden 7 0.38 0.29 CHRR median leiden 6 0.38 0.32 Correlation CHRR agglomerative 5 0.11 0.28 112