Integrated Bioinformatic Analysis of Key Biomarkers and Signaling Pathways in Psoriasis


 Background: Psoriasis is a relatively common autoimmune inflammatory skin disease with a chronic etiology. The present study was designed to detect novel biomarkers and pathways associated with psoriasis incidence. Methods: Differentially expressed genes (DEGs) associated with psoriasis in the Gene Expression Omnibus (GEO) database were identified, and their functional roles and interactions were then annotated and evaluated through GO, KEGG, and gene set variation (GSVA) analyses. In addition, the STRING database was leveraged to construct a protein-protein interaction (PPI) network, and key hub genes from this network were validated as being relevant through receiver operating characteristic (ROC) curve analyses of three additional GEO datasets. The CIBERSORT database was additionally used to assess the relationship between these gene expression-related findings and immune cell infiltration. Results: In total 197 psoriasis-related DEGs were identified and found to primarily be associated with the NOD-like receptor, IL-17, and cytokine-cytokine receptor interaction signaling pathways. GSVA revealed significant differences between normal and lesional groups (P < 0.05), while PPI network analyses identified CXCL10 as the hub gene with the highest degree value, whereas IRF7, IFIT3, OAS1, GBP1, and ISG15 were promising candidate genes for the therapeutic treatment of psoriasis. ROC analyses confirmed that these 6 hub genes exhibited good diagnostic efficacy (AUC > 70%), and were predicted to be associated with increased sensitivity to 10 drugs (P < 0.01). The CIBERSORT database further predicted that these hub genes were associated with infiltration by 22 different immune cell types. Conclusion: These results offer a robust foundation for future studies of the molecular basis for psoriasis, potentially guiding efforts to treat this common and disruptive disease.


Introduction
Psoriasis is a chronic autoimmune-mediated in ammatory skin condition that arises as a consequence of both environmental and genetic factors, leading affected patients to suffer from scaly erythema [1]. An estimated 2% of people in Europe and 0.47% of people in China are thought to suffer from psoriasis [2], and many of these patients suffer from poor quality of life [3]. While psoriasis is not a life-threatening condition, it can impose signi cant psychological and nancial burdens on affected individuals. A number of therapies that can effectively treat psoriasis have been developed in recent years [4][5], including systemic retinoic acid delivery, immunosuppressant or cyclosporine treatment, and external glucocorticoid or retinoic acid application. However, some patients experience signi cant adverse reactions to these treatments. Owing to its complex and multifactorial etiology, psoriasis remains impossible to cure. It is thus vital that the molecular mechanism governing psoriasis development be clari ed in order to guide more effective patient treatment.
Psoriasis exhibits polygenic patterns of inheritance and is primarily driven by T cells, with a range of factors stimulating excessive autoimmune activity in affected individuals. Recent studies of the immunological basis for psoriasis have led to the development and commercialization of novel biological therapeutics targeting TNF-α, PDE-4, T cells, small cell signal transduction molecules, and VEGF. These inhibitors exhibit excellent clinical e cacy, but they are extremely expensive and have the potential to increase the risk of certain infections in treated patients, making them imperfect therapeutic agents.
Several different chemokine signaling pathways have been linked to the incidence of psoriasis [6][7].
These include several CXC family chemokines that can facilitate in ammatory T cell responses (CXCL10, CXCL16) and neutrophil responses (CXCL1, CXCL5) in psoriatic contexts, with these factors signaling through receptors including CXCR1 and CXCR2 in T cells, and CXCR3 and CXC6 in neutrophils [8].
Keratinocytes and other cells in the dermis of plaque psoriasis patients express CXCL10 and its receptor CXCR3, and CXCL10 levels in patient serum fall following patient treatment. While patients with newly diagnosed early-stage psoriasis often exhibit high CXCL10 levels, these levels tend to be lower in the serum of patients who have been suffering from this disease for longer periods of time [9]. Antibodies targeting CCL2 and CXCL10 can suppress the expression of these chemokines in psoriatic lesions, thereby reducing the severity of psoriasis as measured using PASI scores [10]. More research, however, is needed to fully understand the role of CXCL10 in psoriasis patients.
In recent decades, a range of microarray and sequencing technologies have been used to detect genetic and transcriptomic changes associated with particular disease states. These analyses have supported the identi cation of differentially expressed genes (DEGs) and functional pathways associated with the development and progression of psoriasis. However, individual microarray analyses have the potential to yield false-positive results that cannot be replicated. In addition, microarrays conducted using different platforms may exhibit result discordance owing to small sample sizes and to differences in the underlying platforms and statistical methods employed in individual studies. As such, further comprehensive analyses are necessary to ensure the reliable identi cation of robust diagnostic biomarkers and therapeutic targets capable of aiding in the treatment of psoriasis.
The present study was therefore designed to conduct an integrated analysis of three microarray datasets from the Gene Expression Omnibus (GEO) so as to identify key genes that were differentially expressed when comparing normal and psoriatic skin samples. Additional enrichment and network analyses were then performed to clarify the relationships between these DEGs and the progression of psoriasis.
Together, our data provide a comprehensive framework that highlights the molecular mechanisms governing psoriasis pathogenesis, with CXCL10 in particular being identi ed as a promising biomarker of this debilitating disease.

