Proteo-metabolomics and patient tumor slice functional studies point to the centrality of amino acids for re-wired mitochondria in fibrolamellar carcinoma
Creators: name: Donald Long Jr., Don Long
Publisher: name: Zenodo
Year: 2024
Subjects: Machine Learning
Concordance Analysis
Correlation Analysis
Differential Protein Expression
Proteomics
Metabolomics
Data Imputation
Data processing
Figure Generation

description: Short description of the uploaded folders and files

(1) SourceData_TablesS1-S16_DataS1-S3_061024.zip--All source tables and source data that accompany manuscript

(2) CRM_ComplimentaryFigs&Legends_061024.pdf--Complimentary figures that accompany responses to peer reviewers of manuscript

(3) DLJr_PS_Scripts_CoreAnalyses_061024.zip--Scripts for performing core analyses detailed in our manuscript

(4) DLJr_PS_Scripts_FigureGeneration_061024.zip--Scripts for generating all of the main figures in our manuscript

(5) Data_For_Input.zip--Data that serves as input for the scripts

(6) MachineLearning_Schematic_1.png--General overview of machine learning pipeline

(7) MachineLearning_Schematic_2.png--Depiction of machine learning algorithm's core principles

 

After download of scripts and input data, be sure (1) to create a folder entitled "the_results" in your present working directory, and (2) that all files in 'Data_For_Input.zip' are located in a folder entitled 'Data_for_Input' once unzipped.
description: Fibrolamellar carcinoma (FLC) is a rare, early-onset liver cancer.  The five-year survival rate is low, and there is a critical need for new therapeutics.  The primary driver in FLC is the fusion oncoprotein, DNAJ-PKAc, which remains challenging to target therapeutically.  It is critical, therefore, to expand understanding of the FLC molecular landscape to identify druggable pathways/targets. We have performed the most comprehensive characterization of FLC in both proteomics (n = 23 patient samples; n = 8485 proteins) and metabolomics (n = 26 patient samples; n = 135 metabolites). Additionally, we have conducted respirometry analyses as well as RNAi- and small molecule inhibitor-mediated loss of function assays in FLC tumor tissue slices from patients.  We propose a model of cellular energetics in FLC that centers on amino acids.  Our model points to proline anabolism that is very likely mediated by ornithine aminotransferase hyperactivity and ornithine transcarbamylase hypoactivity with serine and glutamine catabolism fueling the process. We also highlight a dependency of FLC tumors on VDAC, which functions as a gatekeeper for mitochondrial metabolites including pyruvate. The metabolic rewiring in FLC proposed by our model, with a particular emphasis on mitochondria, can be exploited for therapeutic vulnerabilities.
description: Data Processing and Differential Abundance Analysis (Proteomics)

Samples were run through LC MS/MS in two consecutive batches (16 samples per batch). Protein quantification by TMT16plex reporter ions was calculated using TMT PSM ratios to the entire sample set (all 16 TMT-channels) and normalized to the sample median. The median PSM TMT reporter ratio from peptides unique to a gene symbol was used for quantification. Protein false discovery rates were calculated using the picked-FDR method using gene symbols as protein groups and limited to 1% FDR [1]. Quantified proteins (n = 10065) displayed 1) values for all samples across both batches, 2) values missing for all samples in one of the batches, or 3) values missing for all samples across both batches. For the full and omics-matched datasets, only proteins containing values for all samples across both batches were retained (n = 8485). Covariate adjustment using R was performed to remove batch effect. Python and R were used to generate hierarchal clustering with heatmap and PCA-based clustering, respectively. Data profiles from metastatic (MET) and primary (PRI) tissue samples were aggregated and deduplicated (random removal of either MET or PRI that came from the same patient) to create the FLC cohort. DeqMS [2] was used for differential protein expression analysis.

Data Processing and Differential Abundance Analysis (Metabolomics)

