Integrative Bioinformatics Analysis and Machine Learning Reveal TSPAN7 Regulates the Involution of Infantile Hemangioma via NK Cells

Introduction

Hemangiomas are developmental vascular anomalies that may present as hamartomas, benign tumors, or vascular malformations, and are classified as infantile or congenital types. Infantile hemangiomas (IH) are GLUT1-positive lesions, unlike congenital types, that typically appear within the first two months of life, grow rapidly during the first year, and gradually regress, whereas congenital hemangiomas are present at birth and may become symptomatic due to growth, infection, hormonal changes, or trauma.^1–3 While most IH lesions undergo spontaneous regression without functional sequelae, some may cause ulceration, disfigurement, or organ dysfunction, requiring medical intervention. Propranolol remains the first-line therapy, but some patients show an incomplete response or relapse after discontinuation, highlighting the need for a deeper understanding of the molecular mechanisms underlying IH involution.^4–6

The involution phase is histologically defined by the gradual replacement of proliferation endothelial cells with fibrofatty tissue, a process driven by the differentiation of hemangioma stem cells (HemSCs) into adipocytes.^7–9 Emerging evidence suggests that lipid metabolism plays a pivotal role in this transition. Transcriptomic analyses have revealed significant alterations in lipid-associated gene expression between proliferation and involution IH tissues.^10,11 The activation of key adipogenic regulators such as peroxisome proliferator-activated receptor gamma (PPARγ) has been shown to promote hemangioma-derived mesenchymal stem cell (Hem-MSC) differentiation toward an adipocytic lineage while suppressing angiogenesis, thereby facilitating IH involution.^12,13 Moreover, metabolic pathway analyses have indicated substantial metabolic reprogramming during IH involution, but the integration of lipid-related changes with the immune microenvironment remains poorly characterized.

Recent bioinformatics studies have used the GSE127487 dataset to explore gene expression dynamics across IH phases. Weighted gene co-expression network analysis (WGCNA) and machine learning approaches have identified candidate genes such as SH3KBP1 associated with disease progression, and immune infiltration analyses have suggested dynamic shifts in immune cell populations during proliferative and involutive stages.¹⁴ However, systematic screening for lipid metabolism-related genes that drive involution and elucidation of their interplay with immune cells have not been performed.

In this study, we integrated time-series expression analysis, differential expression profiling, and machine learning algorithms to identify key lipid metabolism-related genes associated with IH involution using the GSE127487 dataset. Through WGCNA co-expression analysis, immune infiltration assessment, and mediation analysis, we explored the potential relationships between hub genes and immune cells. These findings may provide clues for further investigations into the roles of lipid-related genes and immune cells in IH involution.

Methods

Data Processing and Differential Expression Analysis

Figure 1 illustrates the technical workflow of this study. The GSE127487 microarray dataset was downloaded from the Gene Expression Omnibus (GEO) database using the R package GEOquery. The dataset was generated on a GPL10558 platform (Illumina HumanHT-12 V4.0 expression beadchip). A total of 23 IH tissue samples and five normal skin tissue samples from Homo sapiens were included. Among these, five IH samples were obtained from patients who received propranolol treatment.¹⁵ For subsequent analyses, we selected six untreated proliferation-phase IH (PIH) samples and six untreated involution-phase IH (IIH) samples to avoid confounding factors caused by treatment.

Figure 1 Diagram of the study work flow.

Time-series gene expression patterns across 6, 12, and 24 months were analyzed using the Mfuzz package in R, which employs fuzzy c-means clustering to identify clusters with distinct temporal trajectories.

Differentially expressed genes (DEGs) between PIH and IIH samples were identified using the limma package with thresholds of |log₂ fold change (logFC)| > 2 and P-value < 0.05. Key involution-related genes were identified from the overlap between involution-associated temporal genes and DEGs.

Lipid metabolism-related genes were retrieved from the GeneCards database using the keyword “lipid” and genes with a relevance scores > 2 were retained. Overlapping these lipid-associated genes with key involution-related genes generated lipid-related involution genes.

