Skip to main content

Identification of the susceptible genes and mechanism underlying the comorbid presence of coronary artery disease and rheumatoid arthritis: a network modularization analysis



The comorbidities of coronary artery disease (CAD) and rheumatoid arthritis (RA) are mutual risk factors, which lead to higher mortality, but the biological mechanisms connecting the two remain unclear. Here, we aimed to identify the risk genes for the comorbid presence of these two complex diseases using a network modularization approach, to offer insights into clinical therapy and drug development for these diseases.


The expression profile data of patients CAD with and without RA were obtained from the GEO database (GSE110008). Based on the differentially expressed genes (DEGs), weighted gene co-expression network analysis (WGCNA) was used to construct a gene network, detect co-expression modules, and explore their relation to clinical traits. The Zsummary index, gene significance (GS), and module membership (MM) were utilized to screen the important differentiated modules and hub genes. The GO and KEGG pathway enrichment analysis were applied to analyze potential mechanisms.


Based on the 278 DEGs obtained, 41 modules were identified, of which 17 and 24 modules were positively and negatively correlated with the comorbid occurrence of CAD and RA (CAD&RA), respectively. Thirteen modules with Zsummary < 2 were found to be the underlying modules, which may be related to CAD&RA. With GS ≥ 0.5 and MM ≥ 0.8, 49 hub genes were identified, such as ADO, ABCA11P, POT1, ZNF141, GPATCH8, ATF6 and MIA3, etc. The area under the curve values of the representative seven hub genes under the three models (LR, KNN, SVM) were greater than 0.88. Enrichment analysis revealed that the biological functions of the targeted modules were mainly involved in cAMP-dependent protein kinase activity, demethylase activity, regulation of calcium ion import, positive regulation of tyrosine, phosphorylation of STAT protein, and tissue migration, etc.


Thirteen characteristic modules and 49 susceptibility hub genes were identified, and their corresponding molecular functions may reflect the underlying mechanism of CAD&RA, hence providing insights into the development of clinical therapies against these diseases.

Peer Review reports


The comorbidity of coronary artery disease and rheumatoid arthritis (CAD&RA) can lead to higher mortality rates than those of independent diseases [1], but the biological mechanisms connecting the two remain unclear. Patients with rheumatoid arthritis (RA) have a markedly higher incidence and mortality of cardiovascular disease than general population [2, 3]. Patients with RA accompanied with coronary artery disease (CAD) and revascularization by percutaneous coronary intervention (PCI) are significantly correlated with a higher risk of long-term, major adverse cardiac events [4]. Meanwhile, RA has been considered an independent risk factor for CAD development [5, 6]. Patients with RA have a larger coronary plaque and inflammation burden compared to patients without RA [7,8,9]. Some CAD-related risk factors, such as dyslipidemia, type 2 diabetes mellitus, hypertension, could also contribute to the CAD risk for patients with RA [10, 11]. Drugs such as corticosteroids, which are utilized for treating RA, might increase cardiovascular risk factors and aggravate heart diseases [12].

The underlying mechanisms of CAD-associated progression of RA are not fully elucidated. Studies found that the high mortality of CAD&RA is due to endothelial dysfunction and the circulating acute phase reactants such as C-reactive proteins [13, 14]. Inflammation can promote coronary atherosclerosis and induce coronary microvascular dysfunction in patients with RA, leading to an inadequate supply of myocardial oxygen, with the primary incipient procedures for the two changes being endothelial dysfunction and immune system dysregulation [15, 16]. Neutrophil activation-related genes of S100A8 and S100A12 are under investigation as therapeutic targets for both RA and CAD, hinting at the common pathogenic mechanisms of CAD&RA [17].

To discover the complex pathological mechanisms of CAD and RA, the conventional single target paradigm is not enough to illuminate the molecular basis of CAD&RA, and a novel systematic paradigm is urgently required. Rather than a simple view of the disease due to individual genomic variations, it requires network perspectives to understand the complex phenome-genome relationships of diseases and their comorbidities. Network medicine is thought to be capable of uncovering complex disease relationships using disease modules and network-based approaches, which may help to discover the shared biological mechanisms of associated diseases [18, 19]. A complex disease is rarely the direct consequence of a single gene alternation; rather, it is the result of the interaction of multiple molecular processes. Disease genes usually interact with each other and form closely connected subgraphs, i.e. disease modules, which play important roles in disease–disease relationships [20]. The identification of precise disease modules may help us understand the molecular interactions of complex diseases. Understanding comorbidities can also help physicians evaluate disease progression and improve treatment. Disease-related genes have been used to assess the similarity between different diseases [21, 22]. Accordingly, module-based strategies rather than single gene and targeted strategies are becoming increasingly important for revealing the relationship between multiple gene interactions and disease mechanisms [23, 24].

In this study, the gene expression array profile of CAD patients with and without RA was used to construct gene co-expression networks by weighted gene co-expression network analysis (WGCNA). Network modularization analyses were performed to identify the characteristic modules and susceptibility hub genes for CAD&RA to reveal the potential molecular mechanisms of the comorbid presence of CAD and RA. The workflow is shown in Fig. S1.

Materials and methods

Gene expression profile data and differentially expressed genes analysis

