Identification of Key Genes and Pathways in Oxaliplatin-Induced Neuropathic Pain Through Bioinformatic Analysis

Abstract

Background

The mechanism of Chemotherapy-induced neuropathic pain (NP) remains obscure. This study was aimed to uncover the key genes as well as protein networks that contribute to Oxaliplatin-induced NP.

Material/Methods

Oxaliplatin frequently results in a type of Chemotherapy-induced NP that is marked by heightened sensitivity to mechanical and cold stimuli, which can lead to intolerance and discontinuation of medication. We investigated whether these different etiologies lead to similar pathological outcomes by targeting shared genetic targets or signaling pathways. Gene expression data were obtained from the Gene Expression Comprehensive Database (GEO) for GSE38038 (representing differential expression in the spinal nerve ligation model rats) and GSE126773 (representing differential expression among the Oxaliplatin-induced NP model rats). Differential gene expression analysis was performed using GEO2R.

Results

Protein-protein interaction (PPI) analysis identified 260 co-differentially expressed genes (co-DEGs). Subsequently, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed three shared pathways involved in both models: Kaposi sarcoma-associated herpesvirus (KSHV) infection, Epstein-Barr virus (EBV) infection, and AGE-RAGE signaling pathway in diabetic complications. Further bioinformatics analysis highlighted eight significantly up-regulated genes in the NP group: Mapk14, Icam1, Cd44, IL6, Cxcr4, Stat1, Casp3 and Fgf2. Our results suggest that immune dysfunction, inflammation-related factors or regulating inflammation factors may also be related to Oxaliplatin-induced NP. Additionally, we analyzed a dataset (GSE145222) involving chronic compression of DRGs (CCD) and control groups. CCD model is a classic model for studying NP. We assessed these hub genes’ expression levels. In contrast with the control groups, the hub genes were up-regulated in CCD groups, the difference was statistically significant, except Stat1.

Conclusion

Our research significantly contributes to elucidating the mechanisms underlying the occurrence as well as the progression of Oxaliplatin-induced NP. We have identified crucial genes and signaling pathways associated with this condition.

Keywords:

• neuropathic pain
• oxaliplatin
• spinal dorsal horn
• bioinformatics
• differentially expressed genes
• hub genes

Introduction

Chemotherapy-induced neuropathic pain (CINP) is a progressive, persistent, and challenging disease characterized by pain, numbness, tingling sensation, and sensitivity to cold. It affects about 50–90% of individuals undergoing chemotherapy treatment.¹ CIPN is mainly induced by first-line chemotherapy drugs such as oxaliplatin, paclitaxel, and vinblastine. With the accumulation of chemotherapy drug toxicity, peripheral neuropathological symptoms become increasingly severe, compelling patients to reduce the chemotherapy drug dosages, shorten the chemotherapy cycle, or even discontinue medication, ultimately impacting patient survival rates. Chemotherapy-induced neuropathic pain (CINP) remains a therapeutic challenge with no effective drugs and treatments currently available.

Oxaliplatin is a widely used chemotherapy-drug for the treatment of colorectal and advanced ovarian cancer.² However, one of the major side effects associated with oxaliplatin treatment is the development of neuropathic pain (NP). This condition is characterized by severe and painful peripheral neuropathy, which can often lead to dose reduction or early discontinuation of chemotherapy.^3,4 Despite its clinical significance, the exact mechanism underlying Oxaliplatin-induced neurotoxicity and the subsequent development of NP remains unclear.

Inflammation of the dorsal root ganglions (DRGs) is the pivotal physiological basis for the occurrence of NP. During the process of pain generation, the DGRs serve as the primary neuron of pain input in the pain mechanism. So the DGRs are considered to be key participants within the pathogenesis that lead to NP.⁵ Studies have found that the expression of neuroinflammatory markers is significantly increased in both the dorsal root ganglion (DRG) and the spinal cord (L4-L5 region) of rats.^6,7 Other studies have found that the main target organ of oxaliplatin’s action is the DGRs,^8,9 its specific mechanism is still unclear. Investigating the shared targets and signaling pathways associated with different etiologies leading to the same outcome can enhance our comprehension of the molecular mechanisms underlying NP.In this study, GSE38038 and GSE126773 datasets obtained from the GEO database were investigated. Combined bioinformatics analyses with enrichment methods, the objective of this study was to find out shared differentially expressed genes (DEGs) and elucidate the functions of them in NP induced by diverse etiologies. Additionally, the construction of a protein-protein interaction (PPI) network was carried out, and the STRING database and Cytoscape software were employed to uncover gene modules and identify central genes as hub genes. Finally, we collected expression data on central genes from CCD model rat induced by chronic compression of DRGs to validate our research findings against the experimental and control groups.

Materials and Methods

Data Source

We conducted a search of the public GEO database (https://www.ncbi.nlm.nih.gov/geo/) and retrieved two GSE datasets (GSE38038 and GSE126773) that pertain to the microarray changes in the dorsal root ganglia (DRGs) of rats. In total, four Oxaliplatin-induced NP samples (GSM3612644, GSM3612645, GSM3612646, GSM3612647) and three normal samples (GSM3612641, GSM3612642, GSM3612643) from GSE126773; four samples of NP induced by spinal nerve ligation (GSM932859, GSM932860, GSM932561, GSM932862), and four sham operation samples (GSM932855, GSM932856, GSM932857, GSM932858) from GSE38038 were selected.

Identification of DEGs

To differentiate DEGs of the diseased and control groups, we utilized the online analysis tool of GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r)¹⁰ on a basis of the R packages of GEOquery as well as Limma, which facilitate data reading and differential expression calculation, respectively. By comparing gene expression profiles across groups, we determined the DEGs associated with the disease condition. Probe sets lacking annotated gene identifiers were excluded from the analysis to ensure the validity and precision of the results, and genes exhibiting multiple probe sets were either eliminated or their expression levels were averaged. Genes were deemed to be significantly DEGs if they met the criteria of an adjusted P-value less than 0.05 and a |logFC (fold change) | of no less than 1.0. To identify common DEGs, a Venn diagram performed using R software.

Functional and Pathway Enrichment Analysis

The Gene Ontology (GO) database provides succinct annotations regarding the characteristics of gene products, encompassing their functional attributes, involvement in biological pathways, as well as cellular localization. Besides, the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway database focuses on the compilation of gene pathway information across distinct species. Results of enrichment analysis were achieved utilizing the Pathview database (https://www.bioinformatics.com.cn).¹¹ Pathview represents an innovative toolkit designed for integrating and visualizing data based on pathways simplifies the process of mapping user-provided data onto relevant pathway graphs. Users are only required to provide their data and data file, after which the tool seamlessly integrates and maps the user data onto the corresponding pathway. It also generates pathway graphs with the mapped data, including the generation of Venn diagrams. A significance threshold of adjusted P < 0.05 was utilized to indicate statistical significance.

PPI Network Construction and Module Analysis

To explore the relationships within relevant proteins, including direct binding interactions and interconnected regulatory pathways both upstream and downstream, the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) (version 11.0)¹² was utilized. The utilization of the STRING database facilitates the construction of a protein-protein interaction (PPI) network, which enables the representation and analysis of complex regulatory connections among proteins. Besides, interactions with a cumulative score more than 0.7 were deemed to have statistical significance. For visualizing this PPI network, we employed Cytoscape 3.7.2 (http://www.cytoscape.org).¹³ To identify key functional modules, we employed Cytoscape’s molecular complex detection technology (MCODE). Selection criteria were were defined as follows: (1) a K-core value of 2; (2) a degree cutoff of 2; (3) a maximum depth of 100, (4) a node score cutoff of 0.2.

Selection and Analysis of Hub Genes

The cytoHubba plug-in incorporated within the Cytoscape software was utilized for the purpose of identifying hub genes. Eight commonly used algorithms, namely MCC, Degree, MNC, Stress, Closeness, BottleNeck, EPC, and Radiality, were utilized to choose these hub genes, followed by the construction of a co-expression network among the hub genes using GeneMANIA, a powerful tool for elucidating intrinsic relationships within sets of genes. The GeneMANIA website (http://www.genemania.org) was employed for this purpose. Mod Enrichr was utilized for enrichment analysis by introducing hub genes. Mod Enrichr is a suite of tools for gene set enrichment analysis(https://maayanlab.cloud/Enrichr/).¹⁴

Hub Genes Validation Study

The expression of identified hub genes was verified in GSE145222, which contained 5 sham groups as well as 4 chronic compression of DRGs model groups. The comparison between groups was conducted using a T-test, a statistical test that assesses the significance of differences (P <0.05) between means.

Ethical Declaration

All data utilized in this study were sourced from publicly available databases. The study does not involve any investigations pertaining to animals or humans.

Results

Identification of DEGs

To identify genes exhibiting differential expression in the DRGs following NP, we retrieved two microarray datasets from the GEO database: GSE126773 (Figure 1A) and GSE38038 (Figure 1B). These datasets provide information regarding gene alterations specifically observed in Rat DRGs after NP. The heap maps of GSE38038 as well as GSE126773 showed the expression levels of DEGs (Figure 1C and D). After performing differential expression analysis on the GSE126773 and GSE38038 datasets, 260 DEGs were identified. These DEGs were found to be overlapping, as illustrated in Figure 2.

Figure 1 Bioinformatics analysis of the DRGs microarray datasets, GSE126773 and GSE38038, after nucleus pulposus (NP). (A) Volcano plot depicting the differentially expressed genes (DEGs) in GSE126773. (B) Volcano plot illustrating the DEGs in GSE38038. (C) Heatmap displaying the top 1000 DEGs in GSE126773. (D) Heatmap showcasing the top 1000 DEGs in GSE38038.

Figure 2 Venn diagram showing the overlapping DEGs between GSE126773 and GSE38038.

Analysis of the Functional Characteristics of Common DEGs

To investigate the biological functions as well as pathways associated with the 260 common DEGs, enrichment analysis was conducted using GO and KEGG Pathway databases. The results of the GO analysis revealed that the enriched biological processes of these genes were primarily in relation with response to peptide and myeloid cell homeostasis (Figure 3A). Additionally, the analysis indicated enrichment in terms of neuron to neuron synapse, postsynaptic density (Figure 3B), growth factor receptor binding, and hormone receptor binding (Figure 3C). Regarding the KEGG Pathway analysis, results indicate significant enrichment of the 260 common DEGs in three key pathways, namely insulin resistance, EGFR tyrosine kinase inhibitor resistance, and toxoplasmosis (Figure 3D).

Figure 3 GO analysis and KEGG pathway analysis of the overlapping DEGs. (A) Enriched biological processes (BP) identified in GO analysis. (B) Cellular component (CC) terms enriched in GO analysis. (C) Molecular function (MF) terms enriched in GO analysis. (D) KEGG pathway enrichment analysis of the integrated DEGs.

PPI Network Construction and Module Analysis

Using Cytoscape, a PPI network was constructed for the common DEGs. The network was generated by incorporating interactions with combined scores greater than 0.7, indicating a high level of confidence in the reliability of the interactions. Through the utilization of the MCODE plug-in within Cytoscape, six gene clustering modules were identified (Figure 4A). Among these modules, a total of twenty-five genes were found to be present in all of them, indicating their potential functional significance. Notably, the three modules with the highest scores encompassed sixteen genes, suggesting their potential importance in the biological context (Figure 4B-D).The GO analysis revealed that the identified genes are significantly associated with several biological processes and molecular functions, including sterol biosynthetic process, perinuclear endoplasmic reticulum localization, 3-keto sterol reductase activity, type 5 metabotropic glutamate receptor binding, CD4 receptor binding, and oxidoreductase activity. Furthermore, the oxidoreductase activity was characterized by the process of oxidizing a pair of donors, leading to the reduction of molecular oxygen as well as the production of two molecules of water (Figure 5A). According to the KEGG pathway analysis, the identified genes were found to be primarily associated with three key pathways: steroid biosynthesis, Epstein-Barr virus (EBV) infection, and the AGE-RAGE signaling pathway in diabetic complications (Figure 5B).

Figure 4 Significant gene clustering module. (A) Six gene clustering modules are shown. (B-D) The three modules with the highest scores include a total of 16 genes.

Figure 5 Common enrichment analysis results, and KEGG pathway graph. (A) Gene Ontology enrichment analysis results are presented. (B) KEGG pathway enrichment analysis results are displayed.

Hub Gene Selection and Analysis

Through the application of the cytoHubba plug-in and employing six algorithms, the top 10 hub genes were identified, as presented in Table 1. Through the analysis of the Venn diagrams, we have identified ten hub genes that are commonly shared, namely Stat1, II6, Stat3, Cd44, Icam1, Casp3, Fgf2, Jak2, Cxcr4, and Mapk14 (Figure 6A). In order to gain a deeper understanding of the 10 identified hub genes, we utilized the GeneMANIA database to conduct further analysis of their co-expression network and associated functions. The analysis revealed a complex PPI network, with physical interactions accounting for 84.50%, shared protein domains for 6.34%, co-expression for 30.45%, co-localization for 4.01%, and predicted interactions for 5.61% (Figure 6B). Function analysis conducted using GeneMANIA demonstrated that these genes are primarily involved in regulating endothelial cell migration, blood vessel endothelial cell migration, myeloid cell homeostasis, and response to cytokine stimulus regulation, among other functions. Notably, eight genes were found to appear in the three modules with the highest scores, indicating their potential significance (Figure 6C). Additionally, these genes exhibited significant expression differences when compared the control group to the treatment group (Figure 6D). For a comprehensive understanding, Table 2 provides the complete names of these genes along with their associated functions. Furthermore, KEGG pathway analysis of these 10 hub genes has revealed their strong association with Kaposi sarcoma-associated herpesvirus (KSHV) infection, EBV infection, and the AGE-RAGE signaling pathway in diabetic complications (Figure 6E and F).

Table 1 The Top 10 Hub Genes by Topological Analysis Methods of CytoHubba

MCC	MNC	Radiality	Closeness	Degree	EPC
Stat3	Stat3	Stat3	Stat3	Stat3	Stat3
Cd44	Casp3	Stat1	Fgf2	Fgf2	Casp3
Cxcr4	Jak2	Cd44	Stat1	Casp3	Cd44
Fgf2	Fgf2	Mapk14	Cd44	Mapk14	Fgf2
Icam1	Cd44	Fgf2	Mapk14	Stat1	Cxcr4
Mapk14	Mapk14	Jak2	Casp3	Cd44	Jak2
Stat1	Cxcr4	Icam1	Jak2	Jak2	Mapk14
Casp3	Icam1	Casp3	Cxcr4	Cxcr4	Icam1
Jak2	Stat1	Cxcr4	Icam1	Icam1	Stat1
II6	II6	II6	II6	II6	II6

Table 2 Details of the Hub genes

NO.	Gene Symbol	Full Name	Function
1	Mapk14	Mitogen-Activated Protein Kinase 14	Serine/threonine kinase which acts as an essential component of the MAP kinase signal transduction pathway.
2	Icam1	Intercellular Adhesion Molecule 1	ICAM proteins are ligands for the leukocyte adhesion protein LFA-1 (integrin alpha-L/beta-2).
3	Cd44	CD44 Molecule	Cell-surface receptor that plays a role in cell-cell interactions, cell adhesion and migration, helping them to sense and respond to changes in the tissue microenvironment.
4	IL6	Interleukin 6	Cytokine with a wide variety of biological functions in immunity, tissue regeneration, and metabolism.
5	Cxcr4	C-X-C Motif Chemokine Receptor 4	Receptor for the C-X-C chemokine CXCL12/SDF-1 that transduces a signal by increasing intracellular calcium ion levels and enhancing MAPK1/MAPK3 activation.
6	Stat1	Signal Transducer And Activator Of Transcription 1	Signal transducer and transcription activator that mediates cellular responses to interferons (IFNs), cytokine KITLG/SCF and other cytokines and other growth factors.
7	Casp3	Caspase 3	Thiol protease that acts as a major effector caspase involved in the execution phase of apoptosis.
8	Fgf2	Fibroblast Growth Factor 2	Acts as a ligand for FGFR1, FGFR2, FGFR3 and FGFR4.

Figure 6 Venn diagram and co-expression network of hub genes. (A) The Venn diagram illustrates the overlapping hub genes identified by six different algorithms. (B) The co-expression network of hub genes and their associated genes was analyzed using GeneMANIA. (C) Module analysis revealed that eight genes were present in the three modules with the highest scores. (D) Expression patterns of the eight hub genes in the control and treatment groups. (E and F) Pathways associated with the hub genes.

Validation of Hub Genes Expression

To assess the reliability of the expression levels of the identified hub genes (Mapk14, Icam1, Cd44, IL6, Cxcr4, Stat1, Casp3, and Fgf2), we selected a dataset (GSE145222) that includes both sham and chronic compression of DRGs (CCD) groups. The CCD model is a well-established model used to study nucleus pulposus (NP). We conducted an analysis of the expression levels of the aforementioned hub genes in the selected dataset. The results demonstrated that, in comparison to the control group, all of the hub genes showed up-regulation in the CCD groups, with statistically significant differences observed, except for Stat1 (Figure 7). These findings serve as additional evidence to support the potential involvement of these hub genes in the context of NP.

Figure 7 Expression level of hub genes in GSE145222. (A-H) Comparison between data sets using a mean t test. P<0.05 was considered statistically significant. (I) Expression of hub genes in both the control and CCD groups.

Discussion

Chemotherapy-induced peripheral neuropathic pain (CIPN) is a common and debilitating adverse reaction associated with anticancer medications. Oxaliplatin, an anticancer drug derived from platinum and utilized in advanced colorectal cancer, frequently induces a distinct type of CIPN marked by increased sensitivity to mechanical / cold stimuli. However, the mechanisms underlying the development and chronicization of Chemotherapy-induced peripheral neuropathy, as well as effective treatment options, remain poorly understood, highlighting the urgent need for therapeutic advancements in this area. In this study, our findings shed light on potential pathways implicated in the development of CIPN. Specifically, we identified KSHV infection, EBV infection, as well as the AGE-RAGE signaling pathway in diabetic complications as potentially crucial pathways in both animal models of peripheral neuropathy induced by spinal nerve ligation and oxaliplatin treatment. Moreover, we found that the genes Mapk14, Icam1, Cd44, IL6, Cxcr4, Stat1, Casp3, and Fgf2 may play significant roles during the development of CIPN.

Kaposi sarcoma-associated herpesvirus (KSHV) is known to cause various human malignancies and hyperproliferative diseases, particularly in individuals with compromised immune systems like those infected with human immunodeficiency virus (HIV). KSHV has the ability to modulate innate immune pathways, enabling infected cells to survive for extended periods following primary infection and during viral latency.¹⁵ EBV primarily infects B cells, remaining dormant in the body until certain conditions disrupt the EBV-host balance and allow the virus to manifest its pathogenic potential. Studies have suggested that EBV can induce immune dysfunction in susceptible individuals, leading to neuroinflammation through autoimmunity or antiviral immune responses.¹⁶ Diabetes is associated with increased production of advanced glycation end products (AGEs) resulting from reactive dicarbonyl compounds in a hyperglycemic environment. AGEs can trigger the expression of pro-inflammatory cytokines through their main receptor, receptor for advanced glycation end products (RAGE). The involvement of the AGE-RAGE signaling pathway has been suggested in diabetic peripheral neuropathy’s pathogenesis. Several studies have shown a connection between diabetes-related peripheral neuropathy and the AGE-RAGE signaling pathway in diabetic complications.^17–19 In the context of our study, we found that immune dysfunction and the AGE-RAGE signaling pathway in diabetes complications may also play a role in the development of oxaliplatin-induced peripheral neuropathy.

In our study, we observed significant upregulation of genes such as Mapk14, Icam1, CD44, IL6, Cxcr4, Stat1, Casp3, and Fgf2 in both Oxaliplatin-induced peripheral neuropathy (NP) and spinal nerve ligation-induced NP animal models compared to the control group. Mapk14, a member of the MAP kinase family, might be activated by various environmental pressures as well as pro-inflammatory cytokines. Previous research on microglia and NP has shown that Mapk14 is significantly upregulated and plays a role in microglial activation in NP animal models.²⁰ Icam1 is a glycoprotein located on the cell surface and acts as an adhesion receptor. It plays a crucial role in regulating the recruitment of white blood cells to areas of inflammation. It participates in important physiological processes such as cell signal transduction, immune response, and inflammatory response. In NP animal models, Icam1 has been found to be significantly increased compared to the control group.²¹ Additionally, another study found upregulated expression of Icam1 in the spinal cord of rat models of NP compared to rat models of inflammatory pain.²² CD44 is a transmembrane glycoprotein that is widely expressed in various cell types. It is involved in various biological processes including cell adhesion, migration, signal transduction, proliferation, and tumor metastasis. Recent studies focusing on CD44 and CIPN have demonstrated that CD44 acts as a mediator of HMWH-induced pain relief in cases of oxaliplatin and paclitaxel-induced CIPN.²³ These findings emphasize the significant role of CD44 signaling in HMWH-induced antihyperalgesia and establish it as a potential therapeutic target for inflammatory conditions and NP.²⁴ Interleukin-6 (IL-6) is an interleukin with dual functionality, serving as both a pro-inflammatory cytokine and an anti-inflammatory myokine. The inflammatory response, involving the propagation of inflammation and recruitment of neutrophils, is initiated when glial cells are activated by noxious stimuli and inflammation. Pro-inflammatory cytokines, including IL-6, play a crucial role in the initiation and maintenance of NP.²⁵ In a study, it was found that the expression levels of the inflammatory factor IL-6 in the ipsilateral L4-L6 spinal dorsal horn were increased in contrast with the contralateral side at 7 days following nerve crush injury.²⁶

Cxcr4 is a specific receptor for the chemokine stromal cell-derived factor-1 (SDF-1) and exerts a strong effect on the promotion of chemotaxis, particularly for lymphocytes. It has been suggested that Cxcr4 is closely related to the occurrence of NP. Research has shown that Cxcr4 expression increases in spinal glial cells of mice with peripheral nerve injury-induced NP. Blocking Cxcr4 can alleviate pain behavior, while overexpressing Cxcr4 can induce pain hypersensitivity.²⁷ Stat1 belongs to the STAT family, which is critical for regulating gene expression within cells and important in immune responses. In the context of chronic NP, a study identified specific blood biomarkers for chronic NP by comparing patients with chronic pain (neuropathic and nociceptive) to painless controls, and Stat1 was among the identified biomarkers.²⁸ Casp3, a downstream effector protease in the apoptosis cascade, plays a crucial role in cell apoptosis. In a study, IIK-7 was found to significantly alleviate mechanical allodynia and glial activation while inhibiting casp-3 proteins. This suggests that IIK-7 reduces NP by inhibiting glial activation and suppressing proteins associated with inflammation as well as apoptosis.²⁹ Fgf-2 belongs to the fibroblast cytokine family, which is significant for various processes such as tendon-to-bone healing, cartilage repair, bone repair, and nerve regeneration. Research has shown that spinal cord astrocytes upregulate Fgf-2, a neurotrophic and gliogenic factor, in response to the ligation of spinal nerves L5 and L6. Studies have revealed that endogenous astroglial Fgf-2 plays a role in sustaining neuropathic pain (NP) tactile allodynia, which is linked to the reactivity of spinal cord astrocytes. Inhibiting spinal Fgf-2 has been shown to ameliorate NP symptoms.³⁰ These findings suggest that factors involved in inflammation and regulation of inflammation may also play a role in oxaliplatin-induced NP. Understanding the involvement of these factors could provide insights into potential therapeutic targets for managing NP.

The objective of this study was to utilize bioinformatics technology to identify common factors in oxaliplatin-induced NP animal models and spinal nerve ligation-induced NP animal models. By identifying DEGs and hub genes, the study aimed to enhance our understanding of the underlying mechanisms of oxaliplatin-induced NP. However, there were several limitations to consider in this study. Firstly, the study design is retrospective, which means that the findings should be validated using external sources to ensure their reliability and reproducibility. Additional studies with independent datasets are necessary to confirm the identified DEGs and hub genes. Secondly, further investigations are needed to validate the functional roles of the identified hub genes using in vitro models. Experimental studies should be conducted to explore the mechanisms by which these hub genes contribute to oxaliplatin-induced NP.

Conclusions

In conclusion, our study utilized bioinformatics analysis to identify common DEGs in two independent datasets related to oxaliplatin-induced NP and spinal nerve ligation-induced NP. Through enrichment and PPI network analysis, we discovered several hub genes that are potentially involved in the shared mechanisms underlying these two types of NP.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Disclosure

The authors declare that there are no conflicts of interest in this work.

Data Sharing Statement

The data that support the findings of this study are downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi; GEO accession: GSE38038; GSE126773).