Machine learning screening for Parkinson’s disease-related cuproptosis-related typing development and validation and exploration of personalized drugs for cuproptosis genes

Introduction

Parkinson’s disease (PD) is a progressive neurodegenerative disorder with a high prevalence. According to statistics in 2016, over 6 million people had PD worldwide and the number of people with PD is estimated to increase with the age of the population (1). Patients with PD often require the companionship of caregivers, most of whom are under great stress, and this service imposes a heavy socio-economic burden on society (2). The clinical features of the disease include motor symptoms, including bradykinesia, resting tremor, and rigidity, and non-motor symptoms, including loss of smell, autonomic dysfunction, depression, cognitive impairment, and insomnia. The main pathological features of PD are progressive loss of nigrostriatal dopaminergic neurons and abnormal α-synuclein aggregation called Lewy vesicles (3), but the exact etiology of PD is still unclear.

In recent years, the pathological process involved in PD-associated neuronal death has been thought to be apoptosis, and a variety of cell death and mechanisms have been found to be involved in PD pathogenesis, including autophagy, necroptosis, cell death-dependence, iron death, and cell scorching (4-7). It has been shown that corynoxine derivative CB6 inhibits PD progression in mice by inducing PIK3C3 complex-dependent autophagy, whereas metformin ameliorates depression-like behavior in experimental PD by inducing autophagy in the substantia nigra and hippocampus (8,9). There is evidence that excessive copper (Cu) administration leads to a 40–60% reduction in the binding of dopamine (DA) D2 receptors to (3H)-spiperone and results in neuronal cell death and α-synuclein aggregation (10). This confirms the possibility that Cu-induced thiol modifications are involved in an important mechanism of PD progression. Another study found that Cu can be involved in the pathogenesis of PD by catalyzing harmful redox reactions including oxygen derivatives that increase oxidative stress (11). Not only this, but Cu in the brain can be directly responsive to neurotransmitters, leading to the oxidation of DA, which in turn promotes the progression of PD (12). Recently, Tsvetkov et al. revealed that Cu toxicity in human cells occurs through the direct binding of excess Cu to the lipid acylated component of the tricarboxylic acid (TCA) cycle. This leads to the accumulation of lipid acylation-related proteins, loss of iron-sulphur cluster proteins, and ultimately cell death due to intracellular proteotoxic stress, a new form of cell death known as ‘cuproptosis’ (13). Therefore, modulating Cu ions and lipid acylation in the TCA cycle may be a potential therapeutic strategy to improve the prognosis of PD.

Integrated bioinformatics analysis is now widely used to identify potential novel biomarkers and their role in various diseases, A series of bioinformatic analyses using expression profiles of peripheral blood mononuclear cells (PBMC) have been performed to identify differentially expressed genes associated with humoral immunomodulatory mechanisms between PD and healthy controls (14), but a relevant and novel comprehensive bioinformatic analysis of PD patients and healthy controls in the substantia nigra densa (SNpc) is lacking. We aim to screen potential PD patients at the genetic level and explore novel mechanisms of PD pathogenesis and initially explore novel therapeutic targets. We collected a total of 19 genes associated with cuproptosis by reviewing the literature (13,15). The GSE49036 dataset, which contains 8 control brain samples and 20 PD brain samples, was used to look for differentially expressed cuproptosis genes and explore their immune infiltration patterns with immune-related genes and PD, and also further explored these cuproptosis-related differentially expressed genes (CR-DEGs)-targeted small molecule drugs. The PD patients were further clustered, and the clinical and clustering key module genes were obtained using weighted gene co-expression network analysis (WGCNA). The most relevant diagnostic marker genes for PD were screened by machine learning and receiver operating characteristic (ROC) curve analysis was used to distinguish patients with PD from healthy controls. The diagnostic efficacy was finally validated using GSE7621 and GSE20141. In conclusion, this study explored the possible mechanisms of action of cuproptosis-related genes (CRGs) in PD and the combined immune infiltration pattern as a new therapeutic target for them. We present the following article in accordance with the TRIPOD reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-22-5756/rc).

Methods

Data download

