Construction of a three commitment points for S phase entry cell cycle model and immune-related ceRNA network to explore novel therapeutic options for psoriasis

1.   Introduction

Psoriasis is a complex chronic inflammatory disease at the genetic, genomic, and cellular levels, with scaling, erythema, and thickening as characteristics of this condition ^[1]. Patients with psoriasis have a reduced quality of life compared with healthy individuals and are at an increased risk of developing depression, cardiovascular diseases, and psoriatic arthritis ^[2]. It takes approximately 50 days for basal keratinocytes to differentiate into anucleate corneocytes in normal skin but merely approximately five days in psoriatic lesions ^[3]. Moreover, the classical view holds that keratinocyte progression and differentiation are mutually exclusive (keratinocytes stop cell cycle progression before differentiation). However, this hypothesis does not explain the growth of keratinocytes during differentiation and the epidermal acanthosis observed in psoriatic lesions ^[4]. Nevertheless, studies have shown that the differentiation of keratinocytes can occur within any cell cycle phase ^[5,6], which is a breakthrough from the previous view.

Although keratinocytes play a pivotal role in psoriasis development, whether the trigger of psoriasis is keratinocytes or immune cells remains unclear. It has also been revealed that keratinocytes and the immune system have complex interactions ^[7]. Furthermore, while the cell cycle dysregulation of keratinocytes leads to keratinocyte proliferation and differentiation, promoting immune cell infiltration via antimicrobial peptides and chemokines, immune cells also stimulate the formation of a cytokine microenvironment that promotes keratinocyte proliferation and differentiation ^[7]. Specifically, researchers have demonstrated the effects of certain cytokines, such as tumor necrosis factor α (TNF-α) and interferon γ (IFN-γ), in cell cycle-related signaling pathways ^[8,9]. In addition, some cell cycle-regulating proteins, for instance, transforming growth factor-beta (TGF-β) ^[7], p53 protein ^[10], MYC protein, and retinoblastoma protein (Rb) ^[11], play crucial roles in regulating keratinocyte proliferation and differentiation. However, these roles require further elucidation.

In the classic cell cycle model, unphosphorylated Rb accompanies newborn cells into the G1 phase. Here, cyclin D-CDK4/6 progressively phosphorylates Rb and partially activates E2F1, E2F2, and E2F3A (E2Fs) due to the action of protein phosphatase 1 (PP1). Owing to this process, increased transcription of CCNE1 and CCNE2 is observed, and cyclin E-CDK2 further phosphorylates Rb to activate E2Fs in late G1 ^[12]. The hyperphosphorylation of Rb persists throughout the cell cycle ^[13]. Although several studies have observed that most newborn cells emerge with hyperphosphorylated Rb based on this cascade, only a few cells are affected by PP1 ^[14,15,16]. Furthermore, it has also been shown that at the end of the G1 phase, cyclin E-CDK2 does not phosphorylate Rb further, and Rb phosphorylation by cyclin D-CDK4/6 is sufficient to maintain the entire cell cycle ^[17]. These studies present a different perspective from the classic cell cycle model. In 2020, Hume et al.^[18] first proposed the CP1–3 (three commitment points for S phase entry) model of the cell cycle, advocating that mitogens compete with DNA damage to regulate whether cells enter the S phase. They divided the specific mechanism into three regulatory points, CP1, CP2 and CP3, mainly controlled by cyclin D and p21, including p53, DNA damage response (DDR), and E2Fs, respectively, providing a new direction for researchers.

Several deep learning- and artificial intelligence-based methods for cell cycle classification, such as label-free cell cycle analysis based on imaging flow cytometry ^[19] and cell image analysis using deep learning methods ^[20], have been recently developed. Alternatively, applying the traversal theory to the cell cycle, cell cycle markers can be used to infer a cell's current phase of the cell cycle, as well as the lengthening/shortening of a particular phase. Liu et al. constructed a computational framework called DLGene based on convolutional neural networks for cell cycle regulatory gene detection ^[21], and Huang et al. used multiple machine learning approaches to identify new potential cell cycle-related genes based on single-cell RNA sequencing (scRNA-seq) data ^[22]. With the development of transcriptomics techniques, various methods to analyze time-course expression profile data and gene regulatory networks have emerged. For instance, Huang et al. developed GeneReg, a time delay linear regression model-based method for analyzing time-delayed gene regulatory networks ^[23]. Moreover, Liu et al. proposed a network inference method, IMBDANET, based on an improved Markov blanket discovery algorithm to analyze gene regulatory networks ^[24]. Kordmahalleh et al. developed a hierarchical recurrent neural network for identifying time-delayed gene interactions ^[25], while Yang et al. proposed a restricted gene expression programming (RGEP)-based parallel algorithm MPRGEP ^[26]. These studies provide practical methods for constructing dynamic gene regulatory networks and simulating interactions in time-delayed gene regulatory networks, which can be applied to cell cycle and differentiation studies. Moreover, to better reconstruct the cell cycle phases from gene expression data, the number of parameters in scRNA-seq can be increased to thousands, and the data can be processed by bioinformatics methods, such as the CellCycleScoring function in the R language "Seurat" package ^[27].

With respect to skin cells and psoriasis, Grabe et al. simulated the pathological conditions of psoriasis by altering transit amplified cells ^[28], and Zhang et al. constructed a model to study epidermal homeostasis under normal and pathological conditions and simulated keratinocyte proliferation and differentiation in psoriasis, predicting homeostatic patterns ^[29]. Moreover, Ohno et al. constructed a model to simulate the dynamic equilibrium between the epidermis and dermis ^[30]. However, it remains unclear how to systematically analyze the cell cycle of skin cells in patients with psoriasis.

Since the competing endogenous RNA (ceRNA) hypothesis was first deduced by Salmena et al. ^[31] in 2011, many studies have confirmed this hypothesis. For the first time, the ceRNA hypothesis connects the functions of long noncoding RNA (lncRNA), messenger RNA (mRNA), and microRNA (miRNA), conceiving the regulatory network between them.

The increase in omics data allowed computational biologists to develop methods for mining RNA-RNA crosstalk and RNA-disease relationships, especially those involving ceRNAs. Cancerin ^[32], a LASSO regression-based method; Cernia ^[33], a support vector machine-based method; and SPONGE ^[34], a method based on multi-miRNA sensitivity correlation analysis have been developed to analyze RNA–RNA interactions based on statistics. Zhang et al. developed models for the prediction of lncRNA–miRNA association based on network distance analysis ^[35] and a semi-supervised interactome network ^[36]. Liu et al. proposed a deep forest-based approach (DFELMDA) to predict miRNA–disease associations ^[37]. Other machine learning-based models, such as graph neural networks (GNN), have also been used for predicting the association between biomarkers and diseases ^[38,39].

