Screening for candidate diagnostic biomarkers for keloid using a machine learning algorithm
A total of 283 heme metabolism-related genes were included in this study. Differential expression analysis of the GSE44270 dataset comparing keloid and normal skin tissues identified 25 significantly differentially expressed genes (Figure 1A and Supplemental Table S3). To further screen disease-associated biomarkers, LASSO regression identified 9 candidate genes (Figure 1B,C and Supplemental Table S3), while the random forest (RF) algorithm selected 11 genes with high predictive importance (Figure 1D and Supplemental Table S3). The overlap between the LASSO and RF results was visualized using a Venn diagram, yielding six core biomarkers, namely, FLVCR1, TMCC2, EIF2AK1, XK, HPX, and KEL (Figure 1E and Supplemental Table S3). Receiver operating characteristic (ROC) analysis in the GSE44270 cohort demonstrated favorable diagnostic performance for all six biomarkers, with AUC values of 0.8016 for FLVCR1, 0.7063 for TMCC2, 0.7817 for EIF2AK1, 0.7460 for XK, 0.7500 for HPX, and 0.7857 for KEL (Figure 1F). Based on these six biomarkers, a diagnostic nomogram for keloid was subsequently constructed using the rms package in R (Figure 1G).

Figure 1: Identification of candidate heme metabolism-related genes associated with keloid using machine-learning algorithms. (A) Box plot illustrating the differential expression of heme metabolism-related genes between keloid and normal tissues. (B,C) LASSO logistic regression analysis for screening candidate diagnostic markers. (D) Candidate biomarkers selected by the RF algorithm. (E) Venn diagram showing the overlapping genes identified by the two machine-learning algorithms. (F) ROC curve analysis evaluating the diagnostic performance of the candidate biomarkers. (G) Nomogram for keloid prediction based on the six-gene signature. Abbreviations: LASSO, least absolute shrinkage and selection operator; RF, random forest; ROC, receiver operating characteristic. Statistical significance: ns, P > 0.05; *, P < 0.05; **, P < 0.01; ***, P < 0.001; and ****, P < 0.0001. Please click here to view a larger version of this figure.
The predictive performance of the diagnostic nomogram was assessed in both the training cohort (GSE44270) and the validation cohort (GSE7890). The model demonstrated excellent diagnostic accuracy, achieving AUC values of 0.984 (95% CI: 0.950–1.000) and 0.922 (95% CI: 0.806–1.000), respectively (Figure 2A,D). To further evaluate the robustness and potential overfitting of the six-gene diagnostic signature, additional internal validation analyses were performed in the discovery cohort (GSE44270). Fivefold cross-validation demonstrated consistent discriminative performance across subsets, with a mean AUC of 0.925, indicating that the model maintained stable classification performance despite variations in the training samples. Furthermore, bootstrap validation with 1,000 resampling iterations yielded a mean AUC of 0.930 (95% CI: 0.794–1.000). After adjustment for potential optimism caused by the limited sample size, the optimism-corrected AUC remained 0.930, suggesting that the diagnostic performance of the six-gene signature was relatively stable after internal validation. Furthermore, decision curve analysis (DCA) suggested that the nomogram exhibited a higher potential net benefit than alternative diagnostic strategies across a range of threshold probabilities, although these findings should be interpreted with caution given the limited sample size (Figure 2B,E). Additionally, keloid samples exhibited significantly higher risk scores than healthy controls in both the training and validation cohorts (Figure 2C,F), further demonstrating the diagnostic model's stability and reliability.

