The study cohort comprised participants diagnosed with OSA and/or CAD, and healthy controls. We consecutively screened and recruited patients admitted for coronary angiography and polysomnography (PSG) based on clinical indications. Participants diagnosed with CAD were selected based on evaluations by two cardiovascular specialists through CT coronary angiography or conventional coronary angiography. Exclusions included for those with severe infectious diseases, other cardiac conditions inducing angina (such as hypertrophic cardiomyopathy, myocarditis, or severe valvular heart disease), active malignancies or chemotherapy treatment, severe cardiac dysfunction, pacemaker implantation, frequent ectopic beats exceeding >10% of total heartbeats, prior revascularization procedures (coronary angioplasty or bypass grafting), major surgeries conducted within the past six months, or prior continuous positive airway pressure therapy. Patients with OSA were included if diagnosed via polysomnography (PSG) with an apnea-hypopnea index (AHI) ≥5. While excluding those who had been hospitalized for heart failure or respiratory disease in the past six weeks, were unable to undergo PSG, required oxygen therapy, had airflow-limiting conditions, or presented with renal insufficiency or long-term medication usage. Healthy volunteers without relevant diseases were recruited as controls. All subjects granted written informed consent, and the study was approved by the Institutional Review Board of the Second Affiliated Hospital of Fujian Medical University [NO. 597(2023), Quanzhou, Fujian, China]. Totally, 23 OSA, 25 CAD, and 24 CADOSA patients were enrolled, with health controls (n = 24) included for comparative analysis.

Clinical characteristics, including age, gender, body mass index (BMI), and relevant medical history, were documented as comprehensively as possible for each participant (Supplementary Table S1). Blood samples were taken and followed by PBMCs isolation using Ficoll-Paque density gradient centrifugation within two hours. The PBMCs were washed twice with phosphate-buffered saline and resuspended in Trizol reagent for RNA extraction, and then kept at −80℃ for RNA sequencing.

Total RNA was extracted utilizing the TRIzol Reagent (ThermoScientific, cat. #15596026) following the manufacturer's protocol. Subsequently, genomic DNA was removed by subjecting the samples to DNase I to eliminate DNA contamination. RNA quality and quantity were evaluated using the Agilent 2,100 Bioanalyzer and the NanoDrop™ OneC spectrophotometer, respectively. Samples with an RNA Integrity Number (RIN) >7.0 and a 260/280 ratio between 1.8 and 2.1 were preserved for subsequent analysis. RNA sequencing libraries were constructed and sequenced on the DNBSEQ-T7 platform (MGI), producing 150 bp paired-end reads. The raw sequencing data were processed through the following steps: (1) Adapter sequences and low-quality bases were trimmed using fastp (v0.23.2) with the following parameters: –dont_eval_duplication –qualified_quality_phred 20 –unqualified_percent_limit 8. (2) Clean reads were aligned to the human reference genome (GRCh38) using STAR (v2.7.6a, parameters: outFilterMatchNminOverLread 0.66 –outSAMmapqUnique 60 –chimOutJunctionFormat 1 –chimSegmentMin 10 –chimOutType SeparateSAMol). (3) Gene-level counts were quantified utilizing featureCounts (v1.5.1) with parameters of -d 30 -D 1,000 -C -s 0 -t exon -g fcid –primary -O.

Principal component analysis (PCA) was performed on Fragments Per Kilobase of exon model per Million mapped fragments (FPKM)- normalized gene expression data to reduce dimensionality and visualize sample clustering patterns based on overall transcriptomic variation. DEGs between each disease (OSA, CAD, CADOSA) group and the control group based on raw count data were identified using the edgeR and DESeq2 R package (v3.40.2). Initial quality filtering excluded genes detected in fewer than 25% of samples. Genes with an adjusted P-value (Benjamini-Hochberg correction) <0.05 and an absolute fold change >1.5 were considered significant. To ensure robust results, only genes consistently identified by both edgeR and DESeq2 analytical methods were retained as disease-associated DEGs. The Kyoto encyclopedia of genes and genomes (KEGG) enrichment of DEGs were performed utilizing the ClusterProfiler R package. The KEGG pathways with P < 0.05 were considered significantly enriched. The top 20 or significantly (P < 0.05) enriched functional pathways were visualized. Comparative analysis of DEGs from CADvsHC, OSAvsHC, and CADOSAvsHC through Venn diagramming enabled identification of condition-specific DEGs (unique to CAD or OSA) and comorbidity-associated DEGs (shared across all comparisons).

The Mfuzz (v2.58.0) R package was used to examine temporal gene expression patterns across disease progression during “HC→CAD/OSA→CADOSA” development based on FPKM- normalized gene expression values. The cluster number of Mfuzz were set as 6. Genes with statistically significant clusters (P < 0.05) retained for downstream analysis. Core genes for CAD and OSA were defined as the intersection between disease-specific DEGs and genes exhibiting inverted V-shaped expression patterns that peaked in the respective disease state. For CADOSA comorbidity, core genes were identified as those demonstrating both upward-trending expression patterns (maximal in CADOSA) and significant overlap with key comorbidity-associated DEGs.

To elucidate key molecular interactions, we established protein-protein interaction (PPI) networks for disease-specific DEGs or CADOSA-core DEGs using the STRING database. The resulting interaction data were imported into Cytoscape for topological analysis. Gene centrality was systematically evaluated through five distinct algorithms implemented in the CytoHubba plugin: Maximal Clique Centrality (MCC), Edge Percolated Component (EPC), Closeness, Degree, and Betweenness. Hub genes were rigorously defined as those consistently ranked among the top 10 most central nodes across all five-network metrics, ensuring the identification of biologically relevant and highly connected molecular regulators.

To identify optimal diagnostic biomarkers, we implemented a comprehensive machine learning pipeline. First, candidate genes were obtained from the intersection of core and hub genes. To develop predictive model of each disease, transcriptomic data was individually and randomly divided into training (70%) and validation (30%) sets. We then applied three distinct machine learning algorithms to prioritize robust biomarkers: (1) Elastic Net regression (with optimal regularization parameter λ = 0.05), (2) Random Forest, and (3) Support Vector Machine (SVM). Machine learning models were built using 10-fold cross-validation to improve generalizability. Final biomarker selection was determined by consensus between genes with non-zero coefficients form Elastic Net (indicating strong predictive contribution) and top-ranked genes based on variable importance (Mean Decrease Gini) from Random Forest. The performance of combined diagnostic model using all candidate and individual biomarker was systematically evaluated using ROC analysis.

The abundances of immune cells in four group samples were assessed individually using CIBERSORT algorithms. The difference in cell proportion was statistically evaluated. The associations between cells and biomarker genes were calculated and visualized using “ggcor” R package.

Gene expression of CADOSA biomarkers were verified using cDNA samples from 12 CADOSA and 12 HC cohorts. Quantitative PCR was performed in 10 μl reactions containing TB Green Premix Ex Taq II, gene-specific primers (MMP9, S100A12), and cDNA template, using the following cycling conditions: 95 °C for 30 s, then 40 cycles of 95°C for 5 s and 60 °C for 30 s Gene expression was calculated by the 2^(-ΔΔCt) method using β-actin as the endogenous control. The primers were synthesized by Sangon Biotech (Shanghai, China), and their sequences are: S100A12-forward primer 5′-GGGGCATTTTGACACCCTCT-3′ and S100A12-reverse primer 5′-CGCAATGGCTACCAGGGATA-3′; GLUD2-forward primer 5′-CAGCATCGTGGAGGACAAGT-3′ and GLUD2-reverse primer 5′-CGGTAGCCTTCGATGACCTC-3′; β-actin-forward primer 5′-GAGAAAATCTGGCACCACACC-3′, β-actin-reverse primer 5′-GGATAGCACAGCCTGGATAGCAA-3′.

Statistical analyses of clinical variables were performed to compare the demographic and clinical characteristics among the control, OSA, CAD, and CADOSA groups. All statistical analyses in this study were conducted using SPSS 27.0 software. The normality and homogeneity of variance for all factors were evaluated using the Kolmogorov–Smirnov test. Categorical variables, such as gender, were presented as frequencies (%) and analyzed using the Chi-square test. Continuous variables that follow a normal distribution are represented as mean ± standard deviation (SD), while non-normally distributed continuous variables are shown as median (interquartile range, IQR). One-way ANOVA was employed for intergroup comparisons of normally distributed continuous variables with homogeneous variance. Upon detecting significant differences were detected (P < 0.05), post-hoc pairwise comparisons were conducted using either LSD (for homogeneous variance) or Tamhane's T2 method (for heterogeneous variance). For data that is not normally distribute or exhibits non-homogeneous variance, the Kruskal–Wallis H-test was employed for intergroup comparisons, followed by Bonferroni-adjusted pairwise comparisons upon the identification of significant differences. All tests were two-tailed, with a P < 0.05 was considered statistically significant.

We compared demographic and clinical variables across four cohort groups (Table 1): healthy controls (HC, n = 24), CAD (n = 25), OSA (n = 23), and comorbid CADOSA (n = 24). No significant differences were observed in age, BMI, lipid profiles, and liver enzymes among the four groups. Marked disparities in gender distribution were noted (P = 0.002), with a higher male predominance in the OSA and CADOSA cohorts (76% and 91.7%, respectively). Statistical comparisons revealed that the heart rate in the OSA group was significantly higher than in the HC group (P < 0.05), which may be attributed to chronic nocturnal hypoxia-induced sympathetic overactivation. Regarding blood biochemical indices, the CADOSA group exhibited significantly higher levels of total triglycerides (TG) and creatinine (CR) compared to the CAD group (P < 0.05), suggesting that comorbid CAD and OSA exacerbates TG metabolic dysregulation and renal impairment. Among myocardial injury markers, the CAD group demonstrated significantly elevated D-dimer, NT-proBNP, and Hs-TnT levels compared to either the OSA or HC groups (P < 0.05), confirming increased thrombotic risk and myocardial damage in CAD patients. Echocardiographic parameters showed progressive increases in left ventricular end-diastolic diameter (LVD) and left ventricular end-systolic diameter (LVS) from HC→OSA→CAD→CADOSA, with the CADOSA group exhibiting significantly higher values than HC and/or OSA groups (P < 0.05), indicating that comorbid OSA aggravates cardiac systolic dysfunction in CAD patients. Polysomnography results demonstrated significantly higher AHI values and lower SpO2 levels in both OSA and CADOSA groups compared to HC and CAD groups (P < 0.05), reflecting substantially impaired respiratory function in OSA and comorbid CADOSA patients. Furthermore, we observed a progressive, though statistically non-significant, increase in the prevalence of carotid plaque, left ventricular hypertrophy, and valvular regurgitation across the spectrum from HC→OSA→CAD→CADOSA. This trend suggests that the synergistic effects of CAD-induced myocardial ischemia and OSA-mediated intermittent hypoxia may potentiate the development of these pathological manifestations, potentially exacerbating disease severity of CAD and OSA.

The RNA-seq yielded sequencing data size ranging from 5.44 Gbp to 7.05 Gbp for each sample across four groups, resulting in 673.74 Gbp data in total (Supplementary Table S2). After quality filtration, >90% of reads were preserved as clean data for subsequent reads mapping, gene quantification and differential expression analysis. The detailed workflow of this study was shown as Figure 1.

Principal component analysis based on transcriptomic gene expression showed that CAD and HC samples were apparently distinguished, indicating an intrinsic difference in gene regulation between the two groups (Figure 2A). The limma and DESeq2 methods were combined to identify commonly and DEGs in the CAD. As a result, the CAD group exhibited 571 up-regulated and 261 down-regulated DEGs (Supplementary Table S3), with over two-thirds being up-regulated and nearly one-third being down-regulated (Figures 2B,C). These 832 DEGs were significantly enriched in 30 KEGG pathways (level 3 in Figure 2D), which were classified into 15 level 2 categories and five level 1 pathways. Notably, “immune system” and “signaling molecules and interaction” were the top two abundantly enriched categories. The immune-related pathways, including “Cytokine–cytokine receptor interaction”, “B cell receptor signaling pathway,” “Complement and coagulation cascades”, and “Neutrophil extracellular trap formation”, exhibited prominent enrichment, suggesting a robust involvement of immune regulation in the CAD transcriptomic changes. The signaling molecules and interaction pathways such as “Neuroactive ligand–receptor interaction”, “cytokine-cytokine receptor interaction”, and “cell adhesion molecules” were significantly enriched, highlighting the critical role of intercellular communication in CAD development. Metabolic reprogramming was another key feature of the CAD, with DEGs strongly associated with “Tryptophan metabolism”, “Histidine metabolism”, “Glutathione metabolism”, “Ether lipid metabolism”, and “Arachidonic acid metabolism”. These findings underscore the importance of amino acid and lipid metabolism in the biological regulation of the CAD. Furthermore, pathways related to “phagosome” and “efferocytosis” were enriched, indicating heightened immune clearance activity and potential involvement in inflammation, tissue homeostasis. In the context of human diseases, DEGs were enriched in “Hepatocellular carcinoma”, “Basal cell carcinoma”, “Tuberculosis”, and “Malaria”, suggesting a possible link between the identified genes and infection or oncogenic processes. Collectively, these results demonstrate that the DEGs are not only pivotal in immune and signaling pathways but also play a central role in metabolic adaptation, particularly in amino acid and lipid metabolism, which may have broader implications for the CAD mechanisms and therapeutic targeting.

Samples from the OSA and HC groups exhibited clear separation by PCA analysis, confirming distinct transcriptomic profiles (Figure 3A), with limma and DESeq2 identifying 87 up-regulated and 79 down-regulated DEGs (Figures 3B,C; Supplementary Table S4). KEGG enrichment analysis revealed these DEGs were significantly associated with cytoskeleton in muscle cells (level 3) of cell motility (level 2), PI3K-Akt signaling of signal transduction, arrhythmogenic right ventricular cardiomyopathy of cardiovascular disease, nitrogen metabolism of energy metabolism, and cholinergic synapse of nervous system (Figure 2D). While cancer-related pathways like proteoglycans in cancer were enriched, their relevance likely relates to shared oxidative stress mechanisms rather than direct oncogenesis. These multi-systemic pathway alterations collectively demonstrate how chronic intermittent hypoxia in OSA drives molecular changes affecting cardiovascular integrity, neuromuscular function, and metabolic homeostasis, providing a framework for understanding OSA-related complications.

Similarly, the CADOSA samples were clustered and separated from HC samples (Figure 4A). Compared to healthy controls, CADOSA exhibited 269 up-regulated and 107 down-regulated DEGs (Figures 4B,C, Supplementary Table S5), showing more DEGs than OSA alone but fewer than CAD alone. KEGG enrichment analysis revealed these DEGs were significantly associated with human diseases, organismal systems and cellular processes. The dysregulation of immune-related pathways included “NOD-like receptor signaling pathway”, and “neutrophil extracellular trap formation” suggests chronic inflammation. Other consisted of the “PPAR signaling pathway” in endocrine system and “protein digestion and absorption” in digestive system might indicate body function disturbances. Similar to the CAD, the comorbidity also exhibited the same enriched pathways of cell motility (cytoskeleton in muscle cells), and transport and catabolism (phagosome and efferocytosis). The infection-related pathways (such as “Staphylococcus aureus infection” and “malaria”) and primary immunodeficiency in human diseases reflects immune vulnerability. These multi-systemic alterations highlight the synergistic effects of CAD and OSA, where chronic intermittent hypoxia might exacerbate both immune dysregulation and cardiovascular damage.

To investigate the distinct molecular profiles among CAD, OSA, and CADOSA, we identified unique and shared DEGs across the three comparisons. The Venn diagram revealed 555 CAD-specific DEGs and 44 OSA-specific DEGs, and 73 CADOSA-core DEGs which were commonly shared in three comparisons (CAD vs. HC, OSA vs. HC, and CADOSA vs. HC) (Figure 5A). Pathway enrichment revealed distinct profiles of these disease-specific and comorbidity-core DEGs, with CAD exhibited the most complex molecular regulations (Figure 5B). The CAD-specific DEGs were predominantly associated with immune dysregulation (complement and coagulation cascades, hematopoietic cell lineage, and B cell receptor signaling pathway) and metabolic disturbances (glutathione/tryptophan/histidine metabolism). The OSA-specific DEGs showed enrichment in arginine biosynthesis, while the shared DEGs in CADOSA converged on cholinergic synapse.

To further elucidate the molecular interactions, we constructed protein-protein interaction (PPI) networks using STRING and identified hub genes for each disease group. Hub gene selection was performed using five topological algorithms (MCC, EPC, Closeness, Degree, and Betweenness) implemented in the CytoHubba plugin for Cytoscape, and defined as with the intersection of top-ranked 10 genes across all five methods. Analysis results revealed distinct network characteristics of three diseases. The 555 CAD-specific DEGs formed the most extensive and complex PPI network, with three key hub genes identified, including CD68, ICAM1, and IL10 (Figures 5C,D). The 44 OSA-specific DEGs generated a relatively simple interaction network yielding eight hub genes, consisting of HBD, ALAS2, NIP1L4, CPA3, HDC, GLUD2, MYZAP, and PKP2 (Figures 5E,F). While the 73 shared CADOSA DEGs produced an intermediate complexity network with six central hub genes, including HP, TIMD4, S100A12, PLBD1, MMP9, and VCAN (Figures 5G,H).

Through time-series analysis of DEGs, we characterized dynamic gene expression patterns during disease progression from healthy controls (HC) to CAD/OSA and subsequently to CADOSA comorbidity. In the “HC→CAD→CADOSA” transition, we identified four expression clusters (Supplementary Figure S1): an upward-trending cluster (81 genes) peaking in CADOSA, a downward-trending cluster (180 genes), an inverted V-shaped cluster (580 genes) maximal in CAD, and a V-shaped cluster (120 genes). Similar analysis of “HC→OSA→CADOSA” revealed four patterns (Supplementary Figure S2): an upward-trending cluster (260 genes) highest in CADOSA, a V-shaped cluster (73 genes), a downward cluster (84 genes), and an inverted V-shaped cluster (40 genes) peaking in OSA. Integration of these patterns identified 374 CAD-specific core genes (intersecting unique DEGs with CAD-peaking inverted V-shaped patterns), 15 OSA-specific core genes showing similar characteristics, and 8 CADOSA core genes exhibiting upward trends in both progression pathways.

Through integration of core genes with hub genes, we identified candidate biomarkers including three CAD-associated (CD68, ICAM1, IL10), two OSA-associated (GLUD2, PKP2), and two CADOSA-associated (MMP9, S100A12) genes. Following rigorous evaluation using three machine learning approaches, we established the final biomarker panels and developed diagnostic models for each condition. The composite diagnostic model demonstrated excellent discriminatory capacity across all conditions. For CAD, the elasticNet-based model achieved an AUC of 0.92, while individual biomarkers CD68, ICAM1, and IL10 showed robust standalone performance with AUCs of 0.94, 0.84, and 0.75, respectively (Figure 6A). CD68 had the highest expression, while IL10 displayed the lowest expression. Similarly, the random forest-derived OSA model exhibited strong diagnostic accuracy (AUC = 0.93) followed by SVM (AUC = 0.91), with GLUD2 maintaining high predictive value as a single biomarker (AUC = 0.82) (Figure 6B). For CADOSA comorbidity, the random forest model incorporating both biomarkers performed optimally with an AUC of 0.89, while individual analysis revealed significant diagnostic capability for S100A12 (AUC = 0.83) and MMP9 (AUC = 0.78) (Figure 6C), confirming their clinical utility as standalone indicators. The S100A12 showed much higher expression level than the MMP9. Subsequent qRT-PCR validation confirmed their significant upregulation in CADOSA patients (Figure 6D), further supporting their potential clinical utility as independent diagnostic indicators.

Functional enrichment analysis of DEGs highlighted the involvement of immune dysfunctionin the pathogenesis of CAD, OSA, and CADOSA. To further characterize immune dysregulation, we performed immune cell infiltration analysis in PBMCs across the three disease groups. We identified 22 distinct immune cell subtypes, with monocytes being the most abundant, followed by CD4+ memory resting T cells, NK resting cells, and CD8+ T cells, collectively accounting for nearly 80% of PBMCs (Figure 7A). Compared to healthy controls, all three patient groups exhibited a significant reduction in CD4+ memory resting T cells, CD8+ T cells, B memory cells, and B naïve cells, alongside a pronounced increase in monocytes (Figure 7B).

To assess immune cell interactions, we conducted correlation analyses in CAD samples. CD4+ memory resting T cells showed a negative association with regulatory T cells and resting mast cells, while NK resting cells were inversely correlated with macrophages M0 (Figure 7C). Notably, the three most abundant immune subsets—CD4+ memory resting T cells, NK resting cells, and CD8+ T cells—all demonstrated a negative relationship with monocytes. Furthermore, we examined the association between CAD biomarker genes and immune cells. CD68 expression was positively correlated with monocytes and NK resting cells but negatively correlated with CD8+ T cells, whereas IL10 exhibited a positive relationship with monocytes.

In OSA patients, CD4+ memory resting T cells were inversely associated with activated CD4+ memory T cells and macrophages M0. NK resting cells showed a negative correlation with macrophages M0 and monocytes, while CD8+ T cells were positively linked to activated CD4+ memory T cells and macrophages M0 but negatively correlated with monocytes. However, no significant relationship was observed between the OSA biomarker GLUD2 and immune cell proportions (Figure 7D).

For the CADOSA cohort, we found that CD8+ T cell proportions were negatively associated with CD4+ memory resting T cells, and monocytes exhibited an inverse correlation with NK resting cells (Figure 7E). Additionally, we evaluated the relationship between two key CADOSA biomarkers and immune cell infiltration. S100A12 displayed strong positive correlations with monocytes but negative associations with NK resting cells (Figure 7C). Conversely, MMP9 was positively correlated with macrophages M0, CD4+ memory resting T cells, and B memory cells (Figure 7E).