Integrative Bioinformatics Analysis of Pyroptosis-Related Genes and Analysis of Immune Cell Infiltration in Infantile Hemangioma Regression

Introduction

Infantile hemangioma (IH) is the most common benign tumour in infants.¹ IHs have unique clinical characteristics and can be divided into three phases: proliferation, quiescence and involution.² IHs usually arise within the first few weeks of life and proliferate for 4–18 months of age. The involution phase starts at approximately 12 months of age, when the hemangioma tissue is replaced by adipose and fibrous tissue.^1,3 Exploring the mechanism underlying IH progression from proliferation to involution is very important for the development of new therapeutic targets.

The mechanism underlying spontaneous IH regression is unclear. In the involution phase, increased apoptosis and decreased Bcl-2 expression are observed, indicating that programmed cell death (PCD) is involved in this process.^4,5 Pyroptosis is a type of PCD that is distinct from apoptosis, and it is characterized by cell swelling, cytolysis and the release of inflammatory factors. The two pathways of pyroptosis include the caspase-1-dependent classical pathway and the caspase-4/5/11-mediated nonclassical pathway.⁶ These pyroptosis pathways activate gasdermin D (GSDM) proteins to translocate to the cell membrane and form pores. Subsequently, the cell continuously swells until the membrane ruptures, which is accompanied by the release of inflammatory cytokines such as IL-1β and IL-18 into the extracellular space, which recruit more inflammatory cells to aggregate, inducing a severe inflammatory response.⁷ Pyroptosis is widely involved in tumours,⁸ metabolic diseases,⁹ and inflammatory diseases.¹⁰ However, few studies have investigated the role of pyroptosis in hemangiomas, and the mechanism that links pyroptosis and IH regression still needs to be elucidated.

Bioinformatics analysis is a methodology that includes processing, analysing, and interpreting biological data by integrating methods and technologies from multiple fields, such as computer science, statistics, and biology. Previous studies have used bioinformatics to identify potential genes involved in the pathogenesis of infantile hemangiomas.^11–13 In this study, we used various bioinformatic tools to investigate the relationships between pyroptosis-related genes (PRGs) and IHs. We identified 14 PRGs related to IH regression (IH-PRGs) from the PubMed database and performed a functional enrichment analysis of these genes to further understand their role in IH regression. The hub genes from the IH-PRGs were subsequently screened with PPI network. The diagnostic and predictive values of the hub-IH-PRGs were evaluated, and their correlations with IH immune cell infiltration were analysed. This study on the role of pyroptosis in IH regression has increased our understanding of IH pathogenesis and provided new insights for its treatment.

Methods

Data Processing and Differential Expression Gene Analysis

Figure 1 shows the technical roadmap used for this study. The R package GEOquery was used to download GSE127487 array data from the GEO database. The chip was sequenced on the GPL10558 (Illumina HumanHT-12 V4.0 expression beam) platform. The samples were obtained from 23 samples of hemangioma tissue from Homo sapiens and 5 samples of normal skin tissue.¹⁴ Among these samples, 5 samples of hemangioma tissue were from patients who were treated with propranolol. We selected 6 samples of untreated proliferation-phase hemangioma (PIH) tissue and 6 samples of untreated involution-phase hemangioma (IIH) tissue for this study.

Figure 1 Diagram of the study work flow.

To identify the differential expressed genes between PIHs and IIHs, the limma package of R software was used. A box plot was constructed to check the normalization of the samples. We defined differentially expressed genes (DEGs) as those with a |logFC|>1 and a P value<0.05. Through a PubMed search for research related to cell apoptosis, a total of 356 genes (pyroptosis related genes, PRGs) related to cell apoptosis were retrieved from the GeneCards database, AmiGO2 database, and MSigDB database.¹⁵ The differentially expressed genes associated with pyroptosis in IH regression (IH-PRGS) were obtained by intersecting the DEGs with the PRGs.

Functional Enrichment Analysis

To explore the functional and pathway enrichment of IH-PRGs, Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed with the clusterProfiler package in R. GO terms and pathways with P values<0.05 were considered significant.

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) was employed to evaluate the distribution pattern of genes within a predetermined set of a gene table This is accomplished by sorting genes on the basis of their correlation with genotypes to determine their impact on phenotypes. We obtained the gene sets “c5.go.all.v2022.1.Hs.symbols.gmt” and “c2.cp.all.v2022.1.Hs.symbols.gmt” from MsigDB and used the clusterProfiler package in R to perform GSEA, where a p value less than 0.05 was considered statistically significant.

Protein‒protein Interaction Network Construction and Hub Genes Identification

We used the STRING database to construct a protein‒protein interaction network with the IH-PRGs. Cytoscape 3.6.1 was used to visualize the results. We adopted CytoHubba, a Cytoscape plug-in, to screen the hub genes (hub-IH-PRGs).

Construction of Interaction Networks Between IH-PRGs and Related miRNAs and Transcription Factors

The TarBase v8.0 and ENCODE databases were used to identify miRNAs and TFs associated with IH-PRGs. The target IH-PRG-miRNA and IH-PRG-TF networks were subsequently visualized with Cytoscape software.

Expression Analysis and ROC Validation of the Hub-IH-PRGs

We analysed the expression levels of the hub-IH-PRGs in the PIH and IIH groups, and the results are presented as a box plot. The ROC curve is used to evaluate the diagnostic and predictive value of these genes, where an area under the ROC curve (AUC) value greater than 0.7 is considered accurate for prediction purposes.

Immune Cell Infiltration Analysis and Correlation Analysis of the Hub-IH-PRGs

We used single-sample gene set enrichment analysis (ssGSEA) to analyse the proportion of immune cells in each sample to identify significantly differentially enriched immune cells between the PIH and IIH groups. Finally, we calculated the Pearson correlation coefficient between the expression levels of the hub-IH-DEGs and the significantly enriched immune cells.

Statistical Analysis

All the data analyses were conducted with R software version 1.2. For the comparison of two consecutive variables, the statistical significance of normally distributed variables was estimated with an independent Student’s t test, and the differences between nonnormally distributed variables were analysed with the Mann‒Whitney U-test (ie, the Wilcoxon rank sum test). Bonferroni adjustment was used for multiple testing correction.

Results

Differential Expression Gene Analysis

Normalized box plots (Figure 2) revealed that the signal intensity of each sample in the two groups was approximately at the median level, indicating a good degree of sample normalization. On the basis of the screening criteria, 1554 DEGs were identified, of which 731 were upregulated and 823 were downregulated. A volcano plot and heatmaps of the data are displayed in Figure 3A and B. The intersection of PRGs with DEGs was selected after Venn analysis (Figure 4), and 14 overlapping genes were identified as IH-PRGs.

Figure 2 Normalized box plots of the GSE127487 dataset samples. Blue represents the IIH group, and red represents the PIH group. The abscissa in the figure represents the sample number, and the ordinate represents the chip signal intensity.

Figure 3 Differentially expressed genes (DEGs). (A) Volcano plot of the DEGs between IIHs and PIHs. The abscissa is the log2-fold change, the ordinate is -log10 (p value), red nodes represent upregulated DEGs, blue nodes represent downregulated DEGs, and grey nodes represent genes that are not significantly differentially expressed. (B) Heatmaps of the DEGs. The horizontal axis indicates the sample number, the vertical axis indicates the respective DEGs, red represents high gene expression, blue represents low gene expression, red bars indicate PIH tissue, and blue bars indicate IIH tissue.

Figure 4 Venn diagrams of the differentially expressed genes (DEGs) and pyroptosis-related genes (PRGs). The pink circles represent the DEGs between PIHs and IIHs, and the blue circles represent the genes related to pyroptosis. Fourteen genes that were common to the IH regression-associated DEGs and pyroptosis-related genes were identified: IL6, MIR20B, DPEP1, IL13RA2, PCSK9, GBP5, MALAT1, EGFR, CLEC5A, MELK, IL32, IRF1, ADORA3, and TRIM24.

Functional Enrichment Analysis of IH-PRGs

Subsequently, GO functional enrichment analysis was performed on the IH-PRGs (Figure 5A). The differentially expressed IH-PRGs were associated with biological processes, such as the cellular responses to chemical stimuli, response to stress, response to chemical, cellular response to organic substance, apoptotic process, cell promotion, programmed cell death, regulation of apoptotic process, cellular response to oxidative stress, regulation of programmed cell death, etc. (Figure 5B). Additionally, the differentially expressed IH-PRGs are enriched in cellular components, such as the interleukin-6 receptor complex, multiple bodies, internal vessel, extracellular space, pericrystalline fibres, plasma membrane part, extracellular region part, extrinsic component of the external side of the plasma membrane, endolysosome membrane, cell surface, endolysosome, etc. (Figure 5C). In terms of molecular function, these genes are enriched in virus receptor activity, hijacked molecular function, very low-density lipoprotein particle binding, metal dipeptidase activity, signalling receptor binding, oestrogen response element binding, very low-density lipoprotein particle receptor binding, protein tyrosine kinase activity, molecular function regulator, and receptor regulator activity (Figure 5D). The KEGG annotation of the IH-PRGs revealed that they were enriched in the JAK-STAT signalling pathway, pertussis, cytokine–cytokine receptor interaction, EGFR tyrosine kinase inhibitor resistance, C-type lectin receptor signalling pathway, HIF-1 signalling pathway, TNF signalling pathway, and FoxO signalling-related pathways (Figure 5E and F).

Figure 5 GO and KEGG enrichment analyses. (A) GO enrichment analysis. The ordinate indicates the -log(p value), the abscissa represents the enriched GO terms, and the bar colours indicate the BP (red), CC (green) and MF (blue) GO terms. (B–D) Chord plot showing the top 10 enriched BP (B), CC (C) and MF (D) GO terms. (E) KEGG pathway enrichment analysis. The abscissa indicates the size of the -log10(p value), and the ordinate indicates the pathway name. (F) Chord plot showing the top 10 enriched KEGG pathways.

Gene Set Enrichment Analysis

The GSEA results revealed that the genes associated with PIHs and IIHs were primarily involved in nuclear chromosome segregation, mitotic nuclear division, meiotic cell cycle process, mitotic cell cycle phase transition, epidermis development, skin development, regulation of antigen processing and presentation, and antigen processing and presentation (Figure 6A). The biological pathways impacted by gene expression in the dataset included the cornified envelope, keratinization, type I hemidesmosome assembly, apoptotic cleavage of cell adhesion proteins, cell cycle checkpoints, cell cycle mitosis, mitotic metaphase and anaphase, mitotic prometaphase, Rho GTPases activate formins, and ATM signalling in development and disease (Figure 6B–F).

Figure 6 GSEA–GO and KEGG analyses of the GSE127487 dataset. (A) GSEA–GO analysis of the GSE127487 dataset. In the distribution curves, the y-axis represents the gene set, the x-axis represents the log fold change (FC) distribution of the core molecules in each gene set, and the shape of the peak represents the logFC of the core molecules in the gene set. The curve peak height corresponds to the position where the logFC of most molecules in this group is concentrated. If the normalized enrichment score (NES) of the corresponding gene set is negative, the peak of the gene set is generally to the left of zero; if the NES of the corresponding gene set is positive, the peak of the gene set is generally to the right of zero. (B–F) GSEA–KEGG analysis of GSE127487 dataset showing that the enriched pathways were formation of the cornified envelope, keratinization, type I hemidesmosome assembly, apoptotic cleavage of cell adhesion proteins, cell cycle checkpoints, cell cycle mitosis, mitotic metaphase, and anaphase, mitotic prometaphase, Rho GTPases activate formins, and ATM signalling during development and disease.

PPI Network Construction and Hub Gene Identification

There were a total of 14 IH-PRGs and 7 PPI pairs in the constructed PPI networks. The four genes with the strongest cooperative relationships were IL6, EGFR, IRF1, and IL32 (Figure 7).

Figure 7 Protein–protein interaction (PPI) network. (A) PPI network of IH-associated pyroptosis-related genes (IH-PRGs) constructed with the STRING database. Each network node represents a protein, and the lines represent protein–protein associations. (B) Visualization produced in R software of the PPI network depicting the protein interactions of IH-PRGs from the table provided by the STRING database. The 4 genes with the highest number of interactions are shown.

Analysis of IH-PRGs and Related miRNAs and Transcription Factors

We constructed an IH-PRG-miRNA interaction network consisting of 12 genes and 188 miRNAs (Figure 8A). The top five differentially expressed IH-PRGs were IRF1 (targeted by 73 miRNAs), EGFR (targeted by 56 miRNAs), IL6 (targeted by 40 miRNAs), MELK (targeted by 32 miRNAs), and TRIM (targeted by 24 miRNAs). The IH-PRG-TF interaction network consisted of 11 genes and 245 TFs (Figure 8A). The top five IH-PRGs were MALAT1 (regulated by 197 TFs), IRF1 (regulated by 66 TFs), TRIM24 (regulated by 26 TFs), and IL32 (regulated by 24TFs) (Figure 8B).

Figure 8 Correlations between IH-associated pyroptosis-related genes (IH-PRGs) and miRNAs and transcription factors (TFs). (A) IH-PRG-miRNA network. The blue nodes represent miRNAs, and the red nodes represent IH-PRGs. (B) IH-PRG-TF network. The green nodes represent TFs, and the red nodes represent IH-PRGs.

Hub-IH-PRG Expression Analysis and ROC Validation

The expression levels of the hub-IH-PRGs in the GSE127487 dataset are presented as a box plot (Figure 9A). In the GSE127487 dataset, there was a statistically significant difference in the expression levels of IL6, EGFR, IRF1, and IL32 (p<0.05). The expression levels of IL6 and EGFR were lower in PIHs than in IIHs, whereas the expression levels of IRF1 and IL32 were greater in PIHs than in IIHs. The ROC curve revealed that all four genes have diagnostic value for PIH (AUC>0.7) (Figure 9B).

Figure 9 Expression of hub-IH-associated pyroptosis-related genes (hub-IH-PRGs). (A) Expression levels of the 4 hub-IH-PRGs in the GSE127487 dataset. Red represents the PIH group, blue represents the IIH group, the abscissa represents the hub-IH-PRGs, and the ordinate represents the gene expression value. (B) Receiver operating characteristic (ROC) curve prediction of the hub-IH-PRGs in the GSE127487 dataset. *p<0.05, ** p<0.01 compared with control.

Immune Cell Infiltration Analysis

Tdg cells were significantly more abundant in PIH tissue than in IIH tissue, whereas Th17 and CD8 T cells were significantly less abundant in PIH tissue than in IIH tissue (Figure 10A). Further analysis of the correlation between the hub-IH-PRGs and significantly different infiltrating immune cells revealed a significant positive correlation between IL-6 and Th17 cell infiltration and a significant negative correlation with Tgd cell infiltration; EGFR was negatively correlated with Tgd cell infiltration; and IRF1 and IL32 were significantly negatively correlated with Th17 cell infiltration (Figure 10B).

Figure 10 Analysis of immune cell infiltration. (A) Differences in the enrichment of 24 immune cells between the PIH and IIH groups. Red represents the PIH group, blue represents the IIH group, the abscissa represents the immune cells, and the ordinate represents the score. (B) Heatmap of the correlation analysis between the hub-IH-PRGs and infiltrating immune cells. Blue indicates a negative correlation with immune cells, red indicates a positive correlation with immune cells, the horizontal axis indicates hub-IH-PRGs, and the vertical axis indicates immune cells. *p<0.05, ** p<0.01 compared with control.

Discussion

IH is characterized by spontaneous regression, but the mechanism underlying IH regression is unclear.^1,16 PCD is involved in the spontaneous regression of infantile hemangiomas.^4,5 As a new type of programmed cell death, the relationship between pyroptosis and spontaneous regression of hemangiomas has not been documented. This is the first study to report the correlation between pyroptosis-related genes and the regression of infantile hemangiomas.

Fourteen genes that were differentially expressed between hemangiomas in the proliferation and involution stages and related to pyroptosis were identified genes that contribute to the regression of hemangiomas. GO enrichment analysis revealed that these IH-PRGs were enriched mainly in the terms “cell response to stimuli”, “inflammatory response”, “immune response”, and “PCD”. KEGG analysis revealed that these IH-PRGs were enriched mainly in the following signalling pathways: JAK-STAT, C-type lectin receptor signalling pathway, HIF-1 signalling pathway, TNF signalling pathway, FoxO signalling pathway, and other signalling pathways. The JAK-STAT pathway is extensively involved in stem cell maintenance and IH proliferation. The downregulation of the JAK-STAT pathway can lead to the depletion of hemangioma stem cells, which is an important reason for the regression of IHs.^17–19 Propranolol can promote the regression of IHs by inhibiting the JAK-STAT pathway.²⁰ However, JAK-STAT signalling regulates pyroptosis and participates in multiple disease processes in a bidirectional manner, facilitating pyroptosis in certain pathological states^21,22 while inhibiting it in other states.^23,24

We identified IL-6, EGFR, IRF1, and IL32 as hub genes that have high diagnostic value for regression-stage hemangiomas and may be potential targets for the staging, diagnosis, and treatment of hemangiomas. IL-6 is an important cytokine that regulates immune and inflammatory responses.^25,26 Compared with that in normal skin tissues, the expression of IL-6 mRNA in hemangioma tissues is increased, and IL-6 is thought to promote the proliferation of hemangioma endothelial cells.^19,27,28 However, serum IL-6 levels are increased in patients with nonprogressive hemangiomas compared with those with progressive hemangiomas,²⁹ and tissue IL-6 mRNA expression in proliferation-phase IHs is lower than that in the degeneration-phase IHs,²⁷ indicating that IL-6 is involved in hemangioma regression. IL-6 inhibits gasdermin E (GSDME)- and gasdermin D (GSDMD)-mediated pyroptosis.³⁰ EGFR is a member of the ErbB receptor family that is located on the surface of cell membranes. EGFR binds to EGF or TGF-α and is associated with the inhibition of tumour cell proliferation, angiogenesis, invasion, metastasis, and apoptosis.³¹ EGFR-targeted inhibitors can induce strong pyroptosis in lung cancer, which is an important biological mechanism for the treatment of lung cancer.³² IRF-1 is a transcription factor that can be expressed through IFN-γ, which activates the JAK/STAT1 pathway. Mono-(2-ethylhexyl) phthalate increases the expression of IRF-1, which subsequently increases Hippo–YAP signalling-mediated hemangioma cell proliferation.³³ Moreover, IRF1 is an upstream regulator of PANoptosis, and molecules in the IRF1 pathway could be targeted in the treatment of diseases associated with aberrant inflammatory cell death.³⁴ IL-32 is an inflammatory cytokine that promotes the release of various cytokines and plays an important role in the regulation of inflammatory diseases and the tumour immune microenvironment.³⁵ IL-32 can regulate pyroptosis by activating the NOD1/2/TLR4/NLRP3 pathway.³⁶

A comparison of the infiltration of immune cells in PIH and IIH tissues via single-sample GSEA revealed that the numbers of CD8⁺ T cells and Th17 cells in IIH tissues were significantly greater than those in PIH tissues, whereas the number of Tdg cells in IIH tissues was significantly lower. Interestingly, the expression of the hub genes in the regression stage was consistent with that of the differentially abundant infiltrating immune cells. IL-6 was significantly positively correlated with Th17 cell infiltration and significantly negatively correlated with Tgd cell infiltration; EGFR was negatively correlated with Tgd cell infiltration; and IRF1 and IL32 were significantly negatively correlated with Th17 cell infiltration. IH involution might be an immune-mediated process.³⁷ CD8⁺ T cells are specifically recruited to hemangioma tissue and induce apoptosis in tumour cells. Moreover, pyroptosis promotes the infiltration of CD8⁺ T cells into tumour tissue. IL6 can induce initial CD4⁺ T cells to differentiate into Th17 cells through the JAK-STAT signalling pathway, while Th17 cells secrete IL-17 and IFN-γ, which promote CD8⁺ T cell activation. Activated CTLs induce pyroptosis through the granule exocytosis pathway and play an important role in antitumour immunity.^38,39 The main T cells in the skin, γδ T cells, participate in the development of the microenvironment. Cytotoxicity is the main biological effect of γδ T cells. Additionally, γδ T cells can induce pyroptosis.⁴⁰

This study has several limitations. Despite performing a multifaceted bioinformatics analysis, further experimental research is needed to confirm the present findings. In addition, the specific mechanism underlying pyroptosis and immune cell infiltration during the regression of IHs remains unclear. The role of pyroptosis in the regression of IHs needs further exploration.

In summary, four PRGs, namely, IL6, EGFR, IRF1 and IL32, may play significant roles in the regression of IHs. This study provides insights into the molecular mechanisms underlying IH pathogenesis, highlighting the importance of pyroptosis and immune cell infiltration. The identification of diagnostic markers and potential therapeutic targets may improve the management of IHs.