Systems level profiling of chemotherapy-induced stress resolution in cancer cells reveals druggable trade-offs

Significance Cancer therapies often fail to cure patients because a proportion of tumor cells withstand the toxic effects of chemotherapy. How surviving cancer cells recover from sublethal drug-induced stress is not known, but given that cellular resources are finite, stress resolution may come at the expense of less essential systems. Here, we studied the global cellular events of stress buildup and resolution in the bone marrow cancer, multiple myeloma, after proteasome inhibition, a commonly used therapeutic approach. Using a temporal multiomics approach, we delineate the unexpectedly complex and protracted changes myeloma cells undergo during stress resolution and demonstrate that recovering cells are more vulnerable to specific insults than acutely stressed cells. Thus, the findings may provide avenues for optimizing cancer therapies.

Markov time specifies a desired resolution of the graph clustering. The heatmap shows the VI between different clusterings at different Markov times. Blocks of dark areas mean that clusterings are similar and indicate a robust clustering. The same number of clustering across a long Markov time is also a signal of a robust clustering. Bottom: Biological Homogeneity Index (BHI) z-scores of the clusterings at different Markov times. The scores indicate how much biological information in the Gene Ontology (GO) is captured by the clustering. From the results shown in both panels, the optimal number of clusters to be applied was determined to be 6. B) and C) Enrichment of KEGG (B) and GO (C) terms in the 6 transcriptional clusters.

Figure S3. Proteasome inhibitor-induced metabolic pathway changes in RPMI-8226 cells.
A) Heatmap representing relative levels of 19S and 20S proteasome subunit transcripts determined by RNA sequencing (data shown are the mean of five independent experiments relative to day 0). B) GSVA of the indicated pathways from independent databases ("Proteasome pathway", Biocarta; "Proteasome", KEGG; "Ubiquitin mediated proteolysis", KEGG). C) Distribution of NRF2 target genes in
Relative expression (fold change) or late (day 4) after a Cfz pulse. B) GO biological process term enrichment for differentially expressed genes in RPMI-8226 cells treated with GCN2iB (1 µM, 48 h) in the absence of Cfz and early (day 2) or late (day 4) after a Cfz pulse. Only top 20 enriched terms are shown. Cut-off lines drawn at equivalent of P = 0.05 (dotted), P = 0.01 (dashed), P = 0.001 (solid). C) Number of deregulated genes in the 6 most significantly enriched GO Biological Process terms that did not show substantial overlap with other highly ranked terms. D) Functional subcategorisation of the 60 deregulated genes included in the GO term "Carboxylic acid metabolic process" (GO: 0019752), E) Functional subcategorisation of the deregulated genes included in the GO term "Organic cyclic compound metabolic process" (GO: 1901360).

Figure S7. Differential gene expression in GCN2-dependent and GCN2-independent cancer cell lines.
Heatmap showing transcriptome differences between GCN2-dependent and -independent CCLE cell lines, based on DepMap CERES scores for GCN2 (EIF2AK4) dependency. Differential expression analysis was conducted with DESeq2 with "tissue" as covariate. Normalised counts were transformed via variance stabilisation.

Supp Figure 11
Expression Z-score Proteasome inhibition triggers acute proteotoxic stress, which entails a drop in intracellular glucose levels and increased pyruvate availability compatible with heightened glycolytic activity. Tricarboxylic acid (TCA) cycle intermediates and mitochondrial proteins including electron transport chain (ETC) components decrease, while reduced protein degradation leads to amino acid shortages including glutamine scarcity. During stress recovery, glucose uptake becomes yet more suppressed in parallel with reduced GLUT1 expression. Persistent amino acid scarcity triggers GCN2 signalling to regulate fatty acid (FA) metabolism, revealing cellular dependency on GCN2, while mitochondrial respiration and ETC proteins remain suppressed.

