Next Article in Journal
The Role of Epigenetics in Primary Biliary Cholangitis
Next Article in Special Issue
miR-199a: A Tumor Suppressor with Noncoding RNA Network and Therapeutic Candidate in Lung Cancer
Previous Article in Journal
Epigenetic Mechanisms of Epidermal Differentiation
Previous Article in Special Issue
Exosomal and Non-Exosomal MicroRNAs: New Kids on the Block for Cancer Therapy
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 23
Issue 9
10.3390/ijms23094875
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Novel MicroRNA-Regulated Transcript Networks Are Associated with Chemotherapy Response in Ovarian Cancer

by
Danai G. Topouza
1,
Jihoon Choi
1,
Sean Nesdoly
2,
Anastasiya Tarnouskaya
2,
Christopher J. B. Nicol
1,3,4 and
Qing Ling Duan
1,2,*
1
Department of Biomedical and Molecular Sciences, Queen’s University, 18 Stuart St., Kingston, ON K7L 3N6, Canada
2
School of Computing, Queen’s University, 21-25 Union St., Kingston, ON K7L 2N8, Canada
3
Department of Pathology and Molecular Medicine, Queen’s University, 88 Stuart St., Kingston, ON K7L 3N6, Canada
4
Division of Cancer Biology and Genetics, Queen’s University Cancer Research Institute, Queen’s University, 10 Stuart St., Kingston, ON K7L 3N6, Canada
*
Author to whom correspondence should be addressed.
Int. J. Mol. Sci. 2022, 23(9), 4875; https://0-doi-org.brum.beds.ac.uk/10.3390/ijms23094875
Submission received: 30 March 2022 / Revised: 25 April 2022 / Accepted: 26 April 2022 / Published: 28 April 2022
(This article belongs to the Special Issue MicroRNAs in Cancer Therapy)
Browse Figures
Versions Notes

Abstract

:
High-grade serous ovarian cancer (HGSOC) is a highly lethal gynecologic cancer, in part due to resistance to platinum-based chemotherapy reported among 20% of patients. This study aims to generate novel hypotheses of the biological mechanisms underlying chemotherapy resistance, which remain poorly understood. Differential expression analyses of mRNA- and microRNA-sequencing data from HGSOC patients of The Cancer Genome Atlas identified 21 microRNAs associated with angiogenesis and 196 mRNAs enriched for adaptive immunity and translation. Coexpression network analysis identified three microRNA networks associated with chemotherapy response enriched for lipoprotein transport and oncogenic pathways, as well as two mRNA networks enriched for ubiquitination and lipid metabolism. These network modules were replicated in two independent ovarian cancer cohorts. Moreover, integrative analyses of the mRNA/microRNA sequencing and single-nucleotide polymorphisms (SNPs) revealed potential regulation of significant mRNA transcripts by microRNAs and SNPs (expression quantitative trait loci). Thus, we report novel transcriptional networks and biological pathways associated with resistance to platinum-based chemotherapy in HGSOC patients. These results expand our understanding of the effector networks and regulators of chemotherapy response, which will help to improve the management of ovarian cancer.

1. Introduction

High-grade serous ovarian cancer (HGSOC) is a highly lethal gynecologic cancer, in part due to resistance to first-line, platinum-based chemotherapy treatment among 20% of patients [1]. Chemotherapy-resistant patients have a significantly shorter overall survival (OS) than sensitive patients, and many experience tumor recurrence within six months of completing chemotherapy [2]. There is currently no strategy for predicting response to platinum-based chemotherapy, which reflects our limited understanding of the underlying molecular mechanisms of chemotherapy resistance [3].
MicroRNAs (miRNAs) are small single-stranded noncoding RNAs that post-transcriptionally repress mRNA expression and are involved in the regulation of all biological processes and diseases [4], including ovarian cancer pathogenesis and chemotherapy response [5]. MiRNA and gene expression profiling of patient tumors has the potential to identify signatures that determine disease prognosis [6]. Such prognostic signatures have been implemented in clinic for some cancers [7], but no such tests exist for ovarian cancer patients to date. The detection of transcriptomic signatures of progression-free survival and chemotherapy resistance is an area of active interest, and several previous studies have defined mRNA and miRNA signatures of chemotherapy resistance in HGSOC patients.
The majority of earlier studies reporting genes associated with platinum-based chemotherapy resistance in ovarian cancer patients had employed univariate analysis methods on transcriptomics data [8]. These methods assume that chemotherapy response is driven by a single gene. It is well established, however, that chemotherapy response, like other drug-response outcomes, is a complex multifactorial trait modulated by multiple genes contributing to common biological pathways [9,10]. To date, few studies have investigated chemotherapy response in ovarian cancer using multivariate, machine-learning methods to generate novel hypotheses for the underlying biological mechanisms [11,12,13,14,15]. Fewer still have incorporated multiple types of omics datasets such as mRNA and miRNA sequencing in addition to genomics data to further investigate the regulation of associated gene networks [16,17]. Finally, earlier studies that employed multivariate methods had used expression microarray data, which do not allow for discovery of novel transcript isoforms [18,19]. In contrast, RNA-sequencing data using next-generation technology include gene transcripts that may have been missed by traditional microarray profiling.
In this study, we apply both univariate and multivariate analysis methods to high-throughput miRNA- and mRNA-sequencing data from tumors of HGSOC patients to identify novel biological pathways and networks associated with chemotherapy response. We further determine miRNAs and single-nucleotide polymorphisms (i.e., expression quantitative trait loci, eQTLs) correlated with the expression of the associated transcriptional networks. These findings are validated in two independent ovarian cancer cohorts and improve our understanding of the biological mechanisms underlying resistance to platinum-based chemotherapy in HGSOC patients.

2. Results

2.1. Study Design

We obtained miRNA- and mRNA-sequencing data from HGSOC patients of The Cancer Genome Atlas (TCGA). An analysis pipeline for these data was curated as summarized in Figure 1 and detailed in the Methods and Appendix A. Our pipeline applied select software and tools for quality control of raw sequence reads, alignment to a reference genome, quantification of reads into transcript isoform expression values, outlier removal, and normalization. We separated the TCGA patients into chemotherapy-sensitive and chemotherapy-resistant groups based on their platinum-free interval (Methods). After preprocessing, we performed differential expression and network coexpression analyses, as well as integration of miRNAs, mRNAs, and eQTLs. Our results were validated in two independent HGSOC cohorts.

2.2. Differentially Expressed miRNAs and mRNAs between Sensitive and Resistant Patients

Differential miRNA expression analysis revealed 21 differentially expressed miRNA isoforms between chemotherapy-sensitive and -resistant patients (adjusted p < 0.05), which map to 16 unique miRNAs (Figure 2A, Table S1). Pathway enrichment analysis of these miRNA isoforms revealed 16 pathways, such as blood vessel morphogenesis and negative regulation of autophagy (Table S2).
Differential expression analysis identified 196 mRNAs associated with chemotherapy response (adjusted p < 0.05) that map to 190 unique genes (Figure 2B, Table S3). Pathway enrichment analysis of these associated transcripts indicated enrichment of 41 annotation terms, including B-cell receptor regulation, complement activation, and peptide chain elongation (Table S2).

2.3. miRNA Coexpression Networks Involved in Lipid Transport and Oncogenic Pathways Associated with Chemotherapy Response

Weighted gene coexpression network analysis (WGCNA) of the miRNA dataset constructed 100 coexpression modules (Table S4), of which three are associated with chemotherapy response. The ivory (p = 0.0098, log OR = −5.67) and lightcoral (p = 0.042, log OR = −4.36) modules are negatively associated with chemotherapy resistance, while the third plum network (p = 0.045, log OR = 4.29) is positively associated with chemotherapy resistance. The ivory network consisted of 25 miRNA isoforms mapping to 11 unique miRNAs (Figure 3A, Table S5), which are enriched for seven pathways and functions, including regulation of lipoprotein transport and cholesterol efflux (Table S2). The lightcoral module consists of 17 isoforms of miR-187 and the plum network consists of 17 isoforms of miR-221 and miR-222 (Figure 3A, Table S5). While no pathway annotations were derived for the lighcoral module, the plum module is enriched for 18 pathways and oncogenic functions, such as inhibition of the TRAIL-activated apoptotic pathway and inflammatory cytokine production, and upregulation of protein kinase B signaling (Table S2).

2.4. mRNA Coexpression Networks Involved in Protein Ubiquitination and Fatty-Acid Metabolism Associated with Chemotherapy Response

WGCNA of the mRNA transcripts resulted in 58 coexpression modules, of which 2 were associated with platinum-based chemotherapy (Table S6). First, the lavenderblush3 module is negatively associated with chemotherapy resistance (p = 0.016, log OR = −5.40). This module contains 39 transcripts, mapping to 31 unique genes (Figure 3B, Table S7). Pathway analysis indicates enrichment of biological pathways related to protein ubiquitination, and the binding motif for the transcription factor GABP-alpha (Table S2). Second, the darkolivegreen module is significantly upregulated in chemoresistant patients (p = 0.032, log OR = 13.63) and contains 82 transcripts mapping to 80 unique genes (Figure 3B, Table S7). This module is significantly enriched for the protein-containing complex term, as well as for eight pathways involved in fatty-acid metabolism with nominal significance (Table S2).

2.5. Germline eQTLs May Regulate the Expression of Associated miRNAs, mRNAs, and Networks

Integrative analysis with germline SNP data identified 268 unique cis-eQTLs associated with the expression of significant mRNAs and miRNAs (Table S8). A total of 20 SNPs are associated with the expression of 7 significant miRNAs, and 248 SNPs are associated with the expression of 55 significant mRNAs. Of the 268 SNPs, 118 are novel eQTLs, whereas 126 are previously known and 24 are not yet recorded in the annotation database. The majority (227) are predicted to alter regulatory motifs, and 67 are associated with 94 human phenotypes from published genome-wide association studies. The most common phenotypes are related to triglycerides, high-density lipoprotein (HDL), and low-density lipoprotein (LDL) cholesterol.

2.6. Network Integration Reveals miRNA-Mediated Regulation of Chemotherapy Response Mechanisms

Integration of the associated miRNA and mRNA networks determined that the plum miRNA network significantly correlates with the lavenderblush3 mRNA network (Spearman correlation of module eigengenes, ρ = −0.26, p < 0.001), and the ivory miRNA network significantly correlates with the darkolivegreen mRNA network (Spearman correlation of module eigengenes, ρ = −0.17, p = 0.023). Annotations using miRNet and miRGate determined that 20 of these mRNA-miRNA interactions are experimentally validated, while 15 others are supported by in silico predictions (Table 1). The experiments that validated miRNA binding to mRNA molecules involved HITS-CLIP (high-throughput sequencing of RNA isolated by crosslinking immunoprecipitation), AGO-CLIP (Argonaute-crosslinking and immunoprecipitation), CLASH (crosslinking, ligation, and sequencing of hybrids), and PAR-CLIP (photoactivatable ribonucleoside-enhanced crosslinking and immunoprecipitation) assays (Table S9). The computational prediction algorithms inferred miRNA-mRNA binding by assessing the complementarity between the miRNA seed sequence and the mRNA transcript, as well as the mRNA-miRNA duplex energy (Table S9). Combined with the potential cis-eQTL regulation of mRNAs and miRNAs in these networks, these results reveal an integrative, multi-omics view of transcriptional networks associated with chemotherapy response in ovarian cancer (Figure 4).

