A Robust Panel Based on Mitochondrial Localized Proteins for Prognostic Prediction of Lung Adenocarcinoma

Abstract

Due to high energy and material metabolism requirements, mitochondria are frequently active in tumor cells. Our study found that the high energy metabolism status is positively correlated with the poor prognosis of patients with lung adenocarcinoma. We constructed a scoring system (mitoRiskscore) based on the gene expression of specific mitochondrial localized proteins through univariate and LASSO cox regression. It has been shown that high mitoRiskscore was correlated with a shorter survival time after surgery in patients with lung adenocarcinoma. Compared with the typical TNM grading system, the mitoRiskscore gene panel had higher prediction accuracy. A vast number of external verification results ensured its universality. Additionally, the mitoRiskscore could evaluate the metabolic pattern and chemotherapy sensitivity of the tumor samples. Lung adenocarcinoma with higher mitoRiskscore was more active in glycolysis, and oxidative phosphorylation expression of proliferation-related pathway genes was also significantly upregulated. In contrast, patients with low mitoRiskscore had similar metabolic patterns to normal tissues. In order to improve the accuracy of prediction ability and promote clinical usage, we developed a nomogram that combined mitoRiskscore and clinical prognostic factors to predict the 3-year, 5-year, and 10-year survival rates of patients. We also performed in vitro experiments to verify the function of the key genes in the mitoRiskscore panel. In conclusion, the mitoRiskscore scoring system may assist clinicians to judge the postoperative survival rate and chemotherapy of patients with lung adenocarcinoma.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In recent decades, lung cancer is a common occurrence and the leading cause of cancer-related deaths. According to epidemiological data in 2018, there were 2.09 million lung cancer cases and 1.76 million deaths worldwide [1]. Lung adenocarcinoma (LUAD) accounts for approximately 85% of lung cancers and is the most common histological type. At present, the known factors related to the occurrence of LUAD include gender, age, smoking, environment, and heredity [2]. With the rise of medical technology, lung cancer treatment is gradually increasing, from simple surgical resection to surgery combined with adjuvant chemotherapy, immunotherapy, targeted drug therapy, and biological therapy [3]. Although we have made innovations in lung cancer treatment, LUAD’s annual mortality rate has not decreased as expected. It may be due to our lack of assessment of LUAD patients’ survival status after surgery [4]. If we can accurately predict patients’ survival rate after surgery and each patient’s treatment preference, clinicians can provide patients with accurate postoperative treatment. For example, patients with low predicted survival rates should be reviewed more frequently and treated with high-dose drugs; for patients with high predicted survival rates, more conservative treatments should be used to improve patients’ quality of life.

For the prognostic evaluation of LUAD, although there are some methods (such as age, stage, and grade) to predict patients’ survival rate, these prediction methods are not accurate enough. It is generally believed that patients with early lung cancer have a good prognosis, but the five-year survival rate after surgical intervention fluctuates by 34% (10% to 44%) [5]. Besides, age and stage cannot explain the evaluation basis of patients’ resistance to immunotherapy [6]. Therefore, we need to develop more methods to accurately predict the postoperative survival and treatment preferences of LUAD patients.

Mitochondria are the most important places for energy and matter conversion in organisms. When a cell undergoes a malignant transformation, the mitochondrial metabolism mode, transport mode, kinetics, and response to oxidative stress will change significantly [6]. Among them, the intensity of glucose metabolism in cancer cells is significantly increased to produce more intermediate metabolites [7]. Due to the increased energy demand of cancer cells, the mitochondrial oxidative phosphorylation of (OXPHOS) is also significantly improved [8]. In addition, the decomposition of fatty acids and glutamine in cancer cells was also enhanced to provide intermediates for mitochondrial productivity [9]. The high metabolism of mitochondria in cancer cells produces excessive amounts of reactive oxygen species (ROS), cause the death of normal cells, thereby promoting tumor progression [10]. Abnormal mitochondrial function plays an essential role in the occurrence and development of cancer. From another perspective, we can judge the malignant degree and progress trend of the tumor by evaluating the activity of the mitochondria of the tumor.

We analyzed four LUAD mRNA microarrays and reached the same conclusion. Compared with normal tissues, LUAD has a significant increase in mitochondria’s metabolic activity, especially the expression of genes encoding mitochondrial localization protein (mito-protein genes). Besides, we also found that a vast number of mito-protein genes are significantly related to patients’ prognosis. Therefore, we established a mitochondrial risk score (mitoRiskscore) to predict the prognosis and treatment preferences of LUAD patients. The robustness of the risk score has been verified in multiple external data sets. It brings the possibility of precise treatment of LUAD patients after surgery.

2. Materials and Methods

2.1. Data Source and Study Population

The LUAD RNA-seq dataset used in this study came from the Cancer Genome Atlas project (TCGA) [11]. A total of nine microarray datasets of LUAD were downloaded from Gene Expression Omnibus (GEO). Four datasets containing LUAD samples and control samples were used to analyze the abnormal activity of mitochondria in lung cancer (GSE7670, GSE18842, GSE19188, and GSE31210). Six datasets containing prognostic information of patients were used to verify the prediction ability of mitoRiskscore (GSE3141, GSE8894, GSE31210, GSE50081, GSE68465, and GSE72094). The details about the datasets used in this study are presented in Table 1. The self-composed running scripts, together with the processed results of current study, were merged into a repository that is available at https://github.com/XR-Zhang-group/mitoRiskscore.

Table 1

Basic information of datasets used in this study.

Datasets	Platform	Country	No. of patients	No. of controls	Cancer type	Prognostic information
GSE7670	GPL96	China	31	27	LUAD	—
GSE18842	GPL570	Spain	46	45	NSCLC	—
GSE19188	GPL570	Netherland	45	65	LUAD	—
GSE31210	GPL570	Japan	226	20	LUAD	Yes
GSE3141	GPL570	USA	58	—	LUAD	Yes
GSE8894	GPL570	South Korea	63	—	LUAD	Yes
GSE50081	GPL570	Canada	127	—	LUAD	Yes
GSE68465	GPL96	USA	442	—	LUAD	Yes
GSE72094	GPL15048	USA	398	—	LUAD	Yes
TCGA-LUAD	IlluminaHiSeq	USA	513	59	LUAD	Yes

LUAD: lung adenocarcinoma; NSCLC: non-small-cell lung cancer; TCGA: The Cancer Genome Atlas.

2.2. Pathway Enrichment Analysis

First, the differentially expressed genes between the LUAD samples and the control samples were obtained by using the R package “limma.” Then, the R package “clusterProfiler” [12] was used to perform gene set enrichment analysis (GSEA). The five pathways most related to mitochondrial activity were selected for evaluation. The pathway gene sets (MITOCHONDRIAL PART, HALLMARK FATTY ACID METABOLISM, HALLMARK GLYCOLYSIS, HALLMARK HYPOXIA, HALLMARK OXIDATIVE PHOSPHORYLATION, and HALLMARK REACTIVE OXYGEN SPECIES PATHWAY) were downloaded from the molecular signature database (MSigDB) [13]. Mito-protein genes were obtained from the gene set “MITOCHONDRIAL PART” in the MSigDB. The systematic name of this pathway is M18830.

2.3. Preliminary Screening of Prognostic-Related Genes

Univariate Cox proportional hazard regression analysis was performed on each mito-protein gene to screen for genes that were significantly related to the overall survival rate of TCGA-LUAD [14]. Then, the bootstrapping method was used to test the robustness of these genes [15]. 70% of the samples were randomly selected from the cohort, and the survival impact of these genes was evaluated. This process was repeated 1000 times. A gene that had a valid $P$ value at all times was considered to have passed the bootstrapping robustness test.

2.4. Mito-Protein Gene Panel Generation

Because of the large number of genes that passed the bootstrapping robustness test, it is not convenient for follow-up clinical application. The least absolute shrinkage and selection operator (LASSO) Cox regression analysis was performed on the genes that passed the test. While reducing the number, the genes most related to prognosis were also screened out. Through LASSO analysis, the gene represented by the minimum penalty parameter λ will be selected to establish a prognostic risk score formula (mitoRiskscore). We used the “pheatmap” package to visualize the expression of each gene in mitoRiskscore. The Kaplan-Meier method was used to calculate the survival rate, and the log-rank test was performed to evaluate its statistical significance. Univariate and multivariate cox regression analyses were used to test whether the prediction model based on mitoRiskscore was an independent prognostic factor. The time-dependent receiver operating characteristic curve (TDROC) from the “survivalROC” package [16] was applied to test the predictive ability of mitoRiskscore at 1, 3, and 5 years.

2.5. Gene Set Variation Analysis (GSVA)

Gene set variation analysis is one of the most commonly used methods for analyzing biological processes [17]. The gene set file “http://c2.cp.kegg.v7.3.symbols.gm” downloaded from MSigDB was used for GSVA, and the “GSVA” package in the R software was used. The significance threshold was set at $\text{FDR}<0.05$.

2.6. Prediction of Chemotherapy Response

The R package “pRRophetic” [18] and the TCGA-LUAD data were used to predict the drug sensitivity of each group of patients. Among them, ridge regression was used to estimate the maximum half-inhibitory concentration (IC50) of the samples, and ten-fold cross-validation was used to evaluate the accuracy.

2.7. Cell Culture and Treatment

Human lung adenocarcinoma cells, A549, were cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum, 50 U/mL penicillin, and 50 mg/mL streptomycin. Cells were incubated at 37°C under 5% CO₂. To knockdown VDAC1, A549 cells were transfected with VDAC1 siRNAs (Tsingke Biotechnology, China) using Lipofectamine 2000 regent (Thermo Fisher Scientific, USA) according to the manufacturer’s recommendation. The sequences of the two siRNA were as follows:

siRNA #1-F (5${}^{\prime }$-3${}^{\prime }$): GGAUACACUCAGACUCUAATT

siRNA #1-R (5${}^{\prime }$-3${}^{\prime }$): UUAGAGUCUGAGUGUAUCCTT

siRNA #2-F (5${}^{\prime }$-3${}^{\prime }$): GGAUGGCAAGAACGUCAAUTT

siRNA #2-R (5${}^{\prime }$-3${}^{\prime }$): AUUGACGUUCUUGCCAUCCTT

2.8. RNA Extraction and qRT-PCR

Total RNA was extracted using TRIzol (Invitrogen, USA). qRT-PCR was carried out with PrimeScript RT reagent Kit and TB Green Premix Ex Taq II kit (TaKaRa, China) according to the manufacturer’s protocol. The primers used in the study were as follows:

VDAC1-F1: GCAAAATCCCGAGTGACCCAGA

VDAC1-R1: TCCAGGCAAGATTGACAGCGGT

SLC25A42-F: GCAGCTACTATGGCTTCCGT

SLC25A42-R: TTTCCTTCGGGGTTACGGC

ABAT-F: CAAGGAAAGAGGGCAGAGGG

ABAT-R: GGGTATTTCAGCCGTGGGAA

IVD-F: GTGAGTACATCGGAGCCCTG

IVD-R: TGCTAAAGCCAGGCATACCC

AMT-F: CGAGGCTGGAGGCATCTTAG

AMT-R: ACAGCACTGGTCATGAAGGG

2.9. Extracellular Acidification Rate (ECAR) Detection

Seahorse XF glycolysis pressure test kit was used to test the glycolysis ability of cells. On the night before the experiment, A549 cells were transfected into si-VDAC1 or control si-RNA for 48 hours, and the number of cells per well was paved into XF24 test plate for overnight culture. At the same time, the test board will be soaked in the test solution overnight. On the next day, the cells were removed and washed with glycolysis pressure test solution and then cultured in a CO₂-free incubator for one hour. During the period of cell culture, we prepared the concentration of glucose, oligomycin, and 2-DG. Take out the hydrated test board and add each medicine to the hole as required. Set the experimental template on the XF controller, and the default Mix-Wait-Measure time is 3 min-2 min-3 min. Usually, the basic value is measured three times before adding the drug and then three times after each addition. Put it into the test board and wait for the machine to be calibrated, then replace the utility with the cell culture plate for the ECAR test. Each experiment was repeated three times.

2.10. Western Blotting

72 hours after transfection of si-VDAC1, A549 cells were collected and incubated in 200 μ L 1$\ast$SDS-PAGE Sample Loading Buffer (Cat No. P0015L), 98°C for 10 minutes until the cells were completely lysis. Proteins were resolved on 10% gel using PAGE gel quick preparation kit (10%, CAT: PG112) and transferred to Trans-blot Turbo nitrocellulose membranes (Bio-Rad). Then, incubate the membrane in 10% milk and seal it for 2 hours. For primary antibody–protein hybridization, membranes were probed with the following antibodies at 4°C overnight (Beyotime): LDHA Rabbit Polyclonal Antibody, CAT: AF0216; COX IV Rabbit Polyclonal Antibody, CAT: AF6549; Hexokinase II Rabbit Polyclonal Antibody, CAT: AF7080; SDHB Rabbit Polyclonal Antibody, CAT: AF7956; and PGK1 Rabbit Monoclonal Antibody, CAT: AF1825. Thereafter, secondary anti-mouse (Beyotime CAT: LF101) or anti-rabbit IgG antibodies (Beyotime CAT: LF102) were incubated for 1 h at room temperature. Protein bands were developed with chemiluminescent reagents (Beyotime) and imaged with a Tanon 5200.

2.11. Cell Proliferation Experiment (CCK8)

