Identification of key apoptosis-related genes and immune infiltration in the pathogenesis of psoriasis

Background Psoriasis is a condition in which skin cells build up and form itchy scales and dry patches. It is also considered a common lifelong disease with an unclear pathogenesis. Furthermore, an effective cure for psoriasis is still unavailable. Reductive apoptosis of keratinocytes and immune infiltration are common in psoriasis. This study aimed to explore underlying functions of key apoptosis-related genes and the characteristics of immune infiltration in psoriasis. We used GSE13355 and GSE30999 to screen differentially expressed apoptosis related genes (DEARGs) in our study. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and gene set enrichment analysis (GSEA) were performed using clusterProfiler package. Protein–protein interaction (PPI) network was constructed to acquire key DEARGs. Transcription factor (TF)–target and miRNA–mRNA network analyses, drug sensitivity prediction, and immune infiltration were applied. Key DEARGs were validated using real-time quantitative PCR (RT-qPCR). Results We identified 482 and 32 DEARGs from GSE13355 and GSE30999, respectively. GO analysis showed that DEARGs were commonly enriched in cell chemotaxis, receptor ligand activity, and signaling receptor activator activity. KEGG pathway analysis indicated that viral protein interaction with cytokine and cytokine receptor was maximally enriched pathway. The GSEA analysis of GSE13355 and GSE30999 demonstrated a high consistency degree of enriched pathways. Thirteen key DEARGs with upregulation were obtained in the PPI network. Eleven key DEARGs were confirmed using RT-qPCR. Additionally, 5 TFs and 553 miRNAs were acquired, and three novel drugs were predicted. Moreover, Dendritic.cells.activated exhibited high levels of immune infiltration while Mast.cells.resting showed low levels of immune infiltration in psoriasis groups. Conclusion Results of this study may reveal some insights into the underlying molecular mechanism of psoriasis and provide novel targeted drugs. Supplementary Information The online version contains supplementary material available at 10.1186/s41065-022-00233-0.

psoriasis is unclear, abnormal proliferation of keratinocytes and disruption of the immune system are central to its pathogenesis [8]. Therefore, exploring gene functions and underlying characteristics of immune infiltration in psoriasis is necessary.
Keratinocytes manifest abnormal proliferation and a marked inhibition of apoptosis [9]. The inhibition of apoptosis is related to imbalances in epidermal homeostasis and results in psoriatic hyperplasia [9,10]. Thus, apoptosis-related genes (APRs) also play a vital role in the pathogenesis. Although mine differential genes have been extensively investigated with the progress of bioinformatics analysis, studies on the function of APRs are limited [11,12]. For instance, Gao et al. revealed seven hub genes, namely, HERC6, ISG15, MX1, RSAD2, OAS2, OASL, and OAS3. Choudhary et al. reported the top 10 hub genes (CCNB1, CCNA2, CDK1, IL1B, CXCL8, MKI67, ESR1, UBE2C, STAT1, and STAT3) but they failed to distinguish between their properties.
Immune system imbalance plays an important role in the formation of psoriatic lesions [13]. Keratinocytes and innate immune cells, such as macrophages, plasmacytoids, and dendritic cells, are activated to secrete inflammatory factors (TNF-α, IFN-γ, etc.), which activate myeloid dendritic cells and then migrate to lymph nodes, under the stimulation of external factors [13,14]. These dendritic cells can in turn activate T lymphocytes. Activated T cells secrete a variety of cytokines that interact with keratinocytes, neutrophils, and macrophages to induce a local persistent inflammatory response, which finally leads to the formation of psoriatic lesions [15].
In this study, we screened Differentially Expressed Apoptosis-Related Genes (DEARGs) associated with psoriasis and explored functions of DEARGs through a comprehensive bioinformatics analysis.

Data selection and DEG screening
We first illustrated the flowchart of current study (Fig. 1), and then collated data according to the GEO data platform ( Table 1). The comparison of psoriasis and control groups showed that the sample size is relatively balanced, thereby indicating the basis of statistical analysis. We obtained 297 up-regulated and 339 down-regulated genes using R software after normalizing the gene expression matrix of GSE13355, as shown in the volcano diagram in Fig. 2a. Additionally the top 10 genes with maximal significant differences were then annotated. The DEG heatmap of GSE13355 shows NN and PP represent control group and psoriasis group respectively (Fig. 2b). Similarly, 38 up-regulated and 3 down-regulated genes were obtained after normalizing GSE30999, as shown in the volcano diagram in Fig. 2c. The DEG heatmap of  GSE30999 shows NL and LS represent control group and psoriasis group respectively (Fig. 2d).