Cellular recovery from proteasome inhibition.
RPMI-8226 myeloma cells were expanded to a total of 1.6 x 10 7 cells in 175 cm 2 tissue culture flasks and split into 16 subcultures on day -2. On day 0, cell density was adjusted to 1 x 10 6 cells/mL with fresh medium, and cells were allowed to equilibrate. Carfilzomib (750 nM) was then added to cultures for 1 h. Cells were then washed with PBS and cultured in fresh medium without carfilzomib. Medium was replaced every 48 h and cell density adjusted to 1x10 6 cells/mL on days 6 and 8. Samples for RNA sequencing (n = 5 subcultures), liquid chromatography-tandem mass spectroscopy (n = 3), tandem mass tag labelling proteomics (n = 2), nuclear magnetic resonance spectroscopy of cell culture supernatants (n = 3), quantitative gene expression analysis by qRT-PCR, immunoblotting, and assessment of viability by Trypan Blue exclusion (n = 4) and Guava® easyCyte TM quantification (n = 16) were collected prior to carfilzomib treatment on day 0, after 24 h (day 1), and immediately before changing culture medium on days 2, 4, 6, 8, and 10. Recovery experiments with other cell lines were performed concordantly, with samples collected as indicated. For recovery experiments including GCN2iB, cell density was adjusted to 1x10 6 cells/mL in fresh medium on days 2 and 4.

Cell viability assays.
Cells were mixed 1:1 (volume) with Trypan Blue (93595, Sigma-Aldrich) and then counted using a haematocytometer. In addition, 10 μL of cell culture was mixed with 190 μL ViaCount assay reagent (4000-0040, Merck-Millipore) and incubated for 5 min at RT. Viability data were acquired on a Guava® easyCyte TM instrument (Merck-Millipore) using GuavaSoft 3.  Fisher) according to the manufacturer instruction, using a BD LSRFortessa™ (BD Biosciences). Cell cycle was analysed using FITC BrdU Flow Kit (559619, BD Pharmingen TM ). RPMI-8226 cells were incubated for 2 h with 10 µM BrdU. Cells were stained for BrdU and 7-ADD according to the manufacturer's instructions, and analysed using a BD LSRFortessa™ (BD Biosciences). Analysis was performed using FlowJo (BD Pharmingen).

Quantitative real-time polymerase chain reaction (qRT-PCR).
Total RNA was extracted with the ReliaPrep TM RNA Miniprep Systems (Z6012, Promega), cDNA was synthesized by GoScript Reverse Transcription System (A5001, Promega). qRT-PCR was performed using PowerUp SYBR Green PCR Master Mix (A25742, Thermo Fisher Scientific) and analysed by StepOne TM Software v2.3 (Applied Biosystems, Thermo Fisher Scientific). mRNA transcript levels were quantified using the standard curve method, GAPDH was used as an internal control for normalization. Genespecific primers are listed in Supplementary Table 14.

Protein extraction and immunoblot analysis.
Whole cell extracts were prepared by lysing cells in RIPA lysis buffer (89900, Thermo Fisher Scientific), supplemented with protease with phosphatase inhibitors (05 892 970 001 and 04 906 845 001 respectively, Sigma-Aldrich), on ice for 15 min. Lysates were centrifuged at 13000 × g for 10 min at 4 °C and the supernatant was collected and stored at -80 °C. Protein concentration was determined using the Pierce BCA protein assay kit (23225, Thermo Fisher Scientific). Electrophoresis and immunoblot analyses were performed using NuPAGE protein analysis system (Thermo Fisher Scientific). Membranes were blocked with 5% skimmed milk in TBS plus 0.1% TWEEN® 20 (TBS-T) for 1 h at RT and incubated with primary antibodies (diluted 1: 1,000) overnight at 4 °C. Membranes were then washed with TBS-T and incubated with the appropriate secondary antibody (1:2500). Bands were visualised using enhanced chemiluminescence (ECL) reagents (RPN2232, Amersham Pharmacia) with a ChemiDoc image system with Image Lab software, version 5.

RNA-sequencing.
Total RNA from 2 x 10 6 cells per sample was extracted with the ReliaPrep TM RNA Miniprep Systems kit (Promega). RNA samples were assessed for quantity and integrity using a Tapestation High Sensitivity RNA ScreenTape (Agilent), and 100 ng of total RNA were used for sequencing. Libraries were constructed using the NEBNext® Ultra™ II Directional RNA Library Prep Kit for Illumina (E7760L, New England Biolabs) according to manufacturer recommendations. Briefly, after fragmentation to approximately 180 nucleotides fragment size, the samples were converted to double stranded DNA and ligated to Illumina adapters. The second strand cDNA was then degraded by the UNG enzyme providing directional information to the library. Each library was uniquely indexed by PCR. Libraries obtained from each sample were quantified and quality-checked on a Tapestation using the High Sensitivity D1000 ScreenTape system prior to pooling. A HiSeq4000 was used to generate 75 paired-end sequencing reads according to manufacturer recommendations. qRT-PCR of a panel of 12 representative genes was performed to confirm key expression data ( Supplementary Fig. S13).