After 36 hours transfection, the cells were counted and inoculated in a 96-well plate. After 24, 48, and 72 hours of culture, 10 μL CCK8 solution (Beyotime, China) was added to each well. The OD value of each hole was measured.

2.12. Wound Healing Assay

After si-VDAC1 or NC mimic transfection, A549 cells grew for 24 hours. When the cell confluence rate was about 70-80%, the cells were scratch-wounded in a straight line. The cells were photographed after wounding and after 36 hours of incubation. The wound healing rate was calculated as the area of cells that migrate from the wound edge into the wound zone.

2.13. Transwell Invasion Assays

The cell invasion assay was performed in 24-well plates using filters with an 8.0 μm pore size (Corning, USA) coated with 1: 8 dilution Matrigel (BD Bioscience, USA). Six hundred μL medium containing 10% FBS was added to the lower chamber, and A549 cells in serum-free medium were seeded into the upper chambers. After 36 hours of incubation, the cells were fixed with 4% paraformaldehyde and stained with 0.5% crystal violet. Cells that had invaded the lower surface were captured and evaluated in 6 randomly selected areas.

2.14. Statistical Analyses

Unpaired Student’s $t$-test or one-way ANOVA test was used to compare normally distributed data. Nonnormally distributed data were performed using the Mann–Whitney $U$ test or Kruskal-Wallis test. The prognostic nomogram’s construction used R package “rms” and followed by Iasonos’ guide [19]. All statistical tests and visual analyses were performed using the R software (version 4.0.3) or GraphPad Prism 6.0 (GraphPad Software, USA).

3. Results

3.1. GSEA Analysis of Mitochondrial Activity in LUAD

In order to explore the mitochondrial activity in lung cancer, we downloaded several mitochondrial-related pathways from MSigDB, including the complete mito-protein genes, fatty acid metabolism, glycolytic pathways, hypoxia and high expression genes, and oxidative phosphorylation. Next, the changes of these pathways in the four microarrays were analyzed, namely, GSE7670: 31 LUAD vs. 27 normal, GSE18842: 46 LUAD vs. 45 normal, GSE19188: 45 LUAD vs. 65 normal, and GSE31210: 226 LUAD vs. 20 normal. In all four data sets (Figure 1), glycolytic activity in LUAD was the most significantly increased. Mito-protein genes and oxidative phosphorylation are also upregulated. However, fatty acid metabolism and hypoxia are indeed downregulated in LUAD. It shows that lung cancer cells are not hypoxic, and fatty acid metabolism is not vigorous. The results of GSEA indicate that there is a significant correlation between LUAD and mitochondrial activity and provide a theoretical basis for the development of mito-protein genes to predict the prognosis of LUAD patients.

[figure omitted; refer to PDF]

3.2. Relationship between Mito-Protein Genes and Prognosis

We used TCGA-LUAD mRNA-seq data and clinical data to perform univariate Cox regression analysis to examine the relationship between mito-protein genes and patients’ prognosis. As shown in Supplementary Figure S1, we found 199 genes that are significantly related to the prognosis of LUAD patients ($p<0.05$, $\text{HR}<-1$, or $\text{HR}>1$) and passed the bootstrap test. Among them, 90 genes have hazard ratios less than 1, indicating that patients who overexpress these genes have a longer survival time. In comparison, 109 genes with hazard ratios greater than 1 have the opposite meaning.

3.3. Construction of a Prognostic Risk Panel

The 199 robust prognostic genes were subjected to LASSO Cox regression analysis for dimension reduction. The convergence of the regression coefficients is shown in Figure 2(a). Through a random sampling method of 10 cross-validation, we found that the model constructed from 5 genes performed best (Figure 2(b)). Based on the model correlation coefficient, we constructed a mito-protein gene panel to predict the prognosis of LUAD patients.$$\begin{array}{}\text{(1)}& \text{mitoRiskscore}=0.8056\ast {\text{Exp}}_{\text{VDAC}1}-0.1433\ast {\text{Exp}}_{\text{SLC}25\text{A}42}-0.1913\ast {\text{Exp}}_{\text{ABAT}}-0.1083\ast {\text{Exp}}_{\text{IVD}}-0.0255\ast {\text{Exp}}_{\text{AMT}}\end{array}$$

[figures omitted; refer to PDF]

We calculated the mitoRiskscore score for each LUAD patient. Patients were divided into two groups (high-risk group and low-risk group) according to mitoRiskscore by using the median of the cohort as the cut-off value. The distribution of mitoRiskscore and patient survival status is shown in Figure 2(c). Figure 2(d) shows the expression of these 5 genes in the two groups of patients. The Kaplan–Meier curve indicated that the survival rate of patients in the high-risk group was significantly lower (Figure 2(e)).

3.4. Independent Prognostic Value of the MitoRiskscore

Univariate and multivariate Cox regression analyses demonstrated that mitoRiskscore was a strong independent risk factor for the overall survival of LUAD patients (Figure 3(a)). The predictive power of mitoRiskscore ($p<0.001$, $\text{HR}=1.704$, 95% $\text{CI}=1.425$-2.038) was even higher than age and TMN stage. In the ROC analysis of one year, three years, and five years (Figure 3(b)), mitoRiskscore showed higher predictive ability than the TNM stage and age. The accuracy of mitoRiskscore’s judgment of the patient’s survival status reached more than 0.7 in all three ROC analyses. These data indicated that mitoRiskscore is an independent prognostic factor and may predict patient prognosis more accurately than existing clinical parameters.

[figures omitted; refer to PDF]

3.5. External Verification of the MitoRiskscore

A robust prediction panel is universally applicable to other cohorts. In order to verify whether mitoRiskscore applies to other LUAD cohorts, we used external data for verification. We enrolled six cohorts with more than 50 samples each to ensure the reliability of the analysis results (GSE3141, GSE8894, GSE31210, GSE68465, GSE72094, GSE50081). We performed the same procedure to divide patients into high-risk and low-risk groups based on mitoRiskscore. The detailed data of the GEO cohorts performed in the external verification was uploaded in https://github.com/XR-Zhang-group/mitoRiskscore. Surprisingly, in all validation cohorts, patients in the high-risk group showed a higher mortality rate. Except for GSE50081 ($p=0.1004$), the other results were statistically significant (Figure 4). Besides, we also use mitoRiskscore to predict patients’ 1-year, 3-year, and 5-year survival rates. In most ROC analyses, mitoRiskscore’s AUC reached 0.6 or more. The above data showed that mitoRiskscore has specific practical application value in predicting the prognosis of LUAD patients.

[figure omitted; refer to PDF]

3.6. Identify the Metabolic Patterns of Patients with Different MitoRiskscore

