Gene Co-expression Network Analysis Associated with Acupuncture Treatment of Rheumatoid Arthritis: An Animal Model

Article information

J Acupunct Res. 2020;37(2):128-135
1Unit of Epidemiology and Biostatistics, Department of Public Health, University of Southern Denmark, Denmark
2Protein Research Group, Department of Biochemistry and Molecular Biology, University of Southern Denmark, Odense, Denmark
3Department of Internal Medicine, University Hospital Svendborg, Svendborg, Denmark
4Unit of Human Genetics, Department of Clinical Research, University of Southern Denmark
*Corresponding author: Qihua Tan, Epidemiology & Biostatistics, Department of Public Health, Faculty of Health Science, University of Southern Denmark, Denmark, E-mail: qtan@health.sdu.dk
Received 2020 March 17; Revised 2020 May 03; Accepted 2020 May 07.

Abstract

Background

Classical acupuncture is being used in the treatment of rheumatoid arthritis (RA). To explore the biological response to acupuncture, a network-based analysis was performed on gene expression data collected from an animal model of RA treated with acupuncture.

Methods

Gene expression data were obtained from published microarray studies on blood samples from rats with collagen induced arthritis (CIA) and non-CIA rats, both treated with manual acupuncture. The weighted gene co-expression network analysis was performed to identify gene clusters expressed in association with acupuncture treatment time and RA status. Gene ontology and pathway analyses were applied for functional annotation and network visualization.

Results

A cluster of 347 genes were identified that differentially downregulated expression in association with acupuncture treatment over time; specifically in rats with CIA with module-RA correlation at 1 hour after acupuncture (−0.27; p < 0.001) and at 34 days after acupuncture (−0.33; p < 0.001). Functional annotation showed highly significant enrichment of porphyrin-containing compound biosynthetic processes (p < 0.001). The network-based analysis also identified a module of 140 genes differentially expressed between CIA and non-CIA in rats (p < 0.001). This cluster of genes was enriched for antigen processing and presentation of exogenous peptide antigen (p < 0.001). Other functional gene clusters previously reported in earlier studies were also observed.

Conclusion

The identified gene expression networks and their hub-genes could help with the understanding of mechanisms involved in the pathogenesis of RA, as well understanding the effects of acupuncture treatment of RA.

Introduction

Rheumatoid arthritis (RA) is an autoimmune inflammatory joint disease/disorder which affects the global population. Without treatment, RA can lead to joint destruction and severe disability. Although disease-modifying anti-rheumatic drugs and biological treatments are available for managing the condition, there is currently no cure for RA. Acupuncture is now widely accepted as a treatment for RA patients [1], however, knowledge about the underlying mechanism of the effect of acupuncture on RA is sparse.

Gene expression microarray is a platform for investigating gene activities under different physiological conditions e.g. disease status, treatment, and patient prognosis. Using an animal model of RA, Nardini et al [2] analyzed single gene expression patterns in association with acupuncture treatment of RA in rats and reported a reduced inflammatory response. In this study, a network-based analysis was used to detect and examine gene cluster profiles that are co-expressed during acupuncture treatment.

Materials and Methods

Samples

The Gene Expression Omnibus database (accession number GSE48025) was used in this study. The data for gene expression microarray was derived from 6 female Wistar rats with collagen induced arthritis (CIA) and a group of 6 female control Wistar rats. All rats were anesthetized with ether and shaved on acupunctural areas. The shaving was performed on all animals to minimize bias. Whole blood samples (5 mL) were extracted prior to acupuncture treatment in 2 CIA and 2 control rats (Time 0). Acupuncture was performed on the remaining animals by punctuation at 4 points; bilateral on the lower spine, between 2nd and 3rd lumbar vertebra (BL 23) and between the tibia and fibula on the hind legs (ST 36). The acupuncture was performed by inserting a thin needle and twirling for 20 seconds clockwise every other day. The insertion points used were similar to earlier studies on classical Chinese medicine as described by Nardini et al [2].

At the Time 1 timepoint (1 hour after acupuncture treatment), whole blood samples (5 mL) were extracted from 2 CIA rats and 2 control rats. At the Time 2 timepoint (34 days after acupuncture treatment), the 2 remaining CIA rats, and the 2 remaining control rats were anesthetized, then whole blood samples were drawn. All 12 rats (4 rats each from Time 0 to Time 2), were sacrificed after blood sampling at the end of the study.

The whole blood samples were analyzed using Affymetrix rat230 2.0 arrays which comprised of more than 31,000 probe sets, analyzing over 30,000 transcripts and variants from over 28,000 well-substantiated rat genes, which were hybridized and scanned according to the manufacturer’s recommendations. The output files were analyzed using the mas5 function from the R package affy to generate expression intensities for each probe. Details on raw microarray data pre-processing (QC and normalization) are described in Nardini et al [2].

Gene co-expression network analysis

Since multiple probes were mapped to the same gene, the expression for each gene was first summarized by calculating the mean. This resulted in 15,397 unique genes. The WGCNA R package provides functions for weighted correlation network analysis which were applied to identify the co-expression modules. Fig. 1 is the workflow for WGCNA where it is assumed that the genetic network in the system is scale free, converting their adjacency matrix into a matrix of connection of strength by raising it to a power of beta, which is chosen to satisfy the assumption of a scale free network.

Fig. 1

The workflow of WGCNA.

In this study, the power of 10 was chosen based on a scale-free connectivity plot (Fig. 2). Pearson correlation was used for statistical analysis by assuming a linear relationship between gene expression profiles, as well as between gene expression and experimental variables. To identify the modules, the adjacency matrix was used to construct the Topology Overlap Matrix (TOM) and from this the dissimilarity matrix 1-TOM was constructed to perform hierarchical clustering. This returns a number of modules corresponding to clusters in other methods, indicated by colors (Fig. 3). Modules highly correlated by their eigengenes (1st principal component of all gene expression in a module; > 0.75) were merged to reduce the number of clusters.

Fig. 2

Selection of soft thresholding power. The left panel shows the scale free fit index on the y-axis as a function of the soft-thresholding power on the x-axis. The blue line is placed at 0.80. The right panel shows the mean connectivity with degree on the y-axis and soft-thresholding power in the x-axis.

Fig. 3

A cluster dendrogram showing the clustering of genes with colors assigned to detected gene modules and the merged gene modules (bottom).

Linear regression was used to explore interactions between early and late responses to treatment of RA. Each trait was focused on one by one (disease/condition: CIA = 1, control = 0; acupuncture time) and the absolute correlation between a module eigengene and the specific trait was established. The association of the individual genes to each module was then quantified.

In WGCNA, the gene significance (GS) is defined as: GSi = |cor (xi, Clinical Trait)|

Here xi is the expression of gene i. As we perform a linear regression, GS was defined as: GSi = −log10 (p value of Trait for gene i). The interpretation of GS is that the bigger the GSi, the more biologically significant the gene. As a measurement of relative importance for individual genes in a module, a module-membership (MM) was defined and calculated MM(i) = cor (xi, ME). The interpretation of this function was that the higher the absolute value of MM(i), the more important gene i was to the module. A high MM would indicate that the gene was a possible hub-gene in the module, and therefore, is of special interest. By examining the relationship between GS and MM, important hub-genes with high gene significance for the specific traits, as well as high module-membership can be identified.

Functional annotation

To investigate the biological processes for the genes in the modules associated significantly with traits, gene ontology (GO)5 and KEGG by the R package ClusterProfiler [3] was used. For the GO results, the command EnrichGO was used. All 3 ontologies were chosen to achieve a better overview of the different processes, and a p value cut off at 0.05 was selected. For better understanding of the biological processes, the GO results were submitted to AmiGO [4], which is an open access database for ontology and annotation data. For the KEGG analysis, enrich KEGG command was used, with the p value cut-off at 0.05. For further information about the pathways, KEGG results can be viewed on GenomeNet (www.genome.jp) in the KEGG PATHWAY database [5]. Hub-genes are deducted from GO results and they are visualized using the CNet plot in ClusterProfiler, the lowest p value shown in blue. Selected possible hub-genes are inputted online to retrieve possible associations with RA.

Results

RA associated gene expression networks

A cluster of 140 genes were significantly differentially expressed in the disease/condition trait (p < 0.001). The main biological process of the GO terms (Table 1) for this module was antigen processing and presentation of exogenous peptide antigen (r = −0.26, p < 0.001). In the molecular function part of the GO results, the top result was lytic vacuole and lysosome, both with p < 0.001. The cellular component GO terms showed activation in the lysosome, endocytic vesicles, and vacuoles. The top KEGG pathways (Table 2) were lysosomal, cAMP signaling, salivary secretion, and adrenergic signaling in cardiomyocytes, all with p < 0.001. Fig. 4 illustrates the gene expression networks and their linked genes in relation to CIA.

Gene Ontology Dark Orange Module.

KEGG Pathway Results Dark Orange Module.

Fig. 4

CNet-plot from the GO results for differential expression between CIA and non-CIA rats. Biological processes are shown with beige dots, cellular components are shown with orange dots, and molecular functions are shown with green dots.

CIA, collagen induced arthritis.

Time-dependent RA associated gene expression networks in acupuncture

