Single-Cell Sequencing Combined with Transcriptome Sequencing to Explore the Molecular Mechanisms Related to Psoriasis

Introduction

Psoriasis, a chronic inflammatory skin disorder,¹ exhibits a variable incidence and prevalence globally, affecting approximately 2% of the world’s population. The epidemiology of psoriasis reveals notable variations across different regions and populations, indicating a complex interplay of genetic, environmental, and lifestyle factors.² The disease can manifest at any age, with a bimodal peak in onset observed in the late teens to early 30s and mid-50s to 60s,³ highlighting distinct age-related patterns. There is a roughly equal gender distribution, although some studies suggest a slight male predominance.⁴ A significant genetic component is evident, as individuals with a family history of psoriasis are at a higher risk, with the likelihood increasing with the number of affected relatives.⁵ Environmental influences, including climate, stress, smoking, alcohol consumption, and certain medications, also play critical roles in the onset and exacerbation of psoriasis.⁴ Geographically, psoriasis is more prevalent in high latitude areas, such as Northern Europe and North America, suggesting that vitamin D levels, climate, and sunlight exposure may influence disease prevalence.⁶

Treatment modalities for psoriasis are multifaceted, encompassing topical therapies, phototherapy,⁷ systemic medications, and advanced biological agents. The overarching objectives of management strategies are symptom alleviation and enhancement of patient quality of life.⁸ Through judicious selection and application of therapeutic interventions, the majority of psoriasis patients achieve substantial disease control and a favorable long-term outlook.⁹

The advent of bioinformatics and molecular biology technologies, notably single-cell sequencing and microarray (gene chip) analysis, has revolutionized the landscape of psoriasis research and treatment.¹⁰ Bioinformatics applications facilitate the high-throughput analysis of genetic data, enabling the identification of key genes implicated in psoriasis pathogenesis. Single-cell sequencing, in particular, offers unprecedented insights into the cellular heterogeneity of psoriatic lesions, identifying novel cellular subtypes and their gene expression profiles.¹¹ Concurrently, gene chip technology allows for the comprehensive analysis of gene expression patterns, further refining our understanding of the molecular mechanisms driving psoriasis. Together, these cutting-edge techniques provide a powerful toolkit for the discovery of potential therapeutic targets, paving the way for the development of precision medicine approaches in the management of psoriasis.¹² Through the integration of bioinformatics analyses with clinical strategies, there is an optimistic prospect for advancing psoriasis treatment, tailored to the genetic and molecular landscape of individual patients, thereby optimizing therapeutic outcomes and patient care.¹³

Single-cell RNA sequencing (scRNA-seq) technology has revolutionized the study of immune cell distributions and their interactions, illuminating the intricate dynamics within disease mechanisms. Its capacity to pinpoint specific cell subsets and delineate their unique responses to therapeutic agents or external factors positions scRNA-seq as a pivotal instrument in the advancement of precision medicine.¹⁴ Unlike bulk RNA sequencing, which aggregates gene expression data from a mixed cell population and presents it as a singular average, scRNA-seq unlocks the ability to investigate transcriptional activities at the individual cell level. Traditional bulk RNA sequencing, while valuable, can sometimes provide an oversimplified view due to the complex mixture of cell types and states within tissues.¹⁵ In contrast, scRNA-seq overcomes these limitations by enabling detailed exploration of cellular diversity, uncovering novel or rare cell types, and providing insights into cellular transitions with unparalleled precision.

Materials and Methods

Data Acquisition from the GEO Repository

The Gene Expression Omnibus (GEO) database, a comprehensive gene expression resource, is curated and sustained by the National Center for Biotechnology Information (NCBI). Researchers can access the dataset GSE151177, which contains 18 samples for single-cell analysis, directly from the GEO public database.

For a broader scope of research, the Series Matrix File of GSE41664 is also available on GEO. Analysis of this file, denoted as GPL570, includes gene expression profiles from a total of 157 patients, categorized into 53 control individuals and 104 patients with the disease under study.

