Identification of diabetes related phenotype and diagnostic biomarkers in coronary artery disease

Data sources and preparation

The datasets concerning diabetes (2025) and CAD (Chen et al., 2023) predominantly derive from the GEO database. The diabetes datasets include GSE20966 (20 samples with diabetes), GSE25724 (six of 13 total samples have diabetes), and GSE38642 (nine of 63 total samples have diabetes), and the CAD datasets include GSE12288 (222 samples in total, with 110 samples in CAD), GSE20680 (195 samples with CAD), GSE20681 (198 samples in CAD), and GSE42148 (24 samples, with 13 samples in CAD). The single cell dataset was GSE32678. Additionally, to increase the sample size and confidence level, various datasets of patients with diabetes mellitus and coronary heart disease were merged, with batch effects excluded by “combat” packages for optimal analysis. For missing values, we fill them with the mean value. For data cleansing, we correct formatting errors, remove duplicate values, and remove unneeded data to ensure data cleanliness. For data with large variance, we usually use the “normalize” function to normalize the data, making the data comparable. In addition, for outliers, we usually remove them or fill them with mean values. Specific dataset inclusions and workflow are detailed in the supplementary material.

WGCNA analysis and differential analysis to identify diabetes related markers

Weighted correlation network analysis (WGCNA) is a systems biology approach for characterizing gene association patterns between samples (Langfelder & Horvath, 2008). This analysis focuses on identifying co-expressed genes based on gene-set introgression and gene-set-phenotype associations. We utilized WGCNA to pinpoint modules co-expressed with diabetes and the core genes in the network.

Identification of diabetes phenotypes in CAD patients

The intersecting genes from WGCNA and differential expression analysis were employed to identify immune-related diabetes phenotypes. We identified these phenotypes using the unsupervised clustering algorithm (consensus clustering). The main basis for the classification is the CDF value and delta region in the consensus clustering algorithm, when the CDF curve is the most downward and the delta region has the largest change, it means that this algorithm has the best classification effect.

Identification diabetes related coronary artery disease diagnostic system

For more efficient results, we applied the LASSO regression algorithms to reduce dimensionality and identified 16 diabetes-related genes. These genes were used to construct the diagnostic monitoring model for CAD and to measure the area under the curve (AUC) to evaluate the diagnostic capacity. The model was constructed by “glmnet” package (seed(123), alpha = 1, nfold=100).

Single-cell landscape analysis

The “Seurat” package was applied for single-cell sequencing analysis. The criteria for cell screening were defined as nCount_RNA >100, nFeature_RNA <5,000, percent.mt <25, and nFeature_RNA >100. Only the top 3,000 variable genes were then chosen for additional examination. Using principal component analysis (PCA), significant dimensions with p < 0.05 were found. Dimensionality reduction was made possible via the t-distributed stochastic neighbor embedding (tSNE) approach, which was combined with an 11-cluster classification analysis conducted on all cells. The “singleR” package was used to annotate cell clusters based on trends in the composition of marker genes.

Vertebrate animal study methods

This study was approved by Animal Experimental Medicine Ethics Committee of the First Affiliated Hospital of Xinjiang Medical University, Approval number: 20210301-215. Twenty-eight-week-old male C57/B6J wild-type and LDLR−/− mice were obtained from the Xinjiang Medical University experimental animal center’s specialized pathogen-free (SPF) barrier system. A high-fat diet consisting of 78.85% chow diet, 21% lard, and 0. 15% cholesterol was provided to the model group whereas the control group was fed chow diet. An injection of pentobarbital was used to induce anesthesia. We use carbon dioxide for the euthanasia of mice.

Mouse euthanasia operation procedures

We placed the mice to be executed in empty boxes and used carbon dioxide asphyxiation to execute them. Criteria established for euthanizing animals prior to the planned end of the experiment was the State Standard of the People’s Republic of China, GB/T 39760-2021.

Gene silence of KCNQ1 and ITPK1

We knocked down of KCNQ1 and ITPK1 through siRNA. The sequence was as follow: KCNQ1(SiRNA): CCGCCUGAACCGAGUAGAATT; ITPK1(SiRNA): AUCUUGUCCAGCUCCGUCATT.

CAD mouse model construction

Twenty-eight-week-old male C57/B6J wild-type and LDLR−/− mice were grown in the Xinjiang Medical University experimental animal center’s specialized pathogen-free (SPF) barrier system. Following randomization, a high-fat diet consisting of 78.85% chow diet, 21% lard, and 0. 15% cholesterol was provided to the CAD model group (LDLR−/− mice) continuously for 12 weeks, whereas the control group was fed chow diet. To create a model of a myocardial infarction, the anterior descending branch of the mice was clamped for seven days at the conclusion of the twelfth week. Changes in electrocardiogram before and after modeling were observed, and blood lipids and cardiac enzymes were measured 2 weeks after surgery to determine whether coronary artery disease occurred in mice. Sections of mouse heart tissue were obtained for immunohistochemistry (IHC) staining. The IHC staining of tissue slices from the cardiac tissues of the two groups (five mice in each group) was then examined (Fig. 1).

Flow cytometry analysis

Apoptosis rates were determined using the Annexin V/PE Apoptosis Detection Kit (BD Biosciences, Franklin Lakes, NJ, USA). Cells were collected with 1 × binding buffer. Subsequently, five µL of Annexin V-PE and five µL of 7-AAD were added, and the incubation was carried out for 5 min in the absence of light. The analysis of apoptosis was finally conducted using a FACSCalibur ffow cytometer.

Immunohistochemistry staining

Detailed procedures for immunohistochemistry can be found in the classic literature on Tavitian et al. (2025). Dimethylbenzene was used to dewax pathological sections, and gradient alcohols were used to hydrate them. Sections were then subjected to boiling sodium citrate (pH 6.0) for antigen retrieval and 3% H2O2 to deactivate endogenous peroxidase. Sections were then treated in turn with primary antibody, HRP-linked secondary antibody, and 5% goat serum. A DAB substrate was used to further visualize the indicated antibody staining. Sections were then examined under a bright-field microscope after being restained with hematoxylin and sealed with neutral glue. The staining intensity (negative, 0; weak, 1; moderate, 2; strong, 3) and the percentage of positive signal (negative, 0; 1–30%, 1; 31–60%, 2; above 60%, 3) were used to assess the protein expression. The IHC score was obtained by multiplying the two scores mentioned above. The main antibodies we used were KCNQ1 (A2174, 1:400; ABclonal, Waltham, MA, USA), ITPK1 (10568-AP-1, 1:500: Proteintech, Rosemont, IL, USA), CD31 (11265-1-AP, 1:600: Proteintech, Rosemont, IL, USA), CD68 (28058-1-AP, 1:300: Proteintech, Rosemont, IL, USA).

Statistical analysis

R version 4.1.0 was the primary tool for all bioinformatics analyses and related statistical calculations. For biological experiment statistical analysis, GraphPad Prism version 8.1 was used. Regarding the multiple testing methods for the statistical analysis, we have used FDR control methods (the BH method). Regarding cross-validation of diagnostic models. We divided the fused dataset into training and testing datasets according to 1:1, and performed cross-validation and correction in the two datasets,

WGCNA and differential analysis to identify key diabetes-related gene markers

To reduce data analysis inaccuracies due to small sample sizes, data from three GEO platform datasets were merged, carefully excluding batch effects to create a combined dataset of diabetic patients and control samples. Initial analysis showed significant data dimension disparities between datasets (Fig. 2A). However, post-fusion, the data exhibited nearly uniform distribution patterns (Figs. 2B, 2C). Subsequent analysis of variance identified differential genes between diabetic patients and healthy controls (p < 0.05 and abs (log (FC)) >1), revealing 2,250 differential genes (Fig. 2D). Further, WGCNA analysis was conducted to isolate modular genes strongly linked to diabetes (Fig. 2E). In the WGCNA procedure, the soft threshold β was set at 9. For network visualization, 400 of these genes were randomly selected (Fig. S1). The WGCNA network was constructed using the partitioning method, merging and analyzing modules with high similarity, ultimately identifying 22 different gene modules (Fig. 2E). Correlation analysis between gene modules and clinical data of diabetic patients indicated that the yellow module had the strongest correlation with diabetes (r = 0.53, p = 5e−6) (Fig. 2F). To pinpoint hub genes within the yellow module, 368 genes with GS >0.2 and MM >0.6 were classified as central diabetes-related genes. The intersecting genes from differential analysis and WGCNA results yielded 32 diabetes-related biomarkers (Fig. 2G).

Diabetes related phenotype and characteristics in coronary artery disease

We employed the same methodology from the analysis of diabetes to boost sample size and lower error in the bioinformatics analysis of CAD patients. We reduced the batch effect by combining four datasets from the GEO platform. When the batch effect was excluded, the initial examination of the four datasets showed significant differences (Fig. 3A) that eventually became negligible (Fig. 3B). We also looked into the capacity of diabetes-related markers to classify patients with CAD, given the possible influence of diabetes on the development of CAD. We classified the gene expressions of 32 diabetes biomarkers in CAD patients using the consensus clustering method, which allowed us to identify diabetes-related phenotypes in CAD patients. Binary typing proved to be the most efficient method when comparing the delta areas and cumulative distribution function (PDF) curves of various categorization outcomes (Fig. S2). Consequently, CAD patients were divided into phenotype 1 and phenotype 2. Furthermore, we discovered that phenotype 1 had marginally higher expression of genes linked to diabetes than phenotype 2 (Fig. 3C). Additional examination of the microenvironmental scores revealed that phenotype 1 outperformed phenotype 2 in terms of stromal, immune, and estimation scores (Figs. 3D–3F).

Identification of microenvironment characteristics of diabetes related phenotype in CAD

We delved into the microenvironmental differences between diabetes-related phenotypes in CAD. We compared the phenotypes with regard to cytokines and inflammatory factors, which are important determinants of the immunological and inflammatory milieu. Typically, phenotype 1 had higher amounts of cytokines and inflammatory factors such CXCL4, CCL5, CXCR8, CXCR9, and other regulatory molecules, while phenotype 2 had lower levels of these factors (Fig. 4A). Furthermore, phenotype 1 also exhibited significant expression levels of numerous interleukins and associated receptors, such as IL2, IL3, PDGFD, etc., which are known to regulate inflammation (Fig. 4A). We also examined each phenotype’s immune cell composition. The estimation of immune cell infiltration was made using several algorithms. Phenotype 2 showed less immune cell infiltration than phenotype 1, which showed greater infiltration of immune cells such as CD8+ T cells, NK cells, etc. (Fig. 4B). Additionally, we examined the strength of endogenous and exogenous immunomodulation between the phenotypes. While endogenous immunity strength was similar between phenotypes, exogenous immunity was significantly stronger in phenotype 1 compared to phenotype 2 (Fig. 4C). Comparative analysis of classical immune features showed that promotive immune features were significantly more robust in phenotype 1 than in phenotype 2 (Fig. 4D).

Diabetes related phenotypes differ in myocardial lesion signaling and metabolic characteristics

Recognizing the significant distinctions between the two diabetes-related phenotypes, we looked at the differential pathways in the 2,180 differential genes that were found. Figure 5A shows a volcano diagram of these genes. We conducted Gene Set Variation Analysis (GSVA) for 10 changed pathways common to CAD and compared these between the two phenotypes, indicating considerable route variations (Fig. 5B). The subgroups’ different metabolic and myocardial injury-related pathways were revealed using enrichment analysis. Gene Ontology (GO) analysis identified biological process differences primarily in “fatty acid metabolic process,” “response to glucocorticoid,” “heart morphogenesis,” and other metabolism-related processes. Additionally, immune-related processes like “leukocyte cell–cell adhesion” and “regulation of innate immune response” were significantly enriched (Fig. 5C). KEGG analysis further enriched metabolic signals such as “Lipid and atherosclerosis” and “Carbon metabolism,” along with myocardial lesion signals like “Diabetic cardiomyopathy” and “Cardiac muscle contraction” (Fig. 5D). We investigated metabolic differences using the MSigDB database, taking into account the variation in metabolic signals among different phenotypes. The most important route variations are represented by a heatmap (Fig. 5E), which mainly shows the metabolic transport and absorption of proteins, fatty acids, glucose, and heart disease-related pathways. The significant distinctions between these symptoms, as revealed by the GSEA analysis, include alterations in metabolism, chemokines, and diabetic cardiomyopathy (Fig. 5F).

CAD diagnostic model constructed with diabetes-related genetic features shows excellent performance

To develop a more precise CAD diagnostic system, we further selected key genes for constructing a diabetes-related CAD typing using Lasso regression. A total of 32 core genes were included in the screening model (seed= 1, alpha=1, family=“gaussian”), and final model scores were derived using the predict algorithm. The final model incorporated 16 diabetes-related genes: “TRIM38,” “SAPG4,” “ODC1,” “NPY,” “MTDH,” “MT1G,” “LMF2,” “KCNQ1,” “ITPK1,” “GATAD2A,” “CRAT,” “CDK5RAP3,” “BAG3,” “ATP6V1B1,” “ARF4,” “AQP8.” (Figs. 6A, 6B). The modeling algorithm is as follows: model score= −0.0012369*exp (TRIM38) + 0.0017723*exp(SAPG4)+ 0.0014256*exp(ODC1)+ 0.0014256*exp(NPY)+ 0.0014256*exp(MTDH)+ 0.0016259*exp(MT1G) −0.0050854*exp(LMF2)+ 0.0015382*exp(KCNQ1)+ 0.0011388*exp(ITPK1)+ 0.0087858exp(GATAD2A)+ 0.0020088*exp(CRAT) −0.0065762*exp(CDK5RAP3) −0.00023941*exp(BAG3)+ 0.00817318exp(ATP6V1B1)+ 0.0065567*exp(ARF4)+ 0.0075323*exp(AQP8).

The predictive performance for each gene and the total LASSO diagnostic model was calculated. The areas under the receiver operating characteristic (ROC) curves were in Figs. 6C, 6D, 6E, respectively. Notably, when all 16 genes were included in a single diagnostic model, the AUC in the training dataset reached 0.823, and 0.758 in testing dataset (Fig. 6F, Fig. 6G).

Single-cell analysis explored CAD patient landscapes as well as expression patterns of modeled key genes

Further analysis was conducted on the expression landscape of diabetes-related genes using single-cell data from patients with coronary myocardial injury. After calculating and filtering the mitochondrial gene and ribosome ratios, as well as applying downscaling and tSNE clustering analyses (Fig. S3), we classified the patient cells into 11 subpopulations. These subpopulations were annotated using the “SingleR” package (Fig. 7A), identifying six cell types: macrophages, endothelial cells, monocytes, NK cells, muscle cells, and T cells (Fig. 7B). Macrophages and endothelial cells showed closer associations with other cells, indicating their potential roles in pathogenesis (Fig. 7C, Fig. S4). Dot plots were used to illustrate the top five genes differing among individual cells (Fig. 7D). The distribution of key markers in individual cells of the diabetes diagnostic model revealed that LMF2 and CDK5RAP3 were predominantly found in macrophages, while ARF4 was widely distributed among various immune and endothelial cells, as well as in muscle cells. KCNQ1, ATP6V1B1, MTDH, and ITPK1 were mainly distributed in macrophages (Fig. 7E, Fig. S5). Additionally, the heatmap exhibited the differential genes of various immune cells (Fig. S6).

KCNQ1 and ITPK1 were found to influence immune microenvironment of coronary artery disease

Building on the diagnostic model analysis, we further investigated the clinical relevance of these diagnostic markers. A mouse coronary artery disease model and a normal mouse control were established, with myocardial tissues removed, fixed, embedded, and subjected to IHC staining (Fig. 8A). We assessed the differential expression of the 16 gene diagnostic markers and some common immunological indicators. The results indicated upregulated expression of KCNQ1 and ITPK1 in mice with coronary artery disease, along with increased expression of CD31 and CD68. This suggests that KCNQ1 and ITPK1 are significant in the pathogenesis of CAD, potentially influencing disease progression by affecting the immune microenvironment in patients with CAD (Fig. 8B). We further compared mRNA levels in different subgroups and found the results to be consistent with those of histochemical staining (Fig. 8C). We further treated cells with CHD inducing conditions such as high glucose and fat, and then knocked down ITPK1 or KCNQ1. We found that high-glucose and high-fat conditions promoted the apoptosis rate of HUVEC, which was reversed by knockdown of KCNQ1 or ITPK1, suggesting the important role of ITPK1 or KCNQ1 in the development of CAD (Figs. 8D–8F).