CCL4 as a Potential Immune-Metabolic Biomarker in Acute Myocardial Infarction via Machine Learning and Single-Cell Transcriptomics

Podsumowanie

This study integrates bulk and single-cell transcriptomics with machine learning to identify CCL4 and five other genes as potential diagnostic biomarkers for acute myocardial infarction, revealing immune-metabolic crosstalk and cell-type-specific expression patterns.

Streszczenie

Acute myocardial infarction (AMI) is a leading cause of morbidity and mortality worldwide, highlighting the need for novel complementary biomarkers. By integrating bulk and single-cell transcriptomic data with machine learning approaches, calmodulin-related genes associated with AMI and explored their immune-metabolic features were identified. Differential expression and weighted co-expression analyses revealed 60 calmodulin-related genes, from which six key genes (SOCS3, GBP4, ST14, KPNA5, STAB1 and CCL4) were screened using multiple machine learning algorithms and validated in independent datasets. Functional analyses indicated enrichment in immune, inflammatory, and metabolic pathways. Immune infiltration and single-cell transcriptomics showed cell-type-specific expression patterns, with CCL4 predominantly expressed in T and NK cells and markedly reduced in AMI samples. qPCR validation confirmed significant expression changes for four of the six genes in the local cohort. Drug-gene interaction and docking analyses suggested candidate compounds for further investigation. Collectively, the findings suggest that CCL4 may serve as a potential diagnostic biomarker for AMI, and the observed immune-metabolic associations provide a basis for future mechanistic and translational studies.

Wprowadzenie

Acute myocardial infarction (AMI), caused by sustained myocardial ischemia, is among the most lethal cardiovascular diseases worldwide¹˒². Its pathogenesis is multifactorial, with major contributors including atherosclerosis, thrombosis, coronary artery spasm, and arrhythmias³˒⁴. Hypertension, diabetes, hyperlipidemia, smoking, obesity, and psychological stress are also recognized as important risk factors⁵. Although advances in interventional techniques, antiplatelet therapies, and thrombolytic strategies have markedly improved short-term survival, the overall incidence of AMI continues to rise, particularly among younger populations⁶. Cardiac troponin (cTn) is currently the gold-standard biomarker for diagnosing AMI. It is a structural protein released into the circulation following myocardial cell necrosis and exhibits high sensitivity and specificity. However, cTn levels may be influenced by renal function and the timing of sample collection⁷˒⁸. Therefore, there remains a need to develop new complementary biomarkers.

Calmodulin (CaM) is a small, ubiquitously expressed, and highly conserved calcium-binding protein that functions as a key intracellular transducer of calcium signals and participates in a broad range of biological processes⁹˒¹⁰. By sensing changes in intracellular Ca²⁺ concentrations, CaM interacts with numerous target proteins, including kinases, phosphatases, and ion channels, to regulate cell proliferation, apoptosis, metabolism, muscle contraction, and inflammatory responses¹¹. In the cardiovascular system, CaM plays a central role in controlling cardiomyocyte contractility, heart rate, electrical stability, and vascular smooth muscle tone¹²˒¹³. Accumulating evidence indicates that CaM dysfunction or dysregulation of CaM-related signaling pathways contributes to several cardiovascular disorders, including arrhythmias, heart failure, hypertension, and cardiac hypertrophy¹⁴˒¹⁵. Furthermore, CaM participates in cardiomyocyte responses to ischemia, oxidative stress, and inflammatory stimuli, suggesting its potential involvement in the pathological progression of AMI¹⁶. Therefore, investigating CaM in the context of cardiovascular disease is of substantial scientific and clinical interest.

Bulk RNA-seq and single-cell RNA-seq datasets were integrated to systematically identify differentially expressed genes (DEGs) significantly associated with the CaM score. Weighted gene co-expression network analysis (WGCNA) was used to further identify CaM-related candidate genes involved in AMI. Multiple machine learning algorithms were then applied to screen for key genes and to construct predictive models to assess individual AMI risk. Transcriptomic features, immune infiltration profiles, and competing endogenous RNA (ceRNA) regulatory networks were also integrated to investigate how these key genes may contribute to immune regulation and potential drug-target interactions. The findings indicate that several core genes are involved in immunometabolic crosstalk pathways and exhibit predictive value and therapeutic potential. This integrated approach refines the molecular understanding of AMI and provides a rationale for biomarker development and targeted intervention.

Protokół

