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) DLJr_PS_Scripts_CoreAnalyses_071424.zip--Scripts for performing core analyses detailed in our manuscript

(2) DLJr_PS_Scripts_FigureGeneration_061024.zip--Scripts for generating figures in manuscript

(3) Data_For_Input_Zenodo_v3.zip--Data that serves as input for the scripts

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

(5) MachineLearning_Schematic_2_Proteomics.png--Depiction of machine learning algorithm's core principles for proteomics

(6) MachineLearning_Schematic_3_Metabolomics.png--Depiction of machine learning algorithm's core principles for metabolomics

 

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.

 

Important Notes for DLJr_PS_Scripts_FigureGeneration_061024.zip

DLJr_Fig1_Generation_051124.R:

Fig 1b in the script corresponds to Fig S3a in the manuscript

Fig 1c in the script corresponds to Fig S3b in the manuscript

Fig 1d in the script corresponds to Fig S3c in the manuscript

DLJr_Fig2_Generation_051124.R:

Fig 2d in script does not correspond to any figure panels in the maunscript

Fig 2e in the script corresponds to Fig 2d in the manuscript

DLJr_Fig3_Generation_051124.R:

Fig 3g in the script corresponds to Fig S4a in the manuscript

Fig 3h in the script corresponds to Fig S4b in the manuscript

Fig 3i in the script corresponds to Fig S4c in the manuscript

Fig 3j in the script corresponds to Fig S4d in the manuscript

DLJr_Fig4_Generation_051124.R:

Figure panels 4a-g in script do not correspond to any figure panels in the maunscript

Fig 4h in the script corresponds to Fig 4a in the manuscript

Fig 4i in the script corresponds to Fig 4b in the manuscript

Fig 4j in the script corresponds to Fig 4c in the manuscript

Fig 4k in the script corresponds to Fig 4d in the manuscript

Fig 4l in the script corresponds to Fig 4e in the manuscript

Fig 4m in the script corresponds to Fig 4f in the manuscript

Fig 4n in the script corresponds to Fig 4g in the manuscript

DLJr_Fig5_Generation_051124.R:

Fig 5b in script does not correspond to any figure panels in the maunscript

Fig 5d in the script corresponds to Fig 5c in the manuscript

Fig 5e in the script corresponds to Fig 5d in the manuscript

Fig 5f in the script corresponds to Fig 5e in the manuscript

DLJr_Fig6_Generation_051124.R:

Fig 6g in the script corresponds to Fig S5a in the manuscript

Fig 6h in the script corresponds to Fig S5b in the manuscript

Fig 6i in the script corresponds to Fig S5c in the manuscript

Fig 6j in the script corresponds to Fig S5d in the manuscript

Fig 6k in the script corresponds to Fig S5e in the manuscript

DLJr_Fig7_Generation_051124.R:

Fig 7h in the script corresponds to Fig S6b in the manuscript

DLJr_Fig8_Generation_051124.R:

Fig 8 in the script corresponds to Fig S7 in the manuscript

 
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 with nested cross-validation (CV) [MachineLearning_Schematic_1; MachineLearning_Schematic_2_Proteomics; MachineLearning_Schematic_3_Metabolomics] was trained on our omics datasets, and the pipeline was subsequently assessed if it was better than random chance at classification of FLC. Total samples (ntotal_proteomics = 36; ntotal_metabolomics = 36) were split into a training set (ntraining_proteomics = 23 and ntraining_metabolomics = 26. These samples were used for machine learning as well as other analyses detailed in this study) and hold-out set (nholdout_proteomics = 13 and nholdout_metabolomics = 10. These samples were used exclusively for machine learning). Data processing was performed independently on the training and hold-out sets to minimize data leakage [14], with the exception of imputation for metabolomics which was performed on all the samples combined. Training data went through two (proteomics) to three (metabolomics) iterations of inner-loop train /test splitting (ninner-loop train_proteomics = 13 samples and ninner-loop test_proteomics = 10 samples; ninner-loop train_metabolomics = 18 samples and ninner-loop test_metabolomics = 8 samples). With each iteration, the inner training (IT) data was subsetted for the most relevant features using a PCA-based method [15,16] that involved the following steps: (1) performing PCA with centering (proteomics, metabolomics) and scaling (metabolomics), (2) selecting principal components that explained greater than 10% of the variance in the data, (3) performing varimax rotation on raw loadings, and (4) isolating features that were greater than or equal to an absolute value of 0.5 in any of the selected principal components. The subsetted IT data was subsequently used for stratified 2-fold,3-repeat CV where each fold (resample) contained at least six samples (nNML_proteomics_resample = 2 and nFLC_proteomics_resample > 4; nNML_metabolomics_resample = 2 and nFLC_metabolomics_resample = 5) within each partition for training and testing. Hyperparameter tuning (described in more detail below) was incorporated into the 2-fold, 3-repeat CV. The best model from the inner training was assessed for performance against the inner test data (ninner-loop test_proteomics = 10 samples; ninner-loop train_metabolomics = 8 samples). Overall evaluation of training performance was attained by averaging of the two (proteomics) to three (metabolomics) training assessments from the two (proteomics) to three (metabolomics) iterations of inner training. The union of the feature subsets from the two (proteomics) to three (metabolomics) iterations was used to subset the whole training dataset (ntraining_proteomics = 23 and ntraining_metabolomics = 26), and 2-fold,3-repeat CV with hyperparameter tuning was performed. A final evaluation was executed by assessing the performance of the trained model from the whole training dataset against the hold-out set (nholdout_proteomics = 13 and nholdout_metabolomics = 10). 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 (the performance of the classifier compared to classification at random according to class frequency), and learning rate (lower values indicate better learning).

Hyperparameter tuning was performed with the ‘bayesOpt’ function of the ‘ParBayesianOptimization’ package in R [17]. 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 ten optimization steps with ten samplings of the search space (optimums) 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 on whole optimization process. Additionally, we utilized the upsampling method, SMOTE (Synthetic minority over-sampling technique) [18],  during resampling to account for class imbalances (FLC samples >> NML samples).

Eight datasets were used for XgBoostML nested CV. 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 postprocessed) were the remaining four inputs. Each input went through ten consecutive rounds of nested CV, and the following comparisons were made with each performance metric: noise-to-adjusted noise, noise-to-original preprocessed, adjusted noise-to-original postprocessed, and original preprocessed-to-original postprocessed.