Additionally, the Series Matrix File for GSE30999 can be downloaded from GEO. This file, also labeled as GPL570, contains the expression profiles of 170 patients, evenly divided with 85 individuals in the control group and an equal number in the disease group. This balanced cohort facilitates comparative analysis, crucial for identifying disease-specific gene expression patterns.

Single-Cell Transcriptomic Profiling

Utilizing the Seurat software package, we initiated the analysis by importing the expression profiles. Genes exhibiting low expression were filtered out adhering to the parameters: nFeature_RNA greater than 200 and less than 6000, percent mitochondrial genes less than 10%, and nCount_RNA below 50,000. After standardization and normalization of the data, we conducted principal component analysis (PCA) followed by Uniform Manifold Approximation and Projection (UMAP) to delineate the data structure further. The determination of the optimal number of principal components was facilitated by the Elbow Plot method. Positional relationships between the various cell clusters were elucidated through UMAP analysis. To cluster annotation, we cross-referenced the literature and the CellMarker database, which enabled us to label clusters with cellular subtypes pivotal to the understanding of the disease mechanism.

Cell Communication

CellCall emerges as a sophisticated toolkit designed for deciphering the nuances of intercellular communication networks alongside internal regulatory signals, by integrating both intracellular and intercellular cues. It meticulously compiles data sets on the ligand-receptor-transcription factor (L-R-TF) axis, rooted in the KEGG pathway database. Leveraging prior knowledge on L-R-TF interactions, CellCall correlates the expression of ligands and receptors with the downstream activity of transcription factors (TFs) specific to certain L-R pairs, thereby facilitating the inference of intercellular communication dynamics.

WGCNA Analysis

The Weighted Gene Co-expression Network Analysis (WGCNA) approach enables the identification of cohesively expressed gene modules, thereby elucidating the interplay between gene networks and pivotal genes within these networks. Utilizing the WGCNA-R package, a co-expression network encompassing all genes in the dataset was constructed. The top 5000 genes, distinguished by their variance, were selected for in-depth analysis, setting a soft threshold at 9. This analysis transforms the weighted adjacency matrix into a topological overlap matrix (TOM) to assess network connectivity, employing hierarchical clustering to depict the TOM matrix’s clustering tree structure. Distinct branches and color codlings in the clustering tree delineate diverse gene modules, grouping genes by expression patterns into multiple modules.

Model Construction

Selection of candidate genes was followed by the application of lasso regression to construct predictive correlation models. Each gene’s expression value contributed to a risk score formula for patients, weighted by the estimated regression coefficients from lasso regression analysis. Patients were stratified into low-risk and high-risk groups based on median risk score values, employing the ROC curve to evaluate model prediction accuracy.

GSEA Enrichment Analysis

Gene Set Enrichment Analysis (GSEA) was utilized to discern signaling pathway variations between high and low expression groups. Background gene sets, drawn from version 7.0 of the MsigDB database, facilitated differential pathway expression analysis and the identification of significantly enriched gene sets, ranked by consistency scores (adjusted p value < 0.05). GSEA’s application is crucial in research that integrates disease classification with biological significance.

GSVA (Gene Set Difference Analysis)

Gene Set Variation Analysis (GSVA) offers a non-parametric, unsupervised approach for evaluating transcriptome-wide gene set enrichment, translating gene-level changes into pathway-level insights. By scoring selected gene sets, GSVA assesses potential biological function shifts across samples, utilizing gene sets downloaded from the Molecular Signatures Database.

miRNA Network Construction

MicroRNAs (miRNAs), as regulators of gene expression through mRNA degradation or translation inhibition, were analyzed for their roles in key gene regulation. Key gene-associated miRNAs were identified using the miRcode database, with their networks visualized via Cytoscape software.

Regulatory Network Analysis of Key Genes

The RcisTarget R package was leveraged in this investigation to forecast transcription factor activity, relying on motif analysis. RcisTarget’s computations, grounded in motif identification, enable the derivation of a motif’s normalized enrichment score (NES), which is contingent upon the motif’s prevalence within the database. Beyond the motifs annotated in the source data, this research extrapolated additional annotations through motif similarity and gene sequence comparisons. The initial phase of assessing motif influence on gene expression involved calculating the area under the curve (AUC) for each motif and gene set pair, with NES determinations based on AUC distributions across the gene set.

