Integrated bioinformatics analysis of validated and circulating miRNAs in ovarian cancer
Article information
Abstract
Recent studies have focused on the early detection of ovarian cancer (OC) using tumor materials by liquid biopsy. The mechanisms of microRNAs (miRNAs) to impact OC and signaling pathways are still unknown. This study aims to reliably perform functional analysis of previously validated circulating miRNAs' target genes by using pathfindR. Also, overall survival and pathological stage analyses were evaluated with miRNAs' target genes which are common in the The Cancer Genome Atlas and GTEx datasets. Our previous studies have validated three downregulated miRNAs (hsa-miR-885-5p, hsa-miR-1909-5p, and hsa-let7d-3p) having a diagnostic value in OC patients' sera, with high-throughput techniques. The predicted target genes of these miRNAs were retrieved from the miRDB database (v6.0). Active-subnetwork-oriented Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was conducted by pathfindR using the target genes. Enrichment of KEGG pathways assessed by the analysis of pathfindR indicated that 24 pathways were related to the target genes. Ubiquitin-mediated proteolysis, spliceosome and Notch signaling pathway were the top three pathways with the lowest p-values (p < 0.001). Ninety-three common genes were found to be differentially expressed (p < 0.05) in the datasets. No significant genes were found to be significant in the analysis of overall survival analyses, but 24 genes were found to be significant with pathological stages analysis (p < 0.05). The findings of our study provide in-silico evidence that validated circulating miRNAs' target genes and enriched pathways are related to OC and have potential roles in theranostics applications. Further experimental investigations are required to validate our results which will ultimately provide a new perspective for translational applications in OC management.
Introduction
Despite the numerous research and clinical studies, ovarian cancer (OC) still has a high mortality rate among gynecological cancers since lack of effective biomarkers for early detection and prognosis [1-3]. Also, the treatment efficiency of OC is low depending on multiple challenges such as late diagnosis, lack of reliable markers, the development of resistance to current therapeutics, and phenotype of heterogeneity [1,4]. For these challenges of OC, new studies need to reveal underlying molecular mechanisms, and to discover molecular biomarkers for early diagnosis, prevention, and targeted therapy [4].
MicroRNAs (miRNAs) are small non-coding RNAs that are about ~22 nucleotides in length. Their function is transcriptional and post-transcriptional regulation of gene expression by targeting mRNAs [5,6]. According to the effect on cancer, there are two types of miRNAs which are tumor suppressor miRNAs and oncomiRs. Depending on the type of cancer, oncomiRs or tumor suppressor miRNAs are inhibited or stimulated, respectively [7]. Especially, dysregulations of specific miRNAs affect cancer cell proliferation, differentiation, metastasis, and recurrence formation [8,9]. Various miRNAs have been shown to play different roles in OC. The mechanisms of miRNAs to impact OC and signaling pathways are still unknown. [2,10].
Many experimental studies have verified the detected interactions with bioinformatics analysis and proved the accuracy and predictivity with in-silico tools or databases. Although significant progress has been made about in-silico analysis in evaluating miRNAs, the need for new tools/databases is increasing day by day. It is also because of the limitations in existing tools/databases, that has increased with the development of high-throughput miRNAs technologies to analyze miRNA [11,12].
In this study, pathfindR was used for enrichment analysis of target genes. PathfindR uses active subnetworks, where an active subnetwork can be defined as a subnetwork of interconnected genes in a protein-protein interaction network (PIN), predominantly consisting of significantly altered genes. The tool initially maps the input genes with significance values onto the PIN and identifies active subnetworks, then it performs enrichment analysis on the identified subnetwork gene sets. In general, enrichment approaches overlook the relational information captured in the PIN and the genes neighboring the significant genes are not considered. By identifying active subnetworks, pathfindR exploits interaction information to enhance enrichment analysis. Active subnetworks allow the inclusion of possibly relevant genes that are not significant but connect significant genes in the PIN, and, in turn, the identification of phenotype-associated connected significant subnetworks [13]. This aids pathfindR uncover relevant mechanisms underlying the studied disease/phenotype.
In this study, we developed an in-silico approach to evaluate target genes of OC-related circulating and previously validated miRNAs which may be targeted for therapeutic approaches and utilized for OC management and diagnosis (Fig. 1).
The workflow for bioinformatics analysis of miRNAs. GEPIA, Gene Expression Profiling Interactive Analysis; qPCR, quantitative polymerase chain reaction.
Methods
Identification of the targets genes of validated miRNAs
In our functional analyses, we used hsa-miR-885-5p, hsa-miR-1909-5p, and hsa-let7d-3p which were previously defined by our group as the dysregulated miRNAs that can be a candidate biomarker for OC. These candidate miRNAs, which were determined by microarray, were validated by quantitative polymerase chain reaction (qPCR). Both microarray and qPCR results showed that these three miRNAs were downregulated in the OC group compared with healthy individuals [14,15].
Target genes of validated circulating miRNAs were examined by using the miRDB online database (http://www.mirdb.org). As one of the miRNA-target predictions and the functional annotations databases, miRDB provides access to miRNA-target genes and functions of five different species: human, mouse, rat, dog, chicken. All targets in the database were acquired from the MirTarget database, which was a bioinformatics tool developed with analyzing thousands of miRNA-target interactions obtained from high-throughput techniques. By integrating target prediction and gene ontology enrichment analyses, miRDB presents a streamlined pipeline for quickly identifying miRNA functions [16].
Functional and pathway enrichment analyses
Using all (i.e., union of) target genes of the validated miRNAs, active-subnetwork-oriented Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using pathfindR. The tool pathfindR identifies gene sets that form active subnetworks in a PIN using a list of genes. Afterwards, it performs pathway enrichment analysis. An active subnetwork can be defined as a group of interconnected genes in a PIN that predominantly consists of significantly altered genes. Active subnetworks define distinct disease-associated sets of interacting genes. By incorporating interaction information, pathfindR yields more relevant enrichment results.
For assigning a significance value for each target gene (for use with pathfindR), initially, all Homo sapiens miRNA-target gene scores were obtained from miRDB (v6.0). The significance for each target gene was defined as the probability of observing a score greater than or equal to the score of this target gene over all H. sapiens miRNA-target gene scores [i.e., P(x ≥ observed score)]. For genes that are targeted by more than one miRNA, the lowest significance was kept. The final list of target genes-significance values was, then, filtered keeping genes with significance ≤ 0.5 (corresponds to a score of 67). The significance value used for pathfindR indicates what proportion of all scores (across all H. sapiens miRNA-target gene scores) was as high or higher than the observed score for a given miRNA-target gene pair. The threshold proxy significance value of 0.5 was an ad-hoc choice corresponding to a score of 67, which was roughly equivalent to the threshold value suggested by miRDB for obtaining high-confidence miRNA-target gene pairs.
Determination of differentially expressed target genes
The expressions of the target genes were enriched in The Cancer Genome Atlas and Genotype-Tissue Expression (GTEx) datasets by using GEPIA (Gene Expression Profiling Interactive Analysis) (http://gepia.cancer-pku.cn). GEPIA provides customizable functions, such as tumor and/or normal differential expression analysis, profiling according to cancer types or pathological stages, patient survival analysis, detection of gene expression similarities, correlation analysis and dimensionality reduction analysis using RNA sequencing (RNA-seq) data of The Cancer Genome Atlas (TCGA) and GTEx projects [17].
All target genes enriched were analyzed separately by using the Expression DIY feature and performed for the RNA-seq data of TCGA and GTEx datasets by using the ANOVA test. For each gene expression features of box plots were set as |Log2FC| cutoff: 1 and p-value cutoff: 0,01, jitter size: 0,4 and log2 (TPM + 1) for log-scale. Genes with a p-value <0.05 were identified as differentially expressed target genes.
Overall survival and pathological stage analyses
The associations between expression signatures of shared genes and overall survival and pathological stage analyses were performed in the TCGA and GTEx dataset by the GEPIA platform. For overall survival plot analysis used log-rank test, also known as the Mantel-Cox test, for the hypothesis test and cohorts' thresholds adjusted. The Cox proportional hazard ratio and the 95% confidence interval information are included in the survival plot. The method for differentially expressed target gene analysis is one-way ANOVA, using the pathological stage as variable for calculating differential expression. The expression data are first log2(TPM+1) transformed for differential analysis. They included further statistical analysis, Benjamini and Hochberg's false discovery rate (FDR) adjusted p-value (q-value) < 0.05 was identified as statistically significant. Genes with FDR p-value (q-value) < 0.05 were identified as significant association with overall survival or pathological stage.
Results
Target genes analysis of miRNAs
Target genes of validated circulating miRNAs were comprehensively analyzed by miRDB database. Target genes were retrieved separately determined for each miRNA. No specific filters were applied for target gene prediction in miRDB database. Respectively, 422, 230, and 44 target genes were identified with mining target genes of hsa-miR-885-5p, hsa-miR-1909-5p, and hsa-let-7d-3p. Target genes were filtered by a score of 67 for significance values. After filtering, the number of target genes detected for hsa-miR-885-5p, hsa-miR-1909-5p, and hsa-let-7d-3p were respectively 229, 90, and 23 (Table 1). While two common target genes were found in hsa-miR-885-5p and hsa-miR-1909-5p: ERICH3 (glutamate rich 3) and CAPRIN1 (cell cycle associated protein 1), one common target gene were observed in hsa-miR-885-5p and hsa-let-7d-3p: SEC24D (SEC24 homolog D, COPII coat complex component). No common target gene was found for all three miRNAs. The numbers of target and common genes are shown in Fig. 2. Venn diagrams were drawn by using VENNY 2.1 (https://bioinfogp.cnb.csic.es/tools/venny) online web tool.
Target genes determined in the miRDB database (score ≥ 67, for representation purposes)
Venn diagram of common and overlapping target genes identified for the three miRNAs (score ≥ 67).
Pathway enrichment analyses of target genes of miRNAs
Performing active-subnetwork-oriented enrichment analysis via pathfindR using the union of miRNA-target genes (339 genes in total), 24 enriched KEGG pathways were identified (Table 2). The top three pathways were ubiquitin-mediated proteolysis, spliceosome, and Notch signaling pathway according to p-values. The pathway diagrams for these pathways, in which the target genes of our dysregulated miRNAs are indicated by orange color, are presented in Fig. 3. It is found that a total of ten different target genes were involved in these pathways. A part of the target genes which are TRAF6 (TNF receptor associated factor 6), UBE2N (ubiquitin conjugating enzyme E2 N), UBA7 (ubiquitin like modifier activating enzyme 7), UBE2K (ubiquitin conjugating enzyme E2 K), WWP1 (WW domain containing E3 ubiquitin protein ligase 1) are ubiquitin-mediated proteolysis (Fig. 3A) whereas, SRSF2 (serine and arginine rich splicing factor 2), SRSF6 (serine and arginine rich splicing factor 6), and U2SURP (U2 SnRNP associated SURP domain containing) are involved in spliceosome pathway (Fig. 3B). There are only two genes in the Notch pathway: ATXN1L (ataxin 1 like) and KAT2B (lysine acetyltransferase 2B) (Fig. 3C). While nine of these ten genes are (UBE2N, UBE2K, WWP1, TRAF6, U2SURP, SRSF2, SRSF6 KAT2B, and ATXN1L) targeted by hsa-miR-885-5p, only the UBA7 gene is targeted by hsa-miR-1909-5p. There are no target genes of hsa-let-7d-3p associated within three pathways.
Enriched KEGG pathways identified via pathfindR analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway diagrams of ubiquitin-mediated proteolysis (A), spliceosome (B), and Notch signaling pathway (C). Genes targeted by at least one validated miRNA are indicated in orange color.
Analysis of differentially expressed target genes
Differentially expressed target gene analysis was conducted to examine the expressions of identified target genes of miRNAs in the TCGA and GTEx datasets. In union of miRNA-target genes (339 genes in total), 93 genes were found to be differentially expressed (p < 0.05) in TCGA and GTEx datasets of when compared to cancer tissues with paired normal tissues (Table 3). None of the common genes (ERICH, CAPRIN1, and SEC24D) among target genes were significant in TCGA and GTEx datasets. While miR-885-5p targets 64 out of 93 genes, miR-1909-5p regulates 20 and let-7d-3p 9 target genes.
Overall survival and pathological stage analysis of the differentially expressed target genes regulated by circulating validated miRNAs in TCGA and GTEx datasets
Overall survival and pathological stage analyses
The associations with overall survival outcomes and pathological stage analysis for expression signatures of 93 genes in the TCGA and GTEx datasets were performed by the GEPIA (Table 3). After FDR tests, no significant genes were found to be significant in the analysis of overall survival analyses, but 24 genes were found to be significant with pathological stages analyses (p < 0.05). No significant differentially expressed target genes were found to be with overall survival after statistical analyses. However, 24 differentially expressed target genes were significantly associated pathological stage of OC (Fig. 4). Of 24 genes, ZNF407 (zinc finger protein 407) and UBN2 (ubinuclein 2) genes had the lowest p-value for their relationship with pathological stages for the TCGA and GTEx datasets (p = 0.017856) (Table 3).
Pathological stage analysis of differently expressed target genes in the ovarian cancer dataset using the GEPIA (Gene Expression Profiling Interactive Analysis) platform (p < 0.05). Pr(>F) values were adjusted to Benjamini and Hochberg's false discovery rate test.
Discusssion
Recently, a considerable number of research have been made to determine the effect of miRNAs regulate cancer hallmarks and to develop the early diagnosis and prognosis of cancer [18]. Nevertheless, several challenges such as various sampling methods, sample size, detection techniques, gender, and ethnicity or genetic background, affect the reliance on utilizing miRNAs as biomarkers [1,3,10,19]. One of the fundamental research areas for easily detectable, non-invasive, sensitive, and specific miRNA-based biomarker discoveries are of a great value for accurate and effective early diagnosis, risk prediction, prognosis, recurrence, and effective management of OC [8,10,12].
Identification of miRNAs' target genes has been the focus of computational biology for the last few years. Since detecting all possible miRNA targets with high-throughput technologies is laborious and costly, a wide variety of computational resources has been developed. The methodologies used by databases range from evolutionary conservation evaluations of putative miRNA binding sites to machine learning and classification algorithms. Continuous improvement is needed to develop new tools/databases to accurate predictions of miRNA targets [20,21]. One of these databases is miRDB online database and determines miRNA-target prediction and functional annotations. Estimating that 3.5 million target genes were regulated by 7,000 miRNAs in human, mouse, rat, dog, and chicken, in the latest version of miRDB (v6.0) major updated in 2019. In miRDB, miRNA binding and target down-regulation features used to predict miRNA targets with machine learning methods, generates a prediction scores are in the range of 0–100, and candidate genes with scores ≥ 50 are presented as predicted miRNA targets [16].
The pathfindR tool was used for the enrichment analysis of target genes within the scope of the study. While pathfindR allows better identification of disease-related pathways, it should be noted that the tool requires a significance level per each input gene. To overcome this limitation, we followed an approach where we used the scores for miRNA-target gene pairs to calculate a significance level. Additionally, some genes (usually a small proportion) that are not in the PIN (that do not have any curated interactions) have to be discarded, resulting in the loss of possibly relevant genes.
In our previous studies, three downregulated circulating miRNAs (hsa-miR-1909-5p, hsa-miR-885-5p, and hsa-let-7d-3p) in OC patients were consistently validated by comparison with healthy individuals. In these studies, circulating miRNAs, which were determined to be dysregulated by microarray, were then validated with the qPCR [14,15]. In this study, we predicted the validated circulating miRNAs' target genes with miRDB database and performed functional analysis by pathfindR tool to understand the pathogenesis of OC especially the molecular mechanisms of its development. Also, overall survival and pathological stage analyses were evaluated with differentially miRNAs' target genes which are commonly found in the TCGA and GTEx datasets.
Two hundred and twenty-nine, 90, and 23 of target genes of hsa-miR-1909-5p, hsa-miR-885-5p, and hsa-let-7d-3p were separately determined via miRDB (score ≥ 67), respectively (Fig. 2). Enriched KEGG pathways of target genes were performed with pathfindR tool. Totally, 339 genes included in active-subnetwork-oriented enrichment analysis and as a result, 24 enriched KEGG pathways were identified (Table 2). Top three pathways with the lowest p-values are ubiquitin-mediated proteolysis (p < 0.001), spliceosome (p < 0.001), and Notch signaling pathway (p < 0.001), respectively (Fig. 3). TRAF6, UBA7, UBE2K, UBE2N, and WWP1 genes regulate the ubiquitin-mediated proteolysis pathway. While SRSF2, SRSF6, and U2SURP are involved in the spliceosome, ATXN1L and KAT2B are adjusted to the Notch pathway. Although hsa-miR-885-5p targets ATXN1L, KAT2B, SRSF2, SRSF6, TRAF6, UBE2K, UBE2N, U2SURP, and WWP1 genes, hsa-miR-1909-5p just targets UBA7 gene. No target genes of hsa-let-7d-3p were found to be associated with these three pathways.
Our pathway analysis determined different pathways that can be related to the OC development. Of these pathways, only ubiquitin-mediated proteolysis, spliceosome, and Notch signaling pathway showed that their p-value is less than 0.001. Besides the effects of these pathways on OC, they also have roles on other cancer types. Firstly, ubiquitin-mediated proteolysis is an essential mechanism that is responsible for 80%–90% of intracellular protein degradation and is involved in many cellular processes, including tumorigenesis, tumor survival and apoptosis [22,23]. Ubiquitin-mediated pathway modulates BRCA1/2, p53, ERBB2 gene expressions, ERK pathway, cyclin-dependent cell cycle regulation which are related to OC [24]. Bazzaro et al. [25] claimed that upregulation of proteasome subunit levels occurs in OC and proteasome inhibitors may have utility in the treatment of OC. Secondly, splicing mechanism is a crucial process that regulates cellular proliferation, differentiation, and survival [26]. Dysregulation of splicing processes and splicing factor genes contribute to cancer, including OC [27,28]. Especially, splicing machinery mutations especially contribute to tumorigenesis. Additionally, it is critical to understand that the molecular mechanism of RNA splicing is causing the development of drug resistance in cancer treatment [26]. Regulation of the relationship between splicing and cancer can be led to splicing-based therapies for cancer treatment [29]. Thirdly, Notch signaling pathway regulates not only cell self renewal and differentiation but also a cell to cell communication [30]. Upregulation of Notch signaling pathway proteins have been identified in OC [31]. Several studies have revealed that it is related to poor overall and disease free survival time, and more advanced stages [30,32]. The Notch signaling pathway plays a specific role in deregulation of signaling cascade has been associated with OC. Targeted therapy against Notch pathway activation can offer clinical benefit to OC [30,31].
According to these studies, our findings also indicate that our dysregulated miRNAs are significant in OC via pathway related-target genes. Briefly, circulating downregulated miRNAs could not suppress their target genes in these pathways, hence, have activated these pathways, resulting in the development of OC. As we can see from the previous studies, these pathways were found to be related with OC. We have discovered target genes of three downregulated circulating miRNAs (hsa-miR-1909-5p, hsa-miR-885-5p, and hsa-let-7d-3p) and also related-pathways of these genes. We propose novel mechanisms between miRNAs, target genes and OC that have not been elucidated previously using pathfindR. Our findings can be used as a diagnostic tool in OC. Our perspective on ATXN1L, KAT2B, SRSF2, SRSF6, TRAF6, UBA7, UBE2N, UBE2K, U2SURP, WWP1 and in OC will promote more extensive research on the molecular mechanisms of hsa-miR-1909-5p, hsa-miR-885-5p, and hsa-let-7d-3p and provide a reference for improving the clinical outcome of OC.
The development of prognostic multigene classification protocols can benefit the understanding of tumor biology as well as the prediction of cancer progression and treatment strategies. One critical issue is determining the properly combining the genes [33]. However, studies on the overall survival-related profiles in OC patients have progressed, whereas there have been no large-scale studies based on multicenter validation of gene expression profiles for prediction of disease progression or recurrence in OC patients [34]. Furthermore, pathological staging, which can be determined after surgery and examination of the removed tumor tissue, is likely to be more accurate than clinical staging because it provides direct insight into the extent and nature of the disease [35]. Additionally, differentially expressed miRNAs' target genes, which were identified as multiple candidate genes in OC, were integrated our multigene analysis into the overall survival outcomes and pathological stage analysis using TCGA and GTEx datasets (Table 3). The hsa-miR-885-5p, which was the most gene-targeting score (n = 64) and the most associated miRNA with OC, followed by the hsa-miR-1909-5p and hsa-let-7d-3p miRNAs which were gene-targeting scores (n = 20) and (n = 9), respectively. After FDR analysis, no significant differentially expressed target genes were found to be associated with overall survival. However, we found 24 genes that were targeted by three miRNAs identified in the pathological stage of OC (Fig. 4). Notably, ZNF407 and UBN2 had the lowest p-values for their association with OC pathological stages (p = 0.017856) (Table 3). The molecular role and mechanism of ZNF407 and UBN2 in the development and progression of OC is not well understood. Missense mutations in ZNF407 affect tumor progression in gastrointestinal stromal tumors were reported [36]. Moreover, Tan et al. [37] showed that ZNF407 functiones as a promotive factor in colorectal cancer metastasis and accelerates cell proliferation by regulating phosphoinositide 3-kinase/AKT-mediate pathway. UBN2 is widely expressed in tumor tissues and encodes a nuclear protein that interacts with viral and cellular transcription factors [38]. Zhao et al. [39] suggested that high UBN2 protein expression is an independent prognostic marker to identify patients with poor clinical outcomes in colorectal cancer by affecting the Ras/MAPK pathway. Thus, these genes are likely to be facilitated in therapeutic approaches for OC.
In conclusion, we offer in-silico evidence that validated circulating miRNAs' target genes and enriched pathways are related to OC and have potential roles in theranostics applications. Especially, enrichment pathways and pathological stage-related genes can be combined with validated miRNAs, their multiple analysis can further enhance the molecular etiology of the OC and also can be employed in future research for biomarker and drug development related to OC. Further experimental investigations are required to validate our results which will ultimately provide a new perspective for translational applications in OC management. Our study will allow a greater understanding of broader clinical application prospects.
Notes
Authors’ Contribution
Conceptualization: BD, EG, TG. Data curation: BD, EU. Formal analysis: BD, EU. Funding acquisition: BD, TG. Methodology: BD, EG, EU, OUS, TG. Writing – original draft: BD, EU. Writing – review & editing: BD, EG, EU, OUS, TG.
Conflicts of Interest
No potential conflict of interest relevant to this article was reported.
Acknowledgements
This work was supported by the Istanbul University Scientific Research Projects Department under grant (grant number: 24059). We thank TUBITAK for the scholarship E.U receives (grant number: 118C039).