Identification of hub genes and a prognostic model associated with ovarian cancer

doi:10.21203/rs.3.rs-16943/v1

Download PDF

Research article

Identification of hub genes and a prognostic model associated with ovarian cancer

https://doi.org/10.21203/rs.3.rs-16943/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Ovarian cancer (OC) is one of the most lethal gynecological cancers worldwide. The pathogenesis of the disease and outcomes prediction of ovarian cancer patients remain largely unclear. This study aimed to explore the key genes and biological pathways in ovarian carcinoma development, as well as construct a prognostic model to predict patients’ overall survival.

Results: We identified 164 up-regulated and 80 down-regulated differentially-expressed genes (DEGs) associated with ovarian cancer. GO term enrichment showed DEGs mainly correlated with spindle microtubes. For KEGG pathways, cell cycle was mostly enriched for the DEGs. The PPI network yielded 238 nodes and 1284 edges. Top three modules and 10 hub genes were further filtered and analyzed. Three candidiate drugs targeting for therapy were also selected. Thirteen OS-related genes were selected and an eight-mRNA model was present to stratify patients into high- and low-risk groups with significantly different survival.

Conclusions: The identified DEGs and biological pathways may provide new perspective on the pathogenesis and treatments of OC. The identified 8-mRNA signature has significant clinical implication for outcome prediction and tailored therapy guidance for OC patients.

Oncology

Ovarian cancer

bioinformatics

hub genes

biological function

prognostic model

Ovarian cancer (OC) is the most lethal malignant diseases in the female reproductive system, with over 200,000 new cases and 150,000 deaths happen each year worldwide[1]. Epithelial ovarian cancer accounts for 80–95% of ovarian malignancies, listed as the most common histological type. Since the ovaries locate in the deep pelvis, with mere symptoms emerging at the beginning of ovarian morbid change, the early detection for the malignancy is truly difficult. Hence, when ovarian cancer is detected, the patient is usually at an advanced stage, with invasion and metastasis accompanied[2]. For patients in the early stage, the 5-year survival rate can reach 85–90%, whereas, for advanced-stage patients, the number is mere ~ 20%[3]. Therefore, it is imperative to explore the molecular mechanisms of malignant biological behavior of ovarian cancer cells and to develop more reliable molecular markers for predicting recurrence and evaluating prognosis, further guiding clinicians for therapy.

At present, various high-throughput microarrays and next-generation sequence genomic datasets, which were deposited in the Gene Expression Omnibus (GEO)[4] and The Cancer Genome Atlas (TCGA) databases, have been widely analyzed for identifying differentially expressed genes (DEGs), which could serve as candidate diagnostic or prognostic markers, further effectively improving our understanding of the disease from genetic perspective. Whereas, since the existence of tissue or sample heterogeneity in each independent experiment, as well as the discrepancy of different data processing methods and technology platforms, the DEGs identified from a single-cohort study may generate false positives. Herein, the Robust Rank Aggregation (RRA) method, which analyzes the significant probability of all elements by a probabilistic model, is developed to identify statistically significant genes from multiple datasets and provide more accurate and valuable information for clinical use far beyond one gene list[5]. To date, a bunch of novel prognostic markers has been discovered to potentially improve the efficacy of diagnosis and prognosis of OC[6]. However, these identified markers were only effective for partial stages or grades and were difficult to apply widely. Hence, a prognostic model, which is based on signature gene expression level, with high discriminating power to effectively assist prognosis prediction for each patient is required in clinical practice.

In this study, we downloaded six primary microarray datasets from the GEO database which contained a total of 265 samples, with 201 ovarian cancer samples and 64 normal samples. The gene-set and relative clinical information on ovary tissues of OC patients and healthy females from TCGA and GETx portal were also downloaded. Integrated DEGs between cancerous and normal ovarian samples were screened using the “limma” R package and RRA method. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment of DEGs were performed for next-step functional analysis. The STRING and the Connectivity Map (CMap)online database was then used to analyze the association of DEGs and explore the molecular mechanisms as well as drugs involved in tumorigenesis. Through survival analysis, prognostic mRNAs were also selected. By performing Cox regression analysis, we identified an eight-mRNA signature model with the ability to predict the prognosis of OC patients and independent from clinical factors. Our study provides reliable molecular markers and prognostic models for early detection and outcome prediction, as well as effective drug targets for treating ovarian cancer.

The DEGs among six GEO microarray datasets

The distribution of gene expression from each microarray dataset was presented in volcano plots (Figure 1A-F), and the top 100 significantly up- and down-regulated genes were displayed in the heatmaps (Figure 2A-F). Through RRA analysis of 6 expression microarrays, we identified 605 DEGs, which consist of 301 up-regulated and 304 down-regulated genes, and displayed the top 20 dysregulated genes by "pheatmap" R package in Figure 3. Next, we analyzed the expression profiles of TCGA and GETx, getting 2253 dysregulated genes. Intriguingly, when combined these 2253 DEGs with the 605 differentially expressed genes from GEO datasets, we found that 244 genes were commonly dysregulated in these two databases, with 164 up-regulated (Figure 4A) and 80 down-regulated genes (Figure 4B).

GO term and KEGG pathway enrichment analysis of DEGs

To study the potential biological function of the 244 DEGs, we performed biological pathway analysis and identified significantly enriched pathways via Enrichr web tools. In GO term (Figure 5A), for the BP group, the DEGs were mostly enriched in "regulation of attachment of spindle microtubes to kinetochore", "cellular response to laminar fluid shear stress" and "microtubule cytoskeleton organization involved in mitosis". In MF group, the dysregulated genes were highly correlated to "microtubule-binding", "microtubule motor activity" and "tubulin-binding". As for CC group, the DEGs were closely related to "condensed nuclear chromosome kinetochore" and "mitotic spindle". KEGG pathway analysis showed 244 DEGs highly enriched in "cell cycle" and "Alanine, aspartate, and glutamate metabolism" (Figure 5B).