Machine Learning-Based Hub Gene Screening

Three machine learning algorithms were applied to identify hub genes from lipid-related involution genes.

Least absolute shrinkage and selection operator (LASSO) logistic regression was performed using the glmnet R package with 10-fold cross-validation. The optimal regularization parameter λ was selected as the value within one standard error of the minimum cross-validated deviance (λ.1se), balancing model sparsity and predictive performance. Genes with non-zero coefficients at the selected λ were retained as candidates.

Support vector machine with recursive feature elimination (SVM-RFE) was implemented using the e1071 and caret R packages with a linear kernel. Features were recursively eliminated based on the absolute values of the SVM weight vector, and the optimal feature subset was determined as the set that yielded the highest cross-validation accuracy.

Random forest (RF) was performed using the randomForest R package with 500 trees. Feature importance was computed by the mean decrease in Gini index, and the top 10 genes with the highest importance scores were retained.

To prevent overfitting and ensure model stability, 10-fold cross-validation repeated three times was employed for model tuning. Both LASSO and RF models used this cross-validation scheme. For SVM-RFE, feature selection was nested within cross-validation loops to avoid information leakage.

The genes identified by the three algorithms were considered hub genes. This intersection strategy reduces algorithm-specific bias and enhances the robustness of feature selection given the limited sample size.

Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA was performed using the WGCNA package to construct co-expression networks specific to involution phase. Soft-thresholding power was selected based on scale-free topology fit. Modules with high similarity were merged and module eigengenes were calculated. Pearson’s correlation was used to identify modules significantly associated with TSPAN7 expression in involution phase. Modules with a correlation |r| > 0.5 and P < 0.05, were selected for further GO and KEGG enrichment analysis.

Functional Enrichment Analysis

To explore the biological functions and signaling pathways of the identified modules, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the clusterProfiler package in R. Enriched terms with P < 0.05 were considered statistically significant.

Immune Infiltration Analysis

The MCP-counter algorithm was applied to quantify the absolute abundance of ten immune cell populations and stromal components in each sample using the MCP-counter package. Differences in immune cell infiltration between the PIH and IIH groups were assessed using the Wilcoxon rank-sum test.

Mediation Analysis

A mediation analysis was performed to evaluate whether the effect of TSPAN7 on IH involution was mediated by immune cells. Linear regression models were constructed using the mediation package in R with bootstrapping (1000 simulations) to estimate the direct, indirect, and total effects. The significance of the mediation effect was determined by using a 95% confidence interval.

Statistical Analysis

All statistical analyses were performed using the R software (version 4.2.0). For comparisons between two groups, normally distributed variables were analyzed using the Student’s t-test, while non-normally distributed variables were compared using the Mann–Whitney U-test (Wilcoxon rank-sum test). Multiple testing corrections were applied using the Benjamini-Hochberg false discovery rate (FDR) method, where appropriate. A two-sided P < 0.05 was considered statistically significant.

Results

Identification of Involution-Related Lipid-Associated Hub Genes

Time-series expression analysis using Mfuzz clustering across 6, 12, and 24 months identified nine distinct expression clusters (Figure 2). Based on their temporal trajectories, clusters exhibiting gradually increasing expression (cluster 1), gradually decreasing expression (cluster 4), or stable expression during 6 and 12 months followed by specific downregulation during 24months (cluster 6) were selected as involution-associated temporal genes. Differential expression analysis of proliferation-phase IH (PIH) and involution-phase IH (IIH) samples identified 399 DEGs, including 241 upregulated and 158 downregulated genes (|logFC| > 2, P < 0.05). Intersection analysis of involution-associated temporal genes and DEGs identified 281 core genes linked to involution.

Figure 2 Time-series expression patterns of genes across infantile hemangioma phases (6, 12, and 24 months) revealed by Mfuzz clustering. The x-axis represents the time points (proliferation phase at 6 months, stable phase at 12 months, and involution phase at 24 months), and the y-axis represents relative gene expression levels. Nine distinct clusters (Cluster 1–Cluster 9) were identified, each exhibiting a unique temporal trajectory. Colors indicate the membership degree of each gene in the assigned cluster, with red representing high membership and yellow representing low membership.

