Machine Learning and Bioinformatics Unravel Gene Signatures of Coronary Artery Disease Comorbidity with Periodontitis

2.1 Raw data

Gene expression data were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) using the accession numbers GSE10334 [12] and GSE66360 [13]. The GSE10334 dataset includes 247 samples derived from 90 patients. These samples encompass 183 tissues affected by disease and 64 healthy control tissues from the same patients. The GSE66360 dataset consists of isolated circulating endothelial cells from 49 patients diagnosed with acute myocardial infarction and 50 healthy individuals.

2.2 WGCNA network construction and module identification

The WGCNA package in R [14] was employed to construct co-expression networks and identify co-expressed gene modules. This process involved several steps. First, samples underwent hierarchical clustering to identify outliers, which were subsequently removed alongside any unqualified genes using the goodSamplesGenes function. Second, the soft-thresholding power (β) was calculated using the R function pickSoftThreshold and converted into a topological overlap matrix (TOM). Third, based on the TOM, the network’s interconnectivity was established, and the dynamic tree-cutting algorithm was applied to identify co-expressed gene modules. Fourth, modules are associated with clinical features by calculating gene saliency and module affiliation. The genes involved in the corresponding modules were used for subsequent analysis. Finally, the feature gene network was visualized.

2.3 Identification of common DEGs

GEO2R (www.ncbi.nlm.nih.gov/geo/ge2r) is an online deg analysis tool based on the Limma software package [15]. An adjusted P-value < 0.01 and |logFC| ≥ 1 were defined as thresholds for screening DEGs. We extracted the genes most associated with periodontitis and coronary artery disease from the WGCNA module and cross-analyzed the two gene lists for common DEGs.

2.4 Enrichment analyses of common DEGs

To further understand which biological functions these common DEGs are involved in, we performed gene enrichment analysis using the KOBAS 3.0 database, including gene ontology and Kyoto Encyclopedia of Genes and Genomes pathways [16]. Adjusted P-values < 0.05 was considered significant.

2.5 Machine learning algorithms for identifying feature genes

Machine learning models are constructed using the filtered DEGs. We used four machine learning methods to construct the model respectively:

1. Random Forest (RF): This versatile ensemble method combines multiple decision trees, each built independently from random subsets of the data [17, 18]. Individual trees learn and predict independently, with the final prediction being the average of all tree predictions.

2. Support Vector Machine Recursive Feature Elimination (SVM-RFE): This discriminative classifier iteratively removes the least informative features to identify the most relevant ones for classification [19]. The model is trained using labeled samples, and test samples are classified based on the optimal separation hyperplane.

3. XGBoost: This powerful ensemble method leverages gradient boosting to create a robust model by combining multiple weak decision trees [20]. XGBoost facilitates a meticulous comparison of classification errors and model complexity.

4. Generalized Linear Model (GLM): This flexible method extends the linear model by incorporating a link function to model the relationship between the response variable and the linear combination of predictor variables [21]. We compared the performance of these four methods to identify the one with the highest accuracy for screening disease signature genes. The “DALEX” package was used to interpret the models, allowing us to visualize the residual distributions and feature importance for each method. Performance was further evaluated by visualizing the receiver operating characteristic curves using the “qROC” package. A larger area under the curve indicates higher model accuracy.

2.6 Receiver operating characteristics curve and nomogram construction

To assess the diagnostic potential of identified characteristic genes, their expression levels in coronary artery disease and periodontitis groups were compared. The area under the receiver operating characteristic curve was calculated for each gene, along with its 95% confidence interval, to estimate diagnostic accuracy. To mitigate bias, datasets GSE6751 and GSE71226 were employed as validation cohorts. Only genes exhibiting an AUC > 0.7 in both the original dataset and both validation sets were included for further analysis. Subsequently, nomograms were constructed using the RMS R package [22] to facilitate clinical decision-making.

2.7 Immune Infiltration Analyses

CIBERSORT is a bioinformatics tool that estimates the relative proportions of various cell types within complex tissue samples [23]. We utilized CIBERSORT to evaluate immune cell infiltration in the coronary artery disease and periodontitis groups. This algorithm offers a comprehensive assessment of immune cell abundance by profiling 24 distinct immune cell types. The proportion of each immune cell type in each sample was visualized using histograms. Additionally, heatmaps depicting the correlations between different immune cell types within the CAD and PD groups were constructed using the “corrplot” R package [24]. We further employed t-test to compare immune cell infiltration levels between disease and control groups. Finally, Spearman correlation analysis was performed to explore potential associations between immune cell populations and signature genes.

2.8 Clinical validation of RNA-seq analysis

Blood samples and clinical data of 5 patients with both PD and CAD, and another 5 patients with only CAD, were collected at Shiyan People’s Hospital from May 2024 to July 2024. In order to obtain a homogeneous sample with patients enrolled in a specific age group (35-75 years) and equal representation of men and women, stratified random sampling was performed by dividing patients who met the study criteria into three age strata of 35-50, 51-65, and 66-75 years to reduce possible bias related to gender or different ages. The patients included in the study were categorized into CAD group, PD+CAD group.

The exclusion criteria were (1) those who had taken antibiotics or hormonal drugs within six months; (2) women who were in pregnancy or breastfeeding; (3) systemic diseases and history of drug allergy; (4) those who had a history of treatment for periodontal diseases within one year; (5) those who suffered from renal and thyroid diseases; and (6) those who had used non-steroidal anti-inflammatory drugs within three months. The study experimental design and workflow were approved by the Science and Technology Ethics Review Committee of Hubei Medical College (reference number: 2024-RE-26). Written consent was obtained from patients before inclusion in the study.

All patients included in the study had 5 ml of venous blood collected on admission, placed in an anticoagulation tube containing 3.2% sodium citrate, left to stand for 5 min at room temperature, centrifuged at 3500 r/min, the precipitate extracted, 1 ml of Trizol was added and mixed, and the peripheral blood samples was immediately frozen and stored on dry ice. RNA extraction of the samples was performed using Tianmo#TR205-200 kit. Extracted total RNA was tested for RNA integrity using an Agilent Bioanalyzer 2100 (Agilent technologies, Santa Clara, CA, US) test with a Qubit® 3.0 Fluorometer (Life Technologies, CA, USA) and a Qubit® 3.0 Fluorometer (Life Technologies, CA, USA) and Nanodrop One spectrophotometer (Thermo Fisher Scientific Inc, USA) were used to determine the concentration and purity of total RNA.

Total RNA was extracted for microarray analysis from venous blood of 10 patients who met the study inclusion criteria. The microarray study was performed using Affymetrix GeneChip®Human Genome U133 Plus 2.0 Array for standard operating procedures and quality control of mRNA. Samples from different groups were randomized to prevent batch effects and sample processing was performed simultaneously to minimize variation. Quality control of microarray data was performed using the Transcriptome Analysis Console (TAC). For each mRNA expression probe set, it was then normalized using the TMM (trimmed mean of M values) algorithm (https://www.ncbi.nlm.nih.gov/COG/). The median polish after quantile normalization was used to remove probe affinity effects when calculating probe set summaries. Expression values were log2 of signal intensity. low expression probe sets were removed if they were not expressed in more than 25% of the samples. The expression level of each gene was estimated by calculating the median expression level of all probes within the coding region of that gene. Residual normalized gene expression was fitted to a regression model using the R package “limma”. A mitigated t-test was used to determine the deg between the CAD and PD+CAD patient groups. Genes with P-value < 0.05 and at least a 1-fold change (FC) were determined to be significant DEGs.

Total RNA sequencing libraries were prepared using an Illumina TruSeq Stranded Total RNA Library Prep kit to deplet cytoplasmic rRNA according to the manufacturer’s instructions. All cDNA libraries were QC using the Agilent DNA1000 Kit (Agilent) prior to next-generation sequencing. samples were sequenced on an Illumina NovaSeq 6000 to a depth of approximately 50 million pairs of end-reads per biological sample. FastQ was utilized for quality control of raw reads. Pairwise end-reads were aligned to the human reference genome GRCh38 using the STAR comparator with default parameters. Gene counts were calculated using featureCounts [25]. Gene counts were normalized for “limma” default parameters. And differential gene expression (DE) analysis was performed using “limma” package. The significant deg called with P-value < 0.05 and |logFC| >1. Graphs were generated using R (version 4.3.1).

3.1 WGCNA network construction and module identification

WGCNA was employed to identify co-expressed gene modules potentially associated with disease. Before module identification, samples were clustered to detect and remove outliers. No significant outliers were identified in the GSE10334 (Fig. 2A) or GSE66360 (Fig. 2C) datasets. Within the WGCNA framework, an optimal soft-threshold power (β) of 11 was determined for the GSE10334 dataset (Fig. 2B), while a β of 8 was optimal for the GSE66360 dataset (Fig. 2D). This analysis resulted in the identification of nine modules in the GSE10334 dataset and seven modules in the GSE66360 dataset. Correlation analysis revealed a strong negative correlation between the brown module and periodontitis (r = −0.55) (Fig. 2E). Similarly, the grey module exhibited a strong negative correlation with coronary artery disease (r = −0.65), while the brown was positively correlated with coronary artery disease (r = 0.61) (Fig. 2G). Moreover, the correlation of module membership in brown (r = 0.83) and gene significance for PD samples was observed, as well as grey module membership (r = 0.66) (Figure. 2F, H).

• Download figure
• Open in new tab

Figure 1.

Research design flow chart. In this study, we used WCGNA to explore the correlation between modular signature genes and the occurrence of periodontitis and coronary artery disease, while DEGs for periodontitis and coronary artery disease were obtained using GEO2R. The overlapping genes between modular genes and DEGs were used for further analyses. Subsequently, the optimal machine learning model was selected to identify the signature genes of the two diseases and validate the accuracy of the model, and the expression of the signature genes was verified in the validation dataset. Immune infiltration analysis and gene-immune cell correlation analysis were then performed to demonstrate that they were significantly associated with immune cell dysregulation. Finally, RNA-sequencing of clinical samples was performed to validate the identified genes and explore their potential as biomarkers for early diagnosis and prevention of co-morbid periodontitis and CAD.

• Download figure
• Open in new tab

Figure 2.

(A, C) Clustering dendrogram of samples based on their Euclidean distance in GSE10334 and GSE66360. (B, D) Determination of soft-thresholding power for GSE10334 and GSE66360. (E, G) Heatmap of the correlation between module genes and the occurrence of periodontitis and coronary artery disease. Modules with different colors represent different gene modules, and the numbers in the modules represent the correlation between the module and the phenotype. (F, H) The scatter plots of module membership and gene significance for PD and CAD respectively.

3.2 Identification of common DEGs

Differential expression analysis using GEO2R identified 517 and 687 DEGs in the GSE10334 and GSE66360 datasets, respectively (Figure. 3A, B). We then focused on genes within WGCNA modules that exhibited strong associations with periodontitis and coronary artery disease. Genes from these disease-associated modules were intersected with the GEO2R-identified DEGs to identify commonly dysregulated genes in both diseases. Utilizing Venn diagram analysis, this approach revealed 48 DEGs shared between periodontitis and CAD (Figure. 3C).

• Download figure
• Open in new tab

Figure 3.

(A) The volcano map of GSE10334. (B) The volcano map of GSE66360. Red dots represent up-regulated genes, blue dots represent down-regulated genes and grey dots represent genes with no significant difference. (C) Venn diagram shows that 48 common DEGs in GSE10334 and GSE66360. (D) Enrichment result of common DEGs GO term; (E) enrichment result of common DEGs KEGG pathway. Adjusted P-value < 0.05 was considered significant. The ordinate represents the enriched term, and the abscissa represents the proportion of genes involved in the term. The size of the dots represents the number of genes, and the color of the dots represents the P-value.

3.3 Analysis of the functional characteristics

Gene Ontology enrichment analysis revealed that commonly dysregulated genes were significantly enriched in biological processes including response to lipopolysaccharide, response to bacterial-derived molecules, bone marrow leukocyte migration, and neutrophil migration (Fig. 3D). Cellular component enrichment highlighted secretory granule membranes, secretory granule lumen, and plasma membrane. Furthermore, enriched molecular functions included CXCR chemokine receptor binding, pattern recognition receptor activity, and phosphoprotein binding. KEGG pathway analysis identified significant enrichment of DEGs in pathways associated with lipids and atherosclerosis, the NF-κB signaling pathway with IL-17, and the TNF signaling pathway (Fig. 3E).

3.4 Candidate signature genes of CAD and PD using machine learning

We utilized the identified intersection genes to train four machine-learning models. The “DALEX” package in R was employed to visualize the residual distributions of each model within the training set. Subsequently, feature importance analysis was performed to rank the top ten genes in each model based on their Root Mean Squared Error (RMSE). Finally, the diagnostic performance of all four models in the training set was evaluated using ROC curves generated through 5-fold cross-validation.

Analysis of the residual distributions (Fig.4A, B) revealed that the Support Vector Machine and XGBoost models exhibited lower sample residuals, indicating better model fitting. Furthermore, the XGB model demonstrated the largest area under the ROC curve (AUC) in Fig. 4D, signifying superior prediction accuracy. Consequently, the top five most important genes (FOS, VNN2, PECAM1, MME, RGS1) identified by the XGB model were selected as predictor variables for further analysis (Fig. 4C). To assess generalizability, the XGB model was validated on independent datasets (GSE106090, GSE179789) from the validation set. The results depicted in Fig. 4E and Fig. 4F confirm that the XGB model retains good predictive ability in these external datasets.

• Download figure
• Open in new tab

Figure 4.

Construction and evaluation of SVM, RF, XGB, and GLM machine learning models. (A) Cumulative residual distributions of samples computed by SVM, RF, XGB, and GLM machine learning models respectively. (B) Box line plots of sample residuals based on SVM, RF, XGB, and GLM machine learning models, with the red dots representing the root-mean-square (RMS) of the residuals. (C) The importance scores of the top ten genes obtained by SVM, RF, XGB, and GLM. (D) ROC analyses of the four machine learning models based on the 5-fold cross-validation in the training set. (E, F) ROC analysis of the XGB model in the 2 validation sets.

3.5 Validation of biomarker diagnostic value and nomogram construction

In the external validation set (GSE106090), all five identified signature genes exhibited significantly higher expression levels in periodontitis patients compared to healthy controls (Fig. 5A-E), suggesting their potential involvement in periodontitis pathogenesis. The AUC values of the receiver operating characteristic (ROC) curves for these genes ranged from 0.704 (FOS) to 0.891 (PECAM1) (Fig. 5F-J), all exceeding 0.7, which indicates good diagnostic potential.

• Download figure
• Open in new tab

Figure 5.

Expression of feature genes in GSE106090. (A-E) Expression of signature genes in periodontitis patients and healthy population. (F-J) ROC showing diagnostic performance of signature genes.

We further evaluated the diagnostic efficacy of these genes for coronary artery disease prediction in the validation set (GSE179789). Consistent with the periodontitis findings, all genes displayed higher expression in CAD patients compared to controls (Fig. 6A-E). The AUC values of the ROC curves for these genes ranged from 0.710 (PECAM1) to 0.861 (RGS1) (Fig. 6F-J), again exceeding 0.7, suggesting good diagnostic performance for CAD as well. Collectively, these results indicate that the identified signature genes hold promise as diagnostic biomarkers for both periodontitis and coronary artery disease.

• Download figure
• Open in new tab

Figure 6.

Performance of signature genes in GSE179789. (A-E) Expression of signature genes in patients with coronary artery disease and healthy population. (F-J) ROC showing diagnostic performance of signature genes.

To gain insights into the XGB model’s predictions, we employed a nomogram (Fig. 7A). This graph depicts the estimated risk of developing coronary artery disease compared to periodontitis for a given individual. The expression level of each signature gene translates to a corresponding score, and the sum of these five scores generates a composite score. This composite score serves as an indicator of disease prevalence for a specific sample. The calibration curve (Fig. 7B) demonstrates good agreement between the actual risk of developing the disease and the model’s predicted risk, highlighting minimal error. Additionally, decision curves (Fig. 7C) support the high accuracy of the nomogram for risk prediction.

• Download figure
• Open in new tab

Figure 7.

Nomogram construction and diagnostic efficacy validation. (A) A nomogram was constructed based on five selected feature genes, each corresponding to a score, and the total score of the five feature genes was used to predict the risk of periodontitis and coronary artery disease. (B-C) Calibration curves (B) and decision curves (C) were used to assess the predictive efficiency of the XGB model.

3.6 Immune cell infiltration analysis

Evaluation of the immune cell composition in periodontitis and coronary artery disease showed significant differences in immune cell profiles between the diseased and control groups (Fig. 8A, C). In patients with periodontitis and coronary artery disease, eosinophils, activated memory CD4+ T cells, and resting dendritic cells were highly positively correlated with B cells (Fig. 8B, D).

• Download figure
• Open in new tab

Figure 8.

(A, C) Stacked bar chart of the immune cell. The different colors of the rectangular bars in the diagram represent different immune cells, and the length represents the proportion of immune cells. (B, D) The correlation matrix of immune cell proportions. The numbers in the squares represent the correlation coefficients between the corresponding immune cells. (E) Comparison of immune cell scores in patients with periodontitis and a healthy population. (F) Comparison of immune cell fractions in patients with coronary artery disease and a healthy population. The horizontal axis indicates the different immune cells and the vertical axis indicates the proportion of immune cells. T-test was used to compare the disease group with the control group. (G) Correlation between characterized genes and immune cell fractions in patients with periodontitis. (H) Correlation between characteristic genes and immune cell components in patients with coronary artery disease.

In patients with periodontitis and coronary artery disease, there was a significant increase in monocytes, and neutrophils and a significant decrease in activated memory CD4+ T cells (Fig. 8E, F). We next investigated the relationship between the identified signature genes and immune cell components. In periodontitis, VNN2, PECAM1, MME, and RGS1 displayed strong negative correlations with CD8+ T cells, while VNN2 and MME showed strong positive correlations with neutrophils (Fig. 8G). Similarly, in CAD, VNN2, PECAM1, MME, and FOS exhibited strong positive correlations with neutrophils, whereas RGS1 showed a negative correlation (Fig. 8H). These findings suggest that the signature genes may influence immune responses in both periodontitis and CAD.

3.7 Clinical validation of RNA-seq analysis

Differential gene expression analyses were performed using “limma” packages. Of the five characterized genes we screened using the machine learning algorithm, four of them were expressed differently in two groups. Specially, the FOS gene was confirmed to be highly expressed (logFC > 1.5, P-value = 0.01) in the PC group compared to the C group (Fig. 9). Overall, the expression pattern of the FOS gene was highly consistent in RNA sequencing of clinical samples as well as in the analyses of the present study, confirming the potential of the FOS gene as a signature gene for coronary heart disease and periodontitis.

• Download figure
• Open in new tab

Figure 9.

Peripheral blood RNA sequencing analysis of the coronary artery disease group (C group) and the periodontitis combined with coronary artery disease group (PC group). The box plot indicated that FOS, MME, PECAM1 and VNN2 genes were up-regulated in the PC group compared with C group (|logFC| > 1, P-value < 0.05).