Core Differentially Expressed Genes in Psoriasis Lesions: An Integrated Analysis of Four GEO Datasets

Introduction

Psoriasis is a chronic, immune mediated inflammatory skin disorder that affects approximately 2% of the global population,¹ with marked variability in prevalence across different geographic groups. A study from Denmark² and the US³ reports a prevalence rate of psoriasis of almost 3%. Guinot et al⁴ estimated the prevalence in France at 4.7%.⁴ By contrast, the condition appears to be less common in North African countries, with rates ranging between 2.3% and 3%,^5,6 including a documented prevalence of 2.3% in Morocco.⁶

Microarray and high throughput sequencing technologies are robust and reliable approaches for the fast and accurate identification of differentially expressed genes (DEGs) in human patients and animal models. Several platforms exist to generate large scale transcriptomic datasets, which are publicly available in data repositories such as the Gene Expression Omnibus (GEO). These publicly available data sets provide resources for secondary data analysis and hypothesis generation.⁷

In recent years, numerous transcriptomic studies have investigated the gene expression landscape of psoriasis using microarray technologies.⁸ These efforts have contributed to a deeper understanding of the molecular mechanisms underlying the disease by highlighting distinct gene expression signatures associated with psoriatic lesions.^8,9

Psoriasis pathogenesis is increasingly recognized as the result of a dynamic interplay between keratinocytes, dermal fibroblasts, melanocytes, and immune cells. Previous large scale analyses of Affymetrix transcriptomic data¹⁰ highlighted not only keratinocyte proliferation but also immune mediated pathways, including IL-17/IL-23, TNF, and JAK-STAT signaling. However, Affymetrix platforms have limitations, as probe coverage is incomplete and novel genes or isoforms cannot be identified.¹⁰

More recent reviews have emphasized the central role of inflammatory triggers and cytokine signaling networks.^11,12 Topical corticosteroids, a cornerstone of treatment of psoriasis, exert broad anti-inflammatory effects primarily by binding to cytosolic glucocorticoid receptors, leading to the transrepression of key pro-inflammatory transcription factors like NF-κB and AP-1, thereby suppressing the production of cytokines (eg, IL-17, IL-23, TNF-α) and mediating vasoconstriction. Vitamin D analogues such as calcipotriene also widely used in the psoriasis treatment, act by binding to the vitamin D receptor, modulating keratinocyte differentiation and proliferation, and also exhibit immunomodulatory properties by inhibiting T-cell activation and the Th17 pathway.¹³ Novel biologic and small molecule therapies represent a paradigm shift towards precision targeting, with agents designed to neutralize specific cytokines central to the IL-23/Th17 axis or inhibit intracellular signaling pathways (PDE4 inhibitors, JAK-STAT inhibitors), thereby disrupting the inflammatory cascade at a more upstream and specific point compared to conventional treatments.¹⁴

For this purpose, we performed an analysis of four transcriptomic datasets to determine commonalities across different datasets to underscore the relevant pathways in the molecular landscape of psoriasis and offer potential targets for future functional validation and therapeutic interventions.

Methods

Psoriasis Gene Expression Data Sources

A comprehensive search for psoriasis vulgaris gene expression data was conducted on the NCBI GEO database, focusing exclusively on publicly available human skin biopsy samples. This study was exempted from ethical approval by Sfax University in Tunisia, as it was based exclusively on the analysis of publicly available data and did not involve direct contact with human participants or access to identifiable personal information. The search utilised keywords such as “Psoriasis vulgaris” and “RNA”, yielding 125 datasets, of which only 14 met the inclusion criteria: mRNA expression data from naïve patients, obtained from lesional skin biopsies, and compared with healthy control skin. Datasets were excluded if they involved miRNA, non-coding RNA, pre-treated patients, RNA sequencing from cultured keratinocytes, or perilesional skin samples. Ultimately, seven datasets were initially selected, and four were retained for final analysis as they specifically compared lesional vs. non-lesional skin. These included GSE34248 and GSE41662,¹⁵ GSE50790,¹⁶ and GSE6710.¹⁷ The excluded datasets were omitted due to their inclusion of control samples.

Identification of Differentially Expressed Genes