Identification and functional analysis of DEARGs
The ARG list was downloaded from the GeneCards database, and then DEARGs were screened from DEGs of GSE13355 and GSE30999. The results showed that 482 and 32 DEARGs existed in GSE13355 (Fig. 3a) and GSE30999 (Fig. 3b) respectively. GO (Table 2) and KEGG (Table 3) functional enrichment analyses were performed on DEARGs in the two datasets respectively.
The results of GO showed that DEARGs in GSE13355 were mainly related to leukocyte chemotaxis, cell chemotaxis, leukocyte chemotaxis, and signaling receptor activator activity (Fig. 4a). Meanwhile, the results of KEGG demonstrated that DEARGs in GSE13355 were mainly enriched in viral protein interaction with cytokine and cytokine receptor, cytokine − cytokine receptor interaction, and chemokine signaling pathway (Fig. 4b). Similarly, the GO enrichment of GSE30999 showed that DEARGs was generally enriched in cell chemotaxis, receptor ligand activity, and signaling receptor activator activity (Fig. 4c). The KEGG enrichment of GSE30999 demonstrated that DEARGs were The differential analysis of GSE13355 and GSE30999. a Volcano plot of GSE13355 (Red, green and blue represent up-regulated, down-regulated and no differential genes, respectively). b Heatmap of GSE13355 (Blue indicates psoriasis group, and red indicates control group). c Volcano plot of GSE30999(Red, green and blue represent up-regulated, down-regulated and no differential genes, respectively). d Heatmap of GSE30999 (Blue indicates psoriasis group, and red indicates control group). NN, normal skin from controls; PP, involved skin from psoriatic patients; NL, non-lesional skin; LS, psoriasis lesions; NS, no significance mainly enriched in viral protein interaction with cytokine and cytokine receptor, rheumatoid arthritis, and IL − 17 signaling pathway (Figs. 4d and 5).

Gene set enrichment analysis (GSEA)
GSEA was conducted for the two datasets. "c2.cp.kegg. v7.0.symbols.gmt" was selected as the reference gene set and FDR was set to < 0.25, with p < 0.05 indicating the significant enrichment pathway ( Table 4). The GSEA results showed that GSE13355 mainly existed in up-regulated pathways of allograft_rejection, E2F_targets, and IL6_JAK_STAT3_signaling ( Fig. 5a) and was also significantly enriched in down-regulated pathways of androgen_response, bile_acid_metabolism, epithe-lial_mesenchymal_transition (Fig. 5b). In the same way, the GSEA results of GSE30999 indicated a high similarity with the results of GSE13355 ( Fig. 5c and d).
TF-target, miRNA-mRNA network analysis and drug sensitivity prediction JUN, ATF4, CEBPD, NFKB1 and RELA were obtained as TFs, with corresponding targets of FOSL1, CHAC1, CCL20, CXCL8, CXCL1, and TNIP3 through the TRRUST database for TF prediction of 19 DEARGs. Their regulatory relationships were visualized as a regulatory network chart (Fig. 7a). The IC50 of 138 drugs was then predicted using the ridge regression model, and the three final classes of drugs with p < 0.05 were A.770041, GNF.2, and WO2009093972 (Fig. 7b). Finally, the miRNA of DEARGs was predicted using the TargetScan database. A total of 553 miRNAs were predicted to exhibit regulatory relationships with 18 DEARGs. A network visualization graph was established according to the regulatory relationships (Fig. 7c).

Relationship of immune infiltration with psoriasis and control samples
Immune infiltration analysis was performed between psoriasis and control samples in GSE13355 and GSE30999 datasets on the basis of CIBERSORT algorithm. An immune infiltration heatmap was illustrated for GSE13355. Only Macrophages.M2, Mast.cells.resting, T.cells.gamma.delta and other cells expressed in    (Fig. 8c). The combined results of the comparison between groups presented high levels of infiltration of B.cells.memory, T.cells.CD4.memory.resting, NK.cells.activated, and Dendritic.cells.resting in control groups as well as high levels of infiltration of T.cells. follicular.helper, Dendritic.cells.activated, and Mast.cells. resting in psoriasis groups (Fig. 8d).
Overall, Dendritic.cells.activated demonstrated high levels of immune infiltration and Mast.cells.resting exhibited low levels of immune infiltration in psoriasis groups.