Consequently, five studies have constructed the lncRNA–miRNA–mRNA regulatory network for psoriasis based on the ceRNA theory ^{[40,41,42,43,44]}. In addition, Deng et al. analyzed gene regulatory networks in psoriasis using a tree-based machine learning approach ^[45]. These studies have served as a foundation for our study.

However, we discovered that most of these studies focused on the genetic level. Moreover, how the cell cycle of keratinocytes is regulated in psoriasis, how the CP1-3 regulatory points of the cell cycle function, and how the lncRNA-miRNA-mRNA regulatory network contributes to the overall pathogenesis of psoriasis, especially its relationship with the cell cycle, immune cells, immune functions, and immune pathways, remain unclear. In addition, potential therapeutic drug screening and molecular docking of small-molecule drugs with protein targets have also not been explored. Therefore, this study started from the above aspects to further develop previous research results.

First, we constructed a lncRNA–miRNA–mRNA regulatory model of psoriasis based on the ceRNA theory and microarray datasets different from those in previous studies. Then, we identified a hub gene in the ceRNA network and verified the diagnostic value of the hub gene for psoriasis by validating datasets and constructing a ceRNA subnetwork. Next, we followed the "cell cycle-antimicrobial peptide and chemokine-immune cells-cytokine-cell cycle" trail and the CP1-3 cell cycle model to analyze CP1-3 cell cycle regulatory point changes, including those of immune cells, immune functions, and immune pathways in psoriatic lesions, respectively. We also analyzed the relationship between the ceRNA subnetwork and CP1-3 cell cycle regulatory points, immune cells, and immune functions and pathways, deepening knowledge of the cell cycle of keratinocytes in psoriasis, including the interaction between keratinocytes and the immune system, to reveal the possible role of the ceRNA regulatory mechanism in the above processes. Finally, using the Connectivity Map (CMap), we screened potential therapeutic drugs for psoriasis, after which we used molecular docking to visualize the interaction patterns between a small-molecule drug and its protein targets.

2.   Methods and materials

2.1.   Data acquisition

Datasets that met the following criteria were included in our study: 1) studies that focused on patients with psoriasis, 2) studies whose data and platform information were complete and available for our research, and 3) studies that included psoriatic lesions and normal skin tissue samples. Nevertheless, those that included lesions and uninvolved skin tissue samples from patients with psoriasis were also selected. As a result, five datasets (GSE181318, GSE145054, GSE14905, GSE106992 and GSE117239) were downloaded from the GEO database. Specifically, GSE181318 contains lncRNA, circRNA, and the mRNA expression data from three psoriatic lesions and three normal skin tissue samples; GSE145054 contains miRNA expression data from four psoriatic lesions and four normal skin tissue samples; and GSE14905 contains mRNA expression data from 33 psoriatic lesions and 21 normal skin tissue samples. Subsequently, we used the three above datasets to construct a ceRNA network. While GSE106992 contains mRNA expression data of 67 lesions and 67 uninvolved skin tissue samples from untreated psoriasis patients, GSE117239 contains similar data from 83 lesions and 84 uninvolved skin tissue samples from untreated psoriasis patients, which we used to validate the accuracy of our results. Additionally, since all data were downloaded from a publicly available database (Gene Expression Omnibus), ethics committee approval was unnecessary.

2.2.   Data pre-processing and identification of differentially expressed lncRNAs, miRNAs and mRNAs

First, array probes in the three datasets were converted to matching gene symbols according to platform annotation information. Then, while the Perl language was used for gene biotype re-annotation of the lncRNA and mRNA expression data in GSE181318, the R software's "limma" package was used to normalize the raw data, after which differentially expressed lncRNAs (DElncRNAs) were identified, following the criteria of the adjusted P-value < 0.05 and |log2 (fold change)| > 1. Subsequently, the raw mRNA expression data in GSE181318 and GSE14905 were merged, and batch effects were removed using the R software "SVA" package. Next, we used the R software's "limma" package to normalize raw miRNA expression data in GSE145054 and the merged mRNA expression data. After that, we identified differentially expressed miRNAs (DEmiRNAs) and mRNAs (DEmRNAs), using the criteria of the adjusted P-value < 0.05 and |log2 (fold change)| > 1. Finally, the raw data in GSE106992 and GSE117239 were normalized using the R software's "limma" package for subsequent validation.

2.3.   Construction of the ceRNA regulatory network

With the file "highly conserved microRNA families" downloaded from the miRcode database, we predicted miRNAs binding to DElncRNAs. These miRNAs were then intersected with DEmiRNAs. Subsequently, we constructed a ceRNA network using the miRNAs at the intersection. First, the target mRNAs of miRNAs at the intersection were predicted using the miRDB and TargetScan databases. Only mRNAs with the same result in both databases could be exported. Then, the exported mRNAs were intersected with DEmRNAs. The mRNAs in the intersection were also used to construct a ceRNA network, after which we followed two principles to construct a ceRNA network: 1. each miRNA in the network must have a regulatory relationship with one or more lncRNAs and mRNAs simultaneously and 2. in the regulatory axis, the upregulation or downregulation trend of lncRNA should be similar to that of mRNA but opposite to that of miRNA. Finally, we used Cytoscape (v.3.9.0) to construct the lncRNA-miRNA-mRNA regulatory network of psoriasis.

2.4.   Pathway analysis

We performed gene enrichment using gene set enrichment analysis (GSEA, v.4.1.0). First, the merged and normalized mRNA expression data and groups were converted into "gct" and "cls" files, respectively, for submission to the GSEA software. Accordig to the criteria of the nominal P-value < 0.05 and a false discovery rate (FDR) < 0.025, the signaling pathways that were significantly enriched with differentially expressed genes between psoriasis and normal groups were uncovered. Then, GO annotation and KEGG pathway enrichment analysis of mRNAs in the ceRNA network were conducted using the R software "clusterProfiler" package. Finally, signaling pathways that were significantly enriched were identified, using P-value < 0.05 as the criterion.

2.5.   Construction of a protein–protein interaction (PPI) network and ceRNA subnetwork