The study was conducted in accordance with the Declaration of Helsinki. The protocol was approved by the Ethics Committee of Anhui Public Health Clinical Center on September 19, 2025 (approval ID: PJ-YX2025-062). Written informed consent was obtained from all participants before blood collection. The local cohort included eight patients with acute myocardial infarction (AMI) and eight healthy controls. The research tools used in the protocol are listed in the Table of Materials.

1. Data sources and processing

Bulk RNA sequencing datasets related to myocardial infarction were obtained from the Gene Expression Omnibus (GEO). The GSE59867 and GSE48060 datasets were used for bulk transcriptomic analyses, and the GSE269269 single-cell RNA sequencing dataset was used for cellular-level analyses (Table 1). A set of 255 calmodulin-related genes was obtained from the Human Protein Atlas for subsequent gene-set analyses.

Dataset	Sample type	Sample (controls)	Sample (patients)	Sequencing platform
GSE59867	Bulk RNA-seq	46	111	GPL6244
GSE48060	Bulk RNA-seq	21	31	GPL570
GSE269269	scRNA-seq (peripheral blood)		10	GPL24676

Table 1: Characteristics of the datasets used in the study. The table lists the dataset accession numbers, sample types, numbers of control and patient samples, and sequencing platforms for the bulk RNA and single-cell RNA sequencing datasets. RNA-seq, RNA sequencing; scRNA-seq, single-cell RNA sequencing.

Inter-sample variability in the bulk transcriptomic datasets was corrected using the normalizeBetweenArrays function in the limma package, version 3.60.6. Differential gene expression analysis was then performed using limma. Differentially expressed genes (DEGs) were defined using the thresholds P < 0.05 and |log₂ fold change| > 0.5. The resulting DEGs were visualized using volcano plots and heatmaps and classified as significantly upregulated, significantly downregulated, or not significantly changed.

2. Single-sample gene set enrichment analysis and weighted gene co-expression network analysis

Single-sample gene set enrichment analysis (ssGSEA) was performed using the 255 calmodulin-related genes. The GSVA package, was used to calculate a calmodulin-related gene score, designated Calmodulin_score, for each AMI and control sample. Differences in Calmodulin_score between the AMI and control groups were evaluated using the Wilcoxon rank-sum test.

Weighted gene co-expression network analysis (WGCNA) was performed using the bulk transcriptomic data from patients with AMI. Genes with a mean fragments per kilobase of transcript per million mapped reads value of ≤0.5 were excluded. Samples were clustered to identify and remove outliers.

A soft-thresholding power that achieved a scale-free topology fit of R² > 0.8 was selected. A topological overlap matrix was then constructed. Gene modules were identified using the dynamic tree-cut algorithm with a minimum module size of 200. Modules with highly similar eigengenes were merged using a correlation threshold of >0.75, corresponding to a module-merging threshold of 0.25.

Relationships between module eigengenes and clinical traits, including Calmodulin_score, were assessed using Pearson correlation analysis. The resulting module-trait relationships were displayed in a heatmap annotated with correlation coefficients and corresponding P values. Module membership and gene significance were calculated for each gene. Scatter plots of module membership against gene significance were generated to identify genes with high intramodular connectivity and trait relevance.

3. Identification of AMI-associated calmodulin-related genes

AMI-associated calmodulin-related genes were identified by intersecting the DEGs with genes from WGCNA modules significantly correlated with Calmodulin_score. The overlapping genes were retained for downstream analyses.

Functional enrichment analysis was performed using the clusterProfiler package. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes annotations were used to identify biological processes, molecular functions, cellular components, and signaling pathways associated with the overlapping genes.

4. Identification and validation of key genes using machine learning

Univariate logistic regression analysis was performed using the overlapping genes. Three machine learning algorithms were then applied independently using the following R packages and parameters: random forest, XGBoost and support vector machine.

Each algorithm was used to prioritize genes with predictive value for AMI. Candidate key genes were defined as the genes identified by all three algorithms. Genes that also exhibited significant and directionally consistent differential expression in both the GSE59867 training dataset and the GSE48060 external validation dataset were retained as final key genes.

5. Construction and evaluation of diagnostic models

A logistic regression model was constructed using the identified key genes and the lrm function. A nomogram was generated using the regplot function to display the contribution of each feature gene to the predicted probability of AMI.

Model discrimination was evaluated by receiver operating characteristic analysis using the pROC package. The area under the receiver operating characteristic curve was calculated to assess the model's ability to distinguish AMI from control samples.

Calibration curves were generated to compare predicted probabilities with observed outcomes. Decision curve analysis was performed to estimate the net clinical benefit of the model across a range of threshold probabilities.

6. Gene set enrichment analysis and competing endogenous RNA network construction

Gene set enrichment analysis was performed separately for each key gene using its gene-correlation matrix and the clusterProfiler package. Kyoto Encyclopedia of Genes and Genomes pathway enrichment results were ranked by the absolute normalized enrichment score. The five highest-ranking pathways were displayed for each gene.

Functional associations and gene-interaction networks were examined using GeneMANIA. Potential microRNA regulators of the key genes were predicted using miRanda, miRTarBase, TargetScan, and miRDB. Candidate microRNA-mRNA interactions were identified by intersecting the predictions from the four databases.

Long noncoding RNA-microRNA interactions were obtained from spongeScan. The long noncoding RNA-microRNA and microRNA-mRNA relationships were then integrated to construct a competing endogenous RNA regulatory network. The network was visualized as a Sankey diagram using the ggalluvial package.

7. Drug prediction and molecular docking

Drug-gene interactions were predicted using the Drug-Gene Interaction Database. The resulting interaction network was visualized using network-analysis software.

The UniProt protein identifier for CCL4 was retrieved as P13236. The corresponding three-dimensional protein structure was obtained in Protein Data Bank format (PDB) under the accession number 1HUM (human MIP-1β, X-ray diffraction structure), which was selected for docking. Chain A, representing the biologically relevant monomer, was selected for docking. Protein preparation was performed using the Prepare Protein module in CB-Dock2, which includes removal of water molecules, addition of polar hydrogens, and assignment of Gasteiger charges. The three-dimensional chemical structures of the candidate compounds (clodronic acid and epoetin alfa) were retrieved from the PubChem database in Structure-Data File (SDF) format. Docking simulations were performed using the CB-Dock2 online platform, which employs the AutoDock Vina algorithm for blind docking. The docking site was set to cover the entire protein surface to allow unbiased identification of potential binding pockets. Binding affinity was calculated as the predicted free energy of binding (ΔG) in kcal/mol. Final docking poses and protein-ligand interactions (e.g., hydrogen bonds, hydrophobic contacts) were visualized using PyMOL and CB-Dock2's built-in interaction viewer.

8. Single-cell RNA sequencing data preprocessing

Quality control was performed before downstream single-cell RNA sequencing analysis. Cells were retained when the number of detected genes was between 200 and 10,000, the total unique molecular identifier count was ≥1,000, and the proportion of mitochondrial transcripts was ≤20%.
Cells expressing fewer than 200 genes and genes detected in fewer than three cells were excluded. These filters were applied to reduce the inclusion of low-quality cells and technical noise. Gene-expression values were normalized using the NormalizeData function in the Seurat package. Highly variable genes were identified using the FindVariableFeatures function. Expression values for the highly variable genes were centered and standardized using the ScaleData function.

Batch effects associated with experimental or sequencing variation were corrected using the RunHarmony function from the Harmony integration framework¹⁷.

9. Single-cell dimensionality reduction, clustering, and annotation

Principal component analysis was first applied to reduce the dimensionality of the single-cell RNA sequencing dataset. Uniform manifold approximation and projection, and t-distributed stochastic neighbor embedding were subsequently used to visualize cellular heterogeneity.

Transcriptionally similar cells were grouped using the FindNeighbors and FindClusters functions in Seurat. Differentially expressed marker genes for each cluster were identified using the FindAllMarkers function by comparing each cluster with all remaining clusters.

Cell types were assigned using canonical marker genes obtained from published literature and established cell-marker databases. The spatial distribution and expression levels of key genes were visualized using the FeaturePlot function¹⁸.

10. Quantitative polymerase chain reaction analysis

Peripheral blood samples were obtained from 8 patients with AMI and 8 healthy controls at the Anhui Public Health Clinical Center. The AMI group included patients diagnosed according to the Fourth Universal Definition of Myocardial Infarction, with symptoms consistent with myocardial ischemia and elevated cardiac troponin I levels above the 99th percentile upper reference limit. The control group comprised age- and sex-matched healthy individuals with no history of cardiovascular disease, normal electrocardiograms, and no abnormalities in routine blood tests, liver function, or renal function. For patients with AMI, 3 mL of ethylenediaminetetraacetic acid-anticoagulated blood was collected within 24 h of hospital admission. The same volume was collected from healthy controls during the corresponding study period.

Total RNA was isolated from peripheral blood according to the protocol supplied with the blood RNA isolation kit. RNA concentration and purity were assessed using a NanoDrop spectrophotometer, and RNA integrity was verified by agarose gel electrophoresis. Only samples with an A260/A280 ratio between 1.8 and 2.1 were used for subsequent analyses. A total of 500 ng of RNA was reverse-transcribed into complementary DNA using a first-strand complementary DNA synthesis reagent. The resulting complementary DNA was diluted to a final concentration of 150 ng/mL. Quantitative polymerase chain reaction amplification was performed in a total reaction volume of 10 µL using a SYBR Green-based master mix without passive reference dye. All qPCR reactions were performed in technical duplicates twice, and the subsequent calculations were based on the average Ct values.

Amplification was performed using a real-time polymerase chain reaction instrument. The cycling conditions consisted of initial denaturation at 95 °C for 5 min, followed by 40 cycles of denaturation at 95 °C for 10 s, annealing at 60 °C for 30 s, and extension at 72 °C for 30 s. Melt-curve analysis was performed after amplification.

Gene-expression levels were normalized to β-actin. Relative expression was calculated using the 2^−ΔΔCt method. 

11. Statistical analysis

Statistical analyses were performed in R. Network visualizations were generated using network analysis software. Differences between two groups were evaluated using the Wilcoxon test unless otherwise specified. Continuous variables with a normal distribution were compared using Student’s t-test. Non-normally distributed continuous variables were compared using the Mann-Whitney U test, also referred to as the Wilcoxon rank-sum test. All statistical tests were two-tailed. A P value of <0.05 was considered statistically significant¹⁹.

Wyniki

Identyfikacja i analiza wzbogacenia genów związanych z kalmoduliną powiązanych z AMI
W pierwszej kolejności w zbiorze danych GSE59867 przeprowadzono korekcję efektu serii, aby zapewnić porównywalność między próbkami (Rycyna 1A). Następnie zidentyfikowano łącznie 168 genów różnicowo wyrażonych (DEG), w tym 77 genów o zwiększonej ekspresji i 91 genów o zmniejszonej ekspresji w próbkach AMI w porównaniu z grupą kontrolną (Rycyna 1B). Analiza wz...

Dyskusja