GWAS Analysis

Utilizing the Gene Atlas database, this study accessed a vast compendium of trait-variant associations, underpinned by data from 452,264 participants in the UK Biobank, spanning 778 phenotypes and 30 million loci. This comprehensive resource facilitates the exploration of genetic correlations across a broad spectrum of traits.

Immune Infiltration Analysis

The CIBERSORT approach, esteemed for its application in discerning immune cell compositions within tissue microenvironments, employs support vector regression for deconvolution of immune cell subtype expression matrices. With 547 unique biomarkers identifying 22 immune cell phenotypes, this methodology enabled the quantification of immune cell distributions within patient samples, offering insights into the immune landscape’s role in disease processes.

Cmap Drug Prediction

The Connectivity Map (CMap), a pioneering gene expression database developed by the Broad Institute, focuses on the interplay between small molecule interventions, gene expression changes, and disease states. Incorporating microarray data for 1309 small molecules across five human cell lines under varied¹ experimental conditions, this study harnessed CMap to identify potential therapeutic compounds based on disease-specific gene expression alterations.

Statistical Analysis

Statistical analyses for this research were conducted utilizing the R programming language (version 4.3), with a significance threshold established at p < 0.05. This rigorous statistical framework underpins the validity of the study’s findings, ensuring the reliability of the conclusions drawn from the data analysis.

Results

Preliminary Processing of Single Cell Expression Profile Data

In the analysis of 18 psoriasis-related tissue samples, expression profiles were carefully evaluated. The analysis retained only those cells meeting specific criteria: nFeature_RNA counts greater than 200 but less than 6000, mitochondrial gene content (percent.mt) under 10%, and nCount_RNA less than 50,000. Through this selection process, expression levels from a total of 21,369 cells were included for subsequent examination (as shown in Figure 1A, B).

Figure 1 (A) and (B): This figure presents the distribution of nFeature_RNA, highlighting the criteria applied to filter cells for further analysis. This figure illustrates the mitochondrial gene content (percent.mt) distribution, showing the threshold used to exclude cells with high mitochondrial content. (C): Displays the ten genes with the highest standard deviation in expression across the psoriasis samples, emphasizing variability within the dataset. (D), (E) and (F): Shows the results of principal component analysis (PCA), depicting how the primary components contribute to cell variability. Demonstrates the application of the Harmony algorithm, used for batch effect removal and integration of the data. Visualizes how PCA dimensions were adjusted after applying batch correction. (G): Presents the UMAP plot, highlighting the identification of 18 distinct cellular subtypes within the psoriasis tissue samples.

The analysis highlighted the ten genes with the highest standard deviations in expression, indicating significant variability within the dataset (Figure 1C). The data underwent a series of preprocessing steps including standardization, homogenization for batch effect removal, principal component analysis (PCA), and data integration using the Harmony algorithm, detailed in Figure 1D–F.

The final step involved using Uniform Manifold Approximation and Projection (UMAP) to visualize the relationships among cell clusters, successfully identifying 18 distinct cellular subtypes within the psoriasis samples (Figure 1G).

Subtype Annotation and Receptor-Ligand Interaction Analysis

In a deeper exploration of cellular heterogeneity within psoriasis samples, this study proceeded to annotate identified cellular subtypes. A comprehensive annotation categorized all cells into 12 distinct cell types: CD161_T_cell, CD4_T_cell, CD8_T_cell, KC_Basale, KC_Corneum, KC_Granulosum, KC_Spinosum, Macrophage, Mature_DC, Melanocyte, NK_cell, and SEMIMature_DC, as illustrated in Figure 2A. Additionally, a bubble chart provided a visual representation of 12 classic cellular markers associated with these cell types (Figure 2B). Notably, the frequency of CD4_T_cells within disease samples was significantly elevated in comparison to control samples, as depicted in Figure 2C.

