Construction of miRNA‐lncRNA‐mRNA co‐expression network affecting EMT‐mediated cisplatin resistance in ovarian cancer

Abstract Platinum resistance is one of the major concerns in ovarian cancer treatment. Recent evidence shows the critical role of epithelial–mesenchymal transition (EMT) in this resistance. Epithelial‐like ovarian cancer cells show decreased sensitivity to cisplatin after cisplatin treatment. Our study prospected the association between epithelial phenotype and response to cisplatin in ovarian cancer. Microarray dataset GSE47856 was acquired from the GEO database. After identifying differentially expressed genes (DEGs) between epithelial‐like and mesenchymal‐like cells, the module identification analysis was performed using weighted gene co‐expression network analysis (WGCNA). The gene ontology (GO) and pathway analyses of the most considerable modules were performed. The protein–protein interaction network was also constructed. The hub genes were specified using Cytoscape plugins MCODE and cytoHubba, followed by the survival analysis and data validation. Finally, the co‐expression of miRNA‐lncRNA‐TF with the hub genes was reconstructed. The co‐expression network analysis suggests 20 modules relating to the Epithelial phenotype. The antiquewhite4, brown and darkmagenta modules are the most significant non‐preserved modules in the Epithelial phenotype and contain the most differentially expressed genes. GO, and KEGG pathway enrichment analyses on these modules divulge that these genes were primarily enriched in the focal adhesion, DNA replication pathways and stress response processes. ROC curve and overall survival rate analysis show that the co‐expression pattern of the brown module's hub genes could be a potential prognostic biomarker for ovarian cancer cisplatin resistance.


| INTRODUC TI ON
Ovarian cancer is one of the deadliest gynaecological malignancies with a low survival rate (Near 15% 5-year survival for stage IV). 1 Ovarian cancers are classified as non-epithelial and epithelial ovarian cancers (EOC). EOC is associated with ovarian cancer-related deaths. 2 About 25% of patients with ovarian cancer are resistant to platinum-based therapy. 3 Furthermore, about 80% of patients suffer from recurrence of ovarian cancer, and these tumours are typically platinum-resistant, which leads to chemotherapy failure. 3,4 So, it is crucial to overcome this resistance in ovarian cancer cells. Cisdiamminedichloroplatinum (cisplatin) is a platinum-based chemotherapy drug that treats various cancers like bladder, ovarian, lung, breast and brain cancers. 5,6 Cisplatin forms DNA adducts involved in activating DNA damage recognition, DNA repair and apoptosis signalling pathways. 7 Two mechanisms have been suggested for platinum resistance in cancer cells. In the first mechanism, cisplatin uptake decreases while its detoxification increases. The second mechanism is the activation of anti-apoptotic pathways like NF-κB and MAPK pathways after treatment. [8][9][10][11] Recent evidence shows the association of EMT with cellular resistance to cisplatin. During the EMT process, epithelial cells lose the polarized epithelial structure and transform into moving mesenchymal cells, increasing their invasiveness. It has been shown that EMT is associated with drug resistance and apoptosis scape in a variety of cancer types. [12][13][14][15] Many observations indicate a link between drug resistance and EMT in various cancers like colorectal, 16 breast 17,18 and ovarian. [19][20][21][22][23][24][25] Furthermore, Miow et al. 7 discovered that epithelial-like ovarian cancer cell lines exhibit resistance to cisplatin treatment, along with NF-κB activation and apoptotic impairment. Identification of molecular mechanisms which are involved in this process would be helpful.
Today, co-expression network analysis is used for underlying the regulatory mechanisms relating to the specific biological processes. The WGCNA is a powerful method for analysing gene expression data, discovering modules of highly related genes and connecting each module to sample traits. 26,27 This tool constructs a co-expression network based on the expression profile similarities in samples. 28 WGCNA hierarchical clustering methods use holistic gene expression information to discover gene network signatures in a phenotype, which helps us to reduce bias. 29 In this study, we applied WGCNA to the expression profile of cisplatin-treated ovarian cancer cell lines to identify critical modules in treated cell lines with epithelial status compared to mesenchymal cell lines. These modules were closely associated with cisplatin resistance in ovarian cancer cell lines. The analysis of co-expression networks may decipher new insights into molecular mechanisms and signalling pathways of drug resistance in ovarian cancer to improve its prognosis and treatment.