In order to determine the metabolism and mitochondrial activity in patients with high mitoRiskscore and patients with low scores, we performed a GSEA analysis (Figure 5(a)). The results showed that fatty acid metabolism, glycolysis, hypoxia, oxidative phosphorylation, ROS pathway, and mito-protein genes were significantly activated in patients with high mitoRiskscore compared with patients with low scores. We performed GSVA enrichment analysis to determine the biological behavior characteristics of patients with different mitoRiskscore (Figure 5(b)). The results showed that LUAD with high mitoRiskscore scores has higher glycometabolism and energy generation, as well as higher proliferation activity and nucleotide repairment. In addition, we found that the energy metabolism pattern of low mitoRiskscore LUAD was similar to that of the adjacent cancer controls. They are not active in glycometabolism and oxidative phosphorylation, but they had more lipid metabolism than LUAD with high mitoRiskscore. It confirmed that mitoRiskscore could distinguish LUAD patients with different metabolic patterns.

[figures omitted; refer to PDF]

3.7. MitoRiskscore Predicts Therapeutic Benefit of Chemotherapy

Adjuvant chemotherapy is the main treatment strategy for LUAD patients after surgery. The effectiveness of chemotherapy drugs is often related to the metabolism of tumor cells. We selected six commonly used chemotherapy drugs for LUAD patients and evaluated the chemotherapy sensitivity of patients with different mitoRiskscore. All LUAD samples were divided into high-risk and low-risk groups based on the mitoRiskscore. 59 paracancerous tissues in the TCGA-LUAD dataset were used as a control group. The differences in chemotherapy sensitivity between the three groups were then compared (Figure 6(a)). For cisplatin, the IC50s of the LUAD groups were higher than that of the control group. For vinorelbine and paclitaxel, there were no significant difference in IC50s between the LUAD groups and the control group. The IC50s of etoposide in the LUAD groups were lower than that in the control group. But at the same time, many scattered points were displayed higher than the control group’s average. It means that etoposide’s effectiveness varies greatly, and some patients are not sensitive to it. Furthermore, the LUAD groups were more sensitive to docetaxel and gemcitabine than the control group. Then, we did correlation analysis on the relationship between the IC50 of these two drugs and mitoRiskscore (Figure 6(b)). The results showed that the higher the mitoRiskscore, the lower the IC50 ($R=-0.39$ and $R=-0.42$). It showed that docetaxel and gemcitabine may be two effective drugs for high mitoRiskscore patients.

[figures omitted; refer to PDF]

3.8. Develop a Prognostic Nomogram Based on MitoRiskscore

In order to improve the accuracy of prognosis and facilitate clinical usage, we developed a nomogram that integrates mitoRiskscore and clinical prognostic factors to predict the 3, 5, and 10-year survival rate of patients (Figure 7(a)). The sum of each factor’s contribution scores can be used to determine the prognosis of the patient. The 3-year and 5-year calibration charts (Figure 7(b)) showed that our nomogram performs well compared to an ideal model. The decision curve analysis shows that the clinical utility of our nomogram greatly overwhelmed the clinical features (Figure 7(c)). It showed that by using mitoRiskscore plus clinical features to predict prognosis, more patients could benefit from it.

[figures omitted; refer to PDF]

3.9. In Vitro Verification of MitoRiskscore

To further verify the authenticity of mitoRiskscore, we conducted in vitro experiments on the five genes in mitoRiskscore. First, we used two lung adenocarcinoma cell lines (A549 and PC-9) and a bronchial epithelial cell line (BEAS-2B) to detect the expression of these five genes (Figure 8(a)). As expected, VDAC1 was highly expressed in the two LUAD cell lines and low in normal lung epithelial cells. The expression patterns of the other four genes were reversed. Given that the expression of VDAC1 accounted for the largest proportion in mitoRiskscore, we further conducted in vitro experiments to verify the function of VDAC1.

[figures omitted; refer to PDF]

In order to explore the effect of VDAC1 on the metabolic ability of lung cancer cells, we measured the extracellular acidification rate (ECAR) to quantify the glycolysis ability in A549 cells knocking down VDAC1 in A549 cell line (Figure 8(b) and 8(c)). The results showed that when the expression of VDAC1 in lung cancer cells was downregulated, both the basic glycolysis level and the maximum glycolysis level decreased significantly. The glycolysis level of tumor cells is increased, suggesting that VDAC1 knockdown may inhibit tumor growth by inhibiting cell metabolism. Next, we tested the effects of knocking down VDAC1 on several key enzymes of glycolysis and oxidative phosphorylation in A549 cells (Figure 8(d)). Our results showed that the protein expression of cytochrome c oxidase IV (COXIV) and lactate dehydrogenase (LDHA) decreased after knocking down VDAC1. In contrast, there was no significant change in the protein expression of other enzymes such as succinate dehydrogenase complex subunits B (SDHB), hexokinase II (HKII), and creatine phosphoglycerol kinase 1 (PGK1).

Based on the expression of VDAC1, we divided the TCGA-LUAD samples into high expression (top 20%) and low expression groups (bottom 20%). GSEA was used to analyze the pathways activated in patients with high expression of VDAC1 compared to patients with low expression (Figure 8(e)). The results showed that in patients with high expression of VDAC1, the pathways that promote cell proliferation, such as MYC, MTOR, and energy metabolism, were significantly activated. CCK-8 assays indicated that the downregulation of VADC1 significantly reduced the proliferation ability of A549 cells (Figure 8(f)). Then, a transwell assay was performed to examine the effect of VDAC1 on the invasion ability of LUAD cells. The results showed that knockdown of VDAC1 significantly inhibited the invasion ability of A549 cells (Figure 8(h)). In contrast, we found through wound healing assay that the silencing of VDAC1 did not affect the migration of A549 cells (Figure 8(g)). Taken together, in vitro experiments confirmed that knocking down VDAC1 attenuated the ability of LUAD cells to glycolysis. At the same time, it can reduce the proliferation and invasion ability of LUAD cells.

4. Discussion

Abnormal metabolism is a very typical feature of cancer cells. Due to the central position of mitochondria in cell material metabolism, it is necessary for us to analyze the mitochondrial state of cancer. The most famous theory about the changes in mitochondria’s energy metabolism in cancer was the Warburg effect put forward in the 20th century. That was, most cancer cells produced energy through glycolysis [20]. With the deepening of research, the Warburg effect was gradually overturned. It has been found that the primary way for cancer cells to produce energy was by aerobic oxidation, and the increase in glycolysis was to create more intermediate metabolites [7]. For example, it has been found that the colon cancer cell line SW620 had increased OXPHOS but glycolysis [8]. Ovarian cancer cells gave priority to OXPHOS to address energy requirements [21]. Besides, increased glycolysis does not always promote the development of cancer. In lung cancer cells, knockout of AIF increased glycolysis and reduced oxidative phosphorylation level but inhibited cancer cells’ growth [22]. In addition to OXPHOS, cancer cells also showed an increased breakdown of fatty acids and amino acids. For example, α-ketoglutarate (an intermediate produced by the decomposition of glutamine) could supplement mitochondrial TCA intermediates [9]. The high metabolism of cancer cell mitochondria produced excessive ROS, promoting the death of normal cells, thereby promoting tumor progression [10]. In addition, the dynamics of cancer mitochondria were also changed. Mitochondria maintain their functions through autophagy, division, and fusion. The activity of mitochondria in cancer increased, and there were more mitochondria in a single tumor cell [23]. Therefore, our research focused on the changes in LUAD mitochondria.