Acute myocardial infarction (AMI) results from acute coronary artery occlusion, leading to myocardial ischemic necrosis and high morbidity and mortality worldwide. Calmodulin, a calcium-dependent regulatory protein, plays a pivotal role in intracellular signal transduction. Previous studies have demonstrated that calmodulin regulates cardiomyocyte apoptosis, inflammatory responses, and calcium homeostasis during AMI, suggesting that it may be an important molecule in AMI pathogenesis²⁰. In the pre...

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
β-actin primer	Tsingke	N/A	Forward: 5′-CATGTACGTTGCTATCCAGGC-3′ Reverse: 5′-CTCCTTAATGTCACGCACGAT-3′
calibrate	rms	N/A	Parameters: lrmModel, method = "boot", B = 1000
CB-Dock2	CB-Dock2 Server	Online tool	https://cadd.labshare.cn/cb-dock2/
CCL4 primer	Tsingke	N/A	Forward: 5′-CTGTGCTGATCCCAGTGAATC-3′ Reverse: 5′-TCAGTTCAGTTCCAGGTCATACA-3′
CIBERSORT	Stanford University	Online tool	Parameters: function(sig_matrix, mixture_file, perm = 0, QN = TRUE)Spearman correlation analysis was performed between immune-cell abundance and gene expression.
clusterProfiler	Bioconductor	4.12.6	https://bioconductor.org/packages/clusterProfiler
Cytoscape	Cytoscape Consortium	3.9.1	https://cytoscape.org/
Drug-Gene Interaction Database	DGIdb	Online database	https://www.dgidb.org/
Gene Expression Omnibus	NCBI	Online database	https://www.ncbi.nlm.nih.gov/geo/
Gene Ontology database	Gene Ontology Consortium	Online database	http://geneontology.org/
GeneMANIA	University of Toronto	Online tool	https://genemania.org/
GBP4 primer	Tsingke	N/A	Forward: 5′-AGGCTGCTAAAACACAAGCTG-3′ Reverse: 5′-CCCCAGGTAGAGTGACAATCAT-3′
ggalluvial	CRAN	0.12.5	Sankey diagram generation
GSVA	Bioconductor	1.52.3	Parameter: P < 0.05
Harmony integration framework	CRAN	1.2.0	https://github.com/immunogenomics/harmony
Hifair III 1st Strand cDNA Synthesis SuperMix	YEASEN	11141ES10	Reverse-transcription reagent
Hieff qPCR SYBR Green Master Mix	YEASEN	11201ES	qPCR reagent
Human Protein Atlas	Human Protein Atlas Consortium	Online database	https://www.proteinatlas.org/
KPNA5 primer	Tsingke	N/A	Forward: 5′-TCAAGAGATGCGTAGACGAAGA-3′ Reverse: 5′-ACATTTCTGCGTTTGAACAACTG-3′
Kyoto Encyclopedia of Genes and Genomes database	Kanehisa Laboratories	Online database	https://www.kegg.jp/
LightCycler 480 Instrument II	Roche	LightCycler 480 II	Real-time PCR instrument
limma	Bioconductor	3.60.6	Differentially expressed gene screening
lrm	rms	N/A	Parameters: Type ~ SOCS3 + GBP4 + ST14 + KPNA5 + STAB1 + CCL4, data = aSAH, maxit = 100
miRanda	MicroRNA.org	Online tool	http://www.microrna.org/microrna/home.do
miRDB	miRDB	Online tool	http://mirdb.org/
miRTarBase	National Chiao Tung University	Online database	https://mirtarbase.cuhk.edu.cn/
Protein Data Bank	RCSB	Online database	https://www.rcsb.org/
PubChem	NCBI	Online database	https://pubchem.ncbi.nlm.nih.gov/
PyMOL	Schrödinger, LLC	2.5.4	Visualization of molecular docking results
pROC	CRAN	1.18.5	Receiver operating characteristic curve generation and visualization
Random forest package	CRAN	3.3.1	Parameters: ntree = 100, seed = 200
R software	R Foundation	4.2.2	Bioinformatics and statistical analyses
regplot	CRAN	1.1	Parameters: glm(Status ~ gene, family = binomial(), data = merged_data)lrm(Type ~ SOCS3 + GBP4 + ST14 + KPNA5 + STAB1 + CCL4, data = aSAH, maxit = 100)
rmda	CRAN	1.6	Parameters: decision_curve(Type ~ SOCS3 + GBP4 + ST14 + KPNA5 + STAB1 + CCL4, data = aSAH, family = binomial(link = "logit"), thresholds = seq(0, 1, by = 0.01), confidence.intervals = 0.95)
RNA isolater MolPure Blood RNA Kit	YEASEN	19241ES50	RNA isolation kit
Seurat	Satija Lab	5.1.0	Marker genes used for cell annotation:Monocytes: FCN1, S100A9, S100A8B cells: CD79B, MS4A1T cells: LDHB, CD4, CD3DNK cells: FGFBP2, GZMB, NKG7, KLRD1Megakaryocytes: CD163, CD36, PF4, PPBPBlood cells: ALAS2, AHSP, CA1Plasma cells: IGHA1, CD79A, MZB1, JCHAINDendritic cells: HLA-DMB, HLA-DPA1, HLA-DQB1Mast cells: FCER1A, MS4A2
SOCS3 primer	Tsingke	N/A	Forward: 5′-CCTGCGCCTCAAGACCTTC-3′ Reverse: 5′-GTCACTGCGCTCCAGTAGAA-3′
spongeScan	spongeScan	Online tool	http://spongescan.rc.ufl.edu/
ST14 primer	Tsingke	N/A	Forward: 5′-TTCCTGCCAGTCAACAACGTC-3′ Reverse: 5′-GGTACTGCAAATGCCACACC-3′
STAB1 primer	Tsingke	N/A	Forward: 5′-CCGGGAAATCCTTACCACAGC-3′ Reverse: 5′-ACCTTCGTGTTTGTTGGGTCC-3′
Support vector machine package	CRAN	1.7-13	Parameters: input, k = 10, halve.above = 100
survminer	CRAN	0.5.0	Survival curve plotting
TargetScan	TargetScan	Online tool	http://www.targetscan.org/
UniProt	UniProt Consortium	Online database	https://www.uniprot.org/
Weighted gene co-expression network analysis package	CRAN	1.73	Parameters: R² > 0.8, minModuleSize = 200
XGBoost package	CRAN	1.7.8.1	Parameters: nrounds = c(50, 200), max_depth = c(3, 8), eta = c(0.01, 0.3)