Potential Biomarkers and Signaling Pathways Associated with the Pathogenesis of Primary Ameloblastoma: A Systems Biology Approach

Objective Ameloblastoma is a benign odontogenic tumor that may lead to ameloblastic carcinoma. This study aimed to determine potential signaling pathways and biological processes, critical genes and their regulating transcription factors (TFs), and miRNAs, as well as protein kinases involved in the etiology of primary ameloblastoma. Methods The dataset GSE132472 was obtained from the GEO database, and multivariate statistical analyses were applied to identify differentially expressed genes (DEGs) in primary ameloblastoma tissues compared to the corresponding normal gingiva samples. A protein-protein interaction (PPI) map was built using the STRING database. The Cytoscape software identified significant modules and the hub genes within the PPI network. Gene Ontology annotation and signaling pathway analyses were executed by employing the DAVID and Reactome databases, respectively. Significant TFs and miRNAs acting on the hub genes were identified using the iRegulon plugin and MiRWalk 2.0 database, respectively. A protein kinase enrichment analysis was conducted using the online Kinase Enrichment Analysis 2 (KEA2) web server. The approved drugs acting on the hub genes were also found. Results A total of 1,629 genes were differentially expressed in primary ameloblastoma (P value <0.01 and |Log2FC| > 1). HRAS, CDK1, MAPK3, ERBB2, COL1A1, CYCS, and BRCA1 demonstrated high degree and betweenness centralities in the PPI network. E2F4 was the most significant TF acting on the hub genes. BTK was the protein kinase significantly enriched by the TFs. Cholesterol biosynthesis was considerably involved in primary ameloblastoma. Conclusions This study provides an intuition into the potential mechanisms involved in the etiology of ameloblastoma.


Introduction
Odontogenic tumors are a class of heterogeneous lesions arising from the tooth-forming tissue, including ectomesenchyme and epithelium remnants associated with the teeth' formation [1,2]. ese tumors could affect individuals of different ages, with peripheral or central tumor locations in the maxillary or mandibular region, leading to facial swelling [3,4]. Ameloblastoma is an odontogenic tumor originating from the cells close to the tooth-root derived from the ectoderm germ layer. Although ameloblastoma is the pathology in 1% of all tumors in the jaw region, they are known as the second commonest odontogenic tumor. ey occur equally in men and women with ages ranging from 30 to 50, with approximately 80% of the cases in the mandible. ey are principally benign tumors; however, they show aggressive behavior and could lead to ameloblastic carcinoma (also known as malignant ameloblastoma). e recurrence rate of ameloblastoma is high if the lesions are not cut off perfectly during surgery [5]. Although several types of research have been executed to elucidate molecular changes in odontogenic tumors such as ameloblastoma, their exact underlying mechanisms, including cellular differentiation and tumorigenesis, are unclear and need more studies in this field [6].
Inside human cells, thousands of genes, transcriptomes, and proteins function in various complicated biological networks, including protein-protein interaction (PPI), gene regulatory, metabolic, and signaling networks [7]. erefore, the conventional gene-by-gene method is not as vigorous to achieve a comprehensive view of cellular action mechanisms.
e microarray technique has been developed to simultaneously measure the complete genome's expression in biosystems; this provides scientists an excellent opportunity to disclose signaling pathways and biological processes dysregulated in several human disorders. By using microarray technology, functional genomics is also achievable, unraveling gene functions. Biomarker discovery is another application of this practical approach, leading to identifying disease-specific markers including diagnostic, predictive, and prognostic molecular markers. All these applications could potentially lead to a better diagnosis and medication of the disease. As a high-throughput technique, analyzing the large amount of data generated by microarray technology needs new and vigorous methods to comprehend biological complexity and therapeutic development accurately. To this end, bioinformatics approaches are recommended, which use computational, statistical, and mathematical algorithms to analyze massive biological datasets [8,9].
In this study, we hypothesized that the significant change in the expression of numerous transcripts in primary ameloblastoma tissues compared to the adjacent healthy gingiva samples leads to an abnormal expression of several proteins, which may lead to the misregulation of miscellaneous pathways and biological processes (BPs) involved in the etiology of ameloblastoma. Moreover, the most critical genes (hub genes), transcription factors (TFs), and protein kinases associated with the pathogenesis of primary ameloblastoma could be illustrated by the protein-protein interaction (PPI) network, gene regulatory network (GRN), and protein kinase enrichment analyses, respectively, therefore, may be assigned as potential biomarkers of primary ameloblastoma. e present study aimed to identify (1) differentially expressed genes (DEGs) in primary ameloblastoma tissues compared to the corresponding normal gingiva samples, (2) hub genes and clusters in the PPI network associated with the primary ameloblastoma, (3) TFs and protein kinases regulating hub genes and TFs, respectively, (4) the most important signaling pathways and BPs enriched in primary ameloblastoma, and (5) approved drugs for possible inhibition or activating the upregulated and downregulated hub genes, respectively.