The PPI network of mRNAs in the ceRNA network was constructed using the STRING database, setting the confidence score to > 0.500 and excluding noninteracting genes to simplify the network. Then, the simplified network was imported into Cytoscape (v.3.9.0). Subsequently, referring to the methods of previous studies ^[46,47], we calculated the degree, betweenness, and bridging values of each node in the network using CentiScaPe v.2.2. A top-ranking gene in terms of degree, betweenness, and bridging values was then identified simultaneously as the hub gene, after which the expression of the hub gene was verified in the GSE106992 and GSE117239 datasets. Finally, we obtained its upstream miRNA and lncRNAs from the previous ceRNA network and then constructed a ceRNA subnetwork based on this hub gene.

2.6.   Analysis of lncRNAs and miRNA in the ceRNA subnetwork

We used the ImmLnc database ^[48] to analyze the relationship between lncRNAs in the ceRNA subnetwork and various immune pathways and cells. Significant signaling pathways associated with the miRNA of the ceRNA subnetwork were first obtained using DIANA-miRPath v.3.0, based on TarBase (P < 0.05). Then, we used the ncRI database (http://www.jianglab.cn/ncRI/) ^[49] to finally retrieve the miRNA expression from the ceRNA subnetwork during inflammation.

2.7.   Immune infiltration algorithm

We derived marker genes for five CP1, four CP2, and three CP3 cell cycle regulatory point-associated pathways. Among them, the marker genes of the DDR pathway were obtained from a previous paper ^[50], and the marker genes of the remaining 11 pathways were derived from 11 clusters. We also derived marker genes for 18 immune functions and pathways from 18 clusters, after which we defined the above marker genes as metagenes using the GSVA algorithm and then encompassed the marker genes of 28 types of immune cells ^[51]. Subsequently, we used the R software "GSVA" package to calculate normalized ssGSEA scores for each sample in the previously merged and normalized mRNA expression data to evaluate the enrichment levels of cell cycle-related pathways, immune cells, immune functions, and pathways in each sample. We also compared the normalized ssGSEA scores of the psoriasis group with those of the normal group. Finally, Spearman's rank correlation coefficients and P-values were calculated for the mRNA in the ceRNA subnetwork with each cell cycle-related and immune pathway, including those for immune cell and function. Details of these metagenes are listed in Supplementary File 1. We used TIMER ^[52] to investigate the relationship between the expression level of the hub gene and immune infiltration in pan-cancer. Because the pathogenesis of melanoma is similar to psoriasis ^[53], we focused on the relationship between the expression level of the hub gene in melanoma with the abundance of six types of infiltrating immune cells: CD8+ T cells, CD4+ T cells, B cells, dendritic cells, macrophages, and neutrophils.

2.8.   Acquisition and analysis of scRNA-seq data

A scRNA-seq dataset (GSE176509) containing skin biopsy details of two patients with psoriasis was downloaded from the GEO database. First, we used the R language "Seurat" package to normalize the scRNA-seq data and perform quality control, excluding cells with a detected gene number of ≤ 50 and a proportion of mitochondrial genes per cell of ≥ 5%. Then, we normalized the gene expression of the cells using the LogNormalize function and identified the 1500 most highly variable genes. Subsequently, we used principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) methods to cluster cell samples based on the 1500 most variable genes. We used the CellCycleScoring function of the "Seurat" package to infer the S phase and G2M phase scores of every single cell and subsequently classified the cell cycle state of each cell into G1, S or G2M phases. We also used the R language "SingleR" and "celldex" packages for cell type annotation and the R "monocle" package for pseudotime analysis.

2.9.   Weighted gene correlation network analysis (WGCNA) and single-gene batch correlation analysis-based GSEA

WGCNA was constructed using the R language "WGCNA" package. First, the soft thresholding power (β) was set to 13, and the reliability of genes belonging to modules was described using module membership (MM). Next, the correlation between modules and clinical traits was calculated. We also calculated gene significance (GS) to assess the correlation between each gene and clinical traits. As a result, the hub genes of each module were obtained using the criteria of GS > 0.5 and MM > 0.8, including the correlation of each module with the standardized ssGSEA score of each type of immune cell. Finally, the expression levels of DEmRNAs were analyzed for batch correlation and ranked using Spearman's rank correlation coefficients with the expression level of the hub gene, followed by GSEA using the R language "clusterProfiler" package to assess the functions of the hub gene. The code used during single-gene batch correlation analysis-based GSEA is shown in Supplementary File 2.

2.10.   Screening of potential therapeutic drugs and molecular docking

CMap was screened for potential therapeutic agents. With the Perl language, gene symbols of up-regulated and downregulated mRNAs in the ceRNA network were translated into matching microarray probes and submitted to CMap. Drugs with negative mean values, enrichment correlation coefficients, and P < 0.05 were identified as promising therapeutic drugs. Subsequently, we retrieved the SMILES of the small-molecule drug from the PubChem database and used the SwissTargetPrediction to predict the target genes of the drug. First, the target genes were intersected with the ceRNA network's mRNAs. Then, the STITCH database was used to explore the interactions between the proteins encoded by the intersection/hub genes and the drug. Next, we searched the UniProt database for receptor proteins encoded by the intersection genes and downloaded the 3D structures of these proteins from the RCSB PDB database. Simultaneously, the 2D structure of the small molecular ligand was downloaded from the PubChem database, after which its energy-minimized 3D structure was calculated and derived using ChemBio3D (v.14.0). We also used PyMOL (v.2.5.2) to remove water molecules and separate ligands from these protein structures. However, AutoDockTools (v.1.5.7) was used to hydrogenate protein structures, after which we saved them as "PDBQT" files and constructed docking grid boxes for the protein structures. The 3D structure of the small molecular ligand was also saved as a "PDBQT" file using AutoDockTools (v.1.5.7). Finally, while AutoDock Vina (v.1.1.2) was used to dock the receptor proteins with the small-molecule ligand, PyMOL (v.2.5.2) was used to visualize the conformations with the lowest binding energy.

3.   Results

3.1.   CCND1 was significantly downregulated in psoriatic lesions

