The Gene Expression Omnibus (GEO) is a database created by the NCBI that stores a large amount of gene expression data. We obtained the dataset GSE181540 from GEO (https://www.ncbi.nlm.nih.gov/geo/), which includes 3 HS and 3 NS tissues, as well as 6 NS tissues from GSE158395, 3 NS tissues from GSE190626, and 3 HS tissues from GSE210434. In addition, we obtained dataset OEP002674 from The National Omics Data Encyclopedia (NODE, https://www.biosino.org/), which includes 5 HS and 5 NS samples. Based on previous studies on cuproptosis-related genes (CRGs), we identified a set of 56 genes related to cuproptosis (10, 13–15).

We annotated, normalized, and performed log2 transformation on different datasets using R. Then, we used the “sva” R package to merge the aforementioned datasets, correct batch effects, and extract CRGs from the merged dataset. We performed differential gene analysis on the processed dataset using the “limma” package and screened for differentially expressed genes. The threshold was set to |log2 Fold change (FC)| > 1 and p-value <0.05.

We performed Gene Ontology (GO) enrichment analysis on the genes from three aspects: biological process, cellular component, and molecular function, in order to identify the biological functions of the genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a primary public database related to pathways. We used the hypergeometric test to identify significantly enriched pathways in the genes. We performed functional and pathway enrichment analysis on the differentially expressed genes using these two methods, and set the P-value threshold to 0.05.

We used the “randomForest” R package to perform Random Forest (RF) feature selection on the dataset. The RF algorithm is an ensemble learning method that uses multiple decision trees to form a regressor. It is used to screen for the importance of genes in diseases. Its main advantages are its good performance in processing high-dimensional data and insensitivity to outliers. Then, we used the “e1071” R package to perform support vector machine-recursive feature elimination (SVM-RFE) feature selection on the gene set. SVM-RFE scores and ranks gene features using binary classification and selects the top few genes with the lowest error rate. We took the intersection of the results from these two machine learning methods and the differentially expressed genes that we previously screened for, and used the 11 intersected genes as hub genes for further analysis.

We selected 5 HS samples and 3 NS samples from the dataset GSE188952 in GEO to perform external differential expression validation of the hub genes. We calculate the differences and significance of the genes between the disease group and the control group, and plotted the gene differential boxplot. We selected the genes that were validated and further analyzed them.

The three selected genes were subjected to multivariable logistic regression modeling. Using the generalized linear model with binary classification, the entire dataset was used as the training set to construct the model. ROC curve was plotted and the area under the curve (AUC) was calculated. Finally, bootstrap validation was performed with 1000 resampling iterations to calculate the sensitivity and specificity of the model. Further evaluation of the stability and reliability of the model was conducted through construction of nomogram and decision curves.

Previous studies have demonstrated a significant association between cuproptosis and cellular immunity (16–18). In this study, we performed a Pearson correlation analysis using the corrplot R package to investigate the relationship between the mRNA expression of the 3 CRGs. Additionally, we used the “IOBR” R package to apply the CIBERSORT and MCPcounter algorithms for quantifying the degree of immune cell infiltration in the samples (19). We also analyzed the correlation between the 3 CRGs and immune cell infiltration as well as the expression of inflammatory factors.

We utilized the “ConsensusClusterPlus” R package to perform a k-means consensus clustering analysis, aiming to identify distinct subtypes that are associated with the expression of CRGs.

The differentially expressed genes(|log2 FC|>1, p-value<0.05) identified between two subtypes was used for GO biological function enrichment analysis and KEGG pathway enrichment analysis by using R package “clusterprofiler” and “enrichplot”

Single-cell sequencing dataset analysis has been used to explore a variety of skin diseases with promising results (20–22). In our study, we used the “Seurat” R package to create a Seurat object and normalize and scale cells. By excluding cells with less than 200 expressed genes, more than 2500 expressed genes, or mitochondrial gene content higher than 10%, we retained 45095 cells. We selected the top 3000 highly variable genes for principal component analysis with a setting of 15 principal components for subsequent dimensionality reduction and clustering. We integrated samples and removed batch effects using the “harmony” function, obtained unsupervised cell clusters using a graph-based clustering method (with a resolution of 0.5), and visualized them using t-SNE plots. We used the “FindAllMarkers” function and Wilcoxon rank-sum test algorithm to identify marker genes for each cell cluster, with conditions including |log2 FC|>0.25, P<0.05, and a minimum percentage>0.1. To better identify fibroblast clusters, we selected clusters of fibroblast cell types for further t-SNE analysis, graph-based clustering, and marker gene analysis. We calculated the score of the cuproptosis in each cell using the “AddModuleScore” function and divided cells into high and low groups based on the median score.

We used the recently developed tool “CellChat” R package to perform comprehensive analysis of cell-cell communication molecules. This tool can generate and plot the probability and interaction strength of cell-cell communication from single-cell transcriptome data, allowing for in-depth analysis of cell-cell communication. The normalized count and cell types by Seurat were used for this analysis.

The SCENIC transcription factor inference was performed using the R SCENIC package following the recommended pipeline steps (23).

We collected 3 NS tissues and 3 HS tissues. Prior to surgery, we provided patients with a detailed explanation of the research purpose and process, and obtained their consent to participate in the study. We also obtained written informed consent from all patients or their legal guardians and received approval from the Ethics Committee of the Xijing Hospital, Fourth Military Medical University.

The tissues were subjected to RNA extraction using the TRIzol reagent (Invitrogen, Camarillo, CA, United States). The extracted RNA was then reverse transcribed into cDNA using the PrimeScript RT reagent Kit with gDNA Remover (Takara, Shiga, Japan) at a concentration of 1,000 ng. For qRT-PCR, the TB Green Premix Ex Taq II (Takara, Shiga, Japan) was utilized, and the BIO-RAD CFX Connect Real-Time System (Bio-Rad, Munich, Germany) was used for the analysis. The expression levels of target genes were normalized to GAPDH. Forward 5′-GCACCGTCAAGCTGAGAAC-3′ and reverse 5′-TGGTGAAGACGCCAGTGGA-3′ for human GAPDH, forward 5′-GGCAAGTTCGAGTTCTCCCG-3′ and reverse 5′-CGACCTCCAAATCGTGCTTCT-3′ for human ULK1, forward 5′-CAGTGCGGAGAACACTTGC-3′ and reverse 5′-TGCACATAACCCTGGGACATT-3′ for human MTF1, forward 5′-TGGGACCATATCCAAAGCACA-3′ and reverse 5′-CAAGCTCCTGACTTTTTCTCCAG-3′ for human ATP7A.

The research process of this study can be seen from the graphic flow chart ( Figure 1 ). The “Limma” R package was used to identify differentially expressed CRGs in the dataset. The filtering criteria were set as |log2 FC| > 1 and p-value <0.05, and boxplots and volcano plots of the differentially expressed genes were generated ( Figures 2A, B ). Figure 2C shows the correlation between differentially expressed genes. Among the 56 CRGs, 15 genes were found to be significantly differentially expressed, including 14 downregulated genes and 1 upregulated gene.

The differentially expressed genes were subjected to enrichment analysis using the GO and KEGG methods, and enrichment circle plots were generated for both ( Figures 2D, E ). In GO analysis, the genes were mainly enriched in cellular copper ion homeostasis, copper ion transport, and copper ion homeostasis, while in KEGG analysis, the genes were mainly enriched in the HIF-1 signaling pathway, VEGF signaling pathway, and Mineral absorption pathways.

Machine learning methods are widely used in dermatological diseases (24–26). The SVM-RFE and RF machine learning methods were used to further screen candidate hub genes. First, the RF algorithm was used to rank the importance of each gene ( Figures 3A, B ), and the top 20 genes with the highest importance ranking were selected: ULK1, LOXL2, VEGFA, MAP2K1, ATP7A, PDE3B, LIAS, DBT, CDKN2A, PDHA1, COA6, COX11, SLC25A3, DLAT, SOD1, AOC3, MTF1, TYR, DBH, ATOX1. Then, the top 20 genes with the highest accuracy in the SVM-RFE results were selected: LOXL2, ULK1, DBT, PDE3B, LIAS, COX11, TYR, CCS, ATP7A, ATOX1, VEGFA, DLAT, SCO1, PDK1, MAP2K1, SLC25A3, PDHA1, DLST, DBH, MTF1 ( Figure 3C ). The intersection of the two machine learning methods’ results was taken as the candidate hub genes for further analysis.

The candidate hub genes selected using machine learning and the 15 significantly differentially expressed genes obtained from differential analysis were intersected to obtain 11 genes ( Figure 3D ): MTF1, ULK1, ATP7A, TYR, DBT, SLC25A3, LIAS, COX11, DLAT, VEGFA, PDE3B. The external dataset GSE188952 was used for validation of differential gene expression ( Figure 3E ), and 3 CRGs with significant differential expression in the external dataset were identified: MTF1, ULK1, ATP7A.

To further verify the expression levels of the three CRGs, we collected specimens from 3 patients with NS and 3 HS with patients. The mRNA expression level of MTF1, ULK1 and ATP7A in NS was significantly higher than that in HS, which also consistent with the expression of bulk RNA-seq datasets ( Figure 3F ).

Firstly, a correlation graph between 3 CRGs was drawn ( Figure 4A ), and ROC curve was plotted to evaluate the diagnostic specificity and sensitivity of each gene, and the area under the curve (AUC) was calculated. ATPA (AUC=0.79), ULK1 (AUC=0.92), and MTF1 (AUC=0.76) were obtained ( Figure 4B ). We performed multivariable logistic regression modeling on three hub genes, and the model was evaluated using ROC curve ( Figure 4C ). The AUC was calculated to be 0.914. Further verification was performed using the bootstrap algorithm with R=1000 resampling times, and the average AUC value was0.879893, with a 95% confidence interval (CI) of 0.7303-0.9189 ( Figure 4D ). These results indicate that the diagnostic model has high diagnostic value. A nomogram of the model was drawn, and the total score of the sample was calculated ( Figure 4E ). Finally, the Decision Curve Analysis (DCA) curve was plotted ( Figure 4F ), and it showed that the diagnostic model composed of MTF1, ULK1, and ATP7A had the highest diagnostic performance.

The “IOBR” package in R was used to conduct MCPcounter and CIBERSORT immune infiltration analysis on the 3 CRGs ( Figures 5A, B ). In MCPcounter, almost all immune cells showed a significant positive correlation with 3 CRGs, which were upregulated in the NS group. In CIBERSORT, dendritic cells resting showed a significant positive correlation with CRGs, T regulatory cells showed a positive correlation with ULK1 and ATP7A, and macrophages M0 showed a negative correlation with MTF1 and ATP7A. In addition, we generated correlation plots of the top 5 associations between CRGs and the corresponding immune cells using the CIBERSORT and MCPcounter algorithms ( Figure 5C ). In the analysis of inflammatory factors ( Figure 5D ), IL5, TNF, IL7, PDGFA, CD4, CSF1, HLA-DRB5, HLA-DRB1, HLA-DRA, and other factors showed significant positive correlations with 3 CRGs, while CSF3, IL6, and other factors showed significant negative correlations. Through the analysis of immune infiltration and inflammatory factors, it can be concluded that the CRGs have an extremely important influence on the immune and inflammatory reactions in the development of the disease.

We performed consensus clustering analysis to group the dataset of HS, setting the maximum number of clusters to 5, and scoring the 2-5 clusters ( Figures 6A, B ). The results showed that when the number of clusters was k=2, the inter-cluster correlation was the lowest and the intra-cluster correlation was the highest, and the score was the highest. Therefore, we chose k=2 to group the dataset of HS and analyzed the gene expression of the two groups, and generated heatmaps and differential boxplots ( Figures 6C, D ). Finally, we further analyzed the expression differences of inflammatory factors and immune cells in the two groups ( Figures 6E, F ).

To investigate functional and pathway disparities between subtypes and unravel potential mechanisms of disease progression, we conducted GO and KEGG enrichment analysis using the differentially expressed genes identified between the two subtypes. In the BP category, the genes were mainly enriched in functions such as epidermis development, skin development, axon development, and regulation of neuron projection development. In the CC category, the genes were mainly enriched in functions such as cell-substrate junction, cell leading edge, focal adhesion, and chromosomal region. In the MF category, the genes were mainly enriched in functions such as GTPase regulator activity, nucleoside-triphosphatase regulator activity, and cadherin binding ( Figures 7A, B ). In the KEGG analysis, the genes were mainly enriched in metabolic pathways such as Human papillomavirus infection, Rap1 signaling pathway, Ras signaling pathway, and Focal adhesion ( Figures 7C, D ).

We utilized scRNA-seq datasets containing 3 NS and 3 HS samples from the GEO database to unveil the intrinsic cellular heterogeneity within the skin tissue. After strict quality control, low-quality cells were excluded, and a total of 45095 cells were selected for subsequent analysis. Using dimensionality reduction clustering, we identified 16 cell clusters, which displayed high heterogeneity between different cell populations ( Figure 8A ). Based on marker genes identified in previous studies, we annotated 8 cell types, including Melanocytes (cluster 14, marked by MITF), Langerhans cells (clusters 8, 9, marked by Langerhans cells), T cells (clusters 7, 12, marked by PTPRC), Sweat gland cells (cluster 18, marked by SCGB1B2P, SCGB1D2), Fibroblasts (clusters 0, 2, 4, 5, 10, 15, marked by FBLN1, COL1A1, APOE, and APCDD1), Keratinocytes (clusters 11, 13, marked by KRT1 and KRT14), Smooth muscle cells (cluster 3, marked by ACTA2 and RGS5), and Endothelial cells (clusters 1, 6, marked by THBD and SELE) ( Figure 8B ). We then explored the TOP5 essential genes in seven previously annotated cell clusters and visualized the results ( Figure 8C ). Figure 8D also shows the proportion of fibroblasts in NS and HS. Further analysis of fibroblasts was performed by extracting them and reanalyzing them using tSNE and UMAP, which resulted in their re-clustering into 8 clusters. We scored the activity of cuproptosis in each cell using the AddModuleScore function, and found that the cuproptosis score of fibroblasts in NS was higher than that in HS ( Figure 8E ). Specifically, clusters 1 and 6 had the highest activity of cuproptosis ( Figure 8F ).

To decipher cell-cell signaling, we performed cell-cell communication analysis between different cell types using the “cellchat” R package. We classified fibroblasts into two types based on their previous cuproptosis activity scores, dividing them into high and low cuproptosis activity groups using the median value. Aggregated cell-cell communication networks were constructed based on interaction numbers ( Figure 9A ) and interaction weights ( Figure 9B ). The interaction strengths of cells transmitting and receiving signals were plotted in Figure 9C , which showed that fibroblasts played a crucial role as information senders in cell-cell communication. Low cuproptosis activity fibroblasts had weaker signal transmission strengths than high cuproptosis activity fibroblasts. We compared the interactions between the two types of fibroblasts and other cells. Compared to fibroblasts with low pathway activity, those with high pathway activity can engage in additional cell communication through the MIF pathway and with T cells and Langerhans cells. Additionally, they can communicate with smooth muscle cells via PROS1-AXL interactions.

Figures 10A, B show the top 5 transcription factors with the highest activity, including NFIA, HES1, CEBPB, NFE2L2, and STAT3 in NS and CREB3L1, FOSB, JUND, KLF6, and KLF4 in HS. Figure 10C displays the differential expression of transcription factors between the two groups. We further analyzed the correlation between these transcription factors and CRGs at the bulk transcriptome level, and found that REL, TGIF1, NFIA, IRF1, BCLAF1, NFE2L2, STAT3, KLF9, MAFF, NFIC, MAFB, KLF10, KLF6, CEBPB, KLF4, and other transcription factors are positively correlated with most CRGs, suggesting a potential network of positive regulation. The remaining transcription factors are negatively correlated with most CRGs.