Highlights
-
• Multi-algorithm machine learning prioritizes vesicle transport genes TRAPPC13 and COPS5 for Parkinson’s disease diagnosis.
-
• A two-gene model shows consistent accuracy across three GEO cohorts with good calibration and net benefit.
-
• Non-negative matrix factorization reveals two Parkinson’s disease subtypes linked to vesicle trafficking and immune signatures.
Introduction
Parkinson’s disease (PD) is a common age-related neurodegenerative disorder characterized by progressive loss of dopaminergic neurons in the substantia nigra and widespread presence of Lewy bodies containing α-synuclein aggregates. Reference Ng and Cao1 Clinically, PD presents with heterogeneous motor symptoms (resting tremor, rigidity, bradykinesia, postural instability) and a range of non-motor features such as cognitive impairment and autonomic dysfunction. Reference Irmady, Hale and Qadri2,Reference Su, Hou and Xu3 The incidence of PD increases sharply with age, affecting ∼ 1–3% of the population over 65. Reference Kong, Li and Ma4 Despite its growing prevalence and substantial socio-economic burden, the precise etiopathogenesis of PD remains incompletely understood, and there are currently no disease-modifying therapies. Reference Wei, Zhang and Zuo5,Reference Cai, Tang and Liu6 Diagnosis still relies primarily on clinical evaluation, as affordable and minimally invasive biomarkers with high sensitivity and specificity are lacking. Reference Munhoz, Tumas and Pedroso7 This gap highlights an urgent need for improved molecular understanding and novel biomarkers to facilitate earlier diagnosis and targeted interventions in PD. Reference Zarkali, Thomas and Zetterberg8
Converging evidence suggests that synaptic and vesicle trafficking dysfunction plays a central role in PD pathogenesis. Notably, multiple PD-causing gene mutations (e.g., SNCA, LRRK2, PARK2/Parkin, SYNJ1, DNAJC6) implicate disruptions in synaptic vesicle recycling and endocytic pathways. Reference Bu, Follett and Deng9,Reference Song, Yuan and Cui10 Impaired synaptic endocytosis and vesicle trafficking are increasingly recognized as early pathological events that precede overt neurodegeneration. Reference Vidyadhara, Bäckström and Chakraborty11 For example, a recent neuron model of PD showed that a Parkin (PARK2) mutation disrupts the recycling of synaptic vesicles, leading to dopamine release deficits before neuronal death occurs. Reference Song, Peng and Sauve12 Likewise, loss-of-function mutations in vesicle-associated proteins such as auxilin (DNAJC6) or synaptojanin-1 (SYNJ1) cause early-onset parkinsonism marked by synaptic vesicle endocytosis defects. Reference Köroğlu, Baysal and Cetinkaya13,Reference Lesage, Mangone and Tesson14 Even idiopathic PD brains exhibit evidence of presynaptic dysfunction – for instance, Golgi apparatus fragmentation and vesicular transport abnormalities have been observed, linking PD to broader intracellular trafficking deficits. Reference Mohan, Calenic and Ghiurau15,Reference Tong and Shen16 Moreover, genome-wide studies have begun identifying novel PD risk genes related to vesicle dynamics; a notable example is SV2C, a synaptic vesicle glycoprotein whose variants were highlighted by a recent genome-wide association study as a PD risk factor affecting dopamine neuron vesicle function. Reference Chang, Lim and Foo17 Together, these findings reinforce the hypothesis that vesicle-mediated transport and synaptic organization are critically perturbed in PD, motivating a focused investigation of vesicle transport-related genes in the disease.
In this study, we employed an integrative bioinformatics strategy to elucidate the dysregulation of vesicle-mediated transport-related genes (VMTRGs) in PD and to identify novel diagnostic biomarkers. Rather than using a hypothesis-agnostic feature space, we specifically centered our analysis on VMTRGs, thereby anchoring biomarker discovery to a vesicle-mediated transport–focused biological hypothesis. By combining differential expression analysis, co-expression network construction and a multi-algorithm machine learning framework (least absolute shrinkage and selection operator [LASSO], random forest [RF] and extreme gradient boosting [XGBoost]), we aimed to uncover key molecular drivers linking intracellular trafficking defects to PD pathology. Only genes consistently prioritized across algorithms were retained, enhancing robustness compared with single-model approaches. Our analyses revealed two pivotal genes, TRAPPC13 and COPS5, representing components of the Transport Protein Particle (TRAPP) vesicle tethering complex and the COP9 signalosome (CSN), respectively – machineries essential for vesicle transport and protein homeostasis. Notably, TRAPPC13 has not been prominently reported as a PD-related biomarker in prior transcriptome-based machine learning studies, highlighting its potential novelty in this disease context. We further validated their diagnostic utility across independent cohorts and demonstrated that their expression profiles could stratify PD into distinct molecular subtypes. Collectively, this study provides new insights into PD pathogenesis and offers potential biomarkers and therapeutic targets for precision diagnosis and treatment.
Materials and Methods
Data collection
Transcriptomic datasets of PD were retrieved from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). The GSE20141 dataset, comprising 18 samples (10 PD and 8 healthy controls), was designated as the training cohort. Two independent datasets – GSE20163 (17 samples; 9 PD and 8 controls) and GSE7621 (25 samples; 16 PD and 9 controls) – were used as external validation cohorts. All datasets were based on microarray transcriptome profiling. In addition, a list of 724 VMTRGs was obtained from the Reactome pathway REACTOME_VESICLE_MEDIATED_TRANSPORT. These genes were used as the reference set for subsequent analyses.
Differential expression analysis of vesicle transport-related genes
Differential expression analysis between PD and control samples in the training dataset (GSE20141) was performed using the “limma” package in R, with significance thresholds set at p < 0.05 and |log2fold change| >0.585. This threshold preserved more biologically relevant genes when the sample size was limited. The resulting differentially expressed genes (DEGs) were then intersected with the VMTRGs obtained from the Reactome database to identify differentially expressed vesicle-mediated transport-related genes (DEVMTRGs). Gene Ontology (GO) enrichment analysis of these DEVMTRGs was conducted using the “clusterProfiler” package.
Weighted gene co-expression network analysis (WGCNA) and identification of key modules
WGCNA was performed based on the expression profiles of DEGs identified in the training dataset (GSE20141) using the “WGCNA” package in R. This DEG-based input strategy was used to enhance network stability and reduce noise in the setting of a limited sample size. VMTRGs were incorporated downstream by intersecting PD-associated modules with the curated VMTRG set for subsequent prioritization and interpretation. The optimal soft-thresholding power (β = 7) was selected to ensure scale-free topology of the network. A hierarchical clustering tree was then constructed to identify distinct gene co-expression modules, which were color-coded for visualization. Pearson correlation analysis was used to evaluate the relationships between module eigengenes and clinical traits (PD vs. control) to identify modules most strongly associated with PD. Genes within these key modules were subsequently subjected to GO enrichment analysis to explore their potential biological functions.
Machine learning-based identification of key feature genes
To further identify disease-associated feature genes, an integrative machine learning approach was applied. The intersection between the DEVMTRGs and the genes from the key WGCNA modules was first obtained to generate a candidate gene set. Subsequently, three machine learning algorithms – LASSO regression, RF and XGBoost – were employed to identify robust diagnostic feature genes. These analyses were implemented using the R packages “glmnet,” “randomForest” and “xgboost,” respectively. LASSO regression employed binary logistic regression (family = “binomial”) to generate 1000 lambda (λ) sequences (nlambda = 1000, alpha = 1). It was optimized via 10-fold cross-validation, with the optimal number of features determined by the minimum cross-validation error (lambda.min = 0.001427). RF employed a two-step strategy incorporating 10-fold cross-validation: By constructing ntree = 800 decision tree models, the optimal number of trees (ntree = 23) was determined via Out-of-Bag error rate. Key feature genes were filtered based on MeanDecreaseGini > 0.5. XGBoost employed a binary logistic objective function (max_depth = 6, eta = 0.1, subsample = 0.8, colsample_bytree = 0.8). Optimal iteration count was determined via 5-fold cross-validation (nfold = 5, nrounds = 500), with feature genes contributing to classification selected based on Gain importance metrics. All analyses were performed with a random seed to ensure reproducible results. The overlapping genes identified by all three algorithms were considered key feature genes with potential diagnostic significance. Given the small training cohort, we prioritized parsimonious feature selection by applying multistep dimensionality reduction and retaining only genes consistently identified across three machine learning algorithms, with subsequent evaluation in independent external datasets to reduce overfitting risk.
Receiver operating characteristic (ROC) analysis and diagnostic model construction
ROC curve analysis was performed in the training dataset (GSE20141) to evaluate the diagnostic performance of the feature genes identified through machine learning. The analysis was conducted using the R package “pROC,” and the area under the curve (AUC) was used as a measure of diagnostic accuracy. A diagnostic model was subsequently constructed based on the selected feature genes. A nomogram was generated using the “rms” and “rmda” packages to visualize the model, while decision curve analysis (DCA) and calibration curves were used to assess its clinical utility and calibration performance. The model was further validated in two independent testing datasets (GSE20163 and GSE7621). In addition, expression levels of the diagnostic genes were compared between the PD and control groups using boxplots.
Functional enrichment analysis, regulatory network construction and drug prediction
Single-gene enrichment analysis was performed for each diagnostic feature gene to explore its potential biological functions and associated pathways. Upstream transcription factors (TFs) regulating the diagnostic genes were predicted using NetworkAnalyst (https://www.networkanalyst.ca/). Experimentally validated miRNA–mRNA interactions were retrieved from miRTarBase (https://mirtarbase.cuhk.edu.cn/) to identify potential regulatory miRNAs targeting the diagnostic genes. Based on these miRNAs, corresponding lncRNA–miRNA interactions were further obtained to construct lncRNA–miRNA–mRNA competing endogenous RNA (ceRNA) networks. All regulatory networks were visualized using Cytoscape software. In addition, potential small-molecule drugs targeting the diagnostic genes were predicted using the DSigDB database (https://dsigdb.tanlab.org/DSigDBv1.0/). Specifically, the diagnostic genes were uploaded to the Enrichr platform (https://maayanlab.cloud/modEnrichr/) to perform enrichment analysis against DSigDB, and compounds with significant enrichment scores were considered potential therapeutic candidates.
Immune infiltration analysis
Immune cell infiltration in the training cohort (GSE20141) was evaluated using the single-sample gene set enrichment analysis (ssGSEA: MSigDB C7 immune signature gene sets, v2023.1) algorithm to estimate immune cell populations and immune-related functional scores. In addition, the CIBERSORT algorithm (Standard LM22 Reference Matrix: 22 Immune Cell Types) was applied to further quantify the relative proportions of 22 immune cell types. Spearman correlation analysis was then performed to assess the associations between the diagnostic genes and infiltrating immune cell populations. The reference immune gene sets used for ssGSEA (MSigDB C7 immunologic signature collection, v2023.1) are provided in Supplementary Table S1A. The LM22 signature matrix used by CIBERSORT is provided in Supplementary Table S1B.
Molecular subtype identification based on diagnostic genes
To identify potential molecular subtypes of PD, PD samples from the three datasets (GSE20141, GSE20163 and GSE7621) were merged into a single expression matrix. Batch effects across datasets were corrected using the “removeBatchEffect” function in the “limma” package. Non-negative matrix factorization (NMF) clustering was then performed based on the expression profiles of the diagnostic feature genes. The number of clusters k was tested from 2 to 10, and the optimal clustering number (k = 2) was determined using the cophenetic coefficient, Silhouette score and dispersion coefficient. PD samples were thus classified into distinct molecular subtypes. Principal component analysis (PCA) was used to visualize the clustering results. Boxplots were generated to compare the expression levels of diagnostic genes among the identified subtypes. The ssGSEA algorithm was applied to assess differences in immune cell infiltration and immune-related functional scores across subtypes. In addition, the MCP-counter algorithm was used to quantify immune cell abundance for each subtype. Differential expression analysis between subtypes (p < 0.05, |log2FC| >0.585) was conducted using “limma,” and the resulting genes were subjected to GO and KEGG pathway enrichment analyses using the “clusterProfiler” package.
Statistical analysis
All statistical analyses were conducted using R software (version 4.2.2) and its associated packages. Differential expression between groups was assessed using the Wilcoxon rank-sum test or Student’s t-test as appropriate. Pearson or Spearman correlation analysis was applied to evaluate associations between gene expression and clinical or immune features. ROC curves and AUC values were calculated using the “pROC” package to assess diagnostic performance. Survival and classification accuracy were evaluated with nomogram calibration curves and DCA. All statistical tests were two-sided, and a p-value < 0.05 was considered statistically significant unless otherwise specified.
Results
Differential expression analysis of VMTRGs in PD
Differential expression analysis between PD and control samples in the training dataset (GSE20141) identified numerous DEGs (Figure 1A). By intersecting these DEGs with the VMTRGs from Reactome, 34 DEVMTRGs were obtained (Figure 1B). A heatmap illustrated distinct expression patterns of these 34 DEVMTRGs between the PD and control groups, indicating that vesicle transport-related transcriptional dysregulation is a prominent feature of PD (Figure 1C). GO enrichment analysis revealed that these genes were primarily associated with vesicle organization, Golgi-associated vesicle membranes, TRAPP and vesicle tethering complexes and membrane fission, suggesting their involvement in intracellular trafficking and synaptic vesicle dynamics (Figure 1D).
Identification of differentially expressed vesicle-mediated transport-related genes (DEVMTRGs) in Parkinson’s disease (PD). (A) Volcano plot showing differentially expressed genes (DEGs) between PD and control samples (|log2FC| > 0.585, p < 0.05). (B) Venn diagram illustrating the overlap between DEGs and VMTRGs from Reactome, yielding 34 intersecting genes. (C) Heatmap displaying the expression profiles of the 34 DEVMTRGs across the PD and control groups. (D) Gene Ontology (GO) enrichment analysis of DEVMTRGs, categorized into biological process (BP), cellular component (CC) and molecular function (MF) terms.
WGCNA analysis and identification of key modules
WGCNA was performed using the DEGs from the training dataset to identify gene modules associated with PD. The optimal soft-thresholding power was determined to be 7, ensuring a scale-free topology of the network (Figure 2A). Hierarchical clustering of genes resulted in multiple co-expression modules, each assigned a unique color for visualization (Figure 2B). Correlation analysis between module eigengenes and clinical traits revealed that the green-yellow module showed a strong positive correlation with PD (r = 0.58, p = 0.01), while the blue module was negatively correlated (r = –0.47, p = 0.05) (Figure 2C–E). These two modules were therefore considered PD-related key modules, and the genes within them were defined as hub genes for subsequent analyses. Functional enrichment analysis of the hub genes using “clusterProfiler” demonstrated significant enrichment in biological processes related to vesicle organization, Golgi vesicle transport, endoplasmic reticulum organization and protein localization to the membrane, indicating their involvement in intracellular trafficking and membrane dynamics (Figure 2F).
Weighted gene co-expression network analysis and identification of key Parkinson’s disease (PD)-associated modules. (A) Analysis of scale-free topology fit index and mean connectivity to determine the optimal soft-thresholding power (β = 7). (B) Hierarchical clustering dendrogram of genes, with each color representing a distinct co-expression module. (C) Module–trait relationships showing correlations between module eigengenes and PD status. (D–E) Scatter plots showing the relationships between module membership and gene significance in the green-yellow (positively correlated) (D) and blue (negatively correlated) modules (E). (F) Gene Ontology enrichment analysis of hub genes from PD-related modules, highlighting vesicle transport- and membrane organization-related processes.
Machine learning-based identification of key feature genes
To identify robust diagnostic biomarkers associated with PD, an integrative machine learning approach was applied. The intersection between the DEVMTRGs and the hub genes from key WGCNA modules yielded a set of 12 candidate genes (Figure 3A). Three machine learning algorithms – LASSO regression, RF and XGBoost – were then employed to further refine these candidates and identify the most informative diagnostic features. The LASSO model was tuned based on the optimal λ value (Figure 3B). For RF, model error stabilized with increasing trees (Figure 3C), and feature ranking was quantified using MeanDecreaseAccuracy and MeanDecreaseGini (Figure 3D). XGBoost feature importance was derived from the model’s importance metric (Gain). Genes consistently identified across all three machine learning algorithms were considered robust biomarkers. Two overlapping genes, TRAPPC13 and COPS5, were ultimately selected as key diagnostic feature genes for subsequent analyses (Figure 3E, F).
Identification of key diagnostic feature genes using multiple machine learning algorithms. (A) Venn diagram showing the intersection of differentially expressed vesicle-mediated transport-related genes (DEVMTRGs) and hub genes derived from weighted gene co-expression network analysis (WGCNA) modules. (B) Least absolute shrinkage and selection operator (LASSO) regression analysis and cross-validation curve for selecting the optimal regularization parameter (λ). (C) Error rate of the random forest (RF) model as a function of the number of trees. (D) Importance ranking of the top variables identified by the RF algorithm based on MeanDecreaseAccuracy and MeanDecreaseGini values. (E) Relative feature importance of candidate genes identified by the XGBoost algorithm. (F) Intersection analysis of the genes identified by LASSO, RF and XGBoost models.
ROC curve analysis and diagnostic model validation
ROC curve analysis was performed in the training dataset (GSE20141) to evaluate the diagnostic performance of the two machine learning-identified feature genes, TRAPPC13 and COPS5. Both genes demonstrated strong discriminative ability between PD and control samples (Figure 4A). A diagnostic model incorporating these two genes was subsequently constructed and validated across both the training and independent testing datasets (GSE20163 and GSE7621), showing consistent diagnostic performance (Figure 4B). A nomogram integrating TRAPPC13 and COPS5 expression was developed to facilitate individual risk prediction (Figure 4C). The model’s potential decision-analytic utility (net benefit) and calibration were assessed using DCA and calibration plots, respectively, both indicating satisfactory predictive accuracy (Figure 4D,E). These results support the model’s performance in the analyzed brain-tissue cohorts; its role as an auxiliary tool and its clinical translatability require further validation. Boxplots further confirmed that the expression levels of TRAPPC13 and COPS5 were significantly lower in PD samples compared with controls in the training dataset (GSE20141) (Figure 4F).
Diagnostic performance and model validation of TRAPPC13 and COPS5 in Parkinson’s disease (PD). (A) Receiver operating characteristic (ROC) curves of TRAPPC13 and COPS5 in the training dataset (GSE20141). (B) ROC curves of the two-gene diagnostic model validated in GSE20141, GSE20163 and GSE7621 datasets. (C) Nomogram based on TRAPPC13 and COPS5 expression for individualized diagnostic prediction. (D) Decision curve analysis evaluating the clinical utility of the diagnostic model. (E) Calibration plot assessing the model’s predictive performance. (F) Boxplots showing differential expression of TRAPPC13 and COPS5 between the PD and control groups.
Functional enrichment, regulatory network construction and drug prediction
Single-gene enrichment analysis was conducted for the two diagnostic biomarkers, COPS5 and TRAPPC13, to investigate their potential biological functions. The enrichment plots indicated that both genes were mainly involved in protein localization and post-translational protein modification, highlighting their roles in intracellular trafficking (Figure 5A). To explore potential upstream regulatory mechanisms, a ceRNA network was constructed based on experimentally validated interactions retrieved from miRTarBase. The resulting network included multiple miRNAs (e.g., hsa-miR-26b-5p, hsa-miR-16-5p) and corresponding lncRNAs (e.g., UCA1, PCOTH, LINC00632), suggesting multilayered post-transcriptional regulation of COPS5 and TRAPPC13 (Figure 5B). Additionally, TF-gene regulatory interactions predicted by NetworkAnalyst identified several key TFs potentially regulating these diagnostic genes, including TP53, STAT1, GATA2 and JUN, indicating coordinated transcriptional control (Figure 5C). Potential small-molecule compounds targeting COPS5 and TRAPPC13 were further screened using the DSigDB database via the Enrichr platform, and the top five candidates are summarized in Table 1. These drug candidates were generated in silico and should be interpreted as hypothesis-generating, warranting further experimental validation.
Functional enrichment and regulatory network analysis of diagnostic genes COPS5 and TRAPPC13. (A) Single-gene gene set enrichment analysis plots showing biological processes associated with COPS5 and TRAPPC13. (B) Competing endogenous RNA regulatory network (lncRNA–miRNA–mRNA) constructed based on experimentally validated interactions from the miRTarBase database. (C) TF-gene regulatory network of COPS5 and TRAPPC13 predicted by NetworkAnalyst.
Predicted candidate drugs targeting the diagnostic genes
Immune infiltration analysis
The immune landscape of PD was assessed using multiple computational approaches. The ssGSEA algorithm was used to evaluate immune cell infiltration and immune-related functional activities between the PD and control groups, and significant differences were observed in several specific immune features, including CD8+ T-cell fraction, cytolytic activity and major histocompatibility complex (MHC) class I–related signaling (Figure 6A,B). To further characterize the immune microenvironment, the CIBERSORT algorithm was applied to estimate the relative abundance of 22 immune cell subsets. Immune populations, including plasma cells and resting mast cells, exhibited distinct infiltration patterns between groups (Figure 6C). Spearman correlation analysis revealed strong associations between the diagnostic biomarkers (TRAPPC13 and COPS5) and immune infiltration features, particularly with antigen presentation-related pathways and lymphocyte populations (Figure 6D). Specifically, TRAPPC13 expression was negatively correlated with MHC class I infiltration, while COPS5 expression showed an inverse relationship with B-cell abundance (Figure 6E,F). These findings suggest that vesicle transport-related biomarkers may influence PD pathogenesis through immune regulatory mechanisms.
Immune infiltration analysis in Parkinson’s disease (PD). (A) Single-sample gene set enrichment analysis (ssGSEA)-based comparison of immune cell infiltration scores between PD and control samples. (B) ssGSEA-derived immune function scores in the PD and control groups. (C) Immune cell composition estimated by the CIBERSORT algorithm. (D) Correlation heatmap showing associations between diagnostic genes (TRAPPC13 and COPS5) and immune-related features. (E–F) Scatter plots displaying negative correlations between TRAPPC13 and MHC class I infiltration (E) and between COPS5 and B-cell infiltration (F).
Molecular subtype identification based on diagnostic genes
To explore potential molecular heterogeneity among PD samples, three datasets (GSE20141, GSE20163 and GSE7621) were integrated after batch effect correction using the removeBatchEffect function in “limma.” PCA demonstrated that batch-related variation was effectively eliminated (Figure 7A). Based on the expression profiles of the two diagnostic genes (TRAPPC13 and COPS5), NMF clustering identified two stable molecular subtypes (k = 2) of PD (Figure 7B,C). When k = 2, the cophenetic coefficient reached its maximum, the silhouette score was optimal, the dispersion coefficient decreased significantly and the consensus matrix exhibited a clearly defined block structure. PCA visualization confirmed clear separation between the two subtypes (Figure 7D). Expression analysis revealed that COPS5 expression was significantly higher in subtype 1, whereas TRAPPC13 showed no significant difference (Figure 7E). Immune infiltration analysis by ssGSEA indicated significant differences in type II IFN response score between the two subtypes (Figure 7F). The MCP-counter algorithm further showed distinct proportions of immune cell lineages, particularly endothelial cells and monocytic lineage (Figure 7G). Differential expression analysis between the two subtypes identified genes mainly enriched in immune regulation pathways, including regulation of lymphocyte proliferation, negative regulation of acute inflammatory response and cytokine–cytokine receptor interaction (Figure 7H,I). These findings suggest that the two PD subtypes exhibit distinct immune microenvironmental characteristics and regulatory signaling profiles.
Molecular subtype identification based on diagnostic genes in Parkinson’s disease (PD). (A) Principal component analysis (PCA) plots showing sample distribution before and after batch correction across GSE20141, GSE20163 and GSE7621 datasets. (B–C) Non-negative matrix factorization (NMF) clustering results determining the optimal rank (k = 2) and the corresponding consensus matrix. (D) PCA plot confirming distinct separation between the two identified molecular subtypes. (E) Boxplots of TRAPPC13 and COPS5 expression across subtypes. (F) Single-sample gene set enrichment analysis-based comparison of immune cell infiltration and immune function scores between PD subtypes. (G) Immune cell abundance estimated by MCP-counter across subtypes. (H) Gene Ontology (GO) enrichment analysis of subtype-related differentially expressed genes. (I) KEGG pathway enrichment highlighting cytokine–cytokine receptor interaction and neuroactive ligand–receptor interaction.
Discussion
In this study, we employed an integrated bioinformatics approach to prioritize key vesicle transport-related biomarkers and transcriptomic signatures in PD. Through differential expression analysis of PD transcriptomic data and subsequent network and machine learning prioritization, we identified TRAPPC13 and COPS5 as two robust genes with diagnostic relevance in PD and constructed a parsimonious two-gene model that was evaluated in the training set and independent validation datasets. Our findings highlight the dysregulation of intracellular vesicle trafficking and protein homeostasis pathways in PD and suggest a potential association between these pathways and the neuroinflammatory environment of the disease. Here, we discuss the implications of these results in the context of existing literature, the potential roles of TRAPPC13 and COPS5 in PD pathogenesis, the performance and generalizability of our PD diagnostic model and the broader biological insights gained from transcriptome-based computational analyses.
Our analysis confirmed that genes associated with vesicle trafficking are significantly altered in PD. We identified 34 VMTRGs with differential expression in PD versus control samples, underscoring that vesicle transport processes are perturbed in the diseased nigrostriatal system. Enrichment analysis of these genes revealed overrepresentation of biological processes such as synaptic vesicle organization, Golgi-associated vesicle transport, vesicle tethering complexes and membrane fission – functions critical for neurotransmitter release and intracellular trafficking. Reference Wu, Zheng and Liu18,Reference Brooker, Naylor and Krainc19 These results align with accumulating evidence that synaptic and endomembrane system defects are an early and central feature of PD. Reference Imbriani, Schirinzi and Meringolo20,Reference Vidyadhara, Lee and Chandra21 Impaired synaptic vesicle endocytosis has been shown to lead to dopamine handling deficits and axonal degeneration prior to neuron loss in PD models. Reference Vidyadhara, Somayaji and Wade22 For example, α-synuclein pathology can block endoplasmic reticulum (ER)-to-Golgi vesicular transport and disrupt Golgi organization, contributing to dopaminergic neuron stress. Reference Coune, Bensadoun and Aebischer23 Likewise, inactivation of the PD gene LRRK2 causes Golgi fragmentation and misrouting of vesicles, Reference Mohan, Calenic and Ghiurau15 while pathogenic mutations in Parkin (PARK2) were recently found to hinder synaptic vesicle recycling in human neurons, leading to synaptic dysfunction before cell death. Reference Stern, Lau and Manole24 Our findings provide transcriptome-level support that vesicle transport gene dysregulation is a robust feature of PD, consistent with a state of synaptic dysfunction and organelle trafficking stress in affected neurons. This is compatible with the concept of a “dying-back” neuropathological process in PD, wherein synaptic and axonal disturbances (e.g., failed vesicle recycling, organelle accumulation) may precede and contribute to the eventual degeneration of cell bodies. Reference Burke and O’Malley25 Targeting these early vesicle trafficking deficits – for instance, via therapies that bolster endocytic recycling or Golgi function – may represent potential therapeutic entry points to slow neurodegeneration in PD; Reference Jennings, Huntwork-Rodriguez and Vissers26 however, this therapeutic implication remains hypothesis-generating and requires direct experimental validation.
Through WGCNA and multi-algorithm machine learning, we prioritized TRAPPC13 and COPS5 as the most consistently PD-associated genes among the vesicle transport candidates. Both genes were downregulated in PD samples and together yielded a diagnostic model with good and reproducible performance (AUC >0.9 in training and ∼ 0.8–0.9 in validation cohorts). Notably, TRAPPC13 and COPS5 were selected by three distinct feature selection methods (LASSO, RF, XGBoost) – supporting their robust signal amid the complex gene expression changes in PD. Reference Wu, Qin and Cai27 To our knowledge, neither TRAPPC13 nor COPS5 has been widely reported as a PD biomarker in prior studies, making them novel candidates emerging from our unbiased approach. To contextualize these candidates, we queried disease–gene resources (GeneCards and text-mined databases such as DISEASES) and performed targeted literature mining. TRAPPC13 is not prominently annotated to neurodegenerative disorders in these resources, suggesting that its implication in PD in our study may represent a relatively new disease context. In contrast, COPS5 (JAB1/CSN5) has prior experimental links to neurodegenerative biology, including evidence that oligodendroglial loss of JAB1/CSN5 induces inflammation and neurodegeneration Reference Rivellini, Porrello and Dina28 and that COPS5 can promote amyloidogenic APP processing/Aβ generation. Reference Wang, Dey and Carrera29 Together, these external data support the biological plausibility of COPS5 while underscoring the relative scarcity of TRAPPC13 neurodegenerative annotations.
Across training and external validation cohorts, this parsimonious two-gene signature showed reproducible discrimination with acceptable calibration and potential net benefit in DCA; the nomogram supports individualized probability estimation. Given its brain-tissue derivation and retrospective evaluation, it should be viewed as an auxiliary, hypothesis-generating tool for research stratification rather than a stand-alone diagnostic test. Prospective studies are needed to test peripheral-biospecimen generalizability, early-stage performance and incremental value beyond existing clinical variables/biomarkers.
In our study, TRAPPC13 was consistently downregulated in PD samples and was repeatedly prioritized by the machine learning pipeline, supporting its robustness as a transcriptomic candidate biomarker. TRAPPC13 encodes a subunit of the TRAPP complex, specifically part of the mammalian TRAPPIII complex involved in tethering vesicles for ER-to-Golgi transport and endosome-Golgi trafficking. Reference Ramírez-Peinado, Ignashkova and van Raam30,Reference Hall, Sawant and Gu31 TRAPPC13 has also been implicated in autophagy: it participates in a TRAPPIII-dependent mechanism required for autophagic flux under stress conditions. Reference Papaioannou, Wallace and Malhotra32 Consistent with our differential-expression results, TRAPPC13 was reduced in PD; the following mechanistic interpretations are derived from prior studies and should be viewed as hypotheses. In addition, TRAPPC13 depletion has been shown to reduce Rab1 activity, cause ER stress and impair autophagy in cellular models. Reference Ramírez-Peinado, Ignashkova and van Raam30 In neurons, such effects could lead to accumulation of misfolded proteins and defective organelles – phenomena observed in PD pathology (e.g., aggregation of misfolded α-synuclein and dysfunctional mitochondria due to autophagy failure). Reference Dai, Liu, Ke, Chen, Fang and Zhang33 Although α-synuclein burden and mitochondrial function were not directly measured in our transcriptomic datasets, TRAPPC13 downregulation may plausibly link vesicle trafficking/autophagic impairment to these canonical PD hallmarks. Reduced TRAPPIII-related trafficking and autophagy could compromise the clearance of aggregation-prone proteins (including α-synuclein) and weaken mitochondrial quality control (e.g., mitophagy), thereby favoring protein/organellar stress. Future studies integrating α-synuclein pathology and mitochondrial functional assays are warranted. Interestingly, TRAPPC13 is absent in yeast and appears to represent a human-specific adaptation of TRAPP with possibly specialized functions in certain tissues. Reference Sun and Sui34 Overall, our results are consistent with the growing recognition that ER-Golgi and endolysosomal traffic disruptions are central to PD. It raises the possibility that TRAPPC13 or its associated TRAPPIII complex could be explored as a candidate therapeutic entry point to enhance autophagy or protein trafficking in PD, although this therapeutic implication remains speculative and requires direct experimental testing.
In our study, COPS5 was consistently downregulated in PD tissues and was repeatedly prioritized by the machine learning pipeline, supporting its robustness as a transcriptomic candidate biomarker. COPS5 (also known as JAB1/CSN5) encodes the fifth subunit of the CSN, a conserved multiprotein complex that regulates ubiquitin-mediated protein degradation by de-neddylating cullin-RING ligases. Reference Faull, Lau and Martens35 COPS5 is the catalytic heart of the CSN complex and controls the turnover of numerous proteins involved in cell cycle, DNA repair and stress responses. Reference Schulze-Niemand and Naumann36 The mechanistic links discussed below are inferred from prior literature rather than directly demonstrated in the current transcriptomic study. Diminished COPS5/CSN function could have wide-ranging consequences in neurons: it may impair proteostasis (by dysregulating ubiquitin ligases), exacerbate oxidative DNA damage and promote cellular senescence – all processes linked to neurodegeneration. Reference Zhang, Yu and Li37,Reference Pan, Yang and Claret38 Remarkably, a recent study demonstrated that cell type–specific loss of COPS5 in oligodendrocytes induced DNA damage, cellular senescence, chronic inflammation and secondary neurodegeneration in mice. Reference Rivellini, Porrello and Dina28 Oligodendroglial deletion of COPS5 led to senescent glia that secreted inflammatory factors, activated microglia and ultimately caused demyelination and neuron loss. Reference Rivellini, Porrello and Dina28 This finding underscores COPS5’s role in maintaining central nervous system cellular homeostasis and suppressing pathological inflammation. Reference Tian, Milic and Monasor39 In the context of PD, reduced COPS5 in dopaminergic neurons or glia could potentially contribute to an environment of genomic instability and inflammation that may accelerate neuron death. Reference Trainor, MacDonald and Penney40,Reference Pfeifer41 Moreover, the CSN has been shown to exert anti-inflammatory effects in the brain – for instance, pharmacological CSN activation reduced neuroinflammation and neuronal stress in an ischemic stroke model. Reference Tian, Milic and Monasor39 A deficiency in COPS5 might therefore remove a brake on immune activation, aligning with our observation that COPS5 expression inversely correlated with markers of B-cell and antigen-presenting immune infiltrates in PD. Taken together, our results support COPS5 as a potential nexus between impaired protein clearance and neuroinflammation in PD. Restoring CSN/COPS5 activity might be explored as a strategy to simultaneously bolster degradation of misfolded proteins (like toxic α-synuclein) and temper inflammatory cascades in the PD brain, although this therapeutic implication remains speculative. In addition to shaping inflammatory tone, reduced COPS5/CSN activity may also intersect with classic PD pathology through impaired ubiquitin-dependent proteostasis. Dysregulation of cullin-RING ligases could limit the clearance of aggregation-prone proteins and amplify cellular stress responses that secondarily impair mitochondrial homeostasis. These possibilities warrant targeted validation using assays of α-synuclein clearance and mitochondrial quality control.
Beyond prior biological plausibility, our results suggest that TRAPPC13 and COPS5 may jointly link vesicle trafficking/proteostasis dysregulation to an altered immune milieu in PD. Both genes were downregulated in PD and correlated with immune infiltration and immune functional signatures, consistent with a more activated neuroinflammatory context. Moreover, NMF subtyping based on the two-gene signature delineated subgroups with distinct immune and pathway profiles, suggesting PD heterogeneity in trafficking/protein homeostasis stress coupled to immune activation. Because these results are inferred from bulk transcriptomic data and remain correlative, they should be interpreted as hypothesis-generating. We therefore propose a testable hypothesis that reduced TRAPPC13/COPS5-linked trafficking–proteostasis capacity may exacerbate protein/organellar stress and amplify inflammatory signaling in PD.
Notably, our immune deconvolution suggests that TRAPPC13/COPS5 expression correlates with immune signatures, but this remains associative and not causal. PD is increasingly recognized to involve immunopathology, including microglial activation and antigen presentation programs. Reference Tansey, Wallings and Houser42 Trafficking/autophagy–lysosomal dysfunction may couple proteostasis to immune activation via altered endolysosomal processing and extracellular vesicle signaling. Reference Guo, Wang and Zhao43 These links should be tested using single-cell/spatial profiling and TRAPPC13/COPS5 perturbation with readouts of α-synuclein handling, antigen presentation and inflammatory outputs.
To explore why these genes might be relevant in PD, we performed extensive functional enrichment and regulatory network analyses. GO results for TRAPPC13 and COPS5 consistently pointed to roles in protein localization, vesicle-mediated transport and post-translational protein modification, consistent with their involvement in trafficking and proteostasis. Reference Zhao, Li and Luo44 By constructing a ceRNA network, we identified several non-coding RNAs that may modulate TRAPPC13/COPS5 expression. Notably, our network included microRNAs such as miR-16-5p and miR-26b-5p, which are known to target genes in autophagy and inflammation pathways, and lncRNAs such as NEAT1 and UCA1, which have been previously implicated in PD or other neurological diseases. Reference He, Qiu and Zhang45,Reference Xu, Xie and Jing46 These connections hint that TRAPPC13 and COPS5 are embedded in post-transcriptional regulatory circuits that respond to cellular stress. For instance, NEAT1 is a stress-responsive lncRNA elevated in neurodegeneration and can sponge miR-124 to influence microglial activation. Reference Asadi, Abed and Kouchakali47 TFs regulating TRAPPC13 and COPS5 were inferred in silico using the NetworkAnalyst platform, revealing TP53, STAT1, GATA2 and JUN as candidate regulators. These TFs are biologically plausible in PD context: TP53 (p53) is activated by cellular stress and has been shown to mediate dopamine neuron apoptosis in PD models; Reference Hernández Borrero and El-Deiry48 STAT1 is induced by inflammatory cytokines and could link neuroinflammation to gene expression changes in PD; Reference Pei, Hou and Ma49 JUN (c-Jun) is an AP-1 factor activated in neurotoxic conditions, and GATA2 is expressed in midbrain dopamine neurons and involved in their development. Reference Crocker, Khan and Cameron50,Reference Oliveira, Balling and Smidt51 The enrichment of p53 and STAT/JUN signaling upstream of our genes is consistent with the hypothesis that oxidative DNA damage and pro-inflammatory signaling in PD activate these stress-related TFs, which then may modulate TRAPPC13/COPS5 expression. This hypothesis is conceptually consistent with our earlier reasoning that loss of TRAPPC13/COPS5 could exacerbate stress, creating a feed-forward loop. It is noteworthy that p53 can interact with Parkin and α-synuclein pathways, linking to PD at multiple levels. Reference Liu, Liu and Li52 Future experimental studies could examine whether p53 or cytokine exposure in neuronal cells leads to reduced TRAPPC13 or COPS5 expression, shedding light on this regulatory mechanism.
Regarding drug repurposing, our DSigDB/Enrichr results should be interpreted as hypothesis-generating. Among the predicted candidates, α-lipoic acid (dl-thioctic acid) has been reported to have neuroprotective effects in experimental PD models. Reference Li, Wang and Zhou53,Reference Kulikova, Berezhnoy and Stvolinsky54 The identification of sphingosine further points to potential involvement of sphingolipid/S1P signaling, as S1P-axis modulation (e.g., fingolimod/FTY720) has shown neuroprotective effects in PD mouse models, including evidence implicating mitochondrial quality control/mitophagy. Reference Zhao, Yang and Yang55,Reference Rajan, Sood and Jain56 For the remaining predicted compounds, direct neuroprotection-focused evidence in PD models is limited, and targeted experimental validation is required.
While this study offers novel insights, several limitations should be noted. First, the overall sample size was modest (n = 60 across three datasets), and the training cohort was particularly small (n = 18), which increases the risk of overfitting. To mitigate this risk, we used conservative feature reduction with consensus selection across algorithms and validated the final two-gene signature in two independent external datasets; nevertheless, larger multicenter cohorts and prospective studies are needed to confirm robustness and establish clinically actionable thresholds. Second, analyses were based on bulk microarray data, preventing assessment of cell type–specific expression; single-cell RNA-seq is needed to clarify the localization of TRAPPC13 and COPS5 and to validate immune cell signals inferred from bulk deconvolution. Moreover, microarrays have limited transcript coverage compared with RNA-seq and may miss low-abundance transcripts or isoform-level changes. Third, the findings are correlative and require experimental validation in neuronal or glial models to confirm causal roles in PD, including the observed associations with immune signatures and NMF-defined subgroups. Importantly, because the NMF subtypes were derived from only two genes (TRAPPC13 and COPS5), they should be interpreted as exploratory/hypothesis-generating. Validation using broader gene sets and independent cohorts will be required. In addition, because the present work is based on transcriptomic data, direct pathological or functional measures (e.g., α-synuclein aggregation burden and mitochondrial dysfunction/mitophagy readouts) were not available; integrating these endpoints will be essential in future validation studies. Fourth, the diagnostic model was derived from brain tissue and should be tested in peripheral samples and early-stage PD cohorts. Furthermore, despite validation in public datasets, independent prospective clinical-cohort validation – ideally using peripheral biospecimens – remains lacking. Finally, drug prediction results were obtained in silico and need experimental confirmation to determine whether candidate compounds can modulate TRAPPC13/COPS5 expression or alleviate PD-related pathology.
Conclusion
In summary, our integrative transcriptomic and machine learning analyses identified TRAPPC13 and COPS5 as robust downregulated vesicle transport-related biomarkers in Parkinson’s disease. A parsimonious two-gene model showed reproducible discriminative performance across the training cohort and two independent external validation cohorts, while nomogram, calibration, and decision curve analyses supported its potential value as an auxiliary tool for PD classification and stratification in retrospective datasets. Biologically, these genes were associated with vesicle trafficking/autophagy proteostasis programs and enabled an exploratory two-subtype stratification with distinct immune and pathway features, providing additional insight into PD heterogeneity and its immunopathological context. Further studies in larger multicenter prospective cohorts, clinically accessible biospecimens, and single-cell/spatial and functional models are needed to confirm clinical applicability and clarify causal links with á-synuclein pathology, mitochondrial dysfunction/mitophagy, and immune activation.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/cjn.2026.10587.
Data availability statement
The data and materials in the current study are available from the corresponding author on reasonable request.
Acknowledgments
Not applicable.
Author contributions
XQX and YBC contributed to the study design. XQX and YBC conducted the literature search. XQX acquired the data. XQX and YBC wrote the article and performed data analysis. XBX and YBC revised the article and gave the final approval of the version to be submitted. All authors read and approved the final manuscript.
Funding statement
None.
Competing interests
The authors declare that they have no potential conflicts of interest.
Ethical standards
Not applicable.
Consent for publication
Not applicable.