We first analyzed the changes in mito-protein genes and metabolism in LUAD and normal lung tissues. Not surprisingly, both OXPHOS and glycolysis were very active in LUAD. But we also found that lipid metabolism was not active. Because the transcription of mito-protein genes was abundant in LUAD, we tried to find the relationship between them and patients’ prognosis. Among the 1136 mito-protein genes, we found that 199 are significantly related to LUAD patients’ prognosis, which was a large proportion. Through LASSO cox regression, we further screened five genes (VDAC1, SLC25A42, ABAT, IVD, and AMT) and constructed a prognostic panel named mitoRiskscore.

Among the genes of mitoRiskscore, VDAC1 is the most studied. VDAC1 is a transporter protein on the outer mitochondrial membrane (OMM) [24]. The functions of VDAC1 in metabolism and energy homeostasis are reflected by its facilitation of the transport of ions, nucleotides, and other metabolites across the OMM [25–27]. Specifically, VDAC1 is a transporter of cellular metabolites [25], cholesterol [28], lipids [29], and Ca2+ [30]. VDAC1 allows two-way traffic, mediating the entry of metabolites including pyruvate, malate, succinate, and NADH into mitochondria and the exit of newly formed molecules. The importance of VDAC1 in cell energy and metabolism homeostasis is reflected in the findings that downregulation of VDAC1 expression resulted in reduced metabolite exchange between mitochondria and the cytosol and inhibited cell growth, showing VDAC1 as essential for energy production and cell growth [31]. In our current research, we performed GSEA between the VDAC1 high expression group and the low expression group. The results showed that glycolysis, oxidative phosphorylation, MYC targets, E2F targets, and mTOR signaling pathways were all significantly activated in the high expression group. These activated pathways are related to tumor proliferation and metastasis. This result was consistent with our subsequent in vivo experiments. Because VDAC1 plays an important role in material transport in OMM, we believe that knocking down VDAC1 affects cell metabolism by restricting substance transport. At the same time, the expression of some metabolic enzymes decreased reflexively [32, 33]. We also found that knocking down VDAC1 could suppress the glycolysis metabolism, proliferation, and invasion of LUAD cells, which proved its momentous position in the progression of LUAD.

ABAT (4-aminobutyric acid aminotransferase) is responsible for the catabolism of γ-aminobutyric acid (GABA). GABA is an important inhibitory neurotransmitter in the central nervous system. ABAT deficiency phenotypes include psychomotor retardation, hypotension, hyperreflexia, lethargy, refractory epilepsy, and abnormal EEG. The ABAT gene is closely related to neurological diseases, but it is also associated with cancer. Our survival analysis showed that ABAT was a protective gene. The low expression of ABAT in basal cell-like breast cancer (BLBC) promoted the increase of GABA, thereby enabling the development and metastasis of cancer cells by activating the Ca2 ⁺ -NFAT1 axis [34]. In renal clear cell carcinoma, low-expressed ABAT promoted tumor development by boosting cancer cells’ metabolism, which the transcription factor HNF4A may regulate [35]. ABAT was related to the prognosis of patients with liver cancer, pancreatic cancer, and breast cancer [14, 15, 36]. SLC25A42, AMT, and IVD genes have also been found in other prognostic models of lung cancer [37, 38]. However, no studies have found a relationship between AMT and cancers. Our research found that these three genes’ high expression was related to LUAD patients’ prognosis and was a protective factor. Further study of the specific mechanisms of these three genes inhibiting LUAD may provide new directions for LUAD treatment.

MitoRiskscore contains the above five genes, and its ability to predict the prognosis of LUAD patients is higher than that of clinical features (Figure 3). What is more commendable is that we have used a lot of external datasets to verify this conclusion (Figure 4). MitoRiskscore predicted LUAD patients’ prognosis well in most external verifications, indicating that it is a reliable and effective prognostic evaluation tool. Next, we explored the specific biological significance of mitoRiskscore. High mitoRiskscore represents relatively active metabolism and high expression of mitochondrial localization protein genes (Figure 5). Specifically, the glycometabolism and OXPHOS pathways are significantly activated in LUAD with high mitoRiskscore, which indicates that cells have higher energy production and proliferation capabilities. However, we have also observed activation of the p53 pathway and nucleotide repair in high mitoRiskscore tumors, which may be feedback to the cell’s high DNA replication and proliferation state. On the other hand, we found that lipid metabolism is relatively weak in high mitoRiskscore tumors. Fatty acid synthesis occurs in the cytoplasm. They are derived from acetyl-coenzyme A and are mainly provided by the tricarboxylic acid (TCA) cycle [39]. Based on the above results, we speculate that there is an upregulation of carbohydrate metabolism in lung cancer with highly activated mitochondria, which generates much energy through OXPHOS. The energy and acetyl-coenzyme A produced by these TCAs are not used to synthesize lipids. They are used by cells to maintain a high proliferation state and are associated with poor postoperative survival rates for patients. In short, mitoRiskscore can evaluate the energy metabolism patterns of the samples.

In this study, the R software package “pRRophetic” and mRNA expression profile were used to predict patients’ six drug sensitivity and controls. This method required more clinical trial verification and only gave us a preliminary hint (Figure 6). The results showed no difference in the sensitivity of each group to vinorelbine, and even the control group was more sensitive to cisplatin. Different samples in the same group had different sensitivity to etoposide and paclitaxel. The high mitoRiskscore group’s sensitivity to docetaxel and gemcitabine was higher than that of the low mitoRiskscore group, and both groups were higher than the control group. Furthermore, the sensitivity of patients to docetaxel and gemcitabine was significantly negatively correlated with the mitoRiskscore. It showed that the use of docetaxel and gemcitabine to treat LUAD with active mitochondrial function might achieve better results.

We constructed a nomogram containing mitoRiskscore and clinical features, making the results of this study easier to apply to the clinic (Figure 7). Nomogram is a popular tool for predicting the prognosis of cancer. It uses statistical methods to combine patients’ parameters to predict their prognosis. Due to a combination of factors, the accuracy of nomogram is more accurate than that of a single stage of cancer [19, 40]. The calibration chart analysis and the decision curve analysis showed that the nomogram had better prediction accuracy and can benefit more patients.

One of the innovations of our research is that we developed a LUAD prognosis predictive model based on mitochondrial localized genes for the first time. Since mitochondria are involved in most of the metabolic pathways, our model is more representative than the general metabolic model. In addition, our model screened several genes that have not been studied in the field of lung cancer but are closely related to the prognosis of LUAD patients, such as SLC25A42, AMT, and IVD, which provides a basis for future research on the mechanism of lung cancer.

There were still some limitations in this study. Firstly, the transcriptome data used for panel construction comes from sequencing. Although we had verified the function of mitoRiskscore using chip data, PCR data has not been used yet. Secondly, since this study used gene expression data as categorical variables to be input into Cox regression, it is necessary to further determine the optimal cut-off value. Thirdly, because the study was retrospective, the patient population is heterogeneous, and the results may be biased. The mitoRiskscore and nomogram obtained in this study require more clinical studies to verify their effectiveness.

5. Conclusions

In conclusion, the mitoRiskscore gene panel is an effective tool for the survival prediction of LUAD patients. It can also assess the metabolic status of tumors and assist with clinical chemotherapy. However, further clinical trials are needed to verify our findings.

Authors’ Contributions

Weifeng Chen, Jingyao Wang, and Qiumei Zhao contributed equally to this work.

References

[1] F. Bray, J. Ferlay, I. Soerjomataram, R. L. Siegel, L. A. Torre, A. Jemal, "Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries," CA: A Cancer Journal for Clinicians, vol. 68 no. 6, pp. 394-424, DOI: 10.3322/caac.21492, 2018.

[2] B. C. Bade, C. S. Dela Cruz, "Lung cancer 2020: epidemiology, etiology, and prevention," Clinics in Chest Medicine, vol. 41 no. 1,DOI: 10.1016/j.ccm.2019.10.001, 2020.

[3] N. Duma, R. Santana-Davila, J. R. Molina, "Non-small cell lung cancer: epidemiology, screening, diagnosis, and treatment," Mayo Clinic Proceedings, vol. 94 no. 8, pp. 1623-1640, DOI: 10.1016/j.mayocp.2019.01.013, 2019.

[4] J. A. Barta, C. A. Powell, J. P. Wisnivesky, "Global epidemiology of lung cancer," Annals of Medicine and Surgery, vol. 85 no. 1,DOI: 10.5334/aogh.2419, 2019.

[5] P. Goldstraw, K. Chansky, J. Crowley, R. Rami-Porta, H. Asamura, W. E. Eberhardt, A. G. Nicholson, P. Groome, A. Mitchell, V. Bolejack, P. Goldstraw, R. Rami-Porta, H. Asamura, D. Ball, D. G. Beer, R. Beyruti, V. Bolejack, K. Chansky, J. Crowley, F. Detterbeck, W. E. Erich Eberhardt, J. Edwards, F. Galateau-Sallé, D. Giroux, F. Gleeson, P. Groome, J. Huang, C. Kennedy, J. Kim, Y. T. Kim, L. Kingsbury, H. Kondo, M. Krasnik, K. Kubota, A. Lerut, G. Lyons, M. Marino, E. M. Marom, J. van Meerbeeck, A. Mitchell, T. Nakano, A. G. Nicholson, A. Nowak, M. Peake, T. Rice, K. Rosenzweig, E. Ruffini, V. Rusch, N. Saijo, P. van Schil, J. P. Sculier, L. Shemanski, K. Stratton, K. Suzuki, Y. Tachimori, C. F. Thomas, W. Travis, M. S. Tsao, A. Turrisi, J. Vansteenkiste, H. Watanabe, Y. L. Wu, P. Baas, J. Erasmus, S. Hasegawa, K. Inai, K. Kernstine, H. Kindler, L. Krug, K. Nackaerts, H. Pass, D. Rice, C. Falkson, P. L. Filosso, G. Giaccone, K. Kondo, M. Lucchi, M. Okumura, E. Blackstone, F. Abad Cavaco, E. Ansótegui Barrera, I. Parente Lamelas, A. Arnau Obrer, R. Guijarro Jorge, D. Ball, G. K. Bascom, M. A. González Castro, M. G. Blum, D. Chimondeguy, V. Cvijanovic, S. Defranchi, B. de Olaiz Navarro, I. Escobar Campuzano, I. Macía Vidueira, E. Fernández Araujo, F. Andreo García, K. M. Fong, G. Francisco Corral, S. Cerezo González, J. Freixinet Gilart, L. García Arangüena, S. García Barajas, P. Girard, T. Goksel, M. T. González Budiño, G. González Casaurrán, J. A. Gullón Blanco, J. Hernández Hernández, H. Hernández Rodríguez, J. Herrero Collantes, M. Iglesias Heras, J. M. Izquierdo Elena, E. Jakobsen, S. Kostas, P. León Atance, A. Núñez Ares, M. Liao, M. Losanovscky, G. Lyons, R. Magaroles, L. de Esteban Júlvez, M. Mariñán Gorospe, B. McCaughan, C. Kennedy, R. Melchor Íñiguez, L. Miravet Sorribes, S. Naranjo Gozalo, C. Álvarez de Arriba, M. Núñez Delgado, J. Padilla Alarcón, J. C. Peñalver Cuesta, J. S. Park, H. Pass, M. J. Pavón Fernández, M. Rosenberg, E. Ruffini, V. Rusch, J. Sánchez de Cos Escuín, A. Saura Vinuesa, M. Serra Mitjans, T. E. Strand, D. Subotic, S. Swisher, R. Terra, C. Thomas, K. Tournoy, P. van Schil, M. Velasquez, Y. L. Wu, K. Yokoi, "The IASLC lung cancer staging project: proposals for revision of the TNM stage groupings in the forthcoming (eighth) edition of the TNM classification for lung cancer," Journal of Thoracic Oncology, vol. 11 no. 1, pp. 39-51, DOI: 10.1016/j.jtho.2015.09.009, 2016.

[6] P. Ghosh, C. Vidal, S. Dey, L. Zhang, "Mitochondria targeting as an effective strategy for cancer therapy," International Journal of Molecular Sciences, vol. 21 no. 9,DOI: 10.3390/ijms21093363, 2020.

[7] T. M. Ashton, W. G. McKenna, L. A. Kunz-Schughart, G. S. Higgins, "Oxidative phosphorylation as an emerging target in cancer therapy," Clinical Cancer Research, vol. 24 no. 11, pp. 2482-2490, DOI: 10.1158/1078-0432.CCR-17-3070, 2018.

[8] C. S. Lin, L. T. Liu, L. H. Ou, S. C. Pan, C. I. Lin, Y. H. Wei, "Role of mitochondrial function in the invasiveness of human colon cancer cells," Oncology Reports, vol. 39 no. 1, pp. 316-330, DOI: 10.3892/or.2017.6087, 2018.