Figure 2 (A): Displays the categorization of cells into 12 distinct cell types based on the analysis of psoriasis samples. (B): A bubble chart illustrating the expression levels of 12 classical cellular markers across the identified cell types. (C): Highlights the significantly higher frequency of CD4_T_cells in psoriasis disease samples compared to control samples. (D): Shows a bubble diagram detailing the intercellular communication pathways, emphasizing Th17 cell differentiation and Cellular senescence pathways. (E): Depicts ring interaction plots that visualize ligand-receptor interactions, revealing the direction and intensity of communication between CD4+ Treg cells and other cell types.

Focusing on CD4_T_cells, a targeted gene analysis was conducted, identifying 118 genes exhibiting significant differential expression (|avg_log2FC| > 1 and p_val_adj < 0.05), thereby establishing a candidate gene set for further investigation. Leveraging the CellCall (v1.0.7) R package, an in-depth cell communication analysis was performed. This analysis uncovered pronounced interactions between CD4+ Treg cells and other cellular entities across key biological pathways. A bubble diagram highlighted these communicative pathways, with Th17 cell differentiation and Cellular senescence pathways being particularly noteworthy for their activity levels. Additionally, ring interaction plots provided insights into ligand-receptor dynamics, elucidating the directionality and intensity of the intercellular communication network (Figure 2D and E).

WGCNA Analysis

The study expanded its analytical scope by employing Weighted Gene Co-expression Network Analysis (WGCNA) on the dataset GSE30999, aiming to delineate the underlying regulatory networks. Setting the soft threshold β to 9 facilitated the identification of gene modules within the dataset, as depicted in Figure 3A. This process revealed a total of 10 gene modules within GSE30999, designated as black (227 genes), brown (2966 genes), cyan (77 genes), greenyellow (319 genes), grey (293 genes), lightcyan (67 genes), midnight blue (72 genes), purple (799 genes), salmon (83 genes), and tan (97 genes), shown in Figure 3B.

Figure 3 (A): Network Analysis (WGCNA) to identify gene modules in the dataset GSE30999. (B): Depicts the identification of 10 distinct gene modules, with each module represented by a unique color, such as black, brown, cyan, and others. (C): Shows the correlation analysis, highlighting the brown module, which displayed the strongest association with the studied traits. (D): Demonstrates the intersection of the brown module genes with the 118 marker genes of CD4_T_cells, leading to the identification of 14 candidate genes for further analysis.

Subsequent investigations focused on elucidating the associations between these gene modules and specific traits. The analysis highlighted the brown module, which exhibited the most significant correlation with the traits under study (correlation coefficient = 0.91, p-value = 2e−67), as illustrated in Figure 3C. Further exploration involved intersecting the genes from the brown module with the previously identified 118 marker genes of CD4_T_cells, yielding 14 candidate genes. These genes were subsequently prioritized for in-depth analysis, as demonstrated in Figure 3D.

Feature Selection via Lasso Regression and SVM

In this phase of the study, datasets GSE30999 and GSE41664 served distinct roles as the training and validation sets, respectively, for the purpose of candidate gene identification through feature selection techniques. Lasso regression, applied to the training set, successfully pinpointed six genes as significant biomarkers for psoriasis, as illustrated in Figure 4A–C. The derived model equation for calculating the RiskScore is as follows: RiskScore = S100A4 × (−0.1112228) + JAML × (−0.073860598) + TRAF3IP3 × 0.006690637 + MIAT × 0.072600089 + IL7R × 0.11280755 + KLRB1 × 0.226251628.

Figure 4 (A–C): Displays the results of Lasso regression applied to the GSE30999 dataset, identifying six significant genes as biomarkers for psoriasis. Further details the feature selection process that highlighted these six genes as potential biomarkers. (D): Demonstrates the diagnostic performance of the 6-gene model, with an AUC of 0.9767, indicating high accuracy. (E): Shows the validation of the model using the GSE41664 dataset, confirming robustness with an AUC of 0.9459. (F): Illustrates the pathways associated with S100A4, such as the TGF-beta signaling pathway, PPAR signaling pathway, and cGMP-PKG signaling pathway. (G): Displays JAML’s involvement in key pathways, including the IL-17 signaling pathway, NF-kappa B signaling pathway, and RIG-I-like receptor signaling pathway. (H): Summarizes TRAF3IP3’s involvement in key pathways, including Cytosolic DNA-sensing, NF-kappa B, and Toll-like receptor signaling. (I): Shows the pathways linked to MIAT, such as the Cytosolic DNA-sensing pathway, NOD-like receptor signaling pathway, and TNF signaling pathway.