Data collection and preprocessing
Raw gene expression pro le data (cel les) for the GSE13355 (58 lesional skin samples, 64 normal skin samples), GSE30999 (82 lesional skin samples, 82 normal skin samples), and GSE106992 (125 lesional skin samples, 67 normal skin samples) datasets were downloaded from the GEO database [11]. Hallmark gene sets were downloaded from the Molecular Signatures Database (MSigDB) [12]. Drug sensitivity data were downloaded from the cellminer Database [13].

DEG identi cation
DEGs were identi ed and plotted using R 3.6.2 with the Affy [14] and Limma packages [15]. DEGs were those meeting the following criteria: Adjusted P < 0.05, log |FC| > 1.0.

Principal component analysis
A multivariate principal component analysis (PCA) [16] was used to differentiate between samples. For this analysis, DEGs were treated as variables, and differences between healthy and psoriatic skin samples were identi ed.

Functional enrichment analyses
Gene Ontology (GO) [17] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [18] enrichment analyses were conducted with the R clusterPro ler package [19] in an effort to clarify the functional roles of identi ed DEGs in psoriasis. Enrichment was considered signi cant at a threshold of P < 0.05.
Gene set variation analysis (GSVA) GSVA analyses evaluate pathway and biological process activity over a sample population in an unsupervised manner [20]. The "c2.cp.kegg.v6.2.symbols" and "h.all.v6.2.symbols" gene set les from the Molecular Signatures Database were used to conduct the GSVA analysis with the GSVA R package. P < 0.05 was the signi cance threshold for this analysis.