[9] E. Currie, A. Schulze, R. Zechner, T. C. Walther, R. V. Farese, "Cellular fatty acid metabolism and cancer," Cell Metabolism, vol. 18 no. 2, pp. 153-161, DOI: 10.1016/j.cmet.2013.05.017, 2013.

[10] M. Valko, C. J. Rhodes, J. Moncol, M. Izakovic, M. Mazur, "Free radicals, metals and antioxidants in oxidative stress-induced cancer," Chemico-Biological Interactions, vol. 160 no. 1,DOI: 10.1016/j.cbi.2005.12.009, 2006.

[11] F. Sanchez-Vega, M. Mina, J. Armenia, W. K. Chatila, A. Luna, K. C. La, S. Dimitriadoy, D. L. Liu, H. S. Kantheti, S. Saghafinia, D. Chakravarty, F. Daian, Q. Gao, M. H. Bailey, W. W. Liang, S. M. Foltz, I. Shmulevich, L. Ding, Z. Heins, A. Ochoa, B. Gross, J. Gao, H. Zhang, R. Kundra, C. Kandoth, I. Bahceci, L. Dervishi, U. Dogrusoz, W. Zhou, H. Shen, P. W. Laird, G. P. Way, C. S. Greene, H. Liang, Y. Xiao, C. Wang, A. Iavarone, A. H. Berger, T. G. Bivona, A. J. Lazar, G. D. Hammer, T. Giordano, L. N. Kwong, G. McArthur, C. Huang, A. D. Tward, M. J. Frederick, F. McCormick, M. Meyerson, E. M. van Allen, A. D. Cherniack, G. Ciriello, C. Sander, N. Schultz, S. J. Caesar-Johnson, J. A. Demchok, I. Felau, M. Kasapi, M. L. Ferguson, C. M. Hutter, H. J. Sofia, R. Tarnuzzer, Z. Wang, L. Yang, J. C. Zenklusen, J. Zhang, S. Chudamani, J. Liu, L. Lolla, R. Naresh, T. Pihl, Q. Sun, Y. Wan, Y. Wu, J. Cho, T. DeFreitas, S. Frazer, N. Gehlenborg, G. Getz, D. I. Heiman, J. Kim, M. S. Lawrence, P. Lin, S. Meier, M. S. Noble, G. Saksena, D. Voet, H. Zhang, B. Bernard, N. Chambwe, V. Dhankani, T. Knijnenburg, R. Kramer, K. Leinonen, Y. Liu, M. Miller, S. Reynolds, I. Shmulevich, V. Thorsson, W. Zhang, R. Akbani, B. M. Broom, A. M. Hegde, Z. Ju, R. S. Kanchi, A. Korkut, J. Li, H. Liang, S. Ling, W. Liu, Y. Lu, G. B. Mills, K. S. Ng, A. Rao, M. Ryan, J. Wang, J. N. Weinstein, J. Zhang, A. Abeshouse, J. Armenia, D. Chakravarty, W. K. Chatila, I. de Bruijn, J. Gao, B. E. Gross, "Oncogenic Signaling Pathways in The Cancer Genome Atlas," Cell, vol. 173 no. 2, pp. 321-337.e10, DOI: 10.1016/j.cell.2018.03.035, 2018.

[12] G. Yu, L. G. Wang, Y. Han, Q. Y. He, "clusterProfiler: an R package for comparing biological themes among gene clusters," OMICS, vol. 16 no. 5, pp. 284-287, DOI: 10.1089/omi.2011.0118, 2012.

[13] A. Subramanian, P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A. Gillette, A. Paulovich, S. L. Pomeroy, T. R. Golub, E. S. Lander, J. P. Mesirov, "Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles," Proceedings of the National Academy of Sciences of the United States of America, vol. 102 no. 43, pp. 15545-15550, DOI: 10.1073/pnas.0506580102, 2005.

[14] K. Nagashima, Y. Sato, "Information criteria for Firth’s penalized partial likelihood approach in Cox regression models," Statistics in Medicine, vol. 36 no. 21, pp. 3422-3436, DOI: 10.1002/sim.7368, 2017.

[15] R. Zhou, D. Zeng, J. Zhang, H. Sun, J. Wu, N. Li, L. Liang, M. Shi, J. Bin, Y. Liao, N. Huang, W. Liao, "A robust panel based on tumour microenvironment genes for prognostic prediction and tailoring therapies in stage I-III colon cancer," eBioMedicine, vol. 42, pp. 420-430, DOI: 10.1016/j.ebiom.2019.03.043, 2019.

[16] P. J. Heagerty, T. Lumley, M. S. Pepe, "Time-dependent ROC curves for censored survival data and a diagnostic marker," Biometrics, vol. 56 no. 2, pp. 337-344, DOI: 10.1111/j.0006-341X.2000.00337.x, 2000.

[17] S. Hanzelmann, R. Castelo, J. Guinney, "GSVA: gene set variation analysis for microarray and RNA-seq data," BMC Bioinformatics, vol. 14 no. 1,DOI: 10.1186/1471-2105-14-7, 2013.

[18] P. Geeleher, N. J. Cox, R. S. Huang, "Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines," Genome Biology, vol. 15 no. 3,DOI: 10.1186/gb-2014-15-3-r47, 2014.

[19] A. Iasonos, D. Schrag, G. V. Raj, K. S. Panageas, "How to build and interpret a nomogram for cancer prognosis," Journal of Clinical Oncology, vol. 26 no. 8, pp. 1364-1370, DOI: 10.1200/JCO.2007.12.9791, 2008.

[20] W. H. Koppenol, P. L. Bounds, C. V. Dang, "Otto Warburg's contributions to current concepts of cancer metabolism," Nature Reviews. Cancer, vol. 11 no. 5, pp. 325-337, DOI: 10.1038/nrc3038, 2011.

[21] A. P. Nayak, A. Kapur, L. Barroilhet, M. S. Patankar, "Oxidative phosphorylation: a target for novel therapeutic strategies against ovarian cancer," Cancers, vol. 10 no. 9,DOI: 10.3390/cancers10090337, 2018.

[22] S. Rao, L. Mondragón, B. Pranjic, T. Hanada, G. Stoll, T. Köcher, P. Zhang, A. Jais, A. Lercher, A. Bergthaler, D. Schramek, K. Haigh, V. Sica, M. Leduc, N. Modjtahedi, T. P. Pai, M. Onji, I. Uribesalgo, R. Hanada, I. Kozieradzki, R. Koglgruber, S. J. Cronin, Z. She, F. Quehenberger, H. Popper, L. Kenner, J. J. Haigh, O. Kepp, M. Rak, K. Cai, G. Kroemer, J. M. Penninger, "AIF-regulated oxidative phosphorylation supports lung cancer development," Cell Research, vol. 29 no. 7, pp. 579-591, DOI: 10.1038/s41422-019-0181-4, 2019.

