Novel diagnostic biomarkers related to immune infiltration in Parkinson's disease by bioinformatics analysis

Abstract

Background: Parkinson's disease (PD) is Pengfei Zhang Liwen Zhao Pengfei Zhang Liwen Zhao a common neurological disorder involving a complex relationship with immune infiltration. Therefore, we aimed to explore PD immune infiltration patterns and identify novel immune-related diagnostic biomarkers.

Materials and methods: Three substantia nigra expression microarray datasets were integrated with elimination of batch effects. Differentially expressed genes (DEGs) were screened using the "limma" package, and functional enrichment was analyzed. Weighted gene co-expression network analysis (WGCNA) was performed to explore the key module most significantly associated with PD; the intersection of DEGs and the key module in WGCNA were considered common genes (CGs). The CG protein-protein interaction (PPI) network was constructed to identify candidate hub genes by cytoscape. Candidate hub genes were verified by another two datasets. Receiver operating characteristic curve analysis was used to evaluate the hub gene diagnostic ability, with further gene set enrichment analysis (GSEA). The immune infiltration level was evaluated by ssGSEA and CIBERSORT methods. Spearman correlation analysis was used to evaluate the hub genes association with immune cells. Finally, a nomogram model and microRNA-TF-mRNA network were constructed based on immune-related biomarkers.

Results: A total of 263 CGs were identified by the intersection of 319 DEGs and 1539 genes in the key turquoise module. Eleven candidate hub genes were screened by the R package "UpSet." We verified the candidate hub genes based on two validation sets and identified six (SYT1, NEFM, NEFL, SNAP25, GAP43, and GRIA1) that distinguish the PD group from healthy controls. Both CIBERSORT and ssGSEA revealed a significantly increased proportion of neutrophils in the PD group. Correlation between immune cells and hub genes showed SYT1, NEFM, GAP43, and GRIA1 to be significantly related to immune cells. Moreover, the microRNA-TFs-mRNA network revealed that the microRNA-92a family targets all four immune-related genes in PD pathogenesis. Finally, a nomogram exhibited a reliable capability of predicting PD based on the four immune-related genes (AUC = 0.905).

Conclusion: By affecting immune infiltration, SYT1, NEFM, GAP43, and GRIA1, which are regulated by the microRNA-92a family, were identified as diagnostic biomarkers of PD. The correlation of these four genes with neutrophils and the microRNA-92a family in PD needs further investigation.

Keywords: Parkinson’s disease; bioinformatics analysis; gene set enrichment analysis; hub genes; immune infiltration; weighted gene co-expression network analysis (WGCNA).

FIGURE 1

Flowchart of the study. GEO, gene expression omnibus; DEGs, differentially expressed genes; WGCNA, weighted gene co-expression network analysis; ssGSEA, single-sample gene set enrichment analysis; GSEA, gene set enrichment analysis; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein–protein interaction; PD, Parkinson’s disease; TF, transcriptional factor.

FIGURE 2

Identification of DEGs in the integrated dataset. (A) Volcano plot of all DEGs. The tomato nodes represent upregulated DEGs with p-value < 0.05 and logFC > 0.5; the cyan nodes represent downregulated DEGs with p-value < 0.05 and logFC < –0.5. (B) Heatmap of DEGs in PD samples vs. normal samples. Each row of the heatmap represents one gene, and each column represents one sample. The red and blue colors represent gene expression levels corresponding to upregulated and downregulated expression. DEGs, differentially expressed genes; FC, fold change; PD, Parkinson’s disease.

FIGURE 3

Gene ontology and KEGG pathway enrichment analyses of DEGs. (A) Barplot of KEGG analysis based on the obtained 263 genes. (B) The bubble diagrams show the top ten significantly enriched terms in KEGG analysis. The X-axis is the GeneRatio (gene count/gene size) of the term, and the Y-axis denotes the name of the term. The darker the color is, the smaller the adjusted p-value is. (C) Subnetwork showing the top five KEGG pathways and related genes. (D) The top 5 terms for BP, CC, and MF with p < 0.05 are shown. DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, gene ontology; BP, biological processes; CC, cell component; MF, molecular function.

FIGURE 4

