Development and assessment of diabetic nephropathy prediction model using hub genes identified by weighted correlation network analysis

GEO datasets

We downloaded datasets GSE142153 and GSE26168 (both expression profiling by array) from the GEO database (https://www.ncbi.nlm.nih.gov/gds). The GSE142153 based on the GPL6480 platform includes 3 condition experiments: Healthy control (n=10), DN (n=23), and end-stage renal disease (n=7). The tissue type was from peripheral blood. The GSE26168 based on the GPL6883 platform includes healthy controls (n=8), impaired fasting glucose group (n=7), and T2DM group (n=9). We used the DN group of GSE142153 and the T2DM group of GSE26168 for further analyses.

Patients

This is a retrospective cohort study design. The patients included here had T2DM in the hospital. We collected the baseline data from November 2013 to June 2014, and acquired the follow-data from September 2017 to June 2018. We identified the T2DM according to one or more of the following criteria: HbA1c ≥ 6.5%, fasting blood glucose ≥ 7.0 mmol/L, oral glucose tolerance test ≥ 11.1 mmol/L, and random blood glucose ≥ 11.1 mmol/L [9]. The criteria for inclusion were: age under 75 years old; no severe disease including malignant tumor, uremia dialysis, chronic inflammatory disease, autoimmune diseases, and cerebra-cardiovascular disease. Patients are conscious and can communicate with the investigator. Patients with >30 mg/g (albumin /creatinine ratio) at baseline, or incomplete data or information were excluded. This study was approved by the Ethics Committee of China-Japan Friendship Hospital, and all subjects provided written informed consent.

Bioinformatics analysis

With R packages “limma” and “sva”, we first merged the DN group of GSE142153 and the T2DM group of GSE26168. We used the “removeBatchEffect” function in the R package limma. Differential expression analyses were performed between the two groups. The differentially expressed genes (DEGs) were identified using the following criteria: |log2 fold change|>2 and adjusted P <0.05.

WGCNA is a system biology method used to describe gene association patterns between different samples [10, 11]. It can be used to identify highly covarying gene sets and to find out candidate biomarker genes or therapeutic targets based on the interconnectedness of gene sets and the association between gene sets and phenotypes. The WGCNA consists of gene co-expression network construction, and identification of modules, related modules to external information, relationships between modules, and key drives in interesting modules. We used the following criteria to identify the related modules: the module size of 10-20000, the cut height equal to 0.25 and verbose equal to 5. Highly related genes were identified with thresholds greater than 0.1 in the topological overlap matrix (TOM).

We made an overlap between DEGs and related module genes identified by WGCNA. The intersection genes were used for further analyses. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genome (KEGG) pathway analyses were achieved using R packages “clusterProfiler”, “enrichplot” and “ggplot2”, and the results were presented in bar plots [12, 13]. We further established protein-protein interaction (PPI) using String (https://www.string-db.org/). Any gene without nodes with other genes was excluded. The PPI was inputted into Cytoscape [14]. Using the application of cytoHubba, we identified the top 10 hub genes using five algorithms: maximal clique centrality (MCC), maximum neighborhood (MN), maximum neighborhood component (MNC), Degree, and edge percolated component (EPC). Intersections among the five algorithms were made, which were considered as the key genes.

Clinical data collection

The following information was collected for each patient: age, gender, smoking status, drinking status, body mass index (BMI, weight/height², kg/m²), systolic blood pressure, and diastolic blood pressure. The venous blood was extracted before receiving any treatment. The fasting blood glucose (FBG), HbA1c, total cholesterol (TC), triglyceride (TG), high-density lipoprotein (HDL), low-density lipoprotein (LDL), blood urea nitrogen (BUN), and serum creatinine (SCr) were measured using the Jaffe rate-blanked compensated creatinine assay. Urine acid (UA) was detected using an automated biochemical analyzer. The normal and abnormal standards were from the Guidelines for Prevention and Treatment of Dyslipidemia in Chinese Adults [15]. Patients with urinary albumin creatinine ratio >30 mg/g were diagnosed with DN [16]. According the World Health Organization, people who have smoked continuously or cumulatively for 6 months or more were defined as smokers. Men and women who drank more than 60 g and 40 g of pure alcohol per day respectively were considered as drinking. BMI was calculated as BMI = weight (kg)/height (m)². Hypertension was defined as systolic blood pressure ≥ 140 mmHg and/or diastolic blood pressure ≥ 90 mmHg.

For evaluating the expression levels of seven hub genes, we extracted the total RNA from the peripheral blood of each T2DM patient using a modified RiboPure™-Blood kit (Ambion, Austin, USA) according to the manufacturer’s protocol. The cDNA was synthesized using a reverse transcription kit (Guangzhou Ribobio, Co., Ltd., China). Quantitative real-time polymerase chain reaction (qRT-PCR) was carried out using a Cobas Z480PCR system. The relative expressions of the hub genes were calculated using the 2^-ΔΔCT method. The forward and reverse primer sequences were as follows: TTCTCAGAGGCTGACATGCGCT and CTCGTCCAGAAGGATGTTGGCT (ADRBK1); CATCGAGCAGAGTGTCTACAGC and TGTCGTCGCTATCCAGG TCATG (LMX1A); CACCACGCTCTCCAATGCCTTT and CATCACCAGGATGG ACACAAGC (HTR1B); GACAACCACTGTGAACATCACGC and TGACACTTC CAATGGCTGTGCC (CDH10); TTACTGGTCCGACTGGTACGAC and CAAAGA ACTGCTGAGGCTTGGG (TAC1); CTCCTCTTTGTCATCACGCTTCC and GGA TGAGGACACTGCTGTAGAG (NPY); CTGGAACGTGACCAACGCCATC and TCATTCTCCTCGTACAGGCACG (SLC32A1).

Statistical analysis

All continuous variables were transformed into categorical variables, and expressed using counts and percents like other original categorical variables. We first screened relevant genes from the seven hub genes using the LASSO method, and then calculated the risk score for each patient as follows: risk score =gene₁ expression*β₁+…+gene_n expression*β_n. Then the risk score was inputted as a continuous variable into logistic regression. We established three models for the prediction of DN: model 1 with only clinical parameters, model 2 with only identified gene expression, and model 3 with combined clinical parameters and gene expression. Multivariate logistic regression was used to validate if the risk score was an independent predictor of DN among T2DM patients or not. The receiver’s operating characteristic (ROC) curves were used to assess the discrimination ability of the models. Nomograms were adopted to quantify the predictive ability of the models. Calibration curves were plotted to evaluate the calibration of DN risk, which was accompanied by a Hosmer-Lemeshow test. Decision curve analysis was conducted to evaluate the clinical practicability of the models based on the net benefit according to the different threshold probabilities among T2DM patients.

Data availability

The expression profiling of genes from diabetics and diabetics nephropathy can be available from the GEO (GSE142153 and GSE26168). The other data can be obtained from the corresponding author (Bo Zhang).

Differential expression analysis of genes

Using 23 DN samples and 9 diabetic samples, we performed differential expression analysis to obtain the DEGs that can be associated with DN. First, we merged the expression data of the DN group of GSE142153 and the T2DM group of GSE26168 to exclude the bath effects. The boxplot showed the whole gene expression levels of the merged dataset were approximately the same, which indicated no batch effect (Figure 1A). Using the criteria of |log2-fold changes (FC)|>1.0 and adjusted P <0.05, we obtained 110 DEGs, including 9 upregulated DEGs and 101 downregulated DEGs between the two groups (Figure 1B). The top 60 DEGs were presented in the hierarchical clustering heatmap (Figure 1C), including the significantly upregulated or downregulated DEGs.

Weighted gene co-expression network analysis (WGCNA)

Using WGCNA that describes the expression correlation among multiple genes, we identified the genes of interest. Two analytical methods were applied to identify the target genes. First, we performed the cluster analysis to assess sample similarity. The outlier’s detection identified three samples that were removed (GSM4221586, GSM532842, GSM4221594, Figure 2A). To construct a scale-free network, we performed WGCNA, and the scale-free topological fitting index was 0.88 and the mean connectivity was 90 by setting the soft threshold power at 10 (Figure 2B). Furthermore, we carried out hierarchical clustering to identify different kinds of gene modules using weighted gene co-expression correlation, which can be found through the branches of clustering tree and different colors. We identified 18 modules in size from 20 to 10000 in the network, and merged cut highs of 0.25 (Figure 2C). The 18 modules can be separated into three clusters according to the correlations among the modules (Figure 2D). Besides, the weighted coexpression correlations of all genes were presented in a heatmap plot (Figure 2E). The module MEmidnightblue showed the highest correlations, the correlation coefficients between the control group and the DN group were 0.99. Finally, we obtained 136 highly-related genes in the TOM matrix with a threshold greater than 0.1. The results of the two analyses can be overlapped, and a list of 110 target genes were identified, which may be involved in DN (Figure 3A).

Function enrichment analysis of targeted genes

The occurrence of DN may involve complex pathogenesis and biological processes, and understanding the potential function of target genes will help with DN treatment. The GO enrichment analysis indicates these target genes are mainly related to membrane adhesion, alpha-amino acid biosynthesis, metabolism, and binding, terminus, inhibitory synapse, clathrinid-sculpted vesicle, kinase activity, hormone binding, receptor activity, and transporter activity (Figure 3B). KEGG analysis indicates synaptic vesicle cycle, cysteine and methionine metabolism, N-Glycan biosynthesis, osteoclast differentiation, and cAMP signaling pathway can be involved in the process of DN in diabetic patients (Figure 3C, 3D). We also built the PPI based on these target genes (Figure 4A), and identified the top 10 hub genes: LMX1A, TYRP1, NANOG, ADRBK1, HTR1B, TAC1, NPY, SLC32A1, GAD1, and CDH10 (Figure 4B). Finally, we identified the hub genes using five algorithms (Figure 4C). We overlapped the results of the five algorithms and obtained the final hub genes that were used for model building: ADRBK1, LMX1A, HTR1B, CDH10, TAC1, NPY, and SLC32A1.

Development of DN prediction model based on target genes

Based on 7 target genes, we performed LASSO regression and identified 6 genes (ADRBK1, LMX1A, HTR1B, TAC1, NPY, SLC32A1) included in the prediction model (Figure 4D). We calculated the risk score as follows: risk score = 1.466 * ADRBK1 - 0.627 * LMX1A + 0.199 * HTR1B + 0.130 * TAC1 + 0.477 * NPY + 1.005 * SLC32A1.

The general characteristics between the DN group and the non-DN group were compared. The risk scores were higher in the DN group than in the non-DN group (10.54 vs. 9.03, P<0.001; Table 1). The DN group had higher proportions of males (65.4% vs. 55.9%, P=0.042) and smokers (42.6% vs. 23.1%, P<0.001). The SBP (P<0.001) and FBG (P=0.016) levels were higher in the DN group than in the non-DN group. The TC, TG and BUN levels were also significantly elevated in the DN group (P<0.001). No significant differences were observed in age >60, drinking status, DBP>=90 rate, and HbA1c, HDL, LDL, SCR, or UC levels (P>0.05). The univariate logistic regression indicated the risk score was positively associated with DN (OR (95%CI): 2.72 (2.26-3.26), P<0.001). Besides, male (1.49 (1.03-2.17), P=0.034), FBG>6.9 (1.62 (1.11-2.38), P=0.013), overweight or obesity (1.84 (1.21-2.78), P=0.004), DBP>90 (1.55 (1.02-2.37), P=0.042), and smoking (2.47 (1.69-3.60)), SBP>140 (2.23 (1.55-3.22)), abnormal TC (2.68 (1.84-3.91)) and TG (2.16 (1.48-3.14)), and elevated BUN (2.45 (1.60-3.77)) (all P<0.001) were also risk factors of DN in diabetic patients (Table 2). Using these significant variables in the univariate logistic regression, we performed multiple logistic regression. Our results indicated that risk score was an independent risk factor of DN for diabetic patients (OR:2.76, 9%CI: 2.27-3.41, P<0.001). Besides, smoking, overweight, SBP>=140, abnormal TC, TG and BUN were also risk factors of DN (Figure 5).

Assessment of clinical application of prediction model

To validate the model, we built the following three models: model 1: risk score, mode 2: clinical parameters; model 3: risk score and clinical parameters. For models 1 and 2, the AUCs were 0.763 (95%CI: 0.721-0.805, Figure 6A) and 0.813 (95%CI: 0.775-0.851, Figure 6B) respectively. Model 3 showed the highest predicting ability (AUC:0.870, 9%CI: 0.837-0.902, Figure 6C). Based on model 3, we built a nomogram plot to evaluate the risk probability of DN for individuals (Figure 6D). We only included the significant variables in the univariate logistic regression into the nomogram. Specifically, we can get a point for each variable. We summarized the points of all variables in the nomogram. The summarized points can be projected onto the risk of nonadherence. Then we can know the risk probability of individuals. These descriptions were added in the manuscript.

Furthermore, we evaluated the agreement between actual diagnosed probability and nomogram-predicted probability. Results showed relatively high agreement (Figure 7A). The calculated C-index was 0.870, suggesting that the model has high prediction accuracy. Moreover, decision curve analysis for DN risk indicated that T2DM patients will benefit from this DN risk prediction model when the threshold probability of a patient and clinician was 25%, meaning the net benefit is relatively high based on the DN risk nomogram (Figure 7B). To further assess the model, we randomly extracted half of the whole study sample to validate the prediction model. The internal validation showed that the AUC was 0.821 (95%CI: 0.756-0.886, Figure 7C), and the calibration curve of this sample also indicated high performance (Figure 7D). These results suggest that the DN prediction model can help clinicians make decision for T2DM patients.