A total of 3773 lipid-related genes with relevance scores > 2 were retrieved from the GeneCards database. Overlapping these genes with 281 involution-related key genes identified 23 lipid-related involution genes (Figure 3A). Three machine learning algorithms, LASSO, SVM-RFE, and RF were applied to further refine these candidates. LASSO regression identified 12 genes with non-zero coefficients (Figure 3B and C). SVM-RFE selected seven genes with optimal feature importance (Figure 3D), and random forest analysis ranked the genes based on the mean decrease in accuracy and the Gini coefficient (Figure 3E). By intersecting candidate genes screened via the three algorithms, we identified three hub genes: CFH, TSPAN7, and MIR324 (Figure 3F). Among them, TSPAN7 and CFH displayed a gradually rising expression trend from the proliferation stage to the involution stage (cluster 4). In contrast, MIR324 maintained stable expression at 6 and 12 months, followed by distinct downregulation at 24 months.

Figure 3 Identification of involution-related lipid-associated hub genes. (A) Venn diagram showing the intersection of involution-related key genes and lipid-associated genes, yielding 23 lipid-related involution genes. (B and C) LASSO regression analysis identified 12 genes with non-zero coefficients (B: Coefficient profile; (C) Binomial deviance). (D) SVM-RFE selected 7 genes with optimal feature importance. (E) Random forest analysis ranking genes based on Mean Decrease Gini (>0.2). (F) Venn diagram illustrating the overlap of genes selected by the three machine learning algorithms (LASSO, SVM-RFE, and random forest), identifying three hub genes (CFH, TSPAN7, MIR324).

Weighted Gene Co-Expression Network Analysis Identifies TSPAN7-Associated Modules

To explore functional networks co-expressed with TSPAN7, WGCNA was performed. A soft-thresholding power of 8 was selected to achieve a scale-free topology fit (R² = 0.85) (Figure 4A). Hierarchical clustering identified 20 co-expression modules (Figure 4B). Among these, four modules (blue, brown, turquoise, and yellow) were significantly correlated with TSPAN7 expression (correlation |r| > 0.5, P < 0.05) (Figure 5). These modules were selected for the functional enrichment analysis.

Figure 4 Weighted gene co-expression network analysis (WGCNA). (A) Analysis of scale-free topology fit. A soft-thresholding power of 8 was selected to achieve a scale-free topology fit (R2 = 0.85). (B) Hierarchical clustering dendrogram showing the identification of 20 co-expression modules, each labeled with a distinct color.

Figure 5 Correlation between TSPAN7 and module eigengenes (MEs) in PIH and IIH. Red indicates positive correlation, blue indicates negative correlation. *P < 0.05.

GO and KEGG enrichment analyses of the four TSPAN7-associated modules revealed a striking convergence in immune-related functions. The most significantly enriched biological processes included “natural killer cell-mediated cytotoxicity”, “leukocyte-mediated cell killing”, “regulation of natural killer cell-mediated immunity”, and “adaptive immune response”. These findings suggest a potential association between TSPAN7 and NK cells during IH involution.

Immune Infiltration Landscape During IH Involution

To characterize immune cell dynamics, MCP-counter analysis was performed to estimate immune cell abundance in PIH and IIH samples. The results demonstrated a significant enrichment of NK cells in the proliferation phase compared to that in the involution phase. In contrast, myeloid dendritic cells were significantly more abundant during the involution phase (Figure 6).

Figure 6 Immune cell infiltration landscape in proliferation and involution phases of IH assessed by MCP-counter. The x-axis represents immune cell populations, and the y-axis represents relative gene expression value. *P < 0.05; **P < 0.01.

TSPAN7 Promotes Involution Through NK Cells

