The high‐confidence protein‐protein interaction (PPI) data with score >10000 were downloaded from the Human Protein Reference Database (HPRD, http://www.hprd.org/). The 2747 immune-related genes (HIRGs) were downloaded from the Amigo 2 database (http://amigo.geneontology.org/amigo) with immune as the key word. The 2078 PD-related genes (PDRGs) were downloaded from the DisGeNET database (https://www.disgenet.org/search). The FPKM expression profiles of GSE49036 and GSE20292 datasets were downloaded from the Gene Expression Omnibus (GEO) database. In the GSE49036 dataset, 8 normal and 14 PD samples were selected for data analysis, and 15 normal and 11 PD samples of GSE20292 dataset were selected for validation analysis. The clinical characteristics of GSE20292 and GSE49036 datasets were shown in Supplementary Table S1.

The fragments per kilobase of transcript per million fragments mapped (FPKM) values of gene expression in the GSE49036 dataset were log2-transformed, and the Pearson correlation of the two genes was calculated using the R package psych (version 2.1.9). Then, according to a threshold of |Pearson coefficient value| > 0.7 and FDR <0.05 to obtain the correlation among genes. Furthermore, based on the correlation between the filtered genes, mapped into the protein interaction network of the HPRD database, the common network was selected. The common network was drawn using Cytoscape software (version 3.8.2) (Shannon et al., 2003). According to the high connectivity score of genes in the common network, the number of four types of genes (PD, immune-PD, immune, and others) was counted. The Veen online tool (http://www.bioinformatics.com.cn/static/others/jvenn/example.html) was used to draw a Venn diagram of the correlation of protein-interacting genes, PDRGs, and immune genes. The R package ggplot (version 2.3.3.2) (Villanueva, 2019) was applied to draw a histogram of the four gene categories.

Based on the above-mentioned high connectivity score of genes in the common network, the core genes of PD-related genes were extracted, including PD and immune-PD genes, and their connected genes. Next, only the PD genes and their direct-acting genes were extracted as core genes. Subsequently, the number of four types of genes (PD, immune-PD, immune, other) were counted, and visualized by ggplot (Version 2.3.3.2) package (Villanueva, 2019). Notably, the immune-PD genes were both immune-related genes and PD-related genes. In the interaction network where the core gene was only PD gene, the correlation of different types of genes was calculated. Finally, the expression data of all genes in the network (the core genes are only the PD genes) were extracted and compared with the core genes. Wilcoxon rank-sum test was applied to compare the coexpression correlation coefficients between different gene groups (immune, immune-PD, PD, and other genes). Subsequently, to investigate the level of interaction between different gene groups with neighbors, cumulative distribution function (CDF) was utilized to assess the degree of the expression correlation for each gene group. The Pearson correlations of genes and the genes that related to the corresponding core genes were calculated, and the R package pheatmap was used to draw correlation heatmaps.

Modular clustering analysis was performed on the genes in the IOPDNC to obtain the central core genes using the GraphWeb online website (https://biit.cs.ut.ee/graphweb/). Subsequently, the number of core genes (i.e. PD-related genes) was adjusted to six in the GraphWeb database, with the rest set to default, and the modules with top five functional score values were selected as PD-related gene modules. Then, the R package ConsensusClusterPlus (version 1.54.0) (Wilkerson and Hayes, 2010) was used to perform consistent clustering analysis on the genes with the top 5 functional scores in the module, and the appropriate K value was selected based on the clustering results. In the external validation set GSE20292 dataset, the expression levels of modular genes were extracted in the same way and the accuracy of our screening of modular genes was validated against the same consistent clustering criteria.

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of different module genes were analyzed using the R package Clusterprofiler package (version 4.0.2) (Wu et al., 2021). According to the significance threshold p < 0.05 and count value, the enrichment analysis of each module was carried out separately, and the ggplot (version 2.3.3.2) was used for plotting. According to the website of Pathview (https://pathview.uncc.edu/), the immune-related pathway hsa04650 was selected to visualize the most immune-related pathways in the module.

Based on the 24 immune cell sets, the single sample Gene Set Enrichment Analysis (ssGSEA) algorithm was used to calculate the immune cell infiltration of the PD and normal samples, and the rank sum test was used to analyze the immune differences between the PD and normal samples of cell infiltration.

Peripheral blood mononuclear cells (PBMC) samples from eight normal samples and eight PD patients was collected using vacuum blood tubes containing EDTA anticoagulant in accordance with clinical blood collection techniques. Each PBMC sample was gently shaken repeatedly and loaded into a 4°C thermostat and transferred to the laboratory for subsequent manipulation according to biosafety requirements.

The total RNA of 16 PBMC samples (8 normal samples and 8 PD samples) was extracted to verify the results of the bioinformatics analysis. The top 1 gene of each module (module 1, module 2, module 3, module 4, module 5, and module all) was selected for qRT-PCR experiments. The total RNA of 16 samples was extracted with TRIzol Reagent (Life Technologies-Invitrogen, Carlsbad, CA, United States). Then, these total RNA were reverse transcription into cDNA with the SureScript-First-strand-cDNA-synthesis-kit (Genecopoeia, Guangzhou, China) prior to qRT-PCR. The primers of these genes for qPCR were as follows:

PSMB7-For:CATGGGTTCTGGCTCCTTGG; PSMB7-Rev:CTGGTCCCCTTCTTGTTGGG; GRIN1-For:CAAGAAGGAGTGGAATGGGATG; GRIN1-Rev:GCTCGTTGTT TATGGTTAGCGG; NME1-For:CAACCCTGCAGACTCCAAGC; NME1-Rev:GGTGAAACCACAAGCCGATC; SIN3B-For:ACCCTGCCACCTACAACGG; SIN3B-Rev:TTGTCAGAGGCGAC TGTATGTTTA; HABP4-For:GAGGCAGGCAGACTTCACAG HABP4-Rev:CGAACTCCACATCCACCCAT; STX1A-For:CAATGTGGAACACGCGGTAG; STX1A-Rev: ACA​GTG​GAG​GCG​ATG​ACG​AT.

The expression was uniformized to the internal reference GAPDH and computed employing the 2^−ΔΔCt method.

RIPA Lysis Bufferb (Servucebui) containing a protease inhibitor (Servucebio) was utilized to obtain the protein from tissues, subsequently, immunoblotting was performed. The bicinchoninic acid (BCA) quantification kit was applied to determine protein concentration of the cell lysates. The protein samples were loaded and separated by SDS-PAGE and shifted to PVDF membranes (Millpore, Sigma). Membranes were incubated with specific primary antibodies against PSMB7 (cst), GRIN1 (Affinity), NME1 (Affinity), SIN3B (Affinity), HABP4 (Proteintech), STX1A (BOSTER), and β-Actin (Proteintech) after blocked with 5% nonfat dry milk at 4 °C. Furthermore, secondary antibodies (IgG) were incubated at room temperature for 60 min and visualized using an ECL system. The dilution factors of the primary and secondary antibodies were shown in Supplementary Table S2.

A total of 1022077 gene relationship pairs were selected from the GSE49036 dataset. Then, the common network with 416 nodes and 281 edges was obtained from the GSE49036 dataset and HPRD database (Figure 1A, Supplementary Table S3). Among the common network genes of Figure 1A, the number of immune genes was 63; the number of PD genes was 77; the number of immune-PD genes was 39; the number of other genes were 237 (Figure 1B). Genes with high protein interaction connectivity scores were selected, and a Venn diagram of protein-interacting genes associated with PD and immune genes was drawn. There were 39 intersection genes between 2078 PD genes, 2740 immune genes, and 416 common network genes (Figure 1C). The results indicated that immune‐related genes played a vital role in the IOPDNC network. Together, these results suggested that immune-related genes might be important contributors for PD.

The core genes of PD, immune-PD genes, and their linked genes were extracted from the network of 3.1. Then, a network of these genes was constructed with 130 nodes and 91 nodes (Figure 2A, Supplementary Table S4). The genes whose core genes only were PD and the genes that were directly affected by PD were extracted to construct a network that including 87 nodes and 60 edges (Figure 2B, Supplementary Table S5). Totally five common genes were detected between 2078 PD genes, 2740 immune genes, and 87 core genes (Figure 2C). Furthermore, the number of genes in the four gene categories (PD, immune-PD, immune, and other genes) that the core genes only were PD of the PD gene interaction network were accounted (Figure 2D). Among the protein-interaction network, the number of immune genes was 5; the number of PD genes was 51; the number of immune-PD genes was 5; the number of other genes was 26 (Figure 2D). Moreover, there were significant differently expressed correlations between different gene groups, except immune-PD and PD genes (p = 0.39) (Figure 2E). In the PD genes and its linked genes, the correlation of immune genes, immune-PD genes and other genes was significant (Figure 2F). The correlation between PD genes and its linked genes was significant with |Pearson coefficient value| > 0.7 and FDR <0.05 (Figure 2G). Totally, the results showed that there were topological interactions and expression patterns among the correlations between PD- and immune-related genes.

The modules with top 5 functional scores were selected (Figure 3A, Supplementary Table S6). In the all of 5 modules, 25 PD genes, 11 immune genes, 5 immune-PD genes, and 46 other genes were contained.

In these 6 modules, the genes were divided into 2 clusters when the K = 4, and the expression of these module genes were higher in cluster 1 (Figures 3B–G). In these modules, the expression levels of module genes except PSMB1, KARS, TERT, ZBTB16, NFKB1, MAPKAPK2, MAP3K3, MAP2K4, IRS2, IQGAP1, HDAC1, EPB41L3, C1QBP, BRCA1, BAG2, FYN, MAPK6, MAP3K2, and FNY in the PD group were higher than that in the control group.

The same method was used to validate the accuracy of the screened module genes in the external validation set GSE20292 data set (Figure 4). These genes have better representation and can screen out patients at different stages. These results demonstrated that the model genes could distinguish PD and control samples well.

Module 1 genes were enriched in 104 GO BPs (including 12 GO CCs, 7 GO MFs) and 8 KEGG pathways, and these GO terms and KEGG pathways were mainly related to various metabolic processes, complexes, and disease pathways (Supplementary Figures S1A, B). Module 2 genes were enriched in 303 GO BPs, 55 GO CCs, 12 GO MFs, and 7 pathways, and these GO terms and KEGG pathways were mainly related to various cell migration, protein binding, and disease pathways (Supplementary Figures S1A, B). Module 3 genes were enriched in 60 GO BPs, 17 GOs CC, 51 GO MFs, and 2 KEGG pathways, and these GO terms were mainly related to various telomere. The KEGG pathways were base excision repair and non-homologous end-ioining (Supplementary Figures S1A, B). Module 4 genes were enriched in 198 GO BPs, 35 GO CCs, 32 GO MFs, and 32 KEGG pathways, and these GO terms and KEGG pathways were mainly related to various protein binding, response to stimulus and disease pathways (Supplementary Figures S1A, B). Module 5 genes were enriched in 141 GO BPs, 30 GO CCs, 35 GO MFs, and 12 KEGG pathways, and these GO terms and KEGG pathways were mainly related to various response to stimulus and signaling pathways (Supplementary Figures S1A, B). All of the 5 module genes was enriched in 312 GO BPs, 80 GO CCs, 68 GO MFs, and 49 KEGG pathways, and these GO terms and KEGG pathways were mainly related to various protein binding and disease pathways (Supplementary Figures S1A, B). Additionally, the natural killer (NK) cell mediated cytotoxicity pathway was showed in Figure 5, Fyn, Vav and PKC were significantly enriched.

To explore the differences in immune cell infiltration between the control and PD samples, the ssGSEA algorithm was performed. There were significant differences in aDC, eosinophils, neutrophils, and Th2 cells between the control and PD samples (Figure 6A). In addition, the genes of the 6 modules were significantly associated with these 4 differential immune cells (Figure 6B).

The mRNA and protein expression of 6 module genes were significantly higher expressed in the normal samples than that in the PD samples (Figure 7 and Figure 8), the detailed statistical results for qRT-PCR and western blot were shown in Supplementary Table S7 and Supplementary Table S8, the original bar charts of western blot were shown in Supplementary Figure S2. These results confirm that these module genes could act as potential diagnostic markers for PD.