Three tissue sequencing datasets of RNA microarray data on PD nigrostriatal samples were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/): GSE7621 (9 control brain samples and 16 PD brain samples), GSE20141 (8 control brain samples and 10 PD brain samples), and GSE49036 (8 control brain samples and 20 PD brain samples). These datasets were analyzed using the Affymetrix Human Genome U133 Plus 2.0 array on the GPL570 (Affymetrix GeneChip Human Genome U133 Plus 2.0 Array) platform. The above data were extracted and transformed with gene names using Perl software (https://www.perl.org/). We then screened 19 Cu amination-associated genes from previous literature (16). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Data processing and CR-DEGs screening

To find differences in the expression of cuproptosis genes in the substantia nigra of PD patients compared to normal controls, CRGs were screened from the GSE49036 dataset using the limma package of R software (The R Foundation for Statistical Computing, Vienna, Austria) and CR-DEGs screening was performed on 20 PD samples and 8 control brain tissues. ggpubr package for box plot and heat map presentation of the results.

Gene set enrichment analysis

To understand the impact of CR-DEGs on the mechanisms of PD disease development, we performed a single-sample gene set enrichment analysis (ssGSEA), using the limma package, org.Hs.eg.db, and clusterProfiler package, to an analysis and visualization of the diagnostic genes.

Investigation of immune cell infiltration

To understand the immune microenvironment in PD, we used a previously established method (17) to quantify immune infiltration and associated immune functions by ssGSEA, which calculates an enrichment score representing the level of immune cell infiltration and immune-related pathway activity. The “ggplot2” package was used to create a heat map of the distribution and variation of immune cells.

Competing endogenous RNA network search

To find possible messenger RNA-microRNA-long noncoding RNA (mRNA-miRNA-lncRNA) network relationships for CR-DEGs, we firstly searched MiRanda (http://www.microRNA.org), miRDB (http://www.miRDB.org/), and TargetScan (http://www.TargetScan.org/vert_71/) to find miRNAs for PD-CR-DEGs. The screening criteria were that the miRNAs predicted by all 3 databases were the mRNA-miRNAs we were looking for, and then the lncRNAs corresponding to the miRNAs were found by using the spongeScan website (18) (spongeScan: a web for detecting microRNA binding elements in lncRNA sequences). Finally, the predictions were visualized using Cytoscape 3.8.2 (https://cytoscape.org/).

Drug screening and molecular docking

We used the DGIdb database (https://dgidb.genome.wustl.edu/) for drug prediction of differential genes in order to explore possible targeting small molecule drugs for CR-DEGs. Using the search drug-gene interactions option on the web page, and by entering PD-CR-DEGs into the prediction, results were obtained and downloaded from the webpage by censoring the gene-drugs reported in the literature and visualizing the results of the gene-drug prediction network using Cytoscape 3.8.2. The top 5 gene-drug interaction scores were selected by eliminating chemotherapeutic drugs from the above results and using AutoDock to perform small-molecule drug-large-molecule protein docking, as in our previous studies (19,20). The energy required for docking was heat mapped using the R package pheatmap.

Genotyping and enrichment analysis

First, we found CR-DEGs by the methods described above. We used the R package to classify patients into different molecular subtypes based on the expression of CRGs. The R package ConsensusClusterPlus was used to remove the control samples from the GSE49036 dataset, retaining only the experimental samples. We set the cluster maxK =9, selected the most significant clusters by consistency scoring, and used principal component analysis (PCA) to determine whether the clusters were significant. The R software Limma and pheatmap packages were then used to explore the differential PD CR-DEGs after clustering, and to investigate CR-DEGs with molecular subtype salient features, we performed a gene set variation analysis (GSVA) using the marker gene set (c2. c2.cp.kegg.symbols.gmt, c5.go.symbols.gmt) from the MSigDB database derivation. We then used the Profiler package in R to analyze the enrichment features of the genes. The results were visualized using the “ggplot2” package.

Weighted co-expression network construction

We constructed a co-expression network of DEGs using the WGCNA package in R. Next, Pearson’s correlation matrix was performed on pairs of genes. To further analyze the sample clustering to detect outliers, we calculated the dissimilarity of the module feature genes, selected tangents to the module dendrogram, and merged some modules. Module assignments determined by dynamic tree cuts were shown using the colored rows below the dendrogram. By analyzing the module-trait relationships between module trait genes and clinical traits, two modules associated with a specific trait were identified. The same approach was used to cluster type the two modules associated with a specific trait and their most significant trait modules were Venn diagrammed to find the clinical traits and the common trait genes that were significant after clustering.

Machine learning methods for filtering feature model construction and validation

After clustering, the above obtained clinical traits and significant common trait genes were used to construct a random forest (RF) tree model, support vector machine (SVM) learning model, extreme gradient boosting (XGBoost) model, and general linear model (GLM) using R software, recursive feature elimination (RFE) in Random Forest Algorithm is a supervised machine learning method for sequencing genes associated with copper poisoning in PD. The prediction performance was estimated by ten-fold cross-validation, and the genes whose relative importance was >0.25 were identified as the characteristic genes. SVM is a small sample learning method, which basically bypasses the traditional process of induction to deduction, realizes efficient “reverse reasoning” from training to prediction of samples, and simplifies common classification and regression problems. XGBoost is a representative algorithm based on integrated lift, which complements the overfitting problem of gradient lift model. The GLM is an extension of the traditional linear model, which is an algorithm in which the population mean is passed by a nonlinear join function to better process and obtain non-normally distributed data (21). Using the four methods described above to construct a diagnostic model with characteristics, the data samples were first subjected to residual analysis, and to plot the reverse cumulative distribution of residuals for the four methods. The ROC curves of the four methods were plotted, and the best feature models for machine learning were selected based on the reverse cumulative scores of the residuals and the area under the curve (AUC) values under the ROC curves. The best predictive models were plotted on column line plots, calibration curves, and decision curves using packages such as rms and rmda. Finally, the ROC curves were used to verify the diagnostic performance of the models in the GSE7621 and GSE20141 datasets.

Statistical analysis

All statistical analyses were performed using R version 4.1.0, One-way ANOVA was performed to compare the differences among multiple groups (≥2 groups). The Student’s t-test was used to compare the differences between the two groups. Statistical significance was set at P<0.05.

Results

Study process

Figure 1 show all the processes analyzed. First, PD-brain nigrostriatal sequencing data were retrieved from the GEO database, and the GSE49036 dataset was used to screen CR-DEGs and explore their correlation and single-gene GSEA pathway analysis. Correlations between individual genes and immune cells were then explored, a CR-DEGs mRNA-miRNA-lncRNA network was constructed, a gene-drug network was constructed, and molecular docking was performed to explore the feasibility of their binding. The molecular subtypes of CR-DEGs present in PD were then explored, WGCNs with clinical and clustering typing were constructed, and the clustered typing was crossed over to the clinical module genes. The crossed genes were then used to construct a diagnostic model using machine learning and the diagnostic model was proved using GSE7621 and GSE20141.

Figure 1 Flow chart of the systematic analysis of cuproptosis-related diagnostic genes and medication of PD. In WGCNA selection, red represents a positive association with PD, while blue represents a negative association with PD. In the energy heat map, the redder the color, the greater the energy required for intermolecular binding. In PD-immunocorrelation analysis diagnostic genes, red represents high expression in PD tissue, while blue represents low expression in PD tissue. *, P<0.05; **, P<0.01; ***, P<0.001. GEO, gene expression omnibus; CR-DEGs, cuproptosis-related differentially expressed genes; PD, Parkinson’s disease; WGCNA, weighted gene co-expression network analysis; AUC, area under the curve; ROC, receiver operating characteristic; RF, RF random forest tree; SVM, support vector machine; XGB, Tree Ensemble; GLM, generalized linear model.

Identification of CR-DEGs in PD and correlation analysis between genes

The GSE49036 data were normalized and CRGs were extracted using the limma package of R software and screened for DEG analysis. According to the screening criteria: P<0.001, indicated by “***”, P<0.01, indicated by “**” and P<0.05, indicated by “*”, heat maps (Figure 2A) and box plots (Figure 2B) were drawn using pheatmap and ggpubr software packages to visualize the differential gene expression results.

Figure 2 CR-DEGs in PD. (A) Heat maps of the 3 CR-DEGs in PD; the color of the columns represents the gene expression of the tissue samples, with red indicating high expression and green indicating low expression. (B) The CR-DEGs boxplot in PD. *, P<0.05; **, P<0.01. CR-DEGs, cuproptosis-related differentially expressed genes; PD, Parkinson’s disease.

Correlation analysis of expression of the above-mentioned genes in PD patients was performed using the corrplot package of R software and the results showed that ATP7B was negatively correlated with NFE2L2 and MTF1 and that NFE2L2 was positively correlated with MTF1. The results of expression correlation between genes were also visualized using the circlize package (Figure 3).

Figure 3 CR-DEGs in correlations. ATP7B was negatively correlated with NFE2L2 and MTF1, while NFE2L2 was positively correlated with MTF1. CR-DEGs, cuproptosis-related differentially expressed genes.

GSEA of PD tissue and CR-DEGs

Single-gene GSEA enrichment analysis was performed to explore the potential mechanisms of action of cuproptosis-related genes such as NFE2L2, MTF1 and ATP7B in PD, as shown in Figure 4. Interestingly, all 3 genes were associated with the synaptic vesicle cycle, although they were down-regulated in NFE2L2 and MTF1, down-regulated in ATP7B. In contrast, the ATP7B up-regulated signaling pathway is also related to Citrate cycle (TCA cycle), Nicotine addiction, Retinol metabolism, Starch and sucrose metabolism, etc. whereas NFE2L2 and MTF1 were jointly associated with Proteasome, and MTF1 down-regulation signaling pathway was associated with 2-oxocarboxylic acid metabolism, glycosphingolipid biosynthesis—globo and isoglobo series, and so on. NFE2L2 down-regulation was associated with aminoacyl-tRNA biosynthesis, glycosphingolipid biosynthesis—globo and isoglobo series, and so on.

Figure 4 ATP7B, MTF1 and NFE2L2 pathway correlation exploration. (A) GSEA of ATP7B in the pathway; (B) GSEA of MTF1 in the pathway; (C) GSEA of NFE2L2 in the pathway. GSEA, gene set enrichment analysis.

These results indicate that NFE2L2, MTF1, and ATP7B may be involved in the regulation of the synaptic vesicle cycle, which in turn affects brain neurotransmitter transmission and may be involved in aerobic metabolism of brain tissue, leading to disease progression.

Immunocorrelation analysis of PD tissue and CR-DEGs

To explore the immune cell infiltration of CR-DEGs in PD, we performed GSEA of dataset GSE49036 using the immune gene set “immune.gmt”, as shown in Figure 5A,5B. CD8 T cell was highly expressed in normal controls with PD tissue, whereas plasmacytoid dendritic cell (DC), macrophage, immature B cell, myeloid-derived suppressor cells (MDSCs), mast cell, neutrophil, natural killer (NK) T cell, type 2 T helper cell, T follicular helper cell, and other immune cells were highly expressed in PD substantia nigra versus normal brain tissue. We also explored the relationship between individual genes of CR-DEGs and immune cells (see Figure 6). ATP7B was associated with activated B cell, CD56 dim NK cell, central memory CD4 T cell, central memory CD8 T cell, macrophage, mast cell, MDSC, monocyte, NK cell, NK T cell, plasmacytoid DC, and T follicular helper cell, and other immune cells were negatively correlated, whereas NFE2L2 was negatively correlated with MTF1 and type 17 T helper cell, plasmacytoid DC, NK T cell, central memory CD8 T cell, NK cell, monocyte. Memory B cell, mast cell, and activated B cell were positively correlated, and interestingly, they were negatively correlated with effector memory CD8 T cell expression; it is possible that these immune cells are involved in the midbrain substantia nigra of PD degenerative death of dopaminergic neurons, and CR-DEGs such as NFE2L2, MTF1, and ATP7B may be promising targets for immunotherapy in PD.

Figure 5 Immunogene set analysis of immune infiltration patterns in PD and normal conditions. (A) Heat maps of immunogene set analysis of PD samples and the color of the columns represents the gene expression of the tissue samples, with red indicating high expression and green indicating low expression. (B) Violin plot showing the difference in infiltrating immune cells between the two groups. PD, Parkinson’s disease; MDSC, myeloid-derived suppressor cell.

Figure 6 Immunogene set analysis of immune cell in ATP7B, MTF1, and NFE2L2; heat map of ATP7B, MTF1, and NFE2L2 correlation with immune cells. MDSC, myeloid-derived suppressor cell.

PD-CR-DEGs ceRNA network search

To further explore the complex molecular interaction mechanism between the NFE2L2, MTF1, and ATP7B genes, we created an mRNA-miRNA-lncRNA network relationship diagram, as shown in Figure 7. This complex network has 60 nodes and 99 interactions, including 40 miRNAs and 17 lncRNAs, which may provide us with a new idea to explore the pathogenesis of cuproptosis-related PD.

Figure 7 ceRNA network in ATP7B, MTF1, and NFE2L2. Red circles represent mRNAs, green triangles represent miRNAs, and purple diamonds represent lncRNAs. mRNAs, messenger RNAs; miRNAs, microRNAs; lncRNAs, long noncoding RNAs.

Prospective small molecule drug for PD therapy

In order to explore Prospective small molecule drugs for PD therapy, downloaded the results of NFE2L2, MTF1, and ATP7B targeting drugs from the DGIdb database and selected drugs with documented results (see Table 1) and presented the gene-drug results using Cytoscape (Figure 8A). We finally selected the top 5 NFE2L2 gene-drug interaction scores for small molecule drug-macromolecule protein docking using AutoDock. The energy required for docking is shown in Figure 8B, with the lowest energy required for docking LAGASCATRIOL to NFE2L2. Figure 9A shows a visualisation of the overall picture of the molecular docking results of the NFE2L2 protein in LAGASCATRIOL, Figure 9B shows a local view of the molecular docking results visualising the small molecule drug with the protein, see Figure 9C. LAGASCATRIOL is able to form 1 ionic bond each with arginine at position 456 and 504 of NFE2L2, and 2 ionic bonds with aspartic acid at position 457.

Figure 8 Search for ATP7B, MTF1 and NFE2L2 targeted drugs AND Molecular docking energy heat map. (A) The ATP7B and NFE2L2 targeted drugs; (B) molecular docking of NFE2L2 proteins in andalusol, chembl128729, dimethyl_fumarate, lagascatriol, suforaphane_isoselenocyanate, the ordinate value is the energy of docking (Kcal/mol). The redder the color, the greater the energy required for intermolecular binding.

Figure 9 NFE2L2 protein docking with LAGASCATRIOL. The molecules and drugs with the lowest docking energies are visualized in Figure 8. (A-C) Visualization of molecular docking results for NFE2L2 proteins in LAGASCATRIOL.

Clustering based on the CR-DEGs in the PD sample

In GSE49036, cluster analysis was performed to classify the PD samples into different molecular subtypes and examine the subtypes using PCA based on the expression of the 3 CR-DEGs selected above (Figure 10A,10B). We named the 2 clusters as CRGClusters C1 and C2, respectively. Differential analysis of the 2 CRGClusters showed that ATP7B expression was significantly higher in C1 compared to C2, whereas NFE2L2 and MTF1 in C1 had significantly higher ATP7B expression compared to C2 expression, which was significantly decreased (Figure 10C). In the Kyoto Encyclopedia of Genes and Genomes (KEGG), comparison between the two clusters in the GSEA also revealed significantly different KEGG pathway enrichment profiles, interestingly on PARKINSONS_DISEASE, but also on ubiquitin-mediated proteolysis, oxidative phosphorylation, and citrate cycle tca cycle, which were up-regulated, whereas tgf beta signaling pathway and spliceosome were down-regulated (Figure 10D). In contrast, among the Gene Ontology (GO) biological down-regulated in pathways: immunological synapse formation, positive regulation of hemopoiesis, positive regulation of b cell differentiation, regulation of endodermal cell differentiation, regulation of synaptic vesicle endocytosis, dendrite extension, presynaptic dense core vesicle exocytosis, regulation of postsynaptic vesicle exocytosis, and regulation of postsynaptic cytosolic calcium ion concentration (Figure 10E). These results suggest a possible involvement in the production and transmission of neurotransmitters in the brain, and possibly in dendritic production related to the progression of PD by causing a significant reduction in striatal DA content.

Figure 10 Two different subtypes of PD divided by consensus clustering and GSVA-GO and KEGG enrichment analysis (A) defined two clusters (k=2) and the consensus matrix heat map of their related regions. (B) PCA analysis showed significant transcriptome differences between the two subtypes. C1 is shown in green, C2 in red, and there is no sample overlap between them. (C) Box plots showed that ATP7B, MTF1 and NFE2L2 genes were significantly different between class A and B. Asterisks indicate statistical P values (**P<0.01; ***P<0.001). (D) GSVA-GO enrichment analysis N in the two clusters, where red and blue represent BP, CC, MF in activated GO and BP, CC, MF (E) in inhibited GO, respectively, and GSVA-KEGG enrichment analysis in red and blue represent activation pathway and inhibition pathway. BP, biological process; CC, cell component; MF, molecular function; RIA, ruptured intracranial aneurysm; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PD, Parkinson’s disease; GSVA, gene set variation analysis; PCA, principal component analysis.

WGCNA co-expression analysis of PD samples

We first conducted a WGCNA analysis of clinical trait correlations in order to select the top 25% of genes with the greatest fluctuations for WGCNA analysis and cut off outlier samples, with the remaining samples being included in the analysis (Figure 11A). As shown in Figure 11B, scale independence reached 0.9 when power =8 and mean concordance was high. Therefore, power =8 was used to construct co-expression modules to obtain preliminary module delineation results, and the results of WGCNA showed that different modules were identified with different colors (Figure 11C). To detect outliers, trees were constructed using the eigenvalues of the modules, and then the distances of the trees belonging to the same branch were very close were merged and the intercept value was set to 0.5 (Figure 11D). Co-expression modules were constructed as shown in Figure 11E and the results were obtained after merging similar modules. Based on the eigenvalues of the samples in each module and the characteristics of the samples, correlation analysis was performed to identify modules associated with specific traits, and the genes represented by MEblack were highly positively correlated in PD (Figure 11F).

Figure 11 WGCNA selection of PD disease-related modules co-expressed genes. (A) Outliers were detected in the sample cluster. (B) The cut-off point was set as 0.9, and the soft threshold power was set as β=8. (C) Tree diagram of all DEGs based on the cluster of difference measurement. The colored bands show the results from the automated monolithic analysis. (D) Correlation diagram between modules obtained by clustering according to inter-gene expression levels. (E) Heat map of the correlation between module characteristic genes and phenotypes. We chose the MEblack module for subsequent analysis (the ordinate value is the correlation coefficient of feature module). (F) Module membership vs. gene significance, cor is disease correlation coefficient, P is the P value of correlation coefficient. ME, module eigengene; WGCNA, weighted gene co-expression network analysis; PD, Parkinson’s disease; DEGs, differentially expressed genes.

In the same way that we performed WGCNA analysis of correlations after typing, we cut off outlier samples (Figure 12A) as in Figure 12B when power =7, scale independence reached 0.9, modules were identified by color and eventually found to be characterized after typing (Figure 12C-12F), with the genes represented by MEyellow being highly positively correlated in CRGCluster C1 (Figure 12F).

Figure 12 WGCNA selection of PD cluster-related modules co-expressed genes. (A) Outliers were detected in the cluster 1 and 2-sample cluster. (B) The cut-off point was set as 0.9, and the soft threshold power was set as β=7. (C) Tree diagram of all cluster 1 and 2-sample DEGs based on the cluster of difference measurement. The colored bands show the results from the automated monolithic analysis. (D) Correlation diagram between modules obtained by clustering according to inter-gene expression levels. (E) Heat map of the correlation between module characteristic genes and cluster 1 and 2-sample phenotypes. We chose the MEyellow module for subsequent analysis (the ordinate value is the correlation coefficient of feature module). (F) Module membership vs. gene significance, cor is disease correlation coefficient, P is the P value of correlation coefficient. ME, module eigengene; WGCNA, weighted gene co-expression network analysis; PD, Parkinson’s disease.

Machine learning builds diagnostic models

We crossed the clinical trait and typed common trait gene modules to find a total of 72 trait genes (Figure 13A). We used the GSE49036 dataset to extract the 72 trait genes and used the RF tree, SVM, XGB, and GLM algorithms to build diagnostic models with the traits. We first analyzed the data samples for residuals, the reverse cumulative distribution of the residuals of the four algorithms is plotted (Figure 13B). The XGB model has smaller residuals for most samples, while the GLM has higher residuals for many samples than the XGB model. The box plots of the residuals of the 4 algorithms (Figure 13C), with the red dots indicating the root mean square of the residuals, were ranked as follows: XGB < SVM < RF < GLM, and the 4 models were validated with ROC curves. The AUC values of the four models were as follows: RF =0.833; SVM =0.833; XGB =0.917; GLM =0.833 (Figure 13D).

Figure 13 Machine learning builds models can largely progress between PD and normal samples. (A) The intersection genes of disease feature module and typing feature module were analyzed by WGCNA. (B) RF random forest tree, SVM, XGB, GLM and other algorithms were used to construct the reverse cumulative distribution map of model residuals. The Y-axis value represents the percentile of the outlier. (C) Box plots of sample residuals of the four algorithms. The X-axis value represents the quantile of outliers, and the red dot represents the mean. (D) ROC analysis of RF random forest tree, SVM, XGB, GLM, and other algorithms. The ROC values of RF and SVM were 0.833. (E,F) ROC analysis of XGB model in the GSE7621 and GSE20141 validation set. x-axis: FPR, false positive rate; y-axis: TPR, true positive rate; WGCNA, weighted gene co-expression network analysis; AUC, area under the curve; RF, RF random forest tree; SVM, support vector machine; XGB, Tree Ensemble; GLM, generalized linear model; ROC, receiver operating characteristic; CI, confidence interval.

Finally, we chose the XGB model and used the GSE7621 and GSE20141 gene sets to validate the ROC curves for the XGB model (Figure 13E,13F), with corresponding AUC values of 0.688 and 0.833. Finally, the genes in the model were plotted in a nomogram (Figure 14), which showed that the genes in the model are SLC35D3, FLJ22184, SYT17184, SLC35D3, FLJ22184, FLJ22184, SYT17, and RIMS3.

Figure 14 A nomogram was constructed to predict the probability in PD progress. The values of each variable (SLC35D3, FLJ22184, SYT17, and RIMS3) are summed to obtain a total score. Each variable corresponds to a score, and the total score can be calculated by summing the scores of all variables. PD, Parkinson’s disease.

Discussion

PD is considered a progressive neurodegenerative disease that manifests clinically as resting tremor, rigidity, bradykinesia, or dyskinesia, postural gait abnormalities, and other non-motor symptoms. PD is now understood to be a multisystem disorder accompanied by significant nerve inflammation and immune function dysfunction. Its pathogenesis is currently unclear and includes the presence of oxidative stress, protein aggregation, mitochondrial dysfunction, and inflammation as the cause of neuronal death (22). Among the pathological features of PD are the formation of predominantly intracellular alpha-synuclein-containing inclusion bodies, called Lewy vesicles. The current understanding of neural death in PD is the role of α-synaptic nuclear proteins in brain cell death through the formation and assembly of these protein aggregates. In contrast, Cu¹⁺ can also bind directly to two different regions of the α-synaptic nucleoprotein thereby causing oxidative neurological damage degeneration (23). High concentrations of Cu promote α-synaptic protein misfolding leading to PD progression (24). Alternately, 6-OHDA disrupts copper (Cu) homeostasis through dysregulation of cellular Cu transport, leading to neuronal cell dysfunction and cell death (25). Therefore, Cu accumulation and Cu metabolism are important for the development of PD disease, and in this study we explored potential central genes related to cuproptosis and immune infiltration patterns in the substantia nigra of PD patients through bioinformatics analysis, and we found that plasmacytoid DC was upregulated in tissues with the most significant differences in PD. It has been suggested that DC-mediated T-cell differentiation leads to a reduction in recognized antigenic peptides, resulting in the timely degradation of inflammatory proteins associated with α-synuclein synthesis in PD (26), whereas our immune cell analysis of three CR-DEGs found that they were all associated with NK cells, and it has been shown that increased NK cell numbers and cytotoxicity during PD progression can play a neuroprotective role (27), which may further reveal new targets for PD immunotherapy and provide new potential therapeutic biomarkers and treatment strategies for PD.

Our study used three PD patients’ nigrostriatal sequencing datasets to enrich the exploration of cuproptosis genes, firstly to find DEGs related to cuproptosis in PD patients, namely ATP7B, NFE2L2, and MTF1, and to explore their respective possible pathways of action and relationships with individual immune cells. We also explored their possible target miRNAs and lncRNA, mapped the complex network of ceRNAs, and searched for their corresponding small molecule drugs through the DGIdb database. We also conducted an in-depth analysis of PD patients by clustering and typing to explore the possible mechanisms of tissue typing in PD tissues, and constructed a diagnostic model of PD using WGCNA and machine learning methods.

Our results of GSEA revealed that the genes ATP7B, NFE2L2, and MTF1 were all associated with the synaptic vesicle cycle, which may suggest that these three genes could be involved in the synaptic vesicle cycle to reduce neurotransmitter transmission. ATP7B is a key ATPase in the Golgi complex harbors copper-transporting, which transfers Cu from the cytoplasmic matrix to the Golgi lumen for incorporation into the Cu-dependent enzyme. The Golgi complex also sends these ATPases to appropriate post-Golgi destinations to ensure the correct Cu flux in vivo and to avoid potentially toxic Cu accumulation. Mutations in ATP7A or ATP7B or the proteins that regulate their transport can affect their exit from the Golgi or subsequent fetch back to this organelle (28). α-synuclein assemblies need to exhibit seeding capacity, namely, template capacity for further aggregation, and toxicity mediated at least in part by interference with synaptic vesicle or organelle homeostasis (29). NFE2L2 encodes nuclear factor E2 and related factor 2 (NRF2), which are master regulators of antioxidant enzymes that prevent oxidative stress and inflammation, and perillyl is a key inhibitor of α-synuclein protein synthesis, thereby slowing PD progression. In contrast, galangin attenuates neuronal apoptosis via the Nrf2/Keap1 pathway, thereby slowing PD progression (30,31). MTF1, a key regulator of heavy metal homeostasis and detoxification, provides complementary functions in oxidative stress and other cellular stress responses. MTF1 gene polymorphisms have the potential to influence further exacerbation of PD (32). These studies above suggest that genes such as ATP7B, NFE2L2, and MTF1 may have important implications for the pathogenesis and progression of PD, based on which they are expected to be our next therapeutic targets for the study of PD.

The clinical diagnostic criteria for PD are bradykinesia with resting tremor, rigidity, or both (33). In contrast, the location of the DA transporter has been identified in diagnostic imaging primarily by DA transporter single photon emission computed tomography (DaT SPECT) by showing presynaptic DA neuronal hypofunction in PD and other neurodegenerative diseases. Identification of presynaptic dopamine neuronal dysfunction present in PD and other neurodegenerative diseases is used to differentiate PD from primary tremor, but it is not necessary to routinely perform these scans (34). They cannot distinguish PD from other neurodegenerative diseases that also involve DA transporter dysfunction (e.g., multisystem atrophy, progressive supranuclear palsy). There is growing evidence that PD includes a heterogeneous subtype have implications for diagnosis, prognosis, and expected response to treatment. Initial subtyping focused on motor characteristics (35). One typing approach includes three groups: predominantly mild motor: younger age of onset, mild motor and non-motor symptoms, slow progression and good drug response (36). Therefore, the typing of PD may be a better guide for treatment. We also typed PD and divided it into CRGClusters C1 and C2. Our results showed that ATP7B expression was significantly higher in C1 compared to C2, whereas NFE2L2, MTF1, and ATP7B expression was significantly decreased in C1 compared to C2. Our post-typing GSVA-KEGG enrichment analysis showed that in fractions, C2 could be enriched in PD signal path, as well as in the GO biological pathways: immunological synapse formation, regulation of endodermal cell differentiation, and regulation of synaptic vesicle endocytosis, dendrite extension, presynaptic dense core vesicle exocytosis, regulation of postsynaptic cytosolic calcium ion concentration. We obtained a diagnostic line graph by machine learning methods, so that we can diagnose PD by detecting the expression of genes; our research facilitates us to personalize the diagnosis and guide the treatment of PD by typing.

Currently, DA antagonists remain the first-line drugs used clinically for the treatment of PD, and the use of levodopa improves more function but increases the risk of dyskinesia, especially at higher doses (37). If it is necessary to reduce the dosage of levodopa after the development of dopamine-dysregulation syndrome (DDS), we can in turn try to switch to other drugs, which still have their own side effects, such as hallucinations and delirium (38). In contrast, levodopa equivalent daily dose (LEDD) can usually be reduced by 30–50% after deep thalamic nucleus brain stimulation (DBS) surgery is administered to patients, and one study showed that 7% of 150 patients treated with STN-DBS stopped taking the drug after 3 years (39). Complete cessation after STN-DBS dopaminergic medication is usually not possible because the lack of DA in the limbic system and associated circuits can lead to apathy and depression (39,40), so alternative medications with fewer side effects and effective alternatives need to be chosen. We searched the DGidb database for andalusol, chembl128729, dimethyl_fumarate lagascatriol, and suforaphane_isoselenocyanate, and other small molecules were identified through the DGidb database, and molecular docking was performed to initially explore the amino acid residue linkage sites between these drugs and NFE2L2 and their counterparts. There are studies on the Nrf2-dependent neuroprotective activity of diterpenoids isolated from Sideritis spp., counteracting oxidative damage and their potential role as useful agents for the prevention of those with oxidative stress-mediated dementia. ANDALUSOL was the most active of the diterpenoids studied (41). The neuroprotective effect of dimethyl_fumarate in the MPTP-mouse model of PD (42), the findings of the above-mentioned studies led us to search for these small molecules, which provide a new idea for the next targeted drugs for the treatment of PD.

There were some limitations to our study. Firstly, our findings were derived by computer using in-depth calculations using bioinformatics methods and still require in vitro experiments to validate the results. Secondly, further animal models need to be constructed to collect samples of the specific mechanisms by which the central genes associated with cuproptosis contribute to PD, and initial validation of the animal models of PD using our personalized drugs will be required, and ultimately we will use samples from patients in the early stages of PD for validation, combined with clinical data to provide diagnostic biomarkers and timely drug interventions. In future studies, we will focus on the mechanisms of action of ATP7B, NFE2L2, and MTF1 in PD and explore early diagnostic biomarkers for PD.

Conclusions

Using bioinformatics analysis, CR-DEGs were identified in the substantia nigra of PD patients. ATP7B, together with NFE2L2 and MTF1, were considered as candidate genes for further study. This study provides new genes associated with Cu accumulation leading to cell death in the pathogenesis of PD and explores promising therapeutic targets for related neurodegenerative diseases.

Acknowledgments

Funding: None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://atm.amegroups.com/article/view/10.21037/atm-22-5756/rc

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://atm.amegroups.com/article/view/10.21037/atm-22-5756/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Global, regional, and national burden of neurological disorders, 1990-2016: a systematic analysis for the Global Burden of Disease Study 2016. Lancet Neurol 2019;18:459-80. [Crossref] [PubMed]
2. Zhong M, Peppard R, Velakoulis D, et al. The relationship between specific cognitive defects and burden of care in Parkinson's disease. Int Psychogeriatr 2016;28:275-81. [Crossref] [PubMed]
3. Bloem BR, Okun MS, Klein C. Parkinson's disease. Lancet 2021;397:2284-303. [Crossref] [PubMed]
4. Outeiro TF, Grammatopoulos TN, Altmann S, et al. Pharmacological inhibition of PARP-1 reduces alpha-synuclein- and MPP+-induced cytotoxicity in Parkinson's disease in vitro models. Biochem Biophys Res Commun 2007;357:596-602. [Crossref] [PubMed]
5. Koprich JB, Reske-Nielsen C, Mithal P, et al. Neuroinflammation mediated by IL-1beta increases susceptibility of dopamine neurons to degeneration in an animal model of Parkinson's disease. J Neuroinflammation 2008;5:8. [Crossref] [PubMed]
6. Wu JR, Wang J, Zhou SK, et al. Necrostatin-1 protection of dopaminergic neurons. Neural Regen Res 2015;10:1120-4. [Crossref] [PubMed]
7. Mahoney-Sánchez L, Bouchaoui H, Ayton S, et al. Ferroptosis and its potential role in the physiopathology of Parkinson's Disease. Prog Neurobiol 2021;196:101890. [Crossref] [PubMed]
8. Mendonça IP, de Paiva IHR, Duarte-Silva EP, et al. Metformin improves depressive-like behavior in experimental Parkinson's disease by inducing autophagy in the substantia nigra and hippocampus. Inflammopharmacology 2022;30:1705-16. [Crossref] [PubMed]
9. Zhu Z, Liu LF, Su CF, et al. Corynoxine B derivative CB6 prevents Parkinsonian toxicity in mice by inducing PIK3C3 complex-dependent autophagy. Acta Pharmacol Sin 2022;43:2511-26. [Crossref] [PubMed]
10. Bisaglia M, Bubacco L. Copper Ions and Parkinson's Disease: Why Is Homeostasis So Relevant? Biomolecules 2020;10:195. [Crossref] [PubMed]
11. Filograna R, Beltramini M, Bubacco L, et al. Anti-Oxidants in Parkinson's Disease Therapy: A Critical Point of View. Curr Neuropharmacol 2016;14:260-71. [Crossref] [PubMed]
12. Dodani SC, Domaille DW, Nam CI, et al. Calcium-dependent copper redistributions in neuronal cells revealed by a fluorescent copper sensor and X-ray fluorescence microscopy. Proc Natl Acad Sci U S A 2011;108:5980-5. [Crossref] [PubMed]
13. Tsvetkov P, Coy S, Petrova B, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science 2022;375:1254-61. [Crossref] [PubMed]
14. Xing N, Dong Z, Wu Q, et al. Identification and validation of key molecules associated with humoral immune modulation in Parkinson's disease based on bioinformatics. Front Immunol 2022;13:948615. [Crossref] [PubMed]
15. Pan S, Song C, Meng H, et al. Identification of cuproptosis-related subtypes in lung adenocarcinoma and its potential significance. Front Pharmacol 2022;13:934722. [Crossref] [PubMed]
16. Zhu H, Mao F, Wang K, et al. Cuproptosis-related lncRNAs predict the clinical outcome and immune characteristics of hepatocellular carcinoma. Front Genet 2022;13:972212. [Crossref] [PubMed]
17. Wu J, Chen ZJ, Liang J, et al. Identifying and validating key genes mediating intracranial aneurysm rupture using weighted correlation network analysis and exploration of personalized treatment. Ann Transl Med 2022;10:1057. [Crossref] [PubMed]
18. Furió-Tarí P, Tarazona S, Gabaldón T, et al. spongeScan: A web for detecting microRNA binding elements in lncRNA sequences. Nucleic Acids Res 2016;44:W176-80. [Crossref] [PubMed]
19. Wu J, Li XY, Liang J, et al. Network pharmacological analysis of active components of Xiaoliu decoction in the treatment of glioblastoma multiforme. Front Genet 2022;13:940462. [Crossref] [PubMed]
20. Wu J, Li XY, Liang J, et al. Network pharmacological analysis of active components of Tongqiao Huoxue Decoction in the treatment of intracerebral hemorrhage. Ann Transl Med 2022;10:567. [Crossref] [PubMed]
21. Jin Z, Pei S, Ouyang L, et al. Thy-Wise: An interpretable machine learning model for the evaluation of thyroid nodules. Int J Cancer 2022;151:2229-43. [Crossref] [PubMed]
22. Double KL, Reyes S, Werry EL, et al. Selective cell death in neurodegeneration: why are some neurons spared in vulnerable regions? Prog Neurobiol 2010;92:316-29. [Crossref] [PubMed]
23. Miotto MC, Rodriguez EE, Valiente-Gabioud AA, et al. Site-specific copper-catalyzed oxidation of α-synuclein: tightening the link between metal binding and protein oxidative damage in Parkinson's disease. Inorg Chem 2014;53:4350-8. [Crossref] [PubMed]
24. Rose F, Hodak M, Bernholc J. Mechanism of copper(II)-induced misfolding of Parkinson's disease protein. Sci Rep 2011;1:11. [Crossref] [PubMed]
25. Kondo M, Hara H, Kamijo F, et al. 6-Hydroxydopamine disrupts cellular copper homeostasis in human neuroblastoma SH-SY5Y cells. Metallomics 2021;13:mfab041. [Crossref] [PubMed]
26. Ng L, Wang X, Yang C, et al. Celastrol Downmodulates Alpha-Synuclein-Specific T Cell Responses by Mediating Antigen Trafficking in Dendritic Cells. Front Immunol 2022;13:833515. [Crossref] [PubMed]
27. Zhang L, Zhang Y, Fan D. Pathological Role of Natural Killer Cells in Parkinson's Disease: A Systematic Review. Front Aging Neurosci 2022;14:890816. [Crossref] [PubMed]
28. Polishchuk R, Lutsenko S. Golgi in copper homeostasis: a view from the membrane trafficking field. Histochem Cell Biol 2013;140:285-95. [Crossref] [PubMed]
29. Jan A, Gonçalves NP, Vaegter CB, et al. The Prion-Like Spreading of Alpha-Synuclein in Parkinson's Disease: Update on Models and Hypotheses. Int J Mol Sci 2021;22:8338. [Crossref] [PubMed]
30. Thiruvengadam M, Venkidasamy B, Subramanian U, et al. Bioactive Compounds in Oxidative Stress-Mediated Diseases: Targeting the NRF2/ARE Signaling Pathway and Epigenetic Regulation. Antioxidants (Basel) 2021;10:1859. [Crossref] [PubMed]
31. Chen QX, Zhou L, Long T, et al. Galangin Exhibits Neuroprotective Effects in 6-OHDA-Induced Models of Parkinson's Disease via the Nrf2/Keap1 Pathway. Pharmaceuticals (Basel) 2022;15:1014. [Crossref] [PubMed]
32. Saini N, Georgiev O, Schaffner W. The parkin mutant phenotype in the fly is largely rescued by metal-responsive transcription factor (MTF-1). Mol Cell Biol 2011;31:2151-61. [Crossref] [PubMed]
33. Postuma RB, Berg D, Stern M, et al. MDS clinical diagnostic criteria for Parkinson's disease. Mov Disord 2015;30:1591-601. [Crossref] [PubMed]
34. Suwijn SR, van Boheemen CJ, de Haan RJ, et al. The diagnostic accuracy of dopamine transporter SPECT imaging to detect nigrostriatal cell loss in patients with Parkinson's disease or clinically uncertain parkinsonism: a systematic review. EJNMMI Res 2015;5:12. [Crossref] [PubMed]
35. Fereshtehnejad SM, Postuma RB. Subtypes of Parkinson's Disease: What Do They Tell Us About Disease Progression? Curr Neurol Neurosci Rep 2017;17:34. [Crossref] [PubMed]
36. De Pablo-Fernández E, Lees AJ, Holton JL, et al. Prognosis and Neuropathologic Correlation of Clinical Subtypes of Parkinson Disease. JAMA Neurol 2019;76:470-9. [Crossref] [PubMed]
37. Turcano P, Mielke MM, Bower JH, et al. Levodopa-induced dyskinesia in Parkinson disease: A population-based cohort study. Neurology 2018;91:e2238-43. [Crossref] [PubMed]
38. Ebersbach G, Ip CW, Klebe S, et al. Management of delirium in Parkinson's disease. J Neural Transm (Vienna) 2019;126:905-12. [Crossref] [PubMed]
39. Bertholo AP, França C, Fiorini WS, et al. Medical management after subthalamic stimulation in Parkinson's disease: a phenotype perspective. Arq Neuropsiquiatr 2020;78:230-7. [Crossref] [PubMed]
40. Kleiner-Fisman G, Herzog J, Fisman DN, et al. Subthalamic nucleus deep brain stimulation: summary and meta-analysis of outcomes. Mov Disord 2006;21:S290-304. [Crossref] [PubMed]
41. González-Burgos E, Carretero ME, Gómez-Serranillos MP. Nrf2-dependent neuroprotective activity of diterpenoids isolated from Sideritis spp. J Ethnopharmacol 2013;147:645-52. [Crossref] [PubMed]
42. Campolo M, Casili G, Biundo F, et al. The Neuroprotective Effect of Dimethyl Fumarate in an MPTP-Mouse Model of Parkinson's Disease: Involvement of Reactive Oxygen Species/Nuclear Factor-κB/Nuclear Transcription Factor Related to NF-E2. Antioxid Redox Signal 2017;27:453-71. [Crossref] [PubMed]

(English Language Editor: J. Jones)