Tandem mass tag (TMT) labelling proteomics.
Sample preparation and digestion.
Whole cell extracts were prepared from 100 µg per sample. To reduce cysteine residues, tris (2carboxyethyl) phosphine (TCEP, C4706-2G, Sigma-Aldrich) was added to a final concentration to 10 mM and the samples were incubated at 55 °C for 1 h. Subsequently, the reduced proteins were alkylated with 20 mM iodoacetamide (T7408-100ML, Sigma-Aldrich) for 30 min, protected from light, at RT. Protein precipitation was performed with 1:6 volume of pre-chilled (-20 °C) acetone overnight at -20 °C. After overnight incubation, lysates were centrifuged at 8000 × g, for 10 min, at 4 °C. The supernatant was discarded, and the pellets were air-dried for 10 min to remove residual acetone and then dissolved in 100 µL of 100 mM TEAB. Trypsin (Trypsin Gold, Mass Spectrometry Grade, Promega) was added at a 1:40 (trypsin: protein) ratio and digestion was allowed to proceed at 37 °C for 18 h.

TMT Labelling and Sample Fractionation.
After in-solution digestion the samples were labelled with 0.2 mg of TMT isobaric label reagents reconstituted in anhydrous acetonitrile (Thermo Fisher Scientific). Peptide samples containing 50 µg of protein were labelled with the corresponding TMT 10-plex reagent for 1 h at RT. The reactions were quenched using hydroxylamine added to a final concentration of 0.3% (v/v) for 30 min. TMT-labelled peptides were combined and fractionated using Pierce High pH Reversed-Phase Peptide Fractionation Kit (Thermo Fisher Scientific) following the manufacturer's instructions. All fractions (nine fractions in each set) were dried by vacuum centrifugation and stored at -80 °C until mass spectrometric analysis.

LC-MS/MS analysis.
Peptide mixtures from TMT 10-plex labelled samples were chromatographically resolved on an EASYspray PepMap RSLC C18 column (2 µm, 100Å, 75 µm × 50 cm ID) using an Ultimate 3000 RSLCnano system (Thermo Fisher Scientific) over a 180 min gradient at 40 °C. The peptides were separated using linear gradient of 2-45% solvent B over 145 min. The column was washed with 95% B for 15 min and equilibrated with 2% B for the rest of the acquisition. Solvent A was 0.1% FA and 5% DMSO in HPLCgrade water, and solvent B was 0.1% FA and 5% DMSO in 80% acetonitrile and the flow rate was 250 nL/min. All spectra were acquired using an Orbitrap Fusion Lumos Tribrid mass spectrometer (Thermo Fisher Scientific) and Xcalibur 4.2.28.14 (Thermo Fisher Scientific) was used to control data acquisition. The instrument was operated in data dependent acquisition mode with top scan speed set at 3 s. MS1 spectra were acquired in the Orbitrap at a resolution of 120000 and an ion target of 4E5. The MS2 precursors were isolated using the quadrupole (0.7 Th window) and analysed in the Orbitrap at a resolution of 60000, with an AGC target of 1E5 and a max injection time of 118 ms. Precursors were fragmented by HCD at a normalized collision energy (NCE) of 38%.
Data processing and analysis.
Orbitrap.RAW files were analysed by MaxQuant (version 1.6.0.13), using Andromeda for peptide search. MS2 spectra were searched against the Uniprot database specifying Homo sapiens (Human) taxonomy (Proteome ID: UP000005640, Organism ID: 9606, number of entries: 21039). A list of 247 common laboratory contaminants provided by MaxQuant was also added to the database. All searches were performed using "Reporter ion MS2" with "10-plex TMT" as isobaric labels. For the search the enzyme specificity was set to trypsin with maximum of two missed cleavages. The precursor mass tolerance was set to 20 ppm for the first search (used for mass re-calibration) and to 4.5 ppm for the main search.
Carbamidomethylation of cysteines was specified as fixed modification, oxidized methionines and Nterminal protein acetylation were searched as variable modifications. Reporter mass tolerance was set to 0.003 Da. TMT signals were corrected for isotope impurities based on the manufacturer's instructions. The dataset was filtered on posterior error probability to achieve 1% FDR on peptide and protein level.