Microarray Dataset Retrieval.
A gene expression dataset developed by Kondo et al. [10] was analyzed to study the mRNA expression alterations in primary ameloblastoma tumors compared to the corresponding normal oral tissue derived from patients with ameloblastoma. e diagnoses were completed by pathologists at Aichi Medical University Hospital, based on the standard histological classification of odontogenic tumors by the World Health Organization (WHO) [11]. e scalpel was used to obtain fresh normal and tumor samples (0.5-2 cm in diameter) from patients for further experiments including cDNA in microarray techniques, real-time RT-PCR, and the western blotting analysis. All the ameloblastoma tumors were located in the mandible. e total RNA was extracted with DNase treatment using the TRIzolTM reagent ( ermo Fisher Scientific KK) and NucleoSpin RNA (TaKaRa Bio Inc.). e experimental protocol for the cDNA microarray was executed according to the procedure described by Agilent Technologies [12]. Concisely, the Agilent Low Input Quick Amp Labeling Kit (Agilent Technologies) was used for cDNA synthesis and cRNA labeling with the cyanine 3 (Cy3) dye.
e Gene Expression Hybridization kit (Agilent Technologies) was used for purification, fragmentation, and hybridization of Cy3-labeled cRNA on a Human Gene Expression 8 × 60 K v2 Microarray Chip containing 26,740 Entrez Gene RNAs. Kondo et al. [10] submitted the normalized and raw microarray data to the GEO database, namely, GSE132472. e gene expression dataset of GSE132472 [10] was retrieved for subsequent reanalyzing from the Gene Expression Omnibus (GEO), available at https://www.ncbi. nlm.nih.gov/geo/ [13], which is a public access database of microarray profiles.
is dataset was based on the GPL16699 platform (Agilent-039494 SurePrint G3 Human GE v2 8 × 60 K Microarray). It consisted of eight primary ameloblastoma tumors, eight normal oral tissue samples, one human ameloblastoma cell line, and one human fibroblast cell line derived from ameloblastoma patients. is project has been reviewed and approved by the Ethical Committee of Hamadan University of Medical Sciences, Hamadan, Iran (Ethics no. IR.UMSHA.REC.1400.599).

Multivariate Statistical Analyses.
A new dataset was selected from the GSE132472 TXT file, which consisted of eight ameloblastoma tissue samples and eight normal gingiva tissue specimens collected from patients with primary ameloblastoma. Normalization was executed before statistical analysis, and the unsupervised principal component analysis (PCA) was applied to the dataset to detect outlier sample(s) [14]. Subsequently, the supervised Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA) [14] was used to identify DEGs between earlystage ameloblastoma and healthy tissue samples. All multivariate statistical analyses were performed using the"ropls" package from the R programming environment (version 4.0.0; R Core Team, available at www.R-project.org). Genes having an absolute value of Log2 fold change (|Log2 FC|) > 1 and P value <0.01 were considered statistically significant.

PPI Network Construction and Identification of Significant Clusters and Hub Genes.
e interactions between proteins encoded by the DEGs were identified using the online search tool for the Retrieval of Interacting Genes (STRING) database version 11.5, which is available at https://string-db. org/ [15]. e interaction score cutoff was set to 0.4, and the unconnected nodes were removed from the initial graph. Next, Cytoscape 3.8.0, available at https://www.cytoscape. org [16], was used to visualize the PPI network and execute further structural analyses. Clusters were detected using the Molecular Complex Detection (MCODE) plugin [17], and the modules with the following characteristics were considered significant [18]: MCODE score >3, degree cutoff � 2, node score cutoff � 0.2, k-score � 2, and max depth � 100, and the number of nodes ≥10. MCODE is ordinarily used to identify condensed regions in PPI networks, which are assumed to include proteins involved in common signaling pathways and biological processes [19]. e seed node is also detected by the MCODE, which is known as the vertex of each cluster according to its biological function [20]. Furthermore, the network analyzer tool was utilized to calculate the topological features for each node within the PPI network. e hub genes were detected based on their degree and betweenness centralities [19].

Gene Ontology Annotation and Pathway Enrichment
Analyses. Gene Ontology (GO) annotation analyses, including biological process (BP), cellular component (CC), and molecular function (MF) enrichment analysis, were carried out by using the Database for Annotation, Visualization and Integrated Discovery (DAVID) database version 6.8, which is available at https://david.ncifcrf.gov/ [21]. Moreover, pathway enrichment analysis was conducted using the Reactome database version 77, available at https:// reactome.org/ [22]. Reactome is an open-source and manually curated peer-reviewed pathway database. It provides wise bioinformatics tools to illustrate, interpret, and analyze pathways to support clinical and basic studies, genome analysis, modeling, systems biology, and education [22]. Significant modules were considered for pathway enrichment and BP annotation analysis [18][19][20]23]. However, DEGs were used to identify CCs and MFs significantly dysregulated in primary ameloblastoma. e cutoff conditions were set to false discovery rate (FDR) <0.05 and the number of enriched genes ≥2 [23].

GRN Construction.
e two master regulators control the expression levels of genes within the cell, including TFs and miRNAs [24]. Previous studies have indicated that TFs and miRNAs can either induce or suppress the expression of their target genes [25,26]. e binding of miRNA to the promoter can induce the expression of the target gene while binding to 3′ UTR and 5′ UTR, as well as the coding sequence can result in gene silencing [26,27]. e present study built a GRN consisting of the hub genes, miRNAs, and TFs. e iRegulon plugin was used for the detection of potential TFs for the hub genes [18]. In this regard, the TFs with a normalized enrichment score (NES) >3 were considered statistically significant [28]. Furthermore, the validated upstream miRNAs for the hub genes were identified using the MiRWalk 2.0 database, available at https://zmf. umm.uni-heidelberg.de/apps/zmf/mirwalk2/index.html [29]. Only miRNAs with a number of targets ≥10 were selected for constructing the GRN.

Protein Kinases Enrichment
Analysis. Protein kinases and phosphatases play a critical role in posttranslation modification by adding and removing phosphate group to and from proteins. e activity of several essential proteins in the cell depends on the coordination of these two types of enzymes. erefore, the cell's dysregulation of kinases and phosphatases could cause dysregulation of many crucial signaling pathways and biological processes [30]. Previous studies have linked overexpression of most protein kinases to enhanced cell proliferation, survival, carcinogenesis, metastases, and recurrence of different types of cancer [31][32][33]. us, inhibition of protein kinases has become a promising strategy in cancer therapy, which has already provided satisfying effects [34]. e online Kinase Enrichment Analysis 2 (KEA2) web tool, which is available at https://www.maayanlab.net/KEA2/ [35], was employed to identify protein kinase(s) affecting TFs regulating the hub genes [27]. e KEA2 is being made by the Ma'ayan Lab at the Icahn School of Medicine at Mount Sinai, New York, by manually curating interactions between kinases and their targets in the literature. A list of significant TFs with the criteria of NES >3 was considered as an input list in the KEA2 database.

Identification of Approved Drugs Acting on Hub Proteins.
Following the methods of Mahfuz et al. [27], the online DrugBank database (version 5.1.8, released on 2020-01-03; available at https://go.drugbank.com/) [36] was used to identify experimental/investigational approved drugs for inhibiting and activating the proteins encoded by the upregulated and downregulated hub genes in ameloblastoma, respectively. is may accelerate the development of new therapies for ameloblastoma in the future. However, the therapeutic effects of these drugs in ameloblastoma must be studied in vitro and in vivo. is version of DrugBank provides valuable information about 14,583 drug entries including 2,702 approved small molecules, 1,499 approved biologics (peptides, proteins, vaccines, and allergens), 132 nutraceuticals, and over 6,652 experimental drugs (discovery-phase).

Identification of Differentially Expressed Genes in
Ameloblastoma. Following the normalization process, the initial dataset consisted of 16 observations (normal, 8; ameloblastoma, 8) and 58717 variables (probe IDs). According to the PCA plot (R 2 X � 0.549), one of the observations related to the normal tissue samples was identified as an outlier, and therefore, it was removed before further discriminant analysis (Figure 1(a)). e new dataset consisted of 15 observations and was used for OPLS-DA; this predictive model robustly distinguished primary ameloblastoma from healthy controls (R 2 X � 0.308, R 2 Y � 0.942, Q 2 � 0.63) (Figure 1(b)). For multiple probes related to the unique gene, the average value of their expressions was calculated. Finally, a total of 1629 DEGs, including 541 upregulated and 1,088 downregulated genes, with the criteria of P value <0.01 and |Log2FC| > 1, were identified for primary ameloblastoma tissue compared with the normal International Journal of Dentistry gingiva tissue (Supplementary Table 1). Figure 1(c) shows the volcano plot demonstrating the overview of DEGs, which was achieved using the Shiny apps web-based tool, available at https://huygens.science.uva.nl/ [37].

PPI Network Analysis.
A PPI network was built with a total of 1,500 connected proteins and 9,203 edges. A total of 14 significant clusters were identified within the PPI network using the MCODE plugin ( Figure 2). Number 1 was the most considerable condensed region among all modules, with 51 nodes, 1,111 edges, and an MCODE score of 44.44. e characteristics of all modules are presented in Table 1. In addition to the PPI network analysis, the nodes' topological features, including the degree and betweenness centrality of the proteins, were determined using the network analyzer tool. e nodes with the degree and betweenness centrality criteria above the average value of all vertexes within the PPI network were regarded as hub proteins. e average value of the degree and betweenness centrality for the nodes were determined to be 12.27 and 0.00246, respectively. Accordingly, a total of 106 proteins revealed high levels of centrality and were considered hubs (Supplementary Table 2). e heat-map of hub proteins was generated using the R language (version 4.0.2), which is available at https://cran.rproject.org/bin/windows/base/ [38] (Figure 3(a)). Furthermore, the interactions between the hubs were discovered using the STRING knowledge database (Figure 3(b)). e

GO Annotation and Signaling Pathway
Analyses. e Reactome and DAVID knowledge databases revealed that a total of 529 pathways and 78 BPs were significantly affected in ameloblastoma (FDR <0.05).
e "metabolism of steroids," "extracellular matrix organization," "activation of gene expression by SREBF," "regulation of cholesterol biosynthesis by SREBF," "collagen biosynthesis and modifying enzymes," as well as "cell cycle" were the most significant signaling pathways found to be deregulated in ameloblastoma (Figure 4(a)). Moreover, the most significant BPs deregulated in ameloblastoma were as follows: "cholesterol biosynthesis process," "DNA replication," "extracellular matrix organization," "mitotic nuclear division," "collagen catabolic process," and "cell division" (Figure 4(b)). In addition to the results obtained from the DAVID database, it was found that a total of 20CCs and 9MFs were significantly enriched in ameloblastoma. In this regard, "extracellular exosome" (CC), "extracellular matrix"

Identification of Enriched Protein
Kinases. By using the KEA2 web server, it was found that the tyrosine-protein kinase (BTK) could significantly regulate the general transcription factor II-I (GTF2I) at three different phosphorylation sites, including GTF2I_Y248, GTF2I_Y398, and GTF2I_Y503 (adjusted P value � 3.84E − 07).

Identification of Master Regulators Acting on the Hub
Genes. Master regulators, including TFs and miRNAs, affecting the hub genes were identified, and subsequently, a GRN was built. e TFs with an NES of >3 and miRNAs with a total number of target genes of ≥10 were considered for constructing the GRN. is network included 822 edges and 140 nodes, including 105 hubs, 18 miRNAs, and s17TFs. Of note, the BTK was also included in the GRN ( Figure 5). e list of target genes related to the TFs and miRNAs is presented in Tables 2 and 3, respectively. 3.6. Identification of Approved Drugs. By searching the DrugBank database, a total of 42 and 46 approved drugs (or drug-like compounds) were identified that could act on the upregulated and downregulated hub proteins, respectively. In the case of upregulated proteins, 16 drugs were identified to affect CACNA1C, 13 components with PDGFRB, five compounds with NR3C1, four drugs with FGFR1, two compounds with APOE, one compound with MMP14, one compound with SNAP25, one compound with FYN, one drug with PRKCA, one compound with P4HA1, and one compound with RAC2. Regarding the downregulated proteins, a total of 33 drugs could act on ADRB2, nine components with ANXA1, one compound with MAPK3, one compound with NDUFA9, one drug with GOT2, and one compound with MDH2. Table 4 presents the list of identified approved drugs and their associated targets.

Discussion
Ameloblastoma constitutes 1% of all tumors in the oral cavity, with a prevalence of 0.5 per million individuals each year [39,40]. Although it is a benign lesion, it could lead to malignant ameloblastoma with a dismal prognosis [41]. is study performed integrated bioinformatics analyses to disclose potential biological procedures, genes, TFs, miRNAs, and protein kinases triggering ameloblastoma.
erefore, overexpression of a considerable number of ECM-associated genes in ameloblastoma may be involved in the formation and development of the disease. e GO annotation analysis was conducted using the DAVID database to reveal biological processes affected in ameloblastoma. Accordingly, it was found that the most considerable modules in the PPI network were primarily enriched in the cholesterol biosynthesis process, DNA replication, extracellular matrix organization, mitotic nuclear division, and collagen catabolic process. Moreover, the Reactome pathway analysis demonstrated that the most considerable clusters were principally correlated with the metabolism of steroids, extracellular matrix organization, activation of gene expression by SREBF, regulation of cholesterol biosynthesis by SREBF, and cholesterol biosynthesis, as well as collagen biosynthesis and modifying enzymes. e extracellular matrix term was also associated  International Journal of Dentistry with the most significant cellular components and molecular functions affected by the DEGs in ameloblastoma. Deregulation of GO annotations and signaling pathways related to the extracellular matrix and collagen biosynthesis in ameloblastoma is expected due to the upregulation of several ECM-related genes in this disease. erefore, it is suggested that targeting ECM-related genes may result in therapeutic effects in patients with ameloblastoma, although this needs confirmation in the future. Several flavonoids have shown inhibitory effects on MMP13 and MMP8 in previous studies [45,46]. However, the binding affinity of these natural compounds should be examined on other ECM-related genes such as MMP14. Of note, the approved drug named Marimastat inhibits MMP14 [47].
is may be due to the cellular defense mechanism to reduce the formation of Mst1/Mst2 heterodimers, leading to Hippo pathway activation, although this needs confirmation.
Rotellini et al. [55] designed a study to examine the expression of several cancer-related proteins, including cellular tumor antigen p53 (TP53), serine/threonine-protein kinase B-Raf (BRAF), and EGFR, as well as ERBB2 in maxillary ameloblastoma. is was executed by using immunohistochemistry and fluorescent in situ hybridization analysis. According to the study results by Rotellini et al. [55], the ERBB2 protein was negative in all disease samples including primary and metastatic lesions. According to the present results, ERBB2 was significantly downregulated in ameloblastoma tissue compared to the corresponding normal area (FC � 0.43, P value � 0.002).
In our previous report, a positive correlation was observed between CYCS (cytochrome c protein) overexpression and a dismal prognosis in patients with OSCC; we suggested that this may be due to the enhanced tumorigenesis in OSCC patients with poor prognosis [19]. According to the present results, CYCS was downregulated in ameloblastoma tissue compared with the corresponding healthy controls (FC � 0.314; P value � 0.004); this may be a part of a physiopathological mechanism involved in the initiation and/or progression of ameloblastoma, although confirmation is inevitable. e silenced BRCA1 and/or BRCA2 have been associated with the deficiency in the DNA double-strand repair and genomic instability, leading to cancer predisposition [64]. Our findings showed significant underexpression in ameloblastoma tissue and corresponding healthy controls (FC � 0.44; P value � 0.005).
Accumulating evidence indicates that dysregulation of E2F family members could lead to various malignancies, including breast, ovarian, prostate, bladder, colon, and lung adenocarcinoma [65]. According to the present results, E2F4 was the most significant transcription factor enriched by the hub genes, with an NES of 4.951. Moreover, 28 hub genes were found to be regulated by this gene.

10
International Journal of Dentistry radiotherapy-resistant OSCC tissue array and was significantly associated with poor prognosis in OSCC cohorts. Since the exact role of BTK in the etiology of ameloblastoma has not been elucidated, further research is required to demonstrate the exact function of the gene in the disease. Several bioinformatics studies have previously reported the underlying mechanisms of ameloblastoma and potential biomarkers for the disease. Each method and strategy has its advantages, and by integrating all the outcomes from different studies, more reliable results are obtained. Santos et al. [68] performed a valuable bioinformatics analysis to illustrate possible genes involved in ameloblastoma and keratocystic odontogenic tumor (KCOT). In comparison to our study, some of the methods used by Santos et al. [68] were common, and some others were unique. Santos et al. [68] used the GeneCard, available at https://www.genecards.org/ [69] and STRING databases to extract possible genes involved in the pathogenesis of ameloblastoma and to find interactions between the genes, respectively. Furthermore, Santos et al. [68] expanded the PPI network using the STRING database in their study. e authors calculated the sum of the interaction scores of each node in the PPI network and then adjusted the score by multiplying it by 1,000 [70] to achieve a single score named a weighted number of links (WNL). ey called the most important genes with the highest WNL "leader genes." In our study, the PPI network was constructed based on the DEGs obtained from analysis of the microarray technique, which is a highthroughput approach, and our team did not expand the PPI network. Moreover, our study identified the hub genes based on the degree and betweenness of centralities. Furthermore, Santos et al. [68] executed functional enrichment analysis using the STRING database, but we used the DAVID and Reactome databases to perform GO annotation and pathway enrichment analysis, respectively. According to Santos et al. [68] and the present study, CDK1 was considered a critical gene, playing a significant role in the etiology of ameloblastoma. Zhang et al. [71] recently executed an integrated bioinformatics analysis to demonstrate potential biomarkers and molecular mechanisms associated with the epithelialmesenchymal transition (EMT) and immune infiltration in ameloblastoma.
e authors used the CIBERSORT algorithm to study the immune infiltration in ameloblastoma. Zhang et al. [71] reported that a total of 776 genes were deregulated in ameloblastoma, in which FN1 was upregulated and was linked to the macrophage M2 polarization.  e current study had certain limitations. Only ameloblastoma patients from Aichi Medical University Hospital in Japan were included in the GSE132472 dataset. erefore, our results may not fully translate to patients in other countries. Moreover, only 16 tissue samples, including eight normal and eight ameloblastoma tissues, were included in the GSE132472 dataset, and therefore, the sample size that we reanalyzed was small. Generating and analyzing more datasets with a more significant number of samples is recommended for future studies. In addition, these findings should be confirmed using molecular experiments such as qRT-PCR, western blotting, and immunofluorescence. Moreover, the absence of subtypes of ameloblastoma samples was another limitation of our study. Analyzing different histological types of tumors may lead to significant results associated with the pathogenesis of ameloblastoma subtypes.

Conclusion
e present study suggests that a total of 1,629 genes, including 541 up-and 1088 downregulated genes, are differentially expressed in the primary ameloblastoma tissue compared to the corresponding healthy controls. Moreover, 106 hub genes were identified in the PPI network associated with the pathogenesis of the disease. In this regard, HRAS, CDK1, MAPK3, ERBB2, COL1A1, CYCS, and BRCA1 revealed high levels of centralities based on the degree and betweenness values. Furthermore, 17TFs and 18 miRNAs were also found as master regulators of the hub genes. E2F4 was the most significant TF, with the NES of 4.951. BTK was also found to be significantly enriched by the TFs. e "cholesterol biosynthetic process" and "cholesterol biosynthesis" were demonstrated to be in the top-10 ranked BPs and pathways significantly enriched in ameloblastoma based on their FDR value. Besides, "MAPK1/MAPK3 signaling" (R-HSA-5684996), "MAPK family signaling cascades" (R-HSA-5683057), "MAPK6/MAPK4 signaling" (R-HSA-5687128), and "MAPK cascade" (GO:0000165) were also dysregulated in ameloblastoma (Supplementary Tables 3 and  4). Moreover, translating these results from bioinformatics studies to clinical practice requires more research and concentration on designing new diagnostic kits, which are easy to use in clinical laboratories, and discovering novel therapeutic drugs.
Data Availability e datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.
Ethical Approval e current study was approved by the Ethics Committee of Hamadan University of Medical Sciences, Hamadan, Iran (ethics no. IR.UMSHA.REC.1400.599).