PPI network construction and modules analysis

Using the STRING database and Cytoscape software, a totally 244 DEGs were mapped into the PPI network, which included 238 nodes and 1284 edges (Figure 6A). The PPI enrichment p-value were < 1.0e-16. The top three key modules (Figure 6b-d) within PPI network were then selected (Module 1, MCODE score= 43.391; Module 2, MCODE score= 4.8; Module 3, MCODE score= 3.667) and the biological function of Module 1, which consisted of 47 nodes and 998 edges, was further analyzed. GO analysis indicated that Module1 was mainly associated with "regulation of attachment of spindle microtubules to kinetochore", "condensed nuclear chromosome kinetochore" and "microtubule motor activity". KEGG analysis showed that "cell cycle" and "Oocyte meiosis" were the most highly enriched pathways (Figure S1).

The screen of Hub genes and their characteristics

The top ten hub genes with the highest degree of connectivity were CDC45, CDK1, TOP2A, CDC20, CCNB1, CEP55, UBE2C, HMMR, FOXM1, and TPX2 (Figure 6E). The co-expression analysis results of the hub genes demonstrated that these genes actively interacted with each other (Supplementary 2B). Besides, we established the interaction network of ten hub genes with their related genes and explored the biological role (Figure S2A, C-E) of the network by FunRich. The gene alteration type and frequency, as well as the 50 most frequently altered neighbor genes were also exhibited (Figure 7). Gene alteration frequency of 10 hub genes among 606 TCGA ovarian cancer samples was beyond 50%, with most genes showed amplified and multiple altered (Figure 7A-B). The top 3 most frequent altered genes were FOXM1, CDC20 and CCNB1, with FOXM1 and CDC20 largely amplified while CCNB1 deep deleted. Through analysis of ovarian cancer patients’ gene-set from TCGA, we found that CCNB1, UBE2C, CDK1, CEP55 as well as FOXM1 expressed significantly higher in high-grade tumors and predicted worse outcomes since patients overexpressed above genes owned lower overall survival (OS) and disease-free survival (DFS) rates (Figure 8). The Oncomine database showed results from various studies were consistent to our finding (Figure S3). The HPA website also demonstrate that proteins translated by such 5 hub genes were overexpressed in ovarian cancer tissues (Figure S4). HMMR and TPX2 were also negatively correlated to patients' prognosis while no expression difference was observed in diverse tumor grades and CDC20 was positively associated with tumor grade but not correlated to patients' outcomes.

Related small molecule drugs screening

In total, 244 DEGs were analyzed in CMap to screen small molecule drugs, and the candidate molecules with top ten connectivity score are listed in Table 1. Five of these molecules showed a negative correlation and suggested potential in clinical applications. Among them, Trichostatin A, pyrvinium and 8-azaguanine showed a significantly negative correlation and the stuctures of such candidate molecule drugs was found in the PubChem database and shown in Figure S5.

Construction of prognostic model and evaluation of its predictive ability

Univariate Cox regression analysis revealed that 13 of 240 DEGs were significantly correlated to patients' OS in the training cohort. (Table 2) The 13 OS-related genes were listed as follows: CCND1, SYNE4, CCDC80, TMC4, MCC, FOXQ1, KRTCAP3, CXCR4, IL4I1, DEFB1, CSGALNACT1, KLHL14 and MCUR1. Through LASSO Cox regression, we narrowed the number of 13 prognosis-associated genes into twelve according to the minimum criteria (Figure 9). Next, based on the multivariate Cox model, 8 of 12 candidate mRNAs retained their prognostic significance and were finally selected as independent remarkable prognostic factors, which are TMC4, KLHL14, CXCR4, CCDC80, KRTCAP3, DEFB1, SYNE4 and FOXQ1 (Table 3). To predict patient’s outcomes, we developed an individual’s risk score model as follows: risk score = (0.006809 × expression value of TMC4)+ (0.021258× expression value of KLHL14)+(-0.00839× expression value of CXCR4) + (0.031459 × expression value of CCDC80)+(-0.00903 × expression value of KRTCAP3) + (-0.00156 × expression value of DEFB1) + (0.070689 × expression value of SYNE4) + (0.006726 × expression value of FOXQ1). On the basis of the median risk score, patients were divided into high-risk or low-risk groups. Kaplan-Meier curve analysis showed that the overall survival time of the lower-risk group was significantly longer than the high-risk group (P =1.147e-07) (Figure 10E). ROC curve analysis revealed the area under the ROC curve (AUC) of the prognostic model was 0.815 (Figure 10D). Meanwhile, the risk scores (Figure 10A) of ovarian cancer patients in the training group were ranked for displaying their distribution and the survival status (Figure 10B) was marked on the dot plot. The expression pattern of 8 prognostic mRNAs between high and low-risk groups was also shown in the heatmap (Figure 10C). Univariate and multivariate Cox regression analysis concerning the risk score and clinical factors showed that the prognostic model was able to serve as an independent prognostic indicator (Figure 11A-B). ROC curve analysis also showed that the AUC value of the model was 0.820, much significantly higher than the tumor stage (AUC= 0.542), grade (AUC= 0.574) and patients’ age (AUC= 0.701) (Figure 11C). Interestingly, when combined the risk score with clinical factors, the ROC curve of combination model was much higher than each alone(Figure 11D).

As for the testing cohort, we divided the group into 78 high-risk and 108 low-risk individuals based on the training cohort' cut-off risk score. The outcomes of low- and high-risk group patients of the testing cohort were also measured and the overall survival time of the high-risk group was significantly shorter than the lower-risk group (P =1.721e-02) (Figure 12E). The area under the ROC curve (AUC) of the prognostic model was 0.641 (Figure 12D). The risk scores distribution (Figure 12A) and survival status (Figure 12B) of ovarian cancer patients, as well as the 8-prognostic gene expression heatmap (Figure 12C) in the testing group were also present. Meanwhile, the independency of the prognostic model was confirmed in testing cohort since univariate and multivariate cox regression analysis showed the model correctly predicted high- or low-risk factor group patients' outcomes without relying on any clinical factors (Figure 13A-B). ROC curve analysis showed that the prognostic model exhibited better sensitivity and specificity when compared to tumor stage, grade and patients’ age for the AUC value of the model was much higher than latter (Figure 13C). In accordance with results from training cohort, the combination of risk score and clinical factors showed better OS prediction capability (Figure 13D).

