The datasets analyzed in this study were downloaded from the GEO database. The GSE12720 (20) contains 21 preoperative samples and 21 post reperfusion samples, GSE112713 (21) contains 11 preoperative samples and 11 post reperfusion samples, and GSE14951 (22) contains 5 preoperative samples and 5 post reperfusion samples. They were merged as a training set by the “ComBat” function from “sva” package (23). Principal component analysis (PCA) was performed to verify the batch effect before and after the data operation. PCA is a method of data dimensionality reduction, which extracts the feature vectors (components) of data from high latitude data, converts them into low dimensional data, and displays these features using two-dimensional or three-dimensional graphs. The GSE23649 (33 LIRI and 33 control samples), GSE15480 (6 LIRI and 6 control samples), and GSE151648 (24) (40 LIRI and 40 control samples) datasets were used as validation sets for subsequent biomarker validation.

The PCD-related genes were carefully hand-curated from a variety of reliable sources, such as the KEGG database, the MSigDB database, review articles, and manual curation (8, 25–29). A total of 1,567 genes associated with 19 types of PCD were eventually included, as shown in Supplementary Table S1.

Using the “limma” software, the criteria for detecting differentially expressed genes (DEGs) between LIRI and normal samples were set to |log2 fold change (FC)| > 0.585 and adjust p < 0.05 (18). The expression differences for LIRI were visualized by a volcano plot using the “ggplot2” package in R and a heatmap using the “pheatmap” package in R. The proportions of various immune cell types were evaluated using ssGSEA from the “GSVA” package in order to compare the differences in immunity status between groups (30). Gene Set Enrichment Analysis (GSEA) is a computational method used to determine differences in biological process activity or pathway enrichment across different samples. It evaluates whether predefined sets of genes are statistically significantly associated with a biological state. In this study, all genes were first categorized into two groups based on their positive and negative logFC values. We then performed enrichment analysis on the positive and negative logFC gene groups using the “clusterProfiler” package in R (31).

The Venn map was applied to select the common genes between DEGs and PCD-related genes. The STRING database¹ was used to investigate protein interactions, with the validity of such interactions being determined by a composite score greater than 0.4. The GeneMANIA database² prioritized genes for functional tests. Finally, the Metascape database³ was used to conduct functional enrichment analysis (32).

Different machine learning methods have different characteristics in data processing and pattern recognition, and their combined use can more comprehensively mine the potential information in the data, while avoiding the bias that may be brought by relying on a single method. We used 12 machine learning methods (Lasso, Ridge, Stepglm, XGBoost, Random Forest, Enet, plsRglm, GBM, Naïve Bayes, LDA, glmBoost, and SVM) to build the diagnostic model. 113 combinations of these 12 techniques were investigated through methodical designs for model generation and variable selection utilizing the training dataset (33). Three datasets were then used as external test datasets to validate the model’s performance. We also calculated the AUC values of each model and output the AUC matrix, and then used “ComplexHeatmap” to draw the AUC heatmap for visualizing the performance of each model on the training and testing sets. The best model is the one with the highest average AUC over training and testing cohorts (29).

To identify the hub gene, the cytoHubba plug in Cytoscape was applied (34). GSEA and GSVA were utilized to elucidate the molecular mechanisms between high- and low-hub gene expression samples (35). Meanwhile, the ssGSEA was also used to assess the correlation between significantly different enriched immune cell types and the hub gene. These gene sets used in these analyses were obtained from MSigDB.

We utilize TFTF online (36) to predict the hub gene for TFs based on the integration of five major TF-Target online tools, namely hTFtarget (37), ENCODE (38), JASPAR (39), GTRD (40), and ChIP_Atlas (41). In addition, the TargetScan, miRDB, and miRanda databases were used to anticipate miRNA-mRNA pairs in order to identify the ceRNA network that might be influenced by model genes (42). Genes that were simultaneously listed in three databases were the only ones that were thought to be possible mRNA targets for further research (43). To predict miRNA-lncRNA pairs, the spongeScan database was used. At last, the ceRNA network could be seen using Cytoscape (44). Meanwhile, the Human Protein Atlas (HPA⁴) was utilized to examine the hub gene’ immunofluorescence and single-cell type atlases.

The CMap database⁵ can link diseases, genes, and drugs based on gene expression profiles. The upregulated DEGs of LIRI were entered into the CMap database to identify potential small-molecule compounds for the treatment of LIRI (43). Then, the protein structures of the feature genes were obtained from the PDB database, and the AutoDock tool was utilized to compute the hydrogenation and charge of proteins. PubChemdatabase to download the chemical structure of the drug’s active ingredient. The AutoDock is used to examine the rotatable bonds and charge balance of small molecules. Finally, PyMol software was used to check the docking complex.

The flow chart of the research is presented in Figure 1. Initially, we merged three datasets of LIRI and corrected batch effects using the “sva” software package. As can be seen from Supplementary Figures S1A–D, after data normalization, the difference between batches is effectively eliminated, and the three datasets can be merged. The volcano map and heatmap showed the DEGs for the LIRI (Figure 2A and Supplementary Figure S1E). Given the important role of immune cells in LIRI, differences in immune cell infiltration were subsequently compared between LIRI patients and controls. It was found that LIRI patients had a higher percentage of activated CD4 T cells, dendritic cells, eosinophils, MDSC, mast cells, NK T cells, neutrophils, T helper cells, and CD8 T cells, while central memory CD4 T cells and NK cells were in contrast (Figure 2B). Further GSEA analysis was then performed to explore the differences in biological function between LIRI patients and controls. The results showed that relevant pathways and functions such as apoptotic and cytokines were significantly upregulated in the LIRI group (Figures 2C,D), while metabolism-related pathways and functions were significantly upregulated in the control group (Figures 2E,F).