The CAD and CAD&RA datasets GSE110008 were downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) ( The datasets included eight CAD&RA and eight control CAD samples according to the analysis of biopsies of the ascending aorta, and the platform was Affymetrix Human Genome U133A 2.0 Array (HG-U133A_2). The primary data was annotated to form an expression matrix, each probe was matched to their homologous gene symbols, and the repeated gene symbols in the matrix were excluded.

Differentially expressed genes (DEGs) between CAD&RA and CAD patients were identified using R (version 4.1.1) limma package. Genes with a false discovery rate (FDR) adjusted to p < 0.05 were considered as DEGs. Then, the DEGs were compared with the CAD-related and RA-related genes. To obtain CAD-related and RA-related genes, data were retrieved using the key words “coronary artery disease” and “rheumatoid arthritis” in the HPO (, OMIM ( and dbSNP ( databases during October, 2021.

WGCNA network construction and clinical traits analysis

The R package WGCNA was applied to construct the DEG co-expression network. The DEG dataset was checked through the goodSamplesGenes step in WGCNA [25] to remove unqualified genes which do not qualify for inclusion because of missing values in multiple samples. The co-expression network of DEGs was constructed using appropriate soft-threshold β. Topological overlap measure (TOM) and Dynamic Hybrid Tree Cut algorithm were used to perform hierarchical clustering and partition the branches of dendrogram as a module with the following parameters (minModuleSize = 3, mergeCutHeight = 0.25 and verbose = 3). Then, the correlation coefficient between the expression level of each module and the different disease traits was analyzed.

Analysis of module preservation using Zsummary statistic

To quantitatively analyze whether modules significantly varied between different disease groups, Zsummary [26] statistic was calculated to screen the differentiated modules between CAD and CAD&RA. Modules with a Zsummary ≥ 2 were regarded as preserved common modules, and if a module had a Zsummary score < 2, it was defined as a differentiated characteristic module for CAD&RA. Each identified modules was visualized by the Cytoscape software (version 3.7.2) to display the overall gene relationships that were obtained within a module [27].

Functional enrichment analysis

Genes in the selected modules and all DEGs were respectively uploaded for functional enrichment analysis in the Metascape website ( The website is an open tool that helps the biomedical research community analyze gene/protein lists and make better data-driven decisions. A Gene Ontology (GO) function enrichment analysis and a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis [28] were conducted to identify the function and pathways correlated with these modules (Count = 2; EASE = 0.01; and species = Homo sapiens). P < 0.05 was considered as the cutoff criterion. DiNGO [29] software was used to perform the HPO functional enrichment analysis.

Identification of hub genes

The module membership (MM) was defined as the correlation between a gene and a given module [30]. At the same time, the gene significance (GS) of the gene in the module represented the correlation of the gene with clinical traits [31]. Genes with GS ≥ 0.5 and MM ≥ 0.8 in the clinically related gene module networks were defined as hub genes. Then, the expression level of hub genes with higher GS and MM rank were compared between the CAD and CAD&RA groups. A typical t-test was conducted to compare the difference in the expression with a p < 0.05 to indicate statistical significance. To validate whether these hub genes can classify patients into CAD or CAD&RA, three models—logistic regression (LR) [32], K-nearest neighbor (KNN) [33] and support vector machine (SVM) [34] were applied. Moreover, we used two other data-driven methods to screen featured genes for CAD&RA, i.e. the homogeneity of variance test in machine learning and Chi-square test (χ2) in statistics. Two other methods, i.e. Sab [20] and shortest distance [35], were used to calculate the proximity between the screened genes and CAD/RA-related disease genes from network angel. Smaller Sab and shortest distance values indicate closer proximity between the gene and diseases.


DEGs between CAD and CAD&RA

Among all the 22,215 genes in GSE110008, 278 DEGs (Fig. 1A, Table. S1) were ultimately obtained after the duplicate genes were removed. In the up-regulated genes, XIST, DEFA1, ACTA1, and FAM118A genes’ differential expression were the most significant, and DDX3Y, RPS4Y1, TXLNGY, and KDM5D genes’ differential expression were the most significant among the down-regulated genes.

CAD and RA-related genes were obtained from the HPO, OMIM and dbSNP databases (Table. S2). Nine DEGs overlapped with both CAD and RA-related genes, and 54, 33 and 485 overlapped genes were detected between DEGs and CAD-related genes, DEGs and RA-related genes, CAD-related and RA-related genes, respectively (Fig. 1B).

Fig. 1
figure a

DEGs between CAD and CAD&RA. (A) Volcano plot for the gene’ expression level of GSE110008. Purple dots represent the down-regulated genes (log2Fold Change<-1, p < 0.05), red dots represent up-regulated genes (log2Fold Change > 1, p < 0.05). (B) The overlapped genes between the DEGs and known CAD, RA-related genes

Co-expression modules identification and their correlation to clinical traits

The one-step network construction function in WGCNA was used to construct the gene co-expression network based on the DEGs of CAD&RA. According to the scale-free independence and mean connectivity of the gene matrix, the soft thresholding power β was set at 7. Therefore, 41 modules were obtained, ranging in size from 3 ~ 24 (Table. S3). The cluster dendrogram of module distribution and the coexpression network heatmap are shown in (Fig. 2A-B).

Fig. 2
figure b

Identification of gene co-expression networks, modules and the correlation with clinical traits. (A) Cluster dendrogram of 278 DEGs based on the topological overlap. Each branch of the cluster tree with a certain color represents a co-expression module. (B) Heatmap of the topological overlap matrix (TOM) among all 41 modules of DEGs. (C) Heatmap of module-trait relationships. Each row represents a module, and each column represents a trait. Each cell contains the corresponding correlation coefficient

To confirm the modules’ preservation and reproducibility, 1/2 and 3/4 of the samples were selected as testing sets and the Zsummary value [26] for each module was calculated. In the 1/2 sample testing set, 100% modules had a Zsummary score ≥ 0, 78.05% modules had a Zsummary score ≥ 2; In the 3/4 sample testing set, 100% modules with Zsummary score ≥ 0, and 92.68% modules with Zsummary score ≥ 2. These values demonstrate the robustness of our identified modules (Fig. S2).

To define the modules’ clinical characteristic, the correlation coefficient between modules’ expression and disease clinical traits was calculated. Overall, all the 41 modules, 35 modules (85.37%) had an absolute correlation coefficient to CAD or CAD&RA of over 0.6 (Fig. 2C). Among these modules, 17 modules were positively correlated with CAD&RA and negatively correlated with CAD. Conversely, 24 modules were positively correlated with CAD and negatively correlated with CAD&RA. Among the modules positively correlated with CAD&RA, the turquoise module had the largest corresponding correlation coefficient (0.77, p = 5 e − 04). For the modules negatively correlated with CAD&RA, the orange module (-0.71, p = 0.002), yellow module (-0.71, p = 0.002) and magenta module (-0.7, p = 0.002) had higher correlation coefficients (Fig. 2C).

The characteristic differentiated modules of CAD&RA identification

Judging from the Zsummary, the preserved modules (Zsummary ≥ 2) and differentiated modules (Zsummary < 2) for CAD&RA were identified (Fig. 3A). Finally, 13 modules were selected as the characteristic modules for CAD&RA, i.e., the blue, darkmagenta, lightcyan, lightgreen, lightsteelblue1, mediumpurple3, paleturquoise, plum1, royalblue, saddlebrown, skyblue, skyblue3 and yellowgreen module. In particular, the paleturquoise (OXSR1, ZNF141, CACNA1A, IL19) had a Zsummary less than 0 (Zsummary =-0.21), which represents the obvious differentiation between the two groups. All of the 13 selected modules are shown in Fig. 3(B-N). In addition to Zsummary, the overall expression values of 13 differentially expressed modules between CAD&RA and CAD were assessed using a t-test. We found that six modules have significant differences between the two groups (Fig. S3).

Fig. 3
figure c

The identified differentiated modules for CAD&RA. (A) Preservation analysis of defined modules using Zsummary. The x-axis represents module size; the y-axis represents the Zsummary value. Each labeled color represents a module. The dashed blue line indicates the threshold Zsummary = 2. (B-N) Networks of the 13 characteristic modules. The genes marked in yellow are hub genes

Functional and pathway enrichment analysis

GO function enrichment analysis and KEGG pathway enrichment analysis were performed in 13 differentiated modules and all DEGs respectively. In the GO functional enrichment analysis of all 278 DEGs, the top 20 terms were selected by the p value in each category. Thus, for biological processes (Fig. 4A), the genes were mainly enriched in cation homeostasis, positive regulation of cellular component movement, and nucleosome organization. Regarding the molecular functions (Fig. 4B), the genes were mainly enriched in cell adhesion molecule binding, chromatin binding and cAMP-dependent protein kinase activity. When it comes to cellular components (Fig. 4C), the genes were mainly enriched in lytic vacuole, distal axon and postsynapse. For the 13 differentiated modules, there were 14 enriched GO functions (Fig. 4D), which were mainly cAMP-dependent protein kinase activity, demethylase activity, and regulation of calcium ion import.

Fig. 4
figure d

The enriched GO terms and KEGG pathways. (A-C) Biological process, molecular function, cellular component in GO function for all 278 DEGs. (D) Enriched GO functions of 13 selected modules. On the left of each figure are the on-target numbers of the enriched genes in certain GO terms. (E) Enriched KEGG pathways of 278 DEGs

For the top 20 GO terms, some are common to multiple modules, and some are unique. For example, cAMP-dependent protein kinase activity and heart process are common terms of the blue and skyblue modules; the positive regulation of tyrosine phosphorylation of STAT protein and response to inorganic substance are common terms of yellow, plum3 and royalblue modules. In addition, the term of translation is unique to the blue module (Table. S4).

Two overlapped functions were found between the biological process of DEGs and the 13 modules, tissue migration and response to inorganic substances. Similarly, 2 enriched overlapping molecular functions of DEGs and 13 modules were found, which were cAMP-dependent protein kinase activity and demethylase activity.

In terms of KEGG pathway enrichment analysis, we found that the DEGs could enrich multiple pathways (Fig. 4E). The top three pathways were transcriptional misregulation in cancer, cell adhesion molecules, and mineral absorption. Furthermore, HPO functional enrichment analysis revealed that 13 differentially expressed modules were significantly enriched to 2 HPO terms. One of these terms was Y-linked inheritance and the other was gonosomal inheritance, which contains the genes: DDX3Y, KDM5D, CDKL5, USP9Y, ROM1, KDM6A. Moreover, among the 16 enriched pathways, 12 genes were enriched in transcriptional misregulation in cancer, accounting for 4.4%; 9 genes were enriched in the neuroactive ligand-receptor interaction pathway, accounting for 3.3%; and 8 genes were enriched in cell adhesion molecules pathway, accounting for 2.93% (Table. S5).

The identified hub genes for CAD&RA

The 13 selected modules contained 68 genes. With a GS over 0.5 and a MM over 0.8 as cut-off criteria, 49 genes were identified as hub genes (Table 1). Seven out of these 49 hub genes had a GS greater than 0.6 and MM greater than 0.9, i.e., POT1, ADO, ABCA11P, GALC, ZNF141, GPATCH8 and ATF6. Compared to the known CAD and RA-related genes, 9, 2 and 485 genes were found between CAD-related and hub genes, RA-related and hub genes, CAD-related and RA-related genes, respectively (Fig. 5A). Interestingly, the MIA3 overlapped hub genes were found to be related to both CAD and RA.

Based on their significance, the top five up-regulated genes were XIST, DEFA1, ACTA1, FAM118A and C10orf10, in which the XIST was also a hub gene. The top five down-regulated genes were EIF1AY, KDM5D, TXLNGY, RPS4Y1 and DDX3Y, and all of them were hub genes. The expression level of representative hub genes was significantly different between CAD and CAD&RA (Fig. 5C-I). The results showed that ZNF141 was down-regulated in CAD&RA, while other genes were up-regulated.

Table 1 Hub genes of the differentiated modules for CAD&RA (GS>0.5 and MM>0.8)
Fig. 5
figure e

(A) The overlapped genes among hub gene and known CAD, RA-related genes. (B) Gene expression difference significance ranking. X-axis represents the rank of DEGs, Y-axis represents log2FoldChange. (C-I) The expression level of representative hub genes between CAD and CAD&RA. (J) Classification ability of the three models based on the representative seven hub genes

Moreover, the area under the curve (AUC) values of the 7 hub genes under the three models were greater than 0.88 (Fig. 5J). Compared with the top 7 featured genes based on the homogeneity of variance test (TMX1, TCF7L2, CDC6, ZNF157, HIST3H3, COQ7 and CLDN18) and χ2 test (XIST, DDX3Y, TXLNGY, RPS4Y1, KDM5D, USP9Y and EIF1AY), our 7 susceptible genes (POT1, ADO, ABCA11P, MIA3, ZNF141, GPATCH8 and ATF6) based on GS/MM yielded optimized results in the three models, with AUC of 1.00, 1.00 and 0.88 for LR, KNN, and SVM, respectively (Table. S6), which indicated the excellent classification effect of 7 hub genes. In addition, our 7 hub genes had smaller Sab and shortest distance values than those of genes identified by other two methods, indicating the superiority of the susceptible genes identified by our modular-based analysis (Table S7).


Considering the increased mortality for the comorbid presence of CAD and RA [1,2,3], it is essential to uncover the underlying mechanisms of CAD&RA. For the complexity of CAD and RA, a network modularization approach was used to identify the characteristic module and susceptibility gene for CAD&RA. Thus, 13 modules and 49 hub genes that were related with CAD&RA were screened, and further functional enrichment analysis revealed their potential mechanisms.

Among the identified hub genes, several were reportedly related to CAD or RA. A study demonstrated that the IL19 risk allele was relevant to stroke/MI in SLE and RA, but not in the general population, showing that shared immune pathways may be contained in cardiovascular disease pathogenesis and inflammatory rheumatic diseases [36]. The expression of UTY and PRKY was found associated with the risk of CAD [37, 38]. Studies have proved that the SNP rs17465637 in the MIA3 gene was associated with the risk of CAD and RA [39, 40]. Another study found that high-intensity interval training (HIIT) could improve RA skeletal muscle gene-BCKDHB, which can increase amino acid catabolism and interconversion [41]. FYN is one of genes that is likely to play a significant role in maintenance and functioning of several of the replicated pathways of CAD [42]. Simultaneously, FYN gene is a diagnostic biomarker and one of key driver genes in RA synovial tissue subtypes C1 and C3 [43,44,45]. The ATF6 gene also plays an important role in both CAD and RA [46, 47]. Additionally, a study analyzed that AKAP13 was one the of hub genes unique to CAD [48], a finding consistent with our study. Another study showed that CYP1A2 genotype can modify the risk of RA and CYP1A2*1F allele may relate to leflunomide toxicity in RA patients therapy [49, 50]. A previous study has shown that INSL6 which produced by TNF-polarized macrophages can stimulate bone formation in mice with RA [51]. In addition, the binding of XIST to GATA1 can promote to RA [52]. The genes LIPT1 [53, 54] and TMEM40[55] were reported to be susceptibility genes with RA. Moreover, POT1 expression levels are significantly lower in RA than in the control group in vitro [56]. Mass spectrometry results revealed [57] that GALC expression levels were significantly increased in patients with atherosclerosis. In samples collected from male patients with new-onset heart failure, the RPS4Y1 was overexpression [58]. In addition, coffee intake is correlated with a risk of nonfatal myocardial infarction; this correlation is believed to be influenced by CYP1A2, which is related with the development of RA in Korea [59, 60]. After infarction, the expression of CACNA1A can enhance cardiac differentiation of brown adipose-derived stem cells to regenerate the myocardium after infarction [61]. Besides, using advanced technologies of lncRNA subcellular localization and silencing, lnc-KDM5D-4 expression was shown to be associated with atherosclerosis and CAD in men [62].

A total of 14 GO function terms were enriched by the differentiated modules, the top three terms were cAMP-dependent protein kinase activity, demethylase activity, and regulation of calcium ion import. A related study found that during cardiac preservation, a cAMP pulse could reduce the incidence and severity of transplant-related CAD [63]. Vasoactive intestinal peptide (VIP) may be an effective anti-RA treatment because it leads to the elevation of intracellular cAMP, which can inhibit TNF-α production in macrophages [64]. Another study indicates that a combination of cilostazol and MTX can activate the cAMP-dependent protein kinase pathway in the synovial fibroblasts resulting in the suppression of the inflammation of RA [65]. Fibroblast-like synoviocytes (FLSs) are involved in RA joint destruction, and pathologic process and elevated JMJD3 promotes the proliferation and migration of FLS [66]. A study found that the` Janus kinase-signal transducer and activator of transcription (JAK-STAT) pathway is an emerging target in inflammation, mainly in RA, and it heightens the cardiovascular risk [67]. Overexpression of a histone demethylase KDM4B could boost cell growth, migration and invasion, and inhibit apoptosis of FLS in RA by activating STAT3 signaling [68]. Basal intracellular calcium ion concentrations in patients with inactive RA were significantly higher than in healthy individuals, which in turn were greater than in the active RA group, which showed the important roles of calcium ions in the pathological process of RA [69].

In the top three pathways of 16 KEGG pathways, cell adhesion molecules and mineral absorption were associated with both CAD and RA. For the pathway related to cell adhesion molecules, the expression levels of cell adhesion molecules increased in patients with RA, and were associated with disease activity, oxidative stress, and inflammatory markers targeting the expression of these molecules is an important therapeutic strategy for RA [70, 71]. Moreover, the expression of both CDC42 and microRNA-34a was correlated with that of cell adhesion molecules in patients with CAD [72, 73]. For the mineral absorption pathway, clinical trials revealed that the concentrations of the mineral copper were higher in patients with RA than in healthy people [74], and zinc and selenium levels in patients with CAD admitted for coronary artery bypass grafting were reduced compared to those before surgery [75].

Although we have found several of the related modules and susceptible genes, certain limitations for our study also exist. For limited datasets and samples involving the comorbidities, CAD and RA, cross validation could not be performed. With further clinical sequencing and updated cardiovascular disease and RA-related databases, investigations should continue to validate the modular mechanism of CAD&RA. Besides, the proposed susceptible genes also need further experimental and clinical validation.

In conclusion, thirteen characteristic modules and 49 susceptible hub genes for CAD&RA were identified by network modularization analysis, including ADO, ABCA11P, GALC, ZNF141, GPATCH8, ATF6, MIA3, etc. These hub genes and their corresponding molecular functions may reflect the underlying mechanism of CAD&RA, which can provide novel perspectives for their clinical therapy strategies and precise drug discovery.

Data Availability

The raw data of GSE110008 are freely available from GEO database.


  1. Wang H, Li X, Gong G. Cardiovascular outcomes in patients with co-existing coronary artery disease and rheumatoid arthritis: a systematic review and meta-analysis. Med (Baltim). 2020;99(14):e19658.

    Article  Google Scholar 

  2. Meune C, Touzé E, Trinquart L, Allanore Y. Trends in cardiovascular mortality in patients with rheumatoid arthritis over 50 years: a systematic review and meta-analysis of cohort studies. Rheumatology (Oxford). 2009;48(10):1309–13.

    Article  PubMed  Google Scholar 

  3. Avina-Zubieta JA, Thomas J, Sadatsafavi M, Lehman AJ, Lacaille D. Risk of incident cardiovascular events in patients with rheumatoid arthritis: a meta-analysis of observational studies. Ann Rheum Dis. 2012;71(9):1524–9.

    Article  PubMed  Google Scholar 

  4. Spartera M, Godino C, Baldissera E, et al. Long-term clinical outcomes of patients with rheumatoid arthritis and concomitant coronary artery disease. Am J Cardiovasc Dis. 2017;7(1):9–18. Published 2017 Feb 15.

    PubMed  PubMed Central  Google Scholar 

  5. Kitas GD, Erb N. Tackling ischaemic heart disease in rheumatoid arthritis. Rheumatology (Oxford). 2003;42(5):607–13.

    Article  CAS  PubMed  Google Scholar 

  6. Solomon DH, Karlson EW, Rimm EB, et al. Cardiovascular morbidity and mortality in women diagnosed with rheumatoid arthritis. Circulation. 2003;107(9):1303–7.

    Article  PubMed  Google Scholar 

  7. Karpouzas GA, Malpeso J, Choi TY, Li D, Munoz S, Budoff MJ. Prevalence, extent and composition of coronary plaque in patients with rheumatoid arthritis without symptoms or prior diagnosis of coronary artery disease. Ann Rheum Dis. 2014;73(10):1797–804.

    Article  PubMed  Google Scholar 

  8. Aubry MC, Maradit-Kremers H, Reinalda MS, Crowson CS, Edwards WD, Gabriel SE. Differences in atherosclerotic coronary heart disease between subjects with and without rheumatoid arthritis. J Rheumatol. 2007;34(5):937–42.

    CAS  PubMed  Google Scholar 

  9. Skeoch S, Bruce IN. Atherosclerosis in rheumatoid arthritis: is it all about inflammation? Nat Rev Rheumatol. 2015;11(7):390–400.

    Article  CAS  PubMed  Google Scholar 

  10. Kremers HM, Crowson CS, Therneau TM, Roger VL, Gabriel SE. High ten-year risk of cardiovascular disease in newly diagnosed rheumatoid arthritis patients: a population-based cohort study. Arthritis Rheum. 2008;58(8):2268–74.

    Article  PubMed  Google Scholar 

  11. Solomon DH, Curhan GC, Rimm EB, Cannuscio CC, Karlson EW. Cardiovascular risk factors in women with and without rheumatoid arthritis. Arthritis Rheum. 2004;50(11):3444–9.

    Article  PubMed  Google Scholar 

  12. Hafström I, Rohani M, Deneberg S, Wörnert M, Jogestrand T, Frostegård J. Effects of low-dose prednisolone on endothelial function, atherosclerosis, and traditional risk factors for atherosclerosis in patients with rheumatoid arthritis–a randomized study. J Rheumatol. 2007;34(9):1810–6.

    PubMed  Google Scholar 

  13. Mudau M, Genis A, Lochner A, Strijdom H. Endothelial dysfunction: the early predictor of atherosclerosis. Cardiovasc J Afr. 2012;23(4):222–31.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Crowson CS, Liao KP, Davis JM 3rd, et al. Rheumatoid arthritis and cardiovascular disease. Am Heart J. 2013;166(4):622–628e1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Faccini A, Kaski JC, Camici PG. Coronary microvascular dysfunction in chronic inflammatory rheumatoid diseases. Eur Heart J. 2016;37(23):1799–806.

    Article  PubMed  Google Scholar 

  16. Wang L, Zhang Y, Zhang SY. Immunotherapy for the rheumatoid arthritis-associated coronary artery disease: promise and future. Chin Med J (Engl). 2019;132(24):2972–83.

    Article  CAS  PubMed  Google Scholar 

  17. Jessee R, Peart E, Beineke P, et al. Rheumatoid arthritis complicates noninvasive whole blood gene expression testing for coronary artery disease. Am Heart J. 2017;192:13–8.

    Article  CAS  PubMed  Google Scholar 

  18. Barabási AL, Gulbahce N, Loscalzo J. Network medicine: a network-based approach to human disease. Nat Rev Genet. 2011;12(1):56–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Bloom JS, Kotenko I, Sadhu MJ, Treusch S, Albert FW, Kruglyak L. Genetic interactions contribute less than additive effects to quantitative trait variation in yeast. Nat Commun. 2015;6:8712. Published 2015 Nov 5.

    Article  CAS  PubMed  Google Scholar 

  20. Menche J, Sharma A, Kitsak M, et al. Disease networks. Uncovering disease-disease relationships through the incomplete interactome. Science. 2015;347(6224):1257601.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Sarfati D, Koczwara B, Jackson C. The impact of comorbidity on cancer and its treatment. CA Cancer J Clin. 2016;66(4):337–50.

    Article  PubMed  Google Scholar 

  22. Goh KI, Cusick ME, Valle D, Childs B, Vidal M, Barabási AL. The human disease network. Proc Natl Acad Sci U S A. 2007;104(21):8685–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Beltrao P, Cagney G, Krogan NJ. Quantitative genetic interactions reveal biological modularity. Cell. 2010;141(5):739–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Sales-Pardo M. The importance of being modular. Science. 2017;357(6347):128–9.

    Article  CAS  PubMed  Google Scholar 

  25. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. Published 2008 Dec 29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Langfelder P, Luo R, Oldham MC, Horvath S. Is my Network Module Preserved and Reproducible? Plos Comput Biol 2011;7.

  27. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587–92.

    Article  CAS  PubMed  Google Scholar 

  29. Davidović R, Perovic V, Gemovic B, Veljkovic N. Bioinformatics. 2019;btz836. DiNGO: standalone application for Gene Ontology and Human Phenotype Ontology term enrichment analysis [published online ahead of print, 2019 Nov 8].

  30. Ren Y, van Blitterswijk M, Allen M, et al. TMEM106B haplotypes have distinct gene expression patterns in aged brain. Mol Neurodegener. 2018;13(1):35. Published 2018 Jul 3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Yang Q, Wang R, Wei B et al. Candidate Biomarkers and Molecular Mechanism Investigation for Glioblastoma Multiforme Utilizing WGCNA. Biomed Res Int. 2018;2018:4246703. Published 2018 Sep 26.

  32. Dreiseitl S, Ohno-Machado L. Logistic regression and artificial neural network classification models: a methodology review. J Biomed Inform. 2002;35(5–6):352–9.

    Article  PubMed  Google Scholar 

  33. Zhang S, Li X, Zong M, Zhu X, Wang R. Efficient kNN classification with different numbers of nearest neighbors. IEEE Trans Neural Netw Learn Syst. 2018;29(5):1774–85.

    Article  PubMed  Google Scholar 

  34. Keerthi SS, Shevade SK, Bhattacharyya C. Murthy; improvements to Platt’s SMO algorithm for SVM Classifier Design. Neural Comput. 2001;13(3):637–49.

    Article  Google Scholar 

  35. Wang Y, Yang H, Chen L, Jafari M, Tang J. Network-based modeling of herb combinations in traditional chinese medicine. Brief Bioinform. 2021;22(5):bbab106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Leonard D, Svenungsson E, Dahlqvist J, et al. Novel gene variants associated with cardiovascular disease in systemic lupus erythematosus and rheumatoid arthritis. Ann Rheum Dis. 2018;77(7):1063–9.

    Article  CAS  PubMed  Google Scholar 

  37. Bloomer LD, Nelson CP, Eales J, et al. Male-specific region of the Y chromosome and cardiovascular risk: phylogenetic analysis and gene expression studies. Arterioscler Thromb Vasc Biol. 2013;33(7):1722–7.

    Article  CAS  PubMed  Google Scholar 

  38. Eales JM, Maan AA, Xu X, et al. Human Y chromosome exerts Pleiotropic Effects on susceptibility to atherosclerosis. Arterioscler Thromb Vasc Biol. 2019;39(11):2386–401.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Koch W, Schatke A, Wolferstetter H, Mueller JC, Schömig A, Kastrati A. Extended evidence for association between the melanoma inhibitory activity 3 gene and myocardial infarction. Thromb Haemost. 2011;105(4):670–5.

    Article  CAS  PubMed  Google Scholar 

  40. García-Bermúdez M, López-Mejías R, González-Juanatey C, et al. Association study of MIA3 rs17465637 polymorphism with cardiovascular disease in rheumatoid arthritis patients. DNA Cell Biol. 2012;31(8):1412–7.

    Article  CAS  PubMed  Google Scholar 

  41. Andonian BJ, Johannemann A, Hubal MJ, et al. Altered skeletal muscle metabolic pathways, age, systemic inflammation, and low cardiorespiratory fitness associate with improvements in disease activity following high-intensity interval training in persons with rheumatoid arthritis. Arthritis Res Ther. 2021;23(1):187. Published 2021 Jul 10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Ghosh S, Vivar J, Nelson CP, et al. Systems Genetics analysis of genome-wide Association Study reveals Novel Associations between key biological processes and coronary artery disease. Arterioscler Thromb Vasc Biol. 2015;35(7):1712–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Liu J, Chen N. A 9 mRNAs-based diagnostic signature for rheumatoid arthritis by integrating bioinformatic analysis and machine-learning. J Orthop Surg Res. 2021;16(1):44. Published 2021 Jan 11.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Jung SM, Park KS, Kim KJ. Deep phenotyping of synovial molecular signatures by integrative systems analysis in rheumatoid arthritis. Rheumatology (Oxford). 2021;60(7):3420–31.

    Article  CAS  PubMed  Google Scholar 

  45. Brisslert M, Bian L, Svensson MN, et al. S100A4 regulates the src-tyrosine kinase dependent differentiation of Th17 cells in rheumatoid arthritis. Biochim Biophys Acta. 2014;1842(11):2049–59.

    Article  CAS  PubMed  Google Scholar 

  46. Gonzales K, Feng V, Bikkina P, Landicho MA, Haas MJ, Mooradian AD. The effect of nicotine and dextrose on endoplasmic reticulum stress in human coronary artery endothelial cells. Toxicol Res (Camb). 2021;10(2):284–91. Published 2021 Mar 23.

    Article  PubMed  Google Scholar 

  47. Lee WS, Jeong JH, Lee EG, et al. Tacrolimus regulates endoplasmic reticulum stress-mediated osteoclastogenesis and inflammation: in vitro and collagen-induced arthritis mouse model. Cell Biol Int. 2018;42(4):393–402.

    Article  CAS  PubMed  Google Scholar 

  48. Dong C, Tang L, Liu Z, et al. Landscape of the relationship between type 2 diabetes and coronary heart disease through an integrated gene network analysis. Gene. 2014;539(1):30–6.

    Article  CAS  PubMed  Google Scholar 

  49. Cornelis MC, Bae SC, Kim I, El-Sohemy A. CYP1A2 genotype and rheumatoid arthritis in Koreans. Rheumatol Int. 2010;30(10):1349–54.

    Article  CAS  PubMed  Google Scholar 

  50. Bohanec Grabar P, Rozman B, Tomsic M, Suput D, Logar D, Dolzan V. Genetic polymorphism of CYP1A2 and the toxicity of leflunomide treatment in rheumatoid arthritis patients. Eur J Clin Pharmacol. 2008;64(9):871–6.

    Article  CAS  PubMed  Google Scholar 

  51. Yi X, Liu X, Kenney HM, et al. TNF-Polarized macrophages produce insulin-like 6 peptide to stimulate bone formation in rheumatoid arthritis in mice. J Bone Miner Res. 2021;36(12):2426–39.

    Article  CAS  PubMed  Google Scholar 

  52. Yu B, Chen Y, Chen E, et al. LncRNA RNA XIST binding to GATA1 contributes to rheumatoid arthritis through its effects on proliferation of synovial fibroblasts and angiogenesis via regulation of CCN6. Mol Immunol. 2023;153:200–11.

    Article  CAS  PubMed  Google Scholar 

  53. Zhou Y, Li X, Ng L, et al. Identification of copper death-associated molecular clusters and immunological profiles in rheumatoid arthritis. Front Immunol. 2023;14:1103509. Published 2023 Feb 20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Wang W, Chen Z, Hua Y. Bioinformatics Prediction and Experimental Validation Identify a Novel Cuproptosis-Related Gene Signature in Human Synovial Inflammation during Osteoarthritis Progression. Biomolecules. 2023;13(1):127. Published 2023 Jan 7. doi:

  55. Yu X, Teng H, Marques A, Ashgari F, Ibrahim SM. High resolution mapping of Cia3: a common arthritis quantitative trait loci in different species. J Immunol. 2009;182(5):3016–23.

    Article  CAS  PubMed  Google Scholar 

  56. Qing YF, Zhou JG, Zhao MC, et al. Altered expression of TPP1 in fibroblast-like synovial cells might be involved in the pathogenesis of rheumatoid arthritis. Rheumatol Int. 2012;32(8):2503–10.

    Article  CAS  PubMed  Google Scholar 

  57. Li M, Tong X, Lv P, Feng B, Yang L, Wu Z, Cui X, Bai Y, Huang Y, Liu H. A not-stop-flow online normal-/reversed-phase two-dimensional liquid chromatography-quadrupole time-of-flight mass spectrometry method for comprehensive lipid profiling of human plasma from atherosclerosis patients. J Chromatogr A. 2014 Dec 12;1372 C:110–119. doi: 10.1016/j.chroma.2014.10.094. Epub 2014 Nov 3. PMID: 25465009.

  58. Heidecker B, Lamirault G, Kasper EK, Wittstein IS, Champion HC, Breton E, Russell SD, Hall J, Kittleson MM, Baughman KL, Hare JM. The gene expression profile of patients with new-onset heart failure reveals important gender-specific differences. Eur Heart J. 2010 May;31(10):1188–96. Epub 2009 Dec 22. PMID: 20031959; PMCID: PMC2869442.

  59. Cornelis MC, El-Sohemy A, Kabagambe EK, Campos H, Coffee. CYP1A2 genotype, and risk of myocardial infarction. JAMA. 2006;295(10):1135–41.

    Article  CAS  PubMed  Google Scholar 

  60. Cornelis MC, Bae SC, Kim I, El-Sohemy A. CYP1A2 genotype and rheumatoid arthritis in Koreans. Rheumatol Int. 2010;30(10):1349–54.

    Article  CAS  PubMed  Google Scholar 

  61. Wang H, Shi J, Wang Y, et al. Promotion of cardiac differentiation of brown adipose derived stem cells by chitosan hydrogel for repair after myocardial infarction. Biomaterials. 2014;35(13):3986–98.

    Article  CAS  PubMed  Google Scholar 

  62. Molina E, Chew GS, Myers SA et al. A Novel Y-Specific Long Non-Coding RNA Associated with Cellular Lipid Accumulation in HepG2 cells and Atherosclerosis-related Genes. Sci Rep. 2017;7(1):16710. Published 2017 Dec 1. doi:

  63. Wang CY, Aronson I, Takuma S, et al. cAMP pulse during preservation inhibits the late development of cardiac isograft and allograft vasculopathy. Circ Res. 2000;86(9):982–8.

    Article  CAS  PubMed  Google Scholar 

  64. Foey AD, Field S, Ahmed S, et al. Impact of VIP and cAMP on the regulation of TNF-alpha and IL-10 production: implications for rheumatoid arthritis. Arthritis Res Ther. 2003;5(6):R317–28.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Kim HY, Lee SW, Park SY, et al. Efficacy of concurrent administration of cilostazol and methotrexate in rheumatoid arthritis: pharmacologic and clinical significance. Life Sci. 2012;91(7–8):250–7.

    Article  CAS  PubMed  Google Scholar 

  66. Jia W, Wu W, Yang D, et al. Histone demethylase JMJD3 regulates fibroblast-like synoviocyte-mediated proliferation and joint destruction in rheumatoid arthritis. FASEB J. 2018;32(7):4031–42.

    Article  CAS  PubMed  Google Scholar 

  67. Baldini C, Moriconi FR, Galimberti S, Libby P, De Caterina R. The JAK-STAT pathway: an emerging target for cardiovascular disease in rheumatoid arthritis and myeloproliferative neoplasms. Eur Heart J. 2021;42(42):4389–400.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Zhang X, Nan H, Guo J, Liu J. KDM4B overexpression promotes the growth, Migration, and Invasion of Rheumatoid Arthritis Fibroblast-Like Synoviocytes by activating STAT3 pathway. Biochem Genet. 2021;59(6):1427–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Goulding NJ, Guyre PM. Impairment of neutrophil fc gamma receptor mediated transmembrane signalling in active rheumatoid arthritis. Ann Rheum Dis. 1992;51(5):594–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Salem HR, Zahran ES. Vascular cell adhesion molecule-1 in rheumatoid arthritis patients: relation to disease activity, oxidative stress, and systemic inflammation. Saudi Med J. 2021;42(6):620–8.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Veale DJ, Maple C. Cell adhesion molecules in rheumatoid arthritis. Drugs Aging. 1996;9(2):87–92.

    Article  CAS  PubMed  Google Scholar 

  72. Zhou M, Wu J, Tan G. The relation of circulating cell division cycle 42 expression with Th1, Th2, and Th17 cells, adhesion molecules, and biochemical indexes in coronary heart disease patients. Ir J Med Sci. 2022;191(5):2085–90.

    Article  CAS  PubMed  Google Scholar 

  73. Li H, Chen M, Feng Q, et al. MicroRNA-34a in coronary heart disease: correlation with disease risk, blood lipid, stenosis degree, inflammatory cytokines, and cell adhesion molecules. J Clin Lab Anal. 2022;36(1):e24138.

    Article  CAS  PubMed  Google Scholar 

  74. Yazar M, Sarban S, Kocyigit A, Isikan UE. Synovial fluid and plasma selenium, copper, zinc, and iron concentrations in patients with rheumatoid arthritis and osteoarthritis. Biol Trace Elem Res. 2005;106(2):123–32.

    Article  CAS  PubMed  Google Scholar 

  75. Al-Bader A, Christenson JT, Simonet F, Abul H, Dashti H, Schmuziger M. Inflammatory response and oligo-element alterations following cardiopulmonary bypass in patients undergoing coronary artery bypass grafting. Cardiovasc Surg. 1998;6(4):406–14.

    Article  CAS  PubMed  Google Scholar 

Download references


We acknowledge the support from National Natural Science Foundation of China and thank all other members of Institute of Chinese Materia Medica, China Academy of Chinese Medical Sciences for suggestions during the project.


This work was funded by the Scientific and Technological Innovation Project of China Academy of Chinese Medical Sciences (NO. CI2021A05052 to B.Li), the Fundamental Research Funds for Central public welfare research institutes (ZXKT21024 to B.Li) and the National Natural Science Foundation of China (No. 81873199 to HZ and 81904064 PW).

Author information

Authors and Affiliations



Siqi Zhang wrote the main manuscript text, Qikai Niu prepared the disease’ data, Lin Tong and Sihong Liu checked the associated data, Pengqian Wang and Haiyu Xu gave some guidance and Bing Li and Huamin Zhang guided the writing and revision of the entire article. All authors reviewed the manuscript.

Corresponding authors

Correspondence to Bing Li or Huamin Zhang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1: Table. S1 278 differentially expressed genes.

Supplementary Material 2: Table. S2 CAD-related and RA-related genes in 3 disease databases.

Supplementary Material 3: Table. S3 The relationship between gene and module.

Supplementary Material 4: Table. S4 The common and unique GO items between modules.

Supplementary Material 5: Table. S5 The KEGG pathway of 278 DEGs.

Supplementary Material 6: Table. S6 The comparison of 7 susceptible genes with other featured genes.

Supplementary Material 7: Table. S7 The proximity between genes and CAD&RA-related diseases genes


Supplementary Material 8: Overall workflow of this study. Fig. S1A The source of gene expression profile data of CAD and CAD&RA and differentially expressed genes analysis. B The overlapped genes among CAD, RA and DEGs in disease databases. C Co-expression network construction, module division and clinical traits analysis. D The identification of differentiated modules using Zsummary statistic and hub genes screening. E Functional and pathway enrichment analysis, gene expression difference significance ranking analysis and the expression level of representive hub genes


Supplementary Material 9: Fig. S2 Validation the robustness of modules with Zsummary. The x-axis represents module size; the y-axis represents the Zsummary value. Each labeled color represents a module. The dashed blue line indicates the threshold Zsummary = 2. A 1/2 samples of CAD and CAD&RA test. B 3/4 samples of CAD and CAD&RA test.


Supplementary Material 10: Fig. S3 Evaluation the 13 differentiated modules with t-test. The x-axis represents module; the y-axis represents the p value.

Supplementary Material 11: The code of this study.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, S., Niu, Q., Tong, L. et al. Identification of the susceptible genes and mechanism underlying the comorbid presence of coronary artery disease and rheumatoid arthritis: a network modularization analysis. BMC Genomics 24, 411 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: