Construction of a drug, miRNA, and transcription factor regulatory network in membranous nephropathy based on oxidative stress-related genes

Abstract

In this study, oxidative stress-related genes (OSRGs) were used to discover membranous nephropathy (MN) genes, evaluate their biological activities and molecular pathways, and develop regulatory networks for drugs, miRNAs and transcription factors. All 862 differentially expressed genes (DEGs) in the GSE108113 dataset were evaluated using the gene expression matrix. Cross-testing of 437 OSRg in the gene cards revealed 28 DE-OSRG genes. We performed the KEGG and gene ontology enrichment studies. Response to reactive oxygen species, chemical stress, oxidative stress, bound globin-hemoglobin complex, endocytic vesicle lumen, blood microparticles, conjugated globin binding, peroxidase activity and heme binding are the major GO-enriched active DE-OSRGs. As in the KEGG pathway study. After excluding the disconnected genes, 19 DE-OSRGs were introduced into the PPI network, and the top five CDE-OSRGs were found: G6PD, HBB, TP 53, KEAP 1, and PXDN. The ROC curve study shows the CDE-OS.

KEYWORDS:

• Membranous nephropathy
• oxidative stress
• regulatory network
• core gene
• miRNA
• transcription factor

1. Introduction

Adult nephrotic syndrome is mostly caused by Membranous Nephropathy (MN). Membranous nephropathy is a disease in which immune complexes deposit on the basement membrane of glomerular capillaries, leading to basement membrane lesions. It is characterized by nephrotic syndrome, including proteinuria, hypoalbuminemia, hyperlipidemia, and edema [1, 2]. Idiopathic membranous nephropathy accounts for 75%∼80% of the total incidence rate of membranous nephropathy. Its specific cause has not been determined yet, and its pathogenesis is complex. It is generally believed that the deposition of immune complexes causes damage to the basement membrane and podocytes, resulting in proteinuria [3]. At present, the clinical treatment mainly includes symptomatic support treatment, including rest, sodium restriction, reasonable protein and calorie intake, diuretic detumescence, ACEI/ARB urine protein lowering, statins for lipid regulation, preventive anticoagulation and other measures. Corticosteroids and cytotoxic drugs can also be used for immunosuppressive therapy [4].

Oxidative stress (OS) is a condition in which the body’s oxidative and antioxidant actions are out of balance, leading to oxidation and inflammatory neutrophil infiltration, elevated protease release, and a high concentration of oxidative intermediates [5]. OS is an important mechanism in immune response, mainly formed by intracellular oxidation processes, involving complex gene expression and signal transduction networks. OS can be achieved through the regulation of oxidative stress-related genes (OSRGs) (mainly composed of nuclear factors, proteases, antioxidant factors and detoxification enzymes) [6]. OS works both ways. It is a harmful byproduct of free radicals in the body that is thought to play a major role in ageing and disease. Secondly, it can preserve the body’s regular biological processes and serve as a crucial benchmark for the diagnosis, treatment, and prevention of a wide range of illnesses [7, 8].

Targeted drugs refer to drugs or their formulations that have been endowed with targeting ability. Its purpose is to enable drugs or their carriers to target specific lesion sites and accumulate or release effective ingredients at the target site. Targeted formulations can form relatively high concentrations of drugs in the target area, thereby improving drug efficacy while suppressing toxic side effects and reducing damage to normal tissues and cells. In the preparation of drugs, the physiological structural characteristics of specific tissues and organs can be utilized to enable natural distribution differences of drugs in the body, thereby achieving targeted effects. There are also various methods, including coupling probe molecules that can specifically bind to target molecules such as antibodies, peptides, sugar chains, nucleic acid aptamers, etc., to the surface of drugs or their carriers through chemical or physical methods to achieve targeted effects.

At present, there are few literature reports on the molecular mechanism of MN based on OSRGs. So, in this study, the Gene Expression Omnibus (GEO) differentially expressed oxidative stress-related genes (DE-OSRGs) were used among MN individuals, and the general population was found using the GeneCards dataset. The biological functions and molecular pathways associated with these genes were then investigated, furthermore by building a protein–protein interaction (PPI) network, core DE-OSRGs (CDE-OSRGs) were further discovered. For these genes, correlation and receiver operator characteristic (ROC) curve analyses were carried out. Finally, the outcome explores the interaction between CDE-OSRGs and targeted drugs and predicts CDE-OSRGs binding miRNAs and transcription factors (TFs).

2. Materials and methods

2.1. Acquisition and organization of sample data

The clinical data associated with the dataset GSE108113 was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) [9]. Following this, probe IDs in the original matrix were converted into gene symbols, gene labeling was performed using the platform annotation file, and expression values were log2 transformed [10]. As a result, the gene expression matrix of MN, consisting of 87 MN samples and 11 healthy control (NC) samples, was obtained.

2.2. Differential expression gene analysis in MN

The R packages limma, pheatmap, and ggplot2 were used to conduct a statistical test on the gene expression matrix in order to determine the differentially expressed genes (DEGs) among the MN and NC categories. As a consequence, the DEGs and their expression values were extracted and visualized via a heatmap and volcano map [11, 12].

2.3. Determining the genes linked to oxidative stress that are differently expressed in Minnesota

Getting a list of genes associated with oxidative stress (OSRGs) from the GeneCards database (https://www.genecards.org/) [13] was our initial step. Next, intersections between OSRGs and DEGs were discovered using the R programme VennDiagram. This allowed us to identify the genes that are connected to differentially expressed oxidative stress (DE-OSRGs) [14]. As a final step, these DE-OSRGs were visualized, and their location was determined by their distribution in human chromosomes [15].

2.4. GO and KEGG enrichment analyses of DE-OSRGs

The names of DE-OSRGs were first translated into gene IDs compatible with the R language [16] using the R package org.Hs.eg.db. Then, to clarify the biological significance of the DE-OSRGs, Gene Ontology (GO) enrichment analysis was performed on the DE-OSRGs, which include biological processes (BP), molecular functions (MF), as well as cellular components (CC), using the R packages clusterProfiler, enrichplot, ggplot2, as well as GOplot. Furthermore, Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analysis was carried out on the DE-OSRGs to pinpoint important pathways that are intimately linked to the control of the growth and advancement of MN [17, 18].

2.5. Protein–protein interaction (PPI) network analysis

A comprehensive examination of the interaction between regulatory proteins associated with DE-OSRGs in the MN biological system was conducted via the STRING database (http://www.string-db.org/) [19]. Thus, the cytoHubba plugin in Cytoscape software was used to identify core genes inside the protein–protein interaction (PPI) network. The Degree approach was used to designate the top five hub genes as core DE-OSRGs (CDE-OSRGs). Finally, the PPI network and CDE-OSRGs subnetwork were visualized [20].

2.6. Correlation analysis and receiver operator characteristic (ROC) curve analysis of CDE-OSRGs

Boxplots 21 and 22 were used to display the relative expression levels of CDE-OSRGs in various samples using the R packages limma and ggpubr. These levels were retrieved using the gene expression matrix and compared to the differential expression of these genes within the MN and NC categories [21, 22]. To verify the diagnostic accuracy of CDE-OSRGs, an ROC curve was constructed by referencing R packages glmnet and pROC, using the quantitative accuracy metric of the area under the ROC curve (AUC) [23]. Additionally, the correlation among CDE-OSRGs was assessed using the R packages ggplot2, ggpubr, ggExtra, corrplot, and circlize (correlation coefficient > 0.3 and P value < 0.05 were deemed statistically significant). Additionally, scatter plots and a correlation heat map were created for visualization [24].

2.7. Examination of the relationship among target medications and CDE-OSRGs

The Drug-Gene Interaction database (DGIdb) (https://dgidb.genome.wustl.edu/), which serves as a comprehensive collection of Drug-Gene Interactions, provides genetic data related to both established and potential drugs [25, 26]. In this study, the search field was initially populated with CDE-OSRGs. Subsequently, filtering conditions were applied to specify the drug type as Approved, Antineoplastic and Immunotherapies. The default selections for database source, gene category, and interaction type were maintained. Following this, the resulting gene-drug interaction file was downloaded. Finally, the software Cytoscape was utilized to visualize the interaction between the CDE-OSRGs and the target drugs.

2.8. Prediction of miRNAs and transcription factors (TFs) binding to CDE-OSRGs

An extensive atlas of miRNA (or TF) – target gene interactions is available at https://www.mirnet.ca/, the home of the miRNet database [27]. 14 distinct miRNA databases are integrated (TarBase, miRTarBase, miRecords, miRanda, miR2Disease, HMDD, PhenomiR, SM2miR, PharmacomiR, EpimiR, starBase, TransmiR, ADmiRE, and TAM 2.0) with regard to functional annotation information [28]. First, a list of CDE-OSRGs was added to the search box of the database, with “Organism” selected as “Homo sapiens”, “ID type” selected as “Official gene symbol”, “Tissue” selected as “Kidney”, “Targeted by” selected as “miRNA” and “TF” then the information was submitted and the table of the interaction between miRNAs (TFs) and CDE-OSRGs was downloaded. Finally, Cytoscape software was used to visualize the interactions and regulatory relationships between them.

3. Results

3.1. Identification of DE-OSRGs in MN

The initial step involved compiling dataset GSE108113, which consisted of 87 MN and 11 NC samples, to create an expression matrix encompassing 24,889 genes (refer to Supplementary file Gene expression matrix.xls). Subsequently, a differential gene expression analysis was conducted to identify a total of 862 DEGs that distinguished the MN group from the NC group. This group included 496 genes with low expression and 366 genes with high expression (Figure 1(a)). The top 50 genes with the most notable variations in expression were shown visually in the heat map (Figure 1(b)). Furthermore, 437 OSRGs were collected from the GeneCards database and intersected with these DEGs and OSRGs to produce 28 DE-OSRGs (Figure 2(a)). Figure 2(b) and (c) show the location distribution of these DE-OSRGs in 23 pairs of human chromosomes.

Figure 1. The volcano plot (a) and heatmap (b) of DEGs expressed between the MN and NC groups. Red dots or squares indicate highly expressed genes; Green dots or blue squares indicate lowly expressed genes.

Figure 2. (a) The Venn diagram of overlapped genes (DE-OSRGs). (b, c) The location distribution of these DE-OSRGs in 23 pairs of human chromosomes.

3.2. Function and pathway enrichment of DE-OSRGs in MN

GO and KEGG enrichment analyses were performed to investigate the possible function of the 28 DE-OSRGs. According to the results of the GO enrichment study, the primary roles of BP are in the metabolic processes of reactive oxygen species, chemical stress and oxidative stress response control. Functions of CC were related to haptoglobin−hemoglobin complex, endocytic vesicle lumen, blood microparticle, ribosomal subunit, and heterochromatin. MF functions included haptoglobin binding, peroxidase activity, heme binding, tetrapyrrole binding, nuclear glucocorticoid receptor binding, and molecular carrier activity (Figure 3(a)). In the KEGG pathway analysis, DE-OSRGs mainly affected African trypanosomiasis, Malaria, Proteoglycans in cancer, Basal cell carcinoma, Hepatocellular carcinoma, Bladder cancer, Fluid shear stress, Thyroid cancer, and atherosclerosis, Gastric cancer and Breast cancer, etc. (Figure 3(b)).

Figure 3. Bar graphs displaying the GO results (a) and KEGG (c) analyses of enrichment.

Gene counts that are enriched in various GO items or KEGG pathways are plotted on the x-axis, and descriptions of these GO items or KEGG pathways are plotted on the y-axis. The P-value is represented by the colour of the bar graph, while the gene count is shown by its length.

3.3. PPI network of DE-OSRGs in MN

A total of 19 DE-OSRGs (the disconnected genes were removed) were introduced into the PPI network, followed by an analysis to assess the connectivity of each DE-OSRG within the interaction network. This analysis led to the identification of the top five DE-OSRGs with the highest connectivity, referred to as CDE-OSRGs. The visualization of both DE-OSRGs and CDE-OSRGs within the interaction network was conducted using Cytoscape software, as depicted in Figure 4(a) and (b). Figure 4(b) specifically displayed the five CDE-OSRGs, namely G6PD, HBB, TP53, KEAP1, and PXDN.

Figure 4. The PPI networks consisted of (a) 19 DE-OSRGs and (b)5 CDE-OSRGs. Each gene is represented by a node, and the connections between nodes represent the interactions between genes. 5de osrg is distinguished by nodes in red or orange, and the intensity of the colour reflects their degree of connectivity.

3.4. Verification and correlation of CDE-OSRGs in MN

The differential expression of CDE-OSRG was visualized in the MN and NC groups using boxplots (Figure 5(a)). These five CDE-OSRGs (G6PD, HBB, TP53, KEAP1, and PXDN) showed high expression in the MN group. As Figure 5(b) illustrates, the areas under the Curve standards of G6PD, HBB, TP53, KEAP1, and PXDN were 0.804, 0.888, 0.968, 0.936, and 0.821, respectively, indicating that CDE-OSRGs have high diagnostic accuracy in distinguishing MN patients. In addition, Figure 5(c) shows the expression correlation between CDE-OSRGs and each other. Amongst these, the PXDN expression level had a negative correlation with KEAP1 and a positive correlation with HBB. There was a positive correlation and a negative correlation between the expression level of TP53, PXDN and KEAP1.

Figure 5. Verification and correlation of CDE-OSRGs in MN. (a) Box plots of CDE-OSRGs (G6PD, HBB, TP53, KEAP1, and PXDN) expression values between MN and NC groups. (b) ROC curves of the five CDE-OSRGs. AUC, area under curve. (c) Heat map and scatter plots of the expression correlation between CDE-OSRGs.

3.5. Regulatory network of CDE-OSRGs and targeted drugs

An extensive table comprising CDE-OSRGs and their respective targeted drug interactions was obtained, as outlined in the supplementary file (Drug-gene interactions.xls), utilizing the DGIdb database. A regulatory network that visually depicted the associations between CDE-OSRGs and targeted drugs was constructed through the use of Cytoscape. As illustrated in Figure 6, the CDE-OSRGs were visually represented by red square nodes, while the molecular drugs targeted for regulation were depicted as blue diamond nodes. The lines connecting these nodes symbolized the interactions between them. This drug regulatory network encompassed a total of four CDE-OSRGs, 61 molecular drugs, and 62 pairs of interacting relationships.

Figure 6. Regulatory network between CDE-OSRGs and targeted drugs. CDE-OSRGs are represented by red nodes, targeted drugs are represented by blue nodes, and the lines between them represent interaction pairs.

3.6. Regulatory network of CDE-OSRGs and miRNAs (TF)

The tables of the interactions between miRNAs (TFs) and CDE-OSRGs, as outlined in the supplementary file (Gene-miRNA interactions.xls and Gene-TF interactions.xls), were obtained by utilizing the miRNet database. The two regulatory networks were constructed through the utilization of Cytoscape, which visually depicted the associations between miRNAs (TFs) and CDE-OSRGs. As illustrated in Figures 7 and 8, the CDE-OSRGs were visually represented by red square nodes, while the miRNAs (TFs) that bind to the CDE-OSRGs were depicted as green oval nodes (yellow diamond nodes). The lines connecting these nodes symbolized the interactions between them. The Gene-miRNA regulatory network encompassed a total of five CDE-OSRGs, 154 miRNAs, and 220 pairs of interacting relationships and the Gene-TF regulatory network encompassed a total of three CDE-OSRGs, 60 TFs, and 61 pairs of interacting relationships.

Figure 7. Regulatory network between CDE-OSRGs and miRNAs. CDE-OSRGs are represented by red nodes, miRNAs are represented by green nodes, and the lines between them represent interaction pairs.

Figure 8. Regulatory network between CDE-OSRGs and TFs. CDE-OSRGs are represented by red nodes, TFs are represented by yellow nodes, and the lines between them represent interaction pairs.

4. Discussion

MN is an autoimmune glomerular disease, the primary clinical signs of which are edoema, hypertension, proteinuria, hematuria, and so forth [29]. MN develops gradually, with a peak onset age of 40–60 years. The exact cause of MN is unknown. Studies have found that MN may have a certain genetic tendency, and the phenomenon of family aggregation is more obvious [30, 31]. When biological cells are exposed to different stimuli, the body produces high-level, high-activity substances like nitrogen free radicals as well as reactive oxygen species. This imbalances the body’s oxidation–reduction reaction rate, as well as the body’s antioxidant capacity is exceeded by the oxidation rate of various oxides, which damages human tissues [32]. OS is thought to have a significant role in tumour development and a host of other disorders, in addition to being a significant contributor to ageing and illness [33]. Therefore, this study identified core genes of MN based on OSRGs, explored their key biological functions and molecular pathways, and constructed regulatory networks for drugs, miRNAs, and TFs.

87 MN and 11 NC samples made up the dataset GSE108113, which was initially obtained and aggregated for the current investigation in order to create an expression matrix with 24,889 genes. 862 DEGs identified in differential gene expression analysis intersected with 437 OSRGs collected from the GeneCards database, resulting in 28 DE-OSRGs It is suggested that these DE-OSRGs regulate the occurrence and progression of MN. Secondly, enrichment analyses were performed of GO and KEGG for 28 DE-OSRGs to examine their potential roles. After that, 19 DE-OSRGs (disconnected genes deleted) were added to the network of protein–protein interactions (PPIs) and it was determined G6PD, HBB, TP53, KEAP1, and PXDN had the highest connectedness.

G-6-PD originates from red blood cells and catalyzes glucose-6-phosphate to produce NAD-PH, which is a coenzyme of glutathione reductase. Reduced glutathione (GSH) is a necessary condition for maintaining hemoglobin stability and red blood cell membrane integrity. In individuals with G-6-PD deficiency in red blood cells, the red blood cell membrane loses thiol protection and its function is impaired, ultimately leading to hemolysis. The clinical examination of red blood cell G-6-PD is mainly used to diagnose hemolytic anemia. The HBB gene is one of the most important hemoglobin genes in the human body, encoding HBB, a protein that plays a crucial role in oxygen transport. Defects in the HBB gene can cause changes in the morphology of hemoglobin molecules, causing red blood cells to deform under hypoxia, leading to serious consequences such as vascular blockage, ischemia, and hypoxia. TP53 is a gene on chromosome 17 that encodes a tumour suppressor protein that includes transcriptional activation, DNA binding, and oligomeric domains. This gene encodes a tumour suppressor protein that includes transcriptional activation, DNA binding, and oligomeric domains. Encoding proteins regulate the expression of target genes in various cellular stress responses, thereby inducing cell cycle arrest, apoptosis, aging, DNA repair, or metabolic changes. This genetic mutation is associated with various human cancers, including hereditary cancers such as Li Fraumeni syndrome. The selective splicing of this gene and the use of alternating promoters can lead to multiple transcriptome variations and isomers. Other isomers have also been shown to be due to the use of alternating translation start codons from the same transcript. The protein Keap1 encoded by the Keap1 gene is a key molecule that regulates intracellular redox balance. It regulates the stability and activity of Nrf2 by binding to another protein. Nrf2 is a transcription factor that plays a protective role in cells during oxidative stress. When cells face oxidative stress, Keap1 binds to Nrf2 and degrades it, thereby inhibiting the activity of Nrf2. However, when cells are strongly stimulated by oxidative stress, the binding of Keap1 and Nrf2 decreases, and Nrf2 stabilizes and enters the nucleus, binding to the antioxidant stress response element (ARE), initiating the transcription of a series of antioxidant stress response genes, thereby protecting cells from oxidative stress damage. The PXDN gene belongs to the catalase family and is a gene related to extracellular matrix proteins with peroxidase activity. Many studies have found that the PXDN gene is closely related to the occurrence and development of various malignant tumours. This indicates that there may be a synergistic or other related relationship between the occurrence and development of membranous nephropathy and blood system diseases, oxidative stress reactions, cancer, etc.

Further, ROC curve analysis was done, and the G6PD, HBB, TP53, and KEAP1 AUC values and PXDN were 0.804, 0.888, 0.968, 0.936 and 0.821, respectively, indicating that CDE-OSRGs had high diagnostic accuracy in MN patients. Furthermore, TP53 expression levels were positively connected with PXDN and negatively correlated with KEAP1, whereas PXDN expression levels were favourably correlated with HBB and negatively correlated with KEAP1, according to correlation analysis. Nuclear factor erythroid 2-related factor 2 (Nrf2) signalling was enhanced in genetically modified mice, as studied by Rush et al. [34]. They discovered that Nrf2 inhibitors, Kelch-like ECH-associated protein 1 (KEAP1) enhances glomerular sclerosis, renin destruction and exfoliation, podocyte injury, pedal process regression, and interstitial fibrosis. Pessoa et al. [35] found that Spironolactone, a mineralocorticoid receptor blocker, can improve the nephropathy of spontaneously hypertensive diabetic rats by reestablishing glucose-6-phosphate dehydrogenase (G6PD) activity, and reduce tissue damage caused by OS without affecting blood sugar and blood pressure.

Finally, three regulatory networks were constructed to visually describe the interactions between CDE-OSRGs and targeted drugs, miRNAs, and TFs. Although the pharmacological mechanism of Tripterygium wilfordii (T. wilfordii)’s activity is unknown, it has been frequently utilized to treat idiopathic membranous nephropathy (IMN). Shi et al. [36] used network pharmacology to identify genes that overlapped between IMN target genes and T. wilfordii active chemicals. They also discovered that T. wilfordii may potentially target genes like TP53 through PPI networks for therapeutic purposes.

Nevertheless, this study is subject to certain limitations. Firstly, the absence of clinical or basic experiments to validate the extent of DE-OSRGs expression poses a constraint. Secondly, a comprehensive functional assessment is required in future investigations to elucidate the pathophysiological molecular mechanisms of CDE-OSRGs and ascertain the reliability of these five genes. Additionally, the inclusion of a substantial sample size of MN patients and the collection of more comprehensive clinical follow-up data will contribute to corroborating and refining the findings of this study.

Acknowledgments

The authors thank the data provided by the GEO database.

Authors’ contributions

The study was designed by Lihua Wang. The first draft was written by Lihua Wang. Xinyan Feng gathered unprocessed data. Bioinformatics and statistical analysis were carried out by Xinyan Feng. The work was directed by Lihua Wang.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability

You may request access to the dataset that was used and/or analysed for this work by getting in touch with the relevant author.