We obtained RNA sequencing (RNA-seq) data and corresponding clinical information from the GEO public database (https://www.ncbi.nlm.nih.gov/geo/) for a total of 695 NSCLC samples. The obtained data underwent standardization using the R package “limma”. Among these samples, a training cohort consisting of 539 merged samples from the GSE30219 (16) and GSE32863 (17) series matrices were selected. To mitigate batch effects within the merged samples, we used the R package “SVA” (18). Additionally, we included 156 samples from the GSE19188 (19) series matrix for external validation. Single-cell data from 11 tumor samples were also included from the GSE131907 (20). Furthermore, we retrieved the mitophagy-related gene set from the Reactome database and the inflammatory factor-related gene set from the BioCarta database.

Data from the GSE30219 and GSE32863 series matrices were merged and classified into normal and tumor samples based on clinical information. We extracted MRGs from the Reactome database and generated a heat map using the R package “pheatmap”. DEGs were detected between the normal and tumor samples, with a false discovery rate (FDR) threshold of less than 0.05. Subsequently, DEGs boxplots and volcano plots were constructed using the R packages “ggplot2” and “ggpubr”. In addition, a correlation heatmap of DEGs was visualized using the R package “corrplot”.

We performed analyses of scRNA-seq data using the R packages “Seurat” and “SingleR” (21, 22). To retain high-quality scRNA-seq data, we applied a series of filtering steps on the raw matrix. Specifically, we included only genes expressed in at least three cells, excluded cells expressing fewer than 200 genes, and removed cells with mitochondrial gene expression exceeding 8%. First, we applied the NormalizeData function to normalize the scRNA-seq data and then used the FindVariableFeatures function to identify the top 2,000 highly variable genes. Subsequently, we performed principal component analysis (PCA) using the RunPCA function from the “Seurat” package, reducing the dimensionality of the scRNA-seq data based on the top 2,000 highly variable genes. We used the ElbowPlot analysis to identify significant principal components and performed cell clustering analysis using the top 10 principal components obtained from PCA. The FindNeighbors function was employed to construct a k-nearest neighbor graph based on the Euclidean distance in the PCA space, and the UMAP algorithm was applied for data dimensionality reduction and visualization. We used reference data from the human primary cell atlas (23) for cluster annotation and the CellMarker database (CellMarker (xbio.top)) (24) for reference-based cell annotation.

To evaluate the gene expression scores of MRGs for each cell in the scRNA-seq dataset, we employed the AddModuleScore function to compute specific scores for individual cells. Furthermore, we performed GSVA to assess the enrichment of signature pathways for each cell type using the MSigDB database (http://www.gsea-msigdb.org/) (25). Following this, we performed GSEA on the differentially expressed genes between the high and low groups based on mitophagy scores (26).

Machine learning models have proven to be powerful tools in handling large datasets and facilitating rapid analysis and decision-making in disease diagnosis. In our study, we employed a combined approach using support vector machine (SVM) (27) and random forest (RF) (28) algorithms to select relevant features from a set of 25 genes associated with mitophagy. Subsequently, a prediction model was constructed using logistic regression (29). The sensitivity and specificity of the prediction model were assessed through the analysis of receiver operating characteristic (ROC) curves. The performance of the model was evaluated by calculating the area under the curve (AUC) of the ROC curve. To visualize the relative importance of the feature genes within the model and their contribution to the prediction outcomes, a nomogram model was generated using the “rms” package in R (30). The prediction performance of the model was further evaluated by plotting calibration curves and decision Curve Analysis (DCA) curves. An external dataset was employed to validate the robustness of the results (31).

To validate the prediction model, the HPA database (proteinatlas.org) was used to ascertain if the expression levels of the feature genes in NSCLC differ from those in normal tissues at the protein level.

The bronchial epithelial cell BEAS‐2B and the NSCLC cell lines A549, PC9, and H1299 were from the Cell Bank of Type Culture Collection of Shanghai Institute of Biochemistry & Cell Biology, Chinese Academy of Science. BEAS‐2B cells were cultured in Gibco LHC basal medium (Invitrogen, Carlsbad, CA) supplemented with 10% fetal bovine serum, while the NSCLC cell lines were cultured in Gibco RPMI‐1640 medium (Invitrogen, Carlsbad, CA) supplemented with 10% fetal bovine serum. Total RNA was extracted using TRIZOL reagent (Invitrogen) and then reverse-transcribed into cDNA. β‐actin was used as an internal control. Real‐time PCR was performed using a Roche LightCycler® 480 System. Relative expression levels were calculated as ratios normalized against β‐actin. The QPCR primer sequences are listed in Supplementary Table S1.

To explore the correlation between feature genes and immunity, we employed the “CIBERSORT” package to evaluate the extent of immune cell infiltration in both normal and tumor samples. Furthermore, we performed Spearman’s correlation analysis using the “IOBR” and “psych” packages to investigate the relationships among feature genes, immune infiltration, and inflammatory factors. The results were visually represented using the “ggplot2” package (32, 33).

We performed unsupervised clustering analysis on 351 samples from patients with NSCLC using the “ConsensusClusterPlus” R package based on the expression profiles of six feature genes. The K-means algorithm was employed with a maximum number of subtypes set to 10 (k=10), a sampling rate of 0.8, and 100 resampling iterations. The optimal number of clusters was determined by evaluating the cumulative distribution function (CDF) curve and CDF delta area curves (34).

Following the clustering analysis, we investigated the correlation between the expression levels of MRGs and clinical pathological features within different subtypes of NSCLC. In addition, using the gene set variation analysis (GSVA) package, we calculated the scores of each sample relative to the MRGs. The significance of the differential expression of MRGs between the two groups of NSCLC subtypes was assessed using the Wilcoxon test (35). Furthermore, Kaplan-Meier survival analysis was performed on the different subtypes using the “survival” and “survminer” R packages to evaluate their impact on prognosis.

In this study, we employed the CIBERSORT software code in R to evaluate the proportions of 22 immune cell subtypes within the different subtypes of NSCLC. The LM22 gene expression matrix, which includes 547 genes and defines the composition of 22 immune cell types, served as a standardized reference for characterizing the human immune cell types. By applying the CIBERSORT algorithm, we obtained the profiles of the immune cell infiltrations for each NSCLC subtype. Additionally, we used the “IOBR” package in R to estimate the expression levels of the immune and stromal cells in the samples. A comparison of the abundance of immune cells and stromal cells among different subtypes was performed, and the results were visually presented using box plots (32).

We performed consensus module analysis using the “WGCNA” package in R to identify hub genes within each subtype (36). Initially, we ranked genes based on their median absolute deviation (MAD) values and selected the top 5000 genes. Next, we performed sample and gene filtering to exclude those with excessive missing values. A similarity matrix was constructed, which was then converted into an adjacency matrix using a weight coefficient β = 6. Subsequently, the blockwiseModules function was employed to perform modular analysis on the adjacency matrix, enabling the establishment of a gene co-expression network and the assignment of genes to different modules. For each module, module eigengenes (MEs) were computed, and their Pearson’s correlation coefficients with the phenotype data were calculated. Finally, the correlation between the module genes and phenotype data was visualized using a heat map.

Based on the WGCNA analysis results, we selected genes from the module with the highest correlation to different NSCLC subtypes. Subsequently, we performed enrichment analyses using the Metascape database (https://metascape.org/) (37). In addition, we imported gene sets from the biological processes (BP) category using the “msigdbr” package in R (38) and performed biological process enrichment analysis for different NSCLC subtypes using the “fgsea” package.

First, we analyzed immune checkpoint gene expression in various subtypes of NSCLC. Subsequently, we employed the Tumor Immune Dysfunction and Exclusion (TIDE) database (Tumor Immune Dysfunction and Exclusion (TIDE) (harvard.edu)) to predict the response to immune checkpoint blockade. TIDE can analyze two mechanisms involved in tumor immune evasion: one causing dysfunction of cytotoxic T lymphocyte (CTL) infiltration and the other impeding CTL infiltration into the tumor, known as “exclusion” (39, 40). The use of TIDE facilitated a more effective prediction of the therapeutic efficacy of immunotherapy.

The half-maximal inhibitory concentration (IC50) is a commonly used concentration indicator to assess drug toxicity or therapeutic efficacy. We used the “oncoPredict” package in R and employed expression matrices from the Genomics of Drug Sensitivity in Cancer (GDSC) database, along with drug treatment information, as a training set to evaluate chemotherapy sensitivity in different subtypes of NSCLC (41).

All the statistical analyses were performed by R-4.1.3. Student's t-test or Wilcoxon’s rank sum test was used to detect the significant difference between two independent groups, p < 0.05 was considered statistically significant.

We performed batch correction and merged the GSE30219 and GSE32863 datasets to investigate the expression profiles of MRGs in the combined dataset. Our analysis revealed that the expression levels of MRGs were significantly lower in normal tissues than in tumor tissues ( Figure 1A ). Among the 25 MRGs exhibiting differential expression between the normal and tumor tissues, only MAP1LC3B, PINK1, PARK2, UBB, and UBC were significantly downregulated in the tumor tissues, while expression of the other genes was significantly upregulated in the tumor samples ( Figures 1B–D ). Furthermore, correlation analysis demonstrated the associations among the 25 differentially expressed genes related to mitophagy ( Figure 1E ). Using CIBERSORT analysis, the results showed a significant increase in the infiltration of naive B cells, M0 and M1 macrophages, NK cells, plasma cells, helper T cells, and regulatory T cells in the tumor tissues ( Figure 1F ).

Based on the scRNA-seq data from GSE131907, we performed quality control on 11 tumor samples in the single-cell dataset ( Supplementary Figure 1A ) to ensure the quality of the cells included in the study. Then, the top 2,000 highly variable genes were identified for PCA to reduce dimensionality ( Supplementary Figures 1B, C ). Subsequently, we clustered the cells into 21 clusters ( Supplementary Figure 1D ). We classified them roughly into T cells, NK cells, macrophages, monocytes, mast cells, epithelial cells, endothelial cells, B cells, and fibroblasts based on the results of singleR and CellMarkers database ( Supplementary Figures 1E, F , Figure 2A ). The AddModuleScore function was used to calculate the score of the different cell types based on the MRGs. We found that fibroblasts had the highest MRG score (mitoscore), followed by mast cells ( Figures 2B–D ). We then performed GSVA analysis on different cells and found that the HALLMARK inflammation-related pathways such as apoptosis, inflammatory_response, TGF_BETA_SIGNALING, and TNFA_SIGNALING_VIA_NFKB were upregulated in mast cells ( Figure 2C ). In addition, a GSEA analysis was performed between samples with an Hscore and Lscore for mitophagy-related scores. We found that the differentially expressed genes between the high and low groups were highly enriched in pathways such as KEGG_neuroactive_ligand_receptor_interaction, KEGG_cytokine_cytokine_receptor_interaction, KEGG_natural_killer_cell_mediated_cytotoxicity, and KEGG_ecm_receptor_interaction ( Figure 2E ).

We used SVM and RF algorithms to select feature genes from the 25 MRGs and predict the occurrence of NSCLC. The SVM analysis allowed us to identify the top 10 optimal feature genes associated with NSCLC occurrence. Subsequently, we performed a random forest analysis, which involved selecting the top 10 genes based on the error rate curve. The gene importance scores were obtained from the random forest model ( Figures 3A, B ). The genes selected by the SVM-RFE algorithm and random forest were intersected to obtain the final set of genes (SRC, UBB, PINK1, FUNDC1, MAP1LC3B, CSNK2A1). Multivariable logistic regression analysis demonstrated that these six selected genes were associated with NSCLC ( Figure 3C ). We further evaluated the predictive performance of the model using an ROC curve, and the area under the curve (AUC) was determined to be 0.925. In addition, we validated the model using an external dataset (GSE19188), and the AUC based on this validation was 0.966 ( Figures 3D, E ). To visually illustrate the prediction performance of the feature genes, we created a nomogram model ( Figure 3F ). The model, which integrated the information from all six feature genes, showed superior prediction value compared to individual feature genes. To assess the calibration and clinical utility of the model, we generated calibration plots and DCA curve analysis ( Figures 3G, H ). These analyses demonstrated that the nomogram model had good prediction efficacy for NSCLC. Subsequently, we collected IHC staining images of six feature genes-associated proteins from the HPA database, which were obtained from both NSCLC and normal lung tissue. Significantly higher protein expression levels were observed for three of the feature genes (SRC, CSNK2A1, FUNDC1) in NSCLC samples relative to normal samples, thereby supporting our findings ( Supplementary Figure 2A ). The protein expression levels of the additional three feature genes showed no significant difference between the NSCLC samples and the normal samples. Furthermore, qPCR was used to validate the expression of six feature genes that were screened in both the lung cancer cell lines (A549, PC9, H1299) and the normal bronchial epithelial cell line (BEAS-2B). Results indicated a significant upregulation of the six feature genes in lung cancer cells ( Supplementary Figure 2B ). The expression patterns of CSNK2A1, FUNDC1, and SRC were consistent with our bioinformatics analysis results, therefore, we attempted to construct a prediction model using these three genes. The predictive performance of the model based on the three feature genes was not as good as that of the six-gene model, as evidenced by the ROC curve ( Supplementary Figure 2C ), and the AUC based on this model was only 0.882. In addition, we analyzed the correlation between the six feature genes and immune cells, as well as inflammatory factors. Among the selected feature genes, UBB, PINK1, and MAP1LC3B exhibited a negative correlation with the infiltration of the aforementioned immune cells, while the other feature genes showed a positive correlation ( Supplementary Figure 3A ). UBB, PINK1, and MAP1LC3B were significantly positively correlated with 10 out of 24 inflammatory factors, whereas the remaining feature genes showed the opposite trend ( Supplementary Figure 3B ).

Based on the six feature genes related to mitochondrial autophagy, unsupervised clustering was performed to classify all tumor samples into subtypes. Through comprehensive analysis of the cumulative distribution function (CDF) curve ( Figure 4B ) and delta area ( Figure 4C ), the optimal segmentation efficiency was achieved at k=2, with a clustering number of 2 ( Figure 4A ). The heat map illustrates the differences in expression levels of the feature genes and the distribution of clinical pathological features between the two subtypes ( Figure 4D ). We observed that the prediction model demonstrated good efficacy for both subtypes of NSCLC ( Figure 4E ). Cluster 1 had a significantly higher number of individuals in tumor stages T1–T2 and N0–N1 than Cluster 2 ( Figure 4F ). Survival analysis also demonstrated a significantly higher survival rate in Cluster 1 than in Cluster 2 ( Figure 4G ). GSVA analysis revealed a significantly higher overall expression level of MRGs in Cluster 2 than in Cluster 1 ( Figures 4H, I ).

A weighted gene co-expression network analysis (WGCNA) was performed to identify co-expression modules that showed the highest correlation with each cluster. WGCNA was performed on both datasets of diseased samples after batch normalization. All samples were included in the analysis, and the optimal soft-thresholding value was set to 6. The results revealed a total of nine modules generated ( Figures 5A, B ). Among these modules, the black module and the turquoise module exhibited the strongest correlation with Cluster 1 and Cluster 2, respectively ( Figure 5C ). Subsequently, we performed functional enrichment analysis on the WGCNA modules that exhibited the highest correlation with each cluster. For this analysis, we used the Metascape database to analyze the genes within the black module and the turquoise module. The results revealed significant enrichment in the pathways associated with angiogenesis, endothelial development, and cell adhesion in Cluster 1 ( Figure 6A ). In Cluster 2, there was significant enrichment in the pathways related to cell cycle and metabolism ( Figure 6B ). Furthermore, FGSEA highlighted that Cluster 1 was enriched in the pathway of antigen processing and presentation of exogenous antigens. On the other hand, Cluster 2 exhibited enrichment in pathways such as signal transduction by p53 class mediator and mitochondrial translation ( Figure 6C ).

To investigate the differences in the immune microenvironment among different subtypes of NSCLC, we used the CIBERSORT package and IOBR package to assess the infiltration proportions of distinct immune and stromal cells. The infiltration levels of dendritic cells, M1 macrophages, mast cells, neutrophils, fibroblasts, and T cell CD4 memory cells were significantly elevated in Cluster 2, while monocytes, NK cells, and endothelial cells exhibited higher infiltration levels in Cluster 1 ( Figure 7A , Supplementary Figure 4A ). The expression of most inflammatory factors was higher in Cluster 1 than in Cluster 2 ( Supplementary Figure 4B ).

Immunotherapy checkpoint inhibitors have demonstrated durable efficacy in some NSCLC patients. In this study, we investigated the differences in expression levels of immunotherapy checkpoint genes between two NSCLC subtypes. Our analysis revealed significant differences in the expression levels of several immunotherapy checkpoint genes, including BTN2A1, TNFSF9, CD226, BTN3A1, CD47, CD28, and TNFRSF18 ( Figure 7B ). This suggested that these two NSCLC subtypes might respond differently to immunotherapy checkpoint inhibitors. We further employed the TIDE algorithm to predict the effectiveness of immunotherapy checkpoint inhibition treatment. We observed that Cluster 1 exhibited a significantly lower TIDE score ( Figure 7C ), indicating a better response to immunotherapy. In addition, we found that Cluster 1 had a lower Exclusion score and a higher Dysfunction score compared to other clusters. However, there was no significant difference in the microsatellite instability (MSI) score between the clusters ( Figures 7D–F ). To explore the sensitivity to chemotherapy drugs in different NSCLC subtypes, we used the oncoPredict package to evaluate the IC50 values of several chemotherapy drugs. The results of the drug sensitivities showed that the IC50 values of cisplatin, paclitaxel, vincristine, and gemcitabine were lower in Cluster 2 than in Cluster 1. This meant that patients in Cluster 2 might benefit more from the above chemotherapy drugs ( Figures 7G–J ). Additionally, our correlation analysis between drug response and gene expression revealed that the highly expressed MRGs in Cluster 2 negatively correlate with the IC50 values of the aforementioned chemotherapy drugs ( Supplementary Figures 5A-D ).