Evaluation of the diagnostic performance of this 6-gene model revealed an excellent area under the ROC curve (AUC) of 0.9767, indicating high predictive accuracy (Figure 4D). Subsequent validation using dataset GSE41664 affirmed the model’s robustness, with an AUC of 0.9459 (Figure 4E), underscoring the stability and reliability of the model across different datasets. We selected these six genes S100A4, JAML, TRAF3IP3, MIAT, IL7R, and KLRB1 as key genes for subsequent research.

Pathway Analysis of Key Genes in Psoriasis Progression

The investigation into the specific signaling pathways associated with the six key genes sheds light on the potential molecular mechanisms influencing psoriasis progression. Utilizing Gene Set Enrichment Analysis (GSEA), significant pathways linked to each of the key genes were identified and highlighted for detailed examination.

S100A4 was found to be involved in the TGF-beta signaling pathway, PPAR signaling pathway, and cGMP-PKG signaling pathway (Figure 4F).

JAML showed enrichment in pathways such as the IL-17 signaling pathway, NF-kappa B signaling pathway, and RIG-I-like receptor signaling pathway (Figure 4G).

TRAF3IP3’s associated pathways include the Cytosolic DNA-sensing pathway, NF-kappa B signaling pathway, and Toll-like receptor signaling pathway (Figure 4H).

MIAT plays a role in the Cytosolic DNA-sensing pathway, NOD-like receptor signaling pathway, and TNF signaling pathway (Figure 4I).

IL7R is linked to the Chemokine signaling pathway, p53 signaling pathway, and T cell receptor signaling pathway among others (Figure 5A).

Figure 5 (A): Illustrates IL7R’s association with key pathways such as Chemokine signaling, p53 signaling, and T cell receptor signaling. (B): Highlights KLRB1’s involvement in IL-17, RIG-I-like receptor, and Toll-like receptor signaling pathways. (C): Shows highly expressed JAML enriched in pathways like allograft rejection, complement, and interferon gamma response. (D): Depicts KLRB1’s enrichment in allograft rejection, inflammatory response, and complement pathways. (E): Demonstrates MIAT’s association with E2F targets, allograft rejection, and MYC targets V2. (F): Displays S100A4’s correlation with epithelial-mesenchymal transition, Notch signaling, and myogenesis pathways.

KLRB1 is implicated in the IL-17 signaling pathway, RIG-I-like receptor signaling pathway, and Toll-like receptor signaling pathway, among others (Figure 5B).

Furthermore, Gene Set Variation Analysis (GSVA) revealed that:

Highly expressed JAML is enriched in allograft rejection, complement, interferon gamma response, and other pathways (Figure 5C).

KLRB1 shows enrichment in allograft rejection, inflammatory response, complement, and additional pathways (Figure 5D).

MIAT is associated with E2F targets, allograft rejection, MYC targets V2, and other pathways (Figure 5E).

S100A4 correlates with epithelial-mesenchymal transition, Notch signaling, myogenesis, and more (Figure 5F).

TRAF3IP3 is involved in allograft rejection, E2F targets, complement, and other pathways (Figure 6A).

Figure 6 (A): Shows TRAF33IP’s role in key pathways like allograft rejection, E2F targets, and complement signaling. (B): Highlights IL17R’s involvement in pathways including allograft rejection, inflammatory response, and E2F targets. (C): Visualizes the miRNA-mRNA interaction network involving 81 miRNAs and 189 interactions, showing the regulatory role of miRNAs in key psoriasis-related genes. (D): Displays the enriched transcription factor motifs, particularly highlighting cisbp__M6294 as a key regulator of the six pivotal genes. (E): Presents the Q-Q plot from GWAS data, identifying significant SNPs associated with psoriasis and their connection to the six key genes.