Given the functional enrichment of NK cells in TSPAN7-associated co-expression modules and the observed enrichment of NK cells during proliferation, we hypothesized that NK cells might mediate the effect of TSPAN7 on involution. Mediation analysis was performed to test this hypothesis, with TSPAN7 expression as the independent variable, NK cells abundance as the mediator, and involution status as the outcome variable.

The results revealed the existence of a complete mediation model. The direct effect of TSPAN7 on involution was not significant (Effect= 0.1113, P = 0.4740), whereas the indirect effect mediated by NK cells was highly significant (Effect= 0.1293, P = 0.0120). The total effect of TSPAN7 on involution was significant (Effect= 0.2406, P = 0) (Table 1). Path decomposition showed that TSPAN7 was negatively associated with NK cells abundance (path a = −0.2095, P < 0.01), and that, after controlling for TSPAN7, NK cells abundance was also negatively associated with involution (path b = −0.6174). The product of these two negative paths (a × b) yields a positive indirect effect.

Table 1 Mediation Analysis Results of TSPAN7 and NK Cells

Discussion

IH is distinguished by its unique capacity for spontaneous involution, but the precise molecular mechanisms orchestrating this regression remain unclear.^16,17 In a previous bioinformatics study by our group, we identified IL6, EGFR, IRF1 and IL32 as key pyroptosis-related genes involved in the regression of IH using the same GSE127487 dataset.¹⁸ The present study extends that investigation by focusing on lipid-associated genes. Although lipid metabolism reprogramming has been implicated in the transition from the proliferative to involution phases,^13,19,20 the interplay between lipid-associated genes and the immune microenvironment during IH involution has not been systematically characterized. By integrating time-series clustering and machine learning algorithms (LASSO, SVM-RFE, and random forest), we newly identified CFH, MIR324 and TSPAN7 as key lipid metabolism-related genes that may associated with IH involution. Importantly, WGCNA and mediation analysis revealed that the effect of TSPAN7 may be associated with NK cells during IH involution, a possibility not explored in our prior work.

Three hub lipid-related involution genes were identified: CFH, TSPAN7, and MIR324. CFH belongs to cluster 4, exhibiting a progressively increasing expression pattern during IH involution. CFH (Complement Factor H) is a key regulator of the alternative complement pathway and has been shown to participate in lipid metabolism. In human adipose tissue, high expression of the CFH gene is present in the stromal vascular fraction cells (comprising monocytes, fibroblasts, and endothelial cells), where it mediates the activation of the alternative complement pathway and participates in lipid metabolism reprogramming and insulin resistance.²¹ CFH can also bind to ox-LDL, directly participating in the clearance of oxidized lipids and affecting lipid deposition.²² The stromal vascular fraction cells of the involution phase of IH may highly express CFH, potentially participating in local lipid clearance, transport, or deposition, thereby directly or indirectly influencing adipose tissue formation following vascular regression. The MIR324 gene simultaneously encodes miR-324-5p and miR-324-3p. These two mature microRNA molecules originate from the same precursor hairpin structure but exhibit significant functional differences and have crucial impacts on the pathogenesis of a wide range of disorders.²³ miR-324-3p is downregulated in proliferation-phase IH tissues and functions as a tumor suppressor by targeting the PDRG1 pathway.²⁴ MIR324 plays different roles in different tumor microenvironments, and its role in IH involution warrants further investigation. TSPAN7 exhibited a progressively increasing expression pattern across IH developmental phases, suggesting its potential role in promoting involution. TSPAN7, a member of the tetraspanin family, regulate cell adhesion, migration, and immune synapse formation.²⁵ Previous studies have demonstrated that tetraspanins can modulate immune cell function by organizing membrane microdomains that facilitate receptor clustering and signaling.^26,27 However, its specific role in IH involution has not yet been reported.

The functional enrichment analysis of TSPAN7-associated co-expression modules revealed a striking convergence on immune-related pathways, particularly “natural killer cell-mediated cytotoxicity” and “leukocyte-mediated cell killing”, suggesting that TSPAN7 may participate in IH involution together with NK cells. The immune microenvironment of IH has been increasingly recognized as a critical regulator of disease progression and involution. Macrophages undergo a phenotypic switch from M2-predominant in the proliferative phase to M1-predominant in the involution phase.²⁸ Single-cell RNA sequencing has further resolved the cellular heterogeneity of IH, confirming that macrophages, mast cells, and telocytes dynamically change during involution.¹⁰ Several immune infiltration analyses have revealed alterations in NK cells in IH,^14,29 and accumulating evidence from tumor angiogenesis and vascular biology studies supports the possibility of NK cells-mediated vascular regression.^30–32 However, the role of NK cells in IH has remained largely unexplored. Our immune infiltration analysis revealed that NK cells were significantly enriched in proliferation phase compared to involution phase IH tissues, indicating a phase-specific distribution of NK cells during IH progression. The precise role of NK cells in IH involution remains to be further explored.

Pan-cancer analysis revealed that TSPAN7 expression was positively associated with the infiltration of resting NK cells in renal clear cell carcinoma, supporting a potential link between TSPAN7 and NK cells.³³ This finding prompted us to explore the potential relationship between NK cells and TSPAN7 during IH involution. The complete mediation model demonstrated that the effect of TSPAN7 on IH involution may depend on NK cells, as the direct effect was not significant, whereas the indirect effect mediated by NK cells was highly significant. This finding suggests that TSPAN7 may be associated with IH involution through NK cells. Although direct evidence for TSPAN7-mediated NK cells regulation is lacking, tetraspanins are known to induce cytoskeleton reorganization, which may modulates NK cells activation.³⁴ It should be noted that the finding that TSPAN7 expression correlates negatively with NK cells abundance yet exerts its pro-involution effect through NK cells (complete mediation) may initially appear counterintuitive. We propose that this pattern reflects a quantitative-qualitative trade-off in NK cells biology during IH involution. In the context of IH, TSPAN7 upregulation on hemangioma endothelial cells or infiltrating immune cells may enhance the cytotoxic efficiency of each NK cells, leading to more effective clearance of proliferation endothelium. Following target elimination, activated NK cells may undergo activation-induced cell death (AICD).^35,36 Consequently, the decline in NK cells number at the involution stage may be not a cause of reduced immune surveillance but rather a reflection of successful task completion. The complete mediation model thus indicates that the pro-involution effect of TSPAN7 is conditional upon the presence and activation of NK cells, even if their numerical abundance has started to wane by the time TSPAN7 peaks. However, we acknowledge that causal mediation analysis relies on strong assumptions that are unlikely to be fully satisfied in observational transcriptomic data, including no unmeasured confounding between the exposure and the mediator, no unmeasured confounding between the mediator and the outcome, no exposure-mediator interaction, and correct temporal ordering of variables. Given the cross-sectional nature of the GSE127487 dataset and the limited sample size, these assumptions cannot be rigorously verified. Therefore, the results of the mediation analysis should be interpreted exploratorily, as hypothesis-generating evidence rather than causal inference. We hypothesize that TSPAN7 may modulate NK cells function or turnover, but this requires experimental validation.

This study has several limitations. First, our findings were derived entirely from bioinformatics analyses of public microarray data and require validation through experimental approaches, including in vitro and in vivo models. Second, the precise molecular mechanisms by which TSPAN7 regulates NK cells function in IH remain to be determined. Third, although we excluded propranolol-treated samples to avoid confounding effects, the limited sample size may affect the robustness of our analyses, particularly for WGCNA and machine learning, which typically require larger cohorts. Therefore, our findings should be considered exploratory and hypothesis-generating, requiring validation in independent and larger datasets. Fourth, the roles of other lipid-related hub genes identified in this study, including CFH and MIR324, warrant further investigation.

In conclusion, this study identified three key lipid metabolism-related genes associated with IH involution (CFH, TSPAN7, and MIR324). Among them, TSPAN7 may be associated with IH involution through NK cells. The precise mechanisms of these molecules during IH involution require further experimental investigation. These findings provide clues for further investigations into the roles of lipid-related genes and immune cells in IH involution.