Co-Expression Analysis of Blood Cell Genome Expression to Preliminary Investigation of Regulatory Mechanisms in Uremia.

BACKGROUND Uremia involves a series of clinical manifestations and is a common syndrome that occurs in nearly all end-stage kidney diseases. However, the exact genetic and/or molecular mechanisms that underlie uremia remain poorly understood. MATERIAL AND METHODS In this case-control study, we analyzed whole-genome microarray of 75 uremia patients and 20 healthy controls to investigate changes in gene expression and cellular mechanisms relevant to uremia. Gene co-expression network analysis was performed to construct co-expression networks using differentially expressed genes (DEGs) in uremia. We then determined hub models of co-expressed gene networks by MCODE, and we used miRNA enrichment analysis to detect key miRNAs in each hub module. RESULTS We found nine co-expressed hub modules implicated in uremia. These modules were enriched in specific biological functions, including "proteolysis", "membrane-enclosed lumen", and "apoptosis". Finally, miRNA enrichment analysis to detect key miRNAs in each hub module found 15 miRNAs that were specifically targeted to uremia-related hub modules. Of these, miRNA-21-3p and miRNA-210-3p have been identified in other studies as being important for uremia. CONCLUSIONS In summary, our study connected biological functions, genes, and miRNAs that underpin the network modules that can be used to elucidate the molecular mechanisms involved in uremia.

Background

Chronic kidney disease is a clinical syndrome consisting of a variety of symptoms and metabolic disorders that arise from the progressive and irreversible decline of kidney function. The end stage of chronic kidney disease is uremia [1-3], which is not an independent kidney disease, but rather a series of clinical manifestations that make uremia a syndrome that is common to most end-stage kidney diseases. During end-stage uremia, nearly all vital organ systems are affected. For example, uremia patients can experience severe conditions such as heart failure, psychosis, coma, or other life-threatening complications [4,5]. Uremia can be caused by many factors such as chronic glomerulonephritis, chronic pyelonephritis, renal tuberculosis, renal arteriosclerosis, stones in the urinary tract, enlarged prostate, bladder cancer, lupus, and diabetes [5,6]. The kidneys can be damaged by various factors that severely affect kidney function, leading to defects in nitrogen metabolism and waste excretion. Excess accumulation of wastes can upset internal balances that in turn cause further damage to the kidneys and other organs. Renal failure can also be caused by extra-renal factors such as heart failure, urinary tract obstruction, and urinary reflux [7,8].

Uremia results in the bodily retention of many different compounds that can accumulate to dangerous levels [9]. Current research on uremia has mainly focused on the relationship between metabolomics and renal disease [10], and more than 90 different uremia-related compounds, some of which are toxic, have been identified [11]. However, the exact genetic and/or molecular mechanism of uremia is still not fully understood because of its complex etiology and the difficulty in obtaining tissue samples. Thus, several studies have examined transcriptome expression of blood or urine in uremia, and results from these studies identified many miRNAs and genes, as well as related biological functions and pathways, that are important for uremia [12-15].

For a more comprehensive understanding of the mechanisms that are involved in uremia, we analyzed a whole-genome microarray case-control study of 75 uremia patients and 20 healthy controls included in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) [14], and performed a gene co-expression analysis to construct networks of differentially expressed genes (DEGs) in uremia. We also determined hub co-expression models from the co-expression networks. Finally, miRNA enrichment analysis was used to detect key miRNAs in each hub model. We hope the results from this study may provide a more comprehensive understanding of the molecular mechanisms that are important in uremia.

Material and Methods

Data Preprocessing

For this study we selected the gene expression profile dataset with accession number GSE37171 [14] from the NCBI GEO. The dataset was acquired using an Affymetrix Human Genome U133 Plus 2.0 Array (platform number GPL570) that analyzed peripheral blood from 75 end-stage renal failure (uremia) patients and 20 healthy controls. The raw data were first normalized by Robust Multi-array Average (RMA) background subtraction, and quantile normalized using Affymetrix Expression Console Software; all background noises were removed. Then we used Significance Analysis of Microarrays (SAM) algorithm (|log fold change|>1; FDR <0.05) [16] to screened differentially expressed genes (DEGs) between uremia patients and control subjects.

Gene co-expression analysis and identification of co-expressed hub modules

Gene co-expression networks were performed by Pearson’s correlation with |r| >0.9 and p<0.01. The co-expressed hub module analysis was performed on the co-expression network using the Molecular Complex Detection (MCODE) plugin in Cytoscape with a degree cutoff=2, node score cutoff=0.2, k-core=2, and maximum depth=100 [17]. Significant co-expressed modules with MCODE score >4 and more than six nodes were selected.

We used the Database for Annotation, Visualization, and Integrated Discovery (DAVID) [18] to identified the biological processes and/or functions involving the hub modules. One of DAVID’s features is a functional annotation cluster tool wherein similar gene ontology (GO) categories, based on the parent/child GO term associations and the number of shared genes, can be grouped into a functional cluster. We used this feature to identify the relationship between the GO terms (biological processes in the GO database) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The threshold enrichment score for the enriched clusters was arbitrarily set as 1.

miRNA enrichment analysis

In this study the target genes of all miRNAs were predicted using five prediction algorithms: PicTar [19], DIANA-microT [20], miRanda [21], RNA22 [22], and TargetScan [23]. Genes that were predicted by at least three algorithms were considered as target genes. We collected 249 miRNAs, 5,828 genes, and 60,038 interactions of miRNA to gene.

The miRNA enrichment analysis is a useful tool to reveal modules of target gene sets with expression levels that change significantly relative to the background. These target gene sets can be used to illustrate potential microRNA regulatory mechanisms. We enriched the DEGs in each module to genes targeted by certain miRNA. The p values of certain miRNAs and their target genes was estimated by the cumulative hypergeometric algorithm, and the function is as follows:

In which M is the total number of genes tested, N is the number of genes targeted by a certain miRNA, and K represents the number of DEGs in a given module; x represents the number of genes overlap between DEGs in a given module and genes targeted by a certain miRNA. At last, the threshold of the p value was p<0.01.

Results

Identification of DEGs in uremia patients

We identified the DEGs between uremia patients and controls by the SAM algorithm (|log fold change| >1; FDR <0.05) and found 1,399 DEGs, including 311 upregulated genes and 1,088 downregulated genes.

Gene co-expression networks, co-expressed hub modules, and functional enrichment analysis

To investigate the interrelation between uremia-related DEGs, co-expression networks analysis was used to construct DEGs interaction networks in this study. As the results showed in Figure 1, we found a total of 1,046 genes (179 and 867 upregulated and downregulated genes, respectively) and 59,290 interactions between the genes (Figure 1). The interactions could be grouped into four networks that had more than four nodes, including one network wherein all nodes were downregulated genes (861 nodes, 57,604 interactions) and three networks that had nodes where all genes were upregulated (132 nodes, 1,589 interactions; 32 nodes, 83 interactions; 5 nodes and 5 interactions, respectively).

Figure 1

DEG co-expression network. Red and green nodes represent upregulated and downregulated genes, respectively.

We next performed a co-expressed hub module analysis on the co-expression networks using the Molecular Complex Detection (MCODE) plugin in Cytoscape (MCODE score >4 and nodes >6) [17]. We also examined the biological processes and/or functions involving the modules by performing functional enrichment and clustering analysis tools in DAVID [18]. This analysis identified nine modules (Figure 2, Table 1): M1 (score=88.249) containing 221 downregulated DEGs and 19,503 interactions that were significantly enriched in “proteolysis”, “membrane-enclosed lumen”, “negative regulation of gene expression”, “apoptosis”, and “RNA splicing via transesterification reactions”; M2 (score=16.976) containing 126 downregulated DEGs and 2,139 interactions that were significantly enriched in “membrane-enclosed lumen”, “chromosome organization”, “nucleotide binding”, “endosome”, and “Golgi apparatus”; M3 (score=14.944) containing 36 upregulated DEGs and 538 interactions that were significantly enriched in “membrane-enclosed lumen”;M4(score=6.944) containing 72 down DEGs and 500 interactions that were significantly enriched in “protein localization”, “post-transcriptional regulation of gene expression”, “RNA processing”, and “hormone secretion”; M5 (score=3.778) containing nine upregulated DEGs and 34 interactions that were significantly enriched in “metal ion binding”; M6 (score=3.643) containing 14 upregulated DEGs and 51 interactions that were significantly enriched in “induction of apoptosis”; M7 (score=3.415) containing 41 downregulated DEGs and 141 interactions that were significantly enriched in “mRNA processing” and “GTPase regulator activity”; M8 (score=2.267) containing 30 downregulated DEGs and 68 interactions that were significantly enriched in “protein localization”; M9 (score=2.071) containing 14 upregulated DEGs and 29 interactions that were significantly enriched in “membrane-enclosed lumen”.

Figure 2

Nine co-expressed hub modules from gene co-expression networks. Red and green nodes represent upregulated and downregulated genes, respectively.

Table 1

Functional enrichment cluster analysis for 9 co-expressed hub modules.

Module	Annotation Cluster	Gene number	Enrichment score	P value range
M1	Proteolysis	18	3.14	0.00019–0.072
	Membrane-enclosed lumen	29	1.58	0.018–0.025
	Apoptosis	9	1.35	0.042–0.046
	RNA splicing via transesterification reactions	5	1.05	0.031–0.061
M2	Membrane-enclosed lumen	37	2.06	0.0035–0.029
	Chromosome organization	10	1.73	0.0097–0.052
	Nucleotide binding	31	1.67	0.00083–0.087
	Endosome	8	1.4	0.0046–0.42
	Golgi apparatus	11	1.37	0.03–0.059
M3	Membrane-enclosed lumen	7	1.02	0.034–0.32
M4	Protein localization	10	1.8	0.0046–0.018
	Post-transcriptional regulation of gene expression	7	1.75	0.00012–0.53
	mRNA processing	7	1.54	0.008–0.12
	Hormone secretion	4	1.42	0.0038–0.38
M5	Metal ion binding	5	1.2	0.022–0.094
M6	Apoptosis	2	1.01	0.043–0.092
M7	mRNA processing	4	1.48	0.017–0.086
	GTPase regulator activity	3	1.01	0.044–0.15
M8	Protein localization	4	1.19	0.058–0.082
M9	Membrane-enclosed lumen	3	1.01	0.031–0.33

miRNA enrichment analysis of hub modules

miRNA enrichment analysis was used to predict miRNAs that could significantly target genes in each module. This analysis clearly showed that 15 miRNAs were selectively enriched in three modules (Table 2, Figure 3, p<0.01). Of these, hsa-miR-505-3p,hsa-miR-490-3p, hsa-miR-411-5p, hsa-miR-212-3p, hsa-miR-154-5p, and hsa-miR-132-3p were significantly enriched for DEGs in M1, which may indicate that these miRNAs are involved in regulating “proteolysis”, “membrane-enclosed lumen”, “negative regulation of gene expression”, “negative regulation of gene expression”, “apoptosis”, and “RNA splicing via transesterification reactions” functions by targeting the DEGs in this module. Furthermore, the miRNAs hsa-miR-875-5p, hsa-miR-590-3p, hsa-miR-448, hsa-miR-381-3p, hsa-miR-300, hsa-miR-22-3p, hsa-miR-203a, hsa-miR-154-5p, and hsa-miR-153-3p were significantly enriched in M4 DEGs, suggesting that these miRNAs may regulate “protein localization”, “posttranscriptional regulation of gene expression”, “RNA processing”, and “hormone secretion” functions. The miRNA hsa-miR-139-5p was significantly enriched in M8 DEGs, such that this miRNA may affect “protein localization” functions. Notably, hsa-miR-154-5p was significantly enriched with DEGs in both M1 and M4.

Table 2

miRNA enrichment analysis results for module DEGs.

miRNA	Module target gene	Module	P
hsa-miR-505-3p	10	M1	0.000488
hsa-miR-490-3p	5	M1	0.009921
hsa-miR-411-5p	6	M1	0.000655
hsa-miR-210-3p	13	M1	0.007337
hsa-miR-154-5p	7	M1	0.00217
hsa-miR-132-3p	15	M1	0.00154
hsa-miR-875-5p	3	M4	0.001269
hsa-miR-590-3p	13	M4	0.00063
hsa-miR-448	8	M4	0.003215
hsa-miR-381-3p	8	M4	0.001984
hsa-miR-300	12	M4	0.000338
hsa-miR-21-3p	7	M4	0.000219
hsa-miR-203a	8	M4	0.007956
hsa-miR-154-5p	3	M4	0.006956
hsa-miR-153-3p	7	M4	0.009682
hsa-miR-139-5p	3	M8	0.008584

Figure 3

Key miRNAs that significantly target the hub modules; vee node: miRNA.

Discussion

Gene co-expression network analysis is a useful tool in genomics studies, because of its extraction of genes/miRNAs with similar expression patterns from high throughput molecule expression datasets [24]. The analysis could identify the relationships between genes and facilitated the understanding of the regulatory mechanism of genes [25]. In the present study, we used the DEGs to construct the gene co-expression networks, which not only could identify the gene interactions/networks related to uremia, but also could reduce the false positive rate of the DEGs [26]. Meanwhile, the hub module analysis of differentially expressed genes identified clusters of genes with similar expression profiles that are suggestive of co-regulation. Exploration of these co-regulated gene modules could help find underlying regulation molecules that are pivotal for gene expression in uremia [27]. To explore the molecular mechanisms of DEGs in uremia, miRNA enrichment analysis was used to predict pivotal miRNAs that could significantly target genes in each functional module. Results from these bioinformatics studies might inform future studies of the pathways and mechanisms that could be targeted for uremia treatments.

Blood is an important material for investigation of chronic kidney disease (CKD). Recently, some reports have explored the molecule alterations in the blood of CKD patients [28-31]. Scherer et al. [14] found that over one third of human genes were differential expressed in blood of uremia, and involved in various biological functions and pathways, including many immune and biological mechanisms. Meanwhile, Granata et al. [13] screened genome expression in peripheral blood mononuclear cells of CKD patients, and found some DEGs involved in mitochondrial respiratory system. Moreover, Neal et al. [15] found the mass of miRNAs were reduced in the blood of CKD patients. All these reports indicated a profound magnitude of the changes in CKD, but not enough to define further critical molecular mechanisms in CKD. However, the bioinformatics tools used in our study might provide us with insight into critical biological alterations.

In this study, co-expression analysis and co-expressed hub module analysis found nine hub modules that may play important roles in uremia. For example, module M1 was significantly enriched in “proteolysis”, “membrane-enclosed lumen”, “negative regulation of gene expression”, “apoptosis”, and “RNA splicing via transesterification reactions”. These biological functions have been demonstrated to have significant roles in uremia [32-35]. In accord with the clinical characteristics of uremia, protein fragments in plasma and urine are a hallmark of many renal diseases, and can be an important component of plasma and urine for diagnostic and research purposes. Meanwhile, uremia is typically accompanied by apoptotic erythrocyte death that may lead to anemia. Module M4 was significantly enriched in “protein localization”, “posttranscriptional regulation of gene expression”, “RNA processing”, and “hormone secretion” functions, suggesting that in uremia homeostasis may be severely affected and accompanied by significant disturbances in RNA and protein metabolism [36,37]. Notably, the “membrane-enclosed lumen” function was significantly enriched in modules M1, M2, M3, and M9. Indeed, uremia is typically accompanied by changes in cellular homeostasis that involve metabolic disturbances in membrane-enclosed lumens such as the nucleus and endoplasmic reticulum.

microRNAs are small non-coding, single-stranded RNAs that regulate messenger RNAs at the post-transcriptional level. The role of miRNAs in chronic kidney disease (CKD) has recently gained increased attention [38,39]. In this study, we found 15 key miRNAs that were significantly enriched in three co-expressed modules. Of these miRNAs, miRNA-21-3p and miRNA-210-3p were shown to be significantly related to uremia [12,15]. Moreover, hsa-miR-154-5p was significantly enriched in both M1 and M4 DEGs. However, the other miRNAs identified in this study await further characterization and investigation for their role in uremia-related mechanisms.

Conclusions

Our study provides a comprehensive perspective of gene expression in uremia using a combination of gene co-expression network analyses and DEGs analyses to identify nine hub co-expressed modules that can be assigned to corresponding biological functions. The 15 key miRNAs identified in the miRNA enrichment analysis likely also play important roles in uremia. Further study of the networks and miRNAs identified in this study may allow a better connection of biological functions, genes, and miRNAs that underpin the network modules, and in turn could be used to elucidate the molecular mechanisms involved in uremia. In conclusion, our study provided an orientation for further work, and the bioinformatics analysis pipeline of this study provides a model for others. However, this work is still preliminary and needs more confirmatory research in future.