Observations that were either blank or displayed a zero value were considered missing [3,4]. Subsequently, metabolites (n=180) were filtered using the modified 80% rule where any metabolite that contained less than or equal to 20% missing values in at least one of the cohorts ( tumor or adjacent non-malignant tissue) were retained (n =135). At this point, approximately 9% of observations were missing values. The following five different imputation methods were assessed for optimally imputing missing values in our dataset: imputeBasic, multivariate imputation by chained equations (MICE) [5], missMDA [6], missForest [7], and softImpute [8]. A test dataset (n = 76) was derived from the metabolomics dataset where the modified 80% rule had been applied. The test dataset was void of any missing values. For each imputation method (except for MICE, a very computationally intensive platform) one-hundred simulations were run generating three different missing data patterns (missing completely at random or MCAR, missing at random or MAR, and missing not at random or MNAR) [9] at 3%, 10%, 30%, or 50% of the data missing. Thirty simulations were run for MICE. In total, 1200 simulations each (100 simulations x 3 missing data patterns x 4 different degrees of data missing) were run for imputeBasic, missMDA, missForest, and sofImpute, while 360 simulations (30 simulations x 3 missing data patterns x 4 different degrees of data missing) were run for MICE. Five quality metrics generated from the simulations were used to compare and contrast the imputation methods—raw bias (difference between predicted value of point estimate and it’s true value from the test dataset), percent bias (raw bias expressed as a percent), coverage rate (proportion of confidence intervals generated from simulations that contain the true value), average width (average width of confidence interval), and root mean squared error (represents compromise between bias and variance or accuracy and precision). Given the distribution of missing values in our metabolomics dataset and our unawareness of any bias that would influence the distribution, we presumed the missing values were MCAR. At a missing value rate of 10%, missForest displayed the best combination of quality metrics, and it was used for imputing missing values (~ 9%) in our full and omics-matched metabolomic datasets. Samples were run through LC/MS in four consecutive batches. In the resulting peak intensity file, MET and PRI tissue samples were aggregated and deduplicated (random removal of either MET or PRI that came from the same patient) to create the FLC cohort. Subsequently, the peak intensity files from the full and omics-matched data was processed in the following manner after imputation: (1) batch median normalization (each observation of a metabolite divided by its batch median), (2) log2 transformation, (3) co-variate adjustment for batch, (4) reverse-log2 transformation, (5) metabolite median normalization (each metabolite within a sample divided by the median of all metabolites within that sample), (6) log2 transformation and (7) co-variate adjustment for injection order. Differential abundance analysis was performed using MetaboAnalyst [10,11].

Concordance Analysis (Proteomics vs Metabolomics)

Omics-matched tumor samples (patient tumor samples that were present in both the proteomics and metabolomics datasets; n = 13 out of the 17 omics-matched samples) were used as input for pipeline adapted from Benedetti et al 2023 [12] to assess level of agreement between proteomic and metabolomic abundance in FLC. Importantly, concordance assumes no underlying distribution. We leveraged the ‘concordance’ function from the ‘survival’ package in R to perform the computations (n = 1,009,715 measured associations).

Concordance Analysis (Proteomics vs Mitochondrial Abundance)

Mitochondrial abundance was measured via MitoTracker Deep Red (MTDR), and proteins specific to the mitochondria (identified by MitoCarta3.0 [13]) were retained for downstream analysis (n = 1003). Four of the six NML (non-malignant liver tissue adjacent to tumor) samples measured via MTDR were patient-matched with samples in the mitochondria-specific, omics (MSO) dataset. The two unmatched samples measured via MTDR were randomly paired with two samples from the remaining MSO sample set (n = 3). For FLC, three of the six samples measured via MTDR were patient-matched with samples in the MSO dataset. The three unmatched samples measured via MTDR were randomly paired with three samples from the remaining MSO sample set (n = 13).

The library of mitochondrial pathways contained in MitoCarta3.0 (n= 150) were filtered by the following criteria: (1) The sub-pathways of generalized over-arching pathways were retained to avoid redundancy of protein members. For example, “Gluconeogenesis” was retained, but its overarching pathway, “Carbohydrate Metabolism”, was removed. (2) Retained pathways had at least one protein member identified in the MSO dataset that was significantly concordant or discordant with mitochondrial abundance. (3) Retained pathways had at least ten protein members. Subsequently, mitochondrial pathways were ranked by (1) highest proportion of pathway members that were significantly associated with mitochondrial abundance (pathway members identified in MSO dataset that were significantly associated with mitochondrial abundance divided by all pathway members identified in MSO dataset), (2) highest proportion of pathway members identified in MSO dataset, and (3) absolute value of scaled concordance (2 * concordance – 1). Finally, a threshold was set for attaining the top twenty mitochondrial pathways where at least 10% of a pathway’s protein members that were identified in the MSO dataset were significantly associated with mitochondrial content.

To address the potential issue of reproducibility associated with random pairing, we performed Monte Carlo simulations with the aforementioned pipeline for identifying all possible sets of top twenty mitochondrial pathways from all possible sample combinations of NML (4 matched samples + [3choose2] = 3 possible combinations) or FLC (3 + [13choose3] = 286 possible combinations). Subsequently, we attained the top twenty mitochondrial pathways with the greatest hits. A hit was considered a gene that was significantly associated with mitochondrial abundance in at least one of the simulations.

Machine Learning Pipeline

An extreme gradient boost machine learning (XgBoostML) pipeline (MachineLearning_Schematic_1.png) with Nested 4-fold (proteomics) and 3-fold (metabolomics) cross-validation (CV) was established to assess the predictive power of differentially abundant proteins and metabolites. FLC and NML samples from proteomics and metabolomics data were stratified across four and three outer folds, respectively. All folds contained at least two NML samples, with the exception of one-fold that contained only one NML sample. With each iteration, a different fold served as a hold-out for testing while the remaining folds served to generate a training model (MachineLearning_Schematic_2.png). The outer training folds were further stratified across two inner folds, one training and one testing (MachineLearning_Schematic_2.png). As with the outer folds, successive iterations (2-fold CV) gave both folds a chance to serve as a training and testing fold (MachineLearning_Schematic_2.png).

The Inner folds were used in finding the optimal tuning parameters with the ‘bayesOpt’ function of the ‘ParBayesianOptimization’ package in R [14]. The following tuning parameters were considered: max_depth (maximal depth of decision tree), eta (learning rate, lower rates allow for better granularity of the model), gamma (modulator for tree complexity), colsample_bytree (the proportion of features used for each tree), min_child_weight (additional modulator for tree complexity), nrounds (number of boosting rounds), and subsample (proportion of samples used for each tree). A grid design containing fifty different combinations of tuning parameters was used to initialize the search (generated by using the ‘maximinLHS’ function from the ‘LHS’ package in R). After initialization, the optimization algorithm was programmed to perform twenty optimization steps with three samplings (local optimums) of the search space at each step. The acquisition function type used to assess performance was excepted improvement (ei) accompanied by the tunable parameter, epilson, which balances exploitation of the search space with exploration. Other important parameters considered were the following: convergence threshold (maximum acceptable difference in ei between successive optimization steps), utility (score that represents the potential to find a better parameter set than the previous), acquisition threshold (percentage of global optimal utility required for a local optimum to be included as a candidate parameter set in the successive optimization step), and time limit per optimization step.

The optimal tuning parameters found from the Bayesian optimization strategy were used, along with the inner folds, to generate a 2-fold,5-repeat cross validation model with differentially abundant proteins or metabolites as input (MachineLearning_Schematic_2.png). Additionally, we utilized the upsampling method, SMOTE (Synthetic minority over-sampling technique) [15], to account for class imbalances in our inner folds (FLC samples >> NML samples). Subsequently, the training model was tested for its performance using the outer fold test hold-out. Performance metrics of interest were balanced accuracy (the accuracy of the training model taking class imbalance of the testing set into account), precision-recall area under curve (similar to ROC AUC but takes class imbalance of the testing set into account), Kappa (level of agreement between training and testing sets), and learning rate.

Eight datasets were used for nested K-fold input. Background datasets (noise) were created by permutating the original datasets (proteomics and metabolomics) before data processing (normalization, transformation, or co-variate adjustment). A second set of background datasets (adjusted noise) were generated by processing the former datasets (noise). The original datasets (both preprocessed and processed) were the remaining four inputs. Each input went through five consecutive rounds of nested K-fold CV, and the following comparisons were made with each performance metric: noise-to-adjusted noise, noise-to-original preprocessed, noise-to-original processed, adjusted noise-to-original preprocessed, adjusted noise-to-original processed, original preprocessed-to-original processed.