| Acquisition of microarray datasets
The flow chart of our study is shown in Figure 1. Raw CEL files of Microarray dataset GSE47856 from the NCBI Gene Expression Omnibus (GEO) database were collected. 7 The dataset platform was

| Microarray data preprocessing and DEGs identification
A total of 48 samples in CEL format were simultaneously preprocessed using Robust Multi-chip Analysis (RMA) function for background adjustment, quantile normalization and summarization. 30,31 To come by the highest possible level of data quality and to eliminate mistargeted and nonspecific probes on the microarrays, the Principal Component Analysis (PCA) was used to identify and remove outlier samples from the dataset. 32 The multiple probes measure the expression of a given gene, and it is necessary to collapse the multiple probe sets to the same gene by applying the collapseRows function, which was reported as an effective method previously (MaxMean was used for collapsing rows). 33 Finally, 20,849 genes were used as input in the DEGs and co-expression network analyses. The limma package was used for DEGs analysis through linear modelling and empirical Bayes methods. 34 The criteria considered for DEG extraction were as |log2 fold change| ≥ 0.58 and adjusted p-value <0.05.

| Construction of a signed-hybrid weighted gene co-expression network
To construct the co-expression networks, we used the WGCNA package. A signed-hybrid weighted gene co-expression network was built based on Mesenchymal and Epithelial gene expression. The pickSoftThreshold function of the WGCNA package was used to set soft threshold power β as a tradeoff between scale-free topology and mean connectivity for Mesenchymal and Epithelial.

| Generating adjacency and TOM similarity matrices
Based on the selected soft-power, calculation of the adjacency matrix into a topological overlap matrix (TOM) was accomplished to minimize the effects of noise and spurious associations. Based on the TOM dissimilarity, hierarchical clustering was exerted to classify highly co-expressed genes as dense interconnected branches of the tree (dendrogram) into the same modules and extracted through the dynamic hybrid tree-cutting algorithm. Modules with high eigengene correlation were merged using the mergeCloseModules function (cutHeight = 0.3, corresponding to correlation of 0.7). 35 Application of eigengenes in WGCNA would be in modules summarization and measuring module memberships (kME) to earn suitable target genes.
These genes are recognized as connectivity based on eigengenes, which are calculated by the moduleEigengenes function. Eigengenes would be conducted as the first principal component of each module leading to the weighted average of the module's co-expression profiles by summarizing and comparing them. 36

| Construction of protein-protein interaction Network
The protein-protein interaction (PPI) network of the non-preserved module was constructed using the STRING tool in Cytoscape version 3.7.2. 40 The potential interaction between genes at the protein level predicts protein interactions and the weight of each edge (line) in the PPI network. The criteria used for the network is a combined score greater than 0.4. The molecular complex detection (MCODE) analysis was used to identify the significant clusters of the candidate module PPI network with degree cut-off 2, max depth 100, k-core 2 and node score cut-off 0.2. 41 MCODE's highest-ranked cluster was screened by the cytoHubba Cytoscape plugin, using Maximal Clique Centrality (MCC) parameter to detect hub genes. 42

| Reconstruction of miRNA-lncRNA-TF-hub gene co-expression network
The significantly differentially expressed miRNAs, lncRNAs and transcription factors (TF), which have co-expression correlation with hub genes, were identified. The co-expression network was visualized using Cytoscape version 3.7.2.

| Hub gene screening and validation
Among the candidate hub genes, the Kaplan-Meier survival curve and the receiver operating characteristic (ROC) curve were used to predict the potential ability of each gene to be independent predictors using the KM-plotter 43 (Figure 2).
An optimal soft-thresholding power is primarily required to con-  Figure 3A as Epithelial and Figure 3B as Mesenchymal clustering dendrogram.
The number of genes per module has been displayed in Table S1.
Mesenchymal and Epithelial status were used as independent variables to calculate the module's tendency to ovarian cisplatin resistance initiation and progression.