2.7. Replication in Two Independent Ovarian Cancer Cohorts

Replication of the differentially expressed miRNAs used miRNA-seq data from the MITO cohort (Table S10). The ivory and plum miRNA network modules replicated in the MITO cohort (p = 6.1 × 10−4, log HR = −0.78 and p = 0.022, log HR = −0.52, respectively), while the lightcoral module reached nominal significance (p = 0.057, log HR = −0.43) (Figure 5A–C).
Replication analysis of the differentially expressed mRNAs used RNA microarray data from the AOCS cohort (Table S11). The lavenderblush3 and darkolivegreen mRNA network modules replicated in the AOCS cohort (p = 1.6 × 10−10, log HR = −1.27 and p = 1.3 × 10−10, log HR = −1.27, respectively) (Figure 5D,E).

3. Discussion

In this study, we identified miRNA and mRNA transcripts and networks associated with chemotherapy response in HGSOC patients. Our findings implicate novel and known biological pathways that were further replicated in independent cancer cohorts. In addition, we identified potential interactions among these miRNA and mRNA networks, as well as eQTLs that potentially regulate the associated transcripts. Thus, our results provide an integrative, multi-omics view of biological networks associated with chemotherapy response in HGSOC.
We identified one miRNA coexpression network (ivory), associated with chemotherapy sensitivity in HGSOC, which is involved in the negative regulation of lipid transport. This enrichment is mainly mediated by miR-128-1 and miR-128-2, which play a key role in cholesterol and lipid homeostasis through their suppression of the ABCA1 cholesterol efflux transporter and the low-density lipoprotein receptor (LDLR) [43,44]. MiR-148a, which is significantly downregulated in resistant patients, is also a regulator of these key genes [43]. The overexpression of ABCA1 is associated with reduced survival in OC patients [45], and levels of LDLR are increased in chemoresistant OC cell lines [46]. In addition, overexpression of miR-128 promotes sensitivity to cisplatin in previously resistant OC cells [47]. Our results are consistent with the chemosensitivity-promoting role of miR-128 and its potential activity in cholesterol efflux inhibition alongside miR-148a in this cohort.
We identified a second miRNA coexpression network (plum) consisting of miR-221 and miR-222 isoforms, which have been implicated in the development of chemotherapy resistance in OC. Expression of miR-221/miR-222 transcripts is high in cisplatin-resistant OC cell lines, and their inhibition increases cellular sensitivity [48]. Overexpression of miR-221 and miR-222 promotes proliferation of OC cell lines [49,50] and is associated with reduced disease-free and overall survival [49]. Thus, our findings are consistent with earlier studies showing increased activity of miR-221 and miR-222 in chemoresistant tumors.
The lavenderblush mRNA coexpression network was significantly upregulated in platinum-sensitive patients, which replicated in another independent cohort (AOCS). This module consists of genes involved in ubiquitin-mediated proteolysis in the endoplasmic reticulum (ER). We also detected a significant downregulation of genes responsible for translation initiation in sensitive patients in our differential mRNA expression analysis. These findings suggest that the unfolded protein response (UPR), a cellular process responsible for resolving ER stress, may be increasingly activated in sensitive patients compared to resistant cases. The UPR alleviates ER stress through several pathways, including increased ER-associated protein degradation (ERAD) to remove misfolded proteins, and inhibition of translation to reduce protein load in the ER [51]. ER stress promotes cisplatin resistance in OC cell lines [52] and the upregulation of ERAD genes such as VCP in the lavenderblush3 module is associated with longer OS and platinum sensitivity in HGSOC cohorts [17,53]. Finally, the lavenderblush3 genes VCP, DNAJA1, and TOPORS are overexpressed in platinum-sensitive HGSOC patients as part of a cell cycle and damage-response-associated network [16].
The second mRNA coexpression network associated with chemotherapy resistance in our HGSOC cohort that replicated in the AOCS (darkolivegreen) included genes associated with fatty-acid metabolism (SREBF1, ACAA1, ACADVL), and the protein kinase B oncogene (AKT1), which promotes de novo lipid biosynthesis in cancer [54]. SREBF1 is a key enzyme for cholesterol and fatty-acid synthesis, and an essential gene for OC tumor growth [55]. Specifically, SREBF1 is activated by AKT1, promoting fatty-acid synthesis [56], which favors cell proliferation in OC [57]. Expression of ACADVL, involved in the β-oxidation of long-chain fatty acids, is linked to OC metastasis and cell survival [58]. Our findings indicate the upregulation of these lipid metabolism genes among chemotherapy-resistant patients. Lipid metabolism dysregulation activates the UPR, which triggers lipid metabolism-based adaptations in the cell through several pathways, including SREBF1 regulation [59]. The interaction of these pathways may present a link between our two gene coexpression modules and warrants further study.
Differential expression analysis identified a downregulation of mRNA transcripts involved in the adaptive immune system, which is associated with chemoresistance. Previous studies reported that a high tumor immune score is a strong predictor of chemosensitivity in HGSOC [60]. In addition, there are potential links between this immune response activation, UPR, and lipid metabolism. ER stress can induce proinflammatory cytokine production and UPR activation in tumor cells [61], which can disrupt dendritic cell function in the OC tumor microenvironment [62]. Moreover, dendritic cell function can also be inhibited by increased lipid uptake in various cancers [63].
Integrative analysis of mRNA-seq and miRNA-seq datasets identified potential interactions of the associated transcript coexpression modules. The overexpression of miR-221/222 in resistant patients may be inhibiting the chemosensitivity-associated lavenderblush3 mRNA network, revealing a novel potential mechanism of chemotherapy resistance. This finding, combined with the accumulating evidence of miR-221/miR-222 involvement in chemoresistance, may point to a promising avenue for therapeutic intervention. However, overexpression of miR-221/miR-222 promotes UPR-induced apoptosis in hepatocellular carcinoma (HCC) cells [64]. Additionally, ER stress suppresses miR-221/miR-222 in HCC, promoting resistance to apoptosis. The contribution of this mechanism to chemotherapy response in HGSOC is currently unclear and presents an area for future investigation.
We also identified potential regulation of the darkolivegreen mRNA module by the ivory miRNA network, which may inhibit lipid metabolism in chemotherapy-sensitive patients. As increased lipid metabolism by cancer cells is a known mechanism of chemoresistance in HGSOC, this miRNA-mediated inhibition may present a novel mechanism of chemotherapy sensitivity.
Finally, cis-eQTL analysis identified known and novel genomic variants correlated with the expression of mRNAs and miRNAs, which are associated with lipid-related phenotypes. These SNPs have not been previously associated with platinum-based chemotherapy response in ovarian cancer [65] and may represent novel genomic associations with chemotherapy response. High HDL and triglyceride levels have been correlated with increased cancer stage at diagnosis in OC patients [66]. In addition, advanced-stage OC patients with high LDL levels have a shorter PFS than patients with normal levels [67]. Further investigation of these eQTLs is necessary to further elucidate their role in platinum-based chemotherapy resistance and HGSOC prognosis.
To date, there have been few studies to investigate miRNA and mRNA network associations with chemotherapy response in ovarian cancer using multivariate analysis methods. Bernardini et al. (2005) used unsupervised two-dimensional hierarchical clustering and feature selection to identify genes that were predictive of response to platinum chemotherapy [12]. Spentzos et al. (2005) used pattern recognition and compound covariate predictive algorithms that identified a multigene pattern to classify patient chemotherapy response [11]. Bagnoli et al. (2016) applied a semisupervised principal component analysis method on selected miRNAs, leading to a prognostic miRNA model whose expression is associated with risk of disease progression [15]. Chen et al. (2018) used WGCNA to identify gene networks associated with chemosensitivity [13]. Zhang et al. (2019) also used WGCNA and identified gene networks associated with chemoresistance [14]. The above studies all performed multivariate analysis using a single data modality; two additional studies made use of multiple omics datasets to investigate chemotherapy response. Benvenuto et al. (2020) used a micrographite algorithm to integrate significant mRNA and miRNA expression profiles into a single network that distinguished chemotherapy-sensitive and -resistant patients [16]. Finally, our group previously used WGCNA to identify gene networks associated with chemotherapy sensitivity; the expression of significant genes was integrated with patient germline genomic polymorphisms to identify cis-eQTLs that may regulate the expression of these networks [17]. All of the above studies used expression data from microarrays, which do not allow for the distinction of miRNAs from the same sequence family and of mRNA transcript isoforms [18,19]. Our use of miRNA- and mRNA-sequencing data, as well as our analysis of three different data modalities, results in a more detailed profiling of the transcriptome and generates novel hypotheses for the biological mechanisms underlying variable chemotherapy response in ovarian cancer.
Our study provides novel insight of the underlying mechanisms modulating resistance to platinum-based chemotherapy in HGSOC. Specifically, we conducted whole-transcriptome analysis of miRNA-seq and mRNA-seq data to generate novel mechanistic hypotheses using both univariate and network methods. Moreover, we integrated this data with miRNA-seq and genome-wide SNPs to determine potential regulation of the associated transcripts and networks. Our findings implicate novel and known signaling pathways and networks associated with chemotherapy response in HGSOC as well as regulators, which could become novel drug targets. Further studies are needed to validate these findings in other cancers, and to investigate the contribution of these networks to patients’ overall survival.

4. Materials and Methods

4.1. Chemotherapy Response Classification

Sequencing of miRNA and mRNA was derived from chemotherapy-naïve tumors of 191 and 205 HGSOC patients of TCGA, respectively [68]. Patients who received platinum-based adjuvant chemotherapy were selected and classified for chemotherapy response based on their platinum-free interval. Sensitive patients remained cancer-free for at least 12 months after chemotherapy completion, whereas resistant patients experienced cancer recurrence within 6 months (Table 2; Appendix A).

4.2. Processing of Sequencing Data

The miRNA-sequencing data were downloaded as quantified expression files (level 3 data from TCGA). Sequencing reads from mRNA were downloaded as FASTQ files (level 1 data from TCGA), filtered for base-quality, aligned, and quantified (detailed in Appendix A). Both mRNA and miRNA datasets underwent outlier detection, normalization, and nonspecific filtering, resulting in 49,116 mRNA and 4479 miRNA transcripts for further analyses.

4.3. Differential Expression Analysis