[23] G. R. Anderson, S. E. Wardell, M. Cakir, C. Yip, Y. R. Ahn, M. Ali, A. P. Yllanes, C. A. Chao, D. P. McDonnell, K. C. Wood, "Dysregulation of mitochondrial dynamics proteins are a targetable feature of human tumors," Nature Communications, vol. 9 no. 1,DOI: 10.1038/s41467-018-04033-x, 2018.

[24] R. Benz, "Permeation of hydrophilic solutes through mitochondrial outer membranes: review on mitochondrial porins," Biochimica et Biophysica Acta, vol. 1197 no. 2, pp. 167-196, DOI: 10.1016/0304-4157(94)90004-3, 1994.

[25] V. Shoshan-Barmatz, V. De Pinto, M. Zweckstetter, Z. Raviv, N. Keinan, N. Arbel, "VDAC, a multi-functional mitochondrial protein regulating cell life and death," Molecular Aspects of Medicine, vol. 31 no. 3, pp. 227-285, DOI: 10.1016/j.mam.2010.03.002, 2010.

[26] M. Colombini, "VDAC structure, selectivity, and dynamics," Biochimica et Biophysica Acta (BBA) - Biomembranes, vol. 1818 no. 6, pp. 1457-1465, DOI: 10.1016/j.bbamem.2011.12.026, 2012.

[27] B. Shanmugavadivu, H. J. Apell, T. Meins, K. Zeth, J. H. Kleinschmidt, "Correct Folding of the β -Barrel of the Human Membrane Protein VDAC Requires a Lipid Bilayer," Journal of Molecular Biology, vol. 368 no. 1, pp. 66-78, DOI: 10.1016/j.jmb.2007.01.066, 2007.

[28] M. B. Rone, J. Fan, V. Papadopoulos, "Cholesterol transport in steroid biosynthesis: role of protein-protein interactions and implications in disease states," Biochimica et Biophysica Acta, vol. 1791 no. 7, pp. 646-658, DOI: 10.1016/j.bbalip.2009.03.001, 2009.

[29] C. Martel, M. Allouche, D. D. Esposti, E. Fanelli, C. Boursier, C. Henry, J. Chopineau, G. Calamita, G. Kroemer, A. Lemoine, C. Brenner, "Glycogen synthase kinase 3-mediated voltage-dependent anion channel phosphorylation controls outer mitochondrial membrane permeability during lipid accumulation," Hepatology, vol. 57 no. 1, pp. 93-102, DOI: 10.1002/hep.25967, 2013.

[30] H. Huang, X. Hu, C. O. Eno, G. Zhao, C. Li, C. White, "An Interaction between Bcl-x L and the Voltage-dependent Anion Channel (VDAC) Promotes Mitochondrial Ca 2+ Uptake," The Journal of Biological Chemistry, vol. 288 no. 27, pp. 19870-19881, DOI: 10.1074/jbc.M112.448290, 2013.

[31] S. Abu-Hamad, S. Sivan, V. Shoshan-Barmatz, "The expression level of the voltage-dependent anion channel controls life and death of the cell," Proceedings of the National Academy of Sciences of the United States of America, vol. 103 no. 15, pp. 5787-5792, DOI: 10.1073/pnas.0600103103, 2006.

[32] E. Zerbib, T. Arif, A. Shteinfer-Kuzmine, V. Chalifa-Caspi, V. Shoshan-Barmatz, "VDAC1 silencing in cancer cells leads to metabolic reprogramming that modulates tumor microenvironment," Cancers, vol. 13 no. 11,DOI: 10.3390/cancers13112850, 2021.

[33] T. Arif, Y. Krelin, I. Nakdimon, D. Benharroch, A. Paul, D. Dadon-Klein, V. Shoshan-Barmatz, "VDAC1 is a molecular target in glioblastoma, with its depletion leading to reprogrammed metabolism and reversed oncogenic properties," Neuro-Oncology, vol. 19 no. 7, pp. 951-964, DOI: 10.1093/neuonc/now297, 2017.

[34] X. Chen, Q. Cao, R. Liao, X. Wu, S. Xun, J. Huang, C. Dong, "Loss of ABAT-mediated GABAergic system promotes basal-like breast cancer progression by activating Ca2+-NFAT1 axis," Theranostics., vol. 9 no. 1, pp. 34-47, DOI: 10.7150/thno.29407, 2019.

[35] J. Lu, Z. Chen, H. Zhao, H. Dong, L. Zhu, Y. Zhang, J. Wang, H. Zhu, Q. Cui, C. Qi, S. Wang, S. Chen, J. Shao, "ABAT and ALDH6A1, regulated by transcription factor HNF4A, suppress tumorigenic capability in clear cell renal cell carcinoma," Journal of Translational Medicine, vol. 18 no. 1,DOI: 10.1186/s12967-020-02268-1, 2020.

[36] J. Budczies, S. F. Brockmöller, B. M. Müller, D. K. Barupal, C. Richter-Ehrenstein, A. Kleine-Tebbe, J. L. Griffin, M. Orešič, M. Dietel, C. Denkert, O. Fiehn, "Comparative metabolomics of estrogen receptor positive and estrogen receptor negative breast cancer: alterations in glutamine and beta-alanine metabolism," Journal of Proteomics, vol. 94, pp. 279-288, DOI: 10.1016/j.jprot.2013.10.002, 2013.

[37] M. Zhang, L. Sun, Y. Ru, S. Zhang, J. Miao, P. Guo, J. Lv, F. Guo, B. Liu, "A risk score system based on DNA methylation levels and a nomogram survival model for lung squamous cell carcinoma," International Journal of Molecular Medicine, vol. 46 no. 1, pp. 252-264, DOI: 10.3892/ijmm.2020.4590, 2020.

[38] K. Zhao, Z. Li, H. Tian, "Twenty-gene-based prognostic model predicts lung adenocarcinoma survival," Oncotargets and Therapy, vol. 11, pp. 3415-3424, DOI: 10.2147/OTT.S158638, 2018.

[39] M. Merino Salvador, M. Gómez de Cedrón, J. Moreno Rubio, S. Falagán Martínez, R. Sánchez Martínez, E. Casado, A. Ramírez de Molina, M. Sereno, "Lipid metabolism and lung cancer," Critical Reviews in Oncology/Hematology, vol. 112, pp. 31-40, DOI: 10.1016/j.critrevonc.2017.02.001, 2017.

[40] Y. W. Won, J. Joo, T. Yun, G. K. Lee, J. Y. Han, H. T. Kim, J. S. Lee, M. S. Kim, J. M. Lee, H. S. Lee, J. I. Zo, S. Kim, "A nomogram to predict brain metastasis as the first relapse in curatively resected non-small cell lung cancer patients," Lung Cancer, vol. 88 no. 2, pp. 201-207, DOI: 10.1016/j.lungcan.2015.02.006, 2015.