An online interactive web tool, GEO2R,¹⁸ was used to analyze the raw data of microarrays and identify DEGs between patient groups. GEO2R uses moderated t-statistics to compare gene expression levels in different groups. The p-value < 0.05 and logarithmic fold change |log2FC| ≥ 2 were used as the threshold to obtain statistically significant DEGs. Hence, upregulated genes (p-value < 0.05, log2FC ≥ 2) and downregulated genes (p-value < 0.05, log2FC ≤ −2) were grouped depending on their expression levels in respect to the cut-off values. Importantly, prior to cross-dataset comparison, each platform’s probe identifiers were first mapped to their corresponding official HGNC Gene Symbols using the platform-specific annotation tables available within GEO2R. All subsequent integration and Venn diagram intersections were therefore performed at the Gene Symbol level, not at the raw probe level, ensuring valid cross-platform comparability.

The identification of Differentially Expressed Genes (DEGs) was initiated by analyzing the raw microarray data, with gene expression levels compared between samples. This was achieved using the publicly available online interactive web tool, GEO2R,¹⁸ which simplifies the differential expression analysis for publicly archived GEO datasets. GEO2R was configured to utilize the GEOquery package for parsing the processed data into R structures, and subsequently, the robust limma (Linear Models for Microarray Analysis) package was employed. The analysis relied on moderated t-statistics from limma for testing, which intrinsically handles multiple-testing corrections on p-values to mitigate the occurrence of false positives.

To define the DEGs, stringent cut-off values were applied. The filtering process included the following steps: Statistical Significance Filter: A p-value <0.05 was required. This criterion was used to establish that the change in expression was unlikely to have occurred by chance. Magnitude of Change Filter: A minimum absolute logarithmic fold change of ∣log2FC∣≥2 was set, defining a substantial difference in expression. Gene Classification: Genes were subsequently classified as either upregulated (p-value <0.05 and log2FC≥2) or downregulated (p-value <0.05 and log2FC≤−2). These stringent thresholds were applied deliberately to ensure that only the most prominent probes were identified, confirming that the final list of DEGs reflected the most substantial alterations in the transcriptome.

Experimental Design

The grouping within the datasets was kept unmodified as described in each corresponding dataset. This allowed for the pooling of the probesets from the three datasets in respect to their pathology in two subgroups. The first subgroup corresponded to patients with lesions (lesional - L).

The second group was composed of the same patients’ skin biopsies without lesions and was grouped as non-lesional (NL). Patient numbers in each dataset, general information on GEO datasets, and platforms are provided in Table S1. Lesional and non-lesional samples of each dataset are listed in Table S2.

Analysis of Differential Gene Expression Subgroups

Following the initial data processing, the pooled sample collection was organized into the three defined pathology groups: control (Ctl), lesional (L), and non-lesional (NL). A key focus of the study was the comparison between the two patient derived subgroups, L and NL. Specifically, differential gene expression analysis was conducted by comparing the L group against the NL group. This intra-patient comparison was performed to isolate the gene expression changes specifically associated with the active psoriatic lesion, thereby minimizing variation introduced by differences in genetic background or systemic factors between individuals. Up and downregulated probes were visualized with volcano plots using the bioinfokit tool.¹⁹ Common DEGs for up and downregulated probesets in the four datasets were identified with online tool “Bioinformatics and Evolutionary Genomics” and visualized using Venn diagrams.²⁰ The intersection was performed exclusively at the annotated Gene Symbol level, using the platform specific annotation tables provided within GEO2R. This conservative vote counting strategy maximizes specificity, yielding a high confidence core signature of consistently dysregulated genes across all four independent datasets.

Protein–Protein Interaction Network Construction, Functional Enrichment, and Pathway Analysis

