Quality control and identification of DEGs
Fig. 1 shows the overall flowchart of this present research. The processed data boxplot displayed that the medians of each sample were almost on the same line (Fig. 2A). The PCA and UMAP plots suggested a good distinction between DN and control samples (Fig. 2B, C). After downloading the tubulointerstitial data from GSE30529, we recognized 463 DEGs that consist of 340 up-regulated DEGs and 123 down-regulated DEGs. Fig. 2D and Fig. 2E showed the volcano map and heatmap of the DEGs.
Analysis of the functional characteristics
To evaluate GO analysis of all genes, we performed GSEA analysis. We found that the most enriched gene sets were associated with extracellular matrix structural constituent, complement activation, and positive regulation of T cell proliferation in DN samples (Fig. 3). We subsequently performed the GO and KEGG pathway analyses of DEGs. Compared to the controls, DN patients are significantly enriched in the extracellular matrix and inflammatory response, including collagen-containing extracellular matrix, extracellular structure organization, and regulation of immune effector process (Fig. 4A). Based on the KEGG pathways, the top 10 enriched pathways involved Phagosome, Staphylococcus aureus infection, Pertussis, Leishmaniasis, Tuberculosis, Complement and coagulation cascades, Cell adhesion molecules, Viral myocarditis, Systemic lupus erythematosus, and Toxoplasmosis (Fig. 4B).
PPI network construction and module analysis
The complex PPI network was constructed with 277 nodes and 803 edges in Fig. 5A. Through plugin MCODE of Cytoscape, we recognized 3 significant modules with an MCODE score > 5. Module one had the highest score with 11 (11 nodes and 55 edges), module two had the second-highest score with 5.625 (17 nodes and 45 edges), followed by module three had the third-highest score with 5.60 (6 nodes and 14 edges) (Fig. 5B-E).
Identification of hub genes and prediction of target miRNAs
Four algorithms were used to identify the hub genes via the Venn diagram. As a result, 3 intersected genes were obtained (Fig. 6A). Table 1 lists the hub genes together with their full name and descriptions. We then explored 24 target miRNAs of 3 hub genes based on the regulatory relationship between mRNAs and miRNAs. There are 25 mRNA-miRNA pairs in the co-expressed mRNA-miRNA network (Fig. 6B).
Validation and ROC curve of hub genes
To assess the reliability of these hub gene expressions, we used validation datasets to explore the mRNA levels of the obtained hub genes. The results displayed that 3 hub genes showed up-regulated mRNA levels in DN (Fig. 7A). To further evaluate the sensitivity and specificity of the hub genes, we analyzed the ROC curve with the validation datasets. Among the results, CSF2RB had the highest areas under the curve (AUC) score ( AUC value: 0.897), CD53 had the second-highest AUC score ( AUC value: 0.846), followed by LAPTM5 had the third-highest AUC score ( AUC value: 0.843) (Fig. 7B).
Clinical analysis
Based on the Nephroseq v5 platform, we investigated the possible impact of the hub genes in tubulointerstitial samples in DN. Correlation analyses showed that the mRNA levels of CSF2RB, CD53, and LAPTM5 in human tubulointerstitial samples were negatively correlated with GFR (Fig. 8A-D), indicating that these hub genes may accelerate the tubulointerstitial lesions of DN. Then, the mRNA levels of CD53 and LAPTM5 in human tubulointerstitial samples had positive correlations with proteinuria (Fig. 8E, F), further suggesting that CD53 and LAPTM5 may be involved in the progression of tubulointerstitial injuries of DN. Additionally, we found that, in contrast with the subnephrotic proteinuria group, the mRNA expression level of LAPTM5 in renal tubulointerstitial samples was higher in the nephrotic proteinuria group (Fig. 8G). Furthermore, the mRNA levels of CSF2RB, CD53, and LAPTM5 were conversely related to their age (Fig. 8H, J).
Construction of ceRNA networks
It is well known that miRNAs can inhibit mRNA translation or lead to mRNA degradation by binding to target mRNAs, thus achieving the regulation of gene expression. However, the ceRNA hypothesis interpreted that IncRNAs, circRNAs, and pseudogenes that act as competing endogenous RNAs or natural microRNA sponges can increase gene expression by binding microRNA response elements (MREs) [21]. Based on the Starbase platform, we predicted corresponding IncRNAs and circRNAs of target miRNAs. Finally, we selected the IncRNAs, circRNAs that appeared in most of the forecasted results of target miRNAs to be our final IncRNAs, circRNAs. Of the predicted results in the Starbase tool, we regarded the circRNAs that qualified with the highest score or /and the most samples in the circBase database to be the eventual target circRNAs since one transcript corresponds to multiple shear sites of circRNA.
Just as the forecasted results of circRNAs and IncRNAs (Fig. 9A, B, C), we acquired 22 circRNAs and 1 IncRNA of 3 miRNAs of CD53, which were paired into 66 circRNA-miRNA interactions, 3 IncRNA-miRNA interactions, and 3 miRNA-mRNA interactions; 9 circRNAs and 2 IncRNAs of 18 miRNAs of CSF2RB, which were paired into 135 circRNA-miRNA interactions, 18 IncRNA-miRNA interactions, and 29 miRNA-mRNA interactions; and 5 circRNAs and 1 IncRNA of 4 miRNAs of LAPTM5, which were paired into 20 circRNA-miRNA interactions, 4 IncRNA-miRNA interactions, and 4 miRNA-mRNA interactions. Given the ceRNA hypothesis, we subsequently put two down-regulated miRNAs and two up-regulated IncRNAs into further analyses. According to these data, we can infer that the regulatory pathway of NEAT1/XIST-hsa-miR-155-5p/hsa-miR-486-5p-CSF2RB may play an essential role in DN (Fig. 9D).