We downloaded the transcriptomic and corresponding clinical data of PC patients from the following four omics-related platforms: The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/), ArrayExpress, Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), and International Cancer Genome Consortium (ICGC, https://dcc.icgc.org/) databases. In this analysis, only PC patients who had full follow-up data were included. Data from RNA seq (values expressed as FPKM) were translated into TPM for the TCGA-PC cohort to make them more similar to microarray results. To avoid the batch effect, we use the ComBat function in the SVA package to correct the data from different platforms (41). Overall, we complied the expression profile and clinical information of 930 PC patients consisting of TCGA-PC, GSE62452, GSE28735, GSE57495, ICGC-AU, ICGC-CA, and MTAB-6134 datasets (16–19).

We also downloaded and complied the transcriptomic data of PNET patients from the following three cohorts: GSE73338, GSE98894, and ICGC-PAEN-AU cohorts (20). Similar de-batch technologies were also employed to correct the data from different platforms. Overall, a total of 226 patients with PNET were obtained for the subsequent analysis.

Subsequently, 1136 MTGs were obtained from the MitoCarta3.0 database (http://www.broadinstitute.org/mitocarta), and using the query “oxidative stress,” the GeneCard database (https://www.genecards.org/) provided information on 9512 human genes associated with OS. Finally, 469 MTGs-OS were preserved for further analysis after taking the intersection.

Overall, 178 PC and 171 normal pancreas tissues in the GTEx and TCGA cohorts were obtained for performing differential expression analysis. Subsequently, prognosis-related analysis was conducted on 930 PC samples with complete follow-up records. The ‘limma’ tool in R was employed to evaluate differentially expressed MTGs-OS between PC and normal pancreas tissues from the GTEx and TCGA cohorts. To ascertain the predictive capabilities of MTGs-OS, we applied the “survival” and “survminer” R packages. After the intersection, clustering analyses and model development were carried out using the obtained differentially expressed MTGs-OS with predictive significance.

To summarize the molecular traits of the above MTGs-OS in a variety of human cancers, pan-cancer analysis integrating genomics, transcriptomics, and clinical information was subsequently carried out (21, 22). First, we downloaded and compiled the raw data and clinical records of pan-cancer cohorts from the XENA website. We systematically studied whether the differentially expressed genes in PC and para-cancerous tissues had similar expression characteristics in other malignant tumors. In addition, we investigated the prognostic performances, mutation traits, and methylation levels of MTGs-OS in a series of human cancers. Finally, the ssGSEA algorithm was utilized to derive the MTGs scores of each cancer patient. The GSEA method was then applied to predict the discrepancies in pathway activities between patients with high- and low- MTGs scores (23).

The intrinsic molecular heterogeneity of patients with PNET was identified by a method similar to that of previous cluster analysis. Firstly, the MTGs scores of each PNET patient was evaluated based on the GSVA package, and the MTGs scores were combined with the expression profile of related genes. Based on the MTGs scores and gene expression profile, cluster analysis was carried out in patients with PNET. Similarly, 226 patients with PNET were divided into C2 and C3 subtypes, and the MTGs scores of different subtypes were compared. The two subtypes with no difference in MTGs scores were further merged to get S1 and S2 PNET genotypes, and the MTGs scores between S1 and S2 were compared again. The classical metabolic pathway score and immune pathway score of each PNET patient were calculated based on GSVA package to characterize the inherent molecular characteristics of different PNET subtypes.

The detailed analytic workflow of this research was displayed in Figure 1 . Through the classical bioinformatics difference analysis, we found 469 MTGs-OS in PC ( Figure 2A ). Among 469 MTGs-OS, 182 differentially expressed MTGs-OS were discovered using the lima program in R ( Supplementary Table 1 ). More importantly, 137 MTGs-OS were detected to be remarkably associated with PC prognoses ( Supplementary Table 2 ). After the intersection, 56 differentially expressed MTGs-OS with prognostic values were obtained ( Figure 2B ). The differential expression profiles of above 56 candidate genes were demonstrated in Figure 2C . The prognostic forest plot was also showed in Figure 2D . Notably, most MTGs-OS showed obvious down-regulated expression tendency in PC tissues and were closely associated with favorable prognoses, suggesting their protective roles of MTGs-OS in PC. In addition, the co-expressed relationship among these 56 molecules was intensively explored ( Figure 2E ).

Then, to examine the molecular traits of the 56 MTGs-OS in other human cancers, pan-cancer analyses were carried out based on genomics, transcriptomics, and clinical information. First, we researched the expression characteristics of 56 MTGs-OS in other malignant tumors ( Figure 3A ). Most of the 56 MTGs-OS were expressed differentially between tumors and adjacent normal samples in diverse cancers, including CHOL, COAD, KIRC, UCEC, and LUSC. Besides, the levels of FKBP10, PIF1, and PMAIP1 in nearly all tumor tissues were elevated in contrast with those in paired normal samples; conversely, in almost every case, tumor tissues had lower expression of PDK4, ALDH1L1, ACACB, ACSF2, EPHX2, GATM, ACAT1, and ETFDH than their matched normal tissues. The association of gene expression with patient survival time was then evaluated, and the predictive capabilities of the 56 MTGs-OS in pan-cancers were investigated. As depicted in Figure 3B , we performed a univariate cox regression analysis to determine the MTGs-OS that served as risk factors (HR>1 and p<0.05) and those with protective function (HR<1 and p<0.05). We discovered that certain MTGs-OS are risk factors in part of malignancies, including UVM, ACC, LIHC, KICH, and LGG; however, similar to the protective function of MTGs-OS in PC, these genes also had protective properties in KIRC, MESO, and KIRP. To investigate the mutational characteristics of the 56 MTGs-OS in tumors, the proportion of CNV was studied, and the outcomes showed that CNV happened at high frequency (nearly within the range of 10 to 60%) in the diverse types of cancer ( Figure 3C ). Then, we discovered that BLCA, BRCA, COAD, LUAD, SKCM, STAD, and UCEC, all exhibited a higher frequency of SNVs, especially for UCEC; conversely, the SNV frequency in CHOL, ACC, DLBC, UVM, THYM, THCA, TGCT, PCPG, MESO, LAML, and KICH was low, especially for CHOL and UVM. Interestingly, we found that the SNV of ACACB, ALDH1L1, and GLDC was in significantly high frequencies than that of the other MTGs-OS, as illustrated by the heatmap and waterfall diagram of SNV ( Figures 3D, E ). Notably, the methylation of promoters and abnormal DNA methylation might influence gene expression (27). Therefore, we discovered that MTGs-OS had complicated methylation patterns in different types of cancer. As a whole, the changes in methylation levels of MTGs-OS were not significant in most human cancers, including PC, according to the result ( Figure 3F ). EPHX2, ACSF2, CYP27A1, ACSS1, and ALDH1L1 showed hypermethylation in several types of human cancer; on the contrary, AKR1B10, POLG2, NDUFB8, ACACB, and MRPS22 consistently showed hypomethylation in several types of human cancer.

In consideration of the function and impact of MTGs-OS in PC, the enrichment scores of MTGs-OS were initially evaluated via ssGSEA for a sum of 930 patients. The 930 different PC samples were divided up into three unique categories (Cluster1, Cluster C2, and Cluster C3) as per the concentrations of mRNA that MTGs-OS expresses ( Figure 4A ). C1, C2, and C3 comprised patients with normal, active, and inactive MTGs-OS, respectively. The violin plot illustrates that C2 has the highest enrichment score out of the three clusters, followed by C1 and C3 ( Figure 4B ). Subsequently, by plotting the survival curves for the three categories, we assessed the reliability of the clustering technique. As for the 3 clusters, we found that patients in C3 exhibited the lowest overall survival rates whereas those in C2 exhibited the highest overall survival rates ( Figure 4C ), which revealed that the MTGs-OS score was a protective indicator. The result in the part was in accordance with the results of previous univariate Cox regression analysis. Overall, based on these findings, we conclude that our PC classification technique is accurate, credible, and scientific.

At present, the treatment status of pancreatic cancer is still grim. The new concept of oncology has been instrumental in the rapid advancement of cancer detection and therapy in recent years. Targeted therapy and immunotherapy have advanced the comprehensive diagnostic testing and management of PC beyond the reach of standard therapies like surgery, chemotherapy, and radiation treatments. The foundation of molecular targeted therapy is the knowledge of the molecular biology of neoplasm. Key molecules associated with cancer were identified, which might be used as therapeutic targets. According to the characteristics of the therapeutic targets, using highly precise medications and agents to treat cancer is promising. Based on previous research reports, mitochondrial-related OS participates in the onset and progression of tumors. We also discovered that the MTGs-OS score was a marker that might be protective for those with PC. Therefore, we evaluated the connection between MTGs-OS and the effectiveness of targeted pharmaceuticals already available for treating PC. Twelve drugs were chosen, including mTOR inhibitor (Rapamycin), AMPK activator (AICAR), MEK inhibitors (Selumetinib) plant anticancer drugs (cisplatin), pyrimidine antineoplastic drugs (gemcitabine), PI3Kβ inhibitors (AZD6482), tyrosine kinase inhibitor (ponatinib), IGF-1R/IR inhibitor (BMS.754807, linsitinib), metformin and the other chemotherapeutic drugs (doxorubicin, bleomycin). We applied the “pRRophetic” program to make a comparison of the IC50 values of the 12 drugs across these three molecular subgroups. We compiled a list of MTGs-OS clusters and their associated medication sensitivities: Rapamycin: C2 > C1 > C3; AICAR: C3 > C2 > C1; AZD6244(selumetinib): C1 > C2 > C3; cisplatin: C3> C1 > C2; gemcitabine: C3> C1 > C2; AZD6482(PI3Kβ inhibitors): C3> C1 > C2; AP.24534 (ponatinib): C3> C1 > C2; BMS.754807 (IGF-1R/IR inhibitor): C1> C3 > C2; OSI.906(linsitinib): C3 > C2 > C1; doxorubicin: C3> C1 > C2; bleomycin: C3>C1>C2 ( Figures 5A–L ).

To examine the possible regulation mechanism of the MTGs-OS in PC, we evaluated the expression profiles of tumor suppressor genes and oncogenes in the 3 clusters, as illustrated by the heatmap in Figure 6A . We found that the expression levels of most of the oncogenes were much higher in C3 compared with C1 and C2, such as STAB1, TENM3, PIK3CA, ZFHX4, ADAMTS12, GLI3, HMCN1 and so on. On the contrary, C3 was shown to have a lower expression of tumor suppressor genes including RNF43, in contrast with C1 and C2. The stimulation of oncogenes and inhibition of suppressor genes probably resulted in the worst prognosis in C3. Interestingly, as was the case in the oncogenes, we discovered that the expression patterns of several tumor suppressor genes, such as PTEN, were highest in C3, but elevated in C1 and C2. We also examined the levels of activity in commonly recognized metabolic and immune-related pathways across the three clusters and computed the enrichment scores by ssGSEA ( Figures 6B, C ). Interestingly, various metabolism and immune-related pathways exhibited significant differences in the 3 clusters. We found that most of the protective metabolism pathways including KEGG_GLYCINE_SERINE_AND_THREONINE_METABOLISM were more inactive in C3 than that in C1 and C2. Also, almost all the immune-associated pathways were more active in C3 than in C1 and C2, such as KEGG_P53_SIGNALING_PATHWAY and KEGG_PROGESTERONE_MEDIATED_OOCYTE_MATURATION, which suggested that the degree of malignancy for PC was highest in the C3 leading to the worst prognosis. Additionally, the aforementioned findings established that MTGs-OS were protective factors for PC since they were strongly linked to metabolic reprogramming and immune-related pathways.

It has been well-recognized that the tumor microenvironment (TME) performs an instrumental function in PC development (28, 29). It has been postulated that the TME, particularly the features of tumor-infiltrating immune cells (TIICs), might be used as treatment targets for cancer (30). The levels of ROS in the TME are critical for regulating the levels of immune cells and their functions in the TME, particularly the activation, expansion, and effector function of T lymphocytes (31, 32). Within the scope of this research, we delved deeper into the link between MTGs-OS and the infiltration levels of immune cells as well as the corresponding immune functions ( Figure 6D ). The expression of following genes, NDUFA7, NDUFA3, MMAB, MCRIP2, and COQ4, were negatively correlated with immunocyte infiltration and immune-related functions. Conversely, UCP2, ACAT1, CYP27A1, DLAT, FTH1, MPV17, PMAIP1, PPTC7, and GLDC, were positively linked to the infiltration levels of immune cells. Notably, ICGs contributed to anti-tumor immunity by modulating the function of effector T lymphocytes. Therefore, we studied the relationship between ICGs in the 3 clusters ( Figure 6E ). We discovered that the MTGs-OS-inactive cluster was correlated with overexpression of ICGs, illustrating that MTGs-OS-inactive patients had the weakest anti-tumor immunity and a poor prognosis. Subsequently, we assessed the association of MTGs-OS scores with immune cell infiltration ( Figures 6F–K ). Most infiltrating immune cells, comprising neutrophils, APC co-inhibition, CCR, T cell co-inhibition, and macrophages, were shown to have a negative link to MTGs-OS scores. The “ESTIMATE” R program was also used to examine the TIME of the three subtypes. The immune, stromal, and ESTIMATE scores, as well as the tumor purity, were depicted in Figures 7A–D . We discovered that C2 had considerably decreased stromal, immune, and ESTIMATE scores in contrast to C3. On the contrary, although the tumor purity in C2 was substantially increased in contrast with that in C3, the patients in C2 had a better prognosis. At the same time, we evaluated immune cell responses or cellular components in the three PC subtypes according to the EPIC, MCPCOUNTER, CIBERSORT-ABS, XCELL, CIBERSORT, and TIMER algorithms ( Figure 7E ). Overall, the cellular immune responses or cellular components were significantly different among three clusters, which suggested that there was a close link between the infiltration status of most immune cells and the MTGs-OS score.

Previous 56 differentially expressed MTGs-OS with prognostic values were used for conducting the following LASSO-Cox regression analysis. LASSO regression analysis identified 18 candidate molecules and multivariate Cox regression analysis determined 11 hub MTGs-OS (i.e., COQ4, MCRIP2, SOD1, NDUFB8, MRPL50, BIK, MRPL14, RFK, NMNAT3, BNIP3L, MRPS5) for constructing the prognostic model ( Supplementary Figures 1A–C ). To additionally validate the fundamental function of the selected modeling genes in the onset and progression of PC, we investigated the association of the eleven model genes with PAAD stage and metastasis based on the GEPIA and BEST platforms ( Figure 8 ).

Subsequently, 11 hub MTGs-OS’s expression levels were verified by RT-PCR in four cell lines ( Figure 9 ). We observed that the mRNA expression trend of the 11 hub MTGs-OS was almost consistent with the results predicted by our bioinformatics analysis. As for BxPC-3 cell line, the expression levels of BIK, COQ4, MRPL14, MRPL50, and NMNAT3 were significantly up-regulated. As for CF-PAC1 cell line, the expression level of BIK was significantly up-regulated and the expression level of MRPL50 was significantly down-regulated. As for Panc-1 cell line, there is no significant difference as regards above eleven model genes. As for Mia-Paca-2, the expression levels of BNIP3L and NMNAT3 were significantly up-regulated. Using the HPA database, we searched through immunohistochemical findings of tumor and normal pancreatic tissues to verify the expression pattern of 11 hub MTGs-OS at the protein levels ( Figure 10 ). At the same time, we studied the cellular localization of the 11 MTGs-OS, however, we did not find the information of cellular localization of BIK and MCRIP2 in the HPA database. Therefore, we showed cellular localization of the other 9 MTGs-OS (COQ4, SOD1, NDUFB8, MRPL50, MRPL14, RFK, NMNAT3, BNIP3L, MRPS5) in Figure 10 . The expression products of COQ4, MRPS5, MRPL14, NDUFB8, and NMNAT3 were mainly located in the mitochondria. The expression products of RFK, MRPL50, and SOD1 were located Golgi apparatus, mitochondria and cytosol, nucleoplasm and plasma membrane and cytosol, and mitochondria and nuclear speckles, respectively.

Based on the result above, the evaluation of the prognostic performance of the MTGs-RPS was displayed in Figure 11 . PC patients in the train set were classified into two subgroups by KM survival analysis. The high-risk category recorded a lower overall survival time than the low-risk category, which suggested that the MTGs-RPS could precisely differentiate patients with a poor and a favorable prognosis for PC ( Figure 11A ). As per ROC analysis, we proved the accuracy and robustness of MTGs-RPS as a diagnostic tool. The 5-year survival AUC of the ROC curve AUC was 0.780 ( Figure 11E ). Furthermore, both internal and external validations of survival data found that PC patients with low-risk scores recorded superior overall survival rates ( Figures 11B–D ). In addition, when comparing the test1 with the test2 and external validation sets, we discovered that the AUC values of the 5-year ROC curve were 0.658, 0.742, and 0.673, respectively ( Figures 11F–H ). The results above could suggest the dependability and accuracy of the MTGs-RPS.

During the process of establishing the MTGs-RPS, we classified 447 PC patients into low- and high-risk categories as per the median risk in the train set ( Supplementary Figure 2A ). As per the distributions of risk scores, the mortality rate was remarkably greater among the high-risk subset of PC patients ( Supplementary Figure 2E ). The prognostic signature encompassed 11 different MTGs-OS, and their expression patterns were displayed in a heatmap ( Supplementary Figure 2I ). As we had anticipated, the same technique effectively classified PC patients into low- and high-risk categories across all three cohorts (test1, test 2, and test 3) ( Supplementary Figures 2B–D ). Notably, the median risk score functioned as a benchmark in the train dataset when separating the samples in test1, test2, and test3. In the three cohorts, similarities were observed between the risk score distributions and survival status recorded in the train set and those reported in the internal (test1 and test2 cohorts) and the external validation dataset ( Supplementary Figures 2F–H ). The expression patterns of the MTGs-OS reported by heatmaps in the test1, test2, and test3 cohorts were similar to that in the train cohort ( Supplementary Figures 2G–L ).

As mentioned above, targeted therapies have shown promise in the management of PC, and this is becoming increasingly important. Hence, to identify potentially sensitive medications in both high- and low-risk patients, the “Wilcox.test” tool was used. To be selective, a drug had to show statistical significance across all four categories (train, test1, test2, and test3); we found that 17 distinct drugs had different sensitivities for all the cohorts ( Figures 12A–D ). The 17 drugs were listed as follows: AUY922, CGP.60474, cisplatin, CMK, dasatinib, docetaxel, EHT.1864, FTI.277, imatinib, JNJ.26854165, LFM.A13, metformin, midostaurin, pazopanib PLX4720, salubrinal, and thapsigargin.

The different levels of immune cell infiltration in the low- and high-risk patients were evaluated according to the TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC algorithms. For both derivation and validation cohorts, the distribution of immunocyte infiltration exhibited obvious discrepancy between high- and low-risk category, further implying the crucial role of mitochondria genes associated with oxidative stress in tumor immune microenvironment ( Figures 13A–D ). In addition, we also evaluated the levels of expression of ICGs in high- and low-risk patients, given the significance of ICGs in the anti-tumor activity of immune cells. All of the ICGs in the train set exhibited substantial differences, as our findings showed ( Figure 13E ). In the test1 datasets, there was a substantial difference in four ICGs between low- and high-risk categories ( Figure 13F ). The findings of the immune checkpoints performed on the test2 cohorts were comparable to those performed on the train cohorts ( Figure 13G ). There were nine ICGs that substantially expressed differing amounts in the low- and high-risk patients in the test3 cohorts ( Figure 13H ). Even though the number of genes with differential expression varied across the training and validation sets, the expression pattern of nearly all ICGs was reduced in the low-risk category in contrast to the high-risk category. These findings might explain why the low risk patients had a favorable prognosis.

After performing cluster analysis, 226 patients with PNET were successfully stratified into three clusters (i.e. C1, C2, and C3) ( Supplementary Figure 3A ). Both C1 and C2 clusters demonstrated relatively higher MTGs scores than C3 cluster; however, there is no significant difference of MTGs scores between C1 and C2 clusters ( Supplementary Figure 3B ). Thus, C1 and C2 clusters with similar MTGs scores were subsequently integrated and defined as S1 subtype. The C3 cluster with relatively lower MTGs scores was regarded as S2 subtype. Likewise, significantly higher MTG scores were observed in PNET patients with S1 subtype than S2 subtype ( Supplementary Figure 3C ). In addition, we compared the expression profiles of oncogenes and tumor suppressor genes between S1 and S2 subtypes ( Supplementary Figure 4 ). Of note, the expression levels of GNAS, HRAS, PDZRN3, and SYNE2 showed significantly different between two subtypes. More importantly, we also compared the discrepancies in metabolism and immune-related pathways between two subtypes. Our findings showed that the S1 subtype was accompanied by higher activities of sphingolipid metabolism and lower activities of B cell and T cell receptor signaling pathway ( Supplementary Figure 5 ).