The workflow of this study is shown in Figures 1 and 2. Investigations identified 103 lncRNAs (42 were up-regulated and 61 were downregulated), 50 miRNAs (18 were up-regulated and 32 were downregulated), and 1159 mRNAs (616 were up-regulated and 543 were downregulated) to be differentially expressed between psoriatic lesions and normal skin tissue samples using an adjusted P-value < 0.05 and |log2 (fold change)| > 1 as criteria (Figure 3; Supplementary File 3). Notably, CCND1, encoding cyclin D1, was downregulated in psoriatic lesions (log2 (fold change) = −1.63, adjusted P-value = 9.80×10⁻¹⁵), which contradicts the results of some previous studies ^[54,55,56] and is consistent with only those of a few studies ^[57,58]. Furthermore, the classic cell cycle model suggested that the G1/S transition phase was inhibited. However, the CP1-3 cell cycle model suggested that cells in psoriatic lesions (predominantly keratinocytes) were blocked from passing through CP1 and CP2 ^[18].

3.2.   CeRNA network construction

From the miRcode database, 207 miRNAs binding to DElncRNAs were identified. The 207 miRNAs were intersected with 50 DEmiRNAs, and there were five miRNAs in the intersection. The overlap p-value is 0.01822255 (Figure 4A). We also identified 2311 target mRNAs at five intersecting miRNAs using the TargetScan and miRDB databases. Furthermore, while 2311 target mRNAs intersected with 1159 DEmRNAs, 150 mRNAs were also identified. The overlap p-value is 0.007226333 (Figure 4B). Then, based on the ceRNA hypothesis, Cytoscape (v.3.9.0) was used to construct a ceRNA network. The ceRNA network contained ten lncRNAs, four miRNAs, and 89 mRNAs (Figure 5). Specifically, we also observed that the miRNA that bound CCND1 : hsa-miR-17-5p was associated with the G1 phase cell cycle arrest ^[59], which suggests that cells in psoriatic lesions were prevented from entering the S phase.

3.3.   DNA damage as a significant factor in cell cycle dysregulation of psoriatic lesions

GSEA revealed that 14 pathways were significantly enriched (nominal P < 0.05 and false discovery rate < 0.025) in the psoriasis group (Figure 6A), and among these pathways, the pathway with the highest normalized enrichment score (NES) was base excision repair (BER, NES = 1.80). With a nominal P < 0.05 and false discovery rate < 0.025 as the criteria, no pathway was significantly enriched in the normal group. However, we observed that the normal group was negatively associated with several pathways, such as maturity-onset diabetes of the young (NES = −1.64, P = 0.020) and melanoma (NES = −1.46, P = 0.040). BER is a prominent mechanism of cellular DNA repair that can directly remove damaged bases from DNA, displace them, and regulate many pathological states associated with DNA damage ^[60]. It not only is an isolated DNA damage repair mechanism but also constitutes a regulatory network together with other DNA damage repair pathways ^[61]. This finding suggests that increased DNA damage could occur in psoriatic lesions. Based on these facts, GO annotation and KEGG pathway enrichment analysis were performed to further investigate the mechanisms and pathways associated with the ceRNA network. GO annotations revealed some functions that were related to psoriasis, for instance, the bone morphogenetic protein (BMP) signaling pathway (P = 6.46×10⁻⁵) ^[62] and beta-catenin binding (P = 0.007) ^[63] (Figure 6B, Supplementary File 4). Notably, four of the five mRNAs enriched in the guanyl-nucleotide exchange factor activity were significantly downregulated (DNMBP, SH3BP5, ARHGEF10 and ARHGEF28). Alternatively, the KEGG enrichment analysis indicated that mRNAs in the ceRNA network were significantly enriched in several pathways essential to psoriasis, for instance, the Hippo signaling pathway (P = 0.009) ^[64] and the TGF-β signaling pathway (TGF-β signaling pathway) (P = 0.012) ^[65] (Figure 6C). Among these pathways, the TGF-β signaling pathway is closely related to the keratinocyte's cell cycle ^[7]. Additionally, BMP has been reported to be part of the TGF-β superfamily ^[66]. Hence, investigations revealed that mRNAs enriched in both the TGF-β and BMP signaling pathways were also significantly downregulated (BMP2, SMAD7, RGMB, EGR1 TMEM100 and SORL1), suggesting the downregulation of the TGF-β pathway in psoriasis. This finding indicates that the G1 cell cycle arrest suggested by the significant downregulation of CCND1 and the increased production of β-defensins in psoriasis may be caused by other mechanisms different from the TGF-β pathway, for which we performed a ssGSEA. Complete GSEA, GO annotation, and KEGG pathway enrichment analysis results are presented in Supplementary File 4.

3.4.   CTGF is the PPI network hub and has a high diagnostic value for psoriasis

Subsequently, we explored the mRNA interaction relationships in the ceRNA network at the protein level using STRING (confidence score > 0.500) and removed disconnected nodes. Next, we constructed a PPI network containing 35 nodes (Figure 7A). Agreeing with previous studies ^[46,47], CentiScaPe v.2.2 was used to calculate the degree, betweenness, and bridging values for each node in the PPI network. Table 1 lists the top five nodes in the PPI network for these parameters. Investigations revealed that only CTGF (degree = 5, betweenness = 286, bridging = 49) ranked in the top five for all three values. Therefore, two datasets (GSE106992 and GSE117239) were selected to verify the specific expression of CTGF, which revealed that in both datasets, the expression level of CTGF in psoriatic lesions was significantly lower than that in uninvolved skin tissue samples (P = 3.2×10⁻⁹, P = 6.2×10⁻¹³), which was similar to our previous results (Figure 7B). The receiver operating characteristic (ROC) curves also demonstrated that CTGF had a high diagnostic value in psoriasis samples (Figure 7C–E). Based on these facts, we identified CTGF as the hub gene. We also identified upstream miRNA (hsa-miR-212-3p) and lncRNAs (LINC00173, MIR600HG and MAGI1-IT1) of CTGF from the ceRNA network and then constructed a ceRNA subnetwork based on these data (Figure 7F).

3.5.   Relationships between lncRNAs and miRNA in the ceRNA subnetwork and immune features

Using the ImmLnc database (lncRES > 0.995 and adjusted P-value < 0.05), we obtained immune pathways significantly associated with three lncRNAs in the ceRNA subnetwork (P-value < 0.05 and |correlation coefficient (R)| > 0.3). Complete results are presented in Supplementary File 5.

