Microarray data for CAD were searched in NCBI GEO (http://www.ncbi.nlm.nih.gov/geo/). Inclusion criteria were: i) The dataset belonged to a gene expression profile; ii) samples in the dataset were collected from blood; iii) the samples included patients and healthy controls; iv) the dataset was based on human gene expression profiles; and v) the sample number was ≥20. A total of four datasets, GSE12288 (n=222; CAD, 112; control, 110) (13), GSE20680 (n=139; CAD, 87; control, 52) (14), GSE20681 (n=198; CAD, 99; control, 99) (15) and GSE42148 (n=24; CAD, 13; control, 11; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42148) that met the criteria were included in the present study (Table I).

Raw data in CEL files from GSE12288 were downloaded from the Affy platform (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL96) and subsequently converted into gene symbols under pretreatment. The median method was used to fill in the missing values, the MicroArray Suite method to complete background correction, and the quantiles method to normalize the data, using oligo 1.41.1 in R3.3.1 language (http://www.bioconductor.org/packages/release/bioc/html/oligo.html) (16). GSE20680, GSE20681 and GSE42148 were downloaded from the Agilent platform (Agilent Technologies, Inc., Santa Clara, CA, USA). Microarray raw data (TXT files) of the three datasets were log2-transformed using Limma 3.34.0 (https://bioconductor.org/packages/release/bioc/html/limma.html) to achieve an approximate normal distribution, and were subsequently standardized using the median normalization method (17).

WGCNA is a well-established method for constructing scale-free gene co-expression networks, which is characterized by the use of soft thresholding (18,19). GSE12288 was used as the training dataset, while GSE20680, GSE20681 and GSE42148 were the validations sets. Genes with coefficients of variation <0.1 were discarded. Correlations between the gene expression in the four datasets were evaluated with the verboseScatterplot function of the WGCNA package (https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/). Weighted gene co-expression networks of the genes in the training dataset were constructed using the WGCNA method as previously described (20–22). First, the soft threshold power of β was set as 18 (scale-free R²=0.9), according to the scale free topology criterion. The weighted adjacency matrix was subsequently constructed. Adjacencies and correlations were transformed into a topological overlap matrix (TOM), followed by calculating the corresponding dissimilarity (1-TOM). Subsequently, a hierarchical clustering analysis (23) of genes was performed using 1-TOM as the distance measure. Modules were detected using dynamic tree cut algorithm with a minimum module size of 50 and a minimum cut height of 0.95. Furthermore, module preservation between the training set and the three validation sets was measured using the module preservation function of the WGCNA software package. Possible functions of the preserved module were studied using the user List Enrichment function of the WGCNA package.

DEGs between CAD samples and healthy controls were screened in each of the above four datasets using the metaDE method (10,24) in R language (https://cran.r-project.org/web/packages/MetaDE/). Heterogeneity was examined across the four datasets to assess the consistency of gene expression, by calculating tau2 and Qpval values. A tau2=0 and Qpval >0.05 indicated that the gene was homogeneous and unbiased. The thresholds for DEG identification were tau2=0, Qpval >0.05 and false discovery rate <0.05. Consistency of the DEGs was detected using the MetaDE method, and the DEGs with significant consistency in the four datasets were selected for further analysis.

A PPI network was constructed to evaluate the interactions between genes in the above network. The common genes in the preserved modules that were obtained from WGCNA, and the DEGs with significant consistency, were selected to construct a PPI network based on three databases: The Biological General Repository for Interaction Datasets 3.4.153 (http://thebiogrid.org/) (25); the Human Protein Reference Database Release 9 (http://www.hprd.org/) (26); and STRING 10.5 (https://string-db.org/) (27). The PPIs revealed in at least two of the three databases were extracted for the PPI network, visualized using Cytoscape 3.3 software (http://www.cytoscape.org/) (28). In the PPI network, a node represents a gene; the undirected link between two nodes is an edge, denoting the interaction between two genes; and the degree of a node corresponds to the number of interactions of a gene with other genes in the network.

In order to elucidate the possible biological roles of the genes in the PPI network, GO (29) functional and KEGG (30) pathway enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery 6.8 software (31). GO terms have three categories, including biological process (BP), cellular compartment and molecular function. P<0.05 was considered to indicate a significant difference for GO terms and KEGG pathways.

Following data preprocessing for the four datasets, the present study attempted to identify CAD-associated co-expression modules using WGCNA. Using the aforementioned methods, the correlation coefficients (CCs) of genes in the four datasets were obtained, and the CCs were 0.47–0.92 with P-values <1×10⁻²⁰⁰ between any two datasets (Fig. 1), suggesting a good consistency of the common genes across all datasets.

As a result, 11 co-expression modules were identified for genes in the training set using WGCNA (Fig. 2A). These modules are illustrated in branches of the dendrogram with different colors. To evaluate the robustness of these modules in the training set, the modules were re-constructed in the three validation datasets, GSE20680, GSE20681 and GSE42148, separately (Fig. 2B-D). Genes in the three validation sets were colored, according to the module color in the training set. Multi-dimensional scaling of expression data of all genes in these modules demonstrated that genes in the same module appeared to cluster together and exhibited a similar expression pattern (Fig. 3A). Additionally, hierarchical clustering analysis of these modules in the four datasets revealed that modules of the same branch tended to have similar gene expression patterns (Fig. 3B).

In total, five of the 11 modules, black, brown, magenta, turquoise and yellow modules, with Z-scores >5 were determined to be well preserved across the four datasets, including 873 genes (Table II). This suggested that the five highly preserved modules may be closely associated with CAD. Of the five modules, genes in the black module were significantly linked to ‘response to glucocorticoid stimulus’; those in the brown module were significantly associated with ‘regulation of transcription’; and those in the magenta module may be involved in ‘protein localization’ (Table II). Notably, the turquoise and the yellow module exhibited functions in ‘immune response’ (Table II).

Using the metaDE method, 836 DEGs were identified with significant consistency across the four datasets. A heatmap for these DEGs demonstrated that the expression patterns of these DEGs differed between the CAD and control samples (Fig. 4).

A total of 177 genes were overlapping between WGCNA modules and the DEGs with significant consistency (Fig. 5A). Of them, 92 genes were included in the turquoise module, 26 in the yellow module, 16 in the black module, 31 in the brown module and 12 in the magenta module (Fig. 5B). Based on the three databases mentioned above, a PPI network was built with these overlapping genes. In total, 150 paired PPIs that appeared in at least two of the three databases were included in the PPI network (Fig. 6A), which contained 59 downregulated genes and 40 upregulated genes (Fig. 6B). All genes in the network were ranked in a descending order, according to their degrees. The top five genes were LCK proto-oncogene, Src family tyrosine kinase (LCK; degree=15), euchromatic histone lysine methyltransferase 2 (EHMT2; degree=14), inosine monophosphate dehydrogenase 2 (IMPDH2; degree=12), protein phosphatase 4 catalytic subunit (PPP4C; degree=11) and ζ-chain of T-cell receptor associated protein kinase 70 (ZAP70; degree=11).

The enrichment analysis indicated that genes in the PPI network were significantly associated with numerous GO BP terms, including the ‘transmembrane receptor protein tyrosine kinase signaling pathway’ and ‘peptidyl-tyrosine autophosphorylation’ (Table III). With regard to KEGG pathways, there were five significant pathways, including ‘metabolic pathways’, ‘natural killer cell mediated cytotoxicity’, ‘fructose and mannose metabolism’, ‘primary immunodeficiency’, and ‘Fc gamma R-mediated phagocytosis’ (Table IV). Notably, ‘natural killer cell mediated cytotoxicity’ and ‘primary immunodeficiency’ pathways were significantly enriched with LCK and ZAP70.