Differentially expressed miRNA and mRNA transcripts were detected between sensitive and resistant patients using a negative binomial generalized linear model (GLM) in the DESeq2 R package [27]. DESeq2 uses the median-of-ratios normalization method, which divides transcript counts in each sample by a size factor determined by the ratios of gene counts in the sample to the average gene counts among all samples [69]. This method considers the sequencing depth and RNA composition in each sample and is a recommended normalization method for RNA-sequencing data [70,71]. This analysis controlled for patients’ ages at diagnosis, as resistant patients were significantly older (Table 2). The Benjamini–Hochberg method corrected for multiple testing.

4.4. Weighted Network Correlation Analysis

The weighted gene coexpression network analysis (WGCNA) R package [29] was used to identify modules of coexpressed miRNA and mRNA transcripts using an unsupervised machine-learning approach. We performed multivariate WGCNA to evaluate the association of miRNA and mRNA networks with chemotherapy response. This method groups coexpressed transcripts into modules prior to testing for association to the clinical outcome. This analysis identifies groups of transcripts that individually may have modest effects on chemotherapy response, but collectively contribute to a common biological network or pathway that is significantly associated with the outcome. In addition, this method reduces the sequencing datasets into a smaller number of transcript modules and uses Principal Component Analysis (PCA) to further summarize the information of each cluster into a representative value, referred to as the eigengene, for association testing. This reduces the multiple testing corrections needed. Module eigengenes were used to determine association with chemotherapy response using a GLM, adjusted for patients’ age as a covariate (Appendix A).

4.5. Pathway Enrichment Analysis

Pathway enrichment analysis with g:Profiler [72] was used to determine overrepresentation of biological pathways from lists of differentially expressed transcripts (miRNA and mRNAs) and coexpression networks (Appendix A).

4.6. Expression Quantitative Trait Locus Analysis

Germline single-nucleotide polymorphisms (SNPs) from TCGA–HGSOC patients were imputed as described by Choi et al. [17] before undergoing quality control and linkage disequilibrium-based pruning, retaining 1,722,608 common SNPs for analysis. SNPs were integrated with patient miRNA-seq data (n = 178) and mRNA-seq data (n = 167) to identify correlations with transcript expression (eQTLs) using the MatrixEQTL R package [33] (Appendix A).

4.7. mRNA-microRNA Interaction Analysis

Potential regulation of mRNA networks by miRNAs was tested on a subset of patients with both mRNA and miRNA data from the same tumor (n = 165). We measured the Spearman correlation of module eigengenes from the mRNA and miRNA coexpression networks, as well as the Spearman correlation of individual mRNA and miRNA transcript expression. Results were validated using miRNet [30], a database of experimentally validated mRNA-miRNA interactions, and miRGate [31], a tool that identifies predicted mRNA-miRNA interactions based on the consensus of several algorithms that assess the complementarity between miRNA seed sequences and mRNA transcript sequences, as well as the mRNA-miRNA duplex energies.

4.8. Replication Cohorts and Analysis

Our results were replicated using two independent ovarian cancer cohorts. First, miRNA results were replicated in the Multicenter Italian Trial in Ovarian cancer cohort (MITO; GSE25204; n = 130) [73]. Next, mRNA results were replicated in the Australian Ovarian Cancer Study cohort (AOCS; GSE9891; n = 285) [74]. Replication of differentially expressed transcripts used the auto-cutoff Kaplan–Meier analysis method from the KM Plotter tool [75] to test the association of each transcript with progression-free survival (PFS) in the AOCS and MITO cohorts. Validation of transcript networks used the Prognostic Index estimation method from the SurvExpress tool [32] to test the association of miRNA and mRNA networks with PFS in the above cohorts (Appendix A).

4.9. Software and Statistical Analysis

All statistical analyses were performed using R (v. 3.6.0) [76] in the RStudio environment (v. 1.1.383) [77]. The association tests of continuous patient clinical data, such as age, with chemotherapy-response-employed t-tests, while categorical patient clinical data, such as cancer stage and tumor subtype, were tested with Fisher’s exact tests and Chi-squared tests (stats R package, v. 3.6.0) [76]. Differential expression analyses were performed using DESeq2 (v. 1.26.0) [27] and plotted with ggplot2 (v. 3.3.5) [78] and ggrepel (v. 0.9.1) [79]. Coexpressed transcript networks were detected using WGCNA (v. 1.69) [29] and plotted with Cytoscape (v. 3.7.0) [80]. The association of transcript networks with chemotherapy response was tested using a generalized linear model (stats R package, v. 3.6.0). The correlation of miRNA and mRNA transcript expression was tested using Spearman’s correlation (stats R package, v. 3.6.0). Detection of eQTLs was performed using an additive linear model in the MatrixEQTL R package (v. 2.3) [33]. Validation of differentially expressed transcripts in independent cohorts was performed using Kaplan–Meier analysis with the survival R package (v. 3.2-7) [81] and plotted with survplots from the R package rms (v. 6.1-1) [82]. Validation of network analysis results in independent cohorts was performed using Cox proportional hazards regression models to generate a prognostic index estimator for each network (survival R package, v. 3.2-7).

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms23094875/s1.

Author Contributions

Conceptualization, Q.L.D.; methodology, Q.L.D.; software, D.G.T., J.C. and S.N.; validation, D.G.T.; formal analysis, D.G.T. and J.C.; investigation, D.G.T.; resources, Q.L.D.; data curation, D.G.T., J.C., S.N. and A.T.; writing—original draft preparation, D.G.T.; writing—review and editing, D.G.T., J.C., C.J.B.N. and Q.L.D.; visualization, D.G.T.; supervision, Q.L.D.; project administration, Q.L.D.; funding acquisition, Q.L.D. All authors have read and agreed to the published version of the manuscript.

Funding

D.G.T. is funded by internal awards in the Department of Biomedical and Molecular Sciences at Queen’s University. Q.L.D. receives funding from the Queen’s University National Scholar Award and the Canadian Institutes of Health Research (PJT-178390).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Informed consent was obtained from all subjects involved in the Cancer Genome Atlas (TCGA) research project [68]. Data access and research methods were carried out in accordance with relevant guidelines and regulations set by the National Cancer Institute (NCI) and the National Human Genome Research Institute (NHGRI), in compliance with the HIPAA Privacy Rule, and all protocols were approved by the Human Research Ethics Boards at Queen’s University. The terms of data access are outlined in the Data Use Certification Agreement: https://dbgap.ncbi.nlm.nih.gov/aa/wga.cgi?view_pdf&wlid=10654&tlsid=274 (accessed 22 April 2022). Details of the human subject protection and data access policies implemented by TCGA are as described: https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga/history/policies/tcga-human-subjects-data-policies.pdf (accessed 22 April 2022).

Data Availability Statement

Transcriptomics, genomics, and clinical data from the TCGA cohort supporting the conclusions of this article are available in the Genomic Data Commons (GDC) Data Portal at https://portal.gdc.cancer.gov/, project ID TCGA-OV. Clinical and gene expression data of the AOCS cohort supporting the conclusions of this article are available in the Gene Expression Omnibus (GEO) database at https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/geo/, reference number GSE9899. Clinical and miRNA expression data of the MITO cohort supporting the conclusions of this article are available in the GEO database at https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/geo/, reference number GSE25204.

Acknowledgments

Computations in this manuscript were performed on resources and with support provided by the Centre for Advanced Computing (CAC) at Queen’s University in Kingston, Ontario. The CAC is funded by the Canada Foundation for Innovation, the Government of Ontario, and Queen’s University.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A. Supplementary Methods

Appendix A.1. Chemotherapy Response Classification

Sensitive patients remained cancer-free one year after completion of chemotherapy, whereas resistant patients experienced cancer recurrence within 6 months of chemotherapy completion. We excluded patients who developed new tumors between six months to one year from chemotherapy completion to enrich for genetic differences between the drug response groups. Included patients received platinum-based adjuvant chemotherapy, which consisted of a platinum agent (3.66–3.90% of patients) or a combination of a platinum agent and a taxane (96.10–96.34% of patients) (Table 2). We excluded patients who were treated with an alternative therapy instead of platinum-based adjuvant chemotherapy, patients without a tumor recurrence who passed away within six months of chemotherapy completion, patients who were followed for less than one year after chemotherapy completion, and patients without transcriptomic data.

Appendix A.2. RNA-Sequencing Data Preprocessing

