Differential expression and significance of peripheral blood genes in coronary artery heart disease.

Background: The peripheral blood gene expression profile of patients with coronary artery disease (CAD) has not been fully resolved. The aim of this study was to further analyze the peripheral blood transcriptome information of CAD patients and to uncover key genes and regulatory mechanisms in the pathogenesis and disease progression of CAD.
Methods: The Gene Expression Omnibus (GEO) database was applied to screen out differentially expressed genes (DEGs) in the peripheral blood of CAD patients, and the DEGs were subjected to Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and gene set enrichment analysis (GSEA). The core genes were screened by GO, KEGG, and GSEA, and the gene-gene interaction (GGI) and protein-protein interaction (PPI) networks of DEGs were constructed. The GeneCards database was used to obtain CAD-related genes, and the GEO dataset was used to obtain intersecting genes. The intersecting genes were analyzed for bioenrichment and prediction of potential therapeutic agents, and predictive models were constructed for the intersecting genes. Finally, immune infiltrating cells from the GEO dataset were analyzed.
Results: A total of 79 DEGs were screened in the peripheral blood of CAD patients, of which three were autophagy-related genes. Biological enrichment analysis showed that the DEGs were associated with metabolic pathways, and vascular smooth muscle contraction and were mainly involved the MAPK signaling pathway, metabolic pathways, and the PI3K-Akt signaling pathway. The S100A8, ENTPD1, and MMP9 further screened were screened. A total of 11 CAD crossover genes and 75 potential therapeutic agents were obtained, and the column line graph prediction models constructed for S100A8, HSPB1, F5, MMP9, and PDE9A had good predictive power. There were significant differences in immune cells in CAD patients compared to healthy individuals, especially in T cells regulatory (Tregs) and B cells naïve. Conclusions: The peripheral blood of CAD patients screened by the GEO dataset was significantly different from that of the healthy population, and the DEGs and intersecting genes were involved in numerous key biological processes that may be involved in the development and progression of CAD and could serve as its regulatory sites and therapeutic drug targets. 2022 Journal of Thoracic Disease. All rights reserved.

Introduction

The development of coronary artery disease (CAD) is controlled by multiple risk factors and is one of the most common causes of death (1). Clinical observations as early as the 1950s support the idea that CAD risk is heritable (2), and environmental and genetic factors such as hypertension, obesity, dyslipidemia, and family history may contribute to the development of CAD (3-5). According to the latest global burden of disease report, ischemic heart disease is the tenth highest cause of disease in the middle-aged and elderly population (6), and CAD is the leading cause of death due to cardiovascular disease, accounting for approximately 45% of all cases (7). Currently, percutaneous transluminal coronary angiography, which is the standard for the diagnosis of CAD, is an invasive test that may cause serious adverse effects and, to a certain extent, limit its use in clinical practice. However, at this stage, there is a lack of effective biomarkers to predict the risk of CAD development in clinical practice.

The literature shows that as early as 1967, the Framingham cohort study revealed the association of many risk factors with CAD, established a cardiovascular disease risk assessment model, and calibrated the model to be applicable to populations around the world (8,9). This model is more accurate in predicting short-term disease in middle-aged and older adults, but is less effective in younger patients, and does not predict long-term disease risk. Based on Framingham’s model, subsequent researchers have made improvements, but it still has some shortcomings (10,11). The above models are based on traditional CAD risk factors such as gender, age, blood pressure, and lipids, and the prediction results are more reliable only when the risk factors accumulate to a certain level, thus limiting their generalizability, predictive value, and capacity to predict early-onset CAD.

The transcriptome is a collection of all RNAs produced by a species or a specific cell type. Based on high-throughput analysis and detection technology, thousands of targets and pathways can be screened to obtain common gene variant loci for CAD, which can reveal the differences in gene expression and structure of CAD, elucidate the molecular mechanism, predict and intervene in CAD, and potentially prevent serious adverse cardiovascular events. Seven common disease genes, including CAD, have been identified in previous studies, and a database of related genetic information has been established (12). In recent years, a genome-wide association study has identified more than 160 susceptibility loci for CAD (13), which has greatly facilitated the process of CAD genetic research. The aim of this study was to analyze the RNA expression profiles of CAD patients in anticipation of discovering new biomarkers as predictors of CAD pathogenesis for early diagnosis and specific targeted therapy for CAD. We present the following article in accordance with the STREGA reporting checklist (available at https://jtd.amegroups.com/article/view/10.21037/jtd-22-991/rc).

Methods

Dataset acquisition

The GSE20680 and GSE20681 datasets were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Cases with intraluminal stenosis less than 50% were considered the control group and those greater than 50% were considered the disease group. In GSE20680, there were 87 cases in the disease group and 108 cases in the control group. In GSE20681, there were 99 cases in the disease group and 99 cases in the control group. The datasets were based on the GPL4133 platform, and we used R language (version 3.6.3; The R Foundation for Statistical Computing, Vienna, Austria) to normalize the data using the “limma” package (Figure S1).

Differentially expressed genes and autophagy intersection genes

We used R language to pre-process and analyze the microarray data, and used the “limma” package to screen the differentially expressed genes (DEGs) for CAD. The “pheatmap” package was used to visualize the differential genes. The VennDiagram package was used to analyze the DEGs of the two datasets separately to obtain the intersecting genes. Autophagy-related genes were downloaded from the HADb database (http://www.autophagy.lu/index.html), and the “VennDiagram” package was used to obtain their intersecting genes.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) bioenrichment analysis

GO and KEGG enrichment analysis of DEGs was performed using the clusterProfiler, enrichplot, and ggplot2 packages in R. Significant functional enrichment was considered when Q<0.05. The enrichment analysis of GO was performed from three aspects: biological process (BP), cellular composition (CC), and molecular function (MF). Significant molecules or genes in signaling pathways were identified by KEGG.

Gene set enrichment analysis (GSEA)

GSEA of differential genes was performed using the clusterProfiler, ReactomePA, and enrichplot packages in R language. The top five enrichment results were visualized.

Gene-gene interaction (GGI) and protein-protein interaction (PPI) network establishment and core gene and microRNA screening

Construct A GGI network of differential genes was constructed using the genemania database (http://genemania.org/). The PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes/Genomes database (STRING; https://string-db.org) and the dataset was visualized using cytoscape (version 3.9.0; https://cytoscape.org/index.html) for the top ten differentially linked genes calculated using Cytohubba. The microRNAs (miRNAs) of CAD differential genes were predicted using the FunRich database (version 3.1.3; http://www.funrich.org/). The results were entered into cytoscape for visualization and analysis.

CAD crossover gene function and drug sensitivity analysis

Using the GeneCards database (https://www.genecards.org/), CAD therapeutic targets were retrieved using the keyword “coronary heart disease” and then intersected with the GSE20680 and GSE20681 datasets. Then, GGI, GO, and KEGG enrichment analyses were performed on the intersecting genes. Finally, potential drugs targeting the intersecting genes were predicted using The Drug-Gene Interaction Database (DGIdb; www.dgidb.org) and visualized using Cytoscape for network modules.

CAD intersection gene prediction model construction

The prediction performance of the intersecting genes was analyzed using the “pROC” package in R to compare the prediction performance of the genes intersecting the two datasets. The top five area under the curve (AUC) values of the intersecting genes in the two datasets were included in the subsequent model building, and the “rms” package was used to build the CAD prediction model.

Immuno-infiltration analysis of DEGs

Two normalized datasets, GSE20680 and GSE20681, were imported into the Cell Tpe Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT; https://cibersort.stanford.edu/), and the expression matrices of human immune cell subtypes were deconvoluted to obtain the percentages of 22 immune cell types. The results were visualized using R language.

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Statistical analysis

All bioinformatics analyses in this study were performed using the R software. The PHEATMap package was used to construct the expression heat maps of important genes in CAD patients and healthy controls. Limma package in R language was used for statistical test to compare the expression differences of important genes between CAD patients and healthy controls. The AUC value in the ROC curve was used to evaluate the diagnostic efficacy, and P<0.05 was considered statistically significant.

Results

DEGs and autophagy crossover genes

Based on the screening criteria of adj. P<0.05 and FC absolute value >2,221 DEGs were screened in the GSE20680 dataset, of which 150 were up-regulated and 71 were down-regulated; 297 DEGs were screened in the GSE20681 dataset, of which 212 were up-regulated and 85 were down-regulated. The expression changes in the test and control groups and their gene and sample clustering results are shown in . Some 79 co-expressed DEGs existed in the GSE20680 and GSE20681 datasets, which were further intersected with autophagy-related genes and three of them were found to be autophagy-related genes ().

Figure 1

Differential expression analysis of the GEO dataset. (A,B) Heat map and volcano plot of differential expression analysis of GSE20680 dataset; (C,D) Heat map and volcano plot of differential expression analysis of GSE20681 dataset. FC, fold change; GEO, Gene Expression Omnibus.

Figure 2

Venn diagram. (A) DEGs of GSE20680 and GSE20681 datasets taken as intersection; (B) intersection of DEGs with mitochondrial autophagy-related genes taken as intersection. DEGs, differentially expressed genes.

Biological enrichment analysis of DEGs

The GO analysis revealed that the functions of the co-expressed DEGs were mainly focused on the following aspects: BP was mainly the immune response-regulating signaling pathway, cellular polysaccharide catabolic process and polysaccharide CC is a ficolin-1-rich granule, external side of plasma membrane, and cytoplasmic vesicle lumen; MF was mainly involved in monocarboxylic acid binding, 3,5-cyclic-GMP phosphodiesterase activity and 1-acylglycerol-3-phosphate-acyltransferase activity, nucleotide metabolism, glycerophospholipid metabolism, the NF-kappa B signaling pathway, and cell adhesion molecules ().

Figure 3

GO and KEGG analysis of differentially expressed genes. (A) GO analysis; (B) KEGG analysis. BP, biological process; CC, cellular composition; MF, molecular function; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The GSEA analysis of DEGs in the GSE20680 dataset showed that DEGs were concentrated in amyotrophic lateral sclerosis, metabolic pathways, Parkinson’s disease, phagosome, and vascular smooth muscle contraction. In contrast, the GSE20681 data mainly involved the MAPK signaling pathway, metabolic pathways, pathways in cancer, pathways of neurodegeneration—multiple diseases, and the PI3K-Akt signaling pathway ().

Figure 4

GSEA analysis of DEGs. (A) GSE20680 dataset; (B) GSE20681 dataset. GSEA, gene set enrichment analysis; DEGs, differentially expressed genes.

GGI and PPI network establishment and core gene and miRNA screening

We first analyzed the interactions between the differential genes to find the genes that might share functions with it, then explored the reciprocal proteins of the corresponding proteins of these genes by STRING, and then screened the top 10 core genes: RPL39, RPL19, RPS12, RPL18A, BYSL, ENTPD1, S1OOA8, LY96, MNDA, and MMP9 core genes (). To further explore the potential role of these DEGs in miRNA, we used the funrich database for prediction. The results yielded 29 genes including S100A8, ENTPD1, and AP1S2ADA and their possible target miRNAs ().

Figure 5

DEG network construction and core gene analysis. (A) DEG GGI; (B) core gene PPI; (C) visualization network of top 10 core genes. DEG, differentially expressed gene; GGI, gene-gene interaction; PPI, protein-protein interaction.

Figure 6

Analysis of DEG-miRNA interaction network. DEG, differentially expressed gene.

Functional and drug sensitivity analysis of CAD crossover genes

A total of 2,282 potential pathogenic genes for CAD were obtained from the GeneCards database and were intersected with the GSE20680 and GSE20681 datasets to obtain 11 genes: EPOR, ALOX5AP, PDE9A, PYGL, F5, MMP9, ADA, S100A8, CBS, ENTPD1, and HSPB1 (). These genes are closely related to the regulation of body fluid levels, cytoplasmic vesicle lumen, and monocarboxylic acid binding, and are mainly involved in purine metabolism, nucleotide metabolism, and the VEGF signaling pathway (). Based on the drug sensitivity analysis of the intersecting genes, we also identified 75 potential target therapeutics ().

Figure 7

CAD intersection gene function and drug sensitivity analysis. (A) CAD intersection gene Venn diagram; (B) intersection gene GGI analysis; (C) intersection gene GO analysis; (D) intersection gene KEGG analysis; (E) intersection gene drug sensitivity analysis. CAD, coronary artery disease; GGI, gene-gene interaction; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

To explore the diagnostic efficacy of these intersecting genes, we analyzed the GSE20680 and GSE20681 datasets and constructed prediction models, respectively. The results showed that five genes, S100A8, HSPB1, F5, MMP9, and PDE9A, had good AUC values and the constructed column line graph prediction models had good predictive power in both GSE20680 and GSE20681 datasets ().

Figure 8

Intersection gene column line graph prediction model. (A) GSE20680 dataset column line graph model; (B) GSE20681 dataset column line graph model.

CAD immune infiltration

The inflammatory response is an important factor in the development of CAD, so we analyzed the level and correlation of 22 immune infiltrating cells in the GSE20680 and GSE20681 datasets, respectively. In the GSE20680 dataset, the levels of T cells CD4 naïve, T cells CD4 memory activated, T cells CD8, monocytes, and neutrophils fluctuated significantly between controls and CAD patients, and there were complex positive or negative correlations between these immune cells. The T cells regulatory (Tregs) were significantly reduced in CAD patients compared to controls (). In the GSE20681 dataset, T cells CD8, macrophages M0, B cells naïve, T cells CD4 memory activated, and natural killer (NK) cells resting were the immune cells that changed more in controls and CAD patients, with a significant difference in the level of B cells naïve ().

Figure 9

GSE20680 dataset immune infiltration analysis. (A) Bar graph. (B) heat map; (C) immune cell correlation analysis; (D) violin plot of immune cell infiltration levels in healthy individuals and CAD patients. CAD, coronary artery disease.

Figure 10

GSE20681 dataset immune infiltration analysis. (A) Bar graph; (B) heat map; (C) immune cell correlation analysis; (D) violin plot of immune cell infiltration levels in healthy individuals and CAD patients. CAD, coronary artery disease.

Discussion

In this study, two sets of genes related to CAD were obtained by screening microarray data from the GEO database through bioinformatics analysis, and they were analyzed separately from controls to obtain DEGs and autophagy-related genes. The GO, KEGG, and GSEA analyses of DEGs revealed that some DEGs were involved in purine metabolism, metabolic pathways, the MAPK signaling pathway, and the PI3K-AKT signaling pathway. The main lesion of CAD vessels is coronary atherosclerosis (AS), and the main hallmark is inflammatory cell infiltration of the arterial vessel wall and abnormal proliferation of vascular smooth muscle cells or macrophages (14). To further explore the relationship with CAD pathogenesis-related genes, enrichment functions, transcription factors, and possible immune cells involved, we performed PPI and miRNA analysis in an attempt to identify potential target genes and performed immune infiltration and drug sensitivity analysis.

Mitochondria are the main site of nutrient metabolism and adenosine triphosphate (ATP) production in most cells, and normal mitochondrial structure and function is a prerequisite for various intracellular activities. Myocardial ischemia-reperfusion injury (MIRI) can cause significant mitochondrial damage, and the damaged mitochondria can be encapsulated by intracellular membrane structures to form autophagic vesicles and eventually mitochondrial autophagy (15). When myocardial tissue ischemia occurs, mitochondrial autophagy can help cardiomyocytes adapt to tissue ischemic and hypoxic conditions and improve the chances of cell survival (16). In addition, it has been found that basal level autophagy helps to slow down the rate of AS plaque formation and maintain the stability of AS plaques (17). However, current studies are inconclusive as to whether autophagy actually facilitates cell survival after the onset of MIRI (18,19). Autophagy is involved throughout the development of AS, and in our present study we also observed that autophagy is involved in the regulation of CAD, but its specific mechanism of action is not yet clear.

In functional enrichment analysis, it was shown that the intermingled genes are involved in regulating signaling pathways such as MAPK and PI3K-AKT. Mitogen-activated protein kinase (MAPK) is involved in regulating various physiological functions such as cell growth, division, death, and stress apoptosis, with ERK1/2 and p38 MAPK being key targets of endothelial cell inflammation and apoptosis in the mechanism of AS development (20-22). In addition, the MAPK pathway may be involved in the differentiation of monocytes to macrophages, and blocking this process may contribute to the treatment of AS (23). The phosphatidylinositol 3-kinase (PI3K)/phosphokinase B (AKT) signaling pathway has a wide range of biological effects and is involved in regulating cell division, proliferation, apoptosis, metabolism, and other activities that are closely associated with cardiovascular diseases such as hypertension, ischemic cardiomyopathy, and AS (24-26). These studies confirm the importance of the MAPK and PI3K-AKT signaling pathways in the development of pre-CAD and also remind us to focus on the management of inflammation in CAD.

The S100A8 gene encodes a protein that is a member of the S100 protein family and is involved in regulating many cellular processes. Stress-induced elevation of S100A8/A9 has been evidenced in patients with CAD with impaired cortisol response and is associated with poor prognosis in CAD patients (27). Recent studies have shown that S100A8/A9 is positively correlated with bactericidal/permeability-increasing protein and is a potential biomarker of myocardial infarction in patients with acute myocardial infarction (28). In addition, some researchers found that HSP27 expression was increased in single nucleated cells of human peripheral blood and significantly correlated with the severity of CAD, rendering it a potential prognostic marker (29). Matrix metalloproteinase-9 (MMP9), an important member of MMP, is involved in the whole process of AS and high MMP9 expression is associated with plaque instability (30,31). In a subsequent analysis, these genes were found to be positively or negatively correlated with some immune cells, suggesting that the genes and the immune microenvironment influence each other, but the exact mechanism is unknown. Our study suggests that these genes have an important role in the development of CAD and have good predictive efficacy, and also suggests that these target genes may be potential therapeutic targets.

The instability of atherosclerotic plaques is a common pathological basis of CAD, and the process of atherosclerotic plaque formation involves not only lipid deposition but also a systemic chronic inflammatory response, with circulating neutrophils, lymphocytes, and other inflammatory cells involved in the inflammatory response and immune process during the course of CAD (32,33). A study confirmed that T lymphocytes accounted for 10–20% of all nucleated cells within atherosclerotic plaques (34). Moreover, different phenotypes of T lymphocytes have different roles in the development of CAD, with T helper lymphocytes 17 (Th17) promoting the development of AS, whereas T regulatory lymphocytes (Tregs) play a protective role (35,36). By sequencing peripheral blood from CAD patients, a study found that Tregs were reduced in CAD, and subsequently identified several markers of Tregs (37). In addition, clinical studies have shown that Tregs may be involved in immune regulation early in the atherosclerotic process (38). Macrophages are the “key drivers” of AS development, affecting plaque stability and AS outcome (39). Macrophages are regulated by multiple signaling pathways, among which the PI3K/Akt signaling pathway plays an important role in macrophage survival, proliferation, and migration, and is involved in regulating macrophage polarization, autophagy, and lipid metabolism (40-42). In our study, we found that a variety of immune cells such as Tregs, macrophages, B cells naïve, and T cells CD8 were imbalanced in CAD, which is consistent with the above-mentioned reports. B cells naive and T cells CD8 are significantly elevated in patients with CAD, which may be an indicator of potential disease severity and also suggest clinicians to enhance the treatment of aseptic inflammation, but the mechanisms of immune cells involved in CAD need to be investigated in depth in follow-up. In addition, these results further suggest that the role of the PI3K/Akt signaling pathway in CAD should not be underestimated.

In this study, we analyzed the transcriptomic information of peripheral blood from CAD patients by using a larger sample dataset. A total of 79 DEGs were obtained after differential expression analysis with controls, but only 11 intersected genes were obtained after intersection with potential pathogenic genes of CAD, which we consider to be related to the large number of unknown CAD-related genes. A predictive model consisting of five genes, S100A8, HSPB1, F5, MMP9 and PDE9A, can help identify patients who are likely to develop CAD, as well as select the corresponding targeted drugs for patients with combined mutations in such targets. In addition, peripheral blood, rather than coronary artery tissue, was used as the data set specimen for this study. Lack of research funding and limitations in clinical specimen collection prevented actual clinical gene validation in this study. In addition, we were unable to obtain information on coronary stenosis in these database patients to combine these markers for further correlation analysis.

Conclusions

In summary, this study found that the mRNA and immune infiltrating cells detected in peripheral blood of CAD patients were significantly different from those of the healthy population, and that DEGs and crossover genes were involved in many key biological processes. The DEGs and intersecting genes identified in this study can be used as regulatory and therapeutic targets, and the models constructed in this study have predictive power.

The article’s supplementary files as