Figure 2: Validation of the nomogram for keloid prediction. (A) ROC curve evaluating the predictive performance of the nomogram in the GSE44270 dataset. (B) DCA evaluating the clinical utility of the nomogram in GSE44270. (C) Risk score distribution comparing keloid and healthy samples in GSE44270. (D) ROC curve evaluating the predictive performance of the nomogram in the independent GSE7890 dataset. (E) DCA evaluating the clinical utility of the nomogram in GSE7890. (F) Risk score distribution comparing keloid and healthy samples in GSE7890. Abbreviations: ROC = receiver operating characteristic; DCA = decision curve analysis. Please click here to view a larger version of this figure.
Diagnostic biomarkers are associated with the immune characteristics of keloid
To explore the relationship between the six diagnostic biomarkers and the immune microenvironment, correlation analysis was performed to assess the associations between biomarker expression and immune cell infiltration. The results revealed that all six biomarkers were significantly associated with multiple infiltrating immune cell populations (Figure 3A). Specifically, FLVCR1 expression was negatively associated with T follicular helper cells (Figure 3B). TMCC2 showed positive correlations with natural killer cells and activated dendritic cells, whereas it was negatively correlated with immature dendritic cells and immature B cells (Figure 3C–F). In addition, EIF2AK1 expression was negatively associated with CD56dim natural killer cells (Figure 3G), whereas XK was negatively associated with eosinophils (Figure 3H).

Figure 3: Correlation between candidate heme metabolism-related genes and immune cell infiltration. (A) Heatmap showing correlations between candidate genes and immune cell populations. Red indicates positive correlations, whereas blue indicates negative correlations. (B). Correlation between FLVCR1 expression and T follicular helper cells. (C-F) Correlations between TMCC2 expression and natural killer cells, activated dendritic cells, immature dendritic cells, and immature B cells, respectively. (G) Correlation between EIF2AK1 expression and CD56dim natural killer cells. (H). Correlation between XK expression and eosinophils. Abbreviations: FLVCR1 = feline leukemia virus subgroup C receptor 1; TMCC2 = transmembrane and coiled-coil domains 2; EIF2AK1 = eukaryotic translation initiation factor 2 alpha kinase 1; CD56dim = cluster of differentiation 56 dim. Please click here to view a larger version of this figure.
Single-cell transcriptome data analysis
To characterize the expression patterns of the identified diagnostic biomarkers within the keloid microenvironment, we analyzed the single-cell RNA sequencing dataset GSE163973. After quality control and data integration, 21,488 high-quality cells were retained for downstream analyses. Cells with fewer than 200 or more than 6,000 total unique molecular identifier (UMI) counts were excluded, and potential doublets were identified and removed using the DoubletDetection package. The 2,000 genes exhibiting the greatest expression variability were selected, followed by dimensionality reduction and visualization with Uniform Manifold Approximation and Projection (UMAP). A total of 10 major cell populations were identified, including endothelial cells, fibroblasts, muscle fibers, keratinocytes, immune cells, lymphatic endothelial cells, glandular cells, neural cells, melanocytes, and an unclassified cell population (Figure 4A,B). Expression profiling revealed distinct cell-type-specific distribution patterns of the diagnostic biomarkers. FLVCR1 was predominantly expressed in endothelial cells and melanocytes, whereas EIF2AK1 showed relatively high expression in neural cells, glandular cells, and fibroblasts. HPX was primarily enriched in melanocytes, while KEL exhibited predominant expression in glandular cells (Figure 4C,D).

Figure 4: Distribution of heme metabolism-related diagnostic biomarkers in the keloid single-cell transcriptome. (A) UMAP plot showing 21 cell clusters comprising 21,488 cells from keloid samples. (B) Cell-type annotations based on the annotations reported in the original study. (C) Feature plots showing the expression of heme metabolism-related diagnostic biomarkers across different cell types. (D) Bubble plot showing the average expression levels and percentages of cells expressing the heme metabolism-related diagnostic biomarkers across different cell types. Abbreviation: UMAP = uniform manifold approximation and projection. Please click here to view a larger version of this figure.
Identification and interaction network analysis of candidate diagnostic biomarkers
To explore the regulatory mechanisms underlying the candidate diagnostic biomarkers, a miRNA–mRNA regulatory network was constructed. To improve the reliability of the predicted interactions, overlapping miRNAs targeting the candidate biomarkers were identified. A total of 282 miRNAs interacting with the six diagnostic biomarkers were obtained, and the resulting regulatory network is shown in Figure 5. Notably, hsa-miR-34a-5p, hsa-let-7a-5p, hsa-let-7d-5p, hsa-let-7e-5p, and hsa-miR-26b-5p were predicted to simultaneously regulate all six candidate biomarkers.

Figure 5: miRNA regulatory network of heme metabolism-related diagnostic biomarkers. The network illustrates the regulatory relationships between the six diagnostic biomarker genes (FLVCR1, HPX, TMCC2, KEL, XK, and EIF2AK1) and their associated miRNAs. Gene nodes represent the diagnostic biomarkers, whereas the surrounding nodes represent miRNAs. Edges indicate experimentally supported miRNA–mRNA interactions. Abbreviations: FLVCR1 = feline leukemia virus subgroup C receptor 1; HPX = hemopexin; TMCC2 = transmembrane and coiled-coil domains 2; KEL = Kell metallo-endopeptidase; XK = X-linked Kx blood group; EIF2AK1 = eukaryotic translation initiation factor 2 alpha kinase 1; miRNA = microRNA; mRNA = messenger RNA. Please click here to view a larger version of this figure.
Experimental validation of FLVCR1 expression and molecular docking analysis of potential therapeutic compounds
To validate the bioinformatic findings and confirm the functional relevance of the identified hub gene, we experimentally assessed FLVCR1 expression in PKF and NHDF. Both qRT-PCR and western blot analyses consistently demonstrated that FLVCR1 was significantly upregulated in keloid fibroblasts compared to normal controls (Figure 6A–C, Supplemental Figure S1, and Supplemental Table S4). This elevated cellular expression supports the potential involvement of FLVCR1-associated heme metabolic dysregulation in keloid pathogenesis.
Given the potential involvement of FLVCR1 in heme metabolism-associated immune changes, we next sought to identify potential therapeutic compounds that could directly target FLVCR1 to interrupt this pathogenic axis. High-throughput virtual screening was performed using a traditional Chinese medicine (TCM) compound library and the prepared protein structure. The 20 compounds with the most favorable docking scores were selected for further evaluation (Supplemental Table S5). In general, a lower binding energy indicates a stronger binding affinity and docking energies below −5 kcal/mol are considered indicative of stable ligand–protein interactions. Among the screened compounds, (+)-Gallocatechin, (−)-Epicatechin, (−)-Gallocatechin, and Cyanidin (chloride) exhibited favorable binding affinities toward FLVCR1. Notably, (+)-Gallocatechin showed the strongest interaction with FLVCR1 by forming four hydrogen bonds with GLU214, ASN245, GLN246, and GLN471, suggesting a stable ligand–protein binding mode (Figure 6D–G). These findings highlight (+)-gallocatechin as a promising candidate for mechanism-based therapeutic intervention targeting FLVCR1.

Figure 6: Experimental validation of FLVCR1 expression and molecular docking of potential compounds targeting FLVCR1. (A) Representative western blot images showing FLVCR1 protein expression in CON and keloid. GAPDH served as the loading control. (B) Quantification of FLVCR1 protein levels normalized to GAPDH. (C) Relative mRNA expression levels of FLVCR1 in CON and keloid fibroblasts were determined by qRT-PCR. GAPDH was used as an internal reference. (D–G) Three-dimensional representations of the predicted binding modes between FLVCR1 and selected small-molecule compounds: (D) (+)-Gallocatechin. (E) (-)-Epicatechin. (F) (-)-Gallocatechin. (G) Cyanidin (Chloride). Abbreviations: FLVCR1 = feline leukemia virus subgroup C receptor 1; CON, control; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; qRT-PCR, quantitative reverse transcription polymerase chain reaction; SD, standard deviation. Data are presented as the mean ± SD. Statistical significance: ns, P > 0.05; *, P < 0.05; **, P < 0.01; ***, P < 0.001; and ****, P < 0.0001. Please click here to view a larger version of this figure.
Confirmation of the stability of the FLVCR1–(+)-Gallocatechin complex with molecular dynamics simulation
To examine the reliability of the predicted ligand–protein binding mode, molecular dynamics (MD) simulation was carried out for the FLVCR1–(+)-gallocatechin complex. The analysis focused on whether the complex remained structurally stable over time and whether ligand binding altered the protein’s conformational behavior using RMSD, RMSF, radius of gyration (Rg), SASA, hydrogen bond analysis, and MM/GBSA calculations. RMSD analysis (Figure 7A) showed that both the apo protein and ligand-bound complex underwent initial fluctuations within the first 20 ns, followed by gradual stabilization, indicating that the systems reached equilibrium during the simulation. After equilibration, the RMSD value of the FLVCR1–(+)-Gallocatechin complex remained below 0.2 nm, suggesting that ligand binding contributed to maintaining the structural stability of FLVCR1. RMSF analysis (Figure 7B) demonstrated that most residues exhibited limited fluctuations throughout the simulation, indicating the preservation of overall protein integrity, while several flexible regions may represent loop regions involved in ligand accommodation. Furthermore, stable Rg and SASA profiles (Figure 7C,D) indicated that the complex maintained a compact conformation, with no obvious structural expansion or changes in solvent exposure. Hydrogen bond analysis (Figure 7E) revealed that the FLVCR1–(+)-Gallocatechin complex maintained persistent intermolecular interactions, with approximately 3–4 hydrogen bonds formed during the simulation, supporting the stability of the ligand–protein association. MM/GBSA analysis further showed that the FLVCR1–(+)-Gallocatechin complex exhibited a favorable binding free energy (ΔGtotal = −34.87 ± 4.13 kcal/mol) (Supplemental Table S6). Energy decomposition analysis indicated that van der Waals interactions (ΔVDWAALS = −46.34 ± 2.16 kcal/mol) and electrostatic interactions (ΔEelec = −14.09 ± 3.45 kcal/mol) were the major favorable contributors to binding, despite the unfavorable contribution from polar solvation energy (ΔGsolvation = 25.55 ± 0.74 kcal/mol) (Supplemental Table S6). Collectively, these MD simulation results demonstrated that (+)-Gallocatechin formed a stable complex with FLVCR1 and further supported the reliability of the docking-predicted binding mode.

Figure 7: Molecular dynamics simulation analysis of the FLVCR1–(+)-Gallocatechin complex. (A) RMSD profiles of apo FLVCR1 and the FLVCR1–(+)-Gallocatechin complex during the 100 ns molecular dynamics simulation. (B) RMSF profile showing residue-level fluctuations of FLVCR1 during the simulation. (C) SASA profile showing changes in the solvent-accessible surface area of the FLVCR1–(+)-Gallocatechin complex. (D) Rg profile evaluating the compactness of the FLVCR1–(+)-Gallocatechin complex during the simulation. (E) Hydrogen-bond analysis showing dynamic intermolecular interactions between FLVCR1 and (+)-Gallocatechin throughout the simulation. Abbreviations: FLVCR1 = feline leukemia virus subgroup C receptor 1; RMSD = root mean square deviation; RMSF = root mean square fluctuation; SASA = solvent-accessible surface area; Rg = radius of gyration. Please click here to view a larger version of this figure.
Data Availability:
The publicly available transcriptomic datasets analyzed in this study can be accessed through the Gene Expression Omnibus (GEO) under accession numbers GSE44270, GSE7890, and GSE163973. The source data generated in this study and underlying the experimental validation, including qRT-PCR measurements, original western blot images, and western blot quantification data, are provided as Supplementary Figure S1, Supplementary Table S1, Supplementary Table S2, Supplementary Table S3, and Supplementary Table S4. The molecular docking results and MM/GBSA binding free energy data are also provided in Supplementary Table S5 and Supplementary Table S6.
Supplemental Figure S1: Original data of western blotting.Please click here to download this file.
Supplemental Table S1: Heme metabolism-related genes.Please click here to download this file.
Supplemental Table S2: Primer sequences of selected genes. Please click here to download this file.
Supplemental Table S3: Machine-learning approaches for identifying potential diagnostic biomarkers in keloid. Please click here to download this file.
Supplemental Table S4: Raw data supporting the experimental validation of FLVCR1 expression. Please click here to download this file.
Supplemental Table S5: Top 20 candidate compounds identified by molecular docking with FLVCR1.Please click here to download this file.
Supplemental Table S6: MM/GBSA binding free-energy analysis of the FLVCR1–(+)-Gallocatechin complex.Please click here to download this file.