WGCNA identified a module consisting of 347 genes as the best representation for the interaction between acupuncture time and disease/condition, or in other words, CIA-specific time effect of acupuncture treatment. The module had a highly significant marginal association with CIA (correlation −0.44 and p < 0.001) but also displays a time-dependent pattern with a module-RA correlation at Time 1 (−0.27; p < 0.001) and at Time 2 (−0.33; p < 0.001). The top GO (Table 3) biological processes included the hemoglobin system2 involving porphyrin-containing compound biosynthetic processes, tetrapyrrole biosynthetic processes, porphyrin-containing compound metabolic processes, and heme biosynthetic processes, all with p < 0.001. The GO results for the molecular functions were peroxidase activity, oxidoreductase activity, acting on peroxide as acceptor, oxygen carrier activity and antioxidant activity, all with p < 0.001.

Gene Ontology Purple Module.

The GO results were consistent with the KEGG pathway analysis results (Table 4), as porphyrin and chlorophyll metabolism (p < 0.001) were the most significant pathway. The remaining KEGG results were malaria, collecting duct acid secretion and ferroptosis, all with p < 0.001.

KEGG Pathway Results for Purple Module.

In Fig. 5, there were many possible hub-genes of interest for the interaction between time and disease/condition or treatment effect. The most significant hub-genes for the interaction include Slc25a39, Alas2 and Urod related to the biological processes, Hba-a1, Hbe2, Cat, Prxl2a and Prdx2 related to the molecular function, and genes Prdx2, Hba-a1, Hbe2 and Cat related to both biological processes and molecular functions.

Fig. 5

CNet-plot from the GO result for differential expression over acupuncture time in CIA rats. Biological processes are shown with beige dots and molecular functions are shown with green dots.

CIA, collagen induced arthritis.

Discussion

Network-based analysis identified a gene cluster significantly associated with RA disease/condition status. In the gene cluster, antigen presenting cells (APC) were highly implicated in the biological process part of the GO-terms correlated with RA [6]. The APC in the disease/condition pattern may come from the CIA process. It is unknown if the presence of APC is due to collagen induction of arthritis or T cell activation initiating RA and involved in the joint destruction [7]. In the cellular component part of the GO results, an activation in lysosome, endocytic vesicles and vacuoles was observed which was connected to several genes that have high significance and which could be possible hub-genes. The GO results are consistent with the KEGG results, where the top pathway is lysosomal, which is an organelle in animal cells involved in degradation of macromolecules requiring a series of processes such as endocytosis, phagocytosis, and autophagy. The molecular function part of the GO results suggested an increased activity in the transport through membranes by calcium gated channels, which correlates with the increased activity by the lysosome in various cellular components.

Although every one of the significant genes in the results could be a hub-gene to focus on, a few genes were selected. The genes Unc93b1, Ifi30 and LOC100909593 are connected to both biological processes and cellular components. The genes Anxa11 and Anxa2 are related to molecular processes and cellular components, and the genes Clec4a3 and Clec4a were related to both molecular functions and biological processes.

The Unc93b1 gene encodes a protein placed in the endoplasmic reticulum and is connected to the control of trafficking toll-like receptor involved in autoimmune diseases [8]. Listed in the rat genome database the gene LOC100909593 is a Major Histocompatibility Complex (MHC) Class II histocompatibility antigen, M alpha chain-like [9]. Its encoded protein is predicted to be part of antigen processing and presentation of antigens via MHC Class II and positive regulation of the immune response [10]. MCH Class II has a strong association with the development of CIA in animal models of RA, where MCH Class II presents the immunodominant Collagen II 256-270 peptide to specific T cells which mediate the inflammatory response [11]. This corresponds well with the module’s relation to the disease/condition status. The gene Ifi30 is also called gamma-interferon-inducible lysosomal thiol reductase. It is expressed in APC and is associated with MHC Class II by contributing to the presentation of peptides loaded on the MHC [12]. Anxa11 has been associated with the inflammatory disease sarcoidosis, which is a chronic inflammatory disorder dominantly manifested in the lungs [13]. Interestingly, the disease has been shown to co-occur with RA [14].

The interaction between time and disease i.e. the effect of acupuncture over time is conditional on or specific to the disease. This interaction effect on gene expression was of high interest in this study because it reflected changes in gene activities associated with acupuncture treatment over time, specifically in the disease group. The module significant for this interaction consisted of 347 genes. The top GO biological processes were linked to the hemoglobin system (Table 3). Among the top GO processes associated with porphyrin, the component in the red blood cells which contains the heme group and involved in heme biosynthesis. This could indicate blood formation and improved oxygen transport in the body. On the other hand, this could also be an effect of the CIA-treatment which is causing increased blood flow in the skin in the CIA rats [2]. The GO results are consistent with the KEGG pathway analysis results, where porphyrin metabolism is the most significant KEGG pathway. Porphyrin derivatives have also been suggested as a photodynamic treatment where the porphyrin derivatives are used as a photosensitizer together with oxygen molecules (which produces a single oxygen and possible reactive oxygen species) and external light to interfere with cellular signaling pathways, This leads to necrosis or apoptosis in the target tissues reducing inflammation of the tissue [15].