Protein–protein interaction (PPI) networks were constructed for differentially expressed genes (DEGs), followed by Gene Ontology (GO) enrichment and pathway analysis of up- and downregulated probe sets using the Enrichr classification system.²¹ Enrichr, a user-friendly gene set enrichment analysis tool, facilitated the functional annotation of DEGs across three GO categories: biological processes (BP), molecular functions (MF), and cellular components (CC). Additionally, Enrichr enabled comparative pathway analysis by integrating multiple pathway resources, ensuring a comprehensive interpretation of functional associations. To further investigate protein interactions, a PPI network was generated using the STRING database (https://www.string-db.org/), with the organism parameter set to Homo sapiens, allowing the retrieval and analysis of the minimum required interaction score at high-confidence (0.700) PPI data. Additionally pathways including KLK13 probe, formation of the cornified envelope, keratinization, and developmental biology, were clustered using k-means analysis.

Results

Differential Expression Analysis

The up- and downregulated probe sets for each dataset were visualized using volcano plots (Figure S1). Among the statistically significant differentially expressed probe sets across the four datasets (GSE34248, GSE41662, GSE50790, and GSE6710), a comparative analysis using Venn diagrams identified 32 and 3 major overlapping up and down regulated probe sets respectively (all identified using nominal p-value < 0.05 and |log2FC| ≥ 2 per dataset, the requirement for concordance across four independent datasets served as the primary stringency filter) in the L vs. NL comparison across studied datasets (Figure 1). In Tables 1 and 2 are mentioned the identified upregulated and downregulated genes respectively. Interestingly, among them, several genes are involved in keratinization and cornified envelope formation. Among the 32 upregulated genes, notable examples include TGM1 (present in all 4 datasets), KRT16 (4/4 datasets), DSC2 (4/4 datasets), and S100A7 (4/4 datasets), all well-established markers of keratinocyte activation and epidermal remodeling. KLK13, also consistently upregulated across all four datasets, was identified as one component of the broader protease network. Among the 3 downregulated genes, MMP20 and KLK4 were identified in the 3-dataset overlap (GSE34248, GSE41662, GSE6710).

Table 1 The Total Number and List of Common Probesets of Up-Regulated Probesets in Respect to Venn Analysis

Table 2 The Total Number and List of Common Probesets of Down-Regulated Probesets in Respect to Venn Analysis

Figure 1 Venn’s diagram showing the common (A) upregulated, and (B) downregulated DEGs from across four datasets. The number of probes at each junction of four datasets is indicated within the intersections. Pathway analyses were undertaken for all commonalities; however, the intersection across four datasets for upregulated and downregulated DEGs was evaluated in greater detail (L, lesional; NL, non-lesional).

Pathway Analysis

Pathway enrichment analysis using EnrichR tools to identify key biological pathways associated with commonly upregulated and downregulated probesets in the L vs. NL skin comparison. Among the most significantly upregulated genes (n = 32), the “Formation of Cornified Envelope” pathway ranked as the top pathway, followed by “Keratinization”. Additionally, “Neutrophil Degranulation” was the immune related pathway linked to psoriasis, appearing in the seventh position (Figure 2a). These findings confirm that keratinization and cornified envelope formation represent the primary and most statistically robust biological processes dysregulated in psoriatic lesions (p = 1.53 × 10⁻⁶ and p = 1.93 × 10⁻⁵, respectively). The “Neutrophil Degranulation” pathway, in which KLK13 participates, represents a secondary immune-related finding. Conversely, among the commonly downregulated probe sets, pathways related to “Mitochondrial Uncoupling” were among the top 10 statistically significant pathways (Figure 2b). The combined power of datasets and patients revealed evident downregulation of glucocorticoid pathways “prednisone ADME” and “glucocorticoid biosynthesis”, which are currently a targeted therapy for psoriasis.²²

Figure 2 Enrichr Pathway Enrichment Analysis: Highest-Ranked Significant Pathways for (A) Up and (B) Down Regulated Genes.

Hub Protein Identification

For the L vs. NL skin comparison, hub genes were identified using the STRING database, and the PPI network was visualized (Figure 3). K-means analysis of upregulated genes revealed three interrelated clusters, primarily encompassing genes involved in skin formation, maintenance, and immune system activation (Figure 3a). In contrast, analysis of downregulated genes identified eight clusters, with six exhibiting interrelated functional associations (Figure 3b).

Figure 3 PPI network and K-means clustering were conducted using pathways that included KLK13 (A) upregulated genes identified three major interrelated clusters, (B) downregulated genes identified six clusters. Each circle represents a probe.

Discussion

To date, several genes have been implicated in psoriasis vulgaris described as a common autoinflammatory genetic disease. The known involved genes play critical roles in inflammatory responses, keratinocyte differentiation and skin barrier function.^23–27

Our analysis identified 32 consistently upregulated genes across four datasets, many of which are involved in keratinization and cornified envelope formation, confirming earlier findings.¹⁰ Tables 3 and 4 present all significantly upregulated and downregulated genes, along with annotation in red indicating genes that were previously reported versus newly identified probes.

Table 3 The Most Differentially Up-Regulated Genes in Lesional Psoriatic Skin Compared to Healthy Skin. Red-shaded Boxes Indicate Genes Commonly Identified Across Datasets

Table 4 Significantly Down Regulated Genes in Lesional Psoriatic Skin Compared to Healthy Skin. Red-shaded Boxes Indicate Genes Commonly Identified Across Datasets

This study revealed multiple genes potentially involved in the pathophysiology of the disease. The common 32 genes among 4 datasets revealed numerous pathways, among which the first and second most statistically significant were keratinization (p-value = 1.53E-6) and cornified envelope formation (p-value = 1.93E-5). Statistically significant pathways showed increased upregulation for TGM, KRT16, and DSC2, which were further analyzed to understand their importance in psoriatic lesional skin biopsies.

TGM1, encoding transglutaminase 1, is essential for forming the cornified envelope during terminal keratinocyte differentiation, and its dysfunction compromises skin barrier integrity, contributing to disease onset.³⁵ KRT16, a type I keratin, is markedly overexpressed in psoriatic lesions, where it promotes keratinocyte hyperproliferation and abnormal differentiation,³⁶ inflammation, though its dysregulation can enhance immune responses and tissue damage.³⁷

Importantly, when comparing our results to previous studies, we also detected immune-related mediators such as CXCL8, CCL20, and IL36G, which have established roles in neutrophil recruitment, IL-17–driven inflammation, and epidermal crosstalk.¹¹ Enrichment analysis further revealed activation of IL-17, IL-23, NF-κB, and JAK-STAT pathways, consistent with current therapeutic targets in psoriasis.¹² These results highlight that psoriasis is not limited to keratinocyte hyperproliferation, but rather reflects coordinated deregulation of keratinocyte, immune cell signaling. KLK13 was also upregulated and identified in “Developmental Biology” (p-value = 0.017). Dysregulation of the KLK cascade is responsible for skin inflammatory diseases (Figure 4). KLK are serine proteases encoded by 15 different genes. Abnormal activation of the KLK proteolytic cascade is reported in psoriasis. In atopic dermatitis, an inflammatory skin disease, the upregulation of some KLK may induce activation of NK-κβ pathway and IL-8.³⁸

Figure 4 Dysregulated Kallikrein (KLK) Signaling in Psoriasis Pathogenesis. Left (KLK13 Upregulation): Increased KLK13 promotes neutrophil degranulation and IL-1β activation (amplified by LL-37). Proteolytic cleavage of Desmocollin-2 (DSC2) by KLK13 reduces keratinocyte cohesion. Right (KLK4 Downregulation): Impaired MMP20-KLK4 axis activation leads to reduced cleavage of Corneodesmosin and Desmoglein-1, inhibiting physiological desquamation.

The serine protease inhibitor PI3 (elafin) is also upregulated in psoriasis, acting to regulate protease activity and mitigate inhibitors to protect tissues from damage caused by excessive leukocyte activity.³⁹ Through this mechanism, elafin may influence the activity of kallikrein related peptidases such as KLK13, which is also upregulated in psoriasis. KLK13 may modulate epidermal desmosomal adhesion by proteolytically processing DSC2,⁴⁰ a desmosomal cadherin critical for keratinocyte cohesion.⁴¹ Disruption in these pathways through genetic mutation or altered expression underlies the chronic inflammation, abnormal keratinocyte behavior, and skin barrier dysfunction that characterize psoriasis. Notably, KLK13 is enriched in the “Neutrophil Degranulation” pathway (p = 0.006) and contributes to IL-1β activation via neutrophils,⁴² a response that can be amplified by LL-37 in psoriatic skin.⁴³ Both elafin and KLK13 are upregulated in lesional and non-lesional psoriatic skin,^44,45 with serum KLK13 levels correlating with Psoriasis Area Severity Index scores.⁴⁵

To investigate the molecular interactions among the key genes implicated in psoriasis, a STRING analysis was performed to map their functional network. The results revealed that these genes are interconnected not only with one another but also with genes from the top ten statistically enriched downstream pathways (Table 2 and Table 4), suggesting a broader functional network underlying disease pathology. These pathways were significantly associated with interleukin signaling (p = 6.85 × 10⁻⁴), chemokine receptor activity (p = 0.003), and immune system regulation (p = 0.004), highlighting their roles in mediating inflammatory responses. Notably, CXCL8 (IL-8), STAT1, and LCN2 emerged as central hub genes within these pathways. CXCL8 is a chemokine crucial for neutrophil recruitment and activation, processes that are heightened in psoriatic skin due to increased CXCL8 expression by keratinocytes, thereby amplifying local inflammation.⁴⁶ STAT1, a pivotal transcription factor in the JAK-STAT signaling cascade, orchestrates responses to interferons and other cytokines, with its upregulation and post-translational modification previously confirmed in psoriatic lesions.^47,48 LCN2 (lipocalin-2) contributes to innate immunity and modulates both neutrophil function and keratinocyte behavior, reinforcing its relevance in psoriasis pathogenesis.⁴⁹

Similar regulation is observed with other kallikreins, such as KLK4, which can cleave corneodesmosin, desmocollin-1, and desmoglein-1, key proteins involved in desquamation and epidermal integrity. Furthermore, transglutaminases (TGM1, TGM3, and TGM5) in the epidermis are essential for skin barrier formation and are activated by cathepsin D. Mutations in these enzymes lead to distinct skin disorders, including ichthyosis (TGM1), hair abnormalities (TGM3), and acral peeling skin syndrome (TGM5), emphasizing the broader relevance of protease regulation in skin homeostasis.

While our strict four-dataset intersection yielded a limited number of commonly downregulated genes, analyzing the intersection of three datasets (GSE34248, GSE41662, and GSE6710) revealed additional biologically relevant patterns, notably the coordinated downregulation of the MMP20-KLK4 protease axis. MMP20 is known to activate pro-KLK4, which subsequently activates other kallikrein zymogens crucial for maintaining epidermal homeostasis and desquamation. The concurrent downregulation of MMP20 and KLK4 observed in this three dataset overlap suggests a disrupted proteolytic activation cascade that may impair normal epidermal remodeling. Furthermore, other downregulated genes identified in this subset, such as TGM4, SLC41A1, and SLC45A3, point to altered keratinocyte function and ion transport. However, because these targets did not meet our stringent four dataset core criteria, their role in psoriasis pathogenesis requires cautious interpretation and highlights the need for further functional validation to confirm their biological relevance.

Further analysis of datasets GSE50790 and GSE6710 revealed downregulation of leptin (LEP), a hormone with wide-ranging effects on skin physiology. Leptin receptors have been identified in the epidermis, particularly in basal keratinocyte and in dermal papilla cells of hair follicles.⁵⁰ Leptin promotes keratinocyte and fibroblast proliferation, facilitates epithelialization, and enhances collagen synthesis, all of which are essential for skin regeneration and barrier maintenance.⁵¹ Additionally, leptin upregulates human β-defensin 2, reinforcing the skin’s antimicrobial defence.⁵² However, literature reports on leptin levels in psoriasis are conflicting; patients with severe psoriasis often exhibit elevated leptin levels, whereas those with milder forms may show reduced expression.⁵² In our analysis, leptin was consistently downregulated across both datasets, alongside adiponectin (GSE41662), another adipokine with anti-inflammatory properties. This dual downregulation may contribute to impaired immune regulation and barrier dysfunction in psoriatic skin. However, because disease severity was not stratified in the datasets, direct comparisons across samples remain limited. Additionally, among the downregulated pathways was “carnitine shuffle”, which is essential for the transport of long chain fatty acids through the membrane of mitochondria, which is crucial for skin beta-oxidation supply, pointing to an increased energy demand, that potentially cannot be maintained due to pathology. Collectively, the observed downregulation of key proteases, differentiation regulators, transporters, and adipokines underscores a systemic shift away from normal epidermal homeostasis and immune competence in psoriasis. These findings highlight potential targets for future functional studies aimed at restoring skin integrity and modulating inflammation.

A significant finding of this study is that if prospective trials incorporate evaluations of patient severity, drug targets may be addressed with enhanced specificity and effectiveness.

A notable methodological limitation of this study is our reliance on the direct intersection of independently generated DEG lists (the Venn diagram approach) rather than employing a formal pooled meta-analysis algorithm (such as RankProd or DERGA). While our stringent intersection criteria successfully isolated a highly robust core signature of 32 upregulated genes, this conservative approach is mathematically susceptible to the limitations and size of the smallest dataset, inevitably yielding a higher false-negative rate. Consequently, our analysis represents a highly filtered subset of biomarkers, and it is highly likely that many biologically relevant genes that a standard meta-analysis would detect were excluded. Future studies integrating raw expression matrices to calculate pooled effect sizes are warranted to capture the broader, more comprehensive transcriptomic landscape of psoriasis.

Conclusion

The use of Affymetrix microarrays may impose important limitations, including incomplete probe annotation and the inability to detect novel transcripts. This may explain why several well-established psoriasis-associated genes (eg, IFNAR1, IFNAR2) were not consistently retrieved under our stringent thresholds (|log2FC| ≥ 2, p < 0.05). When applying more permissive criteria (|log2FC| ≥ 1.0, FDR < 0.05), several of these immune mediators reappeared, indicating that their absence in the strict analysis reflects statistical cutoffs rather than biological irrelevance. This reinforces the importance of integrating multiple datasets and thresholds when interpreting transcriptomic data.