Tables 2 and 3 list the immune pathways most significantly associated with each lncRNA (the smallest adjusted P-values are also given). We identified immune cells negatively associated with each lncRNA (|R| > 0.4). The diseases in which these relationships existed were also detected. Our investigations indicated that the three lncRNAs were significantly associated with antimicrobials, cytokines, cytokine receptors, and chemokines, and were negatively correlated with dendritic, macrophage, and CD4(+) T cells. Subsequently, we used DIANA-miRPath v.3. 0, with P < 0.05 as the cutoff point, further revealing that hsa-miR-212-3p was significantly associated with eight signaling pathways, including the cell cycle, TGF-β signaling pathway, and Hippo signaling pathway. Therefore, based on these results, we searched the studies included in the ncRI database and discovered that the expression of hsa-miR-212-3p was up-regulated during inflammation ^[67,68,69], which agreed with our previous results (Table 4).

3.6.   Cells in psoriatic lesions were inhibited from passing CP1 and CP2 but significantly passed CP3 compared with normal skin tissue samples

We further investigated the interaction between the cell cycle, immune cells, and immune dysfunction in psoriasis, including the role played by the ceRNA subnetwork. For this analysis, the R language "GSVA" package was used first to calculate normalized ssGSEA scores for five CP1-related pathways, four CP2-related pathways, three CP3-related pathways, 28 immune cells, and 18 immune functions and pathways in each sample. Our results showed that among the five CP1-related pathways, while Rb was significantly up-regulated in the psoriasis group, mitogen-activated protein kinase (MAPK) signaling pathway, epidermal growth factor (EGF), and fibroblast growth factor (FGF) were significantly downregulated in the psoriasis group. Simultaneously, MYC was not statistically different between psoriasis and normal groups (Figure 8A). This finding suggests a mitogen decrease in psoriatic lesions. Combined with GSEA, an active BER in the psoriasis group was also suggested, indicating that increased levels of DNA damage and reduced mitogens in the psoriasis group inhibited CDK4/6 and triggered Rb dephosphorylation. This cascade is proposed to be accompanied by CCND1 downregulation, where the MAPK pathway might have played a significant role. Furthermore, we observed that among the four CP2-related pathways, TP53 regulates transcription of genes involved in G1 cell cycle arrest, p53 signaling pathway, regulation of DNA damage response by p53 class mediator resulting in transcription of p21 class mediator and DDR pathway were significantly up-regulated in the psoriasis group in comparison with the normal group. (Figure 8B). This finding served as further evidence to prove that DNA damage in psoriasis prevailed in competition with mitogens and activated p21 via p53, inhibiting CDK4/6 and dephosphorylating Rb. Indeed, the above mechanisms also play an essential role in determining whether cells pass CP1 ^[18].

In contrast to those in CP1 and CP2, our results indicated a significant increase in cellular passage through CP3 in psoriatic lesions compared with normal skin tissue samples. Investigations also revealed a key G1/S transcription factor, the E2F transcription factor family (E2Fs), which increased in the psoriasis group more than in the normal group. Furthermore, while an aggregation in E2F4 and E2F5 was observed in the G0 phase, E2F1, E2F2 and E2F3 were unexpressed. Still, in our results, E2F1 (adjusted P = 0.03), E2F2 (adjusted P = 1.75E−07), and E2F3 (adjusted P = 0.004) were also up-regulated in psoriatic lesions, which is another evidence that the cells are in a proliferative state. Specifically, while both the G1/S (P = 0.039) and G2/M phase transitions (P = 1.59×10⁻¹²) were more active in the psoriasis group than in the normal group, the G1/S phase transition in the psoriasis group was closer to that in the normal group relative to the G2/M phase transition (Figure 8C). Besides, while the normalized ssGSEA scores of 16 immune functions and pathways were significantly different between psoriatic lesions and normal skin samples (P < 0.05), 11 of them were significantly up-regulated in psoriatic lesion samples, such as beta-defensins (β-defensins, P = 1.58×10⁻¹⁰) and IL12 (P = 7.24×10⁻¹⁰). Some immune functions and pathways, such as PI3K/AKT signaling pathway (P = 1.38×10⁻⁸) and protein kinase C (PKC, P = 1.56×10⁻⁸), were significantly downregulated in psoriatic lesions compared with normal skin (Figure 9A). In addition, we discovered that the scores of 20 types of immune cells were significantly different between psoriatic lesion samples and normal skin samples. Of these, 17 were significantly increased in psoriatic lesion samples. Also, the three types of immune cells with the most significant differences were neutrophils (P = 1.51×10⁻¹⁴), activated CD4 T cells (P = 4.07×10⁻¹³), and activated CD8 T cells (P = 1.85×10⁻¹⁰), all of which were up-regulated in the psoriatic lesion samples, whereas the plasmacytoid dendritic cell and type 1 T-helper cell (Th1) were insignificantly different (Figure 9B). The heatmap shows the correlation between immune cells (Figure S1A). The principal component analysis results also showed that a normalized ssGSEA score of 28 immune cells could distinguish the psoriasis group from the normal group (Figure S1B).

3.7.   Correlation of CTGF with cell cycle and immune signatures

Among the CP1-3-related pathways, the results showed that although CTGF was strongly negatively correlated with some pathways, such as the DDR pathway (ρ = −0.75, P-value close to 0), TP53 that regulates the transcription of genes and is involved in G1 cell cycle arrest (ρ = −0.72, P-value close to 0), and G2/M phase transition (ρ = −0.68, P = 6.38×10⁻⁹), it was positively correlated with some other pathways such as the MAPK signaling pathway and fibroblast growth factor (Figure 10A, B, C). Furthermore, among the 28 subsets of immune cells, investigations revealed that while CTGF was strongly and inversely correlated with some immune cells, such as neutrophils, activated CD4 T cells, and CD56dim natural killer cells, it was positively correlated with some immune cells, such as mast cells and central memory CD4 T cells (Figure 10D). Moreover, of the 18 immune functions and pathways, while CTGF was negatively correlated with the p38 MAPK signaling pathway, IFN alpha (IFN-α) and β-defensins, it was positively correlated with the PI3K/AKT signaling pathway and PKC (Figure 10E). Complete results of the correlation analysis are presented in Supplementary File 6. Using TIMER to analyze the relationship between CTGF expression and immune infiltration in melanoma, we found a significant positive correlation between CTGF expression and the neutrophil and macrophage abundance in both primary and metastatic melanoma, but not with CD8+ T cells and B cells. Notably, CTGF showed a significant positive correlation with the abundance of CD4+ T and myeloid dendritic cells in metastatic melanoma but not in primary melanoma (Figure S2). The complete results obtained using TIMER are shown in Supplementary File 7.