LC-MS metabolomics.
Sample Preparation. 10 x 10 6 live cells as determined by Trypan Blue exclusion were used per sample. Briefly, samples were prepared using the automated MicroLab STAR® system. Several recovery standards were added prior to the first step in the extraction process for quality control purposes. To remove protein, dissociate small molecules bound to protein or trapped in the precipitated protein matrix, and to recover chemically diverse metabolites, proteins were precipitated with methanol under vigorous shaking for 2 min, followed by centrifugation. The resulting extract was divided into five fractions: two for analysis by two separate reverse phase (RP)/UPLC-MS/MS methods with positive ion mode electrospray ionization (ESI), one for analysis by RP/UPLC-MS/MS with negative ion mode ESI, one for analysis by HILIC/UPLC-MS/MS with negative ion mode ESI, and one sample was reserved for backup. Samples were placed briefly on a TurboVap® (Zymark) to remove the organic solvent. The sample extracts were stored overnight under nitrogen before preparation for analysis.

Ultrahigh Performance Liquid Chromatography-Tandem Mass Spectroscopy (UPLC-MS/MS).
The sample extract was dried and then reconstituted in solvents compatible to each of the four analysis. Each reconstitution solvent contained a series of standards at fixed concentrations to ensure injection and chromatographic consistency. One aliquot was analysed using acidic positive ion conditions (chromatographically optimized for more hydrophilic compounds), and the extract was gradient eluted from a C18 column (Waters UPLC BEH C18 -2.1 x 100 mm, 1.7 µm) using water and methanol, containing 0.05% perfluoropentanoic acid (PFPA) and 0.1% formic acid (FA). Another aliquot was also analysed using acidic positive ion conditions and the extract was gradient eluted from the same afore mentioned C18 column using methanol, acetonitrile, water, 0.05% PFPA and 0.01% FA and was operated at an overall higher organic content. Another aliquot was analysed using basic negative ion optimized conditions using a separate dedicated C18 column. The basic extracts were gradient eluted from the column using methanol and water, with 6.5 mM ammonium bicarbonate at pH 8. The fourth aliquot was analysed via negative ionization following elution from a HILIC column (Waters UPLC BEH Amide 2.1 x 150 mm, 1.7 µm) using a gradient consisting of water and acetonitrile with 10 mM ammonium formate, pH 10.8. The MS analysis alternated between MS and data-dependent MSn scans using dynamic exclusion. The scan range varied slighted between methods but covered 70-1,000 m/z.

Data Extraction and Compound Identification.
Raw data was extracted, peak-identified and quality control-processed using Metabolon's hardware and software. These systems were built on a web-service platform utilizing Microsoft's .NET technologies, which run on high-performance application servers and fiber-channel storage arrays in clusters to provide active failover and load-balancing. Compounds were identified by comparison to library entries of purified standards or recurrent unknown entities. Metabolon maintains a library based on authenticated standards that contains the retention time/index (RI), mass to charge ratio (m/z), and chromatographic data (including MS/MS spectral data) on all molecules present in the library. Furthermore, biochemical identifications are based on three criteria: retention index within a narrow RI window of the proposed identification, accurate mass match to the library +/-10 ppm, and the MS/MS forward and reverse scores between the experimental data and authentic standards. The MS/MS scores are based on a comparison of the ions present in the experimental spectrum to the ions present in the library spectrum. Metabolites were normalised to dsDNA concentration of the samples.