Raw mRNA-sequencing reads from frozen chemotherapy-naïve tumor samples of high-grade serous ovarian cancer (HGSOC) patients were obtained from the TCGA database [68] on 31 August 2017 using the TCGAbiolinks R package [83]. Patient mRNA-seq data were filtered for base quality using the Trimmomatic software package [22]. Bases were evaluated in a window of four bases progressing from the 5′ end of each read to the 3′ end. Each read was trimmed when the average quality within the window dropped below 15 on the PHRED base 33 quality scale, which indicates a 3.1% probability of error [84]. Sequences matching Illumina HiSeq adapters and Illumina PCR primers were removed.
Patient RNA-sequence reads were aligned to the hg19 reference genome [85] using the HISAT2 software [24], with a mean overall alignment rate of 96.38%. The aligned reads were assembled into transcripts by the StringTie software package [25] and quantified using the GRCh37.87 Ensembl transcript list (release 92) [86], resulting in the quantification of 196,464 gene transcripts and isoforms for each patient.
Aligned and quantified frozen chemotherapy-naïve tumor microRNA-sequencing data (Level 3) for the TCGA HGSOC cohort were obtained from the Broad Institute Genome Data Analysis Center Firehose [26] (http://gdac.broadinstitute.org) on 17 September 2018 (TCGA data version 2016_01_28 for OV). Patient data were acquired in the form of quantified microRNA-seq isoform counts (n = 29,382) in the .txt format. The microRNA (miRNA) isoform data had been aligned to GRCh36/hg18 and quantified using TCGA’s modified version of the British Columbia Genome Sciences Centre miRNA profiling pipeline [87].

Appendix A.3. Patient-Level Quality Control

Patient mRNA data were iteratively filtered for outliers based on Spearman’s correlation of gene expression between patients, using Tukey’s outlier labeling threshold (Q1—α * IQR) [88] and fine-tuning the α value based on our sample size (α = 2.3) [89] as described by Panarelli et al. [90]. This process removed four patients. Patient miRNA expression profiles were also filtered for outliers using the above method, with an α value of 2.4 based on the larger number of patients in the miRNA cohort. This process removed nine patients.

Appendix A.4. Transcript-Level Quality Control

Transcript count data were normalized using the median of ratios method from the DESeq2 R software package [27]. Next, we applied nonspecific filtering to the mRNA expression dataset using median absolute deviation (MAD) [28], and retained the top 25% of mRNA transcripts with the highest variance among patients for further analysis (n = 49,116). As for nonspecific filtering of the miRNA dataset, we omitted isoforms with variance below the median value, resulting in 14,691 remaining miRNA transcripts. In addition, we removed miRNA isoforms with low expression among patients (isoform count mean < 1), resulting in 4479 miRNAs for further analysis.

Appendix A.5. Differential Expression Analysis

Volcano plot figures were generated using the R packages ggplot2 (v. 3.3.5) [78] and ggrepel (v 0.9.1) [79].

Appendix A.6. Pathway Enrichment Analysis

Three differentially expressed miRNAs (miR-151a, miR-103a-1, miR-203a) were re-annotated with the newest miRbase [22] identifiers before pathway analysis with g:Profiler [72]. A total of 3 of the 16 miRNAs (miR-320a, miR-1974, and miR-886) are not yet identified in the g:Profiler database and not included in the pathway enrichment analysis.

Appendix A.7. Weighted Correlation Network Analysis

Patient transcript expression data were analyzed using the weighted correlation network analysis (WGCNA) R package [29]. We used this tool to calculate the Pearson correlation between all transcripts, then used an adjacency function to form a network where each node is a gene, and each connecting edge corresponds to the strength of the correlation between nodes. The power adjacency function raised the correlations to a power β, making the connectivity distribution of our gene expression network approach that of a scale-free network. Fine tuning indicated β = 10 to be optimal for the mRNA dataset, and β = 9 to be optimal for the miRNA dataset. We constructed a Topological Overlap Matrix (TOM), where the topological overlap of two transcripts is a function of their adjacency and number of shared connections. The TOM-based dissimilarity was then used in hierarchical clustering to generate gene modules. We used a minimum module size of 30 to encourage larger gene clusters, as recommended by the WGCNA manual [91]. After clustering, we merged modules that were not sufficiently distinct from each other to produce robust networks that better capture biological pathways. We calculated a module membership value for each transcript that shows how well a transcript fits into each of the available modules. A total of 50% of transcripts from each module were tested for membership. If more than 25% of transcripts tested had a higher membership value for a module other than their own, those two modules were merged. This process was performed in rounds until no more modules could be merged. Finally, principal component analysis was conducted for all transcripts within a module, resulting in a value called an eigengene. We used these module eigengene values and the patient age covariate to construct a generalized linear model that revealed modules with significant correlation with chemotherapy response. Networks with a significant association were functionally annotated using pathway analysis. Cytoscape was used for network visualization [80].

Appendix A.8. Expression Quantitative Trait Locus (eQTL) Analysis

Germline single-nucleotide-polymorphism (SNP) data were profiled in normal tissues of HGSOC patients in the TCGA cohort using the Affymetrix SNP Array 6.0. A total of 47,960,330 genetic polymorphisms for 262 patients were obtained after phasing and imputing as described by Choi et al. [17]. The imputed data were then processed for quality control using PLINK 1.9 [92]. We first performed patient-level quality control: 20 patients were removed due to low (F < −0.05) or excessive (F > 0.05) heterozygosity, and 2 patients were removed due to high genetic relatedness (pi-hat > 0.9), leaving 240 patients for further analysis. We then performed variant-level quality control, starting with linkage disequilibrium (LD)-based variant pruning. This process evaluated variants in a window of 50 SNPs, which shifted by 5 SNPs after every iteration, and removed variants with a variance inflation factor (VIF) larger than 2 (r2 > 0.5). A total of 7,075,175 independent SNPs were retained after LD pruning. After checking for allelic independence with the Hardy–Weinberg equilibrium, 112 more variants were removed. Finally, we retained SNPs with a minor allele frequency ≥ 5% and variant missingness < 10%, resulting in 1,722,608 variants to be used for further analysis. A total of 167 patients from the RNA-Seq cohort and 178 patients from the miRNA-Seq cohort had genomic data available following quality control.
Next, we tested the association of common patient polymorphisms (n = 1,722,608) with mRNA transcript expression (n = 196,464) or miRNA isoform expression (n = 29,382), in order to determine novel expression quantitative trait loci (eQTLs) in our dataset. We used the MatrixEQTL R package [33] to test SNPs located within 1 Mb of recorded transcripts [93] (cis-eQTLs) for association with transcript expression. Based on this distance, 1,676,365 common SNPs were tested for association with mRNA expression, while 356,189 common SNPs were tested for association with miRNA expression. miRNA isoform locations were converted to hg19 coordinates using UCSC’s liftOver utility [94].
We identified 121,900 significantly associated SNP-mRNA transcript pairs and 36,995 significantly associated SNP-miRNA isoform pairs after false-discovery-rate correction. Of those, 248 SNPs were associated with 55 mRNA transcripts found to be significant in differential expression or network analysis, while 20 SNPs were associated with 7 miRNAs that were significant in differential expression or network analysis (Table S8). We mapped these SNPs to an rsID using the dbSNP Build150 Human Variation Set [95], and functionally annotated them using HaploReg v4.1 [96].

Appendix A.9. Replication Analysis

The replication of our findings was performed in two independent ovarian cancer cohorts The replication of mRNA results was performed on tumor transcriptomic data collected on the Affymetrix Human Genome U133 Plus 2.0 expression microarray from a cohort of 285 patients from the Australian Ovarian Cancer Study (AOCS; GSE9891) [74]. The replication of miRNA results was performed on normalized tumor miRNA data collected on the Illumina Human v2 MicroRNA expression beadchip microarray from a cohort of 130 patients from the Multicenter Italian Trial in Ovarian cancer cohort (MITO; GSE25204) [73]. The AOCS expression microarray data were retrieved from GEO and processed as described by Choi et al. [17] before analysis. We retained only patients with the serous ovarian cancer subtype that received platinum-based therapy from AOCS (n = 226). Outliers were identified using the arrayQualityMetrics R package [97], leaving 224 patients from the AOCS cohort and 107 patients from the MITO cohort for further analysis.
Differential mRNA and miRNA expression analysis results were replicated using the expression auto-cutoff and Kaplan–Meier analysis method used by the KM-Plotter tool [75]. The replication of our network analysis results was performed using the Prognostic Index (PI) estimation method as described by Aguirre-Gamboa et al. in the SurvExpress biomarker validation tool [32]. We used the AOCS and MITO datasets to validate our mRNA and miRNA networks, respectively. Survival plots in Figure 5 were generated using the R package survminer (v. 0.4.8) [98].

References

  1. Dion, L.; Carton, I.; Jaillard, S.; Nyangoh Timoh, K.; Henno, S.; Sardain, H.; Foucher, F.; Leveque, J.; de la Motte Rouge, T.; Brousse, S.; et al. The Landscape and Therapeutic Implications of Molecular Profiles in Epithelial Ovarian Cancer. J. Clin. Med. 2020, 9, 2239. [Google Scholar] [CrossRef] [PubMed]
  2. Davis, A.; Tinker, A.V.; Friedlander, M. “Platinum resistant” ovarian cancer: What is it, who to treat and how to measure benefit? Gynecol. Oncol. 2014, 133, 624–631. [Google Scholar] [CrossRef] [PubMed]
  3. Freimund, A.E.; Beach, J.A.; Christie, E.L.; Bowtell, D.D.L. Mechanisms of Drug Resistance in High-Grade Serous Ovarian Cancer. Hematol. Oncol. Clin. N. Am. 2018, 32, 983–996. [Google Scholar] [CrossRef] [PubMed]
  4. Bartel, D.P. Metazoan MicroRNAs. Cell 2018, 173, 20–51. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Alshamrani, A.A. Roles of microRNAs in Ovarian Cancer Tumorigenesis: Two Decades Later, What Have We Learned? Front. Oncol. 2020, 10, 1084. [Google Scholar] [CrossRef]
  6. Qian, Y.; Daza, J.; Itzel, T.; Betge, J.; Zhan, T.; Marme, F.; Teufel, A. Prognostic Cancer Gene Expression Signatures: Current Status and Challenges. Cells 2021, 10, 648. [Google Scholar] [CrossRef]
  7. McVeigh, T.P.; Kerin, M.J. Clinical use of the Oncotype DX genomic test to guide treatment decisions for patients with invasive breast cancer. Breast Cancer Targets Ther. 2017, 9, 393–400. [Google Scholar] [CrossRef] [Green Version]
  8. Lloyd, K.L.; Cree, I.A.; Savage, R.S. Prediction of resistance to chemotherapy in ovarian cancer: A systematic review. BMC Cancer 2015, 15, 117. [Google Scholar] [CrossRef] [Green Version]
  9. Kalow, W. Pharmacogenetics and pharmacogenomics: Origin, status, and the hope for personalized medicine. Pharm. J. 2006, 6, 162–165. [Google Scholar] [CrossRef] [Green Version]
  10. Binju, M.; Padilla, M.A.; Singomat, T.; Kaur, P.; Suryo Rahmanto, Y.; Cohen, P.A.; Yu, Y. Mechanisms underlying acquired platinum resistance in high grade serous ovarian cancer—A mini review. Biochim. Biophys. Acta Gen. Subj. 2019, 1863, 371–378. [Google Scholar] [CrossRef]
  11. Spentzos, D.; Levine, D.A.; Kolia, S.; Otu, H.; Boyd, J.; Libermann, T.A.; Cannistra, S.A. Unique gene expression profile based on pathologic response in epithelial ovarian cancer. J. Clin. Oncol. 2005, 23, 7911–7918. [Google Scholar] [CrossRef] [PubMed]
  12. Bernardini, M.; Lee, C.H.; Beheshti, B.; Prasad, M.; Albert, M.; Marrano, P.; Begley, H.; Shaw, P.; Covens, A.; Murphy, J.; et al. High-resolution mapping of genomic imbalance and identification of gene expression profiles associated with differential chemotherapy response in serous epithelial ovarian cancer. Neoplasia 2005, 7, 603–613. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Chen, Y.; Bi, F.; An, Y.; Yang, Q. Coexpression network analysis identified Kruppel-like factor 6 (KLF6) association with chemosensitivity in ovarian cancer. J. Cell. Biochem. 2018, 120, 2607–2615. [Google Scholar] [CrossRef] [PubMed]
  14. Zhang, L.; Zhang, X.; Fan, S.; Zhang, Z. Identification of modules and hub genes associated with platinum-based chemotherapy resistance and treatment response in ovarian cancer by weighted gene co-expression network analysis. Medicine 2019, 98, e17803. [Google Scholar] [CrossRef] [PubMed]
  15. Bagnoli, M.; Canevari, S.; Califano, D.; Losito, S.; Maio, M.D.; Raspagliesi, F.; Carcangiu, M.L.; Toffoli, G.; Cecchin, E.; Sorio, R.; et al. Development and validation of a microRNA-based signature (MiROvaR) to predict early relapse or progression of epithelial ovarian cancer: A cohort study. Lancet Oncol. 2016, 17, 1137–1146. [Google Scholar] [CrossRef]
  16. Benvenuto, G.; Todeschini, P.; Paracchini, L.; Calura, E.; Fruscio, R.; Romani, C.; Beltrame, L.; Martini, P.; Ravaggi, A.; Ceppi, L.; et al. Expression profiles of PRKG1, SDF2L1 and PPP1R12A are predictive and prognostic factors for therapy response and survival in high-grade serous ovarian cancer. Int. J. Cancer 2020, 147, 565–574. [Google Scholar] [CrossRef] [PubMed]
  17. Choi, J.; Topouza, D.G.; Tarnouskaya, A.; Nesdoly, S.; Koti, M.; Duan, Q.L. Gene networks and expression quantitative trait loci associated with adjuvant chemotherapy response in high-grade serous ovarian cancer. BMC Cancer 2020, 20, 413. [Google Scholar] [CrossRef]
  18. Zhang, W.; Yu, Y.; Hertwig, F.; Thierry-Mieg, J.; Zhang, W.; Thierry-Mieg, D.; Wang, J.; Furlanello, C.; Devanarayan, V.; Cheng, J.; et al. Comparison of RNA-seq and microarray-based models for clinical endpoint prediction. Genome Biol. 2015, 16, 133. [Google Scholar] [CrossRef] [Green Version]
  19. Hafner, M.; Renwick, N.; Farazi, T.A.; Mihailovic, A.; Pena, J.T.; Tuschl, T. Barcoded cDNA library preparation for small RNA profiling by next-generation sequencing. Methods 2012, 58, 164–170. [Google Scholar] [CrossRef] [Green Version]
  20. Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data 2010. Available online: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed on 22 April 2022).
  21. Ewels, P.; Magnusson, M.; Lundin, S.; Kaller, M. MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics 2016, 32, 3047–3048. [Google Scholar] [CrossRef] [Green Version]
  22. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef] [Green Version]
  23. Haeussler, M.; Zweig, A.S.; Tyner, C.; Speir, M.L.; Rosenbloom, K.R.; Raney, B.J.; Lee, C.M.; Lee, B.T.; Hinrichs, A.S.; Gonzalez, J.N.; et al. The UCSC Genome Browser database: 2019 update. Nucleic Acids Res. 2019, 47, D853–D858. [Google Scholar] [CrossRef] [Green Version]
  24. Kim, D.; Paggi, J.M.; Park, C.; Bennett, C.; Salzberg, S.L. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 2019, 37, 907–915. [Google Scholar] [CrossRef]
  25. Pertea, M.; Pertea, G.M.; Antonescu, C.M.; Chang, T.C.; Mendell, J.T.; Salzberg, S.L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 2015, 33, 290–295. [Google Scholar] [CrossRef] [Green Version]
  26. Gehlenborg, N.; Noble, M.S.; Getz, G.; Chin, L.; Park, P.J. Nozzle: A report generation toolkit for data analysis pipelines. Bioinformatics 2013, 29, 1089–1091. [Google Scholar] [CrossRef]
  27. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [Green Version]
  28. Gentleman, R.; Carey, V.; Huber, W.; Hahne, F. Genefilter: Methods for Filtering Genes from High-Throughput Experiments, R Package Version 1.64.0. 2019. Available online: https://bioconductor.riken.jp/packages/3.8/bioc/html/genefilter.html (accessed on 22 April 2022).
  29. Langfelder, P.; Horvath, S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinform. 2008, 9, 559. [Google Scholar] [CrossRef] [Green Version]
  30. Chang, L.; Zhou, G.; Soufan, O.; Xia, J. miRNet 2.0: Network-based visual analytics for miRNA functional analysis and systems biology. Nucleic Acids Res. 2020, 48, W244–W251. [Google Scholar] [CrossRef]
  31. Andres-Leon, E.; Gomez-Lopez, G.; Pisano, D.G. Prediction of miRNA-mRNA Interactions Using miRGate. Methods Mol. Biol. 2017, 1580, 225–237. [Google Scholar] [CrossRef]
  32. Aguirre-Gamboa, R.; Gomez-Rueda, H.; Martinez-Ledesma, E.; Martinez-Torteya, A.; Chacolla-Huaringa, R.; Rodriguez-Barrientos, A.; Tamez-Pena, J.G.; Trevino, V. SurvExpress: An online biomarker validation tool and database for cancer gene expression data using survival analysis. PLoS ONE 2013, 8, e74250. [Google Scholar] [CrossRef] [Green Version]
  33. Shabalin, A.A. Matrix eQTL: Ultra fast eQTL analysis via large matrix operations. Bioinformatics 2012, 28, 1353–1358. [Google Scholar] [CrossRef] [Green Version]
  34. Xue, Y.; Ouyang, K.; Huang, J.; Zhou, Y.; Ouyang, H.; Li, H.; Wang, G.; Wu, Q.; Wei, C.; Bi, Y.; et al. Direct conversion of fibroblasts to neurons by reprogramming PTB-regulated microRNA circuits. Cell 2013, 152, 82–96. [Google Scholar] [CrossRef] [Green Version]
  35. Grosswendt, S.; Filipchyk, A.; Manzano, M.; Klironomos, F.; Schilling, M.; Herzog, M.; Gottwein, E.; Rajewsky, N. Unambiguous identification of miRNA:target site interactions by different types of ligation reactions. Mol. Cell 2014, 54, 1042–1054. [Google Scholar] [CrossRef] [Green Version]
  36. Kameswaran, V.; Bramswig, N.C.; McKenna, L.B.; Penn, M.; Schug, J.; Hand, N.J.; Chen, Y.; Choi, I.; Vourekas, A.; Won, K.J.; et al. Epigenetic regulation of the DLK1-MEG3 microRNA cluster in human type 2 diabetic islets. Cell Metab. 2014, 19, 135–145. [Google Scholar] [CrossRef] [Green Version]
  37. Helwak, A.; Kudla, G.; Dudnakova, T.; Tollervey, D. Mapping the human miRNA interactome by CLASH reveals frequent noncanonical binding. Cell 2013, 153, 654–665. [Google Scholar] [CrossRef] [Green Version]
  38. Balakrishnan, I.; Yang, X.; Brown, J.; Ramakrishnan, A.; Torok-Storb, B.; Kabos, P.; Hesselberth, J.R.; Pillai, M.M. Genome-wide analysis of miRNA-mRNA interactions in marrow stromal cells. Stem. Cells 2014, 32, 662–673. [Google Scholar] [CrossRef] [Green Version]
  39. Lipchina, I.; Elkabetz, Y.; Hafner, M.; Sheridan, R.; Mihailovic, A.; Tuschl, T.; Sander, C.; Studer, L.; Betel, D. Genome-wide identification of microRNA targets in human ES cells reveals a role for miR-302 in modulating BMP response. Genes Dev. 2011, 25, 2173–2186. [Google Scholar] [CrossRef] [Green Version]
  40. Karginov, F.V.; Hannon, G.J. Remodeling of Ago2-mRNA interactions upon cellular stress reflects miRNA complementarity and correlates with altered translation rates. Genes Dev. 2013, 27, 1624–1632. [Google Scholar] [CrossRef] [Green Version]
  41. Kishore, S.; Jaskiewicz, L.; Burger, L.; Hausser, J.; Khorshid, M.; Zavolan, M. A quantitative analysis of CLIP methods for identifying binding sites of RNA-binding proteins. Nat. Methods 2011, 8, 559–564. [Google Scholar] [CrossRef]
  42. Nelson, P.T.; Wang, W.X.; Mao, G.; Wilfred, B.R.; Xie, K.; Jennings, M.H.; Gao, Z.; Wang, X. Specific sequence determinants of miR-15/107 microRNA gene group targets. Nucleic Acids Res. 2011, 39, 8163–8172. [Google Scholar] [CrossRef]
  43. Wagschal, A.; Najafi-Shoushtari, S.H.; Wang, L.; Goedeke, L.; Sinha, S.; deLemos, A.S.; Black, J.C.; Ramirez, C.M.; Li, Y.; Tewhey, R.; et al. Genome-wide identification of microRNAs regulating cholesterol and triglyceride homeostasis. Nat. Med. 2015, 21, 1290–1297. [Google Scholar] [CrossRef] [Green Version]
  44. Adlakha, Y.K.; Khanna, S.; Singh, R.; Singh, V.P.; Agrawal, A.; Saini, N. Pro-apoptotic miRNA-128-2 modulates ABCA1, ABCG1 and RXRalpha expression and cholesterol homeostasis. Cell Death Dis. 2013, 4, e780. [Google Scholar] [CrossRef] [Green Version]
  45. Hedditch, E.L.; Gao, B.; Russell, A.J.; Lu, Y.; Emmanuel, C.; Beesley, J.; Johnatty, S.E.; Chen, X.; Harnett, P.; George, J.; et al. ABCA transporter gene expression and poor outcome in epithelial ovarian cancer. J. Natl. Cancer Inst. 2014, 106, dju149. [Google Scholar] [CrossRef]
  46. Zheng, L.; Li, L.; Lu, Y.; Jiang, F.; Yang, X.A. SREBP2 contributes to cisplatin resistance in ovarian cancer cells. Exp. Biol. Med. 2018, 243, 655–662. [Google Scholar] [CrossRef]
  47. Li, B.; Chen, H.; Wu, N.; Zhang, W.J.; Shang, L.X. Deregulation of miR-128 in ovarian cancer promotes cisplatin resistance. Int. J. Gynecol. Cancer 2014, 24, 1381–1388. [Google Scholar] [CrossRef]
  48. Amini-Farsani, Z.; Sangtarash, M.H.; Shamsara, M.; Teimori, H. MiR-221/222 promote chemoresistance to cisplatin in ovarian cancer cells by targeting PTEN/PI3K/AKT signaling pathway. Cytotechnology 2018, 70, 203–213. [Google Scholar] [CrossRef]
  49. Li, J.; Li, Q.; Huang, H.; Li, Y.; Li, L.; Hou, W.; You, Z. Overexpression of miRNA-221 promotes cell proliferation by targeting the apoptotic protease activating factor-1 and indicates a poor prognosis in ovarian cancer. Int. J. Oncol. 2017, 50, 1087–1096. [Google Scholar] [CrossRef] [Green Version]
  50. Sun, C.; Li, N.; Zhou, B.; Yang, Z.; Ding, D.; Weng, D.; Meng, L.; Wang, S.; Zhou, J.; Ma, D.; et al. miR-222 is upregulated in epithelial ovarian cancer and promotes cell proliferation by downregulating P27(kip1.). Oncol. Lett. 2013, 6, 507–512. [Google Scholar] [CrossRef]
  51. Kim, H.; Bhattacharya, A.; Qi, L. Endoplasmic reticulum quality control in cancer: Friend or foe. Semin. Cancer Biol. 2015, 33, 25–33. [Google Scholar] [CrossRef] [Green Version]
  52. Tian, J.; Liu, R.; Qu, Q. Role of endoplasmic reticulum stress on cisplatin resistance in ovarian carcinoma. Oncol. Lett. 2017, 13, 1437–1443. [Google Scholar] [CrossRef] [Green Version]
  53. Bastola, P.; Neums, L.; Schoenen, F.J.; Chien, J. VCP inhibitors induce endoplasmic reticulum stress, cause cell cycle arrest, trigger caspase-mediated cell death and synergistically kill ovarian cancer cells in combination with Salubrinal. Mol. Oncol. 2016, 10, 1559–1574. [Google Scholar] [CrossRef] [Green Version]
  54. Koundouros, N.; Poulogiannis, G. Reprogramming of fatty acid metabolism in cancer. Br. J. Cancer 2020, 122, 4–22. [Google Scholar] [CrossRef] [Green Version]
  55. Nie, L.Y.; Lu, Q.T.; Li, W.H.; Yang, N.; Dongol, S.; Zhang, X.; Jiang, J. Sterol regulatory element-binding protein 1 is required for ovarian tumor growth. Oncol. Rep. 2013, 30, 1346–1354. [Google Scholar] [CrossRef] [Green Version]
  56. Porstmann, T.; Griffiths, B.; Chung, Y.L.; Delpuech, O.; Griffiths, J.R.; Downward, J.; Schulze, A. PKB/Akt induces transcription of enzymes involved in cholesterol and fatty acid biosynthesis via activation of SREBP. Oncogene 2005, 24, 6465–6481. [Google Scholar] [CrossRef] [Green Version]
  57. Wang, H.Q.; Altomare, D.A.; Skele, K.L.; Poulikakos, P.I.; Kuhajda, F.P.; Di Cristofano, A.; Testa, J.R. Positive feedback regulation between AKT activation and fatty acid synthase expression in ovarian carcinoma cells. Oncogene 2005, 24, 3574–3582. [Google Scholar] [CrossRef] [Green Version]
  58. Wheeler, L.J.; Watson, Z.L.; Qamar, L.; Yamamoto, T.M.; Sawyer, B.T.; Sullivan, K.D.; Khanal, S.; Joshi, M.; Ferchaud-Roucher, V.; Smith, H.; et al. Multi-Omic Approaches Identify Metabolic and Autophagy Regulators Important in Ovarian Cancer Dissemination. iScience 2019, 19, 474–491. [Google Scholar] [CrossRef] [Green Version]
  59. Pinto, B.A.S.; Franca, L.M.; Laurindo, F.R.M.; Paes, A.M.A. Unfolded Protein Response: Cause or Consequence of Lipid and Lipoprotein Metabolism Disturbances? Adv. Exp. Med. Biol. 2019, 1127, 67–82. [Google Scholar] [CrossRef]
  60. Hao, D.; Liu, J.; Chen, M.; Li, J.; Wang, L.; Li, X.; Zhao, Q.; Di, L.J. Immunogenomic Analyses of Advanced Serous Ovarian Cancer Reveal Immune Score is a Strong Prognostic Factor and an Indicator of Chemosensitivity. Clin. Cancer Res. 2018, 24, 3560–3571. [Google Scholar] [CrossRef] [Green Version]
  61. Mahadevan, N.R.; Rodvold, J.; Sepulveda, H.; Rossi, S.; Drew, A.F.; Zanetti, M. Transmission of endoplasmic reticulum stress and pro-inflammation from tumor cells to myeloid cells. Proc. Natl. Acad. Sci. USA 2011, 108, 6561–6566. [Google Scholar] [CrossRef] [Green Version]
  62. Cubillos-Ruiz, J.R.; Silberman, P.C.; Rutkowski, M.R.; Chopra, S.; Perales-Puchalt, A.; Song, M.; Zhang, S.; Bettigole, S.E.; Gupta, D.; Holcomb, K.; et al. ER Stress Sensor XBP1 Controls Anti-tumor Immunity by Disrupting Dendritic Cell Homeostasis. Cell 2015, 161, 1527–1538. [Google Scholar] [CrossRef] [Green Version]
  63. Herber, D.L.; Cao, W.; Nefedova, Y.; Novitskiy, S.V.; Nagaraj, S.; Tyurin, V.A.; Corzo, A.; Cho, H.I.; Celis, E.; Lennox, B.; et al. Lipid accumulation and dendritic cell dysfunction in cancer. Nat. Med. 2010, 16, 880–886. [Google Scholar] [CrossRef] [Green Version]
  64. Dai, R.; Li, J.; Liu, Y.; Yan, D.; Chen, S.; Duan, C.; Liu, X.; He, T.; Li, H. miR-221/222 suppression protects against endoplasmic reticulum stress-induced apoptosis via p27(Kip1)- and MEK/ERK-mediated cell cycle regulation. Biol. Chem. 2010, 391, 791–801. [Google Scholar] [CrossRef]
  65. Buniello, A.; MacArthur, J.A.L.; Cerezo, M.; Harris, L.W.; Hayhurst, J.; Malangone, C.; McMahon, A.; Morales, J.; Mountjoy, E.; Sollis, E.; et al. The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic Acids Res. 2019, 47, D1005–D1012. [Google Scholar] [CrossRef] [Green Version]
  66. Zhang, Y.; Wu, J.; Liang, J.Y.; Huang, X.; Xia, L.; Ma, D.W.; Xu, X.Y.; Wu, P.P. Association of serum lipids and severity of epithelial ovarian cancer: An observational cohort study of 349 Chinese patients. J. Biomed. Res. 2018, 32, 336–342. [Google Scholar] [CrossRef]
  67. Li, A.J.; Elmore, R.G.; Chen, I.Y.; Karlan, B.Y. Serum low-density lipoprotein levels correlate with survival in advanced stage epithelial ovarian cancers. Gynecol. Oncol. 2010, 116, 78–81. [Google Scholar] [CrossRef] [Green Version]
  68. Cancer Genome Atlas Research Network. Integrated genomic analyses of ovarian carcinoma. Nature 2011, 474, 609–615. [Google Scholar] [CrossRef]
  69. Anders, S.; Huber, W. Differential expression analysis for sequence count data. Genome Biol. 2010, 11, R106. [Google Scholar] [CrossRef] [Green Version]
  70. Dillies, M.A.; Rau, A.; Aubert, J.; Hennequet-Antier, C.; Jeanmougin, M.; Servant, N.; Keime, C.; Marot, G.; Castel, D.; Estelle, J.; et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief. Bioinform. 2013, 14, 671–683. [Google Scholar] [CrossRef] [Green Version]
  71. Costa-Silva, J.; Domingues, D.; Lopes, F.M. RNA-Seq differential expression analysis: An extended review and a software tool. PLoS ONE 2017, 12, e0190152. [Google Scholar] [CrossRef] [Green Version]
  72. Raudvere, U.; Kolberg, L.; Kuzmin, I.; Arak, T.; Adler, P.; Peterson, H.; Vilo, J. g:Profiler: A web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 2019, 47, W191–W198. [Google Scholar] [CrossRef] [Green Version]
  73. Bagnoli, M.; De Cecco, L.; Granata, A.; Nicoletti, R.; Marchesi, E.; Alberti, P.; Valeri, B.; Libra, M.; Barbareschi, M.; Raspagliesi, F.; et al. Identification of a chrXq27.3 microRNA cluster associated with early relapse in advanced stage ovarian cancer patients. Oncotarget 2011, 2, 1265–1278. [Google Scholar] [CrossRef] [Green Version]
  74. Tothill, R.W.; Tinker, A.V.; George, J.; Brown, R.; Fox, S.B.; Lade, S.; Johnson, D.S.; Trivett, M.K.; Etemadmoghadam, D.; Locandro, B.; et al. Novel molecular subtypes of serous and endometrioid ovarian cancer linked to clinical outcome. Clin. Cancer Res. 2008, 14, 5198–5208. [Google Scholar] [CrossRef] [Green Version]
  75. Gyorffy, B.; Lanczky, A.; Szallasi, Z. Implementing an online tool for genome-wide validation of survival-associated biomarkers in ovarian-cancer using microarray data from 1287 patients. Endocr. Relat. Cancer 2012, 19, 197–208. [Google Scholar] [CrossRef] [Green Version]
  76. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2019; Available online: https://www.R-project.org/ (accessed on 20 April 2022).
  77. RStudio Team. RStudio: Integrated Development Environment for R; RStudio, Inc.: Boston, MA, USA, 2016; Available online: http://www.rstudio.com/ (accessed on 20 April 2022).
  78. Wickham, H. Ggplot2: Elegrant graphics for data analysis. In Use R! 2nd ed.; Springer: Cham, Switzerland, 2016; p. 1, online resource. [Google Scholar]
  79. Slowikowski, K.; Schep, A.; Hughes, S.; Dang, K.T.; Lukauskas, S.; Irisson, J.O.; Kamvar, Z.N.; Ryan, T.; Christophe, D.; Hiroaki, Y.; et al. Automatically Position Non-Overlapping Text Labels with ‘ggplot2’. R Package Version 0.9.1. 2021. Available online: https://CRAN.R-project.org/package=ggrepel (accessed on 29 March 2022).
  80. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [Google Scholar] [CrossRef]
  81. Therneau, T.M. A Package for Survival Analysis in R. R Package Version 3.2-7. 2020. Available online: https://CRAN.R-project.org/package=survival (accessed on 20 April 2022).
  82. Harrell, F.E., Jr. rms: Regression Modeling Strategies. R Package Version 6.1-1. 2021. Available online: https://CRAN.R-project.org/package=rms (accessed on 20 April 2022).
  83. Colaprico, A.; Silva, T.C.; Olsen, C.; Garofano, L.; Cava, C.; Garolini, D.; Sabedot, T.S.; Malta, T.M.; Pagnotta, S.M.; Castiglioni, I.; et al. TCGAbiolinks: An R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016, 44, e71. [Google Scholar] [CrossRef]
  84. Ewing, B.; Green, P. Base-calling of automated sequencer traces using phred. II. Error probabilities. Genome Res. 1998, 8, 186–194. [Google Scholar] [CrossRef] [Green Version]
  85. Van der Auwera, G.A.; Carneiro, M.O.; Hartl, C.; Poplin, R.; Del Angel, G.; Levy-Moonshine, A.; Jordan, T.; Shakir, K.; Roazen, D.; Thibault, J.; et al. From FastQ data to high confidence variant calls: The Genome Analysis Toolkit best practices pipeline. Curr. Protoc. Bioinform. 2013, 43, 11.10.1–11.10.33. [Google Scholar] [CrossRef]
  86. Zerbino, D.R.; Achuthan, P.; Akanni, W.; Amode, M.R.; Barrell, D.; Bhai, J.; Billis, K.; Cummins, C.; Gall, A.; Giron, C.G.; et al. Ensembl 2018. Nucleic Acids Res. 2018, 46, D754–D761. [Google Scholar] [CrossRef]
  87. Chu, A.; Robertson, G.; Brooks, D.; Mungall, A.J.; Birol, I.; Coope, R.; Ma, Y.; Jones, S.; Marra, M.A. Large-scale profiling of microRNAs for The Cancer Genome Atlas. Nucleic Acids Res. 2016, 44, e3. [Google Scholar] [CrossRef]
  88. Tukey, J.W. Exploratory Data Analysis; Addison-Wesley Publishing Company: Boston, MA, USA, 1977. [Google Scholar]
  89. Hoaglin, D.C.; Iglewicz, B. Fine-Tuning Some Resistant Rules for Outlier Labeling. J. Am. Stat. Assoc. 1987, 82, 1147–1149. [Google Scholar] [CrossRef]
  90. Panarelli, N.; Tyryshkin, K.; Wong, J.J.M.; Majewski, A.; Yang, X.; Scognamiglio, T.; Kim, M.K.; Bogardus, K.; Tuschl, T.; Chen, Y.T.; et al. Evaluating gastroenteropancreatic neuroendocrine tumors through microRNA sequencing. Endocr. Relat. Cancer 2019, 26, 47–57. [Google Scholar] [CrossRef] [PubMed]
  91. Langfelder, P.; Horvath, S. Tutorial for the WGCNA Package for R: I. Network Analysis of Liver Expression Data in Female Mice 2.b Step-by-Step Network Construction and Module Detection. Available online: https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/FemaleLiver-02-networkConstr-man.pdf (accessed on 22 April 2022).
  92. Chang, C.C.; Chow, C.C.; Tellier, L.C.; Vattikuti, S.; Purcell, S.M.; Lee, J.J. Second-generation PLINK: Rising to the challenge of larger and richer datasets. Gigascience 2015, 4, 7. [Google Scholar] [CrossRef] [PubMed]
  93. Stranger, B.E.; Montgomery, S.B.; Dimas, A.S.; Parts, L.; Stegle, O.; Ingle, C.E.; Sekowska, M.; Smith, G.D.; Evans, D.; Gutierrez-Arcelus, M.; et al. Patterns of cis regulatory variation in diverse human populations. PLoS Genet. 2012, 8, e1002639. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  94. Hinrichs, A.S.; Karolchik, D.; Baertsch, R.; Barber, G.P.; Bejerano, G.; Clawson, H.; Diekhans, M.; Furey, T.S.; Harte, R.A.; Hsu, F.; et al. The UCSC Genome Browser Database: Update 2006. Nucleic Acids Res. 2006, 34, D590–D598. [Google Scholar] [CrossRef] [Green Version]
  95. Sherry, S.T.; Ward, M.H.; Kholodov, M.; Baker, J.; Phan, L.; Smigielski, E.M.; Sirotkin, K. dbSNP: The NCBI database of genetic variation. Nucleic Acids Res. 2001, 29, 308–311. [Google Scholar] [CrossRef] [Green Version]
  96. Ward, L.D.; Kellis, M. HaploReg v4: Systematic mining of putative causal variants, cell types, regulators and target genes for human complex traits and disease. Nucleic Acids Res. 2016, 44, D877–D881. [Google Scholar] [CrossRef]
  97. Kauffmann, A.; Gentleman, R.; Huber, W. arrayQualityMetrics—A bioconductor package for quality assessment of microarray data. Bioinformatics 2009, 25, 415–416. [Google Scholar] [CrossRef] [Green Version]
  98. Kassambara, A.; Kosinski, M.; Biecek, P. Survminer: Drawing Survival Curves Using ‘ggplot2’, R Package Version 0.4.8. 2020. Available online: https://cran.r-project.org/web/packages/survminer/index.html (accessed on 20 April 2022).