| Four non-preserved modules are identified in epithelial cell lines through network preservation analysis
We used the Epithelial data set as the reference and the  Table S2). The alteration of connectivity patterns in the non-preserved modules may be related to drug resistance status and recurrence of ovarian cancer. We hypothesized that non-preserved modules in Epithelial cell lines might highlight dysregulated pathways in the drug resistance compared to the sensitive network. Figure 3C shows the preservation statistics of Epithelial modules in the Mesenchymal network.

| Antiquewhite4, brown and darkmagenta modules covering the most number of DEGs among non-preserved modules
A total of 679 and 463 genes showed significant up/down-regulation in non-preserved modules, respectively ( Figure 4A). After discrete assessment of overlapped DEGs among four non-preserved modules ( Figure S3), Epithelial and Mesenchymal markers discovered by Miow et al. 7,47 were also captured to choose critical modules in cisplatin resistance status. Figure  and darkmagenta non-preserved modules, respectively.

| Enrichment analysis in antiquewhite4, brown and darkmagenta non-preserved modules implies cancer progression
GO function and KEGG pathway enrichment analyses were performed to assess the function of co-expressed genes in the antiquewhite4, brown and darkmagenta modules. As shown in Figure

| Candidate modules PPI network construction illustrates key genes in cancer and drug resistance
The antiquewhite4, brown and darkmagenta module genes were  resistance ( Figure 8A). Our results show that AUC for OAS1 was 0.9519 (p < 0.05). At the optimal cut-off value of 3.85, both sensitivity and specificity were 100%. Similar results were obtained for SLC38A9, OAS2, IFIT2, IFIT3, IFIH1, RSAD2 and XAF1 ( Figure 8A; Table 1). These results demonstrate that these hub genes possessed a high ability to discriminate between cisplatin resistance and sensitivity. To analyse the expression signature of candidate hub genes as prognostic biomarkers in response to platinum treatment, publicly available data and tools from GEO, EGA and TCGA databases were utilized. As shown in Figure 8B  were significantly related to the epithelial phenotype, which indi-  Our results show that the overexpression of this lncRNA is also associated with a lower survival rate in ovarian cancer. LINC00965

| DISCUSS ION
is a newly discovered lncRNA, and there is no data about it. Abba et al. reported the molecular effects of LINC00885 as a new oncogenic lncRNA associated with early breast cancer progression.
Overexpression of this lncRNA is closely associated with a lower survival rate and poor prognosis in breast and cervical cancer patients. 127,128 Lu et al. 129 demonstrated that high expression of PAX8-AS1 is associated with poor relapse-free survival in thyroid cancer patients. However, a recent study reported the tumour suppressor role of PAX8-AS1. 130 PAX8 has a role in transforming the ovarian cancer cell into mesenchymal phenotype. 131 We hypothesize that PAX8-AS1 may maintain epithelial characteristics in these cells by regulating PAX8.
In summary, we provide a Holistic biological interpretation of gene expression data derived from different ovarian cancer cell lines. WGCNA analysis identified 40 modules, 4 of which were significantly associated with an epithelial phenotype. We suggest that the brown module contains a potential prognostic biomarker of ovarian cancer progression and drug resistance based on the expression pattern, ROC curve and survival analysis. This module includes 10 hub/key genes co-expressed with different TFs, miRNAs, and lncRNAs. Pathway analysis implicates that this coexpression pattern may have a role in developing drug resistance in ovarian cancer. This study also provides several circRNAs for future in vitro and in vivo investigations of their molecular mechanisms regulating cisplatin resistance in ovarian cancer. Our findings add to our understanding of the processes and genes that underpin ovarian drug resistance through the EMT process. However, solid experimental proof is required to confirm our predictions.

ACK N OWLED G EM ENTS
This study was performed at the University of Isfahan (Isfahan, Iran) and was supported by the Graduate Studies Office.

CO N FLI C T O F I NTE R E S T
The authors declare that the research was conducted without any commercial or financial relationships that could be construed as a potential conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in the NCBI Gene Expression Omnibus (GEO) at (https:// www.ncbi.nlm.nih.gov/geo/), reference number (GSE47856, GSE149146).