The Prognostic Signature correlating to immune cells infiltration

Through TIMER web-tool, we analyzed the relative gene expression levels of six types of immune cells for each patient and found that genes concerning macrophage fraction were expressing significantly higher in the high-risk group (P<0.05) compared to the low-risk group in training cohort (Figure 14). Interestingly, same result was also observed in the testing cohort (Figure 15).

In this study, we used the RRA methods to jointly analyzed six GEO ovarian cancer microarrays which contain 201 ovarian cancer and 64 normal samples, identifying 605 DEGs and overlapped them with dysregulated genes of OC cohort from TCGA and GETx portal, finally getting 164 up-regulated and 80 down-regulated genes.

Functional analysis showed that 244 DEGs were significantly enriched in the cell division cycle, to be clear, in the process of the mitotic spindle. Spindle microtubules have been proved to play crucial role in physiological and pathological processes. As for cell division, only when all chromosomes linked to spindle microtubules through kinetochores and the spindle assembly checkpoint is satisfied, this process could step to anaphase[7]. Suraokar et.al found that the mitotic spindle assembly checkpoint and microtubule network were significantly altered in malignant pleural mesothelioma (MPM) while using epothilone-B, a non-taxane small molecule inhibitor targeting the microtubules, could greatly decrease the viability of MPM cell lines [8]. Rogalska et.al compared the anti-proliferative capacity of epothilone B with paclitaxel on ovarian cancer cell line SKOV-3, found that this effect of Epo B was greater than latter[9]. The researches above were consistent with our study that the mitotic spindle process was dysregulated in OC progression, playing important roles in ovarian cancer cell proliferation and tumor development.

PPI network construction of 244 DEGs includes 238 nodes and 1284 edges, among which we identified 3 key modules. Interestingly, the top1 module was also highly associated with spindle microtubules and chromosome kinetochore, confirming the role of cell cycle in OC pathogenesis. The top ten hub genes from the PPI network were also recognized, which are CDC45, CDK1, TOP2A, CDC20, CCNB1, CEP55, UBE2C, HMMR, FOXM1, and TPX2. Among them, CCNB1, UBE2C, CDK1, CEP55 as well as FOXM1 were found overexpressed in high-grade tumors and predicted worse outcomes. Besides, FOXM1, CDC20 and CCNB1were the most frequent altered genes. These genes were reported to closely associated with the BRCA1/2 mutation process of ovarian cancer. It has been reported that females with BRCA1 or BRCA2 mutations were much more susceptible to get ovarian cancer, accounting for the majority of the cohort[10]. Treszezamsky et.al found that BRCA1- and BRCA2-deficient cells are sensitive to Etoposide, which targeting topoisomerase II (TOP2A) and inducing DNA double-strand breaks [11]. High expression of CCNB1 was also observed in BRCA1-mutant cancer and induction of vinblastine targeting CCNB1 could significantly reduce tumor progression[12]. BRCA2 could interact with Filamin A actin-binding protein, further recruiting endosomal sorting complex required for transport (ESCRT)-associated proteins, Alix and Tsg101, and forming CEP55-Alix and CEP55-Tsg101 complexes at the midbody. The disruption of these processes by BRCA2 mutations results in increased cytokinetic defects, in part explain the instability of whole-chromosome in BRCA2-deficient ovarian cancer and propose potential therapeutic target of CEP55[13]. Olaparib, as a PARP inhibitor (PARP-i), has been widely used in BRCA1 or BRCA2 mutated ovarian cancer patients’ treatment. However, Fang et.al found that Olaparib-induced adaptive response could be disrupted by FOXM1 while inhibiting FOXM1 by Thiostrepton could significantly enhance sensitivity to PARP-i. It is noteworthy that other genes were also been reported in previous studies of cancer. Yang et.al reported that CDC45 activated by DNA J heat shock protein family (Hsp40) member A1 (DNAJA1) could be reversed by KNK437 in colorectal cancer. The joint treatment of KNK437 with 5-FU/L-OHP chemotherapy significantly reduced liver metastasis of CRC. Cyclin-dependent kinase 1 (CDK1), a key regulator for cell cycle, was overexpressed in paclitaxel-resistant ovarian cancer and predicted a poor overall survival[14], while miR-490-3P could reduce CDK1 expression, impeding ovarian cancer cell proliferation [15] and Alsterpaullone could effectively reverse the drug-resistant trend[16]. Through searching CMap, we found Trichostatin A, pyrvinium and 8-azaguanine negatively correlated to the genomic-wide changes of OC. Trichostatin A has been proved to enhance the apoptotic potential of Palladium nanoparticles and increased the therapeutic potential in cervical cancer[17]. Likewise, co-treatment with BEZ235 and Trichostatin A enhanced autophagic cell death via up-regulating LC3B-II and Beclin-1 expression, finally exerting anti-tumor activity in breast cancer[18]. Pyrvinium was found to inhibit cell autophagy and promote cancer cell death. The combination of pyrvinium with autophagy stimuli improves its toxicity against cancer cells[19]. 8-azaguanine has also been used for treatment of various carcinoma, sarcoma, osteogenic sarcoma, lymphosarcoma and melanoma[20]. Hence, except for cisplatin or PARPi in treating ovarian cancer, such small molecules may also reverse the malignant phenotypes of OC and serve as potential drugs for therapy.