3.8.   Cell cycle and cell trajectory analyses based on scRNA-seq data

GSE176509 contains 880 cells from skin biopsies of two patients with psoriasis. After quality control, 747 cells were used for subsequent analysis. Subsequently, PCA was used to first classify the cell samples (Figure S3). The distribution of P-values in each PC is shown in Figure S3. Using the t-SNE method, the 747 cells were divided into 6 clusters (Figure 11A). Cluster 4 was annotated as monocytes, and the remaining 5 clusters were annotated as dendritic cells (DCs) (Figure 11B). Subsequently, we performed cell cycle scoring for each cell. Figure 12A shows the cell counts in the G1, G2M, and S phases. Cell types are annotated in Figure 12B. A pseudotime analysis simulated the state of cell differentiation based on gene expression. All cells were sorted along trajectories to construct a pseudotime axis (Figure 13). Figure 13B and Figure 13C show cell clusters and cell types, respectively.

3.9.   Further exploration of CTGF functions

After calculations, no outlier samples were removed, and all 60 samples were used to construct the adjacency matrix (Figure S4A). The different modules' dynamic tree was constructed (Figure S4B), after which the soft thresholding power "β" was set to 13 (Figure S4C). WGCNA identified seven color-coded co-expression modules, with CTGF being one of the hub genes in the brown module (GS = 0.51, P = 3.31×10⁻⁵; MM = 0.83, P = 2.65×10⁻¹⁶). The brown module was significantly correlated with clinical traits. It was also significantly correlated with activated CD4 T cells, plasmacytoid dendritic cells, neutrophils, and mast cells, which play essential roles in psoriasis (Figure 14A). Single-gene batch correlation analysis-based GSEA investigations also revealed that while CTGF might inhibit allograft rejection, IFN-α responses, and IFN-γ responses, suggesting that CTGF plays an anti-inflammatory role in psoriasis. But the results also showed that it might be associated with the activation of TNF-α. Moreover, it could inhibit the mitotic spindle, MYC targets, E2F targets, G2/M checkpoint, and DNA repair, suggesting that CTGF plays a suppressive cell proliferation role in psoriasis (Figure 14B).

3.10.   Screening of potential therapeutic agents and prediction of drug targets

Subsequently, we examined the pragmatic implications of this study using 89 mRNAs in the psoriasis ceRNA network to predict compounds that might affect psoriasis and CMap. First, potential therapeutic drugs were screened based on their mean values and enrichment correlation coefficients. Results showed MS-275 (also known as entinostat) as the most promising potential treatment option. Table 5 lists the observed predictions obtained from CMap with negative mean and enrichment correlation coefficients, including those with non-null values of 100% and p < 0.05. Next, we further explored the role of MS-275 in psoriasis and predicted 109 target MS-275 genes using SwissTargetPrediction.

The results showed that while the target genes of MS-275 intersected with 89 mRNAs in the ceRNA network, three intersected genes were obtained. Notably, the probability value of F3 among them was 1, indicating a 100% probability that F3 is the target gene of MS-275, according to SwissTargetPrediction (Table 6). We also calculated the Spearman rank correlation coefficients between the genes, and the correlation coefficients were 0.68 for CTGF and F3, 0.68 for CTGF and PIK3R1, and 0.42 for CTGF and FKBP5 (Figure S5A). Subsequently, an interaction network between MS-275, CTGF, F3, PIK3R1, and FKBP5 was obtained by entering them into the STITCH database and setting a confidence score greater than 0.500 (Figure S5B). In addition, we observed that genes closely related to psoriasis, such as STAT3, FOXO1, and MTOR, appeared in the network, proving an interaction between the drug, the intersecting genes, and the hub gene in psoriasis. Therefore, we focused on the interaction between the proteins encoded by the three intersecting genes and MS-275. The complete prediction results obtained using CMap and SwissTargetPrediction are presented in Supplementary File 8.

3.10.   Compound-target docking

Molecular docking was used to verify whether MS-275 played an essential role in regulating the proteins encoded by the three intersecting genes. This result demonstrated that MS-275 was closely associated with the proteins encoded by the three intersecting genes. For example, the affinity between MS-275 and tissue factor (encoded by F3) was −8.8 kcal/mol, forming a stabilized complex with residues PRO-40 and THR-112 by two hydrogen bonds (Figure 15A). Investigations also revealed that the affinity between MS-275 and phosphatidylinositol 3-kinase regulatory subunit alpha (encoded by PIK3R1) was −5.9 kcal/mol, forming a complex stabilized with residues TYR-142 and ASP-110 by two hydrogen bonds as well (Figure 15B). However, the affinity between MS-275 and peptidyl-prolyl cis-trans isomerase FKBP5 (encoded by FKBP5) was −7.2 kcal/mol, forming a complex stabilized by five hydrogen bonds with residues TYR-327, VAL-288, SER-296 and GLU-164 (Figure 15C).

4.   Discussion

A complex regulatory relationship exists between DNA damage and the cell cycle. Although Borska et al. found increased levels of DNA damage in patients with psoriasis compared to healthy individuals ^[70], the relationship between DNA damage and the cell cycle in psoriasis has not been fully elucidated. UV can cause DNA damage in keratinocytes during the treatment of psoriasis, which can effectively alleviate psoriasis ^[71]. However, in contradiction to this, some researchers have also reported that UV irradiation can lead to the production of reactive oxygen species (ROS) in keratinocytes, and the DNA damage caused by ROS can exacerbate inflammation and downregulate the TGF-β pathway, thus aggravating the condition of patients with psoriasis ^[72]. Exploring the regulatory relationship between DNA damage and the cell cycle in psoriasis is meaningful for understanding the pathogenesis of psoriasis and thus guiding treatment.

DNA damage exacerbates the condition of psoriasis patients mainly by promoting inflammation. At the same time, a growing number of researchers have suggested that the interaction between keratinocytes and the immune system plays an important role in the pathogenesis of psoriasis ^[7]. However, the role of ceRNA regulatory mechanisms in the interaction of keratinocytes with the immune system has also not been elucidated. Our study attempts to shed light on this enigma.

