Systems-level modelling of the cancer and immune metabolome to improve immunotherapeutic outcomes

Universidade do Minho Escola de Engenharia Sara Manso de Sousa Cardoso Systems-level modelling of the cancer and immune metabolome to improve immunotherapeutic outcomes Doctorate Thesis Doctorate in Biomedical Engeneering Work developed under the supervision of: Miguel Rocha Noel de Miranda Dina Ruano October, 2022 COPYRIGHT AND TERMS OF USE OF THIS WORK BY A THIRD PARTY This is academic work that can be used by third parties as long as internationally accepted rules and good practices regarding copyright and related rights are respected. Accordingly, this work may be used under the license provided below. If the user needs permission to make use of the work under conditions not provided for in the indicated licensing, they should contact the author through the RepositoriUM of Universidade do Minho. License granted to the users of this work Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International CC BY-NC-SA 4.0 ?iiTb,ff+`2iBp2+QKKQMbXQ`;fHB+2Mb2bf#v@M+@bf9Xyf/22/X2M This document was created with the (pdf/Xe/Lua)L A T EX processor and the NOVAthesis template (v6.7.4) [1]. ii Acknowledgements I would like to express my gratitude to everyone who, directly or indirectly, contributed to the completion of this work. First of all, I would like to thank Fundação para a Ciência e Tecnologia for the PhD scholarship I was awarded (SFRH/BD/138951/2018), without which this work would not be possible. I would also like to thank my main host institution, the Centre of Biological Engineering at the University of Minho, and the Leiden University Medical Center (LUMC), my secondary institution. I would also like to thank my supervisors, for offering me this opportunity and all the support and guidance. Miguel Rocha, my main supervisor, for not just the past four years of working together, but also those that came before. Noel de Miranda and Dina Ruano, my co-supervisors, for their feedback on this work and for welcoming me in the group at Leiden. I would also like to thank everyone that has passed through the BISBII research group at the University of Minho while I was here. To everyone at the LUMC, but more specifically those from the pathology group, thank you all. I would like to thank all my friends for supporting me all these years. Last, and by no means least, my family. No words would be enough to thank them for all the years that lead to now. Once again, thank you all. iii STATEMENT OF INTEGRITY I hereby declare having conducted this academic work with integrity. I confirm that I have not used plagiarism or any form of undue use of information or falsification of results along the process leading to its elaboration. I further declare that I have fully acknowledged the Code of Ethical Conduct of the Universidade do Minho. , (Place) (Date) (Sara Manso de Sousa Cardoso) iv Resumo Modelação do metaboloma do cancro e do sistema imunitário de forma a melhorar resultados imunoterapêuticos A medicina de precisão busca fornecer terapias para diversas doenças que sejam adequadas para (grupos de) pacientes específicos. O cancro é uma doença causada por células anormais que se multiplicam descontroladamente e, dada a sua heterogeneidade e base genética, é um dos desafios mais relevantes para a medicina de precisão. As imunoterapias podem ser adaptadas para indivíduos específicos, fornecendo formas interessantes de combater o cancro, induzindo ou aprimorando as respostas naturais do sistema imunológico dos pacientes, o que pode traduzir-se em terapias com menos efeitos colaterais. Neste trabalho, o objectivo foi o de desenvolver abordagens computacionais baseadas em mineração de dados ómicos e modelação metabólica que ajudem os esforços personalizados de descoberta de medicamentos com o mínimo de efeitos secundários. Para isso, foram desenvolvidos modelos metabólicos de células T de tumores e tecidos saudáveis com base em dados ómicos de single-cell .O single-cell RNAseq ( scRNAseq ) é uma ótima ferramenta para a reconstrução de modelos metabólicos específicos para cada paciente e tipo de célula. No entanto, esta abordagem não foi ainda aproveitada no campo da imunoterapia tumoral. Numa fase inicial, construímos um atlas de dados de scRNAseq para cancro colorretal, que foi usado para reconstruir 196 modelos de vários sub-tipos de células T do micro-ambiente desse tumor. Além disso, realizou-se uma análise do desempenho de vários métodos de deconvolução do tumor, para permitir que os dados de bulk RNAseq amplamente disponíveis sejam usados na modelação metabólica de tipos de células presentes no micro-ambiente do cancro colorretal. Palavras-chave: Modelos metabólicos, Single-cell RNAseq , Cancro colorectal, Deconvolução de tumores v Abstract Systems-level modelling of the cancer and immune metabolome to improve immunotherapeutic outcomes Precision medicine seeks to provide therapies for different diseases that are adequate for specific (groups of) patients. Cancer is a a disease caused by abnormal cells that multiply uncontrollably and given its heterogeneity and genetic basis, is one of the most relevant challenges for precision medicine. Immunotherapies can be tailored for specific subjects, providing attractive ways to fight cancer by inducing or enhancing patients’ natural immune system responses, which can lead to therapies with less side effects. In this work, we aimed to develop computational approaches based on omics data mining and metabolic modelling that could help, in the future, personalised drug discovery efforts with minimum off-target effects. For this, metabolic models of T-cells from tumour and healthy tissues based on single-cell omics data were developed. Single-cell RNAseq is a great tool for the reconstruction of metabolic models with patientand cell-typescpecificity. However, this approach has not been exploited in the field of tumour immunotherapy. We constructed an atlas of scRNAseq data for colorectal cancer, which was used to reconstruct 196 models of various T-cell subtypes from the micro-environment of this tumour. Furthermore, the benchmarking of several tumour deconvolution methods was performed to allow the extensively available bulk RNAseq data to be optimally used in cell-type specific modeling of the colorectal cancer. Keywords: Metabolic modeling, Single-cell RNAseq, Colorectal cancer, Tumour deconvolution vi Contents List of Figures x List of Tables xvii Acronyms xix 1 Introduction 1 1.1 Context and Motivation .............................. 1 1.2 Research Objectives ................................ 2 1.3 Thesis Outline .................................. 2 2 Background 4 2.1 Overview of the Immune System .......................... 4 2.1.1 From innate to adaptive immune system ................. 4 2.1.2 T-cells .................................. 5 2.2 Cancer ...................................... 7 2.2.1 Cancer hallmarks ............................ 8 2.2.2 Role of the Immune System in Cancer .................. 10 2.3 Metabolism in T-cells and Cancer ......................... 11 2.3.1 T-cells .................................. 11 2.3.2 Tumour cells ............................... 14 2.3.3 Tumour Micro-Environment ........................ 16 2.3.4 Targeting metabolism for therapy ..................... 18 2.4 Constraint-based Modeling and Human Cancer .................. 19 2.4.1 Stoichiometric modeling . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.2 Metabolic Models in Humans ....................... 24 2.4.3 OmicsData ............................... 26 2.4.4 Applications of Modeling in Human Cancer . . . . . . . . . . . . . . . . 29 vii 44 (A) Overall correlation vs RMSE for all methods, including the methods combining the best methods for each cell-type ( Combined and Combined_norm ). (B) Sample correlation vs RMSE for Scaden , Combined and Combined_norm . A good method has a high correlation and a small RMSE. Grey horizontal and vertical lines mark a correlation of 0.5 and a RMSE of 0.2, respectively. Scatter plots of estimated vs ground-truth proportions for (C) Combined and (D) Combined_norm ................................ 92 45 Pearson correlation (A) and (B) RMSE values of the methods without and with correction. Methods were tested for correction on the reference matrix ( Corrected Before ) and on the estimated proportions ( Corrected After )......................... 93 46 RMSE (A) and (B) pearson correlation values of the methods without and with correction, separated by cell-type. Methods were tested for correction on the reference matrix ( Corrected Before ) and on the estimated proportions ( Corrected After ). .............. 94 47 For each method, RMSE of the samples is mapped against their corresponding (A) total number of cells and (B) total number of read counts. . . . . . . . . . . . . . . . . . . 95 48 For each method, RMSE of the samples is mapped against their corresponding proportion ofcancercells..................................... 96 49 For each method, RMSE of the samples is mapped against their corresponding proportion of other cells. ..................................... 97 50 For each method, RMSE of the samples is mapped against their corresponding proportion of immune cells. ................................... 97 51 Clusterability results from SIGMA. Each dot corresponds to a cell, which are coloured by dataset of origin. The clusterability of the clusters was not dictated by the dataset of origin. 121 52 Heatmap of CNV predictions for the tumour cells of patient KUL21. . . . . . . . . . . . 122 53 Heatmap of CNV predictions for the tumour cells of patient SMC04. . . . . . . . . . . 123 54 Heatmap of CNV predictions for the tumour cells of patient SMC07. .......... 124 55 Heatmap of CNV predictions for the tumour cells of patient SMC10. . . . . . . . . . . 125 56 Similarity between (A) the structure of the models (i.e., reaction presence/absence) and (B) predicted fluxes under normal human blood medium. The smaller Euclidean distance is, the smallerthesimilarityis................................. 126 57 Pathway coverage (%). ................................. 127 58 Biomass (A) and ATP (B) production when biomass was set as the only objective for all models......................................... 128 59 Cumulative fluxes (mmol/gDW/h) of the reactions that produce NADH, from all source pathways, for proliferative CD4 and CD8 T-cell models . . . . . . . . . . . . . . . . . . . 128 60 Distribution of (A) DNA and (B) RNA production of the T-cell types, with and without glutamine inthemedium..................................... 129 61 Distribution of DNA production of the T-cell types, with and without nucleotides in the medium. 130 xiv 62 Distribution of the biomass flux of the T-cell types, with and without glucose in the medium. 130 63 (A) Number of models where biomass flux increases, decreases or suffers no change. This information is further showed by (B) tissue of origin, (C) CMS type, and (D) cell-type. .. 131 64 Scatter plots of estimated vs ground-truth proportions for the remaining methods that are not presentinfigure42. ................................. 132 65 Scatter plots of estimated vs ground-truth proportions of all methods for cancer cells. . . 133 66 Scatter plots of estimated vs ground-truth proportions of all methods for stromal cells. .. 134 67 Scatter plots of estimated vs ground-truth proportions of all methods for macro/mono lineage cells.......................................... 135 68 Scatter plots of estimated vs ground-truth proportions of all methods for B-cells. . . . . . 136 69 Scatter plots of estimated vs ground-truth proportions of all methods for CD4 T-cells. .. 137 70 Scatter plots of estimated vs ground-truth proportions of all methods for regulatory CD4 Tcells.......................................... 138 71 Scatter plots of estimated vs ground-truth proportions of all methods for CD8 T-cells. . . 139 72 Scatter plots of estimated vs ground-truth proportions of all methods for proliferative T-cells. 140 73 Scatter plots of estimated vs ground-truth proportions of all methods for NK cells. .... 141 74 For cancer cells, scatter plots of estimated vs ground-truth proportions of best uncorrected method and the corrected methods whose RMSE improved relative to the best uncorrected method. ....................................... 141 75 For stromal cells, scatter plots of estimated vs ground-truth proportions of best uncorrected method and the corrected methods whose RMSE improved relative to the best uncorrected method. ....................................... 142 76 For macro/mono lineage cells, scatter plots of estimated vs ground-truth proportions of best uncorrected method and the corrected methods whose RMSE improved relative to the best uncorrectedmethod. ................................. 142 77 For B-cells, scatter plots of estimated vs ground-truth proportions of best uncorrected method and the corrected methods whose RMSE improved relative to the best uncorrected method. 143 78 For CD4 T-cells, scatter plots of estimated vs ground-truth proportions of best uncorrected method and the corrected methods whose RMSE improved relative to the best uncorrected method. ....................................... 143 79 For regulatory CD4 T-cells, scatter plots of estimated vs ground-truth proportions of best uncorrected method and the corrected methods whose RMSE improved relative to the best uncorrectedmethod. ................................. 143 80 For CD8 T-cells, scatter plots of estimated vs ground-truth proportions of best uncorrected method and the corrected methods whose RMSE improved relative to the best uncorrected method. ....................................... 144 xv 81 For proliferative T-cells, scatter plots of estimated vs ground-truth proportions of best uncorrected method and the corrected methods whose RMSE improved relative to the best uncorrected method. ................................. 144 82 For NK cells, scatter plots of estimated vs ground-truth proportions of best uncorrected method and the corrected methods whose RMSE improved relative to the best uncorrected method. ....................................... 144 83 Scatter plot of samples’ total read counts vs total cell counts. . . . . . . . . . . . . . . 145 84 For each method, RMSE of the samples is mapped against their corresponding proportion ofstromalcells..................................... 145 85 For each uncorrected method, scatter plot of pearson correlation vs RMSE of the samples. 146 86 For Scaden , scatter plots of pearson correlation vs RMSE of the samples for each correction type.......................................... 146 xvi List of Tables 1 Overview of the datasets used to construct the CRC atlas . . . . . . . . . . . . . . . . 32 2 Number of cells in each cell-subtype present in the CRC atlas. TA : transit-amplifying cells; CAFs : Cancer-associated fibroblasts; VSMCs : Vascular smooth muscle cells; DCs : Dendritic cells; LTi : Lymphoid tissue-inducer cells. . . . . . . . . . . . . . . . . . . . . . . . . 45 3 Calculation of the gene activity scores (GAS) from a sample. 𝑥𝑔: expression, in CPMs, of gene 𝑔;𝑔𝑙𝑜𝑏𝑎𝑙_𝑚𝑎𝑥: 75th percentile of the expression distribution of all genes in all celltypes of a sample; 𝑔𝑙𝑜𝑏𝑎𝑙_𝑚𝑖𝑛: 10th percentile of the expression distribution of all genes in all cell-types of a sample; 𝑙𝑜𝑐𝑎𝑙_𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑: 35th percentile of the expression distribution of a gene in all cell-types of a sample. . . . . . . . . . . . . . . . . . . . . . . . . . 50 4 Examples for calculating RASs. 𝑥𝑎,𝑥𝑏,𝑥𝑐: expression, in CPMs, of genes 𝑎,𝑏and 𝑐, respectively; 𝑚𝑎𝑥: maximum; 𝑚𝑖𝑛: minimum. .................... 51 5 For each cell-type ( Cell Type) , number of reconstructed models (Number of Models) and their distribution regarding tissue of origin (State Distribution) and CMS classification (CMS Distribution) . ..................................... 55 6 Top twenty pathways with the most median percentage of essential genes. ....... 76 7 Potential essential genes from the Eicosanoid metabolism pathway and respective products. 76 8 Top twenty pathways with the most median number of essential genes. ......... 78 9 Bias factors used for RNA content bias correction. ................... 87 10 Best methods for each cell-type. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 11 Metabolites used in the human blood medium, with the corresponding exchange reaction id from the model. Each metabolite has information on the average concentration (mM) in normal human blood, gathered from SMDB database, and the fluxes (mmol/gDW/h) used in normal and tumour human blood media. ...................... 147 xvii 12 Fold changes between tumour and normal blood reported by the studies and calculated average fold change. Study sinalised with a †reported results for both GC-TOFMS and GCMS-QP201. .................................... 158 13 Number of cells present in each sample for each T-cell subtype considered. Those subtypes with 5 or less cells in a sample were not considered for model reconstruction in that sample. Those with a †did not pass the gap-fill and were not analysed. The last line is the total number of models reconstructed for each T-cell subtype. . . . . . . . . . . . . . . . . 160 14 Essential genes that catalase uptake of metabolites. Cell-types codes: Cyto : cytotoxic CD8; Fol : follicular CD4; IL17 : IL17+ CD4; Mem4 : memory CD4; Mem8 : memory CD8; N4 : naive CD4; N8 ; Prol4 : proliferative CD4; Prol8 : proliferative CD8; Regs : regulatory CD4. The essentiality reported by the two CRISPR-Cas9 studies is provided: - : gene not tested in study; Essential : gene tested and reported as essential; Not essential : gene tested and reported as not essential. ............................... 161 15 Differentially covered pathways between normaland tumourderived regulatory CD4 T-cell models......................................... 162 16 Differentially covered pathways between normaland tumourderived cytotoxic CD8 T-cell models......................................... 163 17 Differentially covered pathways between naive and proliferative CD8 T-cell models. . . . 163 18 Differentially covered pathways between IL17+ and regulatory CD4 T-cell models. .... 164 19 Cell-types used for tumour deconvolution, and respective overlap with CRC atlas and groundtruth phenotypes. ................................... 164 xviii Acronyms APC Antigen-presenting cell ATP Adenosine triphosphate CAF Cancer-associated fibroblasts CBM Constraint-based modeling CMS Consensus molecular subtype CNV Copy number variation CPM Counts per million CRC Colorectal cancer EFM Elementary flux mode EMT Epithelial-mesenchymal transition ETC Electron transport chain FA Fatty acid FADH2 Flavin adenine dinucleotide FAO Fatty acid oxidation FAS Fatty acid synthesis FBA Flux balance analysis FVA Flux variability analysis GAS Gene activity score GPR Gene-protein-reaction GSMM Genome-scale metabolic model IDO Indoleamine-2,3-dioxygenase xix ILC Innate lymphoid cell MCC Mathews correlation coefficient MDSC Myeloid-derived suppressor cell MHC Major histocompatibility complex MSC Mesenchymal stem cell NADH Nicotinamide adenine dinucleotide NO Nitric oxide OXPHOS Oxidative phosphorilation PBMC Peripheral blood mononuclear cell pFBA Parsimonious flux balance analysis PPP Pentose phosphate pathway RAS Reaction activity score RMSE Root mean square error RNAseq RNA sequencing ROS Reactive oxygen species SBML Systems biology markup language scRNAseq single-cell RNA sequencing SRC Spare respiratory capacity TAM Tumour-associated macrophage TCA Tricarboxylic acid cycle TCR T-cell receptor TIL Tumour-infiltrating leukocyte TME Tumour micro-environment UMAP Uniform manifold approximation and projection xx R Introduction 1.1 Context and Motivation Cancer is a disease caused by abnormal cells that multiply uncontrollably, which may generate a solid mass called tumour, and possibly invade the organism. Although tumours are usually portrayed as homogeneous cell populations, they actually contain a diverse amount of cells beyond cancer cells, including stromal and immune cells. In fact, the cells in these tumour micro-environments may be shaped by tumour cells or even by each other to aid tumour proliferation and metastasis, as well as affect responsiveness to therapeutics. Although many characteristics span all or most tumours, the heterogeneity in how, when and what leads these characteristics to manifest are specific to each patient or group of patients. As such, precision medicine seeks to provide therapies for different diseases, including cancer, that are adequate for specific (groups of) patients. Besides finding ways of killing only the malignant cells with minimum impact on healthy cells, natural immune system responses can be induced or enhanced to fight cancer (immunotherapy), which can lead to less side effects. Besides the large number of genomic changes that promote uncontrollable proliferation, cancer cells must undergo metabolic changes, when compared to normal cells, to support the acquisition and maintenance of malignant properties. For instance, changes in energy metabolism such as the shift to aerobic glycolysis in cancer cells even under normal oxygen conditions, as observed firstly by Warburg et al [2] in 1958, allows these cells to rapidly proliferate and survive under stressful conditions. Interestingly, many of the changes observed in the tumour micro-environment are also caused at the metabolic level, providing an interesting target for cancer and immune therapies. The discovery and characterisation of the reprogrammed metabolism of cancer cells and of those in the tumour micro-environment may hence help to study tumour tissue non-invasively, predict tumour behaviour, and prevent tumour progression. The recent major advances in high-throughput omics data, such as genomics, transcriptomics and metabolomics, have allowed modelling the metabolism of several species at a large scale through the reconstruction of Genome-Scale Metabolic Models (GSMMs). Human metabolic models, when applied 1 CHAPTER 1. INTRODUCTION together with data specific to (groups of) patients, can be important contributors in gaining knowledge towards a better comprehension of cancer mechanisms. These models can even be essential, in the future, to find therapies specific for each patient, or group of patients. Thus, this work focused on developing system-level computational approaches based on omics data mining and metabolic modelling, so that efficient immunotherapies, complemented with personalised drug discovery efforts, can be developed in the future with minimum off-target effects. For this, models of T-cells from tumour and healthy tissues based on patientand cell-typespecific omics data were developed. To accomplish this, an atlas of single-cell RNAseq data for colorectal cancer was constructed. Furthermore, benchmarking of tumour deconvolution methods was performed to allow the extensively available bulk RNAseq data to be optimally used in cell-type specific modeling of the colorectal cancer. This will enable the study of the effects of metabolic targets and off-targets on each individual celltype/tissue, as well as interactions between the different cell-types in the tumour micro-environment. 1.2 Research Objectives The main focus of this work is to develop system-level approaches that allow personalised immunotherapies, which will be materialised by the development of computational tools based on constraint-based modelling, omics data mining and immunoinformatics. This work provides innovative contributions along the following axes of research and development: !Reviewing the state of the art on the metabolism of the immune system, cancer, and the tumour micro-environment, as well as on constraint-based modelling and its applications in human and cancer metabolism. !Development of a single-cell RNAseq atlas of colorectal cancer micro-environment, using publicly available datasets; !Development of metabolic models of a set of T-cells, based on the atlas constructed; !Benchmark tumour deconvolution methods using bulk RNAseq data provided by the Leiden University Medical Centre (LUMC); !Implement the code in open-source software, allowing the full reproducibility of the studies conducted, as well as bringing important resources for the scientific community. !Write articles with the results of the work in selected international journals and conferences. 1.3 Thesis Outline The present document is divided into six chapters. The first chapter gives a brief introduction to the subject of the work, by providing a context, motivation and main objectives. 2 1.3. THESIS OUTLINE The following chapter gives a background for this work. It starts by introducing the immune system and cancer, and how these two are connected. This is followed by a description of the metabolism in Tand cancer cells, and how it can be taken advantage of for immunotherapy. Lastly, we detail on what constraint-based modelling is and how it can be used in human cancer. The third chapter describes the creation of a single-cell RNAseq (scRNAseq) atlas of colorectal cancer (CRC) tumour-microenvironment. It starts by explaining what datasets were collected and how they were properly integrated. This is followed by a detailed description of how the cells were annotated, with emphasis on T-cells. Lastly, an overview of the annotated cell-types is given. Chapter four details the creation of genome-scale metabolic models of different T-cell subtypes from the micro-environment of CRC and normal matched colon. The data collected in chapter 3 was used to create these models. This is followed by a thorough analysis of the models, including model structure, flux predictions and gene essentiality. Chapter five details the benchmarking of several tumour deconvolution methods that use single-cell RNAseq data as reference to do so. This chapter includes a summary of each method tested and how they were compared, followed by in-depth analysis of what are the best methods to use in different situations. In the final chapter, the present work is reviewed, followed by a discussion on work to do in the future. 3 CHAPTER 2. BACKGROUND may only contribute to the acquisition of a single hallmark in certain tumours, the same event may lead to simultaneous acquisition of several distinct hallmarks in other tumours [17]. The cell’s paths for becoming malignant vary greatly. Mutations in certain oncogenes and tumour suppressor genes can occur early in some tumours and late in others, which causes the particular sequence in which hallmarks are acquired to vary widely. Nonetheless, the hallmarks that are ultimately reached are shared by virtually all types of tumours [17]. 2.2.2 Role of the Immune System in Cancer The human organism is able to undergo a series of stepwise events called the cancer-immunity cycle so that a productive anti-tumour immunity response can occur. Firstly, neoantigens released from tumours are passively transported or captured and delivered by dendritic cells to regional lymph nodes via afferent lymphatic vessels. In tumour-draining lymph nodes, dendritic cells present the captured cancer-specific antigens to T-cells, activating their responses. At this stage, the balance between T effector cells and T regulatory cells is an important key to the final outcome [26]. The newly activated cytotoxic T-cells exit the lymph node and circulate throughout the body via the bloodstream, whose chemokine gradients and adhesion molecules direct circulating T-cells to extravasate through the blood vessels and migrate into the tumour. Once here, they scan the cancer cells and kill them, which releases additional tumour-associated antigens to increase the range and depth of the response in subsequent steps [26, 27]. While cytotoxic T-cells, and NK cells, engage in tumour killing, T𝐻1and, sometimes, T𝐻17 cells provide important help that boosts cytotoxic immunity [28]. Under ideal circumstances, this cycle leads to the eradication of malignant cells by the cytotoxic cells, establishes tumour-specific immunological memory, and prevents further tumour progression [27]. In the vast majority of established tumours, however, the leukocytes that infiltrate the tumour are insufficient to stop tumour growth [28]. For instance, tumour antigens may not be detected, or dendritic cells and T-cells may treat antigens as self rather than foreign. T-cells may not even properly infiltrate tumours due to factors present in the tumour micro-environment [26]. Also, if T regulatory cells exceed effector responses, they will be able to inhibit these responses against tumour cells. It is clear now that tumours are not an homogeneous cell population of only cancer cells, they are also composed by several types of non-malignant cells. This environment, referred to as tumour microenvironment (TME), includes blood and lymphatic endothelial cells, tumour-infiltrating leukocytes (TILs), mesenchymal stem cells (MSCs) and their differentiated progeny, such as cancer-associated fibroblasts (CAFs) and pericytes, accompanied with the extracellular matrix. The different constituents of the TME can interact closely with each other and the tumour cells. These interactions control and shape tumour cell survival, invasiveness and metastatic dissemination, as well as access and responsiveness to therapeutics. This dynamic relationship is set early on in malignant growth and evolves throughout the life history of a tumour [27, 28]. Besides direct contact and contact through cytokine and chemokine production, metabolism is also an important part in this interaction. 10 2.3. METABOLISM IN T-CELLS AND CANCER 2.3 Metabolism in T-cells and Cancer Glycolysis, oxidative phosphorilation (OXPHOS), fatty acid oxidation (FAO), the tricarboxylic acid (TCA) cycle and fatty acid synthesis (FAS) are core cooperative metabolic pathways in any type of cell. Glycolysis is the pathway by which glucose is broken down into pyruvate, which can then be used as substrate in the TCA cycle in the presence of oxygen, or for production of lactate when the mitochondria are damaged or in the absence of oxygen. However, in 1958, Warburg et al [2] found that activated leukocytes and tumour cells, two types of cells that need to rapidly divide, produce lactate from pyruvate even in the presence of oxygen, a process termed aerobic glycolysis. Although it only produces a net of 2 ATP, aerobic glycolysis provides metabolic intermediates necessary for many other metabolic pathways, allowing these cells to rapidly proliferate. In mitochondria, the TCA cycle incorporates Acetyl-CoA from pyruvate produced by glycolysis or FAO, to generate reducing equivalents (NADH and FADH2), which donate electrons to the electron transport chain (ETC). This leads to the generation of a proton gradient across the inner mitochondrial membrane and ultimately to the generation of up to 36 ATP in a process named OXPHOS. NADH and FADH2can also be produced during FAO. During this process, reactive oxygen species (ROS) are also produced [29]. 2.3.1 T-cells Throughout the life of an immune cell, energy and substrate requirements change considerably, as certain metabolic pathways must be engaged, or suppressed, to facilitate development and activation [29]. Naïve T-cells The naïve 𝐶𝐷8+and 𝐶𝐷4+T-cells formed in the development stage leave the thymus and enter the circulation as resting cells. Traveling throughout the organism on immune surveillance requires constant cytoskeletal rearrangements, a process that is ATP expensive but requires only basal replacement biosynthesis [30]. These cells need a metabolic balance that favours energy production over biosynthesis to move through tissues and prevent cell death, without leaving a quiescence state [29]. To do this, naïve T-cells rely greatly on the high-energy yielding processes of FAO, and pyruvate and glutamine oxidation via the TCA cycle (figure 1)[30]. Naïve T-cells Activation Activated T-cells undergo high production of metabolite precursors and ATP to support the cell growth phase prior to the first T-cell division that will enter the cell in the clonal expansion phase (figure 1). Activated T-cells engage into mitochondrial one-carbon metabolism associated with induction of the folate cycle, which coincides with the increased TCA cycle and OXPHOS [31]. Also, lipid synthesis is suppressed and instead lipid oxidation is promoted [32]. ROS are generated in the mitochondria as a consequence of the TCA cycle and OXPHOS [33]. High levels of ROS lead to uncontrolled proliferation of T-cells. Glutathione, generated through one-carbon units derived from serine in the mitochondria, titrates ROS, maintaining moderate levels of ROS [34]. 11 CHAPTER 2. BACKGROUND Rapid cell division requires biosynthesis of intracellular constituents including lipid membranes, nucleic acids and proteins, which need increased glucose and glutamine to satisfy the metabolic requirements, while decreasing lipid oxidation (figure 1)[30]. Figure 1: Overview of the characteristic metabolic phenotype throughout T-cells life cycle. In green are the metabolic pathways characteristic of quiescent/immunosuppressive cells and in yellow those of proliferating cells. Naïve T-cells and their activated form in the initial growing phase are characterised by high activity of fatty acid oxidation (FAO), tricarboxylic acid (TCA) cycle and oxidative phosphorilation (OXPHOS). Activated naïve T-cells in the cell division phase highly express glycolysis, glutaminolysis and fatty acid synthesis (FAS). Effector T-cells have low TCA and OXPHOS activity, and high glycolysis and glutaminolysis. T𝐻17 cells are further characterised by high levels of FAS. Regulatory T-cells have low glycolysis but high OXPHOS and FAO. Finally, memory T-cells are characterised by low glycolysis, high OXPHOS, and the futile cycle between FAS and FAO. Activated T cells also increase methionine uptake. This amino acid is highly important for proliferating cells, as it is the predominant ’start’ amino acid in protein synthesis. The methionine cycle further generates methyl donors required for RNA and histone methylation[35]. Finally, activated T-cells go through asymmetric cell division. Associated with different distribution of metabolic mediators, asymmetric division drives activated T-cells toward either an effector or memory phenotype [36]. Those cells that inherit higher levels of amino acids go through a regulatory process that makes them more glycolytic and more likely to develop into effector cells [37]. They exhibit increased expression of effector molecules and increased expression of the large neutral amino acid transporter CD98, critical for clonal expansion and effector cell differentiation [38]. The other cells have enhanced lipid metabolism, spare respiratory capacity (SRC, the extra capacity that cells have available to produce energy in response to increased stress or work) and survival. All features of the memory phenotype, these cells are primed to become long-lived memory cells [37]. 12 2.3. METABOLISM IN T-CELLS AND CANCER Effector T-cells Activated effector T-cells use glycolysis and glutaminolysis for ATP generation and redox balance. In fact, these cells upregulate Glut1, which mediates the increased glucose uptake [39]. Cells that starve of glucose have increased amounts of glutamine-derived glutamate and pyruvate to allow them maintain TCA activity (figure 1)[40]. This metabolite allocation allows lipids and other amino acids to be redirected to generate biomass, such as nucleotides for DNA and lipids for membranes, to support cell division and effector cell functions. Effector T cells have increased fission and more punctate mitochondria with looser cristae, which leads to a physical dissociation of ETC supercomplexes that may cause electrons to linger in the complexes and imbalance redox reactions. This imbalance can cause an increase in NADH levels that slow the TCA cycle. To restore the redox balance, cells augment glycolysis and shunt pyruvate as excreted lactate, known as aerobic glycolysis. This allows regeneration of NAD+ from cytosolic NADH [41]. Memory T-cells Besides the initial asymmetric division during T-cell activation, some effector CD4+ and CD8+ T-cells are capable of differentiating into long-lived quiescent memory cells after pathogen clearance instead of suffering apoptosis. In this phase, T-cells no longer undergo rapid growth that requires high rates of biosynthesis. Instead, they require efficient energy generation to support basic cellular functions and prevent cell death [30]. Thus, aerobic glycolysis is reduced, while FAO and mitochondrial metabolism is favoured (figure 1). This also allows memory T-cells to increase their capacity to undergo oxidative metabolism under metabolic stress, due to higher SRC [29, 30]. Fatty acids are synthesized by memory T-cells from glucose in internal lysossomal stores, rather than acquiring them from an extracellular source. The fatty acids formed are then broken down by lisossomalacid-lipase-mediated lipolysis to liberate free fatty acids from storage to be used as substrates in FAO. This futile cycle may be engaged to ensure a continuous lipid supply for FAO, regardless of the extracellular lipid content, and maintain enzyme expression to keep cells primed and ready for rapid recall in the event of pathogen re-encounter [42]. Although SRC and this futile cycle are important for the function of central memory T-cells in the lymph nodes, effector memory T-cells that reside in tissues rely on the import of extracellular fatty acids and on glycolysis [43]. Nevertheless, and although not yet clear, the substrates used for FAO in different memory T-cell populations may be due to substrate availability in different tissue environments [29]. Memory T-cells have increased expression of phosphoenolpyruvate carboxykinase (PCK1), which mediates glycogen biosynthesis and subsequent glutathione production through the pentose phosphate pathway, which maintains memory T-cells by reducing ROS levels [44]. Effector 𝐶𝐷4+T-cells subsets The distinct metabolic pathways engaged by the different 𝐶𝐷4+T-cell subsets (figure 1) are not only crucial for their differentiation and survival, but also support important cell-specific functions. Although all effector 𝐶𝐷4+T-cell subsets have higher rates of glycolysis than naïve T-cells, T𝐻1,T 𝐻2and T𝐻17 cells have higher levels of Glut1 compared with regulatory T-cells [45]. Increased glucose uptake is not only sufficient to selectively enhance effector T function, but the inhibition of glucose metabolism 13 CHAPTER 2. BACKGROUND is capable of selectively inhibit effector T-cells, specially T𝐻17 [46]. T𝐻17 cells display lower OXPHOS than other helper T-cells. Regulatory T-cells are only partially dependent on glycolysis, resulting from the simultaneous increase in OXPHOS [47]. While high levels of glycolysis induces their proliferation, it limits their suppressive capacity [48]. Once glycolysis-related genes are inhibited and lipid and oxidative metabolism genes are promoted, regulatory T-cells reach their maximal suppressive ability. T𝐻17 cells rely on de novo FAS, rather than acquisition of extracellular fatty acids, to meet lipid requirements [49]. Cholesterol biosynthesis is required for suppressive function in regulatory T-cells [50]. In spite of exhibiting high FAO levels, Cluxton D. et al [47] observed that regulatory T-cells do not entirely rely on FAO. Glycolysis can serve as an alternative energy source. Raud B. et al [51] even found that deletion of the mitochondrial long-chain fatty acid transporter CPTI, essential for long-chain FAO, did not affect regular T-cell development and function. 2.3.2 Tumour cells Cancer cells go through metabolic changes relative to normal cells that support the acquisition and maintenance of malignant properties. Many of these changes are observed among most types of cancer cells, which supports reprogrammed metabolism as a hallmark of cancer [52]. Nutrients Uptakes The uptake of nutrients from the environment must be increased to fulfil the biosynthetic demands associated with the high proliferation rate of cancer cells. Glucose and glutamine are the two most important components. Their catabolism provides various carbon intermediates that are used for the assembly of various macromolecules. Furthermore, controlled oxidation of carbon skeletons of glucose and glutamine allows generation of NADH, FADH2and NADPH. Glutamine further contributes with reduced nitrogen for the de novo biosynthesis of nitrogen-containing compounds [52]. Normal cells do not import nutrients in a constitutive manner, as this process is strictly regulated by growth factor signalling and interactions with the extracellular matrix. However, accumulated oncogenic mutations in cancer cells allows tumour cells to constantly scavenge for glucose, glutamine and essential amino acids from the extracellular environment [18]. The glucose transporter Glut1 [53] is up-regulated, just like the glutamine transporters ASCT2 and SN2. Glutaminase expression is also promoted [54], which converts glutamine into glutamate, whose accumulation promotes TCA cycle. Increased cysteine uptake is promoted by the glutamate accumulation, as the xCT transporter imports cysteine in exchange of glutamate [55]. Cysteine, a sulfur containing amino acid, can be involved in the biosynthesis of glutathione, iron-sulfur clusters, and hydrogen sulfide (𝐻2𝑆). 𝐻2𝑆is involved in the protection from oxidative stress, increased mitochondrial respiration, protection from apoptosis, and facilitation of angiogenesis [56]. Although the majority of normal proliferating cells require exogenous supply of glutamine despite the existence of a glutamine biosynthetic pathway, some cancer cells over-express glutamine synthetase and are able to produce glutamine de novo [57]. 14 2.3. METABOLISM IN T-CELLS AND CANCER Tumour cells are also able to use opportunistic modes of nutrient acquisition to access normally inaccessible nutrient sources, as well as to recover pre-made molecules when their synthesis within the cell is compromised[52]. The deficit of unsaturated fatty acids in cancer cells can be overcome by the import of ready-made fatty acids. Cancer cells can even induce the neighbouring normal cells to release stored lipids [58]. Bioenergetics Cancer cells exhibit aerobic glycolysis, a robust provider of precursors and reducing equivalents necessary for biosynthesis of macromolecules essential for cell proliferation [52]. The first metabolite produced in the glycolysis pathway, glucose-6-phosphate, enters the pentose phosphate pathway (PPP), generating NADPH and ribose-5-phosphate, a structural component of nucleotides. In fact, PPP utilisation is frequently elevated in tumour cells [59]. 3-Phosphoglycerate can be used as a precursor for the biosynthesis of serine and glycine, and as means to generate methyl donor groups and NADPH. Serine, for example, is a major substrate of the folate cycle, an essential source of precursors for the biosynthesis of purines and thymidine. Methylene tetrahydrofolate dehydrogenase 2 (MTHFD2), a component of this cycle, has been found to be one of the most frequently over-expressed metabolic enzymes in cancer [60]. Furthermore, the final reaction of glycolysis is catalysed by pyruvate kinase (PK) in the form PKM2 in most tissues, including tumours [61]. PKM2 is activated by serine [62]. Nevertheless, most cancer cells still generate the majority of ATP through mitochondrial function, despite the high glycolytic rates [52]. There is no actual shift between TCA and glycolysis, like initially proposed by Warburg et al [2], but rather a considerable reduction of the TCA cycle activity to a state sufficient to maintain mitochondrial integrity and ATP production. In addition to pyruvate derived from glycolysis, fatty acids and amino acids can supply substrates to the TCA cycle. In fact, glutamine can provide acetyl-CoA as a precursor when pyruvate oxidation to acetylCoA is compromised by hipoxia or ETC impairment. Also, most proliferating cells depend on a continuous supply of glutamine to maintain TCA cycle intermediates [52]. Biosynthesis of macromolecules The production of biosynthetic intermediates through metabolic pathways such as glycolysis, PPP, TCA cycle and non-essential amino acid synthesis allows the assembly of larger and more complex molecules, required for replicative cell division and tumour growth. Among these, the most commonly studied in cancer metabolism are proteins, lipids and nucleic acids [52]. While glutamine-derived glutamate works as a nitrogen donor for the production of several non-essential amino acids via transamination, the amide nitrogen of glutamine is used by asparagine synthetase (ASNS) to produce asparagine from aspartate. Notably, ASNS is frequently up-regulated in tumours and is associated with poor prognosis [63, 64]. Essential amino acids, in turn, are acquired from the extracellular space through surface transporters under the influence of growth factor signalling [65]. 15 CHAPTER 2. BACKGROUND Arginine is a non-essential amino acid that can become conditionally essential in some tumourigenic contexts. For example, arginino succinate synthase (ASS1), essential to the de novo biosynthesis of arginine, is frequently epigenetically silenced in pancreatic and renal cancers [66, 67]. Inactivation of this enzyme causes cancer cells to accumulate ornithine, which is then used in the production of polyamines. These compounds have been shown to inhibit apoptosis and promote tumour growth invasion [52]. Proline can be produced from glutamate or from arginine-derived ornithine. Notably, the principal enzyme in proline production, pyrroline-5-carboxylate reductase (PYCR1), is one of the most commonly overexpressed enzymes in tumours [60]. The suppression of ASS1-driven argininosuccinate production can also cause accumulation of its substrate, aspartate, required for nucleotide production [52]. Purine and pyrimidine nucleotides are required for synthesis of RNA and DNA. The expression of phosphoribosyl pyrophosphate synthetase 2 (PRPS2) [68] and carbamoyl phosphate synthetase II (CAD) [69] are up-regulated by c-Myc in tumour cells. These enzymes are involved in purine and pyrimidine biosynthesis, respectively. The capacity to rapidly produce lipids in cancer cells facilitates the formation of membranes, the alteration of membrane composition in favour of oxidative damage-resistant saturated fatty acids, lipidation reactions, and cellular signalling. The activity of several enzymes involved in lipid synthesis, even in lipidreplete conditions, are up-regulated in cancer cells. 2.3.3 Tumour Micro-Environment The most frequently found tumour infiltrating lymphocytes (TILs) within the TME are tumour-associated macrophages (TAMs) and T-cells [28]. TAMs are alternatively activated macrophages reprogrammed to display various tumour-promoting functions [28, 70]. Polymorphonuclear leukocytes are rarely seen in human TMEs [71]. High numbers of cytotoxic T cells and T𝐻1cells are correlated with better survival in some cancers, including invasive colon cancer, melanoma, multiple myeloma, and pancreatic cancer [28]. For example, T𝐻1cells maximize the killing efficiency of macrophages and proliferation of CD8+T-cells [72]. However, during de novo carcinogenesis, anti-tumour T-cells cannot control tumour growth, due to tumour-induced tolerance mechanisms in most cancers [73]. In fact, tumour cells have the ability to actively downregulate all phases of anti-tumour immune responses through metabolism, affecting the recruitment and function of immune cells. The different constituents of the TME can also interact closely with each other to control and shape tumour cell survival, invasiveness and metastatic dissemination. These interactions dictate the ability of the immune system to fight cancer and even responsiveness to therapeutics. This dynamic relationship is set early on in malignant growth and evolves throughout the life history of a tumour [27, 28]. Figure 2 gives an overview of some of the metabolic interactions within the TME, discussed below. When tumour cells are exposed to hypoxic conditions, the production of the hypoxia-inducible factor 1𝛼 (HIF-1𝛼) is up-regulated, which promotes glycolysis and leads to activation of angiogenesis-promoting 16 2.3. METABOLISM IN T-CELLS AND CANCER factors. In the presence of oxygen, HIF-1𝛼is degraded. The acidic microenvironment caused by tumours, that up-regulates glycolysis and increases production of lactic acid, can ’simulate’ the effects of hypoxia, even in the presence of oxygen. Thus, HIF-1𝛼may not be suppressed even in normoxia conditions [72]. Figure 2: Summary of some of the metabolic interactions within the tumour micro-environment (TME). The up-regulation of glycolysis in tumour cells increases the secretion of lactic acid into the micro-environment, causing an acidic microenvironment that ’simulates’ the effects of hypoxia. This leads to an up-regulation of HIF-1𝛼, activating angiogenesis-promoting factors. These factors can stimulate the proliferation of MDSCs, which liberate arginase into the micro-environment that will consume L-arginine. Shortage of this metabolite can inhibit CD8+T-cell function. Constitutively expressed in most human tumours, IDO is involved in the catabolism of tryptophan, which induces immunosuppression through T-cell anergy and depletion. | IDO: Indoleamine-2,3-dioxygenase; HIF-1𝛼: hypoxia-inducible factor 1𝛼; MDSCs: myeloidderived suppressor cells; NO: nitric oxide; ROS: reactive oxygen species Angiogenic factors also stimulate the proliferation of myeloid-derived suppressor cells (MDSCs), which include immature dendritic cells, neutrophils, monocytes, and early myeloid progenitors. When stimulated, these cells up-regulate and liberate arginase into the micro-environment, which consumes L-arginine. Shortage of L-arginine in the environment can inhibit CD8+T-cell function. Stimulated MDSCs also increase production of nitric oxide (NO) and ROS. NO can suppress T-cell function through inhibition of MHC-II expression and T-cell proliferation and apoptosis [72]. TAMs are capable of blocking CD8+T-cells proliferation or infiltration by releasing factors with immunosuppressive potential like ROS [70]. TAMs can further suppress surface proteins on infiltrating T-cells through nitrosylation, inhibiting T-cells’ anti-tumour functions [74]. ROS released by tumour cells can induce cancer-associated fibroblasts to up-regulate aerobic glycolysis and secrete lactate and pyruvate. Tumour cells then consume these two metabolites [75]. Alternatively, the lactate secreted by cancer cells is taken up by cancer-associated fibroblasts and used as fuel to drive tumour-promoting functional activities [76]. Cancer-associated fibroblasts are the most predominant nonhematopoietic stromal cell type in the TME [27]. 17 CHAPTER 2. BACKGROUND Indoleamine-2,3-dioxygenase (IDO) is constitutively expressed in most human tumours [77]. IDO is involved in the catabolism of tryptophan, an essential amino acid for T-cell proliferation and differentiation [71]. The catabolism of tryptophan into kynurenine by this enzyme induces immunosuppression through T-cell anergy and depletion [77]. Plasmacytoid dendritic cells have shown expression of indoleamine-2,3dioxygenase (IDO), as well as defective production of type I interferon [74]. Presence of TNF-𝛼and IFN-𝛾can dramatically increase MSCs’ expression of inducible nitric oxide synthase (iNOS) [78] and production of IDO [71]. iNOS consumes arginase, producing ROS and NO. Mevalonate pathway intermediates produced by tumour cells were shown to activate and promote 𝛾𝛿T-cells’ anti-tumour responses [15]. 2.3.4 Targeting metabolism for therapy The most common therapies applied to cancer patients are surgery, radiation therapy and/or chemotherapy. Surgery is mostly used for non-invasive solid tumours and coupled with other treatments. While radiation therapy uses high doses of radiation to kill or slow the growth of tumour cells, by damaging their DNA beyond repair, chemotherapy uses drugs to this end. However, these two treatments do not only kill tumour cells, they also affect healthy cells. For example, the most common side effects comprise fatigue, hair loss, nausea and vomiting. Through more recent years, other alternatives for cancer therapy have been studied to generate as few harmfull side effects as possible to the patient. These can be achieved by searching for specific targets in tumours cells that do not affect healthy cells, or even specific ways of enhancing the immune response against cancer. As cancer cells go through metabolic changes relative to normal cells that support the acquisition and maintenance of malignant properties, the metabolism of cancer cells has been studied for immunotherapy. Naturally, inhibition of glycolytic and glutaminolytic enzymes has been extensively studied. Diclofenac has been reported to reduce tumour growth, the quantity of regulatory T-cells and lactate in the microenvironment in a glioma model [79]. Neutralisation of the TME’s acidic environment with bicabornate or esomeprazole, for example, improves cytotoxic T-cell and NK cell anti-cancer immune responses [80]. Hexokinase (HK), a glycolytic protein, is overexpressed in many tumour cells and its inhibition was shown to delay tumour progression in pre-clinical mouse models [81]. However, 1-deoxyglucose (2DG), an inhibitor of HK, also leads to impairment of T-cells’ metabolism [82]. Bis-2-(5-phenylacetamido-1,2,4-thiadiazol-2-yl) ethyl sulfide (BPTES) is a glutaminase inhibitor that showed anti-cancer immunity in several tumour models with elevated activity [83]. OXPHOS is also a great potential target to eliminate tumour cells. The anti-diabetic drug metformin can act as an anti-cancer agent by inhibiting the complex I from ETC. This causes a decrease in ATP levels that lead to cancer cell death [80]. However, this drug’s uptake occurs through the organic cation transporters (OCTs), only present in a few tissues, such as liver and 18 2.4. CONSTRAINT-BASED MODELING AND HUMAN CANCER kidney, and in certain tumour cells [65]. Regarding effects on the immune system, this drug can further enhance memory T-cells and regulatory T-cell expansion [80]. Dichloroacetate (DCA) induces a shift from glycolysis to OXPHOS, thus inhibiting tumour cells growth in vitro and in mouse models. However, it also affects T-cells, favouring regulatory T-cell formation [84]. Down-modulation of IDO has been shown to improve anti-tumour responses [74, 80]. Imatinib, for instance, activates effector T-cells and suppresses regulatory T-cells in an IDO-dependet manner [85]. Inhibition of the rate-limiting enzyme in FAO, CPTI, has anticancer effects in vitro and in vivo . However, etomoxir showed hepatotoxicity in patients with congestive heart failure, and other inhibitors are still to be approved for cancer therapy [86]. Resistance to cancer therapies may result from not taking into consideration the TME, as briefly noted above in few examples. Furthermore, CSCs are typically therapy-resistant due to decreased oxidative stress response, increased genomic stability, and expression of multiple drug resistance transporters [78]. TME cells are not subject to mutational and epigenetic changes that result in drug resistance. Thus, targeting the TME cells, specially immune cells, along with tumour cells can be advantageous. Regarding tumour stroma, targeting the tumour extracellular matrix can boost natural anti-tumour immunity and improve immune-therapeutics efficacy. However, it may also enhance regulatory T-cell infiltration and increase angiogenesis [27]. Despite the progresses, clinical responses may be transitory and have limited benefits in long term [74], mostly due to drug resistance caused by the existence of similar pathways that are alternatively upregulated by the cell. Furthermore, investigation of a potential metabolic target is normally only performed in tumour cells, without counting with the possible negative side effects on other cells in the TME and outside. These problems show the value in creating in silico metabolic models of the whole metabolome of immune and tumour cells to study the effects of metabolic targets and off-targets on each individual cell, as well as interactions between the different cells in the TME upon disruption with approved drugs. This could be further fine-tuned to patient-specific cases, personalising each patient’s therapy. Metabolomics can aid in such a way that cancer therapy and cancer immunotherapy act specifically on malignant cells, remarkably reducing the side effects to the patient. 2.4 Constraint-based Modeling and Human Cancer Two popular, but very different, approaches to model a cell’s metabolism are the kinetic (dynamic) and stoichiometric modeling. The kinetic modeling, as the name suggests, relies on the enzyme kinetics information to model the metabolite concentrations and reaction fluxes through time [87, 88]. However, kinetic models require a lot of details for their construction, which must be obtained through experiments that are difficult to perform. Because of this, they usually end up covering only a few pathways or models from small organisms [88, 89]. 19 CHAPTER 2. BACKGROUND Pruning methods, in turn, start with a set of core reactions, known to be present in the desired tissue by going through literature or experimental data, and remove the other reactions of the generic GSMM whose removal does not cause loss of reaction functionality in the core set. To achieve this, these methods establish a trade-off between maintaining the model as concise as possible and including all core reactions in the final model, allowing a core reaction to be removed if it requires too many undesirable reactions to be active [121]. The main advantages of this last type of methods relies on making it possible for the user to define the set of core reactions using multiple different sources and know that reactions with high evidence of being present in a tissue are always included in the final model, besides generating a flexible and functional metabolic model. Examples of such algorithms are MBA [122], mCADRE [123], fastCORE [124], and CORDA [121]. However, pruning algorithms are not deprived of disadvantages. The order in which reactions are removed can affect the outcome of the final model, as well as the possible removal of fundamental reactions to achieve a concise tissue-specific reconstruction can cause physiologically unlikely flux distributions [121]. The fastCORE [124] algorithm, for instance, takes a set of reactions that have strong evidence to be active in the context of interest and searches for a subnetwork that contains all reactions from the core set and a minimal set of additional reactions, by assessing the flux consistency through FVA. This flux consistency is characterized by each reaction in the subnetwork being active in at least one feasible flux distribution. Several human cancer metabolic models have been constructed in the past years. Among these are models of glioblastoma [125], ovarian cancer [126], hepatocellular carcinoma [126], melanoma [127], breast, urothelial, lung and renal cancers [128], and colorectal cancer [129]. Several other studies have focused on constructing models specific for cancer cell lines [130–132], or even from patient-specific samples [119, 133–136]. Regarding the study of the immune system, tissue-specific reconstruction algorithms have also been used to generate metabolic models of immune cells, including CD4 T [113], naive CD4 T [137–139], T𝐻1 [138, 139], T𝐻2[138, 139] and T𝐻17 [138] cells, as well as macrophages [140, 141], monocytes [113], B cells [113] and NK cells [113]. However, none were applied in the context of cancer. 2.4.3 Omics Data The advent of omics technologies have revolutionised the way biological research is conducted. They allow the analysis of a global set of molecules and their interactions simultaneously. Integrating data from different omics helps creating an overall ’snapshot’ of a cell’s metabolism and evaluate how different mutations affect metabolism. Omics data can be used to find the reactions that should be present in a cell-type or tissue in order to reconstruct a context-specific metabolic model from a generic model. In fact, developing large-scale 26 2.4. CONSTRAINT-BASED MODELING AND HUMAN CANCER metabolic models became possible due to the advance of omics technologies. Genomics It consists on the study of the genomic material of an organism. In cancer patients, genomic tests are often performed before treatment. Genetic mutations are identified by comparing the tumour genome with the patient’s normal tissue or a reference genome at a single base pair resolution [142]. The most used techniques have been Sanger sequencing, microarrays, and Next-Generation Sequencing (NGS) [143]. Sanger sequencing consists on a base-by-base sequencing, but it can only capture up to one thousand bases per run. Microarrays, in turn, involve the binding of the different cDNA sequences obtained in a solution to the respective probe in the array, allowing the measurement of the relative concentrations of those sequences [144]. NGS enables the identification and quantification of transcripts without prior knowledge of a particular gene sequence, and can provide information regarding alternative splicing and sequence variation. Although NGS allows whole genome sequencing (WGS), allowing identification of all coding and non-coding variants, it is also possible to only screen variants in the coding region (whole exome sequencing, WES) [143]. The great technical advances achieved throughout the years to obtain genomic data culminated in the sequence of the whole human genome, the Human Genome Project, in 2003. Since then, Sanger sequencing and microarrays have been completely substituted by NGS in sequencing the human cells’ genome. Transcriptomics This technology studies the transcriptome, which represents all RNA transcripts in a cell, both coding (mRNA) and non-coding (ribosomal, transfer, etc) RNAs. Normally, mRNAs are the main focus, as the quantification of each transcript provides and insight into the gene expression levels. This allows a better understanding of the dynamics of metabolism [143]. To assess transcriptomics data, similar techniques to those in genomics are used, namely microarrays and RNA-sequencing (RNA-seq). While microarrays only allow measurement of relative concentrations of the different transcripts, quantification in RNA-seq is done by counting the number of sequence reads assigned to different transcripts [145]. More recently, single-cell techniques like single-cell RNAseq (scRNAseq) have been developed. This technique allows quantitative characterisation of each cell’s transcriptome in a sample, giving valuable information about cell-types and states at high resolution [146]. Especially in cancer, where samples from bulk transcriptomics are often analysed as if they were an homogenous population of tumour cells, not accounting for the diversity of stromal and immune cells present that can be crucial in the understanding of cancer. Proteomics This omics studies the entire set of proteins in a cell, at a precise developmental or cellular phase. Techniques used in proteomics, mainly mass spectrometry (MS), are less scalable than those used to study nucleic acids, like NGS. Still, considering that protein quantities and activities of malignant cells 27 CHAPTER 2. BACKGROUND are affected by distinct replication and metabolic processes, proteomics has the potential to uniquely characterise the malignant cell or tissue and to discover diagnostic markers of a disease [142, 143]. Although MS allows identification and quantification of proteins in a sample, the eventual function of an existing protein in a sample depends on how the protein is folded, i.e., its 3D structure. To derive the proteins’ structures in a sample, techniques such as X-Ray, Nuclear Magnetic Resonance (NMR) and cryo-electron microscopy come into play. These techniques allow visualization of protein domains, deduce protein function, study structural changes following disease associated mutations, and discover and develop drugs [143]. Proteomics can also be used to study protein-protein interactions (PPIs), which consist on either two proteins that physically interact in a complex, or simply share the same location, as interacting proteins are likely to share common tasks or functions [143]. Because the proteome is extremely dynamic, there is an elevated sample heterogeneity, even with the same types of samples and conditions, that complicates the development of a universal and comprehensive human proteome reference and the comparison of different studies [143]. Epigenomics Epigenomic changes include DNA methylation and chromatin modifications like histone acetylation, methylation, phosphorylation and others. These heritable changes have a great impact in the expression patterns of genes, and can be affected by mutations in enzymes in charge of DNA methylation, demethylation and chromatin modification. This may allow tumour progression if the suppressed or activated gene hinders or promotes tumour progression. Epigenomic changes are often observed in many cancers. Identification of DNA methylation status is performed by using bisulfite treatment, which only modifies unmethylated cytosines to uracils, followed by sequencing and comparison between bisulfite-treated and untreated samples [142]. Regarding histone modifications, these are identified by high-throughput DNA sequencing technologies coupled with chromatin immunoprecipitation (ChIP), where modification-specific antibodies immunoisolate DNA-histone complexes with the desired histone modifications. The selected DNA sequences are then identified via microarrays (ChIP-chip) or sequencing (ChIP-seq) [142]. Metabolomics This technique comprises the analysis of the metabolites produced during biochemical reactions, which depend on gene expression to be produced. This omics is thus closely related to the ones described above, as the production of metabolites can reflect particular combinations of individuals’ genetics and environmental exposures. Metabolomics can be used to develop diagnostics and understand relevant molecular pathways under specific conditions [143]. The most used techniques in this field are Mass Spectrometry, coupled with liquid or gas chromatography (LC/GC-MS), and Nuclear Magnetic Resonance (NMR). These methods, despite being able to give important information, are far less common than those used in the previous omics, and there is not a single one that allows the analysis of the whole metabolome. 28 2.4. CONSTRAINT-BASED MODELING AND HUMAN CANCER 2.4.4 Applications of Modeling in Human Cancer At first, the reconstructed models are used to simulate cellular responses under certain conditions and compare the predictions to experimental data or expected behaviour from literature. This model evaluation helps to further refine the model and assess its accuracy. As cancer cells are expected to have reduced TCA cycle activity and increase glycolysis, even in aerobic conditions, many authors start out by evaluating if the cancer models show reduced fluxes in the reactions that are part of the TCA cycle and increased fluxes in those part of glycolysis, in aerobic conditions [125–127, 132], as opposed to normal cell models [126, 127]. Other functions are tested to assess if the models simulate what has been experimentally found. Özcan E. et al [125] glioblastoma metabolic models showed active flux for the pyruvate dehydrogenase reaction, where glucose is metabolised through pyruvate dehydrogenase rather than pyruvate carboxylase in glioblastoma cells. A metastatic melanoma model [127], when compared to the primary melanoma model, showed increased fluxes in purine and pyrimidine synthesis and glutamate, arginine and proline metabolism, while reaction fluxes in coenzyme-A metabolism, fatty acid synthesis and OXPHOS decreased. As for T-cells, Puniya et al [138] showed that T𝐻1,T 𝐻2and T𝐻17 models had more flux through fatty acid biosynthesis and less flux through fatty acid 𝛽oxidation than the naïve CD4 T-cell model. Limiting glucose from the environment resulted in decreased growth rate in all the 4 models. However, there was not a significant effect in growth rate when glutamine was removed from the medium of T𝐻1,T 𝐻2and T𝐻17 models, and growth rate of naïve CD4 T-cell model was more dependent on glucose and glutamine uptake than the other models. Still, gene deletion analysis revealed that more than 70% of the gene essentiality predictions agreed with previous experimental tests. There are several methodologies that make use of a GSMM to gain insights on the capabilities of a cell’s metabolism. This could be of great value in better understanding the underlying metabolic mechanisms of cancer and the tumour micro-environment. Metabolite Biomarkers A challenge in cancer diagnosis is the identification of metabolite biomarkers that are present in biofluids such as plasma, urine and feces. This allows measurement of biomarkers in a non-invasive, cost-effective way for early diagnosis and monitoring treatment efficiency [147]. These biomarkers can distinguish not only between healthy and disease cases, but also between clinical groups such as subtypes of cancer. The study by Nam et al [148] is an example of potential biomarkers that were found using metabolic models. While succinate and fumarate were shown to be biomarkers for gastric cancer, as seen in literature, due to loss of function in sub-unit B of succinate dehydrogenase (SDHB), palmitate, D-glucose, and adenosine are potential biomarkers of leukemia due to fatty acid synthase’s (FASN) loss of function. Pathway Analysis Pathway analysis allows to understand different important aspects of the metabolic network. The more a reaction appears in different EFMs, the more important it is for the structure of the network [149]; when different sets of connected reactions lead to the same output from the same input, 29 CHAPTER 2. BACKGROUND the network has pathway redundancy. A high degree of redundancy might translate into the network being better at tolerating knockouts or drug inhibition than one of low redundancy, thus showing high robustness [91]. Minimal cut sets (MCSs) [90] represent a set of potential drug targets to prevent the functioning of an objective reaction. Drug Targets As mentioned before, the metabolism of cancerous cells is modified during tumour development to meet requirements of fast cellular proliferation or due to genetic alterations. There might be certain reactions whose enzymes are essential for a tumour cell’s viability but not for healthy cells. These enzymes can become interesting drug-targets. Also, as metabolism is evolutionarily more conserved than other biological processes, it is less probable for cancer cells to evolve resistance to these drugs by developing alternative pathways [150]. While the discovery of new drugs presents itself as a challenging task that requires a very long period of research and development before any new compound can be commercialised, exploiting the properties of already available drugs, whose information about the therapeutic and toxicity effects are already known, frequently becomes the focus when using metabolic models [121,127, 131, 133, 151]. The identification of drug targets comes from studying the essentiality of reactions and metabolites. Ghaffari et al [131] tested cancer cell line models for metabolite essentiality and, after testing the essential metabolites in normal tissue models, 85 metabolites were found to be essential only in the cancer models. One of these metabolites was L-carnitine. Knowing that perhexiline malate salt inhibits CPT1 (responsible for the translocation of conjugated L-carnitine and long chain fatty acids from cytosol to the mitochondria), and partly CPT2, the authors treated two different cancer cell lines with this compound. Significant reduction in both cell lines viability was observed. Yizhak et al [132] predicted 17 metabolic enzymes as targets to mitigate cancer cell migration, as their individual knockout from the models reduced the glycolytic to oxidative ATP flux ratio, positively associated with this event. Indeed, most have been found to have significantly higher expression levels in metastatic breast cancer patients than those non-metastatic and lower expression of 9 of those enzymes are associated wit improved long-term survival. Immune cell models were used to find potential drug targets for T cells in rheumatoid arthritis, multiple sclerosis or primary biliary cholangitis [138], and to study metabolic changes during infection of macrophages by a pathogen [141]. However, none were applied in the context of cancer. 30 j A Colorectal Cancer Atlas of scRNAseq Data The type of cancer focused throughout this work was colorectal cancer (CRC). Also known as colorectal adenocarcinoma, this cancer usually emerges from the glandular, epithelial cells of the colon or rectum [152, 153]. Colon or rectal cancers are often merged due to the many biological and clinical common features [152]. Colorectal cancer is often asymptomatic until it substantially grows and spreads, hindering the prognosis and subsequent survival. In fact, the 5-year survival rate is close to 90% when it is diagnosed at an early stage, but of only 13% when the diagnosis is delayed. The limited number of tests that can be used for timely and efficient screening or diagnosis also affects the diagnosis, with up to 90% of the cases being diagnosed after symptoms onset [152]. According to the statistics provided by the International Agency for Research and Cancer (IARC) of the World Health Organization (WHO) in 2020 [154], colorectal cancer is the third most frequent malignant disease around the world, comprising 10% of total malignancies. The number of deaths in 2020 for colorectal cancer was approximately 935 000, representing 9.4% of cancer-related deaths, only preceded by lung cancer. As CRC is a very heterogenous malignancy with different pathological and genetic signatures, it is very difficult to find one single molecular therapy to treat this type of cancer. In fact, surgery remains the primary course of treatment for early diagnosis cases, while in advanced cases cytotoxic therapies are met with rapid evolution of drug resistance and cancer recurrence [153]. However, therapeutic stratification based on the pathological and gene signatures may ultimately lead to improved therapeutic outcomes. The most widely used stratification approach is the Consensus Molecular Subtypes (CMS) of CRC [155]. Briefly, there are 4 CMS subtypes and an additional mixed phenotype when a clear CMS subtype cannot be assigned [156]. CMS1 is characterised by high microsatellite instability (MSI), CpG island methylator phenotype (CIMP+), hypermutation, frequent mutations of BRAF gene, and immune infiltration and activation. CMS2 is known for high somatic copy number alteration (SCNA) and activation of the WNT and MYC signalling pathways. CMS3, in turn, has frequent KRAS mutations and is characterised by metabolic deregulation. Finally, CMS4 is characterised by high SCNA, stromal infiltration, TGF-𝛽activation, and angiogenesis. 31 CHAPTER 3. A COLORECTAL CANCER ATLAS OF SCRNASEQ DATA To construct a wide number of T-cell genome-scale metabolic models that characterised different subtypes of T-cells, we wanted to make use of scRNAseq data, as this type of data allows quantitative characterisation of each cell’s transcriptome in a sample. This lets us construct models of different T-cell types from the same tumour micro-environment. For this, we first constructed a colorectal cancer (CRC) atlas of scRNAseq data from different studies. The following tools were used to read, process, analyse, and save scRNAseq data: R packages Seurat (v.4.0.3) [157], SeuratDisk (v.0.0.0.9019) [158], SeuratObject (v.4.0.2) [159], clustree (v. 0.4.3) [160], SIGMA (v.0.0.0.1) [161], infercnv (v.1.8.0) [162], CMScaller (v.0.99.2) [163]. The scripts were run with the R version 4.1.0. More detailed information, including scripts, is available in the GitHub project ?iiTb, ff;Bi?m#X+QKfb`+`/QbQf*_*nhGa. 3.1 Datasets Collected Raw counts from four publicly available datasets, named CRC_Qian [164], GSE132465 [165], GSE144735 [165] and Colon_smillie [166], were used to construct the atlas. Three correspond to CRC studies, while one ( Colon_smillie ) has data from colon of patients with ulcerative colitis and healthy individuals. As such, the samples related to the healthy individuals were used. Table 1 shows an overview of the datasets collected. Table 1: Overview of the datasets used to construct the CRC atlas CRC_Qian GSE132465 GSE144735 Colon_smillie Nº Patients 7 23 6 12 Nº Samples 21 33 18 48 Nº Cells (before QC) 44 684 63 689 27 414 51 705 Nº Cells (after QC) 29 793 57 804 24 510 50 811 Technology 3’ 10xGenomics 3’ 10xGenomics 3’ 10xGenomics 3’ 10xGenomics Studies CRC_Qian and GSE144735 have tumour (border and core) and normal matched mucosa samples for each patient. For study GSE132465 , 10 of the 23 patients have tumour and normal matched mucosa samples, while the other 13 only have tumour samples. Regarding Colon_smillie dataset, each donor has a sample extracted from two different locations of the colon. 3.2 Quality Control Before merging the datasets together, they were individually checked for quality control. This allowed us to filter low quality cells and lowly expressed genes. Cells were kept if they followed all the following criteria: (1) At least 1000 UMIs in each cell; (2) Between, including, 250 and 6 000 of detected genes (i.e., genes with non-zero UMIs) in each cell; (3) A complexity metric (log10 𝑁𝑢𝑚𝑏𝑒𝑟_𝐺𝑒𝑛𝑒𝑠 log10 𝑁𝑢𝑚𝑏𝑒𝑟_𝑈𝑀𝐼𝑠) of at least 0.8; (4) Percentage 32 3.3. DATASET INTEGRATION of mitochondrial RNA not greater than 20%. Genes were considered lowly expressed, and thus removed, if they were not detected in 10 or more cells. This quality control was performed using the R language. At this point in the pipeline, datasets were already loaded into a Seurat object and stored as a h5Seurat format file, using the R packages Seurat (v.4.0.3) [157] and SeuratDisk (v.0.0.0.9019) [158]. 3.3 Dataset Integration The different datasets were integrated using Seurat (R package, v.4.0.3) [157], so that proper cell-type annotation could be performed. First, datasets were processed. Each dataset was normalised individually by dividing the counts from each cell by the total of the respective cell and multiplied by the scale factor 10000, followed by natural logarithmic transformation. Then, the top 2000 variables were calculated for each dataset and integration features calculated. Each dataset was then scaled with regression of the difference between the S and G2M phases (so that we could still distinguish proliferative cells from those not proliferating) and the percentage of mitochondrial RNA. PCA was finally performed to reduce dimensionality of datasets, in order to perform a faster and less computationally heavy integration. After this, the integration anchors were found and the data integrated. 3.4 Cell Annotation We started by separating the cells into 5 big groups of cell-types: Epithelial, Stromal, Myeloid, B-cells and T-cells. After this, we further annotated the cells for each of these groups individually, so that we could obtain more detailed and better annotations. To perform the initial separation of the cells into the mentioned groups, we used Seurat (R package, v.4.0.3) [157] to find the clusters at a resolution of 0.1. After inspecting the quality of the clustering by mapping into UMAP plot metrics, such as number of UMIs and genes per cell, cell-cycle scores for the S and G2M phases, and percentage of mitochondrial RNA, we explored the expression of gene markers to annotate the cells (figure 5). The following genes were used: EPCAM (epithelial cells); S100B , COL1A1 , VWF and ENG (stromal cells); CD14 , FCGR3A , FCER1A , GZMB , TPSAB1 and TPSB2 (myeloid cells); CD79A and MZB1 (B-cells); and CD3D (T-cells). For each group of cell-types, we found the clusters using a set of 10 different resolutions, ranging from 0.1 to 1 at a step of 0.1. After assessing the quality of the clustering, we found the resolution that best separates the cells into clusters with interesting biological information using clustree (R package, v.0.4.3) [160] and the expression of genes that were markers for cell-types that we wanted to find. The gene expression was evaluated by visualizing their expression mapped into UMAP and violin plots. Regarding Tand epithelial cells, more steps were performed to be able to fully annotate the cells correctly. These steps are mentioned next in their respective sub-sections. 33 CHAPTER 3. A COLORECTAL CANCER ATLAS OF SCRNASEQ DATA Figure 5: Annotation of the big groups of cells. Top-left UMAP is coloured by the annotated cell-types. All other UMAPs are coloured by the RNA expression of the genes used to identify the cells, which goes from blue (low expression) to yellow (high expression). Grey cells have no expression of the respective gene. 3.4.1 Stromal cells Regarding stromal cells, the following genes were used to annotate the cell-types (figure 6): VWF , PLVAP , CDH5 (endothelial cells); RGCC , RAMP3 (tip-like vascular endothelial cells); ACKR1 , SELP (stalk-like vascular endothelial cells); LYVE1 , PROX1 (lymphatic endothelial cells); S100B, PLP1 (enteric glia); SYNPO2 , CNN1 , PDGFRB (vascular smooth muscle-cells); RGS5 , ABCC9 , KCNJ8 (pericytes); TAGLN , ACTA2 , ACTG2 , MYH11 , MYLK , DES (myofibroblasts); COL1A1 , COL1A2 , COL6A1 , COL6A2 , COL3A1 , DCN (fibroblasts); THY1 , FAP (cancer associated fibroblasts). 34 3.4. CELL ANNOTATION Figure 6: Annotation of the stromal cells. Top-left UMAP is coloured by the annotated cell-types. All other UMAPs are coloured by the RNA expression of some of the genes used to identify the cells, which goes from blue (low expression) to yellow (high expression). Grey cells have no expression of the respective gene. 3.4.2 Myeloid cells Regarding myeloid cells, the following genes were used to annotate the cell-types (figure 7): CD14 , FCGR3A , MARCO , ITGAM (monocyte/ macrophage lineage); FCER1A , CST3 (conventional dendritic cells); IL3RA , SERPINF1 , GZMB , ITM2C (plasmacytoid dendritic cells); KIT , TPSB2 , TPSAB1 (mast cells). In our data, we were not able to separate monocytes from macrophages. However, we further separated the monocyte/ macrophage lineage into further groups: IL1B , IL6 , S100A8 , S100A9 (pro-inflammatory macro/mono); CD163 , SEPP1 , APOE , MAF (anti-inflammatory macro/mono); SPP1 (SPP1+ macro/mono); 35 CHAPTER 3. A COLORECTAL CANCER ATLAS OF SCRNASEQ DATA Finally, the cluster with𝛾𝛿T-cells and other non-conventional T-cells (3) (figure 15) contained𝛾𝛿T-cells ( TRDC , TRGC1 , TRGC2 ), NK cells ( KLRD1 , KLRB1 , XCL1 , XCL2 ), lymphoid tissue-inducer cells ( RORC , LTA , LTB ), NKT cells (NK-like signature genes and some T-cell markers). In this cluster, a sub-cluster of cells with high expression of heat-shock proteins ( HSPA6 , HSPA1B , HSPA1A , HSPB1 , HSP90AA1 , HSPD1 , HSPH1 ) was removed from the dataset. Figure 15: Cluster (3) . UMAPS with (A) Final annotation, (B) clusters before the final annotation, and (C) expression of relevant genes (expression goes from blue (low) to yellow (high), while grey cells is the absence of expression). 3.4.5 Epithelial cells First, we wanted to separate tumour from normal cells in tumour samples. For that, we used the R package infercnv (v.1.8.0) [162] to identify somatic large-scale chromosomal copy number variations (CNVs) by comparing the gene expression of the epithelial cells from tumour samples with the gene expression of the epithelial cells from normal matched samples. After processing the raw counts data accordingly, we used the six-state HMM-based CNV prediction method (i6 HMM) to predict, for each gene in each cell, the CNV level. There are 6 different levels: (1) complete loss; (2) loss of one copy; (3) neutral (i.e., no change); (4) addition of one copy; (5) addition of two copies; (6) addition of more than two copies. From these predictions, we considered the existence of a chromosome arm loss if more than a third of the genes from that chromosome arm were predicted to have levels 1 or 2. If more than a third of a 42 3.4. CELL ANNOTATION chromosome arm’s genes were predicted with levels 4 or greater, we considered to have a gain in that chromosome arm. Finally, an epithelial cell was considered a tumour cell if it had at least one chromosome arm alteration (gain or loss), while those with no large-scale alterations were classified as possibly normal. CNV predictions for some of the patients can be observed in the heatmap figures 53 to 55. For those tumour samples having too many epithelial cells classified as normal (figure 55 is an example of that), we further analysed the expression of genes known to be differentially expressed in CRC, to account for those cases where tumour do not have copy number variation. These genes are: MYC , RNF43 , AXIN2 , CTNNB1 , CD44 , MLH1 [167–172]. While MLH1 is under-expressed in tumours, all others are overexpressed, which happens in our atlas (figure 16). For each of those samples, we compared the expression distribution of these genes between normal, putative normal and putative tumour cells. The genes that had similar distribution between putative normal and putative tumour and different distribution between normal and putative normal were used to re-annotate the putative normal cells. For over expressed genes, the cells were re-annotated as tumour if the gene expression was higher than the median gene expression from normal cells. For the under expressed gene, the cells were re-annotated as tumour if the gene expression was smaller than the median gene expression from normal cells. Figure 16: Distribution of the expression of genes MYC , RNF43 , AXIN2 , CTNNB1 , CD44 , MLH1 in normal epithelial cells vs cells classified as tumour after CNV predictions. The following genes were used to annotate the group of normal epithelial cells (figure 17): LGR5 , SMOC2 , ASCL2 (stem cells); SOX9 , CDK6 , MUC4 , FABP5 , PLA2G2A , LCN2 (progenitor cells); MUC2 , ITLN1 , CLCA1 (secretory progenitors); TOP2A , CCNA2 , MCM5 , OLFM4 , SLC12A2 (transit-amplifying cells); POU2F3 , TRPM5 , SPIB , IL17RB , HTR3E (tuft cells); CHGA , CHGB , CPE , NEUROD1 , PYY (enteroendocrine cells); MUC2 , CLCA (goblet cells); FABP1 , SLC26A3 , TMEM37 , BEST4 (enterocyte/colonocyte cells); LYZ , CA7 , CA4 , SPIB , FKBP1A (paneth-like cells). Tumour cells were annotated according to the consensus molecular subtypes (CMS), using the R package CMScaller (v.0.99.2) [163]. Gene expression of cells was aggregated by sample before prediction. As such, a cell was annotated with the CMS type that was predicted to the respective sample. If the tool 43 CHAPTER 3. A COLORECTAL CANCER ATLAS OF SCRNASEQ DATA was not able to confidently classify a sample, that sample was classified with a Mixed type. Figure 17: Annotation of the normal epithelial cells. Top-left UMAP is coloured by the annotated cell-types. All other UMAPs are coloured by the RNA expression of some of the genes used to identify the cells, which goes from blue (low expression) to yellow (high expression). Grey cells have no expression of the respective gene. 3.5 Overview of the atlas Our atlas has a total of 163 810 cells, separated into 51 044 T-cells, 47 462 epithelial cells, 30 187 stromal cells, 17 674 B-cells, and 17 443 myeloid cells. The number of cells for each cell-subtype identified is summarised in table 2. 44 3.5. OVERVIEW OF THE ATLAS Table 2: Number of cells in each cell-subtype present in the CRC atlas. TA : transit-amplifying cells; CAFs : Cancer-associated fibroblasts; VSMCs : Vascular smooth muscle cells; DCs : Dendritic cells; LTi : Lymphoid tissue-inducer cells. Cancer cells 28 208 Normal Epithelial cells Progenitor cells 2818 19 254 Secretory progenitors 2 094 TA cells 3 862 Tuft cells 191 Enteroendocrine cells 74 Goblet cells Not immature 391 Immature 1 505 Enterocyte/colonocyte cells Immature 3 926 BEST4+ 536 Other 3 201 Paneth-like cells 656 Stromal cells Fibroblasts 14 007 30 187 CAFS 4914 Myofibroblasts 693 Pericytes 2014 Enteric glia cells 1 649 VSMCs 791 Endothelial cells tip-like vascular 3 934 stalk-like vascular 1 792 lymphatic 393 Myeloid cells Anti-inflammatory macro/mono lineage Other 4 963 17 443 SPP1+ 4 180 Pro-inflammatory macro/mono lineage Other 2 548 FCN1+ 1 639 conventional DCs 2072 plasmacytoid DCs 335 Mast cells 1 222 Unkown 484 B-cells Naive 2872 17 674 Memory 7 323 Plasma cells IgA+ 5 624 IgG+ 685 Proliferative 1 170 T-cells CD4+ Poliferative 377 51 044 Naive 8 484 Memory 13 484 Regulatory 6 801 Follicular 1 156 IL22+ 201 IL17+ 755 CD8+ Poliferative 439 Naive 414 Memory 6 534 CXCL13+ 1977 Cytotoxic 1 105 Unconventional 𝛾𝛿 2317 NKT 1 759 NK 1 382 45 CHAPTER 3. A COLORECTAL CANCER ATLAS OF SCRNASEQ DATA Table 2 (cont.) CD8𝛼𝛼 2 488 LTi cells 327 Double-Negative 1044 TOTAL 163 810 The consensus molecular subtype (CMS) of colorectal cancer most present in our atlas is CMS2, with 22 samples, followed by CMS1 (10), CMS4 (8) and CMS3 (4). Two samples were not confidently classified with any CMS and were thus labeled with a Mixed type. Gene set enrichment analysis (GSEA) with the assigned groups confirms the CMS classification (figure 18A). The group of samples with a CMS2 type are characterised by a high MSS/MSI ratio and activated MYC and E2F target gene sets. The CMS1 group, as expected, is MSI-like, but also revealed overexpression of genes related with MTORC1 signaling and E2F targets. CMS3 samples are MSS-like and have up-regulation of metabolic processes, namely fatty acid metabolism and oxidative phosphorilation. CMS4, in turn, have the characteristic strong activation of epithelial mesenchymal transition (EMT), TGF𝛽and angiogenesis. GSEA was not performed with the Mixed samples. Figure 18: (A) Heatmap of the results of gene set enrichment analysis using the CMS predictions. Blue and red mean underand overrepresentation, respectively, of the gene set. (B) UMAP visualization of tumour epithelial cells, coloured by CMS type. (C) Hierarchical cluster of samples, coloured by CMS type. (D) Distribution of the proportions of B-, epithelial, myeloid, stromal and Tcells by CMS type. The UMAP visualization of the tumour cells coloured by the respective sample’s CMS (figure 18B) shows that the cells from different samples but same CMS tend to group together. There is a clear separation of the CMS1, CMS2 and CMS4 cells, while CMS3 seems to overlap with the other CMS types. In fact, when performing hierarchical clustering of the samples (by aggregating the genes counts across 46 3.5. OVERVIEW OF THE ATLAS all cells of each sample, figure 18C), CMS3 samples clustered closer to the Mixed type. The remaining types tended to group with other samples of the same CMS type. We further assessed the proportion of the major cell-types across each CMS type (figure 18D). CMS4 showed a higher proportion in stromal cells than other types, as expected, but we did not observe a higher immune infiltration in CMS1 samples than other CMS types. CMS3 samples showed the highest proportions of Band Tcells, although it should be kept in mind that only 4 samples were classified as CMS3. The construction of this CRC atlas of scRNAseq data with cells spanning various T-cell subtypes from several different colorectal patients with tumour and normal matched samples allowed us to construct T-cell metabolic models that characterised not only the T-cells present in the tumour micro-environment, but also those in an unaffected part of the colon (or rectum) of the patient. 47 9 Modeling T-cells from the Colorectal Cancer Environment This chapter discusses the construction of a wide number of genome-scale metabolic models for different subtypes of T-cells from the micro-environment of colorectal cancer (CRC) and normal matched colon. We made use of the scRNAseq atlas for colorectal cancer (CRC) constructed in the previous chapter. A total of 196 models were constructed. Even though the structure of these models showed a lot of common pathways across cell-types and tissues of origin, the models for regulatory T-cells showed interesting differences between tissue of origin, for example. The prediction of models’ fluxes was very close to expected results regarding biomass and energy production. Cell-type prediction showed good results not only on datasets of gene expression but also on datasets of models’ structure (reaction presence) and of flux predictions. Models were further tested for gene essentiality and for the effect of different media conditions. All code was run with the R version 4.1.0 or the python version 3.6.9. More detailed information, including scripts, is available in the GitHub project ?iiTb,ff;Bi?m#X+QKfb`+`/QbQfJ2i#QHB+n JQ/2Hb. 4.1 Methods Starting from a generic human metabolic model, we reconstructed cell-type specific models with the aid of the scRNAseq data from the colorectal cancer (CRC) atlas constructed in the previous chapter. Not all samples from this atlas were used tough. We only used CRC patients with samples from both tumour and normal matched tissues, after filtering the samples with less than 1000 cells. Within the total of 14 patients with normal matched tissue samples, 7 had samples from tumour border and core. The cell-types considered for the construction of the models were: naïve CD8, memory CD8, proliferative CD8, cytotoxic CD8, naïve CD4, memory CD4, proliferative CD4, IL17+ CD4, follicular CD4, and regulatory CD4 T-cells. 48 4.1. METHODS 4.1.1 Generic human model The generic human model used to reconstruct context-specific models was Human-GEM [173], v.1.8.0, retrieved from the GitHub repository ?iiTb,ff;Bi?m#X+QKfavb"BQ*?HK2`bf>mKM@:1J in June 2021 in the SBML format. Human-GEM is a consensus genome-scale metabolic model (GSMM) created by integration of the preceding models HMR2.0 and Recon3D, and curation of different databases. This model contains a total of 13 802 reactions, 8 378 metabolites and 3 625 genes. Reactions span 9 compartments: extracellular, cytosol, inner mitochondria, mitochondria, endoplasmic reticulum, Golgi apparatus, lysosome, peroxisome, and nucleus. Before model reconstruction, we ensured consistency of the generic model chosen by identifying and removing blocked reactions. These are reactions whose maximum and minimum fluxes are null with open medium exchanges when performing flux variability analysis (FVA), meaning that they are not able to carry flux on any condition. Furthermore, we tested the consistent model for its capability to perform metabolic tasks [174] known to occur in human cells, as well as for biomass production under a specific medium (see section 4.1.3 for a detailed information on the medium used). The python module COBRApy [175] was used to handle the SBML file of the Human-GEM model, as well as to process this model into a consistent one. The evaluation of the consistent model regarding metabolic tasks and biomass production was performed through the python module troppo [176]. 4.1.2 Model Reconstruction To construct cell-type specific models from the human generic model, we used the fastCORE algorithm [177], which was available through the python module troppo [176]. fastCORE takes a set of reactions that have strong evidence to be active in the context of interest and searches for a subnetwork that contains all reactions from the core set and a minimal set of additional reactions to allow flux consistency through FVA. The choice of this algorithm was based on Vieira et al [178], where different combinations of data processing and reconstruction algorithms were tested for the reconstruction of tissue-specific models. The models returned by fastCORE were further gap-filled to ensure that flux through the biomass reaction when running the models with a specific medium (see section 4.1.3 for a detailed information on the medium used) is possible. To do so, the gap-filling was performed assuming no limitation (i.e., high upper bounds) of the metabolites that compose the medium. All models whose biomass flux was still null were excluded from the analysis. As such, in downstream analysis, whenever a model is predicted with null biomass, it can only be due to specific conditions (e.g., concentration of metabolites in medium, inhibited internal reactions, etc.). From expression data to core reaction set To construct the models using scRNAseq data, the raw gene counts were aggregated across cells of the same cell-type, generating a pseudo-bulk expression data for each cell-type in each sample. A cell-type in a sample was only considered for aggregation if more 49 CHAPTER 4. MODELING T-CELLS FROM THE COLORECTAL CANCER ENVIRONMENT than 5 cells representing that cell-type were available. A table summarising the cell-types represented in each scRNAseq sample and respective number of cells is present in the supplementary table 13. After aggregation, a gene expression matrix was created for each sample, where the rows corresponded to the genes and the columns to the cell-types considered for that sample. After this, the gene expression matrices were normalised into CPM counts. This was performed with the aid of the R packages Seurat [157] and SeuratDisk [158]. To get the set of core reactions for each cell-type model reconstruction, we used a similar approach to that of Richelle et al [179]. A gene is active in a sample’s cell-type if its expression is greater than a global maximum threshold, defined by the 75𝑡ℎ percentile of the expression distribution of all genes in all cell-types of that sample. On the other hand, a gene is inactive in a sample’s cell-type if its expression is smaller than a global minimum threshold, defined by the 10𝑡ℎ percentile of the expression distribution of all genes in all cell-types of that sample. Genes whose expression falls between these two thresholds, are considered active, or inactive, based on a local threshold. This local threshold is defined by 25𝑡ℎ percentile of the expression distribution of that gene in all cell-types of that sample. The percentiles were defined based on the best combinations found by Vieira et al [178] for the reconstruction of metabolic models. To achieve this, and for further downstream analysis, we first transformed the normalised pseudo-bulk data into gene activity scores (GASs), calculated using the above thresholds (table 3). Table 3: Calculation of the gene activity scores (GAS) from a sample. 𝑥𝑔: expression, in CPMs, of gene 𝑔;𝑔𝑙𝑜𝑏𝑎𝑙_𝑚𝑎𝑥: 75th percentile of the expression distribution of all genes in all cell-types of a sample; 𝑔𝑙𝑜𝑏𝑎𝑙_𝑚𝑖𝑛: 10th percentile of the expression distribution of all genes in all cell-types of a sample; 𝑙𝑜𝑐𝑎𝑙_𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑: 35th percentile of the expression distribution of a gene in all cell-types of a sample. Condition State GAS 𝑥𝑔≥𝑔𝑙𝑜𝑏𝑎𝑙_𝑚𝑎𝑥 Active log2 𝑥𝑔 𝑔𝑙𝑜𝑏𝑎𝑙_𝑚𝑎𝑥 𝑥𝑔≤𝑔𝑙𝑜𝑏𝑎𝑙_𝑚𝑖𝑛 Inactive log2 𝑥𝑔 𝑔𝑙𝑜𝑏𝑎𝑙_𝑚𝑖𝑛 𝑔𝑙𝑜𝑏𝑎𝑙_𝑚𝑖𝑛 ≤𝑥𝑔≥𝑔𝑙𝑜𝑏𝑎𝑙_𝑚𝑎𝑥 Moderate (Active or Inactive) log2 𝑥𝑔 𝑙𝑜𝑐𝑎𝑙_𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑 To obtain the set of core reactions to give as input to the reconstruction algorithm, we then need to translate the GASs into reaction activity scores (RASs). To do this, the gene-protein-reaction (GPR) rules present in the generic human model are used. These rules describe which combinations of genes are involved in the production of the enzyme(s) that catalyse the reactions. If two genes are necessary together (e.g., enzyme complex) for a reaction to occur, its GPR rule will be Gene A AND Gene B . On the other hand, if only one of two genes is necessary for a reaction (e.g., enzyme isoforms), the GPR rule is Gene A OR Gene B . To translate these GPR rules into the continuous RASs, the AND operators are replaced with a minimum function, assuming that the enzyme/reaction is limited by the lowest expressed gene. The OR operators, on the other hand, are replaced by a maximum function, assuming that the reaction activity is the activity of the highest expressed gene. A few examples on how to get the RASs from the GASs are present in table 4. 50 4.1. METHODS Table 4: Examples for calculating RASs. 𝑥𝑎,𝑥𝑏,𝑥𝑐: expression, in CPMs, of genes 𝑎,𝑏and 𝑐, respectively; 𝑚𝑎𝑥: maximum; 𝑚𝑖𝑛: minimum. GPR Replacement RAS (if 𝑥𝑎<𝑥𝑏<𝑥𝑐) 𝑥𝑎𝑂𝑅 𝑥𝑏𝑚𝑎𝑥(𝑥𝑎,𝑥 𝑏)𝑥𝑏 𝑥𝑎𝐴𝑁𝐷 𝑥𝑏𝑚𝑖𝑛(𝑥𝑎,𝑥 𝑏)𝑥𝑎 (𝑥𝑎𝐴𝑁𝐷 𝑥𝑏)𝑂𝑅 𝑥𝑐𝑚𝑎𝑥 (𝑚𝑖𝑛(𝑥𝑎,𝑥 𝑏),𝑥 𝑐)𝑥𝑐 (𝑥𝑎𝑂𝑅 𝑥𝑏)𝐴𝑁𝐷 𝑥𝑐𝑚𝑖𝑛(𝑚𝑎𝑥 (𝑥𝑎,𝑥 𝑏),𝑥 𝑐)𝑥𝑏 Finally, we only need to obtain the set of core reactions that will be used in the reconstruction model. By the way that the GASs and RASs were calculated (base 2 logarithm of a ratio), active reactions are simply those reactions whose RAS is positive. All these calculations were made using the python language. 4.1.3 Media used in the experiments To obtain flux predictions as close as possible to the in vivo reality of the T-cells, we sought to create a medium that best represented the metabolites in healthy human blood, instead of using a cell culture medium. We also wanted to recreate the metabolic medium in the tumour micro-environment, as it would be even closer to reality. However, there isn’t enough information on the concentration of metabolites in this situation, when compared to that found for metabolites under normal conditions. Nevertheless, several studies focused on how much the presence of metabolites change between the blood of normal individuals and CRC patients were found. Normal Medium The concentrations for the metabolites in the external compartment of our models were gathered from the Serum Metabolome Database (SMDB), which is integrated in the Human Metabolome Database (HMDB) [180]. For those model metabolites with more than one ‘normal’ concentration in the database, the final concentration was averaged. Metabolites with no information in the database regarding ‘normal’ concentrations were not included in the medium. Three metabolites (water, oxygen, and H+) were considered as always available, and thus their bounds were opened (maximum uptake set to 1 000). A total of 537 metabolites composed the final normal medium . Metabolic models work with reaction fluxes (mmol/gDW/h) instead of metabolite concentrations. Regarding the metabolites in a medium, they are represented by fluxes that characterise the entry rate of the metabolites in the cell. As such, the concentrations gathered from the database were transformed into fluxes in the following manner [181]: 𝐹𝑙𝑢𝑥𝑀𝑎 = [𝑀𝑎] [𝑐𝑒𝑙𝑙𝑠].𝑐𝑒𝑙𝑙 𝑤𝑒𝑖𝑔ℎ𝑡 .𝑡𝑖𝑚𝑒 (4.1) 51 CHAPTER 4. MODELING T-CELLS FROM THE COLORECTAL CANCER ENVIRONMENT The other pathways that also have little percentage of reactions present in the models are fatty acid recycling in the endoplasmic reticulum , bile acid recycling , protein assembly , drug metabolism , protein degradation , and acyl-CoA hydrolysis . Bile acid recycling , acyl-CoA hydrolysis and drug metabolism , for example, are pathways mostly associated with the liver [193, 194]. Two clear sets of pathways that are very present in all models are the pathways for 𝛽-oxidation of fatty acids in both the mitochondria and peroxisome, and those related to fatty acid metabolism (omega-3 and -6, linoleate), elongation and biosynthesis. Cholesterol biosynthesis (both Kandustch-Russel and Bloch pathways) is also present. Cholesterol has a critical role in T-cell function and signalling [195, 196] because T-cells rely greatly on motility and membrane-membrane interactions with other cells and cholesterol has been shown to be important in maintaining cell membrane stiffness [197]. Inhibiting enzymes from the cholesterol metabolism and transporters leads to changes in T-cell function, activation, and reprogramming. For example, treating a population of naive T-cells with an inhibitor of the enzyme that catalyses the ratelimiting step of cholesterol biosynthesis suppresses progression of their differentiation and cell cycle [197]. Both CD4 and CD8 naive T-cells reprogram their cholesterol metabolism upon activation, promoting both import and biosynthesis [198–200]. SREBP proteins, required for cholesterol biosynthesis, are essential for cytotoxic CD8 T-cells to acquire sufficient cholesterol levels that allow proliferation and acquisition of an effector phenotype [199]. IL17+ CD4 T-cells’ differentiation induced by ROR𝛾was preceded by enhanced cholesterol biosynthesis [200]. In fact, increased cholesterol content in the plasma membrane has been associated with a pro-inflammatory phenotype, even though membrane cholesterol enrichment in FOXP3+ CD4 regulatory T-cells did not alter their suppressogenic function [201]. Other pathways highly present across models include heparan sulfate degradation, keratan sulfate degradation, chondroitin sulfate degradation, aminoacyl-tRNA biosynthesis, terpenoid backbone biosynthesis, and leukotriene metabolism. The interaction of chemokine, integrins and selectins with glycosaminoglycans (GAGs) regulates the recruitment, adhesion, and migration of leukocytes from the circulation to the site of inflammation. GAGs are divided into four main groups: Heparan sulfate, chondroitin sulfate, keratan sulfate, and hyaluronic acid. Apart from hyaluronic acid, GAGs are attached to the core protein of proteoglycans (PGs)[202]. For instance, heparanase, the only mammalian enzyme that directly cleaves heparan sulfate, has also been described to be expressed in T-cells and up-regulated upon activation [203]. Furthermore, lymphocytes from peripheral blood mononuclear cells (PBMCs) of breast cancer patients were found to display higher heparanase expression than those from healthy patients [204]. Even though its expression has been shown to promote leukocytes migration and penetration of the basement membrane and blood vessel entry [203], it has been suggested that it can either be proor antitumorigenic [205]. Leukotrienes are one of the types of eicosanoids that derive from arachidonate. They have been linked to proliferation, apoptosis, cytokine production, differentiation, and chemotaxis in T-cells [206], and recognised as possibly having both proand antiinflammatory roles in the immune response of T-cells [207]. 58 4.4. PATHWAY COVERAGE ALOX5, an enzyme that catalyses the first reaction from this pathway (conversion of arachidonic acid into leukotriene A4, LTA4), was shown to occur in human T cell lines as well as in purified peripheral blood T cells, including naive and memory CD4 T-cells and cytotoxic CD8 T-cells [208]. The production of leukotriene B4 (LTB4) and cysteinyl leukotrienes (LTC4, LTD4 and LTE4) was also detected in various T-cell lines and primary cells [209, 210]. 4.4.2 Models’ structure differs between normal and tumour tissue We wanted to assess if any cell-type was affected by the tumour micro-environment, i.e., if the reaction presence differed between normal and tumour models. For that, we performed pathway differential analysis based only on comparing the pathway coverage between normaland tumourderived models of each cell-type. From all the cell-types, regulatory CD4 T-cells and cytotoxic CD8 T-cells showed the most interesting results. To calculate the most differently covered pathways for each cell-type, we first calculated the fold change on pathway coverage between normaland tumourderived models using the R package gtools [211]. Only those with absolute fold changes higher than 1.5 were subsequently tested using the non-parameteric test Mann-Whitney using the R package stats [212]. p-values were adjusted using the false discovery rate (FDR) method by Benjamin & Hochberg [213]. Pathways with an adjusted p-value smaller than 0.05 were considered differentially covered. Regulatory CD4 T-cells Thirty-two pathways are differentially covered between normal and tumour regulatory T-cell models (supplementary table 15). When plotting these pathways in a heatmap (figure 22), it is possible to see that, even though most models from normal matched tissue group together, there is not a complete separation between normal and tumour derived models. In fact, three main groups can be distinguished: (1) only models from normal matched tissues; (2) all models from tumour tissue classified as CMS4, two from normal tissue and one from tumour tissue classified as CMS3; (3) remaining tumour derived models. This goes in line with the previously shown results (figure 20) regarding the distribution of number of reactions in the models for each CMS type, where it was visible a decreasing trend: CMS1 > CMS2 > CMS3 > CMS4 > Normal Matched. From all the tumour derived models, the number of reactions in the models from CMS4 samples is the closest to that in the normal matched models. Also, one of the models from CMS3 has clearly less reactions than the other two CMS3 models. These 3 groups are still visible when visualizing all metabolic pathways’ coverage in a heatmap (supplementary figure 57). The pathways fatty acid biosynthesis (even-chain) , glycerolipid metabolism , keratan sulfate biosynthesis , fatty acid biosynthesis (odd-chain) give the clearest separation between tumour (heatmap groups (2) and (3)) and normal (heatmap group (1)) derived models. The TCA cycle ( Tricarboxylic acid cycle and glyoxylate/dicarboxylate metabolism ) is also more significantly present in tumour-derived models than normal matched-derived ones. 59 CHAPTER 4. MODELING T-CELLS FROM THE COLORECTAL CANCER ENVIRONMENT Figure 22: Differentially covered pathways between normaland tumourderived regulatory CD4 T-cell models. When comparing mice regulatory T-cells from a tumour site and the spleen, Pacella et al [214] found that those from the tumour site tended to acquire less FAs, despite the up-regulation of CD36 , a transporter of FAs into the cell. Considering that CPT1A , which moves FAs into the mitochondria, was significantly more expressed in tumour regulatory CD4 T-cells, and that these cells accumulated significantly more FAs, FA synthesis might be more present in tumour regulatory CD4 T-cells. Pacella et al [214] further showed higher consumption of the first intermediates of the TCA cycle by the tumour-infiltrating regulatory CD4 T-cells. There are a number of pathways whose distinction between groups (1) and (2) is not as clear. These pathways include: biotin metabolism , folate metabolism , sulfur metabolism , carnithine shuttle (mitochondrial) , ether lipid metabolism , phosphatidylinositol phosphate metabolism , fatty acid biosynthesis , pool reactions , bile acid biosynthesis and glycosphingolipid biosynthesis-lacto and neolacto series . This shows that regulatory CD4 T-cell models from CMS4 tumour seem to be less anti-inflammatory than the other tumour-derived models. Biotin deficiency decreases differentiation toward anti-inflammatory regulatory T-cells [215]. It thus makes sense that this pathway is very present in regulatory CD4 T-cells, even more so in tumour regulatory CD4 T-cells. Although most regulatory CD4 T-cell models have relatively high presence of biotin metabolism , it is visibly higher in group (3). Folate metabolism coverage is lower in groups (1) and (2). Yamaguchi et al [216] reported that folate 60 4.4. PATHWAY COVERAGE receptor 4 (FR4) is crucial for murine regulatory CD4 T-cell expansion in vivo and its blockade enhanced anti-tumour immunity. High levels of hydrogen sulfide (H2S), produced in the sulfur metabolism pathway, are known to limit the release of pro-inflammatory molecules and promote the secretion of anti-inflammatory cytokines [217]. Knockout of genes that encode for H2S-producing enzymes decreases regulatory CD4 T-cell proliferation [218]. The carnithine shuttle (mitochondrial) is in charge of transporting FAs in and out of the mitochondria, crucial for FAO to occur. Cytotoxic CD8 T-cells Twelve pathways are differentially covered between normal and tumour cytotoxic CD8 T-cell models (supplementary table 16). When plotting these pathways in a heatmap (figure 23), it is possible to see a clear separation between tumourand normalderived models. The pathways that are more clearly different between the two groups are oxidative phosphorylation (OXPHOS), biopterin metabolism , pantothenate and CoA biosynthesis , carnitine shuttle (peroxisomal) , folate metabolism , and estrogen metabolism . There are sphingolipid-related pathways ( glycosphingolipid biosynthesis-ganglio series , sphingolipid metabolism and glycosphingolipid metabolism ). There is one model from a CMS2 tumour, however, that seems closer to the normal-derived models, as pathway coverage is very similar apart from folate metabolism and estrogen metabolism . Also, from the 10 tumour-derived models, only 3 have relatively low coverage of estrogen metabolism . Figure 23: Differentially covered pathways between normaland tumourderived cytotoxic CD8 T-cell models. Folate deficiency was shown to reduce CD8+ T-cells capacity to proliferate in response to activation. This sensitivity to the lack of folate is higher in CD8+ T-cells than CD4+ T-cells [219]. We showed previously that the folate metabolism in regulatory CD4+ T-cell models from CMS4 tumours, and normal-matched tissues, was less covered than the other tumour-derived models. Regarding estrogen metabolism , estrogens are usually correlated with an immunoenhancement effect on the immune system, with CD8+ T-cells 61 CHAPTER 4. MODELING T-CELLS FROM THE COLORECTAL CANCER ENVIRONMENT showing a high response [220]. Navarro et al [221] showed that the female human CD8+ T-cells incubated with higher concentrations of the estrogen 17𝛽-estradiol (E2) were significantly more cytotoxic, while the male human CD8+ T-cells showed an increase but not significant. Of note, the 3 tumour-derived models with low estrogen metabolism are not all from male patients, and the other tumour-derived models are not all from female patients. The same occurs in the normal-derived models. Increased levels tetrahydrobiopterin (BH4), by either providing high concentrations of it in the medium or by overexpression of GCH1 , a gene from the biopterin metabolism , were shown to enhance proliferation of stimulated mice CD4 and CD8 T-cells. Cronin et al [222] further showed that stimulated BH4-deficient T-cells hold decreased mitochondrial respiration and oxygen consumption. Indeed, OXPHOS and pantothenate and CoA biosynthesis (essential for the TCA cycle) are two pathways that are also more present in the tumour-derived cytotoxic CD8 T-cell models than those derived from normal-matched tissue. 4.4.3 Differences in models’ structure between cell-types We also wanted to assess if any pair of cell-types’ metabolism was significantly different, i.e., if the pathway presence differed between them. For that, we performed pathway differential analysis based only on comparing the pathway coverage. From all the pairs of cell-types, regulatory vs IL17+ CD4 T-cells, and naïve vs proliferative CD8 T-cells showed the most interesting results. To calculate the most differentiated pathways for each pair of cell-types, we first calculated the fold change on pathway coverage between two cell-types using the R package gtools [211]. Only those with absolute fold changes higher than 1.5 were subsequently tested using the non-parameteric test MannWhitney using the R package stats [212]. p-values were adjusted using the Benjamin & Hochberg (FDR) method [213]. Pathways with an adjusted p-value smaller than 0.05 were considered differentially covered. IL17+ CD4 vs regulatory CD4 T-cells Nine pathways are differentially covered between IL17+ and regulatory CD4 T-cell models (supplementary table 18). When plotting these pathways in a heatmap (figure 24), it is possible to see a separation between these two cell-types. The pathways that are more clearly different between the two cell-types are biotin metabolism , fatty acid biosynthesis , glycerolipid metabolism , and keratan sulfate biosynthesis . We previously showed that biotin metabolism was more present in tumour-derived regulatory CD4 T-cell models than those derived from normal matched mucosa and that biotin deficiency revealed a decreased differentiation towards anti-inflammatory regulatory T-cells [215]. Now, when comparing IL17+ and regulatory CD4 T-cells, we were also able to capture the difference in biotin metabolism presence between these two cell-types, where this pathway is less present in IL17+ CD4 T-cell models. Indeed, biotin deficiency did not only lead to a decreased differentiation towards anti-inflammatory regulatory Tcells, but also induced Th1 and Th17-mediated pro-inflammatory responses [215]. This shows that, like corroborated in the literature, that biotin metabolism is important for an anti-inflammatory function of regulatory T-cells, while it is not necessary for IL17+ CD4 T-cell function. 62 4.4. PATHWAY COVERAGE Figure 24: Differentially covered pathways between IL17+ and regulatory CD4 T-cell models. The normal-derived regulatory CD4 T-cell models tend to cluster closer to IL17+ T-cell models than the other regulatory CD4 T-cell models. Interestingly, the pathways that are clearly different between these two cell-types were also considered differentially present between normaland tumourderived regulatory CD4 T-cell models. These pathways are biotin metabolism , fatty acid biosynthesis (even-chain) , fatty acid biosynthesis (odd-chain) , glycerolipid metabolism , and keratan sulphate biosynthesis . Protein related pathways (protein modification, assembly and degradation, and peptide metabolism) were also differentially more present in regulatory CD4 T-cells but were not considered differentially present between normaland tumourderived regulatory CD4 T-cell models. Naive CD8 vs proliferative CD8 T-cells Ten pathways are differentially covered between naïve and proliferative CD8 T-cell models (supplementary table 17). When plotting these pathways in a heatmap (figure 25), it is possible to see a clear separation between these two cell-types. The pathways that are more clearly different between the two cell-types are carnithine shuttle (endoplasmic reticular and mitochondrial), thiamine metabolism , sulfur metabolism , and O-glycan metabolism . The carnithine shuttle (mitochondrial) is in charge of transporting FAs in and out of the mitochondria. This is known to be crucial for fatty acid oxidation, process in which naïve T-cells greatly rely on [30]. However, our models show higher presence of carnithine shuttle (endoplasmic reticular and mitochondrial) in proliferative CD8 T-cell models than in naïve CD8 T-cell models, except for a few naïve models that also have a high presence of carnithine shuttle (endoplasmic reticular) . Still, Cong-Hui et al [223] showed that transport of fatty acids into the mitochondria by CPT1 may be required for anabolic processes that support healthy mitochondrial function and proliferation independent of fatty acid oxidation in cancer cells. It could be the same case for proliferative T-cells. Endoplasmic reticulum’s carnithine shuttle has been shown to function as an antioxidant to inhibit endoplasmic reticulum stress in proliferative cells other than T-cells [224]. 63 CHAPTER 4. MODELING T-CELLS FROM THE COLORECTAL CANCER ENVIRONMENT Figure 25: Differentially covered pathways between naive and proliferative CD8 T-cell models. Glycans are essential for signaling and cell-cell interactions, and they have been shown to be important in T-cell development, activity, differentiation and proliferation. Although O-glycans are de novo synthesized by cells, naïve T-cells cannot synthesize core 2 O-glycans. Following TCR stimulation, T-cells increase expression of 2 O-glycan synthesis, which allows proliferative T-cells to extravasate into non-lymphoid tissues [225]. 4.5 Predicting cell-types in the different stages of model reconstruction We sought to assess how well the data from different stages of model reconstruction predict the different cell-types. Using pseudo-bulk transcriptomics data leads to the best results out of all stages (figure 26), with an MCC of 0.793. Nevertheless, the model structure (presence/absence of reactions) is still a good classifier for cell-type, with an MCC of 0.526. If we filter the reactions without a GPR, this value increases to 0.583. Figure 26: MCC results when predicting cell-type using test samples from the pseudo-bulk RNAseq (in CPMs), reactions presence, or pFBA predicted fluxes datasets. 64 4.6. FLUX PREDICTIONS The fluxes predicted using the normal human blood medium lead to the worst results, with an MCC of 0.350. This suggests that the structure of a model is a better indicator of cell-type than predicted fluxes. Also, a good structure that correctly resembles the experimental data does not necessarily mean that predicted fluxes will be good, as the model’s objective and constraints like the medium used affects predictions. Nevertheless, the fluxes still have some predictive power, as the MCC is clearly higher than 0. This time, the pFBA fluxes without filtration of the reactions without a GPR give better predictions than those with (MCC 0.350 > 0.293). Even though different, we can thus expect high flux similarities between the models of the different Tcell subtypes. We calculated the Euclidean distances of the structure (supplementary figure 56A) and of the predicted fluxes under normal medium (supplementary figure 56B) between the different models. Indeed, all models were structurally very different from each other, while leading to very similar flux predictions. 4.6 Flux Predictions 4.6.1 Biomass and ATP production We chose different objectives for the models of proliferative T-cells and the remaining, non-proliferating, T-cells. As mentioned, biomass production is not the main objective of non-proliferative T-cells, with the production of energy being as important. With this, the proliferative T-cells’ models were optimised for biomass production while the remaining were optimised for both biomass and ATP production, using pFBA. Having this in mind, prediction results do show that the biomass of proliferative T-cells is higher than their naïve counterparts (figure 27), as expected. Cytotoxic CD8 T-cell follicular CD4 T-cells and IL17+ CD4 T-cell models also have relatively low biomass. Figure 27: Biomass flux prediction using normal human blood. Some models from naïve CD4, follicular CD4, memory CD4, memory CD8 and regulatory CD4 T-cells showed biomasses as big as those observed in the proliferative T-cells’ models, even though the objectives 65 CHAPTER 4. MODELING T-CELLS FROM THE COLORECTAL CANCER ENVIRONMENT are different. While for memory T-cells the flux going through the biomass reaction varies in both tumourand normal matchedderived models (figure 28A), all naive CD4 T-cell models, except one, that have high biomass are from tumour tissue samples. In line with the high biomass fluxes, all naïve CD4 T-cell models, except one, that have low ATP production levels are from tumour tissue samples (figure 28B). As for regulatory CD4 T-cell models, the models with high biomass are from tumour samples of CMS1 and CMS2 types (figure 28A). Interestingly, regulatory CD4 T-cell models from CMS4 and one from CMS3 clustered closer to those from normal tissue when looking into the pathway coverage in section 4.4.2, while models from CMS1 and CMS2 showed higher coverage of pathways connected to pro-inflammatory functions. While the models from tumour tissues of type CMS1 and CMS2 have high biomass flux, unlike the remaining regulatory CD4 T-cell models, only those from CMS1 have low ATP production (figure 28B). Figure 28: Biomass flux (A) and ATP production (B) predictions using normal human blood, separated by CMS type. 66 4.6. FLUX PREDICTIONS Regarding ATP production, it was expected that it would be higher in non-proliferative T-cells, due to the objectives chosen for pFBA. It is visible (figure 28A) that proliferative T-cells do have less ATP production than the remaining T-cell types. However, apart from cytotoxic CD8 T-cells, some of the models of these non-proliferating T-cell types have relatively lower or equal ATP production to those of proliferative T-cells. We further checked what would happen to the biomass and ATP production if all models, irrespective of cell-type, were optimised for biomass production only (supplementary figure 58). Even though we can see some difference in biomass fluxes between naive and proliferative CD8 T-cells, all cell-types show relatively the same biomass and ATP production flux. This shows the importance of choosing the right objective to predict a model’s fluxes that can depict the real fluxes of a cell-type. Previous studies [129, 137, 148], including one for T-cells, have also applied an objective that would combine biomass and ATP production for models characterising normal, non-proliferative, cells, while proliferative cell models (tumour and normal) would only have the maximisation of biomass as their objective. Performing prediction of the cell-types based on the pFBA fluxes predicted using biomass as the only objective for all models (figure 29) lead to worse results, especially when comparing those obtained using all reactions (MCC of 0.158 vs 0.350). Figure 29: MCC results when predicting cell-type using test samples from the pFBA predicted fluxes datasets. pFBA : pFBA predictions with the different objectives and all reactions; pFBA GPRS : pFBA predictions with the different objectives and only reactions with GPRs; pFBA Biomass : pFBA predictions with biomass as the only objective and all reactions; pFBA Biomass GPRS : pFBA predictions with biomass as the only objective and only reactions with GPRs. Ploting the biomass against the ATP production (figure 30) shows that most proliferative models have relatively low ATP production flux and a variable biomass flux that is mostly high (13 out of 19 models). Regarding the remaining models, which are non-proliferative and thus were maximised for both ATP and biomass production, most either have high biomass and low ATP production (28 out of 123 models), or low biomass and high ATP production (73 out of 123 models). This shows that even when maximising for both biomass and ATP production, the tendency of models with high biomass to have low ATP production, and vice versa, is still captured. Furthermore, most models where this is not observed are from tumour 67 CHAPTER 4. MODELING T-CELLS FROM THE COLORECTAL CANCER ENVIRONMENT Thus, all genes related to the uptake of each of these metabolites would have to be deleted to replicate the results obtained upon their removal from the medium. 4.7 Gene Essentiality 4.7.1 Validation with CRISPR-CAS9 studies We compared the gene essentiality predictions with two CRISPR-CAS9 studies that evaluated how essential the genes were for in vitro human CD4 [236] and CD8 [237] T-cell proliferation. Ting et al [236] evaluated a total of 2 658 genes, of which only 274 are present in genes tested in silico (figure 37A). Shifrut et al [237], on the other hand, tested a total of 9 329 genes, of which only 462 are part of the genes tested in silico (figure 37A). 324 genes tested in silico were not tested in any of the studies (figure 37A). It should be kept in mind that some of the genes could be present in two or all the three datasets, but the gene symbol used could be different. As a simple comparison was made between the symbols present in the three datasets, this might be the case for some genes. Figure 37: Venn diagrams of (A) the genes that were tested by the 3 different datasets, (B) the genes tested in silico and the genes reported as essential by the studies, (C) the in silico predictions and the genes tested, and (D) from the common genes between the studies and the pipeline, the predicted essential genes and the essential genes reported by the respective study, for both CD4 T-cells (top diagram) and CD8 T-cells (bottom diagram). 74 4.7. GENE ESSENTIALITY As no CRISPR-CAS9 studies testing cell proliferation were found for any specific CD4 or CD8 T-cell subtype, we joined the genes that were predicted as essential for each CD4 T-cell subtype together and compared them with those from the study of Ting et al [236]. The same thing was done for the CD8 T-cell subtypes, for comparison with Shifrut et al [237]. From the 932 genes tested in our pipeline, 14 were only essential in CD8 T-cell subtypes, while 32 were only essential in CD4 T-cell subtypes (figure 37C). 78 were essential in both CD4 and CD8 T-cell subtypes. From the 2 658 genes tested by Ting et al [236], 48 were considered essential. Of these 48, only 14 are present in the group of genes tested in silico (figure 37B). For the study of CD8 T-cells [237], 454 genes were found to be essential, in the total of 9 329 genes tested, and only 160 are present in the group of genes tested in our pipeline. The CRISPR-CAS9 studies share 19 essential genes, 7 of which are part of our tested genes. For a fair comparison between our predictions and the studies’ essential genes, we only used the in silico genes that were also tested in the studies (figure 37D). For the CD4 T-cells, only 5 were considered essential in both the study [236] and this work. While 9 in silico predictions were not essential in the study, 20 study’s essential genes were not considered so in our pipeline. For the CD8 T-cells, 27 were considered essential by both our pipeline and Shifrut et al [237]. However, 22 in silico essential genes were not considered essential by the study, and 133 of study’s essential genes were not predicted essential in our pipeline. There is a marked difference between the in silico predictions and the results reported by both studies. However, from the 274 common genes, 240 genes are not considered essential in both CD4 T-cells’ datasets, while 280 were not essential in both CD8 T-cells’ datasets, in the total of 462 common genes. It is good to keep in mind, however, three main points: (1) the CRISPR-CAS9 studies were performed with healthy human T-cells in vitro , while our models represent T-cells from a tumour patient (some from tumour tissue, others from normal matched tissue); (2) the medium used in the predictions does not resemble a medium usually used in vitro , but instead resembles the metabolites present in healthy human blood; and (3) the metabolic models can only predict the essentiality of a gene at the metabolism level, they do not account for regulatory and signalling pathways that can affect the essentiality of a gene. 4.7.2 Pathways affected across cell-types We wanted to check the pathways with the most essential genes, i.e., the most affected pathways. For that, from all the potential essential genes involved in a pathway, we calculated the percentage of those that were essential to each model. We focused on the most affected pathways whose number of potential essential genes is bigger than one (table 6). 75 CHAPTER 4. MODELING T-CELLS FROM THE COLORECTAL CANCER ENVIRONMENT Table 6: Top twenty pathways with the most median percentage of essential genes. Pathways % of essential genes Sulfur metabolism 33.3 Butanoate metabolism 25.0 Cholesterol biosynthesis 2 25.0 Glycosphingolipid biosynthesis-globo series 25.0 Phosphatidylinositol phosphate metabolism 25.0 Aminoacyl-tRNA biosynthesis 19.0 Fructose and mannose metabolism 16.7 Propanoate metabolism 16.7 ROS detoxification 16.7 Eicosanoid metabolism 16.7 Glycosphingolipid metabolism 16.7 Metabolism of other amino acids 14.3 Cholesterol metabolism 12.8 Tryptophan metabolism 12.5 Beta-alanine metabolism 12.5 Chondroitin / heparan sulfate biosynthesis 12.5 Amino sugar and nucleotide sugar metabolism 11.8 N-glycan metabolism 11.8 Pyrimidine metabolism 11.5 Pentose and glucuronate interconversions 11.1 Eicosanoid production has been described to normally have low constitutive levels [206], although eicosanoids are recognised as possibly having both proand antiinflammatory roles in the immune response of T-cells [207]. Some eicosanoids are produced from eicosapentaenoic acid (EPA) and dihomo𝛾-linolenic acid (DGLA), but most are derived from arachidonate [206]. These arachidonate-derived compounds include hydroxyeicosatetraenoic acids (HETEs), epoxides, hydroperoxyeicosatetraenoic acids (HPETEs), lipoxins (LXs), leukotrienes (LT), and prostanoids. All the 6 genes that affect at least one reaction in the Eicosanoid metabolism once deleted catalyse reactions that produce several arachidonate-derived eicosanoids (table 7). Table 7: Potential essential genes from the Eicosanoid metabolism pathway and respective products. Gene Product(s) ALOX12B 12-HPETE CYP2F1 12-HETE CYP4F8 18-HETE CYP4F12 10-HETE LTC4S Leukotriene C4 HPGD Prostaglandin F2𝛼/15-Keto-Prostaglandin F2A 76 4.7. GENE ESSENTIALITY As mentioned earlier, arachidonate-derived prostanoids and leukotrienes have been linked to proliferation, apoptosis, cytokine production, differentiation, and chemotaxis in T-cells [206]. CYP4F8 was tested in both CRISPR-CAS9 studies [236, 237] and LTC4S was tested in the CD8 T-cell related study [237]. Only LTC4S was considered essential. Glycosphingolipids (GSLs) are present in the cell membranes and have been related to T-cell activation, differentiation and function [238]. A group of GSLs are (iso)globosides, whose synthesis is represented in the metabolic models through the metabolic pathway Glycosphingolipid biosynthesis-globo series , one of the most affected pathways by gene essentiality. Indeed, the synthesis of this group of GSLs, as well as other types, has been shown to occur in normal human T-cells [238, 239]. From the 8 genes in this pathway that were tested for gene essentiality, 4 were tested by ther CD8 T-cell related study ( A4GALT , B3GALNT1 , ST8SIA1 and GLA )[237], and only 1 by the CD4 T-cell related study ( GLA )[236]. Only Shifrut et al [237] reported essential genes, namely A4GALT and B3GALNT1 . Other highly affected pathways include Butanoate metabolism , which has been shown to play an important role in the cytotoxic capacity of in vitro mouse CD8+ T-cells when these cells were supplemented with low-dose butyrate [240]. In fact, one of the genes tested from this pathway was considered as essential for the cytotoxic CD8+ T-cell type (the gene was essential for more than 90% of this cell-type’s models). As mentioned before, hydrogen sulfide (H2S) is produced in the sulfur metabolism pathway. Knockout of genes that encode for H2S-producing enzymes was shown to decrease regulatory CD4 T-cell proliferation [218]. Some genes catalyse reactions from different metabolic pathways. So, it should be kept in mind that only the inactivity of some of those reactions affected by a gene deletion might have an actual effect on the model’s objective. So, a pathway that is very affected upon gene deletions might not be ‘causing’ the inability of a model to perform its objective when that pathway is not fully working. Evaluating the essentiality of each reaction individually would aid in understanding what metabolic pathways could be seen as more important for a cell-type. Due to having a lot of genes associated to it, the Transport reactions pathway does not come up as one of the most affected pathways when looking into the percentage of genes that affect the models’ objective. However, it is the pathway that has the biggest number of genes considered as essential in each model. This was expected, due to the high number of potentially essential genes (201) that represent the transport of metabolites between compartments of a cell or the uptake of metabolites. Thus, we also checked the pathways with the highest median number of essential genes (table 8), as looking into the percentage of essential genes in the pathways can focus the analysis only on pathways with less genes. The 20 pathways with the most median number of genes that are essential are crucial for cell proliferation. There are several pathways related to amino acids metabolism ( glycine, serine and threonine metabolism , valine, leucine and isoleucine metabolism , phenylalanine, tyrosine and tryptophan metabolism , arginine and proline metabolism , alanine, aspartate and glutamate metabolism , and cysteine and methionine metabolism ), DNA/RNA production ( nucleotide metabolism and pyrimidine metabolism ), and fatty acids ( fatty acid oxidation and fatty acid biosynthesis ). Other important pathways include oxidative 77 CHAPTER 4. MODELING T-CELLS FROM THE COLORECTAL CANCER ENVIRONMENT phosphorylation , cholesterol metabolism , and aminoacyl-tRNA biosynthesis . Table 8: Top twenty pathways with the most median number of essential genes. Pathways Number of essential genes Transport reactions 40.5 Nucleotide metabolism 11.0 Oxidative phosphorylation 8.0 Sphingolipid metabolism 7.0 Fatty acid oxidation 6.0 Cholesterol metabolism 5.0 Aminoacyl-tRNA biosynthesis 4.0 N-glycan metabolism 4.0 Steroid metabolism 4.0 Glycerophospholipid metabolism 4.0 Pyrimidine metabolism 3.0 Glycine, serine and threonine metabolism 3.0 Valine, leucine, and isoleucine metabolism 3.0 Phenylalanine, tyrosine and tryptophan biosynthesis 3.0 Fatty acid biosynthesis 3.0 Bile acid biosynthesis 3.0 Drug metabolism 3.0 Arginine and proline metabolism 2.5 Propanoate metabolism 2.0 Alanine, aspartate and glutamate metabolism 2.0 Cysteine and methionine metabolism 2.0 4.8 Effect of a tumour blood medium We compared the pFBA predictions obtained using the normal human blood medium, constructed using information from the Serum Metabolome DataBase (SMDB), with those obtained using the tumour human blood medium. This tumour human blood medium was created from information on fold changes between the blood of normal and tumour patients. From the total of 537 metabolites in our normal human blood medium, we only found information for 84 of the metabolites (approximately 15.6%). We performed a non-parametric paired statistical test, Wilcoxon signed rank test , to find which models’ reaction fluxes differed significantly between the two types of medium. This was only observed for 53 out of the 116 models that returned feasible pFBA solutions (approximately 46% of the models, figure 38A). More than half of the models of naive CD4 (56.25%), memory CD4 (71.43%), follicular CD4 (55.56%), proliferative CD4 (57.14%), and proliferative CD8 (55.56%) T-cells significantly differed (figure 38B). No naive CD8 T-cell models were affected by the change in medium. We next checked what pathways changed the most when the medium was changed. For this, we 78 4.8. EFFECT OF A TUMOUR BLOOD MEDIUM selected the top 30 pathways with the highest median of ratio of affected reactions and plotted the distribution of this ratio across all models for each pathway (figure 38C). As expected, the most affected pathways are those directly affected by the metabolites whose concentration changed. Figure 38: (A) Number and (B) percentage of models, by cell-type, whose metabolism is significantly different when the medium was changed to a tumour blood-like one. (C) Top 30 pathways that changed the most when the medium was changed. The change in concentration of amino acids (asparagine, aspartate, glutamate, glutamine, glycine, histidine, isoleucine, leucine, lysine, methionine, serine, threonine tryptophan, tyrosine, and valine) affected the flux of a lot of the reactions from two pathways related to the metabolism of this type of metabolites: alanine, aspartate and glutamate metabolism ; cysteine and methionine metabolism . 𝛽 -oxidation of even-chain fatty acids (mitochondrial) and fatty acid activation (cytosolic) pathways were affected. The concentration of four unsaturated FAs (elaidate, linolenate, nervonic acid, oleate) and two saturated FAs (lauric acid, myristic acid) changed in the tumour medium. Moreover, glycolysis / gluconeogenesis was affected by changes in acetate, glucose, lactate, Pi, PEP and pyruvate levels, while the TCA cycle was affected by citrate, fumarate, glycerate, and malate. Other pathways were affected not because the metabolites that enter those pathways suffered a change in concentration in the medium, but because of the change in other pathways. This is the case of ROS detoxification . We finally checked if the flux of the biomass reaction also changed, considering that pathways like amino acids’ metabolism, fatty acid oxidation, glycolysis, TCA cycle, and oxidative phosphorylation were some of the most affected by the change in medium. Only 45 models (not necessarily part of those that had a significantly different metabolism) had different biomass flux. The biomass flux increased in 26 of these models (supplementary figure 63A). No naïve and cytotoxic CD8 T-cell models’ biomass changed (supplementary figure 63D). 79 CHAPTER 4. MODELING T-CELLS FROM THE COLORECTAL CANCER ENVIRONMENT 4.9 Discussion The metabolic models of T-cells structurally resemble well the respective gene expression data, with notso-distant predictive capabilities of the cell-types (figure 26). However, there are still a lot of similarities between the different types of T-cells, whether in gene expression (figure 19) or in the models’ structure (figure 21). Pathways like cholesterol biosynthesis, heparan sulfate degradation, keratan sulfate degradation, chondroitin sulfate degradation, leukotriene metabolism, and fatty acid related pathways are all highly present across all T-cell models and known for being important for T-cell function. Nevertheless, some groups of models show interesting differences in pathway coverage. For example, regulatory CD4 t-cell models from normal tissue samples seem to have less anti-inflammatory function due to less presence of pathways related to anti-inflammatory functions in T-cells (figure 22). Interestingly, we found models from CMS4 tumour samples to have less coverage of anti-inflammatory related pathways than the other tumour models. These regulatory CD4 T-cell CMS4 models also grouped closer to the normal-derived models, suggesting that regulatory CD4 T-cells from CMS4 tumours are not able to be as anti-inflammatory as they could. Further corroborating the less anti-inflammatory function of normalderived regulatory CD4 T-cells is the close clustering of these models with those of IL17+ CD4 T-cells (figure 24). Normal-derived regulatory CD4 T-cell models also have visibly less reactions than their tumourcounterparts (figure 20A), with models from CMS4 tumour being closer to normal-derived ones than the other tumour-derived models (figure 20B). When performing flux prediction with the pFBA approach, the biomass of normal-derived models, as well as of those from CMS4 tumours, of regulatory T-cells was null or very close to null, unlike those from other tumour-derived models (figure 28). All of this suggests that normal-derived and CMS4 regulatory CD4 T-cell models are also generally less metabolically active than the other regulatory CD4 T-cell models. A good structure that correctly resembles the experimental data does not necessarily mean that predicted fluxes will be good, however. Especially when knowing that these models do not take into consideration enzymatic constraints and gene / metabolic regulation, which are important aspects of a cell that greatly affect its metabolism. As such, flux predictions are still known to be quite imprecise [178], which is noticeable in our cell-type prediction results using reaction fluxes vs gene expression and reaction presence (figure 26). Nevertheless, the prediction of major metabolic aspects related to T-cells were correct. It was shown that the median biomass of proliferative T-cell models was higher than the naïve T-cell models, and that most models (114 out of 142 models with feasible pFBA solutions) have relatively high biomass and low ATP production or vice-versa. Fatty acid oxidation (FAO) is the pathway that contributes the most to produce FADH2 in all models. While models from non-proliferative T-cell models rely mostly in FAO for NADH production, proliferative T-cell models’ biggest source of NADH is glycolysis. Proliferative T-cells, characterised by relying in fatty acid synthesis instead of uptake (figure 1), showed to have null or very close to null uptake of fatty acids in the pFBA predictions (figure 32A). 80 4.9. DISCUSSION The individual deletion of tryptophan and glutamine from the medium decreases T-cell’s biomass to or very close to zero, while the deletion of nucleotides did not affect the biomass, as expected [226, 227, 234]. Even though biomass was expected to decrease, the deletion of glucose did not affect the biomass of the models, which could be explained by existence of alternative metabolites in the medium that are not present in in vitro studies and may reduce the dependence on glucose. All flux predictions mentioned above were made using a normal human blood medium. Thus, we tried to construct a medium as close as possible to what a tumour blood medium could be. As there were no quantitative measurements of a tumour blood medium, unlike the normal one, we constructed the tumour medium looking into studies that compared peak intensities from mass spectrometry between normal and tumour blood samples. We were able to change only 84 of the 537 metabolites we have on the normal medium, which might explain why only 46% of the models differentially changed with the new medium. Another reason for a low percentage of models that differed might be the fact that cells always try to adjust so that the metabolism differs as little as possible and this was captured by the models. Nevertheless, the most affected pathways were those directly affected by the metabolites whose concentration changed. A better quantitative characterisation of a tumour blood medium would be very helpful in the future to better estimate the reaction fluxes of these cell-types in the colorectal cancer environment. 81 8 Benchmarking of Tumour Deconvolution Methods Although scRNAseq gives valuable information about cell-types and states at single-cell resolution, this technique is not an accurate representation of cell-type proportions in a sample. Bulk RNA sequencing, on the other hand, is used for diagnostic purposes in routine clinical settings, allowing an unprecedented amount of data describing tumours. Cell-type profiles can be recovered from bulk RNAseq data and used for cell-type metabolic model reconstruction. The proportions can aid in the construction of community models when modelling the tumour micro-environment. Having this in mind, we tested and compared several tumour deconvolution methods that were developed to use scRNAseq data as reference to estimate cell-type proportions. Our CRC atlas of scRNAseq data was used. 5.1 Methods 5.1.1 Bulk RNAseq Deconvolution The mixed gene expression profile of a tumour micro-environment can be seen as the linear combination of the gene expression profiles of each cell type in that micro-environment and the respective proportions of those cell types. This relation can be portrayed by equation 5.1: 𝐴.𝑋=𝐵(5.1) Here, 𝐵represents a matrix of the gene expression profile for a tumour micro-environment, the bulk data, with an expression value for each gene (rows) in each sample (columns). 𝐴is a matrix representing the cell type gene expression profiles, with an expression value for each gene (rows) in each cell type (columns). This matrix is often referred to as reference matrix. Genes in 𝐴are the same as those in 𝐵. Finally, 𝑋represents the proportions of these cells in the micro-environment. This matrix contains the proportion of each cell type (rows) for each sample (columns). 82 5.1. METHODS In this chapter, we focus on deconvoluting bulk RNAseq data (𝐵) from CRC samples to obtain the cell-type proportions (𝑋). 5.1.2 Bulk data and ground-truth proportions We used in-house bulk RNAseq tumour samples from CRC patients, gathered at the Leiden University Medical Center (Leiden cohort) [241] and whose cell counts were calculated by using Hyperion mass cytometry imaging. A total of 21 samples were available. The ground-truth phenotypes were divided into 9 different cell-types plus a group named ’other cells’ (see section 5.1.3). To calculate the proportion of a cell-type in a sample (figure 39B), we aggregated the cell counts of all phenotypes of that cell-type (figure 39A and supplementary table 19) and divided them by the total number of cells in the sample. Figure 39: Ground-truth (A) cell counts and (B) proportions of the cell-types used for deconvolution. 5.1.3 CRC atlas as the reference data The CRC atlas of scRNAseq data constructed was used to build the reference matrix. Usually, scRNAseq data from peripheral blood mononuclear cells (PBMCs) is used as reference for the deconvolution of tumour samples. However, this will always lead to incorrect estimation of immune cell-type proportions: (1) gene expression profile of immune cells from peripheral blood might not represent well the expression of tumour-infiltrating immune cells; (2) several cell-types present in the tumour are not PBMCs and are thus not accounted for when deconvoluting the samples (all methods normalize the proportions estimation to sum to 1). Mapping cell-types between CRC atlas and ground-truth proportions To test the different deconvolution methods, we had to find the best overlap between the atlas’s cell-types and the groundtruth’s cell phenotypes. While doing so, we found the best balance between having appropriate cell-types matched and enough subtype detail to deconvolute the most important cell subtypes. All cell-types / 83 CHAPTER 5. BENCHMARKING OF TUMOUR DECONVOLUTION METHODS overall RMSE of this method (0.542). Regarding correlation, it is evident that no method scores higher than 0.8 in any cell-type, although the three best methods overall have a correlation greater than 0.8. Figure 43: Cell-type correlation vs RMSE for all methods. A good method has a high correlation and a small RMSE. Grey horizontal and vertical lines mark a correlation of 0.5 and a RMSE of 0.2, respectively. We next evaluated which method is the best for predicting each cell-type individually (table 10). For this, we focused on RMSE and correlations metrics (figure 43), but also on visual inspection of the estimated versus ground-truth proportions (supplementary figures 65 to 73). Table 10: Best methods for each cell-type. Cell-type Best Method Cancer cells Scaden Stromal cells BisqueRNA Macro/mono lineage cells BSeqSC B-cells Scaden CD4 T-cells CIBERSORTx Regulatory CD4 T-cells Scaden CD8 T-cells DigitalDLSorter Proliferative T-cells MuSiC_wGrouping NK cells Scaden The best method to estimate cancer cells is Scaden , followed by CIBERSORTx , DigitalDLSorter , AutoGeneS_nusvr and DWLS (figure 43 and supplementary figure 65). These are the methods that are the 90 5.2. RESULTS best overall. As no other cell-type shows these methods as the clear best ones, this shows that the celltypes that are constantly more present in the samples affect more the overall RMSE and correlation of a method. While stromal cells (figure 43 and supplementary figure 66) are best estimated using BisqueRNA , followed by BseqSC , BseqSC is the best one at estimating the proportions for the cells from the macro/mono lineage (figure 43 and supplementary figure 67). Even though MuSiC_woGrouping and CIBERSORTx have better correlation and slight worse RMSE at estimating B-cells (figure 43), the estimations versus the ground-truth plots (supplementary figure 68) show that Scaden is the right choice. The very low RMSE is not due to all samples being estimated as zero and the bad correlation is highly influenced by two samples, without which Scaden ’s correlation would be better and with which MuSiC_woGrouping and CIBERSORTx have better correlation. When assessing the methods’ RMSE and correlation metrics for regulatory CD4 T-cells (figure 43), DWLS , AutoGeneS_nusvr and Scaden are the three best. However, visually comparing estimated and ground-truth proportions (supplementary figure 70) shows that AutoGeneS_nusvr gives negative proportions, while most samples were predicted by DWLS to not have regulatory CD4 T-cells. Scaden underestimates regulatory CD4 T-cells, but is still able to detect them. Regarding CD4 (figure 43 and supplementary figure 69) and CD8 T-cells (figure 43 and supplementary figure 71), the best methods are CIBERSORTx and DigitalDLSorter , respectively. With better RMSE than DigitalDLSorter , DWLS and MOMF , MuSiC_wGrouping is the best method to predict proliferative T-cells (figure 43 and supplementary figure 72). Even though they have better correlation than MuSiC_wGrouping , MOMF only detects the presence of proliferative T-cells in 4 samples, DWLS returns negative proportions for 4 of the samples, and DigitalDLSorter over-estimates proportions a lot more than the previous methods. NK cells is by far the most difficult cell-type to estimate, with no method performing particularly well (figure 43 and supplementary figure 73). Scaden , followed by DigitalDLSorter , seem to be the methods that better capture NK cells proportions, even though the respective correlations are negative. With these results in mind, we decided to see what would happen if a pipeline using the best method for each cell-type was used (figure 44). Combining the best methods for each cell-type creates better overall estimations, with a correlation that now surpasses 0.9 and a RMSE smaller than 0.1 (figure 44A). We tested how different the results would be between a simple combination of the results from the different methods ( Combined ) and normalising the combined results of each sample to sum to 1 ( Combined_norm , figures 44C and 44D). Overall, normalising the estimation leads to better results, although the difference is of only 0.013 and -0.001 for the correlation (0.925 - 0.912) and RMSE (0.085 – 0.086) respectively. We tried to find an independent dataset for colorectal cancer that also provided information on celltype proportions based on mass cytometry to further corroborate these findings. However, we could not find one. Nevertheless, we checked the overall correlation and RMSE by sample, to check how many samples benefit from a combined pipeline as opposed to the Scaden method (figure 44B). Indeed, most samples (76%) were better estimated using the combined approach. The samples that greatly benefited 91 CHAPTER 5. BENCHMARKING OF TUMOUR DECONVOLUTION METHODS were NIC13, NIC24 and NIC4, which showed some of the worst estimations throughout most methods (supplementary figure 85) including Scaden . The samples that seem to have not benefitted from the combined pipeline were NIC16, NIC20, NIC27, NIC6, and NIC7. Still, the metrics did not change much. Figure 44: (A) Overall correlation vs RMSE for all methods, including the methods combining the best methods for each cell-type ( Combined and Combined_norm ). (B) Sample correlation vs RMSE for Scaden , Combined and Combined_norm . A good method has a high correlation and a small RMSE. Grey horizontal and vertical lines mark a correlation of 0.5 and a RMSE of 0.2, respectively. Scatter plots of estimated vs ground-truth proportions for (C) Combined and (D) Combined_norm . 5.2.3 RNA content bias correction does not improve predictions Both RMSE and correlation metrics (figure 45) of predictions with correction for mRNA content bias do not show improvement from those without any correction, even though the decrease in prediction ability is not too steep. There is only one method with improved predictions upon correction on the reference matrix ( Corrected Before ). This method is Scaden , but the improvement is not high and, in both cases, (correction and no correction) Scaden is one of the three best methods. The improvement observed in the Scaden method seems to be mostly due to an improvement in predicting the proportions of the samples NIC3, NIC4, NIC13, NIC15, NIC24 and NIC29 (supplementary figure 86). Furthermore, the methods that already have correction for RNA content ( AutoGeneS_linear , AutoGeneS_nnls , AutoGeneS_nusvr and MuSiC_woGrouping ) have some of the worst predictions (figure 45 and 41), with only AutoGeneS_nusvr showing satisfactory results overall. 92 5.2. RESULTS Figure 45: Pearson correlation (A) and (B) RMSE values of the methods without and with correction. Methods were tested for correction on the reference matrix ( Corrected Before ) and on the estimated proportions ( Corrected After ). This goes against previous findings [252, 253], where correction for mRNA content bias was shown to improve deconvolution. However, it is noteworthy that Zaitsev et al [252], which tested the correction for this bias on the predicted proportions, used a dataset created with only two cell-types whose total RNA content is considerably different and was measured by analysing bulk RNA data of pure samples. The small number of cell-types, the big difference in RNA content and its calculation from the exact same cells used for creation of the samples used in deconvolution might be decisive for a good correction. In our case, there are much more than 2 cell-types (but the exact cell-types present in the bulk samples is not known), and the total RNA content is known through a scRNAseq dataset independent from the bulk samples. Beyond this, Zaitsev et al [252] developed a complete deconvolution method, i.e., a deconvolution method that does not use prior information on cell-type gene expression to estimate the proportions. Interestingly, the predictions with correction for RNA content on the estimated proportions lead to the worst results overall in our work. Sosina et al [253] tested the correction for this bias on the reference matrix using the MuSiC method with default parameters. In that study, the single-cell reference used was obtained using the Fluidigm C1 system, which normalises cDNA libraries to the same concentration prior to sequencing and thus removes potential variability in RNA abundances across cell-types. The datasets collected in this work were constructed using 10XGenomics. If we explore the results obtained without any correction for RNA content for each cell-type individually, it is possible to see that the biggest cell-type, cancer cells, are not over-estimated in any of the methods that does not use correction (supplementary figure 65). In fact, all methods that already come with RNA content bias correction ( AutoGeneS and MuSiC_woGrouping ) under-estimate cancer cells. It is noteworthy, however, that Scaden slightly improves the estimation of cancer cells upon correction of the estimated 93 CHAPTER 5. BENCHMARKING OF TUMOUR DECONVOLUTION METHODS proportions (figure 46 and supplementary figure 74), whose uncorrected version was the best method at predicting cancer cells. The only other two cell-types whose RNA content bias correction seems to help is stromal and macro/mono lineage cells. Figure 46: RMSE (A) and (B) pearson correlation values of the methods without and with correction, separated by cell-type. Methods were tested for correction on the reference matrix ( Corrected Before ) and on the estimated proportions ( Corrected After ). Stromal cells was the cell-type with the most amount of methods where correction for RNA content improved the predictions (7 out of 9 methods – (figure 46A). Of these, only 3 showed better RMSE than the best uncorrected method for this cell-type. Comparing the estimated and ground-truth proportions (supplementary figure 75), however, shows that the uncorrected method ( BisqueRNA ) might still be the best option to predict stromal cells proportions. For the macro/mono lineage cell-type, BsqesSC corrected for RNA content on the reference matrix shows visibly better results than its uncorrected version (figure 46 and supplementary figure 76). Regarding the remaining cell-types, some methods benefited from this correction, but none showed better results than the uncorrected method that was best for the respective cell-type (figure 46 and supplementary figures 74 to 82). 5.2.4 Effect of real cell-types’ proportions in samples estimations Finally, we were interested in assessing how the amount of cells or RNA content in the samples affected the correct estimation of the samples’ proportions. The total number of cells in a sample seems to not affect much its predictions for most methods (figure 47A). For AutoGeneS_nnls , CIBERSORTx , MuSiC_wGrouping , however, the greater the number of 94 5.2. RESULTS cells, the worse the RMSE of the sample. The correlation between the RMSE and the total number of cells is bigger than 0.34 in these methods. Even though there is some positive correlation between the total number of cells in a sample and its total number of reads (0.309, supplementary figure 83), the samples’ RMSEs from AutoGeneS_nnls and CIBERSORTx improve with the increase of total read counts (-0.316 and -0.233, respectively). The total number of reads also positively affects other methods (figure 47B), like AutoGeneS_nusvr (-0.481), MuSiC_woGrouping (-0.478) and SCDC (-0.443). Figure 47: For each method, RMSE of the samples is mapped against their corresponding (A) total number of cells and (B) total number of read counts. 95 CHAPTER 5. BENCHMARKING OF TUMOUR DECONVOLUTION METHODS Even though the total number of cells or reads might not seem to affect much the overall prediction of samples’ proportions, the proportion of certain cell-types in the samples might affect such predictions. This was hinted before, with cancer and stromal cells being the cell-types that most affect the overall predictions of the methods (figures 42 and supplementary figure 64). With different levels of correlation, the best methods to predict cancer cells are positively affected by the increase of the proportions the cancer cells in the samples (figure 48). Apart from AutoGeneS (-0.034) and DWLS (-0.157), these methods have a correlation smaller than -0.3. DigitalDLSorter is the most affected one, with a correlation of -0.729. All other models, which did not perform as well, are negatively affected by the increase of cancer cells proportions in the samples, especially MuSiC_wGrouping (0.98), BisqueRNA (0.661) and MOMF (0.885). Figure 48: For each method, RMSE of the samples is mapped against their corresponding proportion of cancer cells. Regarding stromal cells (supplementary figure 84), CIBERSORTx , AutoGeneS_nusvr and Scaden are not too affected by the proportion that stromal cells hold in a sample, while DigitalDLSorter is negatively affected, with a correlation of 0.651. The best methods overall are negatively affected by the increasing proportion that cells not deconvoluted ( Other cells ) have in the samples (figure 49). As expected, increasing the proportion of a group of cells not deconvoluted by the methods hinders estimations because the methods assume that the deconvoluted cell-types are the only cell-types present, by making proportions of estimated cell-types to have to sum up to 1. Despite that, most of these methods still score RMSEs smaller than 0.2 for the samples with higher proportion of Other cells . Surprinsingly, some methods benefit from a higher presence of Other cells . These methods, however, score an RMSE higher than 0.2 for practically all samples. 96 5.2. RESULTS Figure 49: For each method, RMSE of the samples is mapped against their corresponding proportion of other cells. The same trend seen for the Other cells seems to happen for the immune cells (figure 50). The best methods overall are negatively affected by the increasing proportion of immune cells, while the worst methods are either positively affected or not affected. Figure 50: For each method, RMSE of the samples is mapped against their corresponding proportion of immune cells. 97 CHAPTER 5. BENCHMARKING OF TUMOUR DECONVOLUTION METHODS 5.3 Discussion In this chapter, we wanted to assess the best methods to use for the deconvolution of colorectal cancer samples of bulk RNAseq data. There are two studies [255, 256] that performed benchmarking of tumour deconvolution methods that include those that use scRNAseq data as reference. However, none were focused on colorectal cancer and methods’ performance can change according to the tissue to be deconvoluted. The authors from one of these studies, a DREAM challenge that sought to evaluate a wide range of deconvolution methods and pipelines [256], even concluded that the best method may be problem specific. Although this challenge covered a lot of methods, very few were based on scRNAseq data and those methods were compared using their own default single-cell reference, which does not allow to properly assess if a better performance is due to the reference used or the deconvolution algorithm. Furthermore, they only used correlation metrics (pearson and spearman) to evaluate the results. We found that the best method overall was CIBERSORTx , closely followed by DigitalDLSorter and Scaden . The authors from the DREAM challenge [256] also reported CIBERSORTx as the top-performing method for their datasets. However, when looking into which methods are the best for each cell-type, CIBERSORTx is only the best for CD4 T-cells and DigitalDLSorter for CD8 T-cells. Scaden , on the other hand, was found to be the best method for 4 cell-types (cancer cells, B-cells, regulatory CD4 T-cells, and NK cells). BisqueRNA , BseqSC and MuSiC_wGrouping were some of the worst methods overall but they were the best at predicting stromall cells, macro/mono lineage cells, and proliferative T-cells, respectively. DWLS , all AutoGeneS methods, MOMF , MuSiC_woGrouping , and SCDC were not the best in any cell-type. Combining the best methods for each cell-type seems promising for estimating the proportions of all cell-types in the samples. However, further testing with an independent dataset with known cell-type proportions should be pursued, to assure that this is scalable for other CRC datasets and samples and does not work only for our specific case. The authors from the DREAM challenge [256] also combined the outputs of different methods but their ensemble method found only modest improvement relative to the top-scoring individual methods. We also found that certain cell-types are more easily deconvoluted than others, with NK cells being, by far, the most difficult cell-type to estimate. CD4 and regulatory CD4 T-cells were also rather difficult to estimate. The total number of read counts seems to influence sample estimations, with those having higher total mRNA content showing better results in most methods. We were also able to show that, as expected, the proportion of cells in the samples that are not part of a cell-type considered for deconvolution hinders the estimation of the deconvoluted cell-types. Notwithstanding, the best methods overall still showed low (< 0.2) RMSE for most samples, in a dataset where the proportion of unknown cells ( Other cells ) in the samples goes up to 0.15. Finally, assuming that bigger cell-types would tend to be over-estimated, it would be expected that cancer, stromal and macro/mono lineage cells would suffer over-estimation across most methods without RNA content bias correction. This was not the case, however. All methods that already implemented 98 5.3. DISCUSSION correction ( AutoGeneS and MuSiC_woGrouping ) under-estimated cancer cells and only some of the other methods benefitted (modestly) from this correction. Stromal cells also slightly benefited from RNA content bias correction. Only macro/mono lineage cells showed to visibly benefit from this correction when comparing the best methods with and without correction. These results give a clear indication of the best methods to use when assessing the proportions of different cell-types in a colorectal cancer sample of bulk RNAseq data. Nevertheless, further work is needed. For example, it is necessary to assess cell-types’ spillover in CRC samples, i.e., what other celltypes are being attributed signal that belongs to a certain cell-type. This can only be done by estimating proportions of samples purified for each cell-type and could help understand to what detail we could estimate cell subtypes. 99 BIBLIOGRAPHY [45] R. Michalek, V. Gerriets, S. Jacobs, et al. “Cutting edge: distinct glycolytic and lipid oxidative metabolic programs are essential for effector and regulatory CD4+ T cell subsets”. In: The Journal of Immunology 186.6 (2011), pp. 3299–3303. [46] L. Shi, R. Wang, G. Huang, et al. “HIF1𝛼–dependent glycolytic pathway orchestrates a metabolic checkpoint for the differentiation of TH17 and Treg cells”. In: Journal of Experimental Medicine 208.7 (2011), pp. 1367–1376. [47] D. Cluxton, B. Moran, and J. Fletcher. “Differential regulation of human Treg and Th17 cells by fatty acid synthesis and glycolysis”. In: Frontiers in Immunology 10 (2019), p. 115. [48] V. Gerriets, R. Kishton, M. Johnson, et al. “Foxp3 and Toll-like receptor signaling balance Treg cell anabolic metabolism for suppression”. In: Nature Immunology 17.12 (2016), p. 1459. [49] L. Berod, C. Friedrich, A. Nandan, et al. “De novo fatty acid synthesis controls the fate between regulatory T and T helper 17 cells”. In: Nature Medicine 20.11 (2014), p. 1327. [50] H. Zeng, K. Yang, C. Cloer, et al. “mTORC1 couples immune signals and metabolic programming to establish Treg-cell function”. In: Nature 499.7459 (2013), p. 485. [51] B. Raud, D. Roy, A. Divakaruni, et al. “Etomoxir actions on regulatory and memory T cells are independent of Cpt1a-mediated fatty acid oxidation”. In: Cell Metabolism 28.3 (2018), pp. 504– 515. [52] N. Pavlova and C. Thompson. “The emerging hallmarks of cancer metabolism”. In: Cell Metabolism 23.1 (2016), pp. 27–47. [53] T. Murakami, T. Nishiyama, T. Shirotani, et al. “Identification of two enhancer elements in the gene encoding the type 1 glucose transporter from the mouse which are responsive to serum, growth factor, and oncogenes.” In: Journal of Biological Chemistry 267.13 (1992), pp. 9300–9306. [54] P. Gao, I. Tchernyshyov, T. Chang, et al. “c-Myc suppression of miR-23a/b enhances mitochondrial glutaminase expression and glutamine metabolism”. In: Nature 458.7239 (2009), pp. 762–765. [55] M. Conrad and H. Sato. “The oxidative stress-inducible cystine/glutamate antiporter, system x (c) (-): cystine supplier and beyond”. In: Amino Acids 42.1 (2012), pp. 231–246. [56] D. Wu et al. “Hydrogen sulfide in cancer: friend or foe?” In: Nitric Oxide 50 (2015), pp. 38–45. [57] K. Rajagopalan and R. DeBerardinis. “Role of glutamine in cancer: therapeutic and imaging implications”. In: Journal of Nuclear Medicine 52.7 (2011), pp. 1005–1008. [58] J. Kamphorst, J. Cross, J. Fan, et al. “Hypoxic and Ras-transformed cells support growth by scavenging unsaturated fatty acids from lysophospholipids”. In: Proceedings of the National Academy of Sciences 110.22 (2013), pp. 8882–8887. [59] H. Ying, A. Kimmelman, C. Lyssiotis, et al. “Oncogenic Kras maintains pancreatic tumors through regulation of anabolic glucose metabolism”. In: Cell 149.3 (2012), pp. 656–670. 106 BIBLIOGRAPHY [60] R. Nilsson, M. Jain, N. Madhusudhan, et al. “Metabolic enzyme expression highlights a key role for MTHFD2 and the mitochondrial folate pathway in cancer”. In: Nature Communications 5 (2014), p. 3128. [61] H. Christofk, M. Vander Heiden, M. Harris, et al. “The M2 splice isoform of pyruvate kinase is important for cancer metabolism and tumour growth”. In: Nature 452.7184 (2008), p. 230. [62] B. Chaneton, P. Hillmann, L. Zheng, et al. “Serine is a natural ligand and allosteric activator of pyruvate kinase M2”. In: Nature 491.7424 (2012), p. 458. [63] K. Sircar, H. Huang, L. Hu, et al. “Integrative molecular profiling reveals asparagine synthetase is a target in castration-resistant prostate cancer”. In: The American Journal of Pathology 180.3 (2012), pp. 895–903. [64] J. Zhang, J. Fan, S. Venneti, et al. “Asparagine plays a critical role in regulating cellular adaptation to glutamine depletion”. In: Molecular Cell 56.2 (2014), pp. 205–218. [65] R. DeBerardinis and N. Chandel. “Fundamentals of cancer metabolism”. In: Science Advances 2.5 (2016), e1600200. [66] T. Bowles, R. Kim, J. Galante, et al. “Pancreatic cancer cell lines deficient in argininosuccinate synthetase are sensitive to arginine deprivation by arginine deiminase”. In: International Journal of Cancer 123.8 (2008), pp. 1950–1955. [67] C. Yoon, Y. Shim, E. Kim, et al. “Renal cell carcinoma does not express argininosuccinate synthetase and is highly sensitive to arginine deprivation via arginine deiminase”. In: International Journal of Cancer 120.4 (2007), pp. 897–905. [68] J. Cunningham, M. Moreno, A. Lodi, et al. “Protein and nucleotide biosynthesis are coupled by a single rate-limiting enzyme, PRPS2, to drive cancer”. In: Cell 157.5 (2014), pp. 1088–1103. [69] S. Eberhardy and P. Farnham. “c-Myc mediates activation of the cad promoter via a post-RNA polymerase II recruitment mechanism”. In: Journal of Biological Chemistry 276.51 (2001), pp. 48562– 48571. [70] D. Hanahan and L. Coussens. “Accessories to the crime: functions of cells recruited to the tumor microenvironment”. In: Cancer Cell 21.3 (2012), pp. 309–322. [71] T. Whiteside. “The tumor microenvironment and its role in promoting tumor growth”. In: Oncogene 27.45 (2008), p. 5904. [72] I. Kareva and P. Hahnfeldt. “The emerging ”hallmarks”of metabolic reprogramming and immune evasion: distinct or linked?” In: Cancer Research 73.9 (2013), pp. 2737–2742. [73] M. Binnewies, E. Roberts, K. Kersten, et al. “Understanding the tumor immune microenvironment (TIME) for effective therapy”. In: Nature Medicine 24.5 (2018), p. 541. 107 BIBLIOGRAPHY [74] T. Gajewski, H. Schreiber, and Y. Fu. “Innate and adaptive immune cells in the tumor microenvironment”. In: Nature Immunology 14.10 (2013), p. 1014. [75] U. Martinez-Outschoorn, R. Balliet, D. Rivadeneira, et al. “Oxidative stress in cancer associated fibroblasts drives tumor-stroma co-evolution: A new paradigm for understanding tumor metabolism, the field effect and genomic instability in cancer cells”. In: Cell Cycle 9.16 (2010), pp. 3276–3296. [76] Y. Rattigan, B. Patel, E. Ackerstaff, et al. “Lactate is a mediator of metabolic cooperation between stromal carcinoma associated fibroblasts and glycolytic tumor cells in the tumor microenvironment”. In: Experimental Cell Research 318.4 (2012), pp. 326–335. [77] C. Uyttenhove, L. Pilotte, I. Théate, et al. “Evidence for a tumoral immune resistance mechanism based on tryptophan degradation by indoleamine 2, 3-dioxygenase”. In: Nature Medicine 9.10 (2003), p. 1269. [78] F. Chen, X. Zhuang, L. Lin, et al. “New horizons in tumor microenvironment biology: challenges and opportunities”. In: BMC Medicine 13.1 (2015), p. 45. [79] S. Chirasani, P. Leukel, E. Gottfried, et al. “Diclofenac inhibits lactate formation and efficiently counteracts local immune suppression in a murine glioma model”. In: International Journal of Cancer 132.4 (2013), pp. 843–853. [80] S. Kouidhi, F. Ben Ayed, and A. Benammar Elgaaied. “Targeting tumor metabolism: a new challenge to improve immunotherapy”. In: Frontiers in Immunology 9 (2018), p. 353. [81] K. Patra, Q. Wang, P. Bhaskar, et al. “Hexokinase 2 is required for tumor initiation and maintenance and its systemic deletion is therapeutic in mouse models of cancer”. In: Cancer Cell 24.2 (2013), pp. 213–228. [82] M. Sukumar, J. Liu, Y. Ji, et al. “Inhibiting glycolytic metabolism enhances CD8+ T cell memory and antitumor function”. In: The Journal of Clinical Investigation 123.10 (2013), pp. 4479–4488. [83] Y. Xiang, Z. Stine, J. Xia, et al. “Targeted inhibition of tumor-specific glutaminase diminishes cellautonomous tumorigenesis”. In: The Journal of Clinical Investigation 125.6 (2015), pp. 2293– 2306. [84] T. Eleftheriadis, G. Pissas, A. Karioti, et al. “Dichloroacetate at therapeutic concentration alters glucose metabolism and induces regulatory T-cell differentiation in alloreactive human lymphocytes”. In: Journal of Basic and Clinical Physiology and Pharmacology 24.4 (2013), pp. 271–276. [85] V. Balachandran, M. Cavnar, S. Zeng, et al. “Imatinib potentiates antitumor T cell responses in gastrointestinal stromal tumor through the inhibition of IDO”. In: Nature Medicine 17.9 (2011), p. 1094. [86] U. Martinez-Outschoorn, M. Peiris-Pages, R. Pestell, et al. “Cancer metabolism: a therapeutic perspective”. In: Nature Reviews Clinical Oncology 14.1 (2017), p. 11. 108 BIBLIOGRAPHY [87] E. Klipp, R. Herwig, A. Kowald, et al. Systems biology in practice: concepts, implementation and application . John Wiley & Sons, 2008. [88] E. Stalidzans, A. Seiman, K. Peebo, et al. “Model-based metabolism design: constraints for kinetic and stoichiometric models”. In: Biochemical Society Transactions (2018), BST20170263. [89] F. Llaneras and J. Picó. “Stoichiometric modelling of cell metabolism”. In: Journal of Bioscience and Bioengineering 105.1 (2008), pp. 1–11. [90] S. Klamt and E. Gilles. “Minimal cut sets in biochemical reaction networks”. In: Bioinformatics 20.2 (2004), pp. 226–234. [91] J. Papin, N. Price, S. Wiback, et al. “Metabolic pathways in the post-genome era”. In: Trends in Biochemical Sciences 28.5 (2003), pp. 250–258. [92] C. Wagner and R. Urbanczik. “The geometry of the flux cone of a metabolic network”. In: Biophysical Journal 89.6 (2005), pp. 3837–3845. [93] F. Llaneras and J. Picó. “A procedure for the estimation over time of metabolic fluxes in scenarios where measurements are uncertain and/or insufficient”. In: BMC Bioinformatics 8.1 (2007), p. 421. [94] J. Orth, I. Thiele, and B. Palsson. “What is flux balance analysis?” In: Nature Biotechnology 28.3 (2010), p. 245. [95] N. Lewis, K. Hixson, T. Conrad, et al. “Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models”. In: Molecular Systems Biology 6.1 (2010), p. 390. [96] R. Mahadevan and C. Schilling. “The effects of alternate optimal solutions in constraint-based genome-scale metabolic models”. In: Metabolic Engineering 5.4 (2003), pp. 264–276. [97] Q. Beg, A. Vazquez, J. Ernst, et al. “Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity”. In: Proceedings of the National Academy of Sciences 104.31 (2007), pp. 12663–12668. [98] J. Park, T. Kim, and S. Lee. “Prediction of metabolic fluxes by incorporating genomic context and flux-converging pattern analyses”. In: Proceedings of the National Academy of Sciences 107.33 (2010), pp. 14931–14936. [99] J. Edwards, R. Ramakrishna, and B. Palsson. “Characterizing the metabolic phenotype: a phenotype phase plane analysis”. In: Biotechnology and Bioengineering 77.1 (2002), pp. 27–36. [100] D. Segrè, D. Vitkup, and G. Church. “Analysis of optimality in natural and perturbed metabolic networks”. In: Proceedings of the National Academy of Sciences 99.23 (2002), pp. 15112–15117. [101] T. Shlomi, O. Berkman, and E. Ruppin. “Regulatory on/off minimization of metabolic flux changes after genetic perturbations”. In: Proceedings of the National Academy of Sciences 102.21 (2005), pp. 7695–7700. 109 BIBLIOGRAPHY [102] N. Duarte, S. Becker, N. Jamshidi, et al. “Global reconstruction of the human metabolic network based on genomic and bibliomic data”. In: Proceedings of the National Academy of Sciences 104.6 (2007), pp. 1777–1782. [103] H. Ma, A. Sorokin, A. Mazein, et al. “The Edinburgh human metabolic network reconstruction and its functional analysis”. In: Molecular Systems Biology 3.1 (2007), p. 135. [104] T. Hao et al. “Compartmentalization of the Edinburgh human metabolic network”. In: BMC Bioinformatics 11.1 (2010), p. 393. [105] I. Thiele, N. Swainston, R. Fleming, et al. “A community-driven global reconstruction of human metabolism”. In: Nature Biotechnology 31.5 (2013), p. 419. [106] A. Mardinoglu, R. Agren, C. Kampf, et al. “Integration of clinical data with a genome-scale metabolic model of the human adipocyte”. In: Molecular Systems Biology 9.1 (2013), p. 649. [107] R. Agren, S. Bordel, A. Mardinoglu, et al. “Reconstruction of genome-scale active metabolic networks for 69 human cell types and 16 cancer types using INIT”. In: PLoS Computational Biology 8.5 (2012), e1002518. [108] P. Romero, J. Wagg, M. Green, et al. “Computational prediction of human metabolic pathways from the complete human genome”. In: Genome Biology 6.1 (2005), R2. [109] M. Kanehisa, M. Furumichi, M. Tanabe, et al. “KEGG: new perspectives on genomes, pathways, diseases and drugs”. In: Nucleic Acids Research 45.D1 (2016), pp. D353–D361. [110] A. Mardinoglu, R. Agren, C. Kampf, et al. “Genome-scale metabolic modelling of hepatocytes reveals serine deficiency in patients with non-alcoholic fatty liver disease”. In: Nature Communications 5 (2014), p. 3083. [111] N. Swainston, K. Smallbone, H. Hefzi, et al. “Recon 2.2: from reconstruction to model of human metabolism”. In: Metabolomics 12.7 (2016), p. 109. [112] E. Brunk, S. Sahoo, D. Zielinski, et al. “Recon3D enables a three-dimensional view of gene variation in human metabolism”. In: Nature Biotechnology 36.3 (2018), p. 272. [113] I. Thiele, S. Sahoo, A. Heinken, et al. “When metabolism meets physiology: Harvey and Harvetta”. In: Manuscript submitted for publication (). [114] J. L. Robinson et al. “An atlas of human metabolism”. In: Science Signaling 13.624 (2020). [115] M. Uhlén, L. Fagerberg, B. Hallström, et al. “Tissue-based map of the human proteome”. In: Science 347.6220 (2015), p. 1260419. [116] S. Becker and B. Palsson. “Context-specific metabolic networks are consistent with experiments”. In: PLoS Computational Biology 4.5 (2008), e1000082. [117] T. Shlomi, M. Cabili, M. Herrgård, et al. “Network-based prediction of human tissue-specific metabolism”. In: Nature Biotechnology 26.9 (2008), p. 1003. 110 BIBLIOGRAPHY [118] H. Zur, E. Ruppin, and T. Shlomi. “iMAT: an integrative metabolic analysis tool”. In: Bioinformatics 26.24 (2010), pp. 3140–3142. [119] R. Agren, A. Mardinoglu, A. Asplund, et al. “Identification of anticancer drugs for hepatocellular carcinoma through personalized genome-scale metabolic modeling”. In: Molecular Systems Biology 10.3 (2014), p. 721. [120] K. Yizhak et al. “Phenotype-based cell-specific metabolic modeling reveals metabolic liabilities of cancer”. In: Elife 3 (2014), e03641. [121] A. Schultz and A. Qutub. “Reconstruction of tissue-specific metabolic networks using CORDA”. In: PLoS Computational Biology 12.3 (2016), e1004808. [122] L. Jerby, T. Shlomi, and E. Ruppin. “Computational reconstruction of tissue-specific metabolic models: application to human liver metabolism”. In: Molecular Systems Biology 6.1 (2010), p. 401. [123] Y. Wang, J. Eddy, and N. Price. “Reconstruction of genome-scale metabolic models for 126 human tissues using mCADRE”. In: BMC Systems Biology 6.1 (2012), p. 153. [124] N. Vlassis, M. Pacheco, and T. Sauter. “Fast reconstruction of compact context-specific metabolic network models”. In: PLoS Computational Biology 10.1 (2014), e1003424. [125] E. Özcan and T. Çakır. “Reconstructed metabolic network models predict flux-level metabolic reprogramming in glioblastoma”. In: Frontiers in Neuroscience 10 (2016), p. 156. [126] E. Motamedian, G. Ghavami, and S. Sardari. “Investigation on metabolism of cisplatin resistant ovarian cancer using a genome scale metabolic model and microarray data”. In: Iranian Journal of Basic Medical Sciences 18.3 (2015), p. 267. [127] R. Metri et al. “Modelling metabolic rewiring during melanoma progression using flux balance analysis”. In: 2017 IEEE International Conference on Bioinformatics and Biomedicine (BIBM) . IEEE. 2017, pp. 134–137. [128] F. Shen et al. “Systematic investigation of metabolic reprogramming in different cancers based on tissue-specific metabolic models”. In: Journal of Bioinformatics and Computational Biology 14.05 (2016), p. 1644001. [129] F.-S. Wang et al. “Genome-Scale Metabolic Modeling with Protein Expressions of Normal and Cancerous Colorectal Tissues for Oncogene Inference”. In: Metabolites 10.1 (2020), p. 16. [130] B. ter Braak et al. “Insulin-like growth factor 1 receptor activation promotes mammary gland tumor development by increasing glycolysis and promoting biomass production”. In: Breast Cancer Research 19.1 (2017), p. 14. [131] P. Ghaffari et al. “Identifying anti-growth factors for human cancer cell lines through genome-scale metabolic modeling”. In: Scientific Reports 5.1 (2015), pp. 1–10. 111 BIBLIOGRAPHY [132] K. Yizhak et al. “A computational study of the Warburg effect identifies metabolic targets inhibiting cancer migration”. In: Molecular Systems Biology 10.8 (2014), p. 744. [133] S. Sahoo et al. “Metabolite systems profiling identifies exploitable weaknesses in retinoblastoma”. In: FEBS Letters 593.1 (2019), pp. 23–41. [134] R. Agren et al. “Identification of anticancer drugs for hepatocellular carcinoma through personalized genome-scale metabolic modeling”. In: Molecular Systems Biology 10.3 (2014), p. 721. [135] F. Gatto, R. Ferreira, and J. Nielsen. “Pan-cancer analysis of the metabolic reaction network”. In: Metabolic Engineering 57 (2020), pp. 51–62. [136] M. Uhlen et al. “A pathology atlas of the human cancer transcriptome”. In: Science 357.6352 (2017). [137] F. Han et al. “Genome-wide metabolic model to improve understanding of CD4+ T cell metabolism, immunometabolism and application in drug design”. In: Molecular BioSystems 12.2 (2016), pp. 431– 443. [138] B. L. Puniya et al. “Integrative computational approach identifies drug targets in CD4+ T-cellmediated immune disorders”. In: NPJ Systems Biology and Applications 7.1 (2021), pp. 1–18. [139] A. M. Abdel-Haleem et al. “Model-Driven Analysis of Gene Expression Data: Application to Metabolic Re-programming during T-cell Activation”. In: Proceedings of the International Conference on Bioinformatics & Computational Biology (BIOCOMP) . The Steering Committee of The World Congress in Computer Science. 2015, p. 91. [140] A. Bordbar et al. “Insight into human alveolar macrophage and M. tuberculosis interactions via metabolic reconstructions”. In: Molecular Systems Biology 6.1 (2010). [141] A. Bordbar et al. “Model-driven multi-omic data analysis elucidates metabolic immunomodulators of macrophage activation”. In: Molecular Systems Biology 8.1 (2012). [142] K. Yu and M. Snyder. “Omics profiling in precision oncology”. In: Molecular & Cellular Proteomics 15.8 (2016), pp. 2525–2536. [143] C. Manzoni, D. Kia, J. Vandrovcova, et al. “Genome, transcriptome and proteome: the rise of omics data and their integration in biomedical sciences”. In: Briefings in Bioinformatics 19.2 (2016), pp. 286–302. [144] R. Bumgarner. “Overview of DNA microarrays: types, applications, and their future”. In: Current Protocols in Molecular Biology 101.1 (2013), pp. 22–1. [145] J. Shendure. “The beginning of the end for microarrays?” In: Nature Methods 5.7 (2008), pp. 585– 587. [146] S. Liu and C. Trapnell. “Single-cell transcriptome sequencing: recent advances and remaining challenges”. In: F1000Research 5 (2016). 112 BIBLIOGRAPHY [147] L. Jerby and E. Ruppin. “Predicting drug targets and biomarkers of cancer via genome-scale metabolic modeling”. In: Clinical Cancer Research 18.20 (2012), pp. 5572–5584. [148] H. Nam et al. “A systems approach to predict oncometabolites via context-specific genome-scale metabolic networks”. In: PLoS Computational Biology 10.9 (2014), e1003837. [149] J. Gagneur and S. Klamt. “Computation of elementary modes: a unifying framework and the new binary approach”. In: BMC Bioinformatics 5.1 (2004), p. 175. [150] L. Jerby and E. Ruppin. “Predicting drug targets and biomarkers of cancer via genome-scale metabolic modeling”. In: Clinical Cancer Research 18.20 (2012), pp. 5572–5584. [151] O. Folger et al. “Predicting selective drug targets in cancer through metabolic networks”. In: Molecular Systems Biology 7.1 (2011), p. 501. [152] C. Mattiuzzi, F. Sanchis-Gomar, and G. Lippi. “Concise update on colorectal cancer epidemiology”. In: Annals of Translational Medicine 7.21 (2019). [153] P. Rawla, T. Sunkara, and A. Barsouk. “Epidemiology of colorectal cancer: incidence, mortality, survival, and risk factors”. In: Przeglad Gastroenterologiczny 14.2 (2019), p. 89. [154] I. A. for Research on Cancer. Cancer Today - Cancer Fact Sheets . 2020. url: ?iiT,ff;+QX B`+X7`fiQ/vf7+i@b?22ib@+M+2`b. [155] Y. Hao et al. “Integrated analysis of multimodal single-cell data”. In: Cell 184.13 (2021), pp. 3573– 3587. [156] J. Guinney et al. “The consensus molecular subtypes of colorectal cancer”. In: Nature Medicine 21.11 (2015), pp. 1350–1356. [157] Y. Hao et al. “Integrated analysis of multimodal single-cell data”. In: Cell 184.13 (2021), pp. 3573– 3587. [158] P. Hoffman. SeuratDisk: Interfaces for HDF5-Based Single Cell File Formats . 2021. url: ?iiTb, ff;Bi?m#X+QKfKQDp2xm`2fb2m`i@/BbF. [159] R. Satija et al. SeuratObject: Data Structures for Single Cell Data . R package version 4.0.2. 2021. url: ?iiTb,ff*_LX_@T`QD2+iXQ`;fT+F;24a2m`iP#D2+i. [160] L. Zappia and A. Oshlack. “Clustering trees: a visualization for evaluating clusterings at multiple resolutions”. In: GigaScience 7.7 (2018). doi: RyXRyNjf;B;b+B2M+2f;Bvy3j. url: ?iiTb, ff/QBXQ`;fRyXRyNjf;B;b+B2M+2f;Bvy3j. [161] M. Mircea. SIGMA: A clusterability measure for scRNA-seq data . R package version 0.0.0.1. 2021. [162] T. Tickle et al. inferCNV of the Trinity CTAT Project. Klarman Cell Observatory, Broad Institute of MIT and Harvard. Cambridge, MA, USA, 2019. url: ?iiTb,ff;Bi?m#X+QKf#`Q/BMbiBimi2f BM72`*Lo. 113 BIBLIOGRAPHY [163] P. W. Eide et al. “CMScaller: an R package for consensus molecular subtyping of colorectal cancer pre-clinical models”. In: Scientific Reports 7 (2017), p. 16618. doi: RyXRyj3fb9R8N3@yRd@R ed9d@t. [164] J. Qian et al. “A pan-cancer blueprint of the heterogeneous tumor microenvironment revealed by single-cell profiling”. In: Cell research 30.9 (2020), pp. 745–762. doi: RyXRyj3fb9R9kk@yky @yj88@y. [165] H. Lee et al. “Lineage-dependent gene expression programs influence the immune landscape of colorectal cancer”. In: Nature Genetics 52.6 (2020), pp. 594–603. doi: RyXRyj3fb9R833@yk y@yeje@x. [166] C. Smillie et al. “Intra-and inter-cellular rewiring of the human colon during ulcerative colitis”. In: Cell 178.3 (2019), pp. 714–730. doi: RyXRyRefDX+2HHXkyRNXyeXykN. [167] A. Loregger et al. “The E3 ligase RNF43 inhibits Wnt signaling downstream of mutated 𝛽-catenin by sequestering TCF4 to the nuclear membrane”. In: Science Signaling 8.393 (2015), ra90–ra90. [168] S Salahshor and J. Woodgett. “The links between axin and carcinogenesis”. In: Journal of Clinical Pathology 58.3 (2005), pp. 225–236. [169] S. Razak et al. “Screening and computational analysis of colorectal associated non-synonymous polymorphism in CTNNB1 gene in Pakistani population”. In: BMC Medical Genetics 20.1 (2019), pp. 1–12. [170] C. F. Rochlitz, R. Herrmann, and E. de Kant. “Overexpression and amplification of c-myc during progression of human colorectal cancer”. In: Oncology 53.6 (1996), pp. 448–454. [171] Z. Wang et al. “The prognostic and clinical value of CD44 in colorectal cancer: a meta-analysis”. In: Frontiers in Oncology 9 (2019), p. 309. [172] S.-M. Wang et al. “Clinical significance of MLH1/MSH2 for stage II/III sporadic colorectal cancer”. In: World Journal of Gastrointestinal Oncology 11.11 (2019), p. 1065. [173] J. L. Robinson et al. “An atlas of human metabolism”. In: Science Signaling 13.624 (2020), eaaz1482. [174] R. Agren et al. “Identification of anticancer drugs for hepatocellular carcinoma through personalized genome-scale metabolic modeling”. In: Molecular Systems Biology 10.3 (2014), p. 721. [175] A. Ebrahim et al. “COBRApy: constraints-based reconstruction and analysis for python”. In: BMC Systems Biology 7.1 (2013), pp. 1–6. [176] J. Ferreira et al. “Troppo-A Python framework for the reconstruction of context-specific metabolic models”. In: International Conference on Practical Applications of Computational Biology & Bioinformatics . Springer. 2019, pp. 146–153. 114 BIBLIOGRAPHY [177] N. Vlassis, M. P. Pacheco, and T. Sauter. “Fast reconstruction of compact context-specific metabolic network models”. In: PLoS Computational Biology 10.1 (2014), e1003424. [178] V. Vieira, J. Ferreira, and M. Rocha. “A pipeline for the reconstruction and evaluation of contextspecific human metabolic models at a large-scale”. In: PLOS Computational Biology 18.6 (2022), e1009294. [179] A. Richelle, C. Joshi, and N. E. Lewis. “Assessing key decisions for transcriptomic data integration in biochemical networks”. In: PLoS Computational Biology 15.7 (2019), e1007185. [180] D. Wishart, Y. Feunang, A. Marcu, et al. “HMDB 4.0: the human metabolome database for 2018”. In: Nucleic Acids Research 46.D1 (2017), pp. D608–D617. [181] M. K. Aurich et al. “Prediction of intracellular metabolic states from extracellular metabolomic data”. In: Metabolomics 11.3 (2015), pp. 603–619. [182] M. Mir et al. “Optical measurement of cycle-dependent cell growth”. In: Proceedings of the National Academy of Sciences 108.32 (2011), pp. 13124–13129. [183] M. Beck et al. “The quantitative proteome of a human cell line”. In: Molecular Systems Biology 7.1 (2011), p. 549. [184] E. H. Chapman, A. S. Kurec, and F. Davey. “Cell volumes of normal and malignant mononuclear cells.” In: Journal of Clinical Pathology 34.10 (1981), pp. 1083–1090. [185] Optimization of Human T Cell Expansion Protocol: Effects of Early Cell Dilution . Technical Bulletin. STEMCELL Technologies, 2020. [186] J. Gu et al. “Metabolomics analysis in serum from patients with colorectal polyp and colorectal cancer by 1H-NMR spectrometry”. In: Disease Markers 2019 (2019). [187] C. Zhang et al. “Metabolomic profiling identified serum metabolite biomarkers and related metabolic pathways of colorectal cancer”. In: Disease Markers 2021 (2021). [188] Y. Qiu et al. “Serum metabolite profiling of human colorectal cancer using GCTOFMS and UPLCQTOFMS”. In: Journal of Proteome Research 8.10 (2009), pp. 4844–4850. [189] S. Nishiumi et al. “A novel serum metabolomics-based diagnostic approach for colorectal cancer”. In: PloS One 7.7 (2012), e40459. [190] B. Tan et al. “Metabonomics identifies serum metabolite markers of colorectal cancer”. In: Journal of Proteome Research 12.6 (2013), pp. 3000–3009. [191] J. Zhu et al. “Colorectal cancer detection using targeted serum metabolic profiling”. In: Journal of Proteome Research 13.9 (2014), pp. 4120–4130. [192] M. Kuhn. caret: Classification and Regression Training . R package version 6.0-88. 2021. url: ?iiTb,ff*_LX_@T`QD2+iXQ`;fT+F;24+`2i. 115 APPENDIX A. SUPPLEMENTARY FIGURES Figure 52: Heatmap of CNV predictions for the tumour cells of patient KUL21. 122 A.1. CHAPTER 3 Figure 53: Heatmap of CNV predictions for the tumour cells of patient SMC04. 123 APPENDIX A. SUPPLEMENTARY FIGURES Figure 54: Heatmap of CNV predictions for the tumour cells of patient SMC07. 124 A.1. CHAPTER 3 Figure 55: Heatmap of CNV predictions for the tumour cells of patient SMC10. 125 APPENDIX A. SUPPLEMENTARY FIGURES A.2 Chapter 4 Figure 56: Similarity between (A) the structure of the models (i.e., reaction presence/absence) and (B) predicted fluxes under normal human blood medium. The smaller Euclidean distance is, the smaller the similarity is. 126 A.2. CHAPTER 4 Figure 57: Pathway coverage (%). 127 APPENDIX A. SUPPLEMENTARY FIGURES Figure 58: Biomass (A) and ATP (B) production when biomass was set as the only objective for all models. Figure 59: Cumulative fluxes (mmol/gDW/h) of the reactions that produce NADH, from all source pathways, for proliferative CD4 and CD8 T-cell models 128 A.2. CHAPTER 4 Figure 60: Distribution of (A) DNA and (B) RNA production of the T-cell types, with and without glutamine in the medium. 129 APPENDIX A. SUPPLEMENTARY FIGURES Figure 61: Distribution of DNA production of the T-cell types, with and without nucleotides in the medium. Figure 62: Distribution of the biomass flux of the T-cell types, with and without glucose in the medium. 130 A.2. CHAPTER 4 Figure 63: (A) Number of models where biomass flux increases, decreases or suffers no change. This information is further showed by (B) tissue of origin, (C) CMS type, and (D) cell-type. 131