By performing Univariate and multivariate Cox regression analysis, as well as LASSO regression methods for 244 DEGs, we developed an eight-mRNA model that could classify OC patients into the high- and low-risk group with significantly different overall survival. We explored the regulatory mechanism of eight-mRNAs in the signature by searching the published article, the majority of which were reported to be associated with tumorigenesis and tumor proliferation. DEFB1 is commonly considered as a single copy gene that encodes beta-defensin 1 (BD-1), a member of the host defense peptide group. In human cancers, BD-1 is proposed to inhibit cell growth and promote apoptosis, acting as a tumor suppressor[21, 22]. Zhang et.al found that forkhead box G1 (FOXG1) and miR-422a negatively regulated each other, forming a double-negative feedback loop to modulate the development and metastasis of hepatocellular carcinoma. KRTCAP3 was reported to be overexpressed in gastric cancer and human keratinocytes[23, 24] while KLHL14 participated in the development and metastasis of endometrial cancer[25]. Besides, KLHL14 was also found to mutate in primary central nervous system lymphoma (PCNSL), playing a role in CNS development[26]. Kyle et al proposed that SYNE4, as an outer nuclear membrane protein, could induce kinesin-mediated cell polarization[27]. The mutation of SYNE4 mediated the distinct disease phenotypes, acting as disease-causing behavior[28]. However, contrary to our observation, CXCR4 is observed to highly expressed in high-grade serous epithelial ovarian cancer which positively related to tumor dissemination and metastasis while CCDC80 is down-regulated in papillary thyroid carcinomas and considered as a tumor suppressor role[29, 30]. Treating with CXCR4 antagonists significantly inhibits tumor pro-invasive phenotype and knockdown of CCDC80 is susceptible to developed thyroid adenoma and ovarian cancer [31, 32]. Note that the role of TMC4 in ovarian cancer pathogenesis has not been studied. This may offer a new direction for TMC4 in ovarian cancer research. More recently, the prognostic value of mRNA-related signature has been reported in several studies[33, 34]. However, current traditional clinical risk factors and clinical models have limited success in predicting OC patients’ outcomes due to the molecular heterogeneity and false-positive rate. Our LASSO regression model results with independent validation suggested that the combination of eight mRNA has good robustness and reproducibility in predicting prognosis for OC patients independent from traditional clinical risk factors, with the area under the ROC curve (AUC) marked 0.815, significantly higher than tumor stage, grade, and patients’ age.

To investigate the biological function of various types of immune cells regarding ovarian cancer, we explored the relative gene of each immune cell type and found that macrophage part expressed expressively higher in high-risk group in both training and testing cohort, pointing out the oncogenic-role of macrophages in ovarian cancer development. Macrophage, as a type of immune-related cells, has already been considered to be closely associated with the malignant biological behavior of various cancers. M2 macrophage-like tumor-associated macrophages (TAMs) secreted EGF and then activated EGFR on tumor cells, further upregulating VEGF/VEGF-R signaling in surrounding tumor cells to finally mediate ovarian cancer cell proliferation and migration[35]. The exosomal miR-223 derived from macrophages under hypoxia condition reduced PTEN expression and led to increased PI3K/AKT signal activation, consequently mediated the drug resistance of EOC cells[36]. Hence, the potential therapeutic tools targeting macrophages may provide new perspective into ovarian cancer treatment.

There are some highlights of our study. First of all, from the RRA integrated approach, we jointly explored six ovarian cancer datasets in GEO databases and TCGA OC patients’ data matrix, finding some interesting niche factors and unique modules that were not seen earlier. Second, this study discovered a multitude of differentially expressed genes and hub genes between ovarian cancer and normal tissues, as well as the mutation condition of these genes. This information summarizes the genetic-level changes during the pathophysiological process of ovarian cancer and provides possible target molecules for further research. Third, the prognostic model in our study can effectively predict OC patients’ outcomes, which provide a new method to help gynecologists evaluate patients’ prognosis in clinical practices.

However, some aspects of our study required improvement. First, our research was completely based on public data analysis, additional experimental studies are needed to explore the detailed molecular mechanism regarding DEGs and pathways, as well as the eight-mRNA prognostic model. Second, candidate drugs targeting hub genes and immune-related cells are needed to explore and clinical trials are also needed to verify whether the hub genes can be targeted to truly exert therapeutic effects and whether the prognostic model effectively predicts patients’ outcomes for OC. That said, with the ever-increasing accessibility and volume of genomic data from clinical patients and the continued development of technologies and algorithms, the bioinformatic analysis will further promote the progress of accurate diagnoses and personalized treatment in ovarian cancer.

In summary, in this study, we identified 244 genes commonly dysregulated in ovarian cancer, with 164 up-regulated and 80 down-regulated genes. The most enriched biological pathways regarding DEGs were cell-cycle related processes. From the PPI network concerning all DEGs which comprised 238 nodes and 1284 edges, we seek out top3 hub modules and top10 hub modules. In addition, we found three molecular drugs may target ovarian cancer for therapy. We also performed cox and lasso regression to present and validate a robust prognostic model aggregating eight signature genes. Using this model, we could further distinguish patients with an elevated risk of mortality independent of other clinical factors, which may help us to improve our understanding of underlying mechanisms involved in OC and guide for diagnosing and prognosis prediction, as well as rational therapy in clinical practice.

Data collection

Through searching on the GEO Repository with “ovarian cancer”, we downloaded the gene expression profiles of GSE54388, GSE40595, GSE38666, GSE27651, GSE18520 and GSE14407 and the corresponding annotation files from the GPL570 [HG-U133_Plus_2] Affymetrix Human GenomeU133 Plus 2 Array platform. GSE54388 contains 22 ovarian tissue samples, with 6 normal ovarian surface epithelium and 16 tumor epithelial components from high grade serous ovarian cancer patients. GSE40595 includes 77 ovarian cancer-associated stroma and epithelium samples, which consist of 31 cancer-associated stroma samples and 32 epithelial tissues from high grade serous ovarian cancer patients, along with 8 stromal component and 6 ovarian epitheliums from the normal ovary. GSE38666 comprises 8 stroma and 8 matched ovariana epitheliums from 12 healthy females, as well as 7 cancer stroma and 7 matched cancer epitheliums from ovarian cancer patients. GSE27651 incorporates 6 normal ovarian surfaces epithelial and 8 serous borderline ovarian tumors, 13 low-grade serous ovarian carcinomas, and 22 high-grade serous ovarian carcinomas. GSE18520 covers 53 advanced stage, high-grade primary ovarian cancer specimens and 10 normal ovarian surface epithelium tissues. GSE14407 involves 12 healthy ovarian surface epithelial samples and paired serous ovarian cancer epithelium. Note that all samples from these GEO datasets are classified into the cancerous or normal part, to be clear, the normal stromal and surface epithelium is defined as normal ovarian tissues, whereas the borderline tumors as well as cancerous stromal and epithelial tissues are considered as malignancies. Besides, we also downloaded the gene expression data and relative clinical information of 374 ovarian cancer patients' samples and normal ovarian tissues from TCGA and GETx portal, respectively.

Screening for DEGs and Integration of microarray data

We used the “limma” R package[37] to 1) integrate the expression profiles from TCGA and GETx portal, 2) standardize the data from the integrated TCGA and GETx expression matrix as well as six GEO datasets, and 3) further screen the DEGs between ovarian cancerous and normal samples. The list of DEGs obtained from six GEO microarray datasets by limma analysis was further integrated by the " RRA" method to get prioritized commonly up- or down-regulated gene list. The final overlapped DEGs for subsequent biological function analysis were the combination of prioritized jointly dysregulated genes from six GEO microarrays and the results from TCGA and GETx database. The cut-off criteria were set as FDR < 0.05 and |log2fold change (FC)| > 1.

GO term and KEGG pathway enrichment analysis

GO classified the known genes into three main biological progress: Molecular Function (MF), Cellular Component (CC) and Biological Process (BP)[38]. KEGG provides researchers high-layer functions of the biological system from molecular-level information[39]. The Enrichr online tool （https://amp.pharm.mssm.edu/Enrichr/）allows for GO term annotation and KEGG pathway for a cluster of genes[40, 41]. We explored the biological functions of overlapped DEGs and hub modules from our PPI network using Enrichr website. P-value < 0.05 was considered as significant enrichment. Likewise, the functional biological pathways of the top 10 hub genes from PPI network were also analyzed by the FunRich tool (version: 3.0) [42]and the top 5 enriched pathways of up- and downregulated genes were displayed as bar charts, respectively. We set the p-value < 0.05 as statistically significant.

PPI network construction and analysis

PPI (Protein-protein interaction) networks display the relationships of various proteins according to their physical or biochemical properties. Search Tool for the Retrieval of Interacting Genes Database (STRING) is a database that encompasses the interaction information between known proteins and potentially interacted proteins[43]. In order to explore the correlations between the DEGs, we used the STRING database to construct a PPI network and visualize our results by Cytoscape software[44]. Confidence score > 0.4 was set as significant. Molecular Complex Detection (MCODE) was utilized to select hub modules of PPI networks in Cytoscape[45]. We set the degree cut-off= 2, node score cut-off= 0.2, k-score= 2, and max. Depth= 100 as the criterion. Then, the significant modules were performed by go and KEGG analysis. Top 10 genes were defined according to the high degree of connectivity in string network. [46]. The co-expression analysis of 10 hub genes was performed by STRING, either.

Validation of the hub genes

We downloaded the raw gene-set of ovarian cancer patients from TCGA to explore the expression differences of hub genes in low and high-grade tumor tissues of ovarian cancer and draw the survival plot using Kaplan Meier plotter web-tool[47]. The gene and protein expression level of grade-related hub genes were then confirmed by Oncomine and The Human Protein Atlas（HPA）database[48, 49]. The Meanwhile, we explored the genetic alteration information of the selected 10 hub genes in ovarian cancer patients by plug-in cBioPortal (cBio Cancer Genomics Portal) tool, which deposits the genomics profiles of various cancer types and provides analysis and visualization of the genomics datasets[50].

Identification of candidate small molecule drugs

The Connectivity Map (CMap) database was able to predict potential drugs which might reverse, or induce the biological state encoded in certain gene expression signatures in OC[51]. The 244 DEGs from our study were used to query the CMap database. The enrichment scores which represent the similarity were calculated (ranging from –1 to 1). The positive connectivity score means an inducing influence on the input signature, whereas drugs with negative connectivity score present reversion impact on the characteristic in human cell lines and are considered as candidate therapeutic molecules. After sorting all instances, the connectivity score of various instances was filtered by p value < 0.05. Next, we investigate the structures of these candidate molecular drugs from the Pubchem database (https://pubchem.ncbi.nlm.nih.gov/).

Establishment and evaluation of the prognostic model

The 374 ovarian cancer patients from the TCGA project were randomly classified into the training cohort (n = 188) and the testing cohort (n = 186). OS-related genes were determined by performing univariate Cox regression analysis in the training cohort with the "Survival" R package and further selected for the next-step screen. Least Absolute Shrinkage and Selection Operator (LASSO) is a parameter selection algorithm which shrinks all high-dimensional regression coefficients and generates the penalty regularization parameter λ via the cross-validation routine by “glmnet” R package. To select the optimal prognostic mRNAs, we adopted LASSO regression among the selected candidate genes and further perform multivariate Cox proportional hazards regression to evaluate their independent prognostic values. The risk-score model for predicting outcomes of ovarian cancer patients was the sum of each optimal prognostic mRNA expression level multiplying relative regression coefficient weight calculated from the multivariate Cox regression model.

All training cohort patients were classified into high- and low-risk groups according to the median risk score. The Kaplan–Meier curves of two diverse groups were plotted using “survfit” function and the receiver operating characteristic (ROC) curve was unfolded for OS prediction to estimate the sensitivity and specificity of the prognostic model. Cox multivariate analysis was also performed to exam whether the risk score was independent of the clinical characters, such as age, tumor stage, and grade. Next, we used the testing group to check the efficacy of the prognostic risk model. Each individual in the testing cohort was also categorized as high-risk or low-risk case by comparing the patient’s risk score with the cut-off value calculated from the training cohort. Kaplan–Meier curve analysis, time-dependent ROC analysis, and cox multivariate analysis were performed, either.

Searching tumor‐infiltrating immune cells associated with patients' prognostic signatures

The TIMER web-tool allows for systematical evaluations of the relationship between the six immune cell types in the tumor microenvironment which are B cell, CD4 T cell, CD8 T cell, neutrphil, macrophage as well as dendritic cell and clinical impact in various cancer types via a novel statistical method[52]. To further explore the prognostic signature, we used the TIMER online tool to search the most significant tumor‐infiltrating immune cells according to the TCGA ovarian cancer gene data. To be clear, the relative gene expression levels of six types of immune cells for each patient in high- and low-risk groups from training and testing cohort were measured.

OC: ovarian cancer; DEGs: differentially-expressed genes; RRA: Robust Rank Aggregation; GEO: Gene Expression Omnibus; TCGA: The Cancer Genome Atlas; GO: gene ontology ; MF:Molecular Function ; CC:Cellular Component ; BP:Biological Process; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: protein–protein interaction; STRING: Search Tool for the Retrieval of Interacting Genes Database; MCODE: Molecular Complex Detection; HPA: The Human Protein Atlas; CMap: the Connectivity Map; ROC: Receiver operating characteristic; AUC: area under the curve; LASSO: Least Absolute Shrinkage and Selection Operator ; OS:overall survival; DFS:disease-free survival.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Written informed consent for publication was obtained from all participants.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Competing interest

The authors declare that they have no competing interests.

Funding

This work was supported by the National Nature Science Foundation of China (81872119) and Jiangsu province medical innovation team (CXTDA2017008).

Author’s contribution

Jinhui Liu, Pinping Jiang and Huangyang Meng conceived and designed the experiments; Siyue Li, Sipei Nie and Sun Rui collected the data and perform the experiments; Yang Jing analysed the data; Pinping Jiang and Huangyang Meng wrote the paper.

Acknowledgements

Not applicable.

Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A: Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: a cancer journal for clinicians 2018, 68(6):394-424.
Urban N, Drescher C: Potential and limitations in early diagnosis of ovarian cancer. Advances in experimental medicine and biology 2008, 622:3-14.
Siegel RL, Miller KD, Jemal A: Cancer Statistics, 2017. CA: a cancer journal for clinicians 2017, 67(1):7-30.
Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, Marshall KA, Phillippy KH, Sherman PM, Holko M et al: NCBI GEO: archive for functional genomics data sets--update. Nucleic acids research 2013, 41(Database issue):D991-995.
Vosa U, Kolde R, Vilo J, Metspalu A, Annilo T: Comprehensive meta-analysis of microRNA expression using a robust rank aggregation approach. Methods in molecular biology (Clifton, NJ) 2014, 1182:361-373.
Mok SC, Bonome T, Vathipadiekal V, Bell A, Johnson ME, Wong KK, Park DC, Hao K, Yip DK, Donninger H et al: A gene signature predictive for outcome in advanced ovarian cancer identifies a survival factor: microfibril-associated glycoprotein 2. Cancer cell 2009, 16(6):521-532.
Hiruma Y, Sacristan C, Pachis ST, Adamopoulos A, Kuijt T, Ubbink M, von Castelmur E, Perrakis A, Kops GJ: CELL DIVISION CYCLE. Competition between MPS1 and microtubules at kinetochores regulates spindle checkpoint signaling. Science (New York, NY) 2015, 348(6240):1264-1267.
Suraokar MB, Nunez MI, Diao L, Chow CW, Kim D, Behrens C, Lin H, Lee S, Raso G, Moran C et al: Expression profiling stratifies mesothelioma tumors and signifies deregulation of spindle checkpoint pathway and microtubule network with therapeutic implications. Annals of oncology : official journal of the European Society for Medical Oncology 2014, 25(6):1184-1192.
Rogalska A, Marczak A, Gajek A, Szwed M, Sliwinska A, Drzewoski J, Jozwiak Z: Induction of apoptosis in human ovarian cancer cells by new anticancer compounds, epothilone A and B. Toxicol In Vitro 2013, 27(1):239-249.
Peto J, Collins N, Barfoot R, Seal S, Warren W, Rahman N, Easton DF, Evans C, Deacon J, Stratton MR: Prevalence of BRCA1 and BRCA2 gene mutations in patients with early-onset breast cancer. Journal of the National Cancer Institute 1999, 91(11):943-949.
Treszezamsky AD, Kachnic LA, Feng Z, Zhang J, Tokadjian C, Powell SN: BRCA1- and BRCA2-deficient cells are sensitive to etoposide-induced DNA double-strand breaks via topoisomerase II. Cancer Res 2007, 67(15):7078-7081.
Choi EK, Lim JA, Kim JK, Jang MS, Kim SE, Baek HJ, Park EJ, Kim TH, Deng CX, Wang RH et al: Cyclin B1 stability is increased by interaction with BRCA1, and its overexpression suppresses the progression of BRCA1-associated mammary tumors. Exp Mol Med 2018, 50(10):136.
Mondal G, Rowley M, Guidugli L, Wu J, Pankratz VS, Couch FJ: BRCA2 localization to the midbody by filamin A regulates cep55 signaling and completion of cytokinesis. Dev Cell 2012, 23(1):137-152.
Yang W, Cho H, Shin HY, Chung JY, Kang ES, Lee EJ, Kim JH: Accumulation of cytoplasmic Cdk1 is associated with cancer growth and survival rate in epithelial ovarian cancer. Oncotarget 2016, 7(31):49481-49497.
Chen S, Chen X, Xiu YL, Sun KX, Zhao Y: MicroRNA-490-3P targets CDK1 and inhibits ovarian epithelial carcinoma tumorigenesis and progression. Cancer Lett 2015, 362(1):122-130.
Bae T, Weon KY, Lee JW, Eum KH, Kim S, Choi JW: Restoration of paclitaxel resistance by CDK1 intervention in drug-resistant ovarian cancer. Carcinogenesis 2015, 36(12):1561-1571.
Zhang XF, Yan Q, Shen W, Gurunathan S: Trichostatin A Enhances the Apoptotic Potential of Palladium Nanoparticles in Human Cervical Cancer Cells. Int J Mol Sci 2016, 17(8).
Chen L, Jin T, Zhu K, Piao Y, Quan T, Quan C, Lin Z: PI3K/mTOR dual inhibitor BEZ235 and histone deacetylase inhibitor Trichostatin A synergistically exert anti-tumor activity in breast cancer. Oncotarget 2017, 8(7):11937-11949.
Deng L, Lei Y, Liu R, Li J, Yuan K, Li Y, Chen Y, Liu Y, Lu Y, Edwards CK, 3rd et al: Pyrvinium targets autophagy addiction to promote cancer cell death. Cell death & disease 2013, 4:e614.
Sugiura K, Hitchings GH, et al.: The effect of 8-azaguanine on the growth of carcinoma, sarcoma, osteogenic sarcoma, lymphosarcoma and melanoma in animals. Cancer Res 1950, 10(3):178-185.
Han Q, Wang R, Sun C, Jin X, Liu D, Zhao X, Wang L, Ji N, Li J, Zhou Y et al: Human beta-defensin-1 suppresses tumor migration and invasion and is an independent predictor for survival of oral squamous cell carcinoma patients. PLoS One 2014, 9(3):e91867.
Donald CD, Sun CQ, Lim SD, Macoska J, Cohen C, Amin MB, Young AN, Ganz TA, Marshall FF, Petros JA: Cancer-specific loss of beta-defensin 1 in renal and prostatic carcinomas. Laboratory investigation; a journal of technical methods and pathology 2003, 83(4):501-505.
Oh JH, Yang JO, Hahn Y, Kim MR, Byun SS, Jeon YJ, Kim JM, Song KS, Noh SM, Kim S et al: Transcriptome analysis of human gastric cancer. Mamm Genome 2005, 16(12):942-954.
Bonkobara M, Das A, Takao J, Cruz PD, Ariizumi K: Identification of novel genes for secreted and membrane-anchored proteins in human keratinocytes. The British journal of dermatology 2003, 148(4):654-664.
Wu X, Miao J, Jiang J, Liu F: Analysis of methylation profiling data of hyperplasia and primary and metastatic endometrial cancers. Eur J Obstet Gynecol Reprod Biol 2017, 217:161-166.
Vater I, Montesinos-Rongen M, Schlesner M, Haake A, Purschke F, Sprute R, Mettenmeyer N, Nazzal I, Nagel I, Gutwein J et al: The mutational pattern of primary lymphoma of the central nervous system determined by whole-exome sequencing. Leukemia 2015, 29(3):677-685.
Roux KJ, Crisp ML, Liu Q, Kim D, Kozlov S, Stewart CL, Burke B: Nesprin 4 is an outer nuclear membrane protein that can induce kinesin-mediated cell polarization. Proc Natl Acad Sci U S A 2009, 106(7):2194-2199.
Sahni N, Yi S, Taipale M, Fuxman Bass JI, Coulombe-Huntington J, Yang F, Peng J, Weile J, Karras GI, Wang Y et al: Widespread macromolecular interaction perturbations in human genetic disorders. Cell 2015, 161(3):647-660.
Gil M, Komorowski MP, Seshadri M, Rokita H, McGray AJ, Opyrchal M, Odunsi KO, Kozbor D: CXCL12/CXCR4 blockade by oncolytic virotherapy inhibits ovarian cancer growth by decreasing immunosuppression and targeting cancer-initiating cells. Journal of immunology (Baltimore, Md : 1950) 2014, 193(10):5327-5337.
Ferraro A, Schepis F, Leone V, Federico A, Borbone E, Pallante P, Berlingieri MT, Chiappetta G, Monaco M, Palmieri D et al: Tumor suppressor role of the CL2/DRO1/CCDC80 gene in thyroid carcinogenesis. J Clin Endocrinol Metab 2013, 98(7):2834-2843.
Righi E, Kashiwagi S, Yuan J, Santosuosso M, Leblanc P, Ingraham R, Forbes B, Edelblute B, Collette B, Xing D et al: CXCL12/CXCR4 blockade induces multimodal antitumor effects that prolong survival in an immunocompetent mouse model of ovarian cancer. Cancer Res 2011, 71(16):5522-5534.
Leone V, Ferraro A, Schepis F, Federico A, Sepe R, Arra C, Langella C, Palma G, De Lorenzo C, Troncone G et al: The cl2/dro1/ccdc80 null mice develop thyroid and ovarian neoplasias. Cancer Lett 2015, 357(2):535-541.
Xu Z, Zhou Y, Cao Y, Dinh TL, Wan J, Zhao M: Identification of candidate biomarkers and analysis of prognostic values in ovarian cancer by integrated bioinformatics analysis. Medical oncology (Northwood, London, England) 2016, 33(11):130.
Yang X, Zhu S, Li L, Zhang L, Xian S, Wang Y, Cheng Y: Identification of differentially expressed genes and signaling pathways in ovarian cancer by integrated bioinformatics analysis. OncoTargets and therapy 2018, 11:1457-1474.
Yin M, Li X, Tan S, Zhou HJ, Ji W, Bellone S, Xu X, Zhang H, Santin AD, Lou G et al: Tumor-associated macrophages drive spheroid formation during early transcoelomic metastasis of ovarian cancer. The Journal of clinical investigation 2016, 126(11):4157-4173.
Zhu X, Shen H, Yin X, Yang M, Wei H, Chen Q, Feng F, Liu Y, Xu W, Li Y: Macrophages derived exosomes deliver miR-223 to epithelial ovarian cancer cells to elicit a chemoresistant phenotype. Journal of experimental & clinical cancer research : CR 2019, 38(1):81.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK: limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic acids research 2015, 43(7):e47.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT et al: Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nature genetics 2000, 25(1):25-29.
Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K: KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic acids research 2017, 45(D1):D353-D361.
Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, Clark NR, Ma'ayan A: Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC bioinformatics 2013, 14:128.
Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, Koplev S, Jenkins SL, Jagodnik KM, Lachmann A et al: Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic acids research 2016, 44(W1):W90-97.
Pathan M, Keerthikumar S, Ang CS, Gangoda L, Quek CY, Williamson NA, Mouradov D, Sieber OM, Simpson RJ, Salim A et al: FunRich: An open access standalone functional enrichment and interaction network analysis tool. Proteomics 2015, 15(15):2597-2601.
Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P et al: STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic acids research 2019, 47(D1):D607-D613.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T: Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome research 2003, 13(11):2498-2504.
Bader GD, Hogue CW: An automated method for finding molecular complexes in large protein interaction networks. BMC bioinformatics 2003, 4:2.
Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY: cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC systems biology 2014, 8 Suppl 4:S11.
Nagy A, Lanczky A, Menyhart O, Gyorffy B: Validation of miRNA prognostic power in hepatocellular carcinoma using expression data of independent datasets. Scientific reports 2018, 8(1):9227.
Rhodes DR, Yu J, Shanker K, Deshpande N, Varambally R, Ghosh D, Barrette T, Pandey A, Chinnaiyan AM: ONCOMINE: a cancer microarray database and integrated data-mining platform. Neoplasia (New York, NY) 2004, 6(1):1-6.
Uhlen M, Zhang C, Lee S, Sjostedt E, Fagerberg L, Bidkhori G, Benfeitas R, Arif M, Liu Z, Edfors F et al: A pathology atlas of the human cancer transcriptome. Science (New York, NY) 2017, 357(6352).
Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, Jacobsen A, Byrne CJ, Heuer ML, Larsson E et al: The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer discovery 2012, 2(5):401-404.
Lamb J, Crawford ED, Peck D, Modell JW, Blat IC, Wrobel MJ, Lerner J, Brunet JP, Subramanian A, Ross KN et al: The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease. Science (New York, NY) 2006, 313(5795):1929-1935.
Li B, Severson E, Pignon JC, Zhao H, Li T, Novak J, Jiang P, Shen H, Aster JC, Rodig S et al: Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome biology 2016, 17(1):174.

Table 1 The top ten OC‐related small molecules with highly significant correlations in results of CMap analysis
Rank	Cmap name	Mean	N	Enrichment	P value
1	trichostatin A	-0.468	182	-0.448	0
2	8-azaguanine	-0.749	4	-0.913	0.0001
3	pyrvinium	-0.619	6	-0.772	0.00026
4	isoflupredone	0.7	3	0.937	0.00042
5	quinpirole	0.628	4	0.866	0.00044
6	vorinostat	-0.539	12	-0.558	0.00054
7	genistein	0.3	17	0.47	0.00064
8	antimycin A	-0.632	5	-0.788	0.00082
9	heptaminol	0.586	5	0.797	0.00084
10	midodrine	0.547	5	0.792	0.00086
Abbreviation: OC: ovarian cancer; CMap: connectivity map

Table2 Univariate cox regression identified 13 DEGs correlated to patients' OS
Gene ID	HR	HR.95L	HR.95H	P-value
CCND1	1.006326448	1.002566689	1.010100307	0.000959275
SYNE4	1.049235696	1.019307082	1.080043064	0.001133431
CCDC80	1.027423125	1.008550366	1.046649046	0.004235937
TMC4	1.008946557	1.002423207	1.015512357	0.007117928
MCC	1.158403166	1.035973925	1.29530084	0.0098773
FOXQ1	1.009835231	1.001811636	1.017923088	0.016186926
KRTCAP3	0.990364725	0.982252065	0.998544389	0.021051161
CXCR4	0.993964188	0.988729288	0.999226806	0.024636469
IL4I1	0.994202937	0.989120462	0.999311528	0.026193083
DEFB1	0.998104524	0.996327955	0.999884261	0.036860792
CSGALNACT1	1.430557852	1.014289405	2.017664543	0.041267552
KLHL14	1.020200309	1.000551278	1.040235213	0.043851588
MCUR1	0.953716572	0.910205203	0.999307955	0.046698901
Abbreviation: DEGs: differentially expressed genes; HR: Hazard Ratio

Table3 Multivariate cox regression selected 8 DEGs correlated to patients' OS
Gene ID	HR	HR.95L	HR.95H	P-value
TMC4	1.006831852	0.999904624	1.01380707	0.053250159
KLHL14	1.021485339	1.001857759	1.041497446	0.031756629
CXCR4	0.991645394	0.985647879	0.997679404	0.006716471
CCDC80	1.031959471	1.011590262	1.052738832	0.001982146
KRTCAP3	0.991005724	0.982638972	0.999443715	0.03674419
DEFB1	0.998439683	0.996606595	1.000276142	0.095816854
SYNE4	1.073246946	1.045846994	1.101364745	8.45E-08
FOXQ1	1.006748837	0.998099969	1.015472649	0.126528285
Abbreviation: DEGs: differentially expressed genes; HR: Hazard Ratio

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification of hub genes and a prognostic model associated with ovarian cancer

Status:

Version 1

Abstract

Figures

Background

Results

Discussion

Conclusion

Methods

Abbreviations

Declarations

References

Tables

Supplementary Files

Status:

Version 1