We constructed a ceRNA network with ten lncRNAs (three up-regulated, seven downregulated), four miRNAs (two up-regulated, two downregulated), and 89 mRNAs (13 up-regulated, 76 downregulated). Based on our extensive literature search, the expression trends of 10 lncRNAs in the ceRNA network observed in this study were all reported for the first time in psoriasis. Among the four miRNAs, the upstream miRNA of CTGF (hsa-miR-212-3p) was significantly associated with the cell cycle, the Hippo signaling pathway, and the TGF-β signaling pathway. Specifically, the Hippo signaling pathway also appeared in the results of the KEGG enrichment analysis. This pathway regulates cell proliferation and apoptosis, including heart development and osteoblast differentiation ^[64]. However, whether its dysregulation is associated with psoriasis comorbid with cardiovascular disease and complicated by psoriatic arthritis remains unexplored. Furthermore, studies have reported that hsa-miR-17-5p is closely associated with the cell cycle and the TGF-β signaling pathway ^[59,73], further suggesting the essential roles of cell cycle dysregulation and TGF-β signaling pathway dysregulation in psoriasis pathogenesis. Hsa-miR-193a-3p is downregulated in cutaneous melanoma. Moreover, the overexpression of miR-193a inhibits the expression of the MAPK signaling pathway and PI3K signaling pathway-related genes, such as ERBB2, KRAS, PIK3R3, MTOR and NUSAP1 ^[74], suggesting that it may also be an effective psoriasis treatment. In this study, among the 89 mRNAs, the expression trends of some mRNAs were also reported for the first time in psoriasis, such as F3, FKBP5, PIK3R1, KDM5A, and MYO5C.

Figure 16 illustrates the possible interaction mechanisms between the ceRNA subnetwork, the cell cycle, and the immune system in psoriasis. DNA damage levels have been significantly elevated in patients with psoriasis compared with healthy individuals ^[70]. Several studies have also highlighted that proliferating cells can undergo cell cycle arrest as a response to DNA damage, with a mechanism involved in this process being the downregulation of cyclin D1, usually associated with the activation of GSK3β, phosphorylation of cyclin D1, and degradation of cyclin D1 by proteases ^[75]. A previous genome-wide CRISPR screening revealed that CCND1 was required for the G1/S phase transition in 79.8% of all 769 cell lines ^[18]. Specifically, CCND1 downregulation suggests that cells in psoriatic lesions (mostly keratinocytes) are blocked from entering the S phase, which according to the classical view explains the increased keratinocyte differentiation in psoriasis ^[4] but not the proliferation of keratinocytes. Additionally, the CP1-3 cell cycle model suggests that while cells in psoriatic lesions are inhibited from passing through CP1 and CP2 regulatory points and are significantly more likely to pass through CP3 than the normal group, overall, they still have an increased passage through the S phase compared to the normal group (P = 0.039), explaining keratinocyte proliferation in psoriasis. However, since CP1 and CP2 regulatory sites are located upstream of CP3, more studies should still elucidate this phenomenon. It has also been reported that the activation of the TGF-β signaling pathway can inhibit the DDR pathway, leading to decreased DNA repair capacity and increased mutations ^[76]. Therefore, we speculate that the upregulation of the DDR pathway in psoriasis may be related to the loss of the inhibitory effect of the TGF-β signaling pathway on it, suggesting that the interaction between the two pathways needs further study. Alternatively, DNA damage causes not only a DDR but also a DNA damage differentiation response (DDDR), leading to increased endoreplication ^[77]. This result is consistent with the increased keratinocyte volume observed in psoriasis and the severe DNA damage, decreased guanyl-nucleotide exchange factor activity, cell cycle dysregulation, increased E2Fs, and up-regulated CCNE1 and CCNE2 suggested by the results in this study. Still, this study did not observe the inactivation of cyclin A and cyclin B, including the blockade of cells from the G2 to M phases. Hence, endoreplication in psoriatic lesions should still be confirmed by further studies. Endoreplication may avoid, to some extent, the increase in DNA instability. However, it has been shown that tumor cells complete mitosis and cytokinesis, even with unstable DNA ^[78]. Thus, whether endoreplication is an essential biological process that distinguishes nontumors from tumors (e.g., psoriasis and epidermal melanoma, which have many similarities in their pathogenesis) remains uninvestigated.

TIMER performs deconvolution based on linear least square regression to assess the abundance of six types of immune cells, reducing the correlation between CD8+ T cells and CD4+ T cells by multiple deconvolutions of the CD8+ and CD4+ T cell ratios ^[52]. Various sets of marker genes for immune cells are necessary for ssGSEA. ssGSEA is performed according to ranking genes by their expression in the sample, scoring them according to their ranking weights, and obtaining a normalized enrichment score by calculating the difference between the cumulative density of genes present and not present in the gene set. An accurate assessment of immune cell abundance can be obtained using both methods. Nevertheless, it is sometimes necessary to use both methods simultaneously for increased accuracy. TIMER analysis showed that CTGF expression was significantly positively correlated with neutrophil abundance but not with CD8+ T cells in melanoma. However, ssGSEA showed that in psoriasis, CTGF expression was significantly negatively correlated with the abundance of both neutrophils and CD8+ T cells. Melanoma shares many similarities with psoriasis in terms of pathogenesis, such as elevated TNF-α, IL-12, and IL-17A concentrations and excessive proliferation of keratinocytes ^[79]. However, the relationship between psoriasis and melanoma has shown contradictory results in different studies ^[53]. It has also been shown that patients with psoriasis treated with biologics have an increased risk of developing melanoma ^[80]. However, the risk of melanoma in patients with psoriasis and the role of different treatments in carcinogenesis remain to be further investigated. And it also remains to be explored whether the increased risk of developing melanoma in psoriasis patients treated with biologic agents is associated with changes in gene expression profiles.

CTGF is downstream of the TGF-β pathway, and its downregulation promotes cell proliferation by some mechanisms, including increasing insulin-like growth factor 1 (IGF-1), activating the nonclassical Wnt pathway, and inhibiting vascular endothelial growth factor (VEGF) loss ^[81] (which might be related to increased angiogenesis in psoriasis). This study discovered that CTGF was associated with various cell cycle regulators, immune cells, immune functions, and pathways, which could play an anti-inflammatory and cell proliferation inhibitory role in psoriasis. Some recent studies have also shown that CTGF inhibits cell proliferation and promotes cell apoptosis ^[81]. However, in psoriasis, whether CTGF downregulation leads to cell proliferation or cell proliferation leads to feedback downregulation of CTGF requires further experimental confirmation. For example, it has been studied that the application of anti-CTGF monoclonal antibodies exacerbated skin lesions in mice with imiquimod-induced psoriasis ^[82]. Additionally, the upstream and downstream regulatory relationships of the ceRNA subnetwork show a sequence consistent with the interaction between keratinocytes and the immune system in psoriasis, suggesting that ceRNA regulatory mechanisms not only act at the genetic level but also play a complex role in the overall pathogenesis. Recently, some researchers have observed that under the influence of LL-37, noncoding RNAs can stimulate adhesion molecules in endothelial cells ^[83], emphasizing the vital role of noncoding RNAs. Still, more work should confirm our hypothesis (Figure 9).