To further explore the role of PCD in LIRI, we collected PCD-related genes from 19 types of PCD. Using the Venn diagram, the DEGs and PCD-related genes were crossed, and 47 PCD-related DEGs were obtained (Figure 3A). The PPI network complex was constructed by importing 47 PCD-related DEGs into the STRING database. PPI analysis shows tight junctions between these genes (Figure 3B). Then, we used GeneMANIA to predict correlations between colocalization, pathways, shared protein domains, co-expression, and prediction. The network showed that these genes were mainly enriched in the regulation of apoptotic signaling pathway, response to molecules of bacterial origin, and cell chemotaxis (Figure 3C). Finally, the Metascape database was used to explore the biological activities and pathways of PCD-related DEGs. The result showed that they were involved in the regulation of PCD, the TNF signaling pathway, and the regulation of hydrolase activity (Figure 3D).

A total of 12 machine learning algorithms were combined to identify the most robust diagnostic model based on 47 PCD-related DEGs. Among them, the Lasso + XGBoost machine learning method has the highest average AUC value in the training set and the three validation sets, so the Lasso + XGBoost model is chosen as the final prediction model (Figure 4A and Supplementary Figures S2A–D). We identified 11 model genes: ANGPTL4, ATF3, BAG3, IL1B, IL6, IRF1, KRT18, MYC, PFKFB3, TNFAIP3, and TNFRSF1B. ROC curves of 11 model genes were shown in Figure 4B. The confusion matrix shows that our diagnostic model has high sensitivity, specificity and accuracy (Figure 4C), and it was confirmed in the verification sets (Figures 4D–F).

We then performed gene expression analysis for these 11 model genes in the training set and validation sets, respectively. Analysis in the training set showed that these 11 model genes were all highly expressed in LIRI patients (Figure 5A). The gene expression analysis of the validation set was shown in Supplementary Figures S3A–C, which was consistent with the training set. Figure 5B demonstrated the correlation between these 11 model genes. To further explore the most critical of these genes, we examined the PPI networks of model genes using the STRING tool (Figure 5C) and inserted the results into Cytoscape. MYC was detected as the hub gene for these model genes using the Cytohubba plugin (Figure 5D).

Next, the relationship between MYC and the different immune cells was further studied. It showed that MYC was significantly positively correlated with Type 2 T helper cells, mast cells, macrophages, eosinophils and activated CD4 T cells, while oppositely in monocytes (Figure 5E). Then we performed GSVA analysis to further explore the biological functions of hub genes in LIRI. MYC was significantly positively associated with circadian rhythm, regulation of autophagy, and regulation of immune response-related pathways and functions in LIRI samples (Figures 6A,B). Subsequently, GSEA analysis was utilized to further explore the possible mechanism of action of MYC in LIRI samples. It showed that MYC was significantly positively associated with B cell receptor, Toll like receptor, and transcription regulator complex-related pathways and functions, while oppositely in ECM receptor interaction, complement and coagulation cascades, and basement membrane-related pathways and functions (Figures 6C–F).

TFs are critical in regulating gene expression and shaping cellular and organismal phenotypes. We therefore combined five databases to characterize TFs in MYC and identified a total of 19 reliable TFs: CTCF, BATF, EGR1, MAZ, PAX5, SP1, TCF12, E2F6, FOS, JUND, TEAD4, NFIC, PBX3, STAT3, TFAP2A, TFAP2C, TCF7L2, FOSL2, and RELA (Figure 7A). In addition, the ceRNA network of MYC was built using various public databases. Ultimately, we were able to identify 2 objective miRNAs and 12 objective lncRNAs of MYC (Figure 7B). The network disclosed transcriptional regulatory mechanisms for MYC. The cellular location of MYC was then investigated. It can be found that MYC is mainly expressed in the nucleoplasm, where green represents the target protein and red represents microtubules (Figure 7C). In addition, we evaluated MYC expression at the single-cell data level using the HPA database. It was found that MYC was mainly distributed in smooth muscle cells and fibroblasts (Figure 7D).

In addition, to search for potential drugs to treat LIRI, we entered upregulated DEGs in LIRI into the CMap database for analysis. The top 10 small-molecule compounds are shown in Table 1 and are considered potential therapeutic drugs for LIRI. Then we conducted molecular docking analysis of these 10 small-molecule compounds with MYC. The minimum binding energy between the hub gene and each small-molecule compound is shown in Table 1. The results showed that the minimum binding energy between the ligand and the receptor was less than −7.0 kcal/mol, so the target proteins had a good affinity for the active ingredient, suggesting that these small-molecule drugs could act on MYC to treat LIRI. Among them, BMS-536924 (−9.8 kcal/mol) and PF-431396 (−9.4 kcal/mol) have the lowest free binding energy for MYC and are most likely to be potential therapeutic agents for LIRI (Figures 7E,F).