By using the R package (Version 4.3.0) GEO query (10), we downloaded the Psoriasis GSE30999 (11) and GSE13355 (12) datasets from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The samples in the GSE30999 and GSE13355 datasets were all derived from Homo sapiens, with skin as the tissue source. The chip platform used in these datasets was GPL570, as detailed in Table 1 . Among them, GSE30999 contained 85 Psoriasis samples and 85 Control samples, while GSE13355 contained 58 Psoriasis samples and 122 Control samples. All Psoriasis samples and Control samples were included in this study. We collected cell aging-related genes, including Cellular Senescence-Related Genes (CSRGs), from the GeneCards database (13) (https://www.genecards.org/) and published literature. The GeneCards database provides comprehensive information on human genes. After using the term “Exosome” as a search keyword and keeping only exosome-associated genes (ERGs) with “Protein Coding” and “Relevance Score > 3,” a total of 158 ERGs were obtained. Likewise, we used “Exosome” as keywords on the PubMed website (https://pubmed.ncbi.nlm.nih.gov/) to identify cell aging-related genes in published literature (14), resulting in a total of 121 ERGs. A total of 274 ERGs were obtained after merging and removing duplicates, and detailed information is shown in Supplementary Table S1 . The R package sva (15) was used to remove batch effects from GSE30999 and GSE13355 to obtain the Combined GEO Datasets. The Combined Datasets included 143 Psoriasis samples and 207 Control samples. Finally, the R package limma (16) was used to standardize the Combined GEO Datasets, annotate probes, and normalize them. The expression matrices before and after batch effect removal were subjected to Principal Component Analysis (PCA) (17) to verify the effectiveness of this process. PCA is a data dimensionality reduction method that extracts the feature vectors (components) of data from high-dimensional data, transforming it into low-dimensional data and displaying these features in 2D or 3D graphs.

The Combined GEO Datasets categorized samples into the Psoriasis group and the Control group. The R package limma was used to perform differential analysis of genes between the Psoriasis group and the Control group. For the threshold of DEGs, we set |logFC| > 1 and adjusted p < 0.05. Genes with logFC > 1 and adjusted p < 0.05 were classified as up-regulated genes, while genes with logFC < -1 and adjusted p < 0.05 were classified as down-regulated genes. The method used for adjusting p-values was the Benjamini-Hochberg procedure. The results of the difference analysis were plotted by volcano plot through the R package ggplot2.

We integrate the GEO dataset to identify ERDEGs related to Psoriasis. A dataset is derived from variance analysis, selecting genes with |logFC| > 1 and adjusted p < 0.05. We identify DEGs associated with external secretion (ERGs) and determine their intersection to establish ERDEGs. This is visualized using Venn diagrams and heat maps generated by the R packages pheatmap and RCircos (18).

GO analysis (19)is a widely used method for large-scale studies of functional enrichment, focusing on Biological Processes (BP), Cellular Components (CC), and Molecular Functions (MF). Similarly, the KEGG (20) is a popular database that provides information on genomes, biological pathways, diseases, and drugs. GO and KEGG pathway enrichment analyses were conducted using the R package clusterProfiler (21) for ERDEGs. The criteria for statistical significance included a p-value of less than 0.05 and a false discovery rate (FDR) of less than 0.25.

GSEA (22) evaluates how genes are distributed in a predefined gene set, which is ranked based on their correlation with a phenotype. This analysis helps determine the genes’ contributions to that phenotype. In this study, we first ranked the genes from the Combined GEO Datasets by their logFC values. Then, we used the R package clusterProfiler to perform GSEA on all genes in the integrated GEO Datasets. The parameters used in the GSEA were as follows: the seed was set to 2023, with a minimum of 10 genes and a maximum of 500 genes per gene set. We accessed the c2 gene sets from the Molecular Signatures Database (MSigDB), specifically using the Cp. All. V2022.1. Hs. Symbols. The GMT file containing all Canonical Pathways (3050) was utilized for the GSEA. The GSEA screening criteria were an adjusted p-value of less than 0.05 and a false discovery rate (FDR) value (q-value) of less than 0.25, with p-value correction applied using the Benjamini-Hochberg method.

The PPI network is a crucial tool for understanding how proteins interact with one another and their roles in biological processes. In this study, we utilized the STRING database (23)(https://string-db.org/)to analyze the interactions of ERDEGs. This database offers known protein-protein interactions and predicts potential ones using various evidence sources, such as experimental data and computational predictions.

We constructed a PPI network for the 16 ERDEGs identified in our study. The network was built using a minimum interaction score greater than 0.15, which represents a low confidence level but allows for a broader exploration of potential interactions (24). The resulting network helps to visualize the complex web of interactions among these genes and may reveal molecular complexes or pathways that are relevant to psoriasis. To identify the most significant genes within the PPI network (referred to as Hub Genes), we employed five algorithms from the CytoHubba plug-in of Cytoscape software (25):Maximal Clique Centrality (MCC), Degree, Maximum Neighborhood Component (MNC), Edge Percolated Component (EPC), and Closeness. These algorithms help to score and rank genes based on their connectivity and importance within the network.

The top 10 genes from each algorithm were selected and compared. A Venn diagram was created to identify the intersection genes that were consistently ranked highly across all algorithms. These intersection genes were identified as exosome-related Hub Genes, likely playing central roles in the pathogenesis of psoriasis.

Transcription factors (TFs) control gene expression through interactions with Hub Genes at the post-transcriptional stage. We utilized the ChIPBase database (26) (https://rnasysu.com/chipbase3/index.php) to retrieve transcription factors, applying a filter condition where the number of samples found (both upstream and downstream) is greater than 9. The regulatory effects of transcription factors on Hub Genes were analyzed, and the mRNA-TF regulatory network was visualized using Cytoscape software. Additionally, miRNAs play a crucial regulatory role in biological development and evolution. They can regulate a variety of target genes, and the same target gene can also be regulated by multiple miRNAs. To analyze the Hub Genes and their relationship with miRNAs, we accessed the StarBase v3.0 database (27)(https://rnasysu.com/encori/) to identify miRNAs associated with the Hub Genes. We selected mRNA-miRNA interaction relationships recorded from at least two sources and visualized the mRNA-miRNA regulatory network using Cytoscape software.

To examine the differences in Hub Gene expression between the Psoriasis and Control groups, we created a comparison map based on their expression levels. To further analyze the differences in Hub Gene expression between the Psoriasis and Control groups in the GEO Datasets, we created a comparison map. Finally, we used the R package pROC to plot the ROC curve for the Hub Genes and calculate the Area Under the Curve (AUC). This analysis aimed to evaluate how Hub Gene expression affects the occurrence of Psoriasis. The AUC of the ROC curve typically ranges from 0.5 to 1, with values closer to 1 indicating better diagnostic performance.

Single-Sample Gene-Set Enrichment Analysis (ssGSEA) (28) quantifies the relative abundance of each immune cell infiltrate. First, the types of infiltrating immune cells were labeled, including Activated Dendritic Cells, Gamma Delta T Cells, Natural Killer Cells, and Activated CD8⁺ T Cells. Next, the R package ggplot2 was employed to create comparison maps illustrating the expression differences of immune cells between the Control and Psoriasis groups in the Combined GEO Datasets. Following this, the R package ggplot2 was utilized to generate comparative plots that illustrate the variances in immune cell expression between the Control and Psoriasis groups in the Combined GEO Datasets. The correlation between immune cells was calculated using the Spearman algorithm. The R package pheatmap was then employed to create a heatmap displaying the correlation analysis results of the immune cells. The correlation between Hub Genes and immune cells was also calculated using the Spearman algorithm, and a correlation bubble plot was drawn using the R package ggplot2 to show the correlation analysis results of Hub Genes and immune cells.

All data processing and analysis in this article were based on R software (Version 4.3.0). To compare continuous variables between two groups, we used the independent Student’s T-Test for normally distributed variables, unless stated otherwise. We applied the Mann-Whitney U Test (Wilcoxon Rank Sum Test) to analyze differences in non-normally distributed variables. We used the Kruskal-Wallis test to compare three or more groups. We used Spearman correlation analysis to assess the correlation coefficient between different molecules. All p-values were two-sided unless specified otherwise. A p-value of less than 0.05 was considered statistically significant.

This study combines two datasets: GSE13355, which includes 58 psoriasis patients and 122 controls, and GSE30999, which consists of 85 psoriasis patients and 85 controls. It aims to systematically identify DEGs in psoriasis and explore their biological functions using GSEA. Additionally, a PPI network was constructed to identify key hub genes. Regulatory networks involving mRNA-TF and mRNA-miRNA were also analyzed to reveal potential transcriptional regulatory mechanisms.We conducted functional enrichment analysis (GO and KEGG), ROC curve assessment, and immune infiltration analysis to evaluate the diagnostic value of these genes in psoriasis. This comprehensive evaluation sheds light on their roles in the immune microenvironment and offers a new perspective on psoriasis pathogenesis and potential biomarkers (Figure 1).

The R package sva was used to remove batch effects from the GSE30999 and GSE13355 datasets, resulting in a combined GEO dataset. First, we used a distribution boxplot ( Figures 2A, B ) to compare the expression values of the datasets before and after batch effect removal. Next, we employed a PCA plot ( Figures 2C, D ) to compare the distribution of low-dimensional features before and after the removal of batch effects. The results from the distribution boxplot and PCA plot indicated that the batch effects in the psoriasis dataset were effectively eliminated following the removal process.

The Combined GEO Datasets were split into two groups: psoriasis and control. We used the R package limma to analyze gene expression differences between the psoriasis and control groups in the Combined GEO Datasets. The results are as follows: In the Integrated GEO Dataset (Combined Datasets), a total of 921 differentially expressed genes (DEGs) met the thresholds of |logFC| > 1 and adjusted p < 0.05. According to the variance analysis results shown in the volcano plot (Figure 3A), there are 555 upregulated genes (logFC > 1 and adj. P < 0.05) and 366 downregulated genes (logFC < 1 and adj. P < 0.05). We identified ERDEGs by intersecting DEGs that met the criteria of |logFC| > 1 and adjusted p < 0.05 with known ERGs, as displayed in the Venn diagram ( Figure 3B ).We identified 16 ERDEGs, which are: CD274, SERPINB3, CD24, RAB27A, BIRC5, CXCL13, SDC4, LGALS3BP, POSTN, SMPD3, RAB27B, IVL, PLCB4, EPCAM, FCGBP, and MYOC. According to the intersection results, the expression differences of ERDEGs between different sample groups in the Combined GEO Datasets were analyzed, and the R package pheatmap was used to create a heatmap to display the analysis results ( Figure 3C ). Finally, the location of the 16 ERDEGs on human chromosomes was analyzed using the R package RCircos to draw a chromosome localization map ( Figure 3D ). Chromosome mapping revealed that most ERDEGs were found on chromosomes 1, 17, 18, and 20. Specifically, IVL and MYOC are on chromosome 1, BIRC5 and LGALS3BP are on chromosome 17, RAB27B and SERPINB3 are on chromosome 18, and PLCB4 and SDC4 are on chromosome 20.

We performed GO and KEGG analyses to explore the relationships between BP, CC, MF, and pathways related to 16 ERDEGs in psoriasis. We used the 16 ERDEGs for GO and KEGG, with specific results presented in Table 2 . The results indicated that the 16 differentially expressed early response genes (ERDEGs) in psoriasis were primarily enriched in processes such as positive regulation of exocytosis, exosomal secretion, extracellular exosome biogenesis, and extracellular vesicle biogenesis. The analysis also revealed enrichment in several KEGG biological pathways, including pancreatic secretion, amoebiasis, cell adhesion molecules, chemokine signaling pathways, and apoptosis in various species. The results of the GO and KEGG enrichment analyses are illustrated using bar graphs and bubble plots ( Figures 4A, B ).Additionally, a network diagram illustrating biological processes (BP), cellular components (CC), molecular functions (MF), and KEGG pathways was generated based on the GO and KEGG enrichment analyses ( Figures 4C–F ).The lines illustrate the corresponding molecules and annotations of the respective entries, with larger nodes indicating a greater number of molecules contained within those entries.

We used GSEA to assess how gene expression levels in the Combined GEO Datasets impact Psoriasis. This analysis explored the relationships between gene expression and associated biological processes, affected cellular components, and molecular functions ( Figure 5A ). The results indicate that the integration of the Combined GEO Datasets reveals significant enrichment in pro-inflammatory and profibrotic mediators ( Figure 5B ), the IL-23 pathway ( Figure 5C ), photodynamic therapy-induced NF-κB survival signaling ( Figure 5D ), apoptosis ( Figure 5E ), and other related biological functions and signaling pathways.

The GSEA analysis results reveal several biological mechanisms related to psoriasis, including the interferon signaling pathway, antimicrobial peptides, SARS-CoV-2 signaling, pro-inflammatory and profibrotic mediators, and cell cycle checkpoints.These findings offer valuable insights into the pathogenesis of psoriasis and suggest potential targets for future research and treatment (Table 3).

First, we conducted a protein-protein interaction analysis and constructed a PPI network of 16 ERDEGs using the STRING database ( Figure 6A ). The results of the PPI network indicated that 15 ERDEGs were related: CD274, SERPINB3, CD24, RAB27A, BIRC5, CXCL13, SDC4, LGALS3BP, POSTN, SMPD3, RAB27B, IVL, EPCAM, FCGBP, and MYOC. Next, we used five algorithms from the CytoHubba plug-in of Cytoscape software to calculate the scores of the 15 ERDEGs. The ERDEGs were then ranked based on these scores. The five algorithms used were: Maximal Clique Centrality (MCC), Degree, Maximum Neighborhood Component (MNC), Edge Percolated Component (EPC), and Closeness. Then, the top 10 ERDEGs from the five algorithms were used to draw protein-protein interaction networks, namely MCC ( Figure 6B ), MNC ( Figure 6C ), Degree ( Figure 6D ), EPC ( Figure 6E ), and Closeness ( Figure 6F ). In this representation, the circle colors range from red to yellow, indicating scores from high to low. Finally, the genes from the five algorithms were interfaced, and a Venn diagram ( Figure 6G ) was drawn for analysis. The intersection genes of the algorithms were exosome-related hub genes, and the 9 hub genes were: POSTN, EPCAM, CD274, CD24, SERPINB3, CXCL13, SMPD3, BIRC5, and RAB27A.

First, we obtained TFs associated with Hub Genes from the ChIPBase database, using a screening criterion of more than 9 samples found (upstream + downstream). We constructed and visualized the mRNA-TF Regulatory Network using Cytoscape software ( Figure 7A ). This included 9 Hub Genes and 32 TFs, with specific information provided in Supplementary Table S2 . Next, we identified miRNAs related to Hub Genes from the StarBase database, using screening criteria that required mRNA-miRNA interaction relationships to be recorded from at least two sources. The mRNA-miRNA Regulatory Network was constructed and visualized by Cytoscape software ( Figure 7B ).This analysis identified 7 Hub Genes and 48 miRNAs, with detailed information available in Supplementary Table S3 .

Figure 8A illustrates the differential expression of nine Hub Genes in the Combined GEO Datasets, comparing the Psoriasis group to the Control group. The results indicate that the expression levels of nine Hub Genes in the Psoriasis and Control groups of the Combined GEO Datasets were highly statistically significant (p < 0.001): POSTN, EPCAM, CD274, CD24, SERPINB3, CXCL13, SMPD3, BIRC5, and RAB27A. Finally, the R package pROC was used to generate ROC curves based on the expression levels of Hub Genes in the Combined GEO Datasets. The ROC curves ( Figures 8B–J ) demonstrated that the expression levels of eight Hub Genes—POSTN, CD274, CD24, SERPINB3, CXCL13, SMPD3, BIRC5, and RAB27A—exhibited high accuracy (AUC > 0.9) in classifying the Psoriasis and Control groups. The expression level of EPCAM demonstrated moderate accuracy (0.7 < AUC < 0.9) in distinguishing between the Psoriasis and Control. groups.

The expression matrices from the Combined Datasets were used to calculate the abundance of immune cell infiltration for 28 immune cells using the ssGSEA algorithm. First, we displayed group comparison plots showing differences in immune cell infiltration abundance among various groups. The group comparison diagram ( Figure 9A ) indicated that all 28 immune cells were statistically significant (p< 0.05), including: activated B cell, activated CD4⁺ T cell, activated CD8⁺ T cell, activated dendritic cell, CD56 bright natural killer cell, CD56 dim natural killer cell, central memory CD4⁺ T cell, central memory CD8⁺ T cell, effector memory CD4⁺ T cell, effector memory CD8⁺ T cell, eosinophil, gamma delta T cell, immature B cell, immature dendritic cell, macrophage, mast cell, MDSC, memory B cell, monocyte, natural killer cell, natural killer T cell, neutrophil, plasmacytoid dendritic cell, regulatory T cell (Treg), T follicular helper cell (Tfh), type 1 T helper cell, type 17 T helper cell, and type 2 T helper cell. Next, a correlation heatmap illustrated the results of immune cell infiltration abundance from the immune infiltration analysis in the Combined GEO Datasets ( Figure 9B ). The results indicated that most immune cells exhibited strong positive correlations, with regulatory T cells (Treg) and MDSC showing the most significant correlation (r = 0.877, p< 0.05). Finally, the correlation between Hub Genes and the abundance of immune cell infiltration was shown by a correlation bubble plot ( Figure 9C ).The correlation bubble plot results indicated that most immune cells displayed strong correlations, with CD274 and activated CD4⁺ T cells showing the most significant positive correlation (r = 0.829, p< 0.05).