We searched for the hub gene by calculating each node's degree, betweenness, and bridging values in the PPI network. In addition to CTGF and CCND1, we found some other cell cycle-related genes in Table 1. For example, EGR1 ranks in the top five for both degree and betweenness values. The MAPK cascade can activate EGR1, and overexpressed EGR1 can activate the MAPK signaling pathway, forming a positive feedback loop ^[84]. EGR1 can up-regulate the expression of some cell cycle-related proteins such as cyclin D1 and cyclin D2, to accelerate cell proliferation ^[85]. The downregulation of EGR1 may lead to p53 inactivation and dysfunctional DNA damage response, revealing that EGR1 may play an essential role in regulating the DNA damage response ^[86]. In addition, EGR1 is one of the activators of the TGF-β pathway ^[87]. BMP2 belongs to the TGF-β superfamily. It has been indicated that downregulation of BMP2 is associated with G1 cell cycle arrest ^[88]. In addition, the BMP signaling cascade promotes regulatory T cell production and alleviates inflammation ^[62]. SMAD7 is also a key regulator of the TGF-β pathway, and its downregulation has been found to be associated with decreased apoptosis ^[89]. GPC6 has been reported to promote cell proliferation in nasopharyngeal carcinoma ^[90] and may be a biomarker for melanoma metastasis ^[91]. Downregulation of SDC2 promotes endothelial cell migration and regulates angiogenesis ^[92]. NEDD4L inhibits IL-6 signaling in keratinocytes and promotes the degradation of GP130, thereby suppressing the progression of psoriasis ^[93]. Overexpression of TXNIP can inhibit cell proliferation by promoting the production of reactive oxygen species ^[94]. Meanwhile, 12-myristate 13-acetate, an inducer of keratinocyte differentiation, was reported to inhibit TXNIP expression ^[95]. The role of TXNIP in psoriasis remains to be explored. With increasing in-depth research on the cell cycle, increasing numbers of novel cell cycle-related genes have been identified, such as the proteasome 26S subunit, ATPase (PSMC) family genes ^[96], ABCB1 ^[97], and FMR1 ^[98]. In the future, we aim to develop methods to more accurately identify cell cycle-related information and explore the heterogeneity of gene regulation during the cell cycle in tumor and nontumor diseases.

Drug repurposing has the advantages of shorter development time, lower development costs, and lower risk of failure than traditional drug development strategies. Drug repurposing has been effective in diseases such as acquired immune deficiency syndrome (AIDS), breast cancer, and coronavirus disease 2019 (COVID-19) ^[99,100,101]. However, despite developments in treatments against psoriasis, there are still many limitations to currently approved drugs. Long-term use of high doses of glucocorticoids is associated with adverse effects such as peptic ulcers, osteoporosis, and withdrawal syndrome. In addition, 30%-50% of patients respond poorly to biologics such as etanercept, adalimumab, and secukinumab, especially those with generalized pustular psoriasis ^[102]. Considering the current lack of research in drug repurposing in psoriasis ^[103], our study that identified MS-275 as a potential therapeutic agent for psoriasis is remarkable for its significance and novelty.

MS-275, also known as entinostat, is a histone deacetylase inhibitor (HDAC inhibitor, HDACi) that selectively inhibits classes I and IV histone deacetylases. It has been demonstrated that inhibiting both classes I and IV HDACs results in enhanced regulatory T cell (Treg) function and greater stability of Foxp3 (forkhead box protein P3) ^[104], which could explain its potential therapeutic effect on psoriasis. Moreover, HDACi can act as epigenetic modifiers; regulate chromatin conformation ^[105], acetylation, and DNA methylation; and alter gene and miRNA expression ^[106]. In addition, the binding of hypophosphorylated Rb to E2Fs can recruit histone deacetylases associated with regulating multiple transcription factors ^[107]. Besides, although HDACi is mainly used to treat liquid tumors (e.g., leukemia and lymphoma), it has recently also demonstrated the ability to treat autoimmune diseases. For example, in rheumatoid arthritis, HDACi reduced TNF-α ^[108], and in an animal SLE model, it reduced pro-inflammatory cytokines ^[109]. HDACi can also reduce MIP-1α and MCP-1, but this appears dose-dependent, as in osteoblasts, where high doses of MS-275 are required to exert anti-inflammatory effects ^[110]. Furthermore, while HDACs can aberrantly modify both histones and nonhistone proteins that make up chromatin and further alter cell proliferation and differentiation, they are significantly higher in psoriatic lesions than in healthy skin. Therefore, they are likely associated with clinical subtypes of psoriasis ^[111]. Moreover, our study suggests that MS-275 is expected to bind spontaneously to key proteins in psoriasis. Still, more animal experimentation and clinical trials should verify the efficacy and adjust the dose of MS-275 and other HDAC inhibitors. Also, in-depth studies on HDAC inhibitors will help design next-generation drugs.

This study has the following limitations. In terms of internal validity, this study identified CTGF as the hub gene and analyzed the function of CTGF in psoriasis, concluding that CTGF may inhibit cell proliferation in psoriasis. However, because of the complexity of cell cycle regulation, whether and how CTGF can directly regulate cell proliferation needs to be verified by cell and animal experiments. In terms of external validity, this study identified MS-275 (entinostat) as a potential therapeutic agent for psoriasis. However, its effectiveness in treating psoriasis needs to be further evaluated in animal experiments and clinical trials. In terms of methodology, this study used multiple bioinformatics methods with online databases for analysis, but no new cell cycle and RNA‒RNA crosstalk analysis methods were developed. Deep learning and artificial intelligence require new methods to more accurately resolve cellular information from microscopy, image flow cytometry, time-course gene expression profiling data, and scRNA-seq data, as well as to mine RNA‒RNA interaction patterns from gene sequence and gene expression profiling data. This will be the focus of a future study.

Acknowledgments

This study was supported by the Foundation of The First People's Hospital of Yibin (2022-KYY-14).

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix