Three PD datasets, GSE7621, GSE20163, and GSE49036, were downloaded from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). The GSE7621 and GSE49036 datasets were generated using the GPL570 (HG-U133 Plus 2) Affymetrix Human Genome U133 Plus 2.0 array, while the GSE2016 dataset was generated using GPL96 [HG-U133A] Affymetrix Human Genome U133A array. GSE7621 dataset consisted of 16 brain nigrostriatal samples from Parkinson’s disease patients and 9 normal nigrostriatal samples from controls. GSE20163 dataset, which served as an external validation cohort, consisted of 8 PD brain substantia nigra samples and 9 control samples. GSE49036 was used as a validation cohort for clinical staging and included brain substantia nigra samples from 8 Braak stage 0, 5 Braak stages 1-2, 7 Braak stages 3-4 and 8 Braak stages 5-6 patients.

RNAseq data, mutation data, and clinicopathological characteristics of TCGA-GBM, consisting of 169 glioma samples, were retrieved from the UCSC Xena website (https://xena.ucsc.edu/) Gene expression data for 249 glioma patients were retrieved from the China Glioma Genome Atlas (CGGA) data portal (http://www.cgga.org.cn) and were used to generate a validation model. All expression data were retrieved in TPM format. Batch correction and integration of the two sets of gene expression data were carried out using the “limma” and “sva” (33) packages in R. The detailed flow chart is shown in Figure 1 .

2601 neurotrophic factor-related genes were downloaded from the GeneCard database (https://www.genecards.org/) (34). Differential gene expression analysis was performed on the TCGA cohort using the “limma” package in R, with | log2FC | > 1.0 and FDR (false discovery rate) < 0.05 as the thresholds. The cutoff p-value of the differentially expressed NFRGs (DENFRGs) for the GEO cohort was set to 0.05, which satisfied the condition of |log2FC|>0.5 The “affy” package in R was used to perform background calibration, normalization, and log2 conversion on all GEO raw data sets (35, 36). The expression values of multiple probes that matched the same gene were averaged. Protein interactions and gene enrichment analysis was carried out using the differentially expressed genes identified from the GEO cohort. The hub genes in the network were screened and visualized using “Cytoscape” software following PPI (Protein-Protein Interaction Networks) analysis on the String online platform.

To identify signature genes, we used two machine learning methods: LASSO regression analysis and random forest. LASSO is utilized as a dimensionality reduction approach to perform variable screening and complexity adjustment when fitting a generalized linear model. The LASSO analysis was carried out with a penalty parameter and a 10-fold cross-verification using the glmnet program (37). RF (Random Forest) is a combination of classifiers with a tree-like structure, and a minimum error regression tree was built to select key variables using the software package “randomForest”. After tenfold cross-validation, the eight genes with the highest relative importance were used to measure predictive performance.

The “clusterProfiler” package of R software was used to carry out the functional enrichment analysis, which included KEGG and GO analysis (38). We adjusted the P values using the Benjamini - Hochberg (BH) technique. A computational technique was used in the gene set enrichment analysis to identify genes that exhibited statistically significant and consistent changes between two biological states. 10,000 permutation tests were used to determine the most important and pertinent signaling pathways. Genes with a corrected P-value and false discovery rate (FDR) below 0.05 were considered to be significant. Statistical analysis and ridge mapping were carried out using the “clusterPro” package in R, which is a non-parametric unsupervised analytic method that is widely employed to evaluate gene set enrichment outcomes in microarrays and transcriptomes. It is primarily used to determine if certain metabolic pathways are enriched across samples by transforming the expression matrix of genes across samples into the expression matrix of gene sets (39). From MSigDB, 50 reference gene sets for hallmark genes were chosen. Gene set variation analysis (GSVA) using the ‘GSVA’ package in R was performed to provide insight into the heterogeneity of biological processes between different clusters.

Batch effects between TCGA and CGGA data were removed by creating precise models using the “sva” package in R. Selected NFRGs underwent Minimum Absolute Shrinkage and Selection Operator (LASSO) regression analysis, with the “glmnet” package in R being used to minimize the number of genes in the final risk model. Models were then built using multivariate Cox regression analysis using the following equation: risk score = (Expi), where Expi was the expression value for each NFRG and was the matching regression coefficient (40, 41). The median risk score was used to split all patients into high- and low-risk groups. The “survminer” and “ggrisk” packages in R were used to create survival curves and risk maps to display the disparities in survival and status of each patient. A separate external cohort, the CGGA cohort, was also employed to evaluate the effectiveness of the prognostic model.

A nomogram was created using risk score and clinicopathological features. Calibration charts were internally validated to ensure accuracy of the models. Decision curve analysis (DCA) was carried out using “ggDCA” package in R to evaluate the net clinical benefit of the models (42). We also plotted subject operating characteristic curves using the “timeROC” package in R to evaluate how well risk scores performed in predicting 1-year, 3-year, and 5-year OS in LGG patients (43).

The relative enrichment scores of tumor-infiltrating immune cells (TIICs) were calculated using the R script ssGSEA (single-sample genomic enrichment analysis). We utilized CIBERSORT to calculate and compare the proportion of immune cell types between the low- and high-risk categories, with the sum of all anticipated immune cell type scores in each sample being equal to 1 (44). The TICCs data was downloaded from TIMER 2.0 (http://timer.cistrome.org). The results from TIMER, CIBERSORT, amounts, MCP-counter, xCELL, and EPIC algorithms were also compared between the two groups. The “oncoplot” function in the “maftools” package of the R software was used to create two waterfall plots to compare the specific mutation characteristics between the high- and low-risk groups.

The tumor genomic analysis platform GSCALite (http://bioinfo.life.hust.edu/web/GSCALite/) integrates genomic data for 33 tumor types from the TCGA library, GDSC (Genomics of Drug Sensitivity in Cancer), CTRP (The Cancer Therapeutics Response Portal) medication response data, and normal tissue data from GTEX (Genotype-Tissue Expression) for comprehensive genomic analysis (45).

The Immunological Cell Abundance Identifier (ImmuCellAI) is a computer program launched in 2020 to predict immunological checkpoint reactions based on the abundance of TICCs, particularly certain T cell subpopulations. Comprehensive immunogenomic analysis findings are provided by the Cancer Immunome Atlas (TCIA) online software. Using a scale from 0 to 10, the Immunophenotype Score (IPS) quantifies the immunogenicity of tumors (46). IPS can be used to predict response to immune checkpoint inhibitors. The “prophytic” package in R was used to compute the half-maximal inhibitory concentration (IC50) of samples from the high and low risk score groups in order to test the ability of the risk score to predict sensitivity of samples to chemotherapy and molecular medicines. Zaoqu Liu et al. from the First Affiliated Hospital of Zhengzhou University developed THE BEST website (http://rookieutopia.com/). The database contains sequencing data from a variety of tumors after treatment with immune checkpoint inhibitors

Tumor Immune Single-Cell Hub (TISCH; http://tisch.comp-genomics.org) is an extensive single-cell RNA-seq database dedicated to TME. It enables comprehensive analysis of TME heterogeneity across different datasets and cell types.

Human astrocytes (HA), U87 and A172 glioma cells, obtained from the Center for Experimental Medicine, Southwestern Medical University. All cells were grown in 10% fetal bovine serum-supplemented DMEM. Cells were incubated at 5% CO2 and 37°C. TRIzol reagent was used to isolate RNA, while PrimeScriptTM RT kit was used to perform reverse transcription. Quantitative PCR was carried out using Takara’s SYBR Green PCR Master Mix on the StepOnePlus system. Ploidy changes at the gene level were determined using the 2-ΔΔCT method, with GAPDH as the normalization gene. The primer sequences involved in this study are as follows.

EN1:FORWARD : GAAGAACGAGAAGGAGGACAAGCG,REVERSE: CGTGGTGGTGGAGTGGTTGTAC.

LOXL1:FORWARD : GAAGAACCAGGGCACAGCAGAC,REVERSE ATGTCCGCATTGTAGGTGTCATAGC.

GAPDH:FORWARD : ATGGGGAAGGTGAAGGTCG,REVERSE : GGGGGTCATTGATGGCAACAATA. Each PCR reaction was performed in triplicate.

Transcriptomics and proteomics methods were used to study protein expression at the RNA and protein levels in human tissues and organs, using data found in Human Protein Atlas (HPA, https://www.proteinatlas.org/).

All analyses were conducted using R version 4.1.1, 64-bit6. Prognosis and patient survival in various subgroups were compared using Kaplan-Meier survival analysis and the log-rank test. The nonparametric Wilcoxon rank sum test was used to compare continuous variables between the two groups, while Kruskal-Wallis test was employed for comparisons among more than two groups. Univariate and multivariate Cox regression (R package “survival”) analyses were used to identify clinical traits with prognostic potential in the high- and low-risk groups. Spearman correlation analysis was used to assess correlation coefficients. P < 0.05 was regarded as statistically significant in all statistical investigations. The ROC curves, the nomogram model and the Concordance Index were generated using the “survivalROC”, “rms” and the “pec” (C-index) packages in R, respectively. The changes in gene expression between the two isoforms were determined using principal component analysis (PCA).

The limma package was used for background correction and normalization of expression data from the GSE7621 dataset. Batch effects were removed using the sva package, and the box plot corresponding to the processing results is shown in Figure 2A . The expression of each of the 145 differentially expressed genes (DEGs) in the GEO cohort is shown in Figure 2B . GO analysis revealed that the DEGs were enriched in “positive control of kinase activity” in the biological process (BP) category, “neuronal cell body” in the cellular component (CC) category and the “nuclear glucocorticoid” in the Molecular Function (MF) category ( Figure 2D ). KEGG pathway analysis showed that the DEGs were more closely related to “Neuroactive ligand-receptor interaction” ( Figure 2C ). We also constructed PPI networks to investigate the interaction of the proteins encoded by the DEGs based on betweenness centrality, with the central genes in the network being marked in red. ( Figure 2E ).

We identified 847 DEGs from differential gene expression analysis of NFRGs between tumor and normal tissues of the TCGA-GBM cohort. Out of the total DEGs, 489 genes were down-regulated whereas 358 genes were up-regulated in tumor tissue ( Figure 3A ). Univariate Cox regression analysis identified 104 differentially expressed NFRGs with prognostic potential for GBM in the TCGA cohort ( Figure 3B ). We then determined the NFRGs that overlapped from the univariate cox analysis of the TCGA cohort and the differentially expressed NFRGs obtained from the GEO cohort. A total of 12 NFRGs overlapped between the two cohorts indicating that they were associated with the occurrence of Parkinson’s disease and prognosis of GBM ( Figure 3C ). Figure 3D shows the correlation of the expression of these 12 NFRGs in the TCGA cohort, while Figure 3E shows the localization of these 12 NFRGs on chromosomes, with EN1 being localized on chromosome 2.

Two machine learning algorithms were used to identify key genes among the 12 NRFGs. The best lambda for the LASSO algorithm was 0.138 after ten cross-validations. Due to higher accuracy in comparisons, we used the minimum criterion for the LASSO classifier, and identified 4 key genes, including IRF7, EN1, PLOD3, and LOXL1 ( Figures 4A, B ). The influence of the number of decision trees is shown in Figure 4C . The x-axis shows the number of decision trees, while the y-axis shows the mistake rate. The top 8 key genes with relative relevance scores identified using the random forest technique were PCSK1, S100A4, EN1, CEBPB, IRF7, L1CAM, PLOD3, and LOXL1 ( Figure 4D ). Four key genes, IRF7, EN1, PLOD3, and LOXL1, overlapped from results of lasso regression and random forest algorithm analysis ( Figure 4E ). Figure 4F shows the correlation of these four feature genes in the GEO cohort.

We then estimated the diagnostic performance of the four key genes. The AUC values of the ROC curves were 0.799 for LOXL1 ( Figure 5A ), 0.778 for PLOD3 ( Figure 5B ), 0.861 for IRF7 ( Figure 5C ), and 0.819 for EN1 ( Figure 5D ). GSEA was used to evaluate the signaling pathways associated with the signature genes. Our results showed that LOXL1 ( Figure 5E ), PLOD3 ( Figure 5F ), IRF7 ( Figure 5G ), and EN1 ( Figure 5H ) were mainly associated with functions of the nervous system and the transmission of neurotransmitters. For example, EN1 was strongly correlated with spinal cerebellar ataxia and dopaminergic synapse-related pathways. The GSVA results demonstrated the correlation of the four signature genes with the HALLMARK pathway ( Figure 5I ).

We also quantified the ssGSEA enrichment scores for several immune cell subpopulations, associated PD functions or pathways, and healthy controls. A heat map was used to display the number of TIICs and immunological responses in each sample ( Supplementary Figure 1A ). Supplementary Figures 1B, C show heat maps displaying the relationship between TIICs and immune function, with darker red denoting a stronger correlation between the two. The association between the four key NFRGs and immune-related pathways in the ssGSEA data was also demonstrated using a heat map ( Supplementary Figure 1D ). These results indicated that the four NFRGs play a role in the immune microenvironment of PD.

In the GSE7621 internal validation cohort, the expression of EN1 and LOXL1 was lower, while the expression of IRF7 and PLOD3 was higher in PD tissues than in normal controls ( Figure 6A ). The expression of the four signature genes in the GSE20163 external validation cohort was similar to that in the internal validation cohort, except for IRF7 ( Figure 6B ). The difference in results may be due to the small sample size. In the GSE49036 dataset consisting of patients with Parkinson’s disease at different Braak stages, there was significant difference in the expression of EN1, PLOD3 and LOXL1 in the different Braak stages 0 to 6 of Parkinson’s disease. However, there was no significant difference in expression of IRF7 among the different stages ( Figures 6C–F ). These results suggest that these NFRGs play a role in the pathogenesis and progression of Parkinson’s disease.

A risk-scoring model was developed based on the 12 NRFGs obtained in Figure 3 to identify potential prognostic biomarkers for GBM. The NFRGs with prognostic potential were subjected to LASSO regression analysis to reduce the number of genes in the final risk model. Ten NFRGs were identified from this step ( Figures 7A, B ). Multivariate Cox analysis identified 7 NFRGs, including EN1, TUBB2A, HSPB1, LOXL1, RGS4, L1CAM, and GPR143 as independent prognostic factors. Risk scores were calculated using the following formula: risk score = expression level of EN1* 0.17 + expression level of TUBB2A* 0.09 + expression level of HSPB1*0.14+ expression level of LOXL1*0.19 + expression level of RGS4*0.09 + expression level of L1CAM*0.08 + expression level of GPR143*0.20.

Patients in the TCGA cohort were classified into high-risk and low-risk groups based on the median risk score. Survival curves revealed that patients in the high-risk group had lower overall survival (OS) compared to the low-risk group in the TCGA and CGGA cohorts ( Figures 7C, D , P<0.05). Furthermore, the risk score was effective at predicting OS in the TCGA cohort. (AUCs for 1-, 3-, and 5-year OS were 0.734, 0.823, and 0.942, respectively; Figure 7E ). However, since GBM patients have dismal prognosis, the AUC values in the CGGA sample were not favorable ( Figure 7F ). In both the TCGA and CGGA cohorts, the area under the curve (AUC) for the risk score over three years was greater than the AUC values for the other clinicopathological features ( Figures 7G, H ). Risk maps were used to display survival results from the TCGA and the CGGA cohorts, while heat maps were used to display variations in the expression of the seven NFRGs across the various risk groups ( Figures 7I, J ).

PCA and t-SNE analyses were then performed using the NFRG classification in the expression profiles of the seven models. In the TCGA cohort ( Supplementary Figures 2A, B ) and the CGGA cohort ( Supplementary Figures 2C, D ), our signatures yielded results indicating a different distribution between the high-risk and low-risk groups. These findings imply that prognostic model can distinguish between high and low risk groups.

Univariate and multivariate Cox analyses showed that risk scores were independent prognostic factors for GBM patients compared to other common clinical characteristics ( Figures 8A, B ). In the CGGA cohort, results of both univariate and multivariate cox analyses showed that risk score was a prognostic factor independent of age, IDH mutation status, or MGMTp_methylation status ( Supplementary Table 1 ). To determine the clinical application of the risk models, age, sex, IDH mutation status, and risk score were included in a nomogram used to predict overall survival in patients with GBM based on the TCGA cohort ( Figure 8C ). We found that the risk score had the biggest influence in predicting OS, an indication that prognosis of GBM could be predicted using a risk model based on the seven NFRGs. At 1, 1.5, and 2 years, the calibration curves demonstrated a reasonable agreement between expected and observed values ( Figure 8D ). The three-year DCA curves ( Figure 8E ) and the temporal c-index values ( Figure 8F ) indicated that our model has the highest net benefit and that the risk model constructed based on the 7 NFRGs has more influence in clinical decision-making than the traditional model. The histogram of the chi-square test showed that risk grouping was only associated with whether IDH was mutated ( Figure 8G ). To validate these findings, we evaluated the relationship between risk score and clinical characteristic and found that individuals without IDH mutations were associated with higher risk scores ( Figures 8H–J ).

To determine the relationship between risk scores and immune cells and functions, we measured the enrichment scores of various immune cell subpopulations, associated activities, or pathways using the “cibersort” and “ssGSEA”. The low-risk group was associated with a larger infiltration of monocytes and M2-type macrophages ( Figure 9A ). In addition, the high-risk group showed a higher type 2 interferon response compared to the low-risk group, while the low-risk group had a higher type 1 interferon response ( Figure 9B ). The risk score was associated with the quantity of immune cells in the GBM tumor microenvironment determined by several methods using Spearman correlation analysis ( Figure 9C ). Furthermore, we discovered that a small number of immune checkpoints, namely CD48 and IDO1, were substantially expressed in the low-risk group compared to the high risk group ( Figure 9D ). These results imply that although patients in the high-risk group have a worse prognosis, they may be more responsive to immunotherapy due to their more active immune function. GSEA was used to investigate potential changes in biological function between risk groups based on the various prognoses of patients in the high-risk and low-risk groups. We chose the top 8 enriched signaling pathways based on normalized enrichment scores (NES) and p-values ( Figure 9E ). Surprisingly, lower risk scores were associated with Alzheimer’s disease and Parkinson’s disease, which is in line with the theme of our study.

To determine the molecular mechanisms driving the abnormal expression of these seven NFRGs, we explored the many histological levels, including genomes and copy numbers. Analysis of single nucleotide gene variant (SNV) data revealed that missense mutations in NFRGs were the most frequent variant categorization in the TCGA-GBM cohort, whereas single nucleotide polymorphism was the most common variation type. Among the SNV categories, C>T showed the highest prevalence ( Supplementary Figure 3A ). To summarize the ratio of pure and heterozygous mutations in the sample’s NFRGs, copy number variation was examined ( Supplementary Figure 3B ). We found that 17 GBM patients had mutations, with L1CAM mutations being the most common ( Supplementary Figure 3C ). Additionally, the Spearman’s correlation coefficient analysis of between copy number variations and gene expression showed that L1CAM copy number variations were downregulated in GBM whereas TUBB2A, HSPB1, LOXL1, RGS4, and GPR143 copy number variations were upregulated ( Supplementary Figure 3D ). Heterozygous variants of HSPB1 were present in most samples, and individual analysis showed that LOXL1 and L1CAM were copy number deletions. Whereas the pure-sibling mutation of GPR143 is mainly a copy number reduction ( Supplementary Figures 3E, F ), HSPB1 and L1CAM primarily amplified pure heterozygous mutations, suggesting that abnormal gene expression may be caused by both copy number variation and single nucleotide variation. The relationship between NFRGs expression and the activity of pathways linked to cancer was further examined. The findings demonstrated that NFRGs contributed to the inhibition of hormonal pathways in GBM patients and activation of the EMT, PI-3K-AKT, and TSC-mTOR pathways ( Supplementary Figure 3G ). We further explored the differential expression of NFRGs in the GDSC and Cancer Therapy Response Portal databases, their corresponding drug sensitivity ( Supplementary Figures 3H, I ). This suggested that the expression of the proposed risk profile genes may be exploited to develop agents for sensitizing drugs as well as predict chemotherapeutic drug sensitivity in patients.

In further experiments, we examined the correlation between risk score and tumor mutational load (TMB) ( Supplementary Figure 4B ) as well as differences in TMB among different risk subgroups ( Supplementary Figure 4A ). Results showed that the TMB was higher in the low-risk group. Thus, we generated two waterfall plots to explore the detailed mutational characteristics between high- and low-risk populations. The results indicated that PTEN, TP53, and TTN were the most commonly mutated genes in both risk groups ( Supplementary Figures 4C, D ).

Analysis of the violin plots designed to demonstrate the link between IPSs and risk groups, showed that high IPSs indicate stronger responses to PD-1 and CTLA-4 blockers ( Figures 10A–D ). Using the “pRRophetic” R package, we explored the potential sensitivity of clinical agents in the high-risk and low-risk groups. Agents available for the treatment of gliomas, such as nilotinib ( Figure 10E ), had a higher IC50 in patients of the high-risk group, whereas ABT737 and KU-55933 had a higher IC50 in patients of low-risk groups ( Figures 10F, G ). To understand the association between risk scores of NFRGs and the benefits of immunotherapy, we investigated a cohort of lung cancer patients treated with PD-1 checkpoint inhibitors (GSE135222) using the BEST database. The ROC curve analysis demonstrated that NFRGs was effective in predicting immunotherapy responsiveness, with low NFRGs expression score correlating with higher degree of immune response to anti-PD-L1 ( Figure 10H ).

Using the single-cell dataset GSE141982 from the TISCH database, we investigated the expression of 7 NFRGs in the GBM TME. It was observed that the GSE141982 dataset was enriched with several cell types of 16 cell populations and 4 cell subpopulations ( Supplementary Figure 5A ). Most endothelial cells, monocyte macrophages, and CD8+ T cells expressed HSPB1 and TUBB2A. The expression of the other NFRGs, which are primarily found in tumor cell cells, was low ( Supplementary Figures 5B, C ).

By analyzing the neurotrophic factor-related genes in PD and GBM, we found that EN1 and LOXL1 we important players in both diseases. In the human protein atlas, the protein expression of EN1 ( Figures 11A, B ) and LOXL1 ( Figures 11C, D ) were higher in GBM relative to normal cortical tissue. RT-qPCR results confirmed the higher expression of EN1 and LOXL1 in both GBM cell lines ( Figures 11E, F ).