Ijms 23 04875 g001 550
Figure 1. Flowchart of study analysis pipeline. First, reports of read-sequence quality were generated using FastQC [20] and MultiQC [21], followed by trimming of low-quality mRNA-seq reads with Trimmomatic [22]. Filtered reads were aligned to the hg19 human reference genome [23] using HISAT2 (hierarchical indexing for spliced alignment of transcripts) [24] and quantified by StringTie [25]. miRNA-seq data were obtained as quantified reads from the Cancer Genome Atlas [26]. Next, transcript expression was normalized using the DESeq2 R package [27], and highly variable transcripts were selected with the varFilter function of the genefilter R package [28]. Moreover, transcript expression was used to test for differential expression using DESeq2 and to construct coexpression networks using the weighted correlation network analysis (WGCNA) R package [29]. Interactions between miRNAs and mRNAs from significant WGCNA networks were detected using Sperman’s correlation, and validated using miRNet [30] and miRGate [31]. Validation of differential expression analysis results was performed using Kaplan–Meier analysis on two independent HGSOC cohorts. Validation of network analysis results was performed using the Prognostic Index estimation method from the SurvExpress tool [32] on the same validation cohorts. Transcript expression data were integrated with genomics data from the same patients to determine eQTLs (expression quantitative trait loci) using the MatrixEQTL R package [33].
Figure 1. Flowchart of study analysis pipeline. First, reports of read-sequence quality were generated using FastQC [20] and MultiQC [21], followed by trimming of low-quality mRNA-seq reads with Trimmomatic [22]. Filtered reads were aligned to the hg19 human reference genome [23] using HISAT2 (hierarchical indexing for spliced alignment of transcripts) [24] and quantified by StringTie [25]. miRNA-seq data were obtained as quantified reads from the Cancer Genome Atlas [26]. Next, transcript expression was normalized using the DESeq2 R package [27], and highly variable transcripts were selected with the varFilter function of the genefilter R package [28]. Moreover, transcript expression was used to test for differential expression using DESeq2 and to construct coexpression networks using the weighted correlation network analysis (WGCNA) R package [29]. Interactions between miRNAs and mRNAs from significant WGCNA networks were detected using Sperman’s correlation, and validated using miRNet [30] and miRGate [31]. Validation of differential expression analysis results was performed using Kaplan–Meier analysis on two independent HGSOC cohorts. Validation of network analysis results was performed using the Prognostic Index estimation method from the SurvExpress tool [32] on the same validation cohorts. Transcript expression data were integrated with genomics data from the same patients to determine eQTLs (expression quantitative trait loci) using the MatrixEQTL R package [33].
Ijms 23 04875 g001
Ijms 23 04875 g002 550
Figure 2. Differential expression analysis of miRNA and mRNA transcripts. (A) Significant differentially expressed miRNA isoforms (n = 21) are shown in red. The Benjamini–Hochberg adjusted significance threshold is indicated by the dashed line. Transcripts with a positive fold change are upregulated in resistant patients, whereas transcripts with a negative fold change are downregulated in resistant patients. (B) Significant differentially expressed mRNA transcripts (n = 196) are shown in red. The Benjamini–Hochberg adjusted significance threshold is indicated by the dashed line. Transcripts with a positive fold change are upregulated in resistant patients, whereas transcripts with a negative fold change are downregulated in resistant patients. Gene symbol labels are provided for the top 20 most significant mRNAs, as well as for mRNAs that were replicated in an independent cohort (Results Section 2.7, Table S11).
Figure 2. Differential expression analysis of miRNA and mRNA transcripts. (A) Significant differentially expressed miRNA isoforms (n = 21) are shown in red. The Benjamini–Hochberg adjusted significance threshold is indicated by the dashed line. Transcripts with a positive fold change are upregulated in resistant patients, whereas transcripts with a negative fold change are downregulated in resistant patients. (B) Significant differentially expressed mRNA transcripts (n = 196) are shown in red. The Benjamini–Hochberg adjusted significance threshold is indicated by the dashed line. Transcripts with a positive fold change are upregulated in resistant patients, whereas transcripts with a negative fold change are downregulated in resistant patients. Gene symbol labels are provided for the top 20 most significant mRNAs, as well as for mRNAs that were replicated in an independent cohort (Results Section 2.7, Table S11).
Ijms 23 04875 g002
Ijms 23 04875 g003 550
Figure 3. Weighted correlation network analysis of mRNA and miRNA transcripts. (A) The miRNA modules lightcoral, plum, and ivory are visualized in their respective colors. Each node represents one miRNA isoform, and each edge represents a connection or coexpression. The distance between nodes is the connection weight, where more similarly expressed transcripts are plotted closer. (B) The mRNA modules lavenderblush3 and darkolivegreen are visualized in their respective colors. Each node represents a transcript, and each edge represents a connection or coexpression. The distance between nodes is the connection weight, where more similar transcripts are plotted closer.
Figure 3. Weighted correlation network analysis of mRNA and miRNA transcripts. (A) The miRNA modules lightcoral, plum, and ivory are visualized in their respective colors. Each node represents one miRNA isoform, and each edge represents a connection or coexpression. The distance between nodes is the connection weight, where more similarly expressed transcripts are plotted closer. (B) The mRNA modules lavenderblush3 and darkolivegreen are visualized in their respective colors. Each node represents a transcript, and each edge represents a connection or coexpression. The distance between nodes is the connection weight, where more similar transcripts are plotted closer.
Ijms 23 04875 g003
Ijms 23 04875 g004 550
Figure 4. Integration of significant networks. The mRNA transcripts (oval nodes) from the lavenderblush3 and darkolivegreen modules are arranged based on coexpression similarity (grey edges). Regulatory interactions (black edges) indicate inhibition (bar-headed arrow) or activation (arrow). The colors of miRNA nodes indicate their membership in the plum and ivory coexpression networks. Four miRNA isoforms (rectangular nodes) were identified, either using experimental methods (solid edges) or in silico predictions (dashed edges), to regulate the expression of genes in the lavenderblush3 network module. In addition, 6 eQTLs (hexagonal nodes) may regulate the expression of the NFX1 gene. Genes in the darkolivegreen module have 10 validated or predicted regulatory miRNAs and 97 regulatory eQTLs. Finally, the ivory miRNA miR-140 may be regulated by the miR-QTL rs71397980. eQTL, expression quantitative trait loci; miR-QTL, microRNA expression quantitative trait loci.
Figure 4. Integration of significant networks. The mRNA transcripts (oval nodes) from the lavenderblush3 and darkolivegreen modules are arranged based on coexpression similarity (grey edges). Regulatory interactions (black edges) indicate inhibition (bar-headed arrow) or activation (arrow). The colors of miRNA nodes indicate their membership in the plum and ivory coexpression networks. Four miRNA isoforms (rectangular nodes) were identified, either using experimental methods (solid edges) or in silico predictions (dashed edges), to regulate the expression of genes in the lavenderblush3 network module. In addition, 6 eQTLs (hexagonal nodes) may regulate the expression of the NFX1 gene. Genes in the darkolivegreen module have 10 validated or predicted regulatory miRNAs and 97 regulatory eQTLs. Finally, the ivory miRNA miR-140 may be regulated by the miR-QTL rs71397980. eQTL, expression quantitative trait loci; miR-QTL, microRNA expression quantitative trait loci.
Ijms 23 04875 g004
Ijms 23 04875 g005 550
Figure 5. Replication analysis of significant transcript coexpression networks. (A–C) Kaplan–Meier analysis of patient PFS association with the prognostic index of miRNA networks ivory, plum, and lightcoral in the MITO cohort. The ivory (A) and plum (B) miRNA network modules replicated in the MITO cohort (p = 6.1 × 10−4, log HR = −0.78 and p = 0.022, log HR = −0.52, respectively), while the lightcoral module (C) reached nominal significance (p = 0.057, log HR = −0.43). (D,E) Kaplan–Meier analysis of patient PFS association with the prognostic index of mRNA networks lavenderblush3 and darkolivegreen in the AOCS cohort. The lavenderblush3 (D) and darkolivegreen (E) mRNA network modules replicated in the AOCS cohort (p = 1.6 × 10−10, log HR = −1.27 and p = 1.3 × 10−10, log HR = −1.27, respectively). PFS, progression-free survival; AOCS, Australian Ovarian Cancer Study; MITO, Multicenter Italian Trial in Ovarian cancer.
Figure 5. Replication analysis of significant transcript coexpression networks. (A–C) Kaplan–Meier analysis of patient PFS association with the prognostic index of miRNA networks ivory, plum, and lightcoral in the MITO cohort. The ivory (A) and plum (B) miRNA network modules replicated in the MITO cohort (p = 6.1 × 10−4, log HR = −0.78 and p = 0.022, log HR = −0.52, respectively), while the lightcoral module (C) reached nominal significance (p = 0.057, log HR = −0.43). (D,E) Kaplan–Meier analysis of patient PFS association with the prognostic index of mRNA networks lavenderblush3 and darkolivegreen in the AOCS cohort. The lavenderblush3 (D) and darkolivegreen (E) mRNA network modules replicated in the AOCS cohort (p = 1.6 × 10−10, log HR = −1.27 and p = 1.3 × 10−10, log HR = −1.27, respectively). PFS, progression-free survival; AOCS, Australian Ovarian Cancer Study; MITO, Multicenter Italian Trial in Ovarian cancer.
Ijms 23 04875 g005
Table
Table 1. Predicted and validated mRNA—miRNA network interactions.
Table 1. Predicted and validated mRNA—miRNA network interactions.
Gene IDmRNA IDmiRNA IDSpearman’s ρp ValueInteraction Evidence
lavenderblush3 mRNA—plum miRNA network interactions
PSIP1ENST00000397519hsa-mir-222-3p−0.320.000021Experimental validation [34]
CDC37L1ENST00000381854hsa-mir-222-3p−0.310.000042Computational prediction
CDC37L1ENST00000381854hsa-mir-221-3p−0.270.00034Computational prediction
CDC37L1ENST00000381854hsa-mir-221-5p−0.270.00034Computational prediction
NFIBENST00000380959hsa-mir-221-3p−0.270.00039Experimental validation [35]
SMU1ENST00000397149hsa-mir-221-5p−0.260.00082Experimental validation [36]
VCPENST00000358901hsa-mir-222-3p−0.250.00089Experimental validation [37]
CAAP1ENST00000333916hsa-mir-222-5p−0.220.0041Computational prediction
CAAP1ENST00000333916hsa-mir-221-5p−0.210.0076Computational prediction
TOPORSENST00000360538hsa-mir-221-3p−0.190.014Experimental validation [38]
SLC25A51ENST00000496760hsa-mir-222-3p−0.190.016Experimental validation [37]
PPAPDC2ENST00000381883hsa-mir-222-5p−0.180.021Computational prediction
VCPENST00000358901hsa-mir-221-3p−0.170.023Experimental validation [39]
VCPENST00000358901hsa-mir-221-5p−0.170.023Experimental validation [38]
NFX1ENST00000379540hsa-mir-221-5p−0.170.027Computational prediction
PPAPDC2ENST00000381883hsa-mir-221-5p−0.170.03Computational prediction
DNAJA1ENST00000330899hsa-mir-221-5p−0.160.035Experimental validation [38]
AK3ENST00000381809hsa-mir-222-5p−0.160.037Computational prediction
darkolivegreen mRNA—ivory miRNA interactions
RPL18ENST00000552347hsa-mir-212-3p−0.290.00017Experimental Validation [40]
MRPL55ENST00000411464hsa-miR-1306-5p−0.250.001Computational prediction
FASTKENST00000466855hsa-mir-140-3p−0.250.0013Experimental Validation [41]
RPL18ENST00000552347hsa-mir-361-3p−0.240.0021Experimental Validation [40]
ACBD4ENST00000321854hsa-miR-128-2-5p−0.220.0041Computational prediction
OAZ1ENST00000581150hsa-mir-361-3p−0.220.0047Experimental Validation [36,40]
ACBD4ENST00000321854hsa-miR-128-1-5p−0.210.0059Computational prediction
MYL12AENST00000578611hsa-mir-103a-3p−0.170.023Experimental Validation [38,40]
MYL12AENST00000578611hsa-mir-107−0.170.023Experimental Validation [38,40,42]
HTRA2ENST00000484352hsa-mir-1306-5p−0.170.029Experimental Validation [40]
ACADVLENST00000579425hsa-mir-103a-3p−0.170.03Experimental Validation [38]
ACADVLENST00000579425hsa-mir-107−0.170.03Experimental Validation [38]
SOD1ENST00000470944hsa-mir-140-3p−0.160.037Experimental Validation [34]
RHOT2ENST00000569675hsa-mir-1306-5p−0.160.044Experimental Validation [40]
TSSK6ENST00000360913hsa-miR-212-5p−0.150.048Computational prediction
RCC1ENST00000373832hsa-miR-1306-3p−0.150.048Computational prediction
RCC1ENST00000373832hsa-miR-1306-5p−0.150.048Computational prediction
Both (100%) of the miRNAs in the plum network significantly interact with 12 of the 31 (39%) genes in the lavenderblush mRNA network. A total of 8 of the 11 (73%) miRNAs in the ivory network significantly interact with 12 of the 80 (15%) genes in the darkolivegreen mRNA network.
Table
Table 2. Characteristics of the HGSOC cohorts with mRNA-Seq and miRNA-Seq data from TCGA.
Table 2. Characteristics of the HGSOC cohorts with mRNA-Seq and miRNA-Seq data from TCGA.
mRNA-Seq CohortmiRNA-Seq Cohort
SensitiveResistantp ValueAll CasesSensitiveResistantp ValueAll Cases
Age
Mean57.262.00.0043 a59.257.861.50.018 a59.3
Range30–8738–87 30–8730–8738–87 30–87
Grade
G218 (15.9%)8 (10.3%)0.27 b26 (13.6%)23 (18.7%)8 (9.8%)0.14 b31 (15.1%)
G3/495 (84.1%)69 (88.5%) 164 (85.9%)99 (80.5%)73 (89.0%) 172 (83.9%)
Ungraded0 (0.0%)1 (1.3%) 1 (0.5%)1 (0.8%)1 (1.2%) 2 (1.0%)
Stage
II7 (6.2%)2 (2.6%)0.48 b9 (4.7%)8 (6.5%)2 (2.4%)0.41 b10 (4.9%)
III91 (80.5%)67 (85.9%) 158 (82.7%)99 (80.5%)67 (81.7%) 166 (81.0%)
IV15 (13.3%)9 (11.5%) 24 (12.6%)16 (13.0%)13 (15.9%) 29 (14.2%)
Cytoreductive surgery outcome
Optimal (≤10 mm)73 (64.6%)54 (69.2%)0.23 c127 (66.5%)77 (62.6%)56 (68.3%)0.21 c133 (64.9%)
Suboptimal (>10 mm)26 (23.0%)20 (25.6%) 46 (24.1%)29 (23.6%)21 (25.6%) 50 (24.4%)
Unknown14 (12.4%)4 (5.1%) 18 (9.4%)17 (13.8%)5 (6.1%) 22 (10.7%)
Overall survival
Mortality events48 (42.5%)60 (76.9%)2.35 × 10−6 c108 (56.5%)55 (44.7%)66 (80.5%)0.0002 c121 (59.0%)
Median months52.028.26.97 × 10−10 a39.951.928.21.13 × 10−10 a39.5
95% CI50.1–54.026.2–30.2 37.9–41.949.9–53.926.2–30.2 37.5–41.5
Adjuvant chemotherapy regimen
Platinum agent only4 (3.5%)3 (3.9%)1 b7 (3.7%)4 (3.3%)4 (4.9%)0.71 b8 (3.9%)
Platinum and Taxane combination109 (96.5%)75 (96.2%) 184 (96.3%)119 (96.8%)78 (95.1%) 197 (96.1%)
Primary tumor mRNA subtype
Differentiated37 (32.7%)22 (28.2%)0.49 b59 (30.9%)42 (34.2%)24 (29.3%)0.68 c66 (32.2%)
Immunoreactive23 (20.4%)10 (12.8%) 33 (17.3%)24 (19.5%)13 (15.9%) 37 (18.1%)
Mesenchymal22 (19.5%)18 (23.1%) 40 (20.9%)24 (19.5%)18 (22.0%) 42 (20.5%)
Proliferative30 (26.6%)27 (34.6%) 57 (29.8%)33 (26.8%)27 (32.9%) 60 (29.3%)
Unknown1 (0.9%)1 (1.3%) 2 (1.1%)0 (0.0%)0 (0.0%) 0 (0.0%)
Total patients11378 19112382 205
a—p-value based on t-test; b—p-value based on Fisher’s Exact test; c—p-value based on Chi-Squared test; CI, Confidence Interval.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Topouza, D.G.; Choi, J.; Nesdoly, S.; Tarnouskaya, A.; Nicol, C.J.B.; Duan, Q.L. Novel MicroRNA-Regulated Transcript Networks Are Associated with Chemotherapy Response in Ovarian Cancer. Int. J. Mol. Sci. 2022, 23, 4875. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms23094875

AMA Style

Topouza DG, Choi J, Nesdoly S, Tarnouskaya A, Nicol CJB, Duan QL. Novel MicroRNA-Regulated Transcript Networks Are Associated with Chemotherapy Response in Ovarian Cancer. International Journal of Molecular Sciences. 2022; 23(9):4875. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms23094875

Chicago/Turabian Style

Topouza, Danai G., Jihoon Choi, Sean Nesdoly, Anastasiya Tarnouskaya, Christopher J. B. Nicol, and Qing Ling Duan. 2022. "Novel MicroRNA-Regulated Transcript Networks Are Associated with Chemotherapy Response in Ovarian Cancer" International Journal of Molecular Sciences 23, no. 9: 4875. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms23094875

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop