Weighted gene coexpression network analysis identifies the key role associated with acute coronary syndrome

Data preprocessing

We obtained 47,280 probes from each sample expression profile in GSE56045. After data preprocessing, we obtained 20,918 probes containing gene symbols from 1,202 samples. The gene expression matrix was associated with the sample phenotype matrix for further analysis.

Weighted gene coexpression networks

After calculation, we believe that, when the correlation coefficient is 0.9 (soft threshold β is 4), the coexpression network has a higher correlation and is more suitable to construct different gene modules (Figure 1A). Together with the topological overlap matrix (TOM), we performed the hierarchical average linkage clustering method to identify the gene modules of each gene network (deepsplit = 2, cut height = 0.4) and then showed the heatmap in Figure 1B. Next, we calculated the gene cluster tree and showed the results in Figure 1D. Finally, about eleven gene modules should be handled by the dynamic tree cut (Figure 1C).

Interest module and functional annotation

As shown in Figure 2, the highest association in the module-feature relationship was the dark-red module and T cell (r² = 0.86, P < 0.001), orange module and B cell (r² = 0.87; P < 0.001). The dark-red module contains 121 genes, while the orange module contains 98 genes. All of these genes were showed in Supplementary Table 1. To determine the correlation between gene significance and color module, we conducted an in-depth calculation. As shown in Figure 3A, the correlation between the gene significance and orange module was 0.57 (P = 9E-10) and that of the dark red module was 0.28 (P = 0.0019). After confirming the correlation between the gene significance and modules, we analyzed the functional enrichment of the 219 genes in these two modules. Gene Ontology (GO) function, KEGG pathway enrichment and Disease Ontology analyses were performed by R (Figure 3B). The details of these analyses can be found in Supplementary Table 2. The biological processes of these two modules were found to be associated with GO:0042110-T cell activation (P = 4.89E-19), GO:0002768-immune response-regulating cell surface receptor signaling pathway (P = 4.66E-14), GO:0042113-B cell activation (P = 8.13E-07), and GO:0002699-positive regulation of immune effector process (P = 7.25E-05). However, in KEGG pathway analysis, these two modules were found to be associated with hsa04660-T cell receptor signaling pathway (P = 4.54E-09), hsa04060-Cytokine-cytokine receptor interaction (P = 6.87E-08), hsa04064-NF-kappa B signaling pathway (P = 2.38E-05), hsa04662-B cell receptor signaling pathway, and hsa04514-Cell adhesion molecules (CAMs) (P = 0.002).

Protein-protein interaction (PPI) network construction and identification of hub genes

Approximately 149 nodes and 1,016 protein pairs were obtained when the combined weight score was set at more than 0.25 (Figure 3C-All). Analysis in the submodule revealed four modules with a score > 6 detected by MCODE (Figure 3C Cluster 1-4). After integrating the GO function, KEGG pathway enrichment and PPI analysis results, we identified several genes with a high degree and MCODE scores: PRF1, GZMB, CD27, CD2, CCL5, CXCR5, CD8A, CCR7, IL2RB, IFNG, CD40LG, IL7R, CD226, LAT, CD6, PLCG1, CD22, CCR4, TNFRSF13B, PRKCQ, LTB, LTA, STAT4, TNFRSF13C. Next, we needed to further verify these genes to determine their functional characteristics.

Dataset validation

First, we analyzed the relationship between the expression level of these 24 genes and the coronary artery calcification score by linear correlation in GSE56045. As shown in Figure 4, only 4 gene expression levels were related to the coronary artery calcification score: CCR7 (R = -0.081, P = 0.0065), CD2 (R = -0.075, P = 0.0012), CXCR5 (R = -0.065, P = 0.029) and IL7R (R = -0.06, P = 0.043). Next, we tested the expression level and function of the four genes in different datasets to confirm whether they were closely related to CAD. As shown in Figure 5A, the expression levels of CCR7, CXCR5 and IL7R were different in the CAD and control group (P < 0.05–0.001), but no significant difference was found with CD2 (P = 0.41). In GSE34822 (Figure 5B), the expression of CCR7, CXCR5 and CD2 was different in progressive and stable plaques, while that of IL7R was not significantly different (P < 0.05–0.001). Additionally, the gene expression level of progressive plaques was lower than that of stable plaques. Figure 5C analyzes GSE129935, which revealed that only CCR7 and CXCR5 had statistical significance in the comparison of stability and acute myocardial infarction (AMI), stability and instability (P < 0.05–0.01). Additionally, the subsequent decrease in expression led to the occurrence of plaque instability, which leads to AMI. GSE60993 (Figure 5D) also well verified the previous conclusion. From this dataset, the expression of these four genes in the healthy control group was statistically significant compared with unstable angina (UA), non-ST elevation myocardial infarction (NSTEMI) and ST elevation myocardial infarction (STEMI) (P < 0.001), and the healthy control group had a higher expression level. GSE59867 (Figure 5E) reflected the change in gene expression from stable to the first day to six months after AMI. From this dataset, we found that only CCR7 and CXCR5 expression was statistically significant (P < 0.05–0.001), the expression level was the lowest on the first day after AMI, and the gene level gradually increased with time. By dividing the expression of CCR7 and CXCR5 in GSE59867 by the mean value, expression greater than the mean value was defined as high expression, and the occurrence of heart failure was defined as the end event. We found that, when CCR7 showed high expression and CXCR5 showed low expression, the occurrence of heart failure was 10.796 times higher than that of CCR7 and CXCR5 (P < 0.05) (Figure 5F). Based on the above five datasets, we believe that CCR7 and CXCR5 play an important role in atherosclerosis and plaque vulnerability. Further verification of these data in a future investigation is warranted.

Subject validation

To verify the exact function of CCR7 and CXCR5, we collected the data of some patients hospitalized due to chest pain, measured the relative gene expression after collecting the peripheral blood, and compared the gene expression level according to the coronary angiography data. Table 1 shows the general situation of 1,528 patients with gender and age matching. We considered all the variable data, including the relative expression of CCR7 and CXCR5, gender, age, smoking, drinking, BMI, systolic blood pressure (SBP), diastolic blood pressure (DBP), serum glucose, TC, TG, high-density lipoprotein cholesterol (HDL-C), LDL-C, apolipoprotein (Apo)A1 and ApoB, which were the best subset of risk factors to develop the acute coronary syndrome (ACS) risk score and risk model (nomogram) (Figure 6). We defined the sores as follows: smoking and/or drinking: yes = 2, no = 1; male = 1; female = 2. The nomogram had excellent discriminative power with a C-statistic and was well calibrated with the Hosmer-Lemeshow χ ² statistic. The predicted probabilities of developing ACS ranged from 0.0004 to 99%. The discrimination accuracy of the model was 0.841 (95% CI, 0.809–0.871). At an optimal cutoff value, the sensitivity and specificity were 64.0% and 90.9%, respectively.

After calculation, the relative expression of levels of CXCR5, ApoA1, LDL-C, serum glucose, smoking and Troponin T were significantly related to the risk of ACS, with statistical significance. The relative expression of peripheral blood RT–PCR showed that the expression of CXCR5 was statistically significant in the comparison of cases (including UA, NSTEMI and STEMI) and controls, and the expression decreased gradually with the increase in plaque vulnerability (Figure 7), a finding that was also consistent with the results of previous multidatasets.

Microarray data

All dataset analyses were performed using R software (version 3.60). We downloaded GSE56045 [21] microarray data from GEO (Gene Expression Omnibus, http://www.ncbi.nlm.nih.gov/geo/). These datasets were based on the platform of the GPL10558 Illumina HumanHT-12 V4.0 expression BeadChip. GSE56045 contains 1,202 peripheral blood samples, as well as clinical information, such as age, whole-blood cell count, and the coronary artery calcification plaque score. First, we added CEL files into the R software using the Affy package for transformation into an expression value matrix. The probe information was then transformed into gene names using the Bioconductor package. The mean value should be chosen when a gene had more than one probe [22]. The preprocessing process of the datasets (GSE23561, GSE34822, GSE59867 and GSE60993) [23–27] used for validation was the same as that for GSE56045. GSE129935 was performed by Fragments Per Kilobase per Million mapped reads (FPKM) and quantile normalized using the robust multiarray average (RMA) method. The probes were then annotated using Bioconductor in R. GSE129935 was also used for validation.

Weighted gene coexpression network analysis

We conducted the analyses in strict accordance with the weighted gene coexpression network analysis process [28]. First, we chose the appropriate soft threshold power according to standard scale-free networks, with which adjacencies between all differential genes were calculated by a power function. Next, a topological overlap matrix (TOM) was derived from the adjacencies, and the corresponding dissimilarity (1-TOM) was counted. To complete the module recognition, we used the dynamic tree cutting method to cluster the genes in layers, using 1-TOM as the distance measure, a minimum size cutoff of 30, and a deepSplit value of cutting of 2. Next, we selected the highly similar modules by clustering and merged them, with the height line set as 0.4. To test the stability of each identified module, module preservation and quality statistics were computed using the module preservation function implemented in the WGCNA package [28]. Because the genes in the gray module cannot be attributed to any other module, all the genes in the gray module were removed.

Interest module and hub gene

We selected the module with the highest correlation with clinical features and the genes in this module with the important biological functions. To identify the most biologically significant module, we used Pearson correlation analysis to evaluate the correlation between clinical features and gene modules. We selected the most relevant modules related to clinical features for further analysis and research. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed to detect potential mechanisms by which these module genes affect correlative clinical features. The thresholds for the P value and false discovery rate (FDR) were set as less than 0.01 and 0.05, respectively. These analyses were completed using clusterProfiler and DOSE package in R [29].

We describe the correlation between the gene expression profile and module eigengenes (Mes) defined as the module membership (MM). The gene significance (GS) can be defined as the absolute value of the correlation between external traits and the gene. The genes of interest in the modules with the highest MM and GS scores were selected for subsequent analysis. We defined the GS and MM scores as more than 0.2 and 0.6, respectively, for intramodular hub genes selected by external traits, and the P cutoff value is set at less than 0.05. The STRING database (https://string-db.org/, Version 11.0) [30] was employed to analyze the protein-protein interaction (PPI) network. We used Cytoscape software (version 3.80) to visualize and construct the gene-gene interaction network. Additionally, we used the Molecular Complex Detection (MCODE) app to screen the most notable clustering modules, with an MCODE score greater than 6 set as the threshold for further analysis [31].

Hub gene validation and survival analysis

Our validation analysis was divided into two stages. In the first stage, we analyzed the function of the hub genes in different datasets; in the second stage, we verified the population. After collecting the patients admitted for chest pain, we extracted the peripheral blood and verified the expression of core genes. GSE23561, GSE34822, GSE59867, GSE60993 and GSE129935 were employed for validation. At first, we analyzed the correlation between the gene expression level and calcification score. Next, we compared the expression differences of hub genes among different groups and showed the results using the ggplot2 package in R. Subsequently, the “survival” package [32] in R was used to perform overall survival (heart failure) and disease-free survival analyses for all hub genes. The patients were divided into four groups (high vs. low) based on the hub gene expression level compared with the mean expression level of that hub gene. A Kaplan–Meier survival plot was also constructed. The ‘rms’ package was used for ACS prediction nomogram. The predictive accuracy of the risk model was assessed by discrimination measured by C-statistic and calibration evaluated by Hosmer-Lemeshow χ² statistic. The discriminatory ability of the model was quantified using the area under the receiver operating characteristic curve (AUC). A 95% CI was calculated for each AUC. In general, an AUC > 0.75 was considered to be relatively good discrimination.

Study population

In total, 1,528 patients were recruited from the inpatient department for a complaint of chest discomfort at the Liuzhou People's Hospital from 2018-3-1 to 2019-12-31 and had undergone coronary angiography. CAD, UA and AMI were diagnosed based on the Fourth Universal Definition of Myocardial Infarction (2018) [33]. Exclusion criteria included subjects with poor compliance, incomplete clinical data, contrast agent sensitivity and autoimmune diseases. Additionally, subjects with obvious surgical contraindications were excluded. Clinical data collection, biochemical measurements and diagnostic criteria were performed according to previous studies [34]. The study adhered to the Declaration of Helsinki of 1975 (http://www.wma.net/en/30publications/10policies/b3/) and its revision in 2008 and the Ethics Committee of Liuzhou People's Hospital agreed with the study design (No: Lunshen-2018-KY; Feb. 12, 2018). Informed consent was obtained from all subjects after receiving a full explanation of the study.

RT–qPCR and statistical analysis

The procedures of blood sample collection, RNA isolation, reverse transcription cDNA and RT–qPCR are the same as those in our previous studies and were carried out in strict accordance with the product instructions and laboratory operating procedures [35]. The specific divergent primers were designed to amplify the transcripts and are shown in Supplementary Table 3. The statistical software package SPSS 22.0 (SPSS Inc., Chicago, IL, USA) and R software (version 3.6.0) were used for statistical analysis. Quantitative variables were expressed as means ± standard deviation (TG levels were shown as medians and interquartile ranges and were analyzed by the Wilcoxon–Mann–Whitney test because they were not a normal distribution). Chi-square analysis was used to assess the difference in the percentages between the groups. All tests were two-sided, and P < 0.05 was considered statistically significant.

Ethics approval

This analysis of publicly available data does not require ethical approval.