Metabolite consumption and release rate.
High-resolution 1-dimensional 1H NMR spectroscopy was used to profile metabolite consumption and release (3). Frozen culture supernatant samples were thawed at RT and subsequently vortexed. A volume of 550 μL of each sample was transferred into microcentrifuge tubes, 50 μL of DSS buffer (12 mM deuterated 4, 4-dimethyl-4-silapentane-1-sulfonic acid, 0.5 M sodium phosphate, 0.02% sodium azide in D2O water) were added to each sample before tubes were vortexed and centrifuged at 20000 x g for 1 min. A volume of 550 μL of each sample supernatant was then transferred into 5 mm diameter Wilmad® NMR tubes (Sigma-Aldrich, Missouri, US).
High-resolution 1-dimensional 1H NMR spectroscopy was performed using the 14.1 T Bruker AVANCE 400 MHz spectrometer (Bruker Biospin) at 298 K. NMR spectra were acquired for each sample using a conventional ZGPR solvent pre-saturation method with a single RF pulse, a recycle delay of 4 s, spectral width of 6 402.049 Hz, 32 FIDs and 64000 data points. Data were automatically digitalised, and Fourier transformed before being processed in MATLAB® software (Mathworks) using in-house scripts developed in the Department of Surgery and Cancer at Imperial College London (London, UK). Phase correction, baseline correction, ppm referencing (DSS) and normalisation to the reference peak (DSS) was automatically done before spectral peaks were identified with reference to the Human Metabolome Database (HMDB) and subsequently integrated using in-house scripts in in MATLAB® software Oxygen consumption rate (OCR) was determined using a Seahorse XFe96 Analyzer (Agilent) following manufacturer's recommendations. Briefly, Seahorse cartridges (102601-100, Agilent) were hydrated and incubated at 37 °C in a CO2-free incubator for 24 h before analysis. Seahorse 96-well plates (102601-100, Agilent) were coated with 0.01% poly-lysine solution (4832, Sigma Aldrich). On the day of the analysis, the Seahorse XFe96 analyser was calibrated using the hydrated cartridge. A total of 10 5 live RPMI-8226 cells per well were plated in XF assay medium (102365, Agilent) complemented with 10 mM glucose (G8769, Sigma Aldrich), 1 mM sodium pyruvate (1136007, Gibco) and 2 mM glutamine (25030081, Gibco) and incubated for 1 h at 37 °C without CO2. OCR was quantified after injection of 1.5 µM oligomycin A, 60 µM 2,4-Dinitrophenol (DNP) (34334, Sigma Aldrich), or 1 µM rotenone/antimycin A. Results were acquired and analysed by the Wave 2.6.1 software (Seahorse bioscience Inc.). Data were normalized to viable cell numbers at the time of analysis as determined by Trypan Blue exclusion, and day 0 to 8 samples were analysed in parallel on one plate.

Gene expression time series modelling and clustering.
RNA-sequencing-derived transcript counts for 18062 genes were regularised and log2-transformed using DEseq2 (4). Low counts were filtered out by applying an expression cut-off at the trough between the two modes of the distribution of the log2-transformed sum of normalized counts. This resulted in 11798 transcripts selected for further analysis. To study the dynamics of the recovery with respect to the expression levels in day 0, we removed the expression data on day 0 (timepoint ! ) from analysis and renormalised the remaining timepoints by subtracting the mean of the log2-transformed regularised counts for the 11798 transcripts. Lastly, we filtered for genes that underwent substantial changes in their expression profiles. To this end, we computed the variance: with " being the expression level averaged across all timepoints. We applied a cut-off var( " ) > 0.1, which resulted in 2542 genes selected for further analysis.
We modelled the time series with a combination of Gaussian Process (GP) regression and graph-based clustering (5). For GP regression we employed a squared exponential function as a covariance function and learned the GP hyper-parameters by maximizing the marginal likelihood for all 2542 genes(4). We then used GP likelihood to compute a similarity score between every pair of genes, which were then converted into a gene-to-gene similarity graph. For clustering, we pruned the graph using k-nearest neighbour (kNN) algorithm with = 7, where a gene is connected only to its seven nearest neighbour genes based on the GP similarity. When the kNN graph had disconnected components, we intersected it with the minimum spanning tree of the fully connected graph. The resulting graph was clustered with the Markov Stability (MS) algorithm for community detection (5). MS produces different graph partitions depending on the Markov time, a resolution parameter that allows to examine the graph structure across scales. We determined the number of clusters by three simultaneous criteria: 1) the algorithm identifies the same number of clusters for a wide range of Markov times (blue curve in SI Appendix Fig.  S2A); 2) the identified clusters are stable, as quantified by the variation of information across Markov times (green curve and heatmap SI Appendix Fig. S2A); 3) the clusters are biologically meaningful, as quantified by the biological homogeneity index (BHI) (6). Based on these criteria the filtered genes were clustered into 6 groups (SI Appendix Fig. S2B,C). The average gene expression profile for each cluster was computed with GP regression of the genes in each cluster. For enrichment analysis we employed the R package clusterProfiler to identify the Gene Ontology (GO) terms and KEGG pathways that were over-represented by the clusters of genes (7). A p-value cut-off of 0.05 was used to retrieve relevant terms and pathways.

