We downloaded the normalised gene expression data of cancer and normal tissues from the Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression Project (GTEx) databases in UCSC Xena (http://xena.ucsc.edu/). Data was graphically displayed using the ggplot2 R package. For differences in protein expression between cancer and normal tissues, we used the clinical proteomic tumor analysis consortium (CPTAC) analysis tool of the UALCAN database (http://ualcan.path.uab.edu/index.html), and no LUAD protein expression data were found for PARK6 (PINK1), PARK10 (ELAVL4), and PARK16 (SLC41A1). The expression of Parkinson family genes at different stages of LUAD was explored using the GEPIA2 database (http://gepia2.cancer-pku.cn/#index).

We obtained and downloaded the gene expression and detailed pathological data of 436 LUAD patients and 425 LUSC patients (primary solid tumour samples with detailed prognostic information and a survival time of up to 15 years were included) from the SangerBox database (http://www.sangerbox.com/) to analyse the impact of Parkinson family genes on cancer prognosis. Data was graphically displayed with the help of the survival R package. The forest plot was constructed using GraphPad Prism 8.

Based on the expression and prognostic value of 25 Parkinson family genes in 436 LUAD and 425 LUSC patient samples, we used the survival R package to construct a Least absolute shrinkage and selection operator (LASSO) regression model. Risk factor analysis was performed using the Hiplot online analysis platform (https://hiplot.com.cn/). A prognostic nomogram model for the 436 LUAD patients was constructed using the rms R package, and the bootstrap method was applied to assess consistency. To validate the model in external data sets, we selected three Gene Expression Omnibus (GEO) chips (GSE37745, GSE31210, and GSE30219) of LUAD patients with detailed prognostic information. The receiver operating characteristic (ROC) curve and time-dependent ROC curve of the working characteristics of subjects were established using R in order to evaluate the survival prediction accuracy of the seven-gene signature for the three chips.

Univariate and multivariate analyses were conducted using the IBM SPSS Statistics 19. The tumour mutation burden (TMB) and neoantigen (NEO) data of 436 LUAD patients were obtained from the Cancer Imaging Archive (TCIA) database (https://tcia.at/home) (Van Allen et al., 2015; Hugo et al., 2016), and both data were available only for 363 patients. We explored mutations of the seven genes screened via LASSO regression (TCGA, PanCancer Atlas) using the cBioportal database (https://www.cbioportal.org/). Employing the HitPredict database (http://www.hitpredict.org/) (Patil and Nakamura, 2005; Patil et al., 2011; López et al., 2015), we searched for interaction partners of the seven genes. The protein network interaction map was constructed using Cytoscape_v3.7.1. GO and KEGG pathway enrichment analysis of the seven genes and Gene Set Enrichment Analysis (GSEA) analysis of high- and low-risk LUAD patients were carried out using the OmicShare 6.2.1 online tool.

In the TIMER database (http://cistrome.dfci.harvard.edu/TIMER/), we explored the impact of mutations in the seven genes on the degree of infiltration for six immune cell types within the tumour microenvironment (TME). Based on the ESTIMATE database (https://bioinformatics.mdanderson.org/estimate/), we further analysed the immune scores of 436 LUAD patients in different risk groups. Next, we used seven analysis methods (TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC) to determine the immune cell infiltration status within the TME of patients in the TIMER2 database (http://timer.cistrome.org/). In addition, we explored differences between the high- and low-risk groups in the steps of the cancer immunity cycle for these patients using the TIP database (http://biocc.hrbmu.edu.cn/TIP/) (Xu et al., 2018).

Gene expression and detailed pathological data of anti-PD-1-treated LUAD patients were obtained from the GSE135222 dataset. The correlation between risk scores and various immunosuppressive molecules was graphically displayed utilising the online analysis tool Hiplot. R packages were used to analyse the correlation between risk score and multiple treatment targets.

In order to explore the expression of Parkinson family genes in LUAD, we obtained the normalised cancer and normal tissue gene expression data of LUAD patient samples from the TCGA and GTEx databases in UCSC Xena. All Parkinson family genes exhibited significant differential expression, with the exception of PARK14 (PLA2G6) and PARK20 (SYNJ1) (Figure 1A). Similar results were obtained for LUSC (Supplementary Figure S1). Furthermore, using the UALCAN database, we analysed the protein expression of Parkinson family genes between LUAD and normal tissue. Seventeen genes exhibited significant differences in expression at the protein level (Figure 1B). Finally, we explored the expression of Parkinson family genes at different TNM stages in the GEPIA2 database and found that four genes had significant expression differences during the progression of LUAD (Figure 1C). In summary, all Parkinson family genes exhibited significant differential expression in LUAD, suggestive of their involvement in the occurrence and development of lung cancer.

Next, we selected 436 LUAD patients from TCGA to investigate the impact of Parkinson family genes on the prognosis of LUAD. We found that 12 genes had significant differences based on LUAD prognosis (Figures 2A,B). Similarly, 11 genes had a significant impact on the prognosis of LUSC patients (Supplementary Figures S2A,B). Moreover, we performed GSEA analysis to explore the influence of deregulation of these genes on 50 hallmark gene sets. Results showed increased expression for these genes associated with poor prognosis can obviously activate a variety of cancer-related pathways (WNT, P53 and NOTCH pathway for LUAD; MYC, G2M checkpoint and Oxidative phosphorylation for LUSC). Interesting, whether in LUAD or LUSC, the impact of the increase of protective prognostic genes and risk genes on the 50 hallmark gene sets were obviously different and even opposed to each other, this finding may provide some help for follow-up in-depth study (Supplementary Figure S3, S4).

However, not all differentially expressed Parkinson genes (in Figure 1A) accurately predicted overall survival in LUAD patients. Therefore, we sought to create a gene signature from 25 Parkinson family genes that had the most pronounced impact on LUAD prognosis. To this end, we constructed a LASSO regression model based on the expression and prognosis data of 436 LUAD patients. We obtained a seven-gene prognostic signature (Figures 2C,D). The complete names of the seven genes and their main functions are listed in Table 1. Meanwhile, we further analyzed the correlation between the 7 genes and genes in the non-small cell lung cancer pathway (map05223) from KEGG database. The results showed that the 7 genes were related to multiple cancer progression marker genes respectively, this result suggest these genes were involved in the development of lung adenocarcinoma via different mechanisms (Supplementary Figure S5). We also constructed a LASSO regression model based on the data of LUSC patients. However, limited by data quality, the results indicated that a combination of two genes had the best performance (Supplementary Figures S2C,D). Therefore, we focused on the seven-gene signature for the prognosis and treatment outcome prediction in LUAD.

To better assess the prognostic value of our seven-gene signature in LUAD, we first conducted risk factor analysis based on the genes. The occurrence of death was significantly correlated with a higher risk factor score (Figure 3A). Furthermore, we divided 436 patients into high- and low-risk groups based on risk factor score and analysed prognosis in the two groups. The high-risk group exhibited a significantly poorer prognosis among patients with different TNM stages (Figure 3B). Univariate and multivariate analyses also indicated that the signature-based risk factor score was a superior prognostic indicator compared to TNM staging (Figures 3C,D). Taken together, the seven-gene signature of Parkinson family genes can be used as an independent prognostic marker for LUAD.

In order to validate and further assess the prognostic value of our seven-gene signature for LUAD, we constructed a nomogram based on data from the TCGA database in order to predict patient survival probability by weighing age, gender, stage, T, N, M, and the signature-based risk score (Figure 4A). Further, we applied the bootstrap method to evaluate the nomogram’s predictive performance. The calibration curves indicated that the nomogram-predicted probability matched the actual 3- and 5-years survival (Figures 4B,C). Subsequently, we selected three GEO datasets (GSE37745, GSE31210, and GSE30219) with detailed prognostic information for external testing (Table 2). The seven-gene signature risk scores suggested a significantly better diagnostic performance than TNM staging based on the ROC curve, also showing excellent prognostic performance at different survival time points in the time-dependent ROC curve (Figures 4D–F). Taken together, risk evaluation based on the seven-gene signature has important clinical significance in the diagnosis and treatment of LUAD.

We explored the possible mechanisms underlying the functional association of signature genes with poor prognosis. Utilising the HitPredict database, we searched for interaction partners of their protein products, which are displayed in the protein interaction network diagram (Figure 5A). We then performed GO and KEGG enrichment analyses for these genes. GO enrichment analysis indicated the genes’ involvement in diverse tumourigenesis- and development-associated signalling, such as metabolism-related pathways, ubiquitination pathways, and RNA transcription-related processes (Figures 5B–D). Furthermore, KEGG pathway enrichment analysis also revealed that immune and cell proliferation-related pathways (RIG-I-like receptor, MAPK, mTOR, and ErBb) were enriched by these genes (Figure 5E).

In view of the fact that various immune-related processes were also observed among functional enrichment results (Figure 5E), we sought to explore the influence of the seven signature genes on the TME. Employing the TIMER database, we obtained the effect of various single gene mutations on the infiltration of six immune cell types (B cells, CD8+T cells, CD4+T cells, macrophages, neutrophils, and dendritic cells) into the TME. The results indicated that mutations of signature genes decreased the degree of infiltration for almost all six immune cell types (Figures 6A–G). We also explored the link between risk factor score and TMB (or NEO) based on TCIA data. A high-risk score was clearly related to a high TMB and high NEO (Figure 6H). Mutation data from the cBioportal database indicated that mutations in these seven genes were rare among TCGA-LUAD patients (Figure 6I).

We further explored the effect of the seven-gene signature on the TME. Estimation analysis indicated that the immune score of the high-risk group was significantly lower than that of the low-risk group for the 436 LUAD patients (Figure 7A). We further observed the different infiltration levels of 22 immune cell types between the low- and high-risk groups using CIBERSORT (Figure 7B). Meanwhile, in order to eliminate possible errors caused by different analysis methods, we employed another six methods (Supplementary Figure S6) and summarised the results, revealing the same trend and significant differences (Table 3). However, these results revealed an interesting phenomenon. In the high-risk group, the infiltration of DC cells, B cells, CD4+ T cells, and CD8+ T cells remained lower, while M1 pro-inflammatory macrophages were upregulated, and M2 anti-inflammatory macrophages and Treg cells were downregulated.

The cancer immunity cycle reflects a complex interaction network involving the chemokine system, other immunomodulators, and the various immune cell types. Based on TIP data, we found that high risk was associated with the release of cancer cell antigens (Step 1), reduced priming and activation (Step 3), recruitment of multiple immune cells (Step 4), and the infiltration of immune cells into tumours (Step 5) (Figure 7C). It is worth noting that Th1 cells were significantly downregulated in the high-risk group, while Th2 cells were upregulated, although the difference was not statistically significant. This partly explains why the various effector cells of cellular immunity exhibited low infiltration with an opposite trend of change compared to regulatory cells. The detailed underlying mechanism remains to be further explored. Similarly, GSEA analysis showed that signature-based high risk was associated with a significant inhibition of tumour immune-related processes (Figure 7D). Exploring the mechanism underlying the increase in M2 cells and Treg cells in the high-risk group may help explain the low infiltration of effector cells.

Finally, we sought to explore the utility of our seven-gene signature as a predictor of immunotherapy and targeted therapy response in LUAD. A high-risk score was significantly negatively correlated with a variety of immunosuppressive molecules (Figure 8A). Furthermore, we explored the correlation between risk score and the efficacy of PD-L1 immunotherapy in LUAD. Data from GSE135222 indicated that after PD-1 immune checkpoint inhibitor treatment, the risk score of dying patients was higher than that of surviving patients, and there was a negative correlation between the survival time for surviving patients and the risk score (Figure 8B). Limited by the sample size, the results were not statistically significant, and we did not find other available immunotherapy cohorts with detailed prognostic and gene expression data. Nevertheless, we believe that these results would be useful for treatment selection. We analysed the correlation between risk score and the expression of multiple driver genes for targeted therapy in TCGA patients, observing a clear correlation (Figure 8C). Furthermore, based on the Genomics of Drug Sensitivity in Cancer (GDSC) database, we predicted the impact of 7 genes on the half-inhibitory concentration (IC50) of some common chemotherapeutics (platinum, paclitaxel, etc) and targeted agents for non-small cell lung cancer (NSCLC). Interestingly, the result showed high expression of SLC41A1, RAB29 and PLA2G6 may be positively correlated with resistance of various common chemotherapeutics, whereas high level of CHCHD2, EIF4G1, DNAJC13 and SFXN5 may be involved in resistance of target therapy for drive gene (Supplementary Figure S7). However, we did not find available targeted therapy datasets for validation. In general, our findings indicated that targeted therapy may be more appropriate than immunotherapy for high-risk patients as a precision medicine approach for LUAD.