In the KEGG results, ferroptosis, a regulated form of cell death that is present in some neurodegenerative diseases and tissue damage, may be associated with RA. Ferroptosis can be induced therapeutically by sulfasalazine, which is a drug used to treat RA [16], but as the CIA rats were not given that drug, the acupuncture treatment may have activated some of the same pathways induced by the drug. As a reaction of the metabolism of porphyrin, it could be a very interesting target for RA researchers [17]. Ferroptosis requires iron, while the latter is associated with activated hemes. Ferroptosis is the accumulation of reactive oxygen species [16], which could be interacting with the porphyrin and be activated by an external light source. This could lead to a decrease in inflammation, but further studies need to be conducted for verification.

This study was limited by the sample size and design. As an association study in an animal model, the results presented here represent correlations rather than causal relationships. Well-defined randomized controlled trials are needed to establish deterministic effects of the treatment.

Conclusion

This network-based analysis detected significant gene clusters enriched for functional pathways and biological processes implicated in an RA model in association with acupuncture treatment. The identified gene expression networks and their hub-genes could help with the understanding of mechanisms involved in the pathogenesis of RA, as well as in the effects of acupuncture treatment of RA.

Acknowledgments

This project was supported by the Region of Southern Denmark Research Grant for Alternative Traditional Chinese Medicine Acupuncture #13/26014.

Notes

Conflict of Interests

The authors have no conflicts of interest to declare.

References

1. Chou PC, Chu HY. Clinical Efficacy of Acupuncture on Rheumatoid Arthritis and Associated Mechanisms: A Systemic Review. Evid Based Complement Alternat Med 2018;20188596918.
2. Nardini C, Devescovi V, Liu Y, Zhou X, Lu Y, Dent JE. Systemic Wound Healing Associated with local sub-Cutaneous Mechanical Stimulation. Sci Rep 2016;6:39043.
3. Yu G, Wang LG, Han Y, He QY. clusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS 2012;16:284–287.
4. Carbon S, Ireland A, Mungall CJ, Shu S, Marshall B, Lewis S, et al. AmiGO : online access to ontology and annotation data. Bioinformatics 2009;25:288–289.
5. Kanehisa M, Goto S, Furumichi M, Tanabe M. KEGG for representation and analysis of molecular networks involving diseases and drugs. Nucleic Acids Res 2010;38(Database issue):D355–360.
6. Aarvak T, Natvig JB. Cell-cell interactions in synovitis: Antigen presenting cells and T cell interaction in rheumatoid arthritis. Arthritis Res 2001;3:13–17.
7. Brand DD, Kang AH, Rosloniec EF. Immunopathogenesis of collagen arthritis. Springer Semin Immunopathol 2003;25:3–18.
8. Fukui R, Saitoh S, Kanno A, Onji M, Shibata T, Ito A, et al. Unc93B1 restricts systemic lethal inflammation by orchestrating toll-like receptor 7 and 9 trafficking. Immunity 2011;35:69–81.
9. Shimoyama M, De Pons J, Hayman GT, Laulederkind SJ, Liu W, Nigam R, et al. The Rat Genome Database 2015: genomic, phenotypic and environmental variations and disease. Nucleic Acids Res 2015;43(Database issue):D743–750.
10. Laulederkind SJF, Hayman GT, Wang SJ, Lowry TF, Nigam R, Petri V, et al. Exploring Genetic, Genomic, and Phenotypic Data at the Rat Genome Database. Curr Protoc Bioinformatics 2012;Chapter 1(Unit1.14)
11. Kjellén P, Brunsberg U, Broddefalk J, Hansen B, Vestberg M, Ivarsson I, et al. The structural basis of MHC control of collagen-induced arthritis; binding of the immunodominant type II collagen 256 – 270 glycopeptide to H-2Aq and H-2Ap molecules. Eur J Immunol 1998;28:755–766.
12. Jun-Ichi S, Kino Y, Yanaizu M, Ishida T, Saito Y. Microglia express gamma-interferon-inducible lysosomal thiol reductase in the brains of Alzheimer’s disease and Nasu-Hakola disease. Intractable Rare Dis Res 2018;7:251–257.
13. Hofmann S, Franke A, Fischer A, Jacobs G, Nothnagel M, Gaede KI, et al. Genome-wide association study identifies ANXA11 as a new susceptibility locus for sarcoidosis. Nat Genet 2008;40:1103–1106.
14. Kobak Ş, Karaarslan AA, Yilmaz H, Sever F. Co-occurrence of rheumatoid arthritis and sarcoidosis. BMJ Case Rep 2015;2015bcr2014208803.
15. Zhao C, Ur Rehman F, Yang Y, Zhao C, Rehman FU, Yang Y, et al. Bio-imaging and Photodynamic Therapy with Tetra Sulphonatophenyl Porphyrin (TSPP)-TiO2 Nanowhiskers: New Approaches in Rheumatoid Arthritis Theranostics. Sci Rep 2015;5:11518.
16. Dixon SJ, Stockwell BR. The Hallmarks of Ferroptosis. Annu Rev Cancer Biol 2019;3:35–54.
17. Blouin JM1, Duchartre Y, Costet P, Lalanne M, Ged C, Lain A, et al. Therapeutic potential of proteasome inhibitors in congenital erythropoietic porphyria. Proc Natl Acad Sci U S A 2013;110:18238–18243.

Article information Continued

Fig. 1

The workflow of WGCNA.

Fig. 2

Selection of soft thresholding power. The left panel shows the scale free fit index on the y-axis as a function of the soft-thresholding power on the x-axis. The blue line is placed at 0.80. The right panel shows the mean connectivity with degree on the y-axis and soft-thresholding power in the x-axis.

Fig. 3

A cluster dendrogram showing the clustering of genes with colors assigned to detected gene modules and the merged gene modules (bottom).

Fig. 4

CNet-plot from the GO results for differential expression between CIA and non-CIA rats. Biological processes are shown with beige dots, cellular components are shown with orange dots, and molecular functions are shown with green dots.

CIA, collagen induced arthritis.

Fig. 5

CNet-plot from the GO result for differential expression over acupuncture time in CIA rats. Biological processes are shown with beige dots and molecular functions are shown with green dots.

CIA, collagen induced arthritis.

Table 1

Gene Ontology Dark Orange Module.

Ontology ID Description Gene Ratio Bg Ratio p p adjusted q Gene ID Count
BP GO:0002478 Antigen processing and presentation of exogenous peptide antigen 9/839 19/12332 1,49E-06 7,67E-03 7,56E-03 LOC100909593/B2m/Clec4a/Ifi30/Tapbp/Fcgr1a/Fcgr2b/Clec4a3/Unc93b1 9
BP GO:0019884 Antigen processing and presentation of exogenous antigen 9/839 24/12332 1,54E-05 3,98E-02 3,92E-02 LOC100909593/B2m/Clec4a/Ifi30/Tapbp/Fcgr1a/Fcgr2b/Clec4a3/Unc93b1 9
CC GO:0000323 Lytic vacuole 53/874 346/12790 2,12E-08 1,24E-05 1,18E-05 March3/LOC100909593/Slc15a3/Anxa11/Slc38a9/Vps11/Ccdc115/Hgsnat/March8/RGD1359108/Lipa/Mmp13/Milr1/Acp2/Cxcr4/Siae/Litaf/Naga/Anxa2/Unc13d/Psap/Hck/Tmem106b/Ifi30/Stx3/Ggh/Fuca1/P2ry2/Ctsf/Rab38/Nsg2/Hgs/Vamp8/Pcyox1/Ctss/Chmp2b/Bri3/Pip4p1/Slc35f6/Slc39a11/Rnaset2/Ctbs/Plekhf1/Rab7b/Gm2a/Lgmn/Asah1/Grn/Unc93b1/Vps26a/Ctsz/Plbd1/Cd68 53
CC GO:0005764 Lysosome 52/874 345/12790 4,96E-08 1,45E-05 1,38E-05 March3/LOC100909593/Slc15a3/Anxa11/Slc38a9/Vps11/Ccdc115/Hgsnat/March8/RGD1359108/Lipa/Mmp13/Milr1/Acp2/Cxcr4/Siae/Litaf/Naga/Anxa2/Unc13d/Psap/Hck/Tmem106b/Ifi30/Stx3/Ggh/Fuca1/P2ry2/Ctsf/Rab38/Nsg2/Hgs/Vamp8/Pcyox1/Ctss/Chmp2b/Bri3/Pip4p1/Slc35f6/Rnaset2/Ctbs/Plekhf1/Rab7b/Gm2a/Lgmn/Asah1/Grn/Unc93b1/Vps26a/Ctsz/Plbd1/Cd68 52
CC GO:0005773 Vacuole 60/874 431/12790 8,78E-08 1,71E-05 1,63E-05 Becn1/March3/LOC100909593/Slc15a3/Anxa11/Slc38a9/Vps11/Ccdc115/Hgsnat/Abcc3/March8/RGD1359108/Lipa/Mmp13/Milr1/Acp2/Cxcr4/Siae/Pip4k2b/Litaf/Naga/Anxa2/Unc13d/Psap/Hck/Tmem106b/Mapk15/Ifi30/Stx3/Ggh/Fuca1/P2ry2/Ctsf/Rab38/Slc22a17/Nsg2/Hgs/Vamp8/Pcyox1/Ctss/Chmp2b/Bri3/Pip4p1/Slc35f6/Slc39a11/Rnaset2/Ctbs/Vac14/Plekhf1/Rab7b/Gm2a/Lgmn/Asah1/Grn/Unc93b1/Vps26a/Ctsz/Plbd1/Cd68/Thbd 60
CC GO:0030139 Endocytic vesicle 23/874 152/12790 2,56E-04 3,23E-02 3,08E-02 Becn1/Anxa11/Vps11/Rab5c/Igf2r/Tirap/Cd2ap/Ehd3/Clta/Myo6/Rab38/Capg/Amot/Vamp8/Cybb/Pld4/Vim/Heatr5a/Rab7b/Rab11fip5/Fcho2/Unc93b1/Pld1 23
CC GO:0016323 Basolateral plasma membrane 30/874 223/12790 2,76E-04 3,23E-02 3,08E-02 Ctnna2/Slc22a5/Slc12a6/Cacna1d/Cldn4/Flot1/Abcc3/Tgfbr1/Heph/Atp1a1/Cadm1/Anxa2/C5ar1/Slc22a2/Tacstd2/Atp2b4/Atp2b1/Kcnn4/Cd38/P2ry2/Entpd1/Atp7a/Vamp8/Slc7a7/Lrp1/Clcnka/Gm2a/St14/Fxyd2/Slc8a1 30
MF GO:0005388 Calcium-transporting ATPase activity 6/795 10/11622 1,66E-05 7,66E-03 7,45E-03 Atp2b4/Atp2b1/Atp2a3/Atp2a1/Atp2b3/Anxa5 6
MF GO:0015662 ATPase activity, coupled to transmembrane movement of ions, phosphorylative mechanism 10/795 30/11622 1,81E-05 7,66E-03 7,45E-03 Atp1b3/Atp1a1/Atp6v1e1/Atp2b4/Atp2b1/Atp2a3/Atp2a1/Atp7a/Atp2b3/Anxa5 10
MF GO:0005509 Calcium ion binding 56/795 497/11622 1,38E-04 3,90E-02 3,80E-02 Myl12a/Dnah7/Stim1/S100vp/Clec3b/Syt10/Crtac1/Anxa11/Man1a2/Unc13c/Nucb1/Hpcal1/Casq1/Necab1/Micu1/Pdcd6/Tbc1d8/Ehd3/Stim2/Mmp13/Plcb1/Mmp12/Plcd1/S100a10/Myl7/Sgca/Notch2/Anxa2/Efhd1/Fat1/Clec4a/Cdhr1/S100a11/Atp2a1/Tpd52/Cetn4/Melk/Dgkg/Clec4a3/Pcdha12/Pcdh20/Lrp1/Tbc1d9/Fam20c/Adgre1/Anxa5/Vldlr/S100a6/Cabp4/Pla2g4a/Itgb1bp2/Sulf2/Anxa4/Aif1/Thbd/Slc8a1 56

Table 2

KEGG Pathway Results Dark Orange Module.

ID Description p p adjusted q Count
rno04142 Lysosome 4,12E+07 1,19E-04 1,00E-04 21
rno04970 Salivary secretion 5,26E+09 7,60E-03 6,40E-03 13
rno04972 Pancreatic secretion 2,15E-04 1,67E-02 1,41E-02 14
rno04024 cAMP signaling pathway 2,79E-04 1,67E-02 1,41E-02 22
rno04380 Osteoclast differentiation 2,89E-04 1,67E-02 1,41E-02 16
rno04071 Sphingolipid signaling pathway 6,34E-04 2,37E-02 1,99E-02 15
rno04925 Aldosterone synthesis and secretion 6,73E-04 2,37E-02 1,99E-02 13
rno04060 Cytokine-cytokine receptor interaction 6,85E-04 2,37E-02 1,99E-02 25
rno05152 Tuberculosis 9,13E-04 2,37E-02 1,99E-02 19
rno04022 cGMP-PKG signaling pathway 9,95E-04 2,37E-02 1,99E-02 18
rno04261 Adrenergic signaling in cardiomyocytes 1,07E-03 2,37E-02 1,99E-02 16
rno04978 Mineral absorption 1,09E-03 2,37E-02 1,99E-02 8
rno04961 Endocrine and other factor-regulated calcium reabsorption 1,12E-03 2,37E-02 1,99E-02 9
rno04350 TGF-beta signaling pathway 1,15E-03 2,37E-02 1,99E-02 12
rno04020 Calcium signaling pathway 1,71E-03 3,29E-02 2,77E-02 19
rno04911 Insulin secretion 2,15E-03 3,79E-02 3,19E-02 11
rno05140 Leishmaniasis 2,27E-03 3,79E-02 3,19E-02 10
rno05410 Hypertrophic cardiomyopathy (HCM) 2,36E-03 3,79E-02 3,19E-02 11
rno04933 AGE-RAGE signaling pathway in diabetic complications 2,78E-03 4,23E-02 3,56E-02 12
rno05414 Dilated cardiomyopathy (DCM) 3,36E-03 4,86E-02 4,09E-02 11

Table 3

Gene Ontology Purple Module.

Ontology ID Description Gene ratio Bg ratio p p adjusted q Gene ID Count
BP GO:0006779 Porphyrin-containing compound biosynthetic process 11/651 22/12332 3,38E-09 1,43E-05 1,43E-05 Slc25a39/Alas2/Urod/Fech/Spta1/Ank1/Alad/Hmbs/Ppox/Fxn/Uros 11
BP GO:0033014 Tetrapyrrole biosynthetic process 11/651 23/12332 6,17E-09 1,43E-05 1,43E-05 Slc25a39/Alas2/Urod/Fech/Spta1/Ank1/Alad/Hmbs/Ppox/Fxn/Uros 11
BP GO:0006778 Porphyrin-containing compound metabolic process 12/651 32/12332 3,61E-08 5,58E-05 5,58E-05 Slc25a39/Alas2/Urod/Fech/Spta1/Ank1/Alad/Eif2ak1/Hmbs/Ppox/Fxn/Uros 12
BP GO:0033013 Tetrapyrrole metabolic process 13/651 40/12332 7,04E-08 8,17E-05 8,17E-05 Slc25a39/Alas2/Urod/Fech/Spta1/Ank1/Alad/Eif2ak1/Hmbs/Ppox/Fxn/Mtr/Uros 13
BP GO:0006783 Heme biosynthetic process 9/651 18/12332 9,54E-08 8,86E-05 8,86E-05 Slc25a39/Alas2/Urod/Fech/Alad/Hmbs/Ppox/Fxn/Uros 9
BP GO:0046501 Protoporphyrinogen IX metabolic process 7/651 10/12332 1,16E-07 8,95E-05 8,95E-05 Alas2/Urod/Fech/Alad/Eif2ak1/Hmbs/Ppox 7
BP GO:0042744 Hydrogen peroxide catabolic process 10/651 24/12332 1,57E-07 1,04E-04 1,04E-04 Snca/Prdx2/Cat/Hbe2/Hba-a1/Gpx1/Hbg1/Hbb/Hbe1/Prdx5 10
BP GO:0051186 Cofactor metabolic process 50/651 438/12332 1,93E-07 1,12E-04 1,12E-04 Hagh/Gclm/Snca/Prdx2/Slc25a39/Alas2/Bpgm/Urod/Fech/Spta1/Ank1/Cat/Gstt3/Pck2/Alad/Chac2/Gsta4/Eif2ak1/Hmbs/Acot1/Hbe2/Hba-a1/Acsl1/Gstm3l/Isca1/Tyms/Mlxipl/Dlst/Gstt2/Aspdh/Gstm7/Elovl1/Gch1/Gpx1/Ghr/Ppox/Hbg1/Fxn/Nadk2/Mtr/Hbb/Iscu/Gpx4/Acmsd/Hbe1/Pdss2/Prdx5/Ldhb/Uros/Mif 50
BP GO:0015669 Gas transport 8/651 16/12332 5,10E-07 2,63E-04 2,63E-04 Car2/Aqp1/Hbe2/Hba-a1/LOC100134871/Hbg1/Hbb/Hbe1 8
BP GO:0042168 Heme metabolic process 9/651 25/12332 2,87E-06 1,33E-03 1,33E-03 Slc25a39/Alas2/Urod/Fech/Alad/Hmbs/Ppox/Fxn/Uros 9
BP GO:0051187 Cofactor catabolic process 12/651 47/12332 4,02E-06 1,70E-03 1,70E-03 Snca/Prdx2/Cat/Chac2/Hbe2/Hba-a1/Aspdh/Gpx1/Hbg1/Hbb/Hbe1/Prdx5 12
BP GO:0015671 Oxygen transport 6/651 11/12332 7,78E-06 3,01E-03 3,01E-03 Hbe2/Hba-a1/LOC100134871/Hbg1/Hbb/Hbe1 6
BP GO:0098869 Cellular oxidant detoxification 16/651 88/12332 1,25E-05 4,46E-03 4,46E-03 Prdx2/Prxl2a/Cat/Hbe2/Hba-a1/LOC100134871/Gstt2/Fbln5/Gstm7/Gch1/Gpx1/Hbg1/Hbb/Gpx4/Hbe1/Prdx5 16
BP GO:0098754 Detoxification 17/651 100/12332 1,70E-05 5,65E-03 5,65E-03 Prdx2/Akr7a3/Prxl2a/Cat/Hbe2/Hba-a1/LOC100134871/Gstt2/Fbln5/Gstm7/Gch1/Gpx1/Hbg1/Hbb/Gpx4/Hbe1/Prdx5 17
BP GO:1990748 Cellular detoxification 16/651 91/12332 1,94E-05 5,80E-03 5,80E-03 Prdx2/Prxl2a/Cat/Hbe2/Hba-a1/LOC100134871/Gstt2/Fbln5/Gstm7/Gch1/Gpx1/Hbg1/Hbb/Gpx4/Hbe1/Prdx5 16
BP GO:0017001 Antibiotic catabolic process 11/651 46/12332 2,00E-05 5,80E-03 5,80E-03 Snca/Prdx2/Cat/Pck2/Hbe2/Hba-a1/Gpx1/Hbg1/Hbb/Hbe1/Prdx5 11
BP GO:0048821 Erythrocyte development 9/651 31/12332 2,13E-05 5,81E-03 5,81E-03 Rhd/Bpgm/Sox6/Ank1/Klf1/Ercc2/Tmod3/Dmtn/Hbb 9
BP GO:0006749 Glutathione metabolic process 11/651 48/12332 3,07E-05 7,93E-03 7,93E-03 Hagh/Gclm/Gstt3/Chac2/Gsta4/Gstm3l/Gstt2/Gstm7/Gpx1/Hbb/Gpx4 11
BP GO:0042743 Hydrogen peroxide metabolic process 10/651 44/12332 7,60E-05 1,86E-02 1,86E-02 Snca/Prdx2/Cat/Hbe2/Hba-a1/Gpx1/Hbg1/Hbb/Hbe1/Prdx5 10
BP GO:0020027 Hemoglobin metabolic process 6/651 16/12332 1,07E-04 2,45E-02 2,45E-02 Alas2/Epb42/Snx3/Fech/Cat/Eif2ak1 6
BP GO:0055072 Iron ion homeostasis 12/651 64/12332 1,11E-04 2,45E-02 2,45E-02 Rhd/Alas2/Epb42/Ncoa4/Fech/Ank1/Eif2ak1/Slc25a37/Atp6v1g1/Fxn/Iscu/Cp 12
BP GO:0006790 Sulfur compound metabolic process 27/651 235/12332 1,17E-04 2,46E-02 2,46E-02 Hagh/Gclm/Snca/Mpst/Gstt3/Chac2/Gsta4/Acot1/Acsl1/Gstm3l/Isca1/Dlst/Gstt2/Gstm7/Elovl1/Gpx1/Ghr/Mical2/Fxn/Cbs/Mtr/Hbb/Iscu/Gpx4/Spock3/Idua/Mamdc2 27
BP GO:0046148 Pigment biosynthetic process 9/651 41/12332 2,30E-04 4,65E-02 4,65E-02 Slc25a39/Alas2/Urod/Fech/Alad/Hmbs/Ppox/Fxn/Uros 9
MF GO:0004601 Peroxidase activity 12/604 44/11622 1,59E-06 1,24E-03 1,21E-03 Prdx2/Cat/Hbe2/Hba-a1/Gstt2/Gstm7/Gpx1/Hbg1/Hbb/Gpx4/Hbe1/Prdx5 12
MF GO:0016684 Oxidoreductase activity, acting on peroxide as acceptor 12/604 48/11622 4,33E-06 1,69E-03 1,65E-03 Prdx2/Cat/Hbe2/Hba-a1/Gstt2/Gstm7/Gpx1/Hbg1/Hbb/Gpx4/Hbe1/Prdx5 12
MF GO:0005344 Oxygen carrier activity 5/604 10/11622 7,55E-05 1,97E-02 1,92E-02 Hbe2/Hba-a1/Hbg1/Hbb/Hbe1 5
MF GO:0016209 Antioxidant activity 13/604 76/11622 1,32E-04 2,57E-02 2,51E-02 Prdx2/Prxl2a/Cat/Hbe2/Hba-a1/Gstt2/Gstm7/Gpx1/Hbg1/Hbb/Gpx4/Hbe1/Prdx5 13

Table 4

KEGG Pathway Results for Purple Module.

ID Description p p adjusted q Count
rno00860 Porphyrin and chlorophyll metabolism 1,18E+08 3,17E-04 3,11E-04 10
rno05144 Malaria 2,42E-04 2,72E-02 2,67E-02 9
rno04966 Collecting duct acid secretion 3,03E-04 2,72E-02 2,67E-02 6
rno04216 Ferroptosis 5,33E-04 3,59E-02 3,52E-02 7