Discussion
An effective cure for psoriasis is still unavailable [4] and its pathogenesis remains poorly defined given that multiple factors come into play [16]. The reductive apoptosis of keratinocytes is a common phenomenon in psoriatic lesions [9,17]. Many APRs are involved in the pathogenesis of psoriasis [18]. Hence, exploring the molecular mechanism of APRs is necessary to achieve novel therapeutic targets. We identified and further explored the functions of 482 and 32 DEARGs from GSE13355 and GSE30999 respectively, to predict three novel targeted drugs. Meanwhile, characteristics of immune infiltration in psoriasis were comprehensively analyzed.
A total of 514 DEARGs were screened from the two datasets of GSE13355 and GSE30999. GO annotation and KEGG pathway analyses of genes were performed to investigate their further functions. The results of GO annotation showed that DEARGs were typically enriched in cell chemotaxis, receptor ligand activity, and signaling receptor activator activity in both datasets. The results of KEGG pathway indicated that the viral protein interaction with cytokine and cytokine receptor was the pathway with maximum enrichment. These results are different from the findings of previous studies [11,12,19]. For example, Choudhary et al. reported that keratinocyte differentiation and positive regulation of cytokine production were the respective biological process and molecular function with maximum enrichment and the cytokinescytokine receptor was the pathways with maximum enrichment. GSEA analysis for GSE13355 and GSE30999 presented that the high degree of consistency between the results of the two datasets. Thus, the consistency of GSEA results verified the validity of screened DEARGs. Thirteen key DEARGs with upregulation were derived in the PPI network, including CXCL8, CCL20, CXCL1, CXCL13, S100A12, GZMB, IL19, ATP12A, FOSL1, HYAL4, RHCG, SERPINB4, and TCN1. Eleven key DEARGs, namely, CXCL8, CCL20, CXCL1, CXCL13, S100A12, GZMB, IL19, ATP12A, HYAL4, SERPINB4, and TCN1, were confirmed via RT-qPCR. Among these genes, CXCL8, CCL20, CXCL1, CXCL13, and S100A12 have been extensively investigated [19][20][21]. Meanwhile, studies on GZMB, IL19 and SERPINB4 are lacking. To the best of our knowledge, ATP12A, HYAL4 and TCN1 in psoriasis remain unverified. Based on the literatures, ATP12A is the nongastric form of H + /K + -ATPase and plays an important role in respiratory diseases [22]. HYAL4 is a member of the hyaluronidases (HYAL) family, which is named as such due to their ability to degrade hyaluronan. HYAL4 is produced by mast cells and plays a particular role in maintaining α-granule homeostasis [23]. Transcobalamin (TCN1) is a vitamin B12-binding protein usually highly expressed in tumor tissues and linked to aggressive tumor behavior and poor prognosis [24].
Nineteen DEARGs were used to construct TF-target and miRNA-mRNA network. The action mechanism of upstream transcription factors and downstream miRNA was easily mined when network relationships of 19 DEARGs were built. This complex network relationships also indicated numerous genes involved in the pathogenesis of psoriasis. Three novel targeted drugs, namely, A.770041, GNF.2, and WO2009093972, were  [25]. Some studies have shown that A.770041 can function similar to cyclosporin A to prevent acute transplant rejection [25,26].
GNF-2, an allosteric inhibitor of Bcr-Abl, is a new anticancer drug to treat resistant chronic myelogenous leukemia [27]. At present, the novel drug WO2009093972 still remains unverified. Therefore, A.770041, GNF.2, and WO2009093972 are possible targeted drugs for psoriasis that require further exploration. Furthermore, characteristics of immune infiltration in psoriasis were analyzed. Various immunocytes involved in psoriasis form a complex network centered on T lymphocytes [15]. Compared with normal control groups, our results showed that high levels of immune infiltration of Dendritic.cells.activated and low levels of immune infiltration of Mast.cells.resting in psoriasis groups. Studies have indicated that dendritic cells are the driver of psoriasis that trigger a series of inflammatory responses [28]. Dendritic cells are the source of cytokines (TNF-α, IFN-γ, etc.) that play a crucial role in the pathogenesis of psoriasis [28,29].
However, this study presents the following limitations. First, some datasets such as GSE14905, GSE34248, GSE41745, and GSE40033, were not used  for analysis due to different sample types or platforms. Second, the sample size in this study was insufficient. Future investigations can attempt to integrate multiple databases to increase the sample size. Third, we lack relevant clinical studies and cannot combine clinical information for analysis. Fourth, although experiments with skin tissue samples were not performed in this study, future explorations should include this, to improve the reliability of the results. Therefore, further studies with additional samples and experiments are required to illustrate the role of key DEARGs and underlying mechanism of psoriasis.
The results of our study may reveal some insights into underlying molecular mechanisms of psoriasis and provide novel targeted drugs.

Conclusions
This study explored functions of key APRs and analyzed underlying characteristics of immune infiltration of psoriasis through a comprehensive bioinformatics analysis. We identified 13 key DEARGs associated with psoriasis via GEO datasets analysis. Moreover, dendritic cells play an important role in the initiation of psoriasis. Notably, A.770041, GNF.2, and WO2009093972 are possible novel   The RT-qPCR results of 13 DEARGs. The transcription levels of CXCL8, CCL20, CXCL1, CXCL13, S100A12, GZMB, IL19, ATP12A, HYAL4, SERPINB4, and TCN1 were significantly up-regulated in M5 groups; however, the transcription levels of FOSL1 and RHCG were down-regulated in M5 groups. (Red represents M5 groups, blue represents normal control groups, *, p < 0.05; **, p < 0.01; ***, p < 0.001). DEARGs, Differentially Expressed Apoptosis Related Genes targeted drugs for psoriasis in the future. However, additional experiments are needed to support further our findings.

Data download and preprocessing
Original data of GSE13355 [30] and GSE30999 [31] were downloaded from GEO (https:// www. ncbi. nlm. nih. gov/ geo/) using the GeoQuery package [32]. All samples in the two datasets were from homo sapiens based on GPL570 ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array). The GSE13355 dataset contained 58 lesion skin samples from psoriatic patients and 64 normal skin samples from healthy controls, while the GSE30999 dataset included 85 paired lesion and non-lesion skin samples from psoriatic patients, all of which were included in this study. Original data of GSE13355 and GSE30999 were read using the affy package [33] to obtain their gene expression matrices. Ethical approval is not necessary because this study does not contain any studies with human participants or animals performed by any of the authors.

Screening and functional analysis of differentially expressed apoptosis-related genes (DEARGs)
The limma package [34] was used to screen differentially expressed genes (DEGs) of GSE13355 and GSE30999. The ggplot2 and pheatmap packages were utilized to illustrate a volcano plot and heatmap of DEGs with a cut-off value setting of adj.p value < 0.05 and |log2FC|> 1. The APR list was downloaded from GeneCards database (http:// www. genec ards. org/) [35], and DEARGs were screened from DEGs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of DEARGs were performed using the clusterProfiler package [36]. The gene set enrichment analysis (GSEA) of gene expression matrix was also conducted using the clusterProfiler package, and "c2.cp.kegg.v7.0.symbols. gmt" was selected as the reference gene set, while a false discovery rate (FDR) < 0.25 and p < 0.05 were considered statistically significantly.

Protein-protein interaction network construction
The VennDiagram package [37] was applied to illustrate the Venn diagram of DEARGs in GSE13355 and GSE30999. The STRING (https:// string-db. org/) database [38] was used to construct the PPI network for common DEARGs in the two datasets, and the NetworkAnalyzer in Cytoscape [39] was utilized to analyze related attributes of nodes in the network.

Analysis of immune infiltration
CIBERSORT algorithm [43] is based on linear support vector regression, deconvolves the transcriptome expression matrix, and thereby estimates the composition and abundance of immunocytes in mixed cells. We uploaded the gene expression matrix to CIBERSORT and filtered out the samples with p < 0.05 to obtain the immune infiltration matrix. A heat map was established using the pheatmap package to show the distribution of 22 immunocytes in each sample. Immune infiltration between different subgroups in the two datasets was illustrated using the ggpubr package and visualized at p < 0.05.

Validation with RT-qPCR
Cultured HaCaT cells were treated with 10 ng/ml of M5 (IL-22, TNF-a, IL-17A, IL-1a, and Oncostatin M) (Pepro-Tech) for 48 h. Untreated and treated cells were regarded as normal control (NC) groups and psoriasis cell model (M5) groups respectively. Total RNA was extracted using TRIpure total RNA extraction reagent (#EP013, ELK Biotechnology, China) and reverse transcribed with EntiLink ™ 1st Strand Cdna Synthesis Kit (#EQ003, ELK Biotechnology, China). The RNA expression was detected according to the manual of the StepOne ™ Real-Time PCR System (Life technologies) using EnTurbo ™ SYBR Green PCR SuperMix (#EQ001, ELK Biotechnology, China). Primer sequences were presented in Supplementary Table 1.

Statistical analysis
R software (version 4.0.0, http://r-proje ct. org/) was used to analyzed the data. T-testing was applied to compare the expression levels of 13 key DEARGs between normal control (NC) groups and psoriasis cell model (M5) groups. The ggplot2 package was utilized to perform the T-testing (*, p < 0.05; **, p < 0.01; ***, and p < 0.001).

Abbreviations
DEARGs: Differentially Expressed Apoptosis Related Genes; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; GSEA: Gene Set