Statistical analysis.
Unless stated otherwise, data show means and standard error of the mean (SEM) from at least three independent experiments. Statistical tests are indicated in the figure legends and were performed using GraphPad Prism 8.4.2 software. Differences among groups were considered statistically significant at P < 0.05 unless stated otherwise.

RNA-sequencing.
Read count abundances were generated from raw data FASTQ files by pseudo-aligning these to the Homo sapiens GENCODE 'comprehensive' reference transcriptome (GRCh38.p12 / release 28) using Kallisto v0.44.0, adjusting for GC content bias and bootstrapping 50x. Bootstrapped transcript isoformlevel counts were then imported to R Programming Language v 3.5.0 (R) and summarised to gene level counts (adjusting for gene-length) using the tximport v1.9.8 (8) package in R.
Protein coding genes (total = 19665) were then isolated from the raw counts based on GENCODE biotype keyword 'protein_coding'. A gene with zero counts across all samples was removed (total = 1603). The raw counts of the remaining genes (total = 18062) (as a tximport object) were then converted to a DESeq2 v1.25.8(4) object for normalisation with betaPrior set to FALSE. For downstream analyses, the negative binomial-distributed normalised counts were transformed to regularised log (rlog) and variance-stabilised (vst) expression levels via the rlog and vst functions of DESeq2 in R, with blind set to FALSE.
Differential expression analysis was conducted on the negative binomial-distributed normalised counts with FDR set at 5%. Following differential expression, log2 fold changes (log2FC) were shrunk via the lfcshrink function of DESeq2. A gene was defined as differentially expressed if it passed Benjamini-Hochberg Q ≤ 0.05 and absolute log2FC ≥ 2 unless stated otherwise.

