The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of Wannan Medical College Second Affiliated Hospital (Approval No. WYEFYLS2025033). Anonymized clinical samples were analyzed in this study, and the ethics committee waived the need for informed consent.
Data source and preprocessing
We obtained AMI datasets (Table 1) from Gene Expression Omnibus (GEO), including GSE5986715 (training), GSE12334216 (validation), and scRNA-seq dataset GSE26926917. It was also sourced 50 TMRGs from the three Molecular Signatures Database (MSigDB)18 gene sets: KEGG_TRYPTOPHAN_METABOLISM.v2024.1.Hs.gmt, REACTOME_TRYPTOPHAN_CATABOLISM.v2024.1.Hs.gmt, and WP_TRYPTOPHAN_METABOLISM.v2024.1.Hs.gmt (Supplementary Table 1).
The raw data from the training and validation sets underwent background correction, quantile normalization, and log 2 transformation using the robust multi-array average (RMA) algorithm. Probes were annotated using the platform’s annotation files, and probes that were unannotated or had empty annotations were removed. In cases where multiple probes corresponded to the same gene, the probe with the highest expression level was retained as the representative for that gene. The preprocessed data were used for subsequent analysis.
Weighted gene co-expression network analysis (WGCNA)
To pinpoint key gene modules linked to AMI-related TMRGs, the WGCNA was applied19 to the training set. Initially, single-sample gene set enrichment analysis (ssGSEA) was used to compare TMRG scores between the two groups (AMI vs. non-AMI). Subsequently, the optimal soft threshold met the scale-free criterion (R2 > 0.9). Following this, dynamic tree pruning with a minimum module size of 50 genes identified co-expressed modules. Spearman’s correlation analysis was then used to assess the correlation between module genes and TMRG scores. Finally, the genes were labeled from the strongest-correlation module as WGCNA-TMRGs for follow-up.
Screening of candidate genes
The limma package ran differential expression analysis, identifying differentially expressed genes (DEGs)20 (AMI vs. non-AMI) under thresholds (p < 0.05, |log2FC| > 0.5). Overlapping genes were obtained by intersecting WGCNA-TMRGs with these DEGs.
Identification of key genes and construction of a diagnostic model
To identify AMI-associated key genes, four Machine Learning (ML) algorithms were applied: LASSO21, XGBoost22, random forest (RF)23, and support vector machine recursive feature elimination (SVM-RFE)24. By integrating feature selection results from multiple algorithms, bias introduced by a single model can be reduced and the robustness of candidate gene screening improved. LASSO used the Fit a cross-validated regularized regression model function to perform 10-fold cross-validation to determine the optimal λ (regularization parameter corresponding to the minimum cross-validated prediction error or binomial deviance) and to reduce model complexity through L1 regularization. XGBoost was configured with a Maximum depth of 6 per decision tree, a learning rate of 0.3 to control the contribution of each boosting iteration, and 100 boosting iterations to control overfitting through parameter constraints and to select features based on gain values. RF was configured with 1000 decision trees to reduce the risk of overfitting through ensemble learning, and feature importance was assessed using variable importance, measured by the mean decrease in Gini impurity. SVM-RFE employed 5-fold cross-validation and standard functions to perform recursive feature elimination, improving model generalization. Differences among candidate-gene sets may result from distinct feature-selection strategies, as LASSO focuses on sparse linear selection, whereas XGBoost and RF capture nonlinear relationships and feature interactions, and SVM-RFE identifies features through recursive elimination. These algorithm-specific characteristics, together with sensitivity to sample size and to correlated variables, may contribute to differences in the selected gene sets. Subsequently, the candidate genes common to all four algorithms were extracted. The expression differences of these candidate genes between the AMI and non-AMI groups were compared across the training and validation sets, and genes showing statistically significant (p < 0.05) differential expression in both sets were considered key genes. Based on this, an AMI diagnostic model was constructed using logistic regression. Model performance was assessed using the area under the receiver operating characteristic (ROC) curve (AUC), calibration plots, decision curve analysis (DCA), and confusion matrices.
Gene set enrichment analysis
To explore the potential functions of key genes, gene set enrichment analysis was performed (GSEA) via the clusterProfiler package 25. A threshold of adjusted p < 0.05 was applied, and the top 10 most significantly differentially expressed Gene Ontology (GO) biological processes (BPs) or Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were visualized.
Drug prediction and molecular docking
To screen for potential therapeutics targeting key genes, molecular docking and drug prediction were conducted. The Drug-Gene Interaction Database (DGIdb)26 served to identify drugs interacting with these genes, after which a drug-gene interaction network was constructed. Subsequently, the drugs with the highest interaction scores and their target genes were selected for molecular docking. Specifically, the 3D protein structures of key genes and drug molecular structures were obtained from PDB (https://www.rcsb.org/) and PubChem databases27, respectively. Molecular docking was then performed using a cavity-detection-guided docking server. A binding energy ≤ −5 kcal/mol indicates a stable binding affinity.
Molecular dynamics simulation
To evaluate receptor–ligand binding dynamics, molecular dynamics (MD) simulations were performed. The receptor was obtained from PDB and parameterized with Amber99SB-ILDN28. The ligand was protonated and pre-optimized in a molecular editor; partial charges were derived from quantum-chemical calculation29. Immersion of this complex in a Transferable Intermolecular Potential with Three Points (TIP3P) triangular-prism water box (10 Å buffer)30 was followed by neutralization with counterions (Na+/Cl-). Full-system minimization followed by restrained solvent/ion relaxation for energy minimization. For equilibration, the system underwent 100 ps of constant-number-of-particles, volume, and temperature (NVT) heating (0–300 K), followed by another 100 ps of constant-number-of-particles, pressure, and temperature ensemble (NPT)) density stabilization. Subsequently, a 100 ns NPT production run was performed using a 2-fs integration step, with trajectories recorded at 10-ps intervals. Particle Mesh Ewald method(PME) handled long-range electrostatics, SHAKE-constrained hydrogen bonds31, and standard thermostat/barostat schemes regulated temperature and pressure. Root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) were used to assess system stability and flexibility, respectively. Additionally, radius of gyration (Rg) and solvent-accessible surface area (SASA) were used to evaluate compactness and solvent exposure.
Immune infiltration analysis
To assess changes in the immune landscape of AMI samples, an immune infiltration analysis was performed. Deconvolution methods were used to assess the distribution of 22 immune cell types in AMI samples. Subsequently, deconvolution algorithms were used to compare immune cell infiltration between AMI and non-AMI samples. The association between infiltration scores and key gene expression was evaluated by Spearman’s correlation analysis.
ScRNA-seq data processing
Single-cell transcriptomic analysis32 was primarily performed using the R package. Cells were retained if the number of detected genes per cell was >100 and <5,000, the total unique molecular identifier (UMI) count per cell was >100 and <50,000, and the mitochondrial gene percentage was <20%. After quality control, gene expression matrices were normalized using the normalize the gene-expression count data function, with log-normalized counts, after scaling each cell to the specified total count method and a scale factor of 10,000. Highly variable genes were identified using the identify highly variable genes function, and expression values were scaled using the center and scale gene-expression values. Principal component analysis was performed, and the first 25 principal components were used for downstream analyses. Batch effects were corrected using the correct batch effects using the Harmony algorithm function. Cell clustering was performed using the Construct a nearest-neighbor graph in the selected reduced-dimensional space and the specified resolution functions, and cell populations were visualized using generate a two-dimensional uniform manifold approximation and projection representation. Cell types were annotated using established marker genes33. To prioritize biologically significant cell populations, AUCell scores were calculated to identify those most strongly associated with the target key genes. Furthermore, transcriptional regulatory mechanisms within key cell clusters were explored in depth using the Single-Cell Regulatory Network Inference and Clustering (SCENIC) method.
qPCR validation
At the Second Affiliated Hospital of Wannan Medical College (Approval No. WYEFYLS2025033), blood samples were collected from 8 AMI patients and 8 healthy controls. Due to the retrospective and anonymized nature of the archived samples, detailed demographic and clinical information was unavailable. Total RNA was extracted using an acid guanidinium thiocyanate–phenol–chloroform reagent, and concentration and purity were subsequently assessed using a microvolume UV spectrophotometer. cDNA was then generated from RNA using reverse transcriptase. After diluting all sample concentrations to 150 ng/mL, PCR amplification was performed on a real-time PCR system. The primer sequences for ADM, ECRP, MCEMP1, TSPO, PPARG, and β-actin are shown in Table 2. Relative expression levels were calculated using the 2−ΔΔCt method.
Statistical analysis
All statistical analyses were carried out using R and a commercial statistics and graphing package (see Table of Materials), whereas network visualization was performed in network visualization software. For high-throughput analyses involving multiple comparisons, including GSEA, multiple-testing correction was performed using the Benjamini-Hochberg false discovery rate (FDR) method, and adjusted p-values were used where appropriate. Group comparisons were conducted via Student’s t‑test (normal) or Mann‑Whitney U test (non‑normal). Significance: two‑tailed p < 0.05.