IL7R’s expression enriches allograft rejection, inflammatory response, E2F targets, and further pathways (Figure 6B).

This suggests that key genes may affect the progression of psoriasis through these pathways.

Key Gene miRNA-mRNA Network Construction and Transcriptional Regulation Analysis

To elucidate the regulatory mechanisms of the six key genes implicated in psoriasis progression, an extensive search was conducted within the miRcode database to identify associated non-coding RNA networks. This exploration resulted in the prediction of 81 miRNAs and 189 miRNA-mRNA interaction pairs. The intricate network of these interactions was subsequently visualized using Cytoscape software (Figure 6C), highlighting the complex interplay between miRNAs and key genes.

Focusing on these six pivotal genes, it was observed that they are subject to regulation by a set of common mechanisms, notably involving multiple transcription factors. To delve deeper into the transcriptional regulation governing these genes, an enrichment analysis of the transcription factors was carried out. This analysis employed cumulative recovery curves and Motif-TF annotation techniques. The analysis pinpointed a motif, cisbp__M6294, as having the highest normalized enrichment score (NES: 7.08), underscoring its significance in the transcriptional regulation of the key genes.

The comprehensive display of all enriched motifs and their corresponding transcription factors (Figure 6D) offers a detailed overview of the transcriptional landscape affecting the key genes.

GWAS Analysis of Psoriasis and Key Gene Loci

To delineate the genetic underpinnings of psoriasis, this study ventured into the analysis of Genome-Wide Association Study (GWAS) data, focusing on the identification of causative regions associated with the six key genes. A Q-Q plot was employed to visualize the significant single nucleotide polymorphisms (SNPs) related to the disease, as revealed by the GWAS data (Figure 6E).

The precise mapping of GWAS data facilitated the delineation of key SNP sites located within regions enriched for disease association. This analysis unveiled the SNP pathogenic regions specifically linked to the genes S100A4, TRAF3IP3, MIAT, IL7R, and KLRB1:

TRAF3IP3 and S100A4 were pinpointed within pathogenic regions on chromosome 1.

IL7R was identified in a pathogenic region on chromosome 5.

MIAT and KLRB1 were found on chromosome 12.

The identification of significant SNP sites associated with these five genes within their respective pathogenic regions provides compelling evidence of their involvement in psoriasis. These findings, cataloged in the supplementary Table 1, underscore the genetic predisposition factors contributing to the disease’s pathology.

Analysis of Immune Infiltration in Psoriasis

The tumor microenvironment, a complex milieu comprising immune cells, extracellular matrix components, growth factors, inflammatory cytokines, and distinct physical-chemical properties, plays a pivotal role in disease diagnosis and the efficacy of clinical interventions. This study delves into the interplay between key genes and immune cell infiltration within the psoriasis dataset to elucidate the potential molecular mechanisms through which these genes influence psoriasis progression. It highlights the diversity in immune cell composition across patients and examines the correlation between different immune cell types (Figure 7A and B).

Figure 7 (A) and (B): Depicts the diversity in immune cell composition across patients within the psoriasis dataset. Highlights the correlation between different immune cell types in the psoriasis microenvironment. (C): Shows the increased presence of Macrophages M1 and T follicular helper cells in psoriasis patients compared to controls, indicating a unique immune landscape in psoriatic pathology.

Comparative analysis revealed that, in the disease cohort, there is a notable increase in the presence of Macrophages M1 and T cells follicular helper compared to controls (Figure 7C), suggesting a distinct immune landscape in psoriatic pathology. Further investigation into the associations between key genes and immune cell populations demonstrated significant correlations:

IL7R exhibits a strong positive correlation with T cells gamma delta, among others, and a notable negative correlation with activated NK cells, etc. (Figure 8A).

Figure 8 (A): IL7R shows a strong positive correlation with T cells gamma delta and a negative correlation with activated NK cells. (B): JAML is significantly positively correlated with T cells gamma delta and negatively with activated NK cells. (C): KLRB1 exhibits a positive association with T cells gamma delta and a negative correlation with activated NK cells. (D): MIAT has a positive correlation with activated CD4 memory T cells and a negative correlation with resting Mast cells. (E): S100A4 is positively correlated with resting Mast cells and negatively with activated CD4 memory T cells. (F): TRAF3IP3 shows a positive correlation with activated CD4 memory T cells and activated NK cells.

JAML is significantly positively correlated with T cells gamma delta, among others, and negatively with activated NK cells (Figure 8B).

KLRB1 shows a significant positive association with T cells gamma delta, among others, and a negative correlation with activated NK cells (Figure 8C).

MIAT has a significant positive correlation with activated CD4 memory T cells, etc., and a negative correlation with resting Mast cells, etc. (Figure 8D).

S100A4 is positively correlated with resting Mast cells, etc., and negatively with activated CD4 memory T cells, etc. (Figure 8E).

TRAF3IP3 displays a significant positive correlation with activated CD4 memory T cells, etc., and with activated NK cells (Figure 8F).

Drug Target Prediction Using Cmap

Leveraging the limma package, this investigation quantitatively assessed the differential gene expression within the GSE30999 dataset to pinpoint genes variably expressed in psoriasis. The criteria for differential gene identification were stringent, with a significance threshold set at a P-value of less than 0.05. This analysis culminated in the identification of 3622 differential genes, comprising 1775 up-regulated and 1847 down-regulated genes.

The differentially expressed mRNAs were categorized into two groups of the top 150 for each expression polarity. Utilizing the Connectivity Map (CMap) database, this study embarked on predicting potential drug targets among these differential genes. The analysis revealed that the expression profiles influenced by certain drugs, including methotrexate, budesonide, aminopurvalanol-a, and selumetinib, exhibited negative correlations with the disease-perturbed expression profiles. Notably, the implicated drugs were identified as candidates capable of ameliorating or potentially reversing the pathological state associated with psoriasis (Figure 9A–D).

Figure 9 (A–D): Illustrate the Connectivity Map (CMap) analysis results, revealing drugs like methotrexate, budesonide, aminopurvalanol-a, and selumetinib, which exhibit negative correlations with disease-altered expression profiles, suggesting their potential as therapeutic candidates for psoriasis. (E) and (F): Presents the single-cell expression analysis of key genes such as S100A4, JAML, TRAF3IP3, MIAT, IL7R, and KLRB1 across various cell types, including CD4_T_cells, KC_Basale, macrophages, and NK cells. Provides detailed profiling of the expression patterns of the same key genes (S100A4, JAML, TRAF3IP3, MIAT, IL7R, KLRB1), focusing on their differential distribution in specific keratinocyte layers (KC_Corneum, KC_Granulosum) and immune cells. (G): Explores the connection between key gene expression and immune/metabolic pathways, showing high expression levels of genes in pathways like allograft rejection, oxidative phosphorylation, and reactive oxygen species signaling.

Single-Cell Expression Profiles of Key Genes and Their Association with Immune/Metabolic Pathways

This study embarked on an intricate analysis of key gene expression within single-cell data to uncover the expression patterns of S100A4, JAML, TRAF3IP3, MIAT, IL7R, and KLRB1 across various cell types. These types include CD161_T_cell, CD4_T_cell, CD8_T_cell, KC_Basale, KC_Corneum, KC_Granulosum, KC_Spinosum, Macrophage, Mature_DC, Melanocyte, NK_cell, and SEMIMature_DC (Figure 9E and F). This comprehensive profiling elucidates the differential expression of these genes across the cellular landscape of psoriasis.

Further, the study probed into the correlation between the expression of these key genes and immune/metabolic pathways. It was found that S100A4, JAML, TRAF3IP3, MIAT, IL7R, and KLRB1 are prominently expressed in the immune pathway related to allograft rejection. In the metabolic pathway, oxidative phosphorylation showed high expression levels, while the reactive oxygen species pathway was notably expressed in the signaling pathways domain (Figure 9G).

Additionally, leveraging data from the GeneCards database (https://www.genecards.org/), the study identified co-expression patterns between immune-related genes (STAT5B, RAG1) and the six key genes, providing a nuanced understanding of gene interaction within immune contexts. This co-expression network was visualized across 12 cell markers, offering insights into the complex regulatory mechanisms at play (Figure 10A and B).

Figure 10 (A) and (B): Visualizes the co-expression network between immune-related genes (eg, STAT5B, RAG1) and the six key genes, using data from the GeneCards database, highlighting interactions across 12 cell markers.

Discussion

The epidemiology of psoriasis, with a prevalence of approximately 2% in the global population,¹⁶ underscores the significant public health burden posed by this chronic inflammatory skin condition.¹⁷ The variable incidence across different geographic regions and populations points to the interplay of genetic, environmental, and lifestyle factors in its pathogenesis. The management and prognosis of psoriasis hinge on a comprehensive understanding of its epidemiological landscape, facilitating the development of targeted therapies and personalized care plans.^18–22 The identification of key genes associated with psoriasis is paramount in this context, offering insights into the disease’s molecular underpinnings and opening avenues for precision medicine approaches.²³

Current research into immune infiltration in psoriasis reveals a complex interplay between various immune cell types and the psoriatic microenvironment.²⁴ The relationship between immune infiltration and the pathophysiology of psoriasis is a subject of intensive study, highlighting the roles of T cells, dendritic cells, and other immune cells in perpetuating inflammatory responses. This research area is crucial for unraveling the mechanisms that drive psoriasis and for identifying potential targets for therapeutic intervention.^25–31

The investigation of key genes in psoriasis and their involvement in specific signaling pathways provides valuable insights into the disease’s molecular mechanisms. Studies have identified several critical genes and elucidated their roles in pathways such as IL-17 signaling,³² TNF signaling, and others that are instrumental in the onset and progression of psoriasis. Understanding these pathways not only sheds light on the disease’s pathogenesis but also identifies potential molecular targets for treatment.³³

Non-coding RNAs associated with key psoriasis genes represent an emerging research frontier, with evidence suggesting they play pivotal roles in regulating gene expression.^34–37 Whether these non-coding RNAs act to inhibit or activate core gene expression in psoriasis is a question of significant interest. The potential regulatory functions of non-coding RNAs on psoriasis-related genes offer a deeper understanding of gene expression control mechanisms in psoriasis,³⁸ presenting novel therapeutic targets and biomarkers for disease management and prognosis.³⁹

In summary, the multifaceted approach to psoriasis research—from epidemiology and immune infiltration studies to molecular investigations of key genes and non-coding RNAs—highlights the complexity of the disease and the necessity for integrated research strategies.^40–44 These efforts converge on the ultimate goal of improving psoriasis management and patient outcomes through precision medicine, targeted therapies, and a more profound understanding of the disease’s molecular basis.⁴⁵

It is essential to acknowledge the limitations inherent in our current dataset,⁴⁶ primarily sourced from public databases rather than direct patient data. This approach, while valuable for preliminary investigations, may not capture the full spectrum of genetic and environmental variability present in individual psoriasis cases.⁴⁷ Recognizing this constraint, our future research endeavors are poised to incorporate data derived directly from patients within a clinical setting.⁴⁸ By leveraging patient-specific information, we aim to validate our findings and refine our understanding of psoriasis pathogenesis and treatment response.⁴⁶ This step will not only enhance the robustness of our conclusions but also ensure that our research outcomes are directly applicable to patient care, thereby bridging the gap between theoretical research and practical, clinical applications.^48–50

Our Drug Target Prediction analysis revealed that the expression profiles influenced by certain drugs, including methotrexate, budesonide, aminopurvalanol-a, and selumetinib, exhibited negative correlations with the disease-perturbed expression profiles. It was found that S100A4, JAML, TRAF3IP3, MIAT, IL7R, and KLRB1 are prominently expressed in the immune pathway related to allograft rejection. In the metabolic pathway, oxidative phosphorylation showed high expression levels, while the reactive oxygen species pathway was notably expressed in the signaling pathways domain. Our study reveals the potential medicine and pathogenesis of psoriasis. These potential pathway and hub genes may provide new ideas for further mechanism research.