The WGCNA process for the integrated dataset. (A) Clustering dendrogram of 38 PD substantia nigra tissue and 29 normal substantia nigra tissue gene expression patterns. (B) Analysis of the scale-free fit index (left) and the mean connectivity (right) for various soft-thresholding powers; the power value β was set as 7 for further analysis. (C) Clustering dendrograms of 5000 genes based on a dissimilarity measure (1-TOM). Seventeen co-expression modules were constructed with various colors under the gene tree, and similar modules were merged into twelve modules with a height cutoff of 0.25. Each color represents one module. (D) Heatmap of associations between modules and clinical traits. Correlation coefficients and p-values are shown in each cell, which were obtained by the intersection of rows and columns. The turquoise module correlated significantly with PD. (E) Visualization of the eigengene dendrogram and eigengene adjacency heatmap. Red indicates more similarity, and blue indicates less similarity. (F) Visualization of 400 random genes from the WGCNA network using a heatmap plot to depict the TOM among all modules included in the analysis. A redder background indicates a higher module correlation. (G) Scatter plot of module membership vs. gene significance for PD in the turquoise module. WGCNA, weighted gene co-expression network analysis; TOM, topological overlap matrix; PD, Parkinson’s disease.

FIGURE 5

Protein–protein interaction network and identification of candidate hub genes. (A) Venn plot showing the intersection between the DEGs and genes in the turquoise module, and 263 CGs were obtained. (B) PPI network of CMs by Cytoscape. The size and gradient color of circles are adjusted by the degree value, which reflects the connectivity between nodes. The size of circles has a positive correlation with the degree value. (C) PPI Count, revealing the number of adjacent nodes of the top 30 genes (ranked from low to high) based on the PPI network. (D) UpSet plot showing the intersection of ten algorithms, namely, MCC (top 30), MNC (top 30), Degree (top 30), EPC (top 30), BottleNeck (top 30), Closeness (top 30), Radiality (top 30), Stress (top 30), and Betweenness (top 30), and PPI Count (top 30). (E) Eleven candidate hub genes. All these candidate hub genes were found to be downregulated. PPI, protein–protein interaction; CGs, common genes; MCC, maximal clique centrality; MNC, maximum neighborhood component; EPC, edge percolated component.

FIGURE 6

Validation of candidate hub genes. (A) Validation of candidate hub genes in the GSE20292 dataset. (B) Validation of candidate hub genes in the GSE20681 dataset (*p < 0.05, **p < 0.01). PD, Parkinson’s disease.

FIGURE 7

Diagnostic efficacy of potential biomarkers (GAP43, GRIA1, NEFL, NEFM, SNAP25, and SYT1) for prediction of PD and GSEA based on expression levels of those potential biomarkers. ROC analysis of six hub genes (GAP43, GRIA1, NEFL, NEFM, SNAP25, and SYT1) for diagnosing PD in the integrated dataset (A) and the validation sets GSE20292 (B) and GSE26927 (C). Single-gene GSEA-KEGG pathway analysis of GAP43 (D), GRIA1 (E), NEFL (F), NEFM (G), SNAP25 (H), and SYT1 (I). PD, Parkinson’s disease; GSEA, gene set enrichment analysis; ROC, receiver operator characteristic.

FIGURE 8

Visualization and evaluation of immune infiltration levels based on ssGSEA. (A) Comparison of 16 immune cells between PD samples and control samples. (B) Box diagram of different immune function expression levels in the PD and control groups. (C) Pearson correlation analysis of 29 types of immune cells and immune-related functions with hub genes (*p < 0.05 and **p < 0.01). PD, Parkinson’s disease.

FIGURE 9

Visualization and evaluation of immune infiltration levels based on the CIBERSORT algorithm. (A) Comparison of 22 immune cells between PD samples and control samples. (B) Pearson correlation analysis of immune cell infiltration with hub genes (*p < 0.05 and **p < 0.01). PD, Parkinson’s disease.

FIGURE 10

Construction of the nomogram. (A) Construction of a nomogram for immune-related hub genes (GAP43, GRIA1, NEFM, and SYT1) for predicting the occurrence of PD. (B) Calibration curve estimates the prediction accuracy of the nomogram for PD patients. (C) The area under the curve (AUC) was 0.905. PD, Parkinson’s disease; AUC, area under the curve.

FIGURE 11

The miRNA-TF-mRNA interaction network. The red circles represent marker hub genes, the blue triangles represent miRNAs, and the green inverted cones indicate TFs. mRNAs, messenger RNAs; miRNAs, microRNAs; TFs, transcription factors.