Functional effects of ccRCC are regulatory layer dependent

ccRCC is a heterogeneous cancer with transformation and progression linked to widespread dysregulation of biological processes. To understand how cellular processes are affected, analyses must extend across regulatory layers and incorporate multiple omics’ data types. Recently, Clark et al. and the CPTAC consortium investigated the impact of genomic alterations on protein regulation and created revised subtypes based on integrated analyses [14]. As part of their study, Clark et al. generated protein, RNA and DNA methylation measurements from 110 patients with ccRCC [14], of which we filter to 76 cases, see “Methods” for details.

We sought to investigate the compounding effects of dysregulation by following gene changes from the DNA methylation via the mRNA to protein layers to unravel the regulatory layer of each gene within ccRCC. The original study by Clark et al. did not include normal methylation samples thus our study extends their work by also considering the effect of epigenetic changes in the context of ccRCC (see Table S1 for sample annotations and “Methods” for quality controls). The CPTAC cohort is demographically homogenous but includes cases across four stages of ccRCC and multiple mutational patterns (Table S1, Additional file 1: Fig. S1a). Data types are herein denoted “data layer” as they correspond to distinct layers of regulation (DNA methylation beta values, normalised mRNA expression and normalised protein expression). PCA emphasises the primary source of variation in each data layer is the sample type (tumour versus normal) with neither of the top two components explaining tumour stage (Additional file 1: Fig. S1a).

We identified genes changing significantly between tumour and normal by performing differential analysis tests for each data layer and found widespread changes, most notably on the mRNA layer (“Methods”, Table S1). Whilst there is minimal negative correlation between gene changes on the DNA methylation and mRNA layers, the mRNA and protein layers are correlated (Fig. 1a). We found 1284 significant genes shared across the mRNA and protein layers and 796 significant genes between all three layers, where significant is defined as p < 0.05, |beta_diff|> 0.1, |protein_logFC|> 0.5 and |rna_logFC|> 1.0.

Fig. 1

SiRCle model logical regulatory rules. a Spearman correlation comparing the change between tumour and normal samples for mRNA expression and DNA methylation, and mRNA expression and protein expression. b Numbers of enriched Gene Ontology (GO) terms when comparing tumour versus normal (TvN) for each data type, indicating shared terms by overlapping regions in the Venn diagram. “Down regulated” corresponds to GO terms enriched in genes that are decreased in tumour samples, whilst “Up regulated” is for GO terms enriched in genes increased in tumour samples. c Alluvial plot depicting the flows used to define SiRCle clusters. The plot is read from the left to the right, in line with the flow of information in a biological system. Each data type has been labelled as a layer, with one of three states defined for each layer based on the results for differential analysis between tumour versus normal in that data type. The “No Change” on the translational layer is further subdivided to reflect protein detection and threshold criteria met. d Explanation of the directions of change and the corresponding biological meaning for each layer. e Summary of collapse of regulatory groupings coloured by regulatory grouping 2 which is used throughout the paper. f How the ccRCC data is dispersed into the different regulatory clusters using the SiRCle grouping method with the biological outcome described on the side

Given the small number of significant genes shared across layers, yet comprehensive changes within each layer, we sought to determine whether the affected genes shared biological function. To do this, we performed over representation analysis (ORA) of Gene Ontology (GO) terms on significant genes from each differential analysis, sub-setting genes by direction of change. We observed both unique and shared biological functions between the layers, with more similar functions enriched across layers for upregulated than downregulated genes (Fig. 1b). Hypermethylated genes were most enriched for terms associated with development (methylation layer), whilst supressed transcripts were most associated with transporter activity (mRNA layer) and supressed proteins with mitochondrial processes (protein layer) (Additional file 1: Fig. S1b–c, Table S2). In line with the findings from the original study [14], immune response terms showed the most significant enrichment across all data layers for genes upregulated on the protein and mRNA layers (Additional file 1: Fig. S1b–c, Table S2).

Following the flow of information extracts clusters that drive specific cellular phenotypes in ccRCC patients

Given the lack of correlation between the data layers and the heterogeneity of functional enrichment, we posited that grouping genes by their pattern across layers, prior to performing ORA would facilitate biological interpretation. Based on the differential analysis, we defined three states comparing tumour and normal for DNA methylation, mRNA expression and protein abundance, namely positive, negative or unchanged, for each gene and layer (Fig. 1c, Table 1, “Methods” section). Given that the protein abundance is the data layer closest to the phenotype, yet at the same time has the lowest coverage we aim to account for this. Hence, “No Change”, which includes proteins that do not meet the significance and/or the fold change threshold, is subdivided into four states for the protein layer: 1. “undetected” (genes detected in other data layers, but not in the proteomics data layer), 2. “not significant” (proteins that do not meet the significance threshold), 3. “significant negative” and 4. “significant positive” (proteins that do not meet the negative/positive fold change threshold but meet the significance threshold) (Fig. 1c, labelled, translational layer). Since this is an ordered series of three 3-state-3-state-6-state transitions between the layers, there are 54 possible “flows”. Each gene will be assigned to one of those 54 “flows”, which is required to enable the summary of genes into clusters.

To perform biologically meaningful clustering, we group genes based on the layer where the gene “changed” in the tumour sample compared to the normal sample and termed this “Regulation Grouping” (RG). We do this by following the central dogma of biology [61], i.e. genetic information proceeds from DNA, to mRNA, then to protein (Fig. 1c–d). For example, if a gene is hypermethylated, has a decrease in mRNA expression and displays a decrease in protein expression, we can likely conclude that the dysregulation first occurred on the DNA methylation layer, meaning this gene is supressed via Methylation-Driven Suppression (MDS). Whilst there are 54 possible flows, we use the RG to summarise these 54 flows into SiRCle clusters based on combinations of regulations between DNA methylation and mRNA expression (MR) and/or between mRNA expression and protein expression (RP) (Fig. 1d–e), ultimately defining the SiRCle cluster names that reflect those regulations (“Methods”, Table 1).

We term the process “Signature Regulatory Clustering” (SiRCle), and throughout the paper, we use “RG2_Change” with ten clusters (Fig. 1f). As an input for SiRCle clustering, we include genes that are either detected on the protein layer or on both the mRNA and DNA methylation layers (background P¦R&M). For more information on thresholds, RG and background see “Methods” and the package guidelines.

Next, we performed ORA on the SiRCle clusters (Fig. 1f) and found that each cluster was enriched for a biological signature important in ccRCC and consolidated the functional enrichment results we observed across each layer when independently analysed (Fig. 2a, Table S2). We found clear biological signatures for both SiRCle clusters regulated on one layer and SiRCle clusters regulated on multiple layers. Specifically, hypomethylation (Methylation-Driven Enhancement (MDE)) was associated with the hypoxic response, oxidative stress response and angiogenesis (Fig. 2a), key players in ccRCC tumours [46]. These terms were not revealed in the top terms when we performed ORA on each layer independently (Additional file 1: Fig. S1b–c).

Fig. 2

Phenotypic changes in ccRCC and TF drivers. a Emapplot visualisation of the over representation analysis (ORA) performed on each SiRCle cluster resulted in biological pathways that are altered comparing tumour versus normal. Pathways were plotted if p-adjusted value (p.adj) < 0.05 and the gene ratio was greater than 5%. The dot size corresponds to the number of genes found in the cluster that are part of the biological pathway. The colour of the dot shows the p.adj of the ORA. The connecting lines (grey) show that the biological pathways have genes in common. b Transcription factor (TF) network of TFs from manually curated repositories (DoRothEA) that may drive genes in the cluster Methylation Driven Enhancement (MDE). The dot size corresponds to the number of genes targeted by the TF. The colour of the dot shows the p.adj. The connecting lines (grey) show the number of common genes the connected TFs regulate. c Wordcloud including the TF targets with shared TF binding in the MDE cluster, with the size corresponding to how many different TFs are predicted to regulate the gene’s expression. d HIF1A ChIP-seq peaks (RCC4, GSM3417827) binding at the transcription start site of VEGFA, SLC2A1 and CA9

In line with previous findings [14, 46, 62], we observed an increase in immune response pathways and show these genes are likely regulated by DNA hypomethylation (MDE, Fig. 2a) or enhanced transcription (Transcription and Processing-Driven Enhancement (TPDE), Fig. 2a). Clark et al. described a high correlation between mRNA and protein expression of immune signatures [14], supporting our finding that the regulation occurs at the transcriptional layer. We were also able to distinguish between immune response genes likely regulated by hypomethylation from those likely regulated by enhanced transcription (Fig. 2a).

We observed that genes in two SiRCle clusters are involved in distinct metabolic rewiring, with metabolic processes such as lipid and amino acid metabolism likely regulated by transcriptional suppression (Transcriptional and Processing-Driven Suppression (TPDS)), whilst mitochondrial respiration and nucleotide metabolism are likely downregulated by translational suppression of the proteins (Translation and translational-modification-Driven Suppression (TMDS), Fig. 2a). Genes regulated by translational suppression appear to be involved in mitochondrial metabolism, mitochondrial translation and mitochondrial morphology (TMDS, Fig. 2a). Hypermethylated genes and/or translationally enhanced genes were enriched for kidney development, hinting towards a loss of cellular identity in cancer cells (MDS_TMDE, Fig. 2a). Using SiRCle clustering, we were able to determine the layer at which a gene’s dysregulation occurs and find clusters corresponding to distinct biological processes that may underpin tumour pathology.

Transcription factors as drivers of genes in SiRCle clusters

A subset of SiRCle clusters contain genes that change their state at the transcriptional (TPDE/TPDS) and/or the methylation layer (MDE/MDS), suggesting that they are regulated by transcription or epigenetic factors. Indeed, a specific transcription factor (TF) may drive such changes in gene expression as TFs can act to enhance or repress gene transcription and DNA methylation changes can alter the TF’s ability to bind to the DNA [63]. We used validated TF-to-target interactions from DoRothEA [55] to identify which TFs were statistically associated with each SiRCle cluster. We use Fisher’s exact test (FET) to measure association between a TF and SiRCle clusters by testing if a TF targets more genes in a single cluster relative to the background frequency (all SiRCle clusters). This analysis only recovered TFs significantly associated with SiRCle clusters that are either up- or downregulated at the methylation and transcriptional layers (Fig. 2b–c, Table S3). The stratification provided via SiRCle clustering helps distil the data from which regulatory logic can be inferred.

Given that HIF1 drives angiogenesis in ccRCC tumours [46], we were unsurprised to find HIF1 TFs target genes in MDE as this cluster was enriched for angiogenesis related GO terms (Fig. 2a–b, Additional file 1: Fig. S2a, Table S3). We found that HIF1A had a significant increase in target-mRNAs driven by greater accessibility to binding due to DNA hypomethylation, which also translates to a detectable change on the protein layer (MDE). Given DoRothEA is not tissue specific, we sought to confirm that the HIF TFs bind to the target genes in ccRCC by using ChIP-seq data from HIF1A and EPAS1 (also known as HIF2A) in kidney cancer cell lines (see “Methods”). We found evidence of HIF1A binding at the transcription start site (TSS) of genes in MDE (VEGFA, CA9 and SLC2A1) (Fig. 2d) including many metabolic enzymes (HK2, ALDOC, GAPDH, PKM, ENO1, ENO2 and LDHA, Additional file 1: Fig. S2a). Our findings support that TF analyses coupled with SiRCle can elucidate relationships at both the transcriptional and methylation layers enabling the identification of novel TF targets (Fig. 2b, Additional file 1: Fig. S2c), which can be investigated for changes in downstream products and/or be followed up experimentally.

ccRCC characteristic metabolic changes are detected on distinct layers

Metabolic rewiring plays a crucial role in ccRCC and hence we sought to investigate the dysregulation of metabolic genes to understand the layer at which the metabolic pathways are altered. We performed gene set enrichment analysis (GSEA) using metabolic signatures from Gaude and Frezza [60], herein referred to as metabolic signatures, to identify pathways with coordinated protein changes in tumour versus normal. The protein layer was chosen for GSEA to capture changes in pathways’ enzyme activity, and then coupled with SiRCle cluster annotations to relate gene changes to the initial layer of dysregulation (Table S4). We found the majority of enzymes involved in glycolysis (~ 70%) were increased comparing tumours versus normal, fitting with the previously established upregulation of anaerobic glycolysis in ccRCC [8]. Interestingly we find that 47% of glycolytic enzymes are hypomethylated and upregulated on both the mRNA and protein layers (SiRCle cluster MDE, Fig. 3a and Additional file 1: Fig. S2b–c). From the TF analysis, we noticed that HIF1A targets in MDE included several glycolytic enzymes, and we also found HIF1A binding at the transcription start site of glycolytic genes in kidney cancer cell line ChIP-seq data (Additional file 1: Fig. S2c).

Fig. 3

Metabolic changes in tumour samples are regulated on distinct data landscapes. a Metabolic changes of glycolysis, tricarboxylic acid (TCA) cycle and serine, glycine and cysteine biosynthesis in ccRCC comparing tumour versus normal. Metabolic enzymes are written in bold and if they have been detected, they are labelled with coloured circles depending to the SiRCle cluster they are part of. Metabolite changes are labelled with plus for a positive and minus for a negative log2FC. The three main metabolic sections, categorized based on the SiRCle cluster to which the majority of metabolic enzymes belong, are summarised in pie charts. b–d Volcano plots are based on the metabolic pathways defined by Gaude and Frezza [60] and the number of proteins detected within the pathway is reported. Proteins that are unique for the metabolic pathways are displayed as circles and proteins that are part of multiple metabolic pathways are displayed as diamonds. The colour code is based on the SiRCle cluster the protein is part of and is summarised in the pie chart. GSEA was performed using the protein statistic. Here, we plot oxidative phosphorylation (b), valine, leucine and isoleucine metabolism (c) and tryptophan metabolism (d), all reported a GSEA p.adj of 0.000966 (the lowest p.adj by threshold limits)

We next sought to dissect the role of mitochondrial dysfunction in ccRCC, including the suppression of mitochondrial electron transport chain (ETC), since this has been previously reported [8] to play a role in ccRCC, yet how metabolic enzymes are regulated remains unclear. We found that enzymes in electron transport chain complexes were depleted in the tumour samples and were likely regulated at the translational layer (SiRCle cluster TMDS, Fig. 3b). ETC includes genes encoded by the nuclear genome regulated on the translational layer (TMDS), whilst genes encoded in the mitochondrial DNA were found in the transcriptional regulation cluster (TPDS_TMDE). Yet, given that we do not include any normalisation for mitochondrial copy number, this could change once the mitochondrial copy number is taken into account. Hence, we propose that decreased genomic translation leads to instability of the complexes and protein degradation and vice versa. The translation-driven downregulation affects almost all enzymes of the TCA cycle (SiRCle cluster TMDS, Fig. 3a and Additional file 1: Fig. S2d). In addition to these primary metabolic rewiring steps, there are many secondary pathways altered in ccRCC [8]. By coupling SiRCle annotations with GSEA, we found that the suppression of cysteine methionine and glutathione (GSH) metabolism [7] occur on the transcriptional layer (SiRCle cluster TPDS, Fig. 3a and Additional file 1: Fig. S2e, f). Similarly, serine biosynthesis, which is important to fuel GSH biosynthesis, also occurs on the transcriptional layer. In accordance with the ORA results (Fig. 2a), we found that other amino acids’ metabolic pathways such as tryptophan, valine, leucine and isoleucine are downregulated either by translational or transcriptional suppression (SiRCle clusters TPDS and TMDS, Fig. 3c–d).

To understand the impact of the enzyme expression on the metabolic profile, we used metabolomics data from 84 ccRCC patients with paired tumour and adjacent tissue samples [7] and performed differential metabolomics analysis to identify the metabolites significantly changing between tumour and adjacent tissue. To consolidate the metabolite levels with enzyme information, we assigned significant metabolite changes to the pathways identified from the protein layer (Fig. 3a). We observed a “split” of glycolysis as previously described by Hakimi et al. [7], whereby metabolites upstream of GADPH are accumulated, and the downstream metabolites are depleted, despite the majority of glycolytic enzymes being upregulated on the protein layer (Fig. 3a, Additional file 1: Fig. S2b, Table S4). A similar split is observed in the TCA cycle, where citrate and succinate are increased despite a decrease in the enzyme’s protein expression, whilst fumarate and malate are depleted in line with the enzyme’s decreased protein expression (Fig. 3a). Overall, we show that SiRCle clustering can elucidate the layer where metabolic changes in ccRCC are orchestrated, highlighting points of intervention for targeting metabolism as anticancer strategy.

Functional differences between late and early-stage tumours show limited agreement across layers

The previous sections explored dysregulation in tumour samples across a heterogeneous cohort of ccRCC patients. However, it has been reported that tumour stage is an important determinant of prognosis and treatment response [45]. We thus sought to identify the regulatory differences between patients with early-stage (stage I, N = 30 and stage 2, N = 8) and late-stage (stage III, N = 27 and stage IV, N = 11) cancer by performing differential analysis on each layer, between tumour versus normal, for each patient group independently (Table S5). Using the significant genes from each analysis, we performed ORA on each data layer (Fig. 4a-b). On the protein and mRNA layers, we found that similar biological terms were enriched across the two patient’s groups, whilst on the methylation layer there were far more terms enriched in the late-stage samples (271) (Fig. 4c–d). The main differences between the top associated GO terms were an association with the regulation of lymphocytes (hypomethylated) and DNA-binding transcription activator activity (hypermethylated) in late-stage samples (Additional file 1: Fig. S3a).

Fig. 4

Late vs early tumours data integration. a Overlap of GO terms from ORA performed on each layer separately using the significant GO terms from differential analysis comparing late-stage and early-stage tumours. b As in a, except the gene overlap not the GO overlap. c Overlap of GO terms from ORA on DNA methylation for tumour vs normal on early-stage and late-stage datasets independently. d As in c except using the protein layer. e Overview of the SiRCle approach from the pre-processing (up to user), and SiRCle clustering, to the SiRCle integration for post hoc analyses

It is unsurprising that the tumour versus normal comparisons yielded similar results despite being run on patients with different stage tumours as the primary source of variation in the dataset is accounted for by the sample type. Hence, we next tested for differences in tumour profiles between late (stage III and IV) and early-stage (stage I and II) tumour samples (Table S5). When comparing late versus early-stage tumours, we found few changes on the protein layer, but aberrant changes on the DNA methylation and mRNA layers (Fig. 4a–b, Additional file 1: Fig. S3b). As tumours progress, the dynamic adaptability of metabolism plays a crucial role to ensure cell growth [64], and therefore, we tested for coordinated changes across metabolic signatures using the differences between late and early-stage tumours. We found minimal shared enrichment (Pcorrected < 0.05) of pathways between layers, with only protein modification enriched on all three layers (Table S5). Whilst we find significantly changing metabolic pathways within each layer, functional information is not readily captured across layers rendering the results challenging to interpret. However, we note that given the hierarchy of GO terms, there is a possibility that the layers have higher shared function via overlapping parental terms not captured by a direct comparison.

Integrated statistical test to identify changes between patient cohorts

Given the independently performed differential analyses were unable to extract shared functions when comparing patients with late versus early-stage tumours, we posited that integrating across the data layers prior to performing differential analysis and biological enrichment may better capture biological signal. We opted to use a VAE to learn gene-wise relationships across the three data layers for each patient, resulting in an integrated value for each patient enabling analysis between patient groups on, e.g. stage difference (Fig. 4e). Akin to matrix factorisation methods, a VAE finds a projection of data that centres on variance. However, unlike linear methods, such as PCA, a VAE does this using a neural network and can thus capture non-linear relationships.

In short, there are two parts of a VAE, (1) an encoding function that transforms data to a compressed representation and (2) a decoding function that recreates the input from the compressed representation. The encoding and decoding functions are shared once the parameters have been learnt from the training data, whilst an encoded value is specific to a given data point. For example, the compressed representation of data point \(_\) is given by \(_ \cong q}_\left(z|_\right)\cong _+_\odot _\), where \(_\) is the latent representation produced by the model conditioned on \(_\), \(_\) is the encoding function, \(_\) is stochastic noise and part of regularisation, \(_\) is the encoded mean and \(_\) is the encoded variance. Note, \(_\) refers to a feature vector of a gene, denoted by index \(i\) for a specific patient \(p\). The feature vector is the normalised values across the data layers (Fig. 4e).

The encoding and decoding functions are “learnt” by minimising an objective lower bound over a dataset, by iteratively updating parameters during training. For our purposes, the objective is to minimise the difference between input and output as per mean squared error and maximise the similarity between the latent space and a Gaussian normal distribution via maximum mean discrepancy (see “Methods” for specifics).

Given the relatively low signal produced by stage, extracting regulatory variations may be obstructed by the noise in the dataset when considering all regulatory relationships at once. As such, instead of learning a dataset-wide representation, we used genes from each SiRCle cluster to learn a representation that defines a given regulatory flow. As such, there is an encoding function \(_\left(z|x\right)\) for each SiRCle cluster.

As we can now calculate an integrated (encoded) value for each gene, for each patient, we can define a gene’s integrated difference as a gene’s mean encoding difference between two sets of patients, e.g. patients with late versus early tumours (Fig. 4e).

Mathematically, we define this below, where z refers to the latent encoding given a particular data point, which corresponds to a patient’s gene value across the layers. Gene index \(i\) is held constant as the mean difference is calculated between patient sets, for example patients with a late tumour may be in set S, and those with an early tumour are in set G.

$$D_||_)=\frac_^_-\frac_^_= }_- }_$$

For biological interpretation, we can correlate the integrated difference to a biological reference such as the mean difference between the two patient sets on the protein layer. For example, in the MDS SiRCle cluster, we see that the late versus early stage integrated difference correlates to the difference on the DNA methylation, protein and mRNA layers between these patient sets, whilst TPDE does not seem to be affected by DNA methylation, as expected (Additional file 1: Fig. S4a).

Using our integrated value for each gene, we next perform a Mann–Whitney U test to identify genes with a significant integrated difference between patients with early-stage (stage I, N = 30 and stage 2, N = 8) and late-stage (stage III, N = 27 and stage IV, N = 11) tumours. In the following sections, we demonstrate that VAE integration prior to performing analyses enables us to identify changes across the regulatory layers that were not found when performing differential analyses independently on each data type. Moreover, by benchmarking the VAE’s capacity to capture biological information on the integrated dimensions to six other cancer integration methods (intNMF [33], iCluster [34], JIVE [35], MCIA [36], RGCCA [37], tICA [38]) using the integration benchmarking package MOMIX [17], we find that our VAE approach performs comparatively, if not better at extracting biologically relevant orderings of genes, than existing methods (Additional file 1: Fig. S4b).

Metabolic alterations in tumour stage impact ccRCC’s metabolic fingerprint

Whilst there were few shared metabolic pathways in the layer-specific analyses, an integrative model may improve our ability to understand metabolic differences between ccRCC patient groups. As metabolic profiles alter during tumour progression [64], we tested for integrated differences between patients with late and early tumours by performing GSEA to identify pathways with coordinated integrated differences between the two groups (Fig. 4e). We define regulation as “up” if the integrated value increases and as “down” if the integrated value decreases. The integrated value is positively correlated with the mRNA or protein layers, as such it carries biological meaning however, “up” regulation of the integrated value of a specific gene does not equal the same increase in mRNA or protein (Additional file 1: Fig. S4a). In line with previous observations [14], we find oxidative phosphorylation was upregulated when comparing late with early-stage (Fig. 5a, Table S6), yet we are now additionally able to understand on what layer the enzyme expression is defined. We observed that oxidative phosphorylation enzymes are likely altered on the translational layer. Reassuringly, oxidative phosphorylation was also a significant pathway in the late versus early-stage tumour protein layer GSEA analysis, but not on the mRNA layer, matching our TMDS label and showing that trivial single-layer relationships are captured by the VAE approach (Table S6). The coordinated change observed in oxidative phosphorylation was also observed for enzymes of the mitochondrial TCA cycle, with some significantly (p.val < 0.05) upregulated enzymes such as IDH1-3, FH and SDHB in late-stage compared to early-stage (Fig. 5b). To understand if this is in accordance with the metabolite levels, we use metabolomics data comparing late (III and IV) with early-stage (I and II) patients published by Hakimi et al. [7]. Interestingly, citrate and cis-aconitate metabolite levels were depleted in late-stage compared to early-stage patients, whilst the downstream metabolites remained unchanged over the stages and have decreased protein expression in the tumour tissue (Fig. 5c).

Fig. 5

Metabolic alterations in tumour stage impact ccRCC’s metabolic fingerprint. Metabolic pathways based on gene expression are defined by Gaude and Frezza [60] with p.adj values corresponding to the GSEA results after ranking the genes in each SiRCle cluster using the VAE integrated rank. Metabolite pathways based on metabolites are defined by Hakimi et al. [7]. a Comparison of the VAE integrated rank of late and early-stage patients for oxidative phosphorylation for genes within the Translation and post-transcriptional Modification Driven Suppression (TMDS) SiRCle cluster. b Comparison of the VAE integrated rank of late with early patients for TCA cycle genes colour coded for the SiRCle clusters the gene belongs to. c Comparison of the \(}_FC\) of tumour versus normal (TvN) of late stage (III and IV) with early stage (I and II) patients for TCA cycle metabolites. d Comparison of the VAE integrated value of late with early patients for genes corresponding to the methionine metabolism pathway colour coded for the SiRCle clusters the gene belongs to. e Comparison of the \(}_FC\) (TvN) of late stage (III and IV) with early stage (I and II) patients for “cysteine, methionine SAM, taurine met