Identification of a Diagnostic Gene Signature Associated with Centrosome Amplification in Pressure Injuries: A Cross-Sectional Transcriptome and Machine Learning Study

Introduction

Pressure injuries (PI) are severe skin and soft tissue lesions that result from prolonged mechanical deformation of tissues compressed between rigid internal structures (eg, bones or tendons) and external support surfaces.^1,2 Approximately 70% of PI occur at the sacrum, ischial tuberosity, or greater trochanters, whereas 15–25% occur in the lower extremities, most commonly at the heel or lateral malleolus. Although these sites are most frequently affected, PI can develop in any anatomical region exposed to sustained pressure, including the elbows, ears, nose, chest, and back.³ PI cause substantial harm to patients, often leading to severe pain, reduced quality of life, and increased mortality, and they impose a considerable economic burden on health care systems.⁴ In Australia, the total public-sector cost attributable to Pls exceeded AUD 9.1 billion in 2020.⁵ Similarly, data from the United States indicate that from 2009 to 2019, the annual prevalence of Pls continued to increase, with mean hospitalization costs rising from USD 69,499 to USD 102,939.^6–9 Current clinical management strategies for PI include surgical interventions (eg, debridement), adjunctive therapies (eg, growth factors), nutritional support, palliative wound care aimed at pain relief, and management of bacterial infections.⁷ Furthermore, nutritional supplementation with collagen peptides, vitamin C, and arginine has been shown to promote PI healing.¹⁰ In addition, the perceived level of collaboration and synergy between patients and nurses significantly influences the effectiveness of PI prevention strategies.¹¹ Current research on PI risk factors generally falls into two broad domains: mechanical boundary conditions (eg, mechanical loading, shear forces, and pressure) and individual susceptibility and tissue tolerance (eg, physiological function and tissue geometry). However, no definitive causal relationship has been established between these two domains and PI development. Therefore, early intervention for patients at risk of PI—particularly older inpatients with chronic diseases—is essential to reduce PI incidence.

The centrosome is a non-membranous organelle comprising two centrioles surrounded by pericentriolar material, and it serves as the primary microtubule-organizing center in animal cells.¹² Moreover, the centrosome plays a critical role in spindle assembly and orientation during mitosis and in cilia formation during interphase.¹² When cells contain more than two centrosomes or exhibit abnormally enlarged centrosome volume, including both numerical and structural abnormalities, this condition is termed centrosome amplification (CA).¹³ Currently, the role of CA in cancer cell proliferation is well recognized, and centrosomes are being actively investigated as potential molecular targets for cancer therapy. In contrast, studies examining CA in non-neoplastic diseases remain limited.^14,15 Ischemia–reperfusion (I/R) injury–induced oxidative stress and subsequent inflammatory responses are widely regarded as key contributors to PI development. I/R injury is characterized by exacerbated hypoxic damage following the restoration of blood flow and oxygen supply, which disrupts intracellular redox homeostasis and leads to excessive production of reactive oxygen species (ROS) in affected tissues.^16,17 ROS can induce DNA damage and centrosome overduplication, thereby triggering CA, and the prevalence of CA has been shown to increase with aging and oxidative stress.¹⁸ Furthermore, supernumerary centrosomes generated during unscheduled polyploidization or aberrant centriole biogenesis can activate NF-κB signaling, promoting sterile inflammation.¹⁹ Proinflammatory cytokines can also induce CA, which further amplifies cytokine release, creating a self-perpetuating vicious cycle.²⁰ CA may also exacerbate local cell death by inducing mitotic catastrophe. Collectively, these findings suggest that CA may play a critical role in PI development and progression by disrupting cell cycle progression, impairing normal mitosis, exacerbating genomic instability, and fostering a chronic, nonhealing wound microenvironment. However, the specific role of CA in PI pathogenesis and its underlying regulatory network remain unclear.

Although PI primarily manifest as localized tissue pathology, existing research indicates that systemic inflammatory responses contribute significantly to their development.²¹ For example, previous studies have reported significantly elevated serum FABP4 levels in patients with recurrent PI following spinal cord injury, whereas FABP3 levels were lower in patients without PI. Concurrently, circulating interleukin-13 and VEGF-A levels were markedly reduced in patients with recurrent PI. These findings suggest that alterations in peripheral blood molecular expression may reflect PI-associated pathophysiological processes.^22,23 However, most existing studies have focused on tissue-level mechanisms, and systemic molecular investigations based on peripheral blood remain limited.

Accordingly, the primary objective of this study was to collect peripheral blood samples from patients with PI for transcriptome sequencing, integrate CA-related genes, and apply differential expression analysis, machine learning, and receiver operating characteristic (ROC) curve analysis to identify key CA-associated genes in PI. These genes or their expression signatures may serve as potential circulating biomarkers for auxiliary early diagnosis or risk stratification in high-risk populations. The secondary objective was to conduct multidimensional analyses—including functional enrichment, immune infiltration, regulatory mechanism exploration, and drug prediction—to preliminarily investigate the biological roles of these key genes and to generate new hypotheses regarding the molecular mechanisms underlying PI development.

Materials and Methods

Data Collection

Patients were included in the analysis if they (1) met the diagnostic criteria for PI, (2) were between 18 and 75 years of age, and (3) had PI located at any anatomical site (eg, the sacrococcygeal region, heel, or ischium). The control group was matched to the PI group by age and sex. All participants provided informed consent and voluntarily participated in the study.

Patients were excluded (1) if they had severe respiratory disease, severe cardiovascular or cerebrovascular disease, severe systemic infection, active autoimmune disease, or coagulation disorders; (2) if they were receiving hospice or palliative care; (3) if they were using medications that could influence study outcome measures (eg, glucocorticoids, immunosuppressants, or biologic agents); (4) if they had severe mental illness or cognitive impairment; or (5) if they were concurrently enrolled in another clinical trial. In addition, control participants were excluded if they had poorly healing wounds or superficial trauma or had known risk factors for PI (eg, long-term bed rest, spinal cord injury, or severe malnutrition).

Forty samples were collected from The First Affiliated Hospital of Henan University of Chinese Medicine and Henan Provincial People’s Hospital, comprising 20 patients with PI and 20 healthy control whole-blood samples. Among them, 30 samples (15 PI and 15 control) were used for transcriptomic differential expression analysis, and an additional 10 newly collected samples were used for experimental validation by quantitative reverse transcription polymerase chain reaction (RT-qPCR). The study was conducted in accordance with the Declaration of Helsinki and was approved by the Ethics Committee of The First Affiliated Hospital of Henan University of Chinese Medicine (protocol code 2025HL-511; approval date: August 05, 2025). Written informed consent was obtained from all participants.

In addition, 411 centrosome amplification–related genes (CARGs) were obtained from the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (Table S1).

Transcriptome Sequencing and Data Preprocessing

Total RNA was extracted and purified from 30 whole-blood samples using TRIzol reagent (Invitrogen, CA, USA). RNA concentration and integrity were subsequently assessed. Samples were considered eligible for downstream analyses if they met the following criteria: RNA concentration > 50 ng/μL; RNA integrity number > 7.0; OD 260/280 ratio > 1.8; and total RNA amount > 1 μg.

Poly(A) RNA was then isolated from 1 μg of total RNA through two rounds of purification. Following fragmentation, poly(A) RNA was sheared into short fragments using the Magnesium RNA Fragmentation Module (NEB; cat. no. E6150, USA) at 94°C for 5–7 min. The fragmented RNA was then reverse-transcribed into cDNA using SuperScript™ II Reverse Transcriptase (Invitrogen; cat. no. 1896649, USA), followed by PCR amplification. The mean insert size of the final cDNA library was approximately 300 ± 50 base pairs (bp). Paired-end sequencing (PE150) was subsequently performed on the Illumina NovaSeq 6000 platform.

After sequencing, low-quality reads were filtered using Fast QC (https://directory.fsf.org/wiki/FastQC), and overall sequencing quality was assessed with MultiQC (https://seqera.io/multiqc/). Clean reads were then aligned to the human reference genome (Homo sapiens, GRCh38) using HISAT2 (https://ccb.jhu.edu/software/hisat2).

Gene expression levels were quantified as FPKM using StringTie (https://ccb.jhu.edu/software/hisat2).²⁴ PCA was performed using the plotPCA function in the DESeq2 package (version 1.42.0) to evaluate the global mRNA expression differences between the PI and control groups.

Differential Expression Analysis

Differentially expressed genes (DEGs) between the PI and control groups were identified using DESeq2 with thresholds set at (|log₂ FC| > 1 and adjusted P < 0.05).^25,26 DEG visualization was performed using volcano plots and heatmaps generated with the ggplot2 (version 3.4.4) and ComplexHeatmap (version 2.18.0) packages. The heatmap provided a focused representation of the most pronounced gene expression changes.

Identification and Functional Analysis of Candidate Genes

Candidate genes were identified by intersecting DEGs with CARGs using the ggvenn package (version 1.4.9).²⁷ Subsequently, GO and KEGG analyses were conducted using clusterProfiler (version 4.10.0), with an adjusted P < 0.05. Protein–protein interaction (PPI) networks were constructed using the STRING database (https://string-db.org/), applying a confidence score > 0.4 (Homo sapiens).²⁸ After excluding outlier genes, the PPI network was visualized using Cytoscape (version 3.10.2). This visualization step was crucial for understanding the complex relationships among the candidate genes and identifying potential hubs or bottlenecks in the network, which could serve as promising targets for further investigation and therapeutic intervention.

Machine Learning

Based on the identified candidate genes, the LASSO, SVM-RFE, and Boruta algorithms were applied to screen feature genes from the transcriptome sequencing data.²⁹ LASSO analysis was conducted using glmnet (version 4.1–8) with a binomial model and 10-fold cross-validation, and features corresponding to the minimum lambda value were retained. SVM-RFE was conducted using the caret package (version 6.0–94), with feature selection defined at the point of maximal model prediction accuracy.³⁰ Additionally, the Boruta algorithm was implemented using the Boruta package (version 8.0.0). Common feature genes were identified by intersecting the gene sets obtained from the three algorithms.

Identification of Key Genes and Nomogram Construction

To evaluate the diagnostic performance of the common feature genes for PI, ROC curves were generated using pROC (version 1.18.5), and AUC was calculated. Key genes with an AUC > 0.7 were selected.^31,32 Based on these genes, a nomogram was constructed to predict PI risk using the rms package (version 6.5.0). Afterward, some curves were plotted. Calibration curves were generated with regplot, and goodness of fit was evaluated using the Hosmer–Lemeshow (HL) test, with P > 0.05 indicating acceptable calibration. ROC curves were replotted using pROC, while decision curve analysis and clinical impact curves (CICs) were generated using the rmda package (version 1.6).

Chromosomal Localization and Tissue Distribution of Key Genes

To visualize the chromosomal localization of the identified key genes, the OmicCircos package (version 1.40.0) was used.³³ This approach provided an overview of their genomic distribution across human chromosomes. To further investigate the biological relevance of these genes in PI development, specifically in the context of CA in human skin tissues and cells, the Bgee database (https://Bgee.org/) was queried to examine their expression patterns across various human skin tissues. This analysis helped identify tissue-specific expression profiles, offering additional insights into the potential involvement of these genes in PI pathogenesis.

Functional Analysis of Key Genes

GeneMANIA (https://genemania.org/) was used to identify genes functionally associated with the key genes. Additionally, gene set enrichment analysis (GSEA) was performed using the transcriptomic data. Specifically, Spearman correlation coefficients between each key gene and all other genes were calculated and ranked in descending order using the stats package (version 4.3.1) (https://www.rdocumentation.org/packages/stats). The KEGG gene set collection (c2.cp.kegg.v7.5.1.symbols.gmt) was obtained from MSigDB. GSEA was then conducted using clusterProfiler, with statistical significance defined as an adjusted P < 0.05 and |NES| > 1. To further characterize pathway alterations during PI progression, the Hallmark gene set collection (h.all.v2024.1.Hs.symbols.gmt) was then downloaded from MSigDB and used as the background gene set.³⁴ In the transcriptome dataset, the GSVA package (version 1.50.0) was applied to calculate the GSVA scores for all pathways in both the PI and control groups based on the expression matrix.³⁵ Differential analysis of GSVA scores between groups was then performed using the limma package (version 3.58.1), with thresholds set at adjusted P < 0.05 and |log2FC| > 1. Enrichment of significantly altered Hallmark pathways in each group was visualized using bar plots.

Immunological Feature Analysis

The CIBERSORT algorithm (version 0.1.0) was used to estimate the relative infiltration of 22 immune cell types in each sample.^36,37 Samples with CIBERSORT output P values of >0.05 were excluded from subsequent analyses. Differences in immune cell infiltration between the PI and control groups were assessed using the Wilcoxon rank-sum test, and immune cell types with significantly different infiltration levels were identified based on an adjusted P < 0.05. Correlation analyses were performed using the psych package to examine relationships among differentially infiltrating immune cells and to evaluate associations between these immune cells and key genes (|r| > 0.30, P < 0.05).³⁸ To further explore interactions between key genes and immune-related factors, multiple categories of immune molecules—including chemokines, receptors, immunosuppressive factors, immunostimulatory factors, MHC genes, and 17 inflammatory factors—were curated from the TISIDB database (http://cis.hku.hk/TISIDB/index.php) and relevant literature. Spearman correlation analysis was then conducted using psych to assess associations between key genes and these immune factors (|r| > 0.30, P < 0.05).

Regulatory Network Analysis

NetworkAnalyst was used to predict transcription factors (TFs) targeting the key genes, and a TF–key gene interaction network was constructed. Putative TF binding sites were identified using the JASPAR database (https://jaspar.elixir.no/). In addition, miRDB (https://mirdb.org/) and TargetScan (https://www.targetscan.org/vert_80/) were used to predict miRNAs targeting the key genes, and overlapping results from both databases were defined as key miRNAs. Subsequently, miRNet (https://www.mirnet.ca/) was used to predict lncRNAs interacting with these key miRNAs. The resulting lncRNA–miRNA–mRNA regulatory relationships were integrated and visualized using Cytoscape to construct a comprehensive regulatory network.

Disease Association Analysis

To investigate associations between key genes and PI-related disease genes, PI-related genes were first identified by querying the GeneCards database using the term “pressure injuries”. The top 20 PI-related genes with the highest relevance scores were selected for subsequent analysis. Correlation analysis was then performed using the transcriptomic data to assess associations between the key genes and these top 20 PI-related genes, with thresholds set at |r| > 0.30 and P < 0.05. To further identify diseases associated with the key genes, the DisGeNET database (https://disgenet.com/) was queried. The top 10 diseases showing the strongest interaction scores with the key genes were selected to construct a disease–gene interaction network.

Drug Prediction and Molecular Docking

The DGIdb database (https://www.dgidb.org/) was used to identify drugs associated with the key genes, with significance defined as P < 0.05. The top 15 candidate drugs, ranked by P value, were selected to construct and visualize a biomarker–drug interaction network using the ggsankey package (version 0.0.99999).

To further assess the binding affinity of these drugs to the key genes, we selected the drug that interacted with all key genes among the top 15 candidates for molecular docking studies. The protein crystal structures of the key genes (receptors) were obtained from UniProt, and the 3D molecular structures of the candidate drug (ligand) were retrieved from PubChem. Molecular docking was conducted using AutoDock2 (https://autodock.scripps.edu/), and binding energies were calculated.³⁹ Docking results were visualized using PyMOL (version 2.5). Binding energies below −5 kcal/mol indicated strong binding affinity, suggesting effective molecular interactions between the drug and the target proteins.

RT-qPCR Validation

To validate the expression levels of key genes identified by transcriptomic analysis, RT-qPCR was performed using an independent clinical cohort. Total RNA was extracted from 10 whole-blood samples using TRIzol reagent (Ambion, USA) according to the manufacturer’s instructions. RNA concentration was measured using a NanoPhotometer N50 (Implen, Munich, Germany). cDNA was synthesized using the SureScript First-Strand cDNA Synthesis Kit, and reverse transcription was performed on an S1000™ Thermal Cycler (Bio-Rad, USA). Primer sequences are listed in Table S2. qPCR was conducted using the CFX Connect Real-Time Quantitative Fluorescence PCR Instrument (Bio-Rad). Relative mRNA expression levels were calculated using the 2^−ΔΔCT method. RT-qPCR results were exported to Microsoft Excel and subsequently visualized using Graphpad Prism 5 (https://www.graphpad.com/).

Statistical Analysis

All statistical analyses were performed using R software (version 4.3.1). Differences between the two groups were assessed using the Wilcoxon rank-sum test, with statistical significance defined as P < 0.05. Multiple comparisons were corrected using the Benjamini–Hochberg method. For RT-qPCR data, between-group comparisons were conducted using a t-test, with P < 0.05 considered statistically significant.

Results

Quality Assessment of Transcriptomic Data

Quality control analysis using FastQC demonstrated that the majority of sequencing reads exhibited high-quality scores, with values concentrated near 40, indicating excellent overall sequencing data quality (Figure S1a). FPKM density distributions showed that gene expression values for most samples were predominantly within the range of 0.1–100, reflecting a consistent and reliable expression profile across samples (Figure S1b).

PCA identified one outlier control sample (A2). Before outlier removal, PC1 and PC2 accounted for 33.6% and 23.2% of the total variance, respectively. Following the removal, PC1 and PC2 accounted for 33.9% and 24.1% of the variance, respectively (Figure S1c). Subsequently, samples from the two groups showed tighter clustering, with clearer separation of expression patterns. A total of 15 PI samples and 14 control samples were retained for subsequent analyses.

Acquisition and Function Analysis of Candidate Genes

Transcriptome analysis identified 2431 DEGs, including 988 upregulated and 1443 downregulated genes in the PI samples (Figure 1a and b). Intersection of these DEGs with 411 CARGs yielded 50 candidate genes (Figure 1c). Functional enrichment analysis revealed significant associations with 296 GO terms, comprising 249 biological processes (BPs), 29 cellular components (CCs), and 18 molecular functions (MFs). Prominent terms included “centrosome cycle” (BP), “spindle” (CC), and “microtubule binding” (MF) (Figure 1d and Table S3). KEGG pathway analysis demonstrated enrichment in eight pathways, notably the “cell cycle” and “MAPK signaling pathway” (Figure 1e). After removing 12 outlier genes, a PPI network comprising 53 significant interactions was constructed (Figure 1f). The five genes with the highest degrees of connectivity were CDK5RAP2, GEN1, RAD51AP1, ALMS1, and CEP250.

Figure 1 Identification and functional analysis of candidate genes. (a) Volcano plot showing differentially expressed genes (DEGs) between the pressure injury (PI) and control groups. The x-axis represents log2(PI/Control), and the y-axis represents −log10(adjusted P value). Genes closer to the upper corners indicate greater statistical significance. Dashed lines on the x-axis indicate log2FC thresholds of ±1, and the dashed line on the y-axis denotes −log10(0.05). (b) Heatmap illustrating expression patterns of the top DEGs. Higher expression levels are shown in red, and lower expression levels are shown in light green. (c) Identification of DE-CARGs. (d) Bar plot of significant Gene Ontology terms among the candidate genes. (e) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of candidate genes. The left panel displays genes ranked by log fold change, with pink indicating upregulation and green indicating downregulation; color intensity reflects the magnitude of change. The right panel presents KEGG pathways ranked by statistical significance. Dot size corresponds to adjusted P values, with smaller dots indicating stronger enrichment. (f) Protein–protein interaction network of candidate genes. Nodes represent genes, and edges represent interactions between the encoded proteins.

Identification of GADD45A, LFNG, and DUSP13 as Key Genes in PI

From the 50 candidate genes, three feature genes—GADD45A, LFNG, and DUSP13—were identified through a series of analyses. Specifically, 6 genes were selected using LASSO regression (lambda.min = 0.044; Figure 2a), 4 genes were identified using SVM-RFE (Figure 2b), and 17 genes were confirmed using the Boruta algorithm (Figure 2c). Intersection of the results from all three methods yielded GADD45A, LFNG, and DUSP13 (Figure 2d). ROC analysis showed that all three genes exhibited strong diagnostic performance for PI, with AUC values exceeding 0.8 (Figure 2e). Accordingly, GADD45A, LFNG, and DUSP13 were identified as key genes associated with PI, indicating potential diagnostic and therapeutic relevance. Consistent with the transcriptomic findings, RT-qPCR validation showed significantly higher expression of GADD45A and DUSP13 and significantly lower expression of LFNG in the PI samples compared with the controls (P < 0.05; Figure 2f).

Figure 2 Machine learning–based identification of key genes. (a) LASSO coefficient profiles showing the trajectories of candidate gene coefficients. The x-axis indicates the log(lambda), and the y-axis represents coefficient values; the dashed line marks the minimum cross-validation error (lambda.min). (b) Accuracy variation curve of the SVM-RFE model obtained by cross-validation. The y-axis represents classification accuracy corresponding to different numbers of selected genes. (c) Feature importance analysis using the Boruta algorithm. Box plots show the distribution of importance scores for real features compared with shadow features (randomly permuted). (d) Venn diagram showing the intersection of genes selected by the three machine learning methods. (e) ROC curves evaluating the diagnostic performance of the identified key genes. The x-axis represents the false-positive rate (1 − specificity), and the y-axis represents the true-positive rate (sensitivity). (f) RT-qPCR validation of key gene expression levels, **p < 0.01.

Construction of a High-Performance Nomogram

Based on the three key genes, a nomogram was constructed by assigning weighted point values to each gene. Higher cumulative scores were associated with an increased risk of PI (Figure 3a). Model calibration confirmed the accuracy of the model, with a P value of 0.919 in the HL test (Figure 3b). The nomogram achieved a perfect AUC of 0.99, indicating outstanding predictive performance (Figure 3c). The decision curve analysis demonstrated that the nomogram provided a greater net clinical benefit than any individual predictor (Figure 3d). CIC analysis further demonstrated that when the risk threshold exceeded 0.6, the predicted PI cases closely aligned with the observed outcomes, supporting the high clinical prediction efficiency of the model (Figure 3e). Overall, these findings indicate strong predictive performance in this exploratory analysis; however, external validation in independent cohorts is required.

Figure 3 Nomogram construction and validation. (a) Nomogram integrating key gene expression levels. The vertical axis lists the variables (gene expression levels), each associated with a corresponding score scale. The horizontal axis represents the score range for each variable. Red dots indicate the scores for individual variables, and vertical lines aggregate these values to generate the total score (Total Points). (b) Calibration curve comparing predicted and observed probabilities of PI. The x-axis represents predicted probability, and the y-axis represents observed probability. (c) ROC curve evaluating overall model performance. (d) Decision curve analysis assessing clinical net benefit across threshold probabilities. The x-axis represents the decision threshold, and the y-axis indicates net benefit. Curves represent individual gene models (GADD45A, red; DUSP13, light blue; LFNG, purple) and the combined model (dark blue). (e) Clinical impact curve showing the number of individuals predicted to be at high risk per 1000 individuals across threshold probabilities. Results are based on internal validation and are intended for exploratory purposes; external validation is required before clinical application.

Figure 4 Genomic localization and tissue expression patterns of key genes. (a) Chromosomal localization of key genes. Bands in the circular plot represent individual chromosomes, with gene positions indicated. (b) Expression profiles of key genes across various skin tissues. The vertical axis shows gene locations, and the horizontal axis represents gene expression scores, with higher values indicating higher expression. The color gradient from light to dark reflects −log10(false discovery rate [FDR]), with darker colors indicating a lower FDR and higher statistical significance. (c–e) Correlations between key gene expression and different skin cell types: (c) GADD45A, (d) LENG, and (e) DUSP13.

Genomic Localization and Tissue Distribution of Key Genes

Chromosomal localization analysis revealed that GADD45A is located on chromosome 1, LFNG on chromosome 7, and DUSP13 on chromosome 10 (Figure 4a). Expression profiling across human skin tissues revealed that GADD45A and LFNG are highly expressed in various skin regions, including leg and abdominal skin (Figure 4b). In skin tissues, GADD45A expression was strongly correlated with endothelial cells and squamous epithelial cells (Figure 4c). For LFNG, a notable correlation was observed between keratinocytes and squamous epithelial cells (Figure 4d). Similarly, DUSP13 expression was notably correlated with keratinocytes and squamous epithelial cells in the skin tissues (Figure 4e).

Functional and Pathway Analyses of Key Genes

A total of 20 genes functionally associated with the key genes were identified. Representative interactions included LFNG with B4LGA1 and GADD45A with TAF6, highlighting involvement in processes such as transferase activity (transferring hexosyl groups), the Notch signaling pathway, and cell cycle arrest (Figure 5a). GSEA demonstrated that GADD45A was significantly enriched in 70 pathways, including “pathogenic Escherichia coli infection” and “regulation of actin cytoskeleton” (Figure 5b and Table S4). LFNG was enriched in 22 pathways, notably the “ribosome” and “spliceosome” pathways (Figure 5c and Table S4), whereas DUSP13 was enriched in 63 pathways, including “endocytosis” and the “B cell receptor signaling pathway” (Figure 5d and Table S4). Notably, all three key genes were co-enriched in the “spliceosome”, “Huntington disease”, and “ubiquitin-mediated proteolysis” pathways, suggesting shared involvement in critical cellular regulatory processes. GSVA further revealed that several pathways, including “apoptosis” and “adipogenesis”, were significantly activated in the PI samples, whereas “Notch signaling” and “KRAS signaling DN” were significantly suppressed (Figure 5e).

Figure 5 Functional characterization of key genes. (a) GeneMANIA network illustrating functional associations among the key genes. Each node represents a gene, and colored edges indicate different types of relationships. Different colors within nodes denote the biological functions or pathways associated with each gene. (b–d) GSEA plots for hallmark pathways associated with each key gene: (b) GADD45A, (c) LENG, and (d) DUSP13. Each GSEA plot comprises three components: (1) an enrichment score curve, in which the x-axis represents all genes ranked by correlation with the key gene (positively correlated genes are shown on the left [red region], and negatively correlated genes are shown on the right [blue region]) and the y-axis represents the running enrichment score; the maximum deviation from zero corresponds to the enrichment score, with genes preceding this peak defined as core enrichment genes; (2) vertical bars indicating the positions of genes from the tested pathway within the ranked gene list; and (3) a ranked gene distribution plot illustrating the correlation metric (eg, signal-to-noise ratio) across all genes. (e) Positive and negative values of the enrichment analysis t value of the PI group and the control group indicate the relative difference in the expression level of the gene set between the two clusters. Positive t values indicate higher pathway activity in the PI group, whereas negative t values indicate higher activity in the control group.

Differences in Immune Cell Infiltration Between the PI and Control Samples

Immune infiltration profiles derived from transcriptomic data are shown in Figure 6a. Significant differences in immune cell infiltration between the PI and control samples were observed for seven immune cell types (Figure 6b). Specifically, memory B cells showed reduced infiltration in PI samples, whereas neutrophils exhibited significantly increased infiltration (P < 0.05). Correlation analysis demonstrated a strong negative association between CD8⁺ T cells and neutrophils (r = −0.87, P < 0.05) and a strong positive correlation between CD8⁺ T cells and memory B cells (r = 0.75, P < 0.05) (Figure 6c). Moreover, neutrophil infiltration was significantly positively correlated with GADD45A (r = 0.67, P < 0.01) and DUSP13 (r = 0.72, P < 0.01) and negatively correlated with LFNG (r = −0.66, P < 0.01) (Figure 6c and d). These findings suggest that neutrophils may play a significant role in the immune modulation associated with PI. Furthermore, correlation analysis revealed that GADD45A, DUSP13, and LFNG were significantly associated with various immune-related molecules, including chemokines (eg, XCL2, CCL5), receptors (eg, CXCR1, CCR9), immunoinhibitory molecules (eg, CD96, TIGIT), immunostimulatory molecules (eg, CD2, TNFRSF25), and MHC genes (eg, HLA-DRA, TAPBP). In contrast, only GADD45A and LFNG showed significant correlations with selected inflammatory factors, such as IL10 and IL6. Notably, LFNG often exhibited correlation patterns opposite to those of GADD45A and DUSP13 with immune-related factors, suggesting a complex and potentially antagonistic regulatory relationship between these key genes and the chemokine–immune network (Figure 6e).

Figure 6 Immune infiltration landscape in PI. (a) Stacked bar plot showing the relative composition of immune cell types across samples; different colors represent distinct immune cell populations. (b) Boxplots comparing significantly altered immune cell populations; different colors represent group membership. (c) Correlation analysis between the key gene expression levels and differentially infiltrating immune cells. Correlation coefficients are shown in the lower-left panels, and statistical significance is shown in the upper-right panels; color represents the size, and asterisks represent significance levels, *p < 0.05, **p < 0.01, ***p < 0.001. (d) Correlation bar plots showing the relationship between key genes (GADD45A, LFNG, DUSP13) and differentially infiltrating immune cells. Bar length represents the correlation coefficient (r), with direction indicating positive or negative correlations. Dot color reflects P values, with red dots indicating statistically significant correlations. (e) Heatmap showing correlations between key genes (LFNG, GADD45A, DUSP13) and immune-related factors, including chemokines, receptors, immunoinhibitory molecules, immunostimulatory molecules, MHC genes, and inflammatory factors. Green indicates positive correlations, and blue indicates negative correlations, *p < 0.05, **p < 0.01.

Potential Regulatory Mechanisms Underlying PI

A total of 11, 16, and 13 TFs were predicted to regulate GADD45A, LFNG, and DUSP13, respectively. Among these, JUN was identified as a common TF targeting all three key genes (Figure 7a). Analysis of TF binding sites revealed a higher frequency of T at site 5 (Figure 7b). In addition, 614 regulatory relationships were identified, enabling the construction of a comprehensive lncRNA–miRNA–mRNA network (3 mRNAs, 48 key miRNAs, 563 lncRNAs). Representative regulatory axes included CYTOR–hsa–mir–125a-5p–LFNG and DLEU1–hsa-mir-671-5p–DUSP13, among others (Figure 7c). Correlation analysis further demonstrated that 10 PI-related genes (eg, IL10, TLR4, and TP53) were negatively correlated with LFNG (r < −0.3, P < 0.05) and positively correlated with GADD45A (r > 0.3, P < 0.05). In contrast, DUSP13 exhibited a notable negative correlation with TP53 only (r < −0.3, P < 0.05; Figure 7d). Disease–gene network analysis revealed that several diseases, including pauciraticular juvenile arthritis and ovarian neoplasm, were associated with GADD45A, suggesting potential links between these diseases and GADD45A (Figure 7e).

Figure 7 Regulatory network analysis. (a) TF-mRNA regulatory network, in which red nodes represent key genes and Orange nodes represent TFs. (b) Binding site of the transcription factor JUN is shown, with the x-axis representing base sequence position and the y-axis indicating information entropy (bits). Higher values indicate greater base conservation. The height of each base (A, T, G, C) represents its relative frequency, with different colors for each base. (c) LncRNA–miRNA–mRNA regulatory network, with mRNAs (key gene) shown in red, miRNAs in blue, and lncRNAs in green. (d) Correlation analysis between key genes and PI-related genes. The x-axis lists PI-related genes, and the y-axis lists key genes (LFNG, GADD45A, DUSP13). Color intensity reflects correlation strength (green, positive [close to 1]; blue, negative [close to −1]), and asterisks denote statistically significant correlations, *p < 0.05, **p < 0.01. (e) Gene–disease association network, with key genes shown in red and associated diseases shown in blue.

Exploration of Potential Therapeutic Drugs for PI

A total of 181, 10, and 8 candidate drugs were predicted to target GADD45A, LFNG, and DUSP13, respectively. Notably, decitabine was identified as one of the three common therapeutic agents targeting all three key genes (GADD45A, LFNG, and DUSP13) (Figure 8a). To further investigate these interactions, molecular docking analyses were performed between decitabine and each of the three key genes. Docking results showed favorable binding between DUSP13 and decitabine, with a binding energy of −6.9 kcal/mol and key interactions at residues Asp188 and Pro26 (Figure 8b). GADD45A and decitabine demonstrated a binding energy of −5.5 kcal/mol, involving residues Gln134 and Lys136 (Figure 8c). Similarly, LFNG and decitabine showed a binding energy of −5.5 kcal/mol, with interactions at residues Thr122 and Ala171 (Figure 8d).

Figure 8 Therapeutic drug prediction and molecular docking. (a) Sankey diagram illustrating drug–gene interactions, with key genes on the left and predicted therapeutic drugs on the right. (b–d) Molecular docking models showing binding conformations between decitabine and each key gene protein. Molecular docking represents a preliminary computational analysis, and the results require experimental validation.

Discussion

PI are localized injuries to the skin and underlying soft tissues caused by sustained mechanical pressure and shear force, most commonly affecting individuals with prolonged immobility and substantially impairing quality of life and clinical outcomes.^40,41 CA, characterized by numerical or structural abnormalities of centrosomes, is well established as a driver of tumorigenesis; however, its role in non-neoplastic conditions, including PI, remains unclear. In this study, three CA-related genes—GADD45A, LFNG, and DUSP13—were preliminarily identified through an integrated analytical framework combining transcriptomic sequencing and multiple machine learning algorithms. A risk prediction nomogram constructed based on these genes showed strong predictive performance for PI risk, with an AUC > 0.7. Furthermore, functional enrichment, immune infiltration, and drug prediction analyses provided initial insights into the regulatory mechanisms and potential therapeutic relevance of these genes in PI pathophysiology. Experimental validation using RT-qPCR further confirmed that GADD45A and DUSP13 were significantly upregulated, whereas LFNG was significantly downregulated, in peripheral blood samples from patients with PI compared with healthy controls. Collectively, these findings provide preliminary molecular evidence for further systematic investigation into the role of centrosome-related abnormalities in PI development and progression.

GADD45A (growth arrest and DNA damage–inducible α) is a member of the GADD45 family and functions as a stress-inducible gene that is significantly upregulated under various injury conditions.⁴² It exhibits proapoptotic and cytotoxic effects, particularly in models of I/R and endoplasmic reticulum (ER) stress.^43,44 For example, in an I/R injury–induced cardiac microvascular endothelial cell model, increased GADD45A expression promotes apoptosis by activating the JNK/p38 MAPK pathway while simultaneously inhibiting angiogenesis by suppressing the STAT3/VEGF pathway, thereby exacerbating tissue damage.⁴⁵ From a miRNA regulatory perspective, GADD45A has been identified as a direct target of miR-301a-3p.⁴⁴ Under ER stress conditions, GADD45A is upregulated through the IRE1-mediated degradation of pre–miR-301a, leading to enhanced apoptosis. This mechanism highlights the critical role of GADD45A in maintaining ER homeostasis and regulating cell fate decisions. In the context of PI, our findings suggest that GADD45A may contribute to disease pathogenesis by participating in cell cycle regulation, DNA damage responses, and immune cell infiltration. Collectively, these observations suggest that GADD45A, as a core regulator of cellular stress responses, may be involved in a “maladaptive tissue damage–inflammation–repair imbalance” during PI development. Future research may therefore consider developing therapeutic strategies targeting miRNA-mediated regulation of GADD45A as a potential therapeutic target for PI treatment and other stress-related disorders.

LFNG encodes a glycosyltransferase that plays a critical role in regulating the Notch signaling pathway.⁴⁶ Specifically, LFNG modifies Notch receptors through O-fucosylation, thereby modulating Notch signal transduction and regulating critical processes, such as cell differentiation and proliferation. Dysregulated LFNG expression has been implicated in various disease contexts, underlying the importance of its regulatory function. For example, exosomes derived from bone marrow mesenchymal stem cells can deliver specific tsRNAs (such as tsRNA-15797), which enhance angiogenic capacity by negatively regulating LFNG expression.⁴⁷ This mechanism has shown therapeutic potential in ischemia-related diseases, including osteonecrosis of the femoral head.⁴⁷ Experimental validation has demonstrated that LFNG overexpression significantly weakens the proangiogenic effects of tsRNA-15797–modified exosomes, suggesting that regulating LFNG expression may present a novel therapeutic strategy for treating ischemic diseases. Machine learning–based analyses combined with clinical validation in the present study identified LFNG as a key DEG associated with CA in patients with PI, with significantly reduced expression. This finding highlights both the diagnostic potential and biological relevance of LFNG in disease contexts. Its involvement in PI progression may be mediated through its effects on cell cycle regulation, cellular proliferation, and apoptosis. It is hypothesized that abnormal LFNG expression may affect immune regulation, cell proliferation, and apoptosis through multiple pathways, thereby contributing to disease development. Collectively, these observations suggest that LFNG may serve as a potential diagnostic biomarker, prognostic indicator, and therapeutic target for clinical treatment.

DUSP13, a member of the dual-specificity phosphatase family, is essential in signal regulation across various disease states and has been implicated in cancer biology, muscle differentiation and regeneration, and cellular stress responses.⁴⁸ During the transition of muscle stem cells from proliferation to differentiation, DUSP13 functions as a key regulatory mediator.⁴⁹ Previous studies have shown that DUSP13 is a direct transcriptional target of the myogenic regulator MYOD. Its deletion results in a significant reduction in muscle regenerative capacity, while its overexpression induces premature differentiation of muscle cells. This mechanism further indicates that DUSP13 may modulate muscle regeneration by regulating MAPK signaling activity. Collectively, these findings indicate that DUSP13 is an important regulator of the cellular signaling pathways involved in both tissue regeneration and stress responses.

In this study, GSEA further revealed that the three key genes—GADD45A, LFNG, and DUSP13—were commonly enriched in several critical pathways, including SPLICEOSOME, HUNTINGTON’S DISEASE, and UBIQUITIN_MEDIATED_PROTEOLYSIS.

Abnormalities in the spliceosome can lead to mRNA splicing errors, thereby affecting protein function and intracellular signal transduction.^50,51 Although research on spliceosome function in PI remains limited, extensive evidence from cancer and neurodegenerative diseases demonstrates that spliceosome abnormalities are closely associated with cell cycle dysregulation, apoptosis, and altered immune responses. In the context of PI, aberrant spliceosome activity may influence pressure-induced cellular injury, tissue repair capacity, and immune regulation, suggesting a potentially important role in PI initiation and progression.⁵² During PI pathogenesis, sustained cellular and tissue damage may trigger protein damage, and the ubiquitin–proteasome system preserves intracellular stability by selectively degrading damaged or misfolded proteins.⁵³ During PI repair, ubiquitination processes may help clear damaged cells, repair defective proteins, and regulate immune cells and tissue remodeling. Accordingly, the ubiquitin-mediated proteolysis pathway is likely to play a central role in coordinating cellular stress responses and immune reactions in PI.

Immune infiltration analysis identified seven immune cell types that differed significantly between the PI and control samples, with all three key genes (GADD45A, LFNG, and DUSP13) showing strong correlations with neutrophil infiltration.⁵⁴ Neutrophils play a critical role in PI development and progression. Prolonged pressure leads to local tissue ischemia, and subsequent reperfusion triggers oxidative stress and inflammatory responses.⁵⁵ Neutrophils are central mediators of these responses, releasing cytokines and inflammatory mediators that can support early tissue repair. However, excessive or sustained neutrophil activation may result in secondary tissue injury and chronic inflammation, thereby exacerbating PI severity.

Drug prediction analysis further identified the epigenetic agent decitabine as a common compound targeting all three key genes, with molecular docking showing the strongest binding affinity to DUSP13 (−6.9 kcal/mol).⁵⁶ Although decitabine is primarily recognized for its antitumor activity through DNA methylation inhibition, this study is the first to suggest its potential role in mitigating oxidative stress–related damage in PI by modulating CA-associated genes such as GADD45A. These findings highlight decitabine as a promising candidate for PI. However, it should be emphasized that drug prediction based on molecular docking is inherently computational and hypothesis-generating and cannot directly establish in vivo efficacy and safety. Consequently, further experimental studies are required to determine whether decitabine can effectively modulate the functions of the identified target genes and ultimately influence the pathological progression of PI, thereby providing a robust foundation for clinical translation.

In this study, GADD45A, LFNG, and DUSP13 were preliminarily identified as key CA–associated genes in PI. Pending large-scale validation, their expression profiles in peripheral blood may serve as candidate biomarkers to support the early identification and dynamic risk monitoring of high-risk patients, particularly by nursing and nutritional support teams. These molecular markers could complement existing clinical and nursing assessment frameworks by enabling a more precise and individualized allocation of preventive resources. From a therapeutic perspective, computational analyses suggest that decitabine may target this gene network, and the associations of the key genes provide preliminary directions for developing novel intervention strategies.

Several limitations of this study should be acknowledged. First, the relatively small sample size may limit the external validity and generalizability of the findings to broader populations or diverse clinical settings. Second, although peripheral blood sampling is clinically practical, transcriptomic analyses based on blood may not fully capture the localized molecular changes occurring in skin or wound tissues at PI sites. Additionally, the construction of predictive models using high-dimensional transcriptomic data in a small cohort may increase the risk of overfitting. This underscores the need for larger independent validation cohorts to confirm the robustness and reliability of the results. Furthermore, the cross-sectional study design, which compares gene expression profiles between patients with PI and healthy controls at a single time point, allows for biomarker identification but cannot preclude causal inference or limit the ability to assess dynamic gene expression changes throughout PI development and progression. Finally, the absence of gene-knockout experiments and animal model validation limits a comprehensive understanding of the mechanistic roles of the identified key genes in PI. Future studies will therefore focus on independent external validation using multicenter, large-sample cohorts and on comparative analyses integrating both wound tissue and peripheral blood samples to improve biological relevance. In addition, prospective longitudinal study designs will be implemented to dynamically track gene expression changes, thereby clarifying the temporal roles of key genes during PI onset, progression, and healing. Gene-editing approaches and animal models will be employed to functionally validate the roles of GADD45A, LFNG, and DUSP13 and to experimentally assess the therapeutic potential of candidate drugs, such as decitabine, with the goal of systematically advancing clinical translation. Despite these limitations, this study provides new perspectives on the molecular mechanisms, early diagnosis, and potential therapeutic targets for PI.

Conclusion

By integrating CA-related genes with peripheral blood transcriptomic data from patients with PI, this study preliminarily identified GADD45A, LFNG, and DUSP13 as potential key genes and validated their differential expression using RT-qPCR. The predictive model constructed based on these genes demonstrated moderate classification performance in the validation cohort, suggesting potential utility for risk assessment. Further exploratory analyses, including immune infiltration profiling, regulatory network analysis, and molecular docking, provided initial insights into PI pathogenesis and potential intervention strategies. However, these findings remain exploratory and require confirmation through large-scale clinical trials, systematic longitudinal studies, and functional experiments to fully clarify their clinical applicability and biological significance.