Protein-protein interaction (PPI) network construction
The STRING database (http://string-db.org)(v 10.0) [21] was used to construct a PPI network incorporating psoriasis-related DEGs with a combined interaction score > 0.4. The network was visualized with the open-source Cytoscape program (v 3.4.0) [22].
The plug-in Molecular Complex Detection (MCODE) (version 1.4.2) of Cytoscape is an APP for clustering a given network based on topology to nd densely connected regions. The PPI networks were drawn using Cytoscape and the most signi cant module in the PPI networks was identi ed using MCODE. The criteria for selection were as follows: MCODE scores > 5, degree cut-off = 2, node score cut-off = 0.2, Max depth = 100 and k-score = 2.

Hub genes selection and analysis
The hub genes were selected with degrees ≥ 10. A network of the genes and their co-expression genes was analyzed using Cytoscape. The CytoHubba plugin in Cytoscape and molecular complex detection (MCODE) were employed to identify hub genes. ceRNA network construction Interactions between lncRNAs and miRNAs were predicted using miRcode [23]. StarBase was used to predict mRNA targets of miRNAs of interest [24], along with miRDB [25], and miRWalk [26]. Those interactions that were predicted by all three databases were the focus of interactions between miRNAs and psoriasis-related DEGs. A lncRNA/miRNA/mRNA ceRNA network was then constructed based on these predicted interactions, with the resultant network being visualized with the Cytoscape software [22].

ROC analysis
Receiver operating characteristic (ROC) curve analyses were used to assess the diagnostic performance of particular DEGs as a means of differentiating between psoriatic and normal samples based on area under the ROC curve (AUC) values [27]. These analyses were conducted with the R pROC package, with the value yielding the maximum sum of sensitivity and speci city serving as a cut-off value[28].

Drug Sensitivity Analysis
Drug sensitivity data were downloaded from GDSC, and the relationships between gene expression and putative psoriasis drug sensitivity were analyzed.

Immune in ltration analysis
The CIBERSORT database was used to evaluate the relationship between genes and immune cell in ltration pro les in normal and psoriatic samples in these three datasets [29]. The resultant data were visualized using bar charts, correlation charts, violin plots, and heatmaps. Correlations between hub genes and in ltrating immune cell populations in these three datasets were assessed, and samples were strati ed into those that exhibit low or high levels of a given gene of interest based on the median expression value for that gene. Box plots were then used to visualize differences in the degree of immune cell in ltration by particular cell types between these high and low expression sample cohorts.

Statistical analysis
Data are means ± standard deviation, and were analyzed using R (v.3.6.2), with P < 0.05 as the signi cance threshold.

Principal component analysis and distribution density mapping of psoriasis-related differentially expressed genes
We began by plotting distribution density maps for sample gene expression data from the GSE13355, GSE30999, and GSE106992 datasets ( Fig. 1A-C), revealing similar expression pro les for both control and lesional samples in all cases, suggesting that su cient data were available to permit variance analysis. Subsequent PCA revealed that control and lesional samples from these three datasets could be separated by the chosen DEGs, indicating that screened DEG expression patterns were speci c and associated with psoriasis ( Fig. 1D-F).
Differential gene expression analysis Next, we identi ed DEGs in these three datasets using the R limma package, including 655, 849, and 286 DEGs in the GSE13355, GSE30999, and GSE106992 datasets, respectively. The resultant data were then subjected to heatmap and clustering analyses, which revealed clear and reliable differentiation between the control and psoriasis samples ( Fig. 2A, C, E). Volcano plots additionally highlighted differences between groups with respect to DEG expression (Fig. 2B, D, F).

Functional enrichment analyses
Next, the 197 DEGs that were shared across these three datasets were next identi ed and used to conduct functional enrichment analyses (Fig. 3). GO terms signi cantly enriched for these genes included BPs such as response to virus, defense response to virus, antimicrobial humoral response, type 1 interferon signaling pathway, and cellular response to type 1 interferon (Table 1A). Enriched CCs included corni ed envelope, speci c granule lumen, tertiary granule lumen, secretory granule lumen, and cytoplasmic vesicle lumen (Table 1B). Enriched MFs included chemokine receptor binding, serine-type endopeptidase activity, serine-type peptidase activity, serine hydrolase activity, and chemokine activity (Table 1C). KEGG pathways revealed these DEGs to be enriched in 5 pathways, including the IL-17 signaling pathway (11 genes, hsa04657), the NOD-like receptor signaling pathway (14 genes, hsa04621), the Viral protein interaction with cytokine and cytokine receptor pathway (10genes, hsa04061), the in uenza A pathway(11genes, hsa05164), and the cytokine-cytokine receptor interaction pathway (13 genes, hsa04060) (Table 1D). While the speci c genes enriched in these pathways for each dataset were not identical, all three datasets were nonetheless enriched in the same pathway (Table 2).

GSVA analysis
The hallmark gene set was next selected as a reference set, after which a GSVA approach was used to assess differences in pathway enrichment between control and psoriatic samples in these three datasets. This analysis revealed that there were differences in the relative degree of pathway enrichment for different datasets, although the overall trends remained the same as above (Fig. 4).

Protein-protein interaction network
Next, interactions among these 197 DEGs were assessed by generating a PPI network with the STRING database (Table 3A). An interaction network diagram for these genes was additionally generated (Fig. 5A), and the MCODE plugin was used to select cluster 1 ( Fig. 5B; Table 3B). The cytoHubba tool was then used to identify key hub genes within this cluster (Fig. 5C), leading to the identi cation of 10 key hub genes, with darker coloration corresponding to a higher degree value (Table 3C). A one-step interaction node for these 10 genes was then prepared (Fig. 5D). Those 6 genes that overlapped between the MCODE and cytoHubba hub gene analyses (CXCL10, IRF7, IFIT3, OAS1, GBP1, ISG15) were then identi ed as putative biomarkers of psoriasis (Table 4).

ceRNA network construction
The predictions made by the starBase, miRDB, and miRWalk databases were next integrated to identify putative miRNAs targeting these hub genes. StarBase was then further used to predict lncRNAs interacting with these miRNAs, facilitating the construction of a ceRNA network (Fig. 6).

Drug sensitivity analysis
Next, drug sensitivity analyses were conducted for these six hub genes (Fig. 8). The results of these analyses suggested that certain drugs were negatively and positively correlated with the expression of these hub genes, with positive correlations indicating that a patient is more likely to be sensitive to a given drug if the corresponding gene is upregulated, and negative correlations indicating the opposite. In total, sensitivity to 10 different drugs was predicted to be increased (P < 0.05), including LDK-378, AP-26113, Alectinib, Lenvatinib, Tanespimycin, Pazopanib, Elesclomol, Estramustine, Abiraterone, Idelalisib, and Itraconazole.

Immune cell in ltration analysis
Lastly, we used to CIBERSORT algorithm to assess immune cell in ltration for each of these three datasets, assessing the immune cell composition of each sample, correlations between immune cell populations, and differences in immune cell levels between normal and diseased samples as represented using heatmaps and violin plots (Fig. 9).
We then conducted correlation analyses to examine the relationships between the hub genes identi ed above and immune cell in ltration in these three datasets, revealing signi cant differences in immune cell levels between samples with low and high levels of hub gene expression (Fig. 10).

Discussion
Psoriasis is an autoimmune in ammatory skin disease with a chronic etiology that impacts an estimated 0.96 per 1,000 people in Europe [30]. Psoriatic lesions can manifest as systemic or localized regions of scaly erythema and plaques, primarily in younger adults. There are four primary clinical types of psoriasis: ordinary, arthritis, pustular, and erythroderma, with the ordinary type accounting for over 99% of the total psoriasis disease burden. This disease generally exhibits a recurrent, chronic course that adversely impacts many aspects of individual physical and mental health, interfering with social interactions and other facets of daily life. In one survey, 79% of psoriasis patients reported that the disease had negatively impacted their lives, with 10% having contemplated suicide [31]. Current therapeutic targets associated with psoriasis include TNF-α, PDE-4, T cells, VEGF, and a range of signaling molecules, with many inhibitors targeting these pathways having been developed and marketed to date. Identifying novel biomarkers of psoriasis has the potential to improve patient quality of care and cost of treatment by facilitating the development of highly speci c, inexpensive, and e cient treatments. Herein, we explored psoriasis-related gene expression patterns through a microarray-based bioinformatic analysis of three GEO gene expression pro les, leading to the identi cation of 197 DEGs associated with psoriasis incidence.
Identi ed DEGs were enriched for pathways including NOD-like receptor signaling, IL-17 signaling, cytokine-cytokine receptor interactions, In uenza A, viral protein interaction with cytokine and cytokine receptors, hepatitis C, drug metabolism-other enzymes, and biosynthesis of unsaturated fatty acids pathways. This analysis revealed a close relationship between IL-17 and psoriasis, consistent with prior reports highlighting a key role for this cytokine in psoriatic disease progression. IL-17 mRNA levels are elevated in psoriatic skin samples relative to normal skin [32]. Treatment of keratinocytes with IL-17 can further promote the upregulation of pro-in ammatory cytokines including IL-6 and IL-8, thereby potentially exacerbating psoriasis [33]. The topical application of the TLR7/8 ligand imiquimod can induce IL-17A and IL-17F expression and drive psoriasis-like dermatitis in mice [34]. Treatment with CsA and anti-TNF-a drugs can further decrease IL-17A, IFN-g, IL-23p19, and CCL20 levels in psoriatic lesions while improving psoriatic eruptions, [33.35-37], suggesting that IL-17A and other proin ammatory cytokines play key roles in psoriasis development.
In our PPI network analysis, we identi ed CXCL10 as having the highest degree value, while IRF7, IFIT3, OAS1, GBP1, and ISG15 were identi ed as key targets for the drug-based treatment of psoriasis. Targeting CXCL10 may thus be a viable approach to treating this debilitating disease. The upregulation of CXCL10 has been detected in many Th1 in ammatory diseases in humans [38]. Reductions in serum CXCL10 levels over time have been reported to be related to new-onset psoriatic arthritis in psoriasis patients [39]. We found that CXC family cytokines including CXCL1, CXCL9, and CXCL10 were primarily enriched in skin-and cytokine signaling-related GO and KEGG pathways, suggesting a potential role for these factors as regulators of psoriasis. Luster et al. previously highlighted the relationship between chemokines and in ammation, likely explaining their predicted role in psoriasis patients [40]. IRF7, IFIT3, OAS1, GBP1, and ISG15 were also identi ed as promising targets for the treatment of psoriasis. Raposo et al. detected signi cant upregulation of 16 different antiviral genes in lesional psoriatic skin, with ISG15 expression being increased more than two-fold [41], and others have similarly observed cutaneous upregulation of antiviral proteins in the context of psoriasis [42], potentially explaining why these patients rarely experience cutaneous viral infections. Psoriasis is also associated with keratinocyte expansion, and these cells may express high levels of antiviral genes. Immune cells also express high levels of these antiviral genes [41], potentially indicating robust immune cell and keratinocyte activity in this pathological context. Consistently, analyses of the CIBERSORT database suggested that these 6 hub genes were associated with in ltration by 22 different immune cell types. Together, these results offer new insights into biomarkers and pathways that may be targeted to more effectively treat psoriasis.
There are several limitations to this analysis. For one, integrated analyses of blood samples that were not conducted in the present study are required to more fully clarify the molecular mechanisms governing psoriasis. Secondly, an individual microarray analysis has the potential to yield false-positive results that cannot be replicated. As such, it is essential that multiple individual datasets be integrated together to facilitate a more robust analysis of the genetic basis for psoriasis. Third, additional mechanistic experiments including western blotting, qPCR, and immunohistochemical assays will be required to clarify the functional importance of the hub genes identi ed in this study as regulators of psoriasis development and progression.

Conclusions
In summary, we conducted the present study to identify key genes likely to be involved in psoriasis development or progression This led us to identify 197 key psoriasis-related DEGs, with CXCL10 representing a particularly promising biomarker worthy of future study as a sensitive indicator of psoriasis. However, further research will be needed to understand the biological role of this chemokine in patients with this condition, and to validate the results of this study so as to improve patient care.