Metabolomics.
The raw area count metabolite levels were normalised and imputed by Metabolon as follows: 1, each sample was normalised by dsDNA concentration; 2, each metabolite was rescaled to have median = 1; and 3, missing values were imputed with the minimum overall value. These values were then log2 transformed and converted to the Z-scale for the purposes of all downstream analyses. Differential metabolite analysis was performed using RegParallel (https://github.com/kevinblighe/RegParallel. doi:10.18129/B9.bioc.RegParallel), which fits a linear model independently to each metabolite level, with final nominal p-values adjusted for FDR (Benjamini-Hochberg).

Clustering and heatmaps.
Supervised (filtered) and unsupervised (unfiltered) clustering were performed using the Heatmap function of the ComplexHeatmap v2.1.0 package (12). Data was transformed to the Z scale (scaled by row/gene) and then clustered via 1 minus Pearson correlation distance (or Euclidean distance) and Ward's linkage. A samples box-and-whisker plot of Z-scores was added to the heatmap bottom. In certain cases, k-means clustering was performed on heatmap data in order to identify k clusters, with the variable-to-cluster assignment then being used to split the heatmap and dendrogram into k separate entities.

Fold-change comparisons: RNA-seq vs. Proteomics.
For fold-change comparisons, the log2FC for all RNA-seq and proteomics genes/proteins that matched on Entrez ID were included, irrespective of statistical significance. biomaRt v2.41.7 (13,14) was used to convert gene symbols to Entrez IDs. Log2FCs were then extracted for day 0 comparisons for these genes/proteins and then plotted via geom_point from ggplot2. A linear regression fit with 95% confidence intervals was added via geom_smooth and stat_smooth. Pearson correlation values and the p-values from the correlation test were added as captions. Plots were arranged via grid.arrange from the gridExtra package.

GO term and KEGG pathway enrichment.
For GO term enrichment, topGO v2.37.0 was used. Enrichment was made against "BP" (biological processes), "CC" (Cellular Component), and "MF" (Molecular Function) ontology with nodeSize = 5. The "classic" algorithm with the Kolmogorov-Smirnov (K-S) test, which takes the rank of each target gene into consideration, was used to determine statistically significantly enriched terms (P < 0.001) -only top 10 terms were chosen.
Pathway enrichment for KEGG was performed using KEGGprofile v1.27.0. Initial enriched pathways were identified via the find_enriched_pathway function, with returned_pvalue = 0.01, returned_adjpvalue = 0.05, and and returned_genenumber = 5, or returned_genenumber = 3 for the GCN2-dependency TCGA signature enrichments due to lesser numbers of genes in each signature. This process was repeated 10x to increase reliability of results, randomly selecting 90% total genes during each bootstrap, and selecting the lowest obtained adjusted p-value where enriched pathways were the same across multiple bootstraps. The final list of pathways was then filtered by adjusted P < 0.001. Prior to topGO and KEGGprofile, gene symbols were converted to Entrez IDs via biomaRt, with non-matching symbols being removed. For the GCN2-dependency TCGA signatures, no cross validation was performed, again due to lesser numbers of genes in each signature.
For each comparison, genes passing Benjamini-Hochberg Q ≤ 0.05 were used. Final plots were drawn via geom_bar from ggplot2, as negative log10 of the enrichment p-value.

Gene-set variation analysis.
Unbiased/Unsupervised gene-set variation analysis (GSVA) was performed using GSVA v1.33.1 (15) and using the C2 curated gene sets from Molecular Signatures Database (C2 MsigDB), with parameters min.sz = 5, max.sz = 999999, and abs.ranking = FALSE. Input in each case were rlog expression levels with gene symbols converted to Entrez IDs via biomaRt. Gene symbols that did not match were removed. For comparisons to day 0, limma was used to fit a linear model to the enriched dataset for the purpose of deriving test statistics. Test statistics were modified via the empirical Bayes methods via ebayes. Separately, GSVA enrichments were also made against subsets of C2 MSigDB signatures. These subsets were identified via key terms as indicated in the manuscript. The resulting enrichments were then plot via corrplot v0.84, with enrichment level scaled between -1 and +1.

Mitocarta 2.0 analyses.
Genes/Proteins passing Benjamini-Hochberg Q ≤ 0.05 were included, and these were separated for those that were down-or up-regulated (negative or positive fold-change, respectively), and presence / absence in MitoCarta2.0. A Pearson χ2 test was performed at each day to gauge the statistical significance of MitoCarta2.0 vs non-MitoCarta2.0 statistically significant genes. ggplot2 was then used to plot count frequencies of these as a stacked bar plot, with each bar colour-coded based on MitoCarta2.0 / non-MitoCarta2.0 genes/proteins. Gene-set variation enrichments using GSVA were also performed separately against just MitoCarta2.0 known mitochondrial genes, i.e., MitoCarta2.0 was the signature being enriched.

DepMap GCN2 dependency and TCGA.
Raw RNA-seq counts were downloaded from CCLE (https://portals.broadinstitute.org/ccle) as CCLE_RNAseq_genes_counts_20180929.gct.gz, with accompanying metadata Cell_lines_annotations_20181226.txt. Counts over cell-lines of interest were extracted and then different analyses were conducted using DESeq2: § Normalisation of all samples followed by differential expression between cell-lines with highest DepMap GCN2 dependency (loCERES) and lowest GCN2 dependency (hiCERES), with inclusion of "tissue" as covariate. § Separate normalisations and differential expression (as above) for each main tissues of interest (CNS, liver, and skin) separate.
Downstream heatmaps were generated as previously described in these methods. These were then normalised via DESeq2 -variance-stabilised expression levels were also produced for downstream analyses.
For the construction of tissue-specific predictive models (one each for CNS, liver, and skin) using the CCLE data, we first identified statistically significant genes between loCERES vs. hiCERES at Benjamini-Hochberg Q ≤ 0.05 and |log2FC| ≥ 1 for each tissue. Using glmnet (16), we then fit an elastic-net penalised regression model and cross-validated this x10. The best predictors from this model were identified as those whose coefficients were not shrunk to zero based on the 1 standard error rule (lambda.1se). Enriched GO terms and KEGG pathways were identified from these genes as previously described in these methods.
For each tissue specific CCLE model, we then made predictions on each respective TCGA dataset, i.e.: CNS à GBM + LGG; skin à SKCM; liver à LIHC. Prior to model predictions, the variance-stabilised expression levels of each TCGA dataset was standardised by subtracting, per gene, the mean, and then dividing by the standard deviation. Only TCGA samples with predictions at > 80% probability were retained. Where possible, and numbers per predicted class permitted, a differential expression analysis was conducted between the loCERES and hiCERES predicted TCGA samples. A logistic regression analysis was also performed against the TCGA biotab clinical data to discover statistically significant associations.
For gene-set enrichment analysis, we used the fgsea package in R. Enrichment was made against selected KEGG pathway with 100000 permutations. Input genes were ranked based on negative log of the adjusted p-value for the predicted low versus high GCN2-dependency TCGA samples.