Exploring prognostic genes in ovarian cancer stage-related coexpression network modules.

Identification of meaningful cluster modules of differential genes or representative biomarkers related to the stages of ovarian cancer (OC) is pivotal, which may help to detect mechanisms of OC progression and evaluate OC patients' prognosis.We downloaded gene expression data and the corresponding clinical information of OC patients from The Cancer Genome Atlas (TCGA) database, which included 379 ovarian cancer patients. Differentially expressed genes (DEGs) of OC patients between stages were picked out using R. There were 731 differential genes between ovarian cancer stage II and stage III (DEGs II-III) and 563 differential genes between ovarian cancer stage III and stage IV (DEGs III-IV), then we performed GO analysis and Kyoto Encyclopedia of Gene and Genome (KEGG) pathway analysis using Database for Annotation, Visualization and Integrated Discovery (DAVID). Moreover, CytoHubba was used to detect the top 20 hub genes in DEGs II-III and DEGs III-IV, followed Cytoscape with search tool for the retrieval of interacting genes (STRING) and MCODE plug-in was utilized to construct protein-protein interaction (PPI) modules of these genes. Three important coexpression modules of DEGs II-III and 3 more meaningful modules of DEGs III-IV were detected from PPI network using molecular complex detection (MCODE) tool. In addition, 5 hub genes in these stage-related DEGs modules with worse overall survival were selected, including COL3A1, COL1A1, COL1A2, KRAS, NRAS. This bioinformatics analysis demonstrated that stage-related prognostic DEGs, such as COL3A1, COL1A1, COL1A2, KRAS, and NRAS might play an unfavorable role in the development as well as metastasis of ovarian cancer. Furthermore, they need to be experimentally verified as a new biomarker to predict OC patient prognosis.

Introduction

Ovarian cancer (OC) is the most lethal gynecological cancer and the fifth most common cause of cancer-related death among women in the United States.[ Due to latent symptoms and lack of reliable early screening means, most OC patients are diagnosed at an advanced stage (stage III–IV; International Federation of Gynecology and Obstetrics, FIGO).[ Early detection and diagnosis of OC remain the main target for successful treatment. For advanced-stage OC are much more likely to have a poor prognosis, exploring gene expression characteristics related to OC stage is critical.

Currently, the only biomarker that is widely used in clinical practice is cancer antigen 125 (CA125).[ This high MW glycoprotein CA125 is elevated in 90% of patients with advanced stage disease. However, a number of false positive results could also occur since levels of CA125 could naturally be elevated during ovulation and may also be elevated due to a range of benign gynecologic causes such as fibroids, endometriosis, and pelvic inflammatory disease among others. CA125 can also be elevated in a variety of cancers other than ovarian such as pancreatic, lung, and breast cancer.[ In addition to CA125, other biomarkers are used routinely in medical practice and these include CA 19-9, CA 15-3, CA 72-4, and CEA. CA 19-9 has the advantage of high sensitivity for mucinous ovarian cancers that fail to express CA125.[ Serum levels of CA 19-9 are elevated in 68% to 83% of mucinous ovarian cancers but in only 28% to 29% of nonmucinous types. CA 15–3, CA 72–4, and CEA levels are found to be elevated in 50% to 56%, 63% to 71%, and 25% to 50% of patients with ovarian cancer.[ However, serial measurement of these tumor markers still plays a vital role in the management of patients with a CA125 negative tumor.[ Panels of biomarkers are thought to offer the potential for higher discriminatory power. Recent studies that constructed putative biomarker panels with samples from the prostate, lung, colorectal, and ovarian (PLCO) cancer trial found no improvement in diagnostic power in preclinical samples.[ Therefore, researchers still try hard to discover new biomarker to assist diagnose ovarian cancer earlier and can accurately predict the patients’ prognosis.

High-throughput sequencing is increasingly used and it has been used as a very significant tool for life sciences, such as cancer early diagnosis, cancer stage, and prognosis prediction.[ In this analysis, we downloaded data from TCGA database and used Edger R package to detect the stage-related differentially expressed genes (DEGs). Followed by, we selected the top 20 hub genes in the 2 groups of DEGs and established PPI network of the stage-related DEGs modules. Moreover, analysis of biological process (BP), molecular function (MF), cellular component (CC), and KEGG pathways of the DEGs and 6 meaningful coexpression modules significantly related to tumor stage were performed. Finally, overall survival (OS) analysis of these genes in the 6 modules was conducted using the Kaplan–Meier plotter online database (http://kmplot.com/analysis/). Then, COL3A1, COL1A1, COL1A2, KRAS, and NRAS were selected as the prognostic genes that could be used as a new potential biomarker for diagnosis and to predict the prognosis of OC patients.

Materials and methods

Microarray data

Gene expression data of OC patients were downloaded from TCGA database (The Cancer Genome Atlas, http://cancergenome.nih.gov/. The National Cancer Institute and National Human Genome Research Institute work with physicians who collect tissue for TCGA to gain approval with local Institutional Review Boards (IRBs). An IRB is a group of scientists, doctors, clergy, and consumers who review and approve the research proposal for every research project that involves human subjects. These boards ensure that the research is well designed, legal and ethical, and does not involve unnecessary risks to patients.). Meanwhile, the relevant clinical information was also obtained. We download 379 tumor samples of stage II to stage IV. The data is level 3 data which has been normalized.

Screening of stage-related DEGs and hub genes

We use R Language (Edger R package) to screen the differential genes between ovarian cancer stage II and stage III, and next obtained the DEGs between stage III and stage IV. The adjusted P value < .05 and FDR < .05 were set as the cut-off criterion. Then, 731DEGs II-III and 563DEGs III-IV were detected. CytoHubba is a tool that provides 11 topological analysis methods including Degree, Edge Percolated Component, Maximum Neighborhood Component, Density of Maximum Neighborhood Component, Maximal Clique Centrality, and 6 centralities (Bottleneck, EcCentricity, Closeness, Radiality, Betweenness, and Stress) based on shortest paths to detect the hub genes. Therefore, the top 20 DEGs II-III hub genes and the top 20 DEGs III-IV hub genes were screened by the cytoHubba according to the high degree of connectivity, then we selected these DEGs and hub genes between stages for the further analysis.

Functional and pathway enrichment analysis of DEGs

DAVID is the Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/) which aims to provide online bioinformatics tools for the functional interpretation of lists of genes or proteins.[ We could visualize the main pathways and biological processes, molecular functions, cellular components among those DEGs through DAVID.

PPI network, module analysis, and hub genes

Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org) is an online tool designed to evaluate the protein–protein interaction (PPI) information.[ To look for the potential interaction among those DEGs, we used STRING and input those DEGs into STRING, set confidence score ≥ 0.4, the maximum number of interactors = 0 as the cut-off criterion. Then we used the MCOED to construct 6 correlation modules. The pathway analysis of genes in each module was performed using DAVID. Also, some of the top 20 DEGs II-III hub genes and the top 20 DEGs III-IV hub genes which were screened by the cytoHubba according to the high degree of connectivity were found in these modules.[

Survival analysis of hub genes in stage-related modules

The Kaplan–Meier plotter (http://kmplot.com/analysis/) is capable of assessing the effect of 54,675 genes on survival using 10,461 cancer samples. These include 5143 breast, 1816 ovarian, 2437 lung, and 1065 gastric cancer patients with a mean follow-up of 69/ 40/ 49/ 33 months. The primary purpose of the tool is a meta-analysis-based biomarker assessment.[ We used this online tool to plot the OS of the hub genes in the stage-related modules.

Visualization of the prognostic genes expression level

The Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/index.html) is an interactive web server for analyzing the RNA expression sequencing data of 9736 normal and 8587 tumors samples from the GTEx and TCGA the projects, based on a standard processing pipeline.[ We displayed the expression level of prognostic genes among stages through the stage plot.

Results

Screening stage-related DEGs and the top 20 hub genes

There were 379 ovarian cancer samples obtained from TCGA in this study. We aimed to obtain the differential genes between stages and found the prognostic genes for ovarian cancer. The R analysis (Edger R package) was applied to detect the DEGs between stage II and stage III, stage III and stage IV, using adjusted P value < .05 and FDR < 0.05 as cut-off criteria. Finally, 731 DEGs II-III that include 328 upregulated and 403 dysregulated genes and 563DEGs III-IV that include 222 upregulated and 341 dysregulated genes were found. Meanwhile, there are 62 overlap genes between DEGs II-III and DEG III-IV. Moreover, the top 20 DEGs II-III hub genes and the top 20 DEGs III-IV hub genes were screened by the cytoHubba according to the high degree of connectivity (Table 1). We selected these DEGs and hub genes between stages for further analysis.

Table 1

Top 20 hub genes in differential genes between stage II and stage III as well as top 20 hub genes in differential genes between stage III and stage IV.

GO function and KEGG pathway enrichment analysis

For a more profound exploring of the screened DEGs, we used DAVID tools to perform GO function and KEGG pathway enrichment analysis. Above all, we mapped 731 DEGs II-III and 563 DEGs III-IV to DAVID web server and conducted GO analysis, the results demonstrated the enriched KEGG pathway of the DEGs II-III and DEGs III-IV. The DEGs II-III were most significantly enriched in protein digestion and absorption, ribosome, melanoma, TGF-beta signaling pathway, vascular smooth muscle contraction, signaling pathways regulating pluripotency of stem cells, focal adhesion, platelet activation, malaria, oxytocin signaling pathway. The DEGs III-IV were mainly enriched in protein digestion and absorption, maturity onset diabetes of the young, ABC transporters, neuroactive ligand-receptor interaction, adrenergic signaling in cardiomyocytes, bile secretion, gap junction, nicotine addiction, pancreatic secretion, melanogenesis. In addition, the BP of DEGs II-III includes regulation of cell proliferation, tube development, homeostatic process, blood vessel development, response to estrogen stimulus, and cell–cell signaling. For MF, these genes were enriched in growth factor activity, carbohydrate binding, structural molecule activity, heparin binding, and cytokine activity. GO CC analysis also revealed that the DEGs II-III were significantly enriched in the extracellular region, extracellular region part, proteinaceous extracellular matrix, ribonucleoprotein complex, and cytoplasmic membrane-bounded vesicle (Fig. 1A). Meanwhile, the BP of DEGs III-IV includes chemical synaptic transmission, cell–cell signaling, positive regulation of cell proliferation, transmembrane transport, cell adhesion. For MF, these genes were enriched in sequence-specific DNA binding, structural constituent of cytoskeleton, heparin binding, transporter activity, and receptor binding. At last GO CC analysis also showed that the DEGs III-IV were significantly enriched in the extracellular space, proteinaceous extracellular matrix, extracellular region, anchored component of membrane (Fig. 1B).

Figure 1

The BP, CC, MF of DEGs II––III (A), and DEGs III-IV (B). BP = biological process, CC = cellular component, MF = molecular function, DEGs II–III = differential genes between ovarian cancer stage II and stage III, DEGs III–IV = differential genes between ovarian cancer stage III and stage IV.

Construct stage-related EDGs modules from PPI co-expression network

We made the PPI network of these DEGs that based on the information in the STRING protein query from public databases. Aim to detect meaningful modules in this PPI coexpression network, we used MCODE plug-in. The modules of DEGs II-III and DEGs III-IV were conducted separately. KEGG pathway enrichment analysis demonstrated that these 3 DEGs II-III modules were mainly associated with RAS signaling pathway, MAPK signaling pathway, PI3K-Akt signaling pathway, and ECM-receptor interaction, focal adhesion (Fig. 2). Meanwhile, these 3 DEGs III-IV modules were mainly associated with neuroactive ligand-receptor interaction, pathways in cancer, retinol metabolism, WNT signaling pathway (Fig. 3).

Figure 2

Top 3 modules from the DEGs II–III protein-protein interaction network and the enriched pathways of DEGs II–III modules. DEGs II–III = differential genes between ovarian cancer stage II and stage III, DEGs III–IV = differential genes between ovarian cancer stage III and stage IV.

Figure 3

Top 3 modules from the DEGs III-IV protein-protein interaction network and the enriched pathways of DEGs III-IV modules. DEGs II–III = differential genes between ovarian cancer stage II and stage III, DEGs III–IV = differential genes between ovarian cancer stage III and stage IV.

The Kaplan–Meier plotter and expression level of prognostic genes

There are 10 hub genes in the 3 stage-related modules of DEGs II-III, they are COL1A1, COL1A2, COL3A1, EGF, IGF1, KRAS, MAPK1, NRAS, TIMP1, VEGFA. Besides, there are also 10 hub genes in the 3 stage-related DEGs III-IV modules, they are ADCY1, ALB, FGF8, GAL, KRAS, PENK, PF4, POMC, PYY, WNT3A. The prognostic information of the DEGs in these modules was obtained from Kaplan–Meier plotter software (http://kmplot.com/analysis/). It revealed that the expression of COL3A1 (HR 1.5 [1.3–1.72], P = 7.7 × 10–9) was associated with worse OS for ovarian cancer patients in these DEGs II-III module as well as COL1A2 (HR [1.39 1.21–1.6], P = 3.4 × 10–6), COL1A1 (HR [1.33 1.16–1.53], P = 5.8 × 10–5), KRAS (HR 1.18 [1.03–1.34], P = .014), NRAS (HR 1.49 [1.2–1.86] P = .00024). And in DEGs III-IV modules the expression of KARS (HR 1.18 [1.03–1.34], P = .014) was associated with worse OS (Fig. 4). Furthermore, we used GEPIA to detect the COL3A1, COL1A1, COL1A2, KRAS, NRAS expression level in different stages (Fig. 5).

Figure 4

Prognostic value of 5 genes (COL3A1, COL1A2, COL1A1, KRAS, NRAS) in ovarian cancer patients. HR = hazard ratio, CI = confidence interval.

Figure 5

GEPIA database displayed that COL3A1, COL1A1, COL1A2, KRAS, NRAS had a strong correlation with the progression of OC based on TCGA data. GEPIA = gene expression profiling interactive analysis, OC = ovarian cancer.

Discussion

In this study, a total of 731 DEGs II-III and 563 DEGs III-IV were screened, then we conducted 3 DEGs II-III coexpression modules and 3 DEGs III-IV coexpression modules. The DEGs II-III modules were mainly associated with RAS signaling pathway, MAPK signaling pathway, PI3K-Akt signaling pathway, ECM-receptor interaction, and focal adhesion. These pathways might associate with the ovarian cancer progression to the advanced stage. Moreover, COL3A1, COL1A1, COL1A2, KRAS, NRAS were selected as the prognostic genes because they are differential genes between the early and advanced cancer stages and might play an unfavorable role in the development of even the metastasis of ovarian cancer. The DEGs III-IV modules were mainly associated with the neuroactive ligand-receptor interaction, pathways in cancer, retinol metabolism, and WNT signaling pathway, and these pathways might contribute to metastasis of ovarian cancer. All these pathways might function in the progression and metastasis of ovarian cancer. For KRAS is the differential gene in both DEGs II-III module and DEGs III-IV module, then we select it as a prognostic gene for its worse overall survival of ovarian cancer patients.

Analysis of the 3 selected DEGs II-III modules from the PPI network showed that advanced ovarian cancer was associated with focal adhesion, PI3K-Akt signaling pathway and ECM- receptor interaction, RAS signaling pathway, MAPK signaling pathway. COL1A1, COL1A2, COL3A1were enriched in focal adhesion (bta04510) and ECM-receptor interaction (bta04512) KEGG pathway. It has been shown that ECM containing a large amount of collagen increases the invasiveness and the progression of tumors.[ It was reported that high expression of focal adhesion kinase activity was associated with elevated level of fibrosis and poor CD8+T cell infiltration. Focal adhesion kinase inhibition could substantially limit tumor progression and extend the survival time of cancer patients.[ In contrast to normal tissues where collagen is organized as thin, long wavy fibrils parallel to the epithelial boundary, collagen fibrils in tumor stroma are thicker and shorter.[ In epithelial ovarian cancer, collagen tracts that are perpendicular to the epithelial boundary have been observed.[ Thus far, the expression of COL1A1 and COL1A2 has been noted in gastric cancer and was positively correlated with the degree of invasion, metastasis, and advanced stages.[COL3A1 is frequently in association with type III collagen.[ In our study, COL3A1 is the differential stage-II-III hub gene and with the worst OS of OC. There was a study demonstrating that COL3A1 mRNA and protein was upregulated in CRC which was associated with clinicopathologic factors and poor survival, and COL3A1 was also increased in plasma of CRC patient. COL3A1 could be a potential diagnostic biomarker for colon cancer.[ A report indicates that COL3A1 gene had a prognostic implication in brain tumor.[ It also was associated with breast cancer development and progression.[ Moreover, it has been suggested that progressive ovarian carcinoma can induce expression of type III procollagen both in the tumor tissue and peritoneal cavity. In addition, in poorly differentiated serous ovarian carcinoma, the formation of type III procollagen may occur in the neoplastic cells.[ In serous ovarian carcinoma, production of type III procollagen has been found to be related to an increased degree of malignancy.[

KRAS, NRAS are predominantly expressed in various malignancies. RAS has 4 major isoforms, HRAS, NRAS and KRAS splice variants, KRas4A and KRas4B.[ RAS protein family members (KRAS4A, KRAS4B, HRAS, and NRAS) function as GDP-GTP-regulated on-off switches, which regulate cytoplasmic-nuclear signaling networks ruling diverse normal cellular processes. Constitutive activating mutations in RAS genes are found in up to 30% of human cancers,[ most frequently in KRAS (85%), then NRAS (15%), then HRAS (1%).[ And remarkably, the oncogenic RAS mutations and mutations in other components of RAS/MAPK signaling pathways seem to be mutually exclusive in most tumors, pointing out that deregulation of RAS-dependent signaling is an essential requirement for tumorigenesis.[ The mutations of 3 oncogenes KRAS and NRAS are very important in the development, spread, as well as in diagnostics and therapy of colorectal cancer and melanoma. In addition, KRAS and NRAS mutations are very well known to be mutated and play important role in pancreatic cancer, lung cancer, bladder cancer, and acute myeloid leukemia.[ Mucinous ovarian cancer tumors have prevalent KRAS mutations. KRAS and NRAS mutations have been shown to have transforming activity, so it showed that these mutations are rare but important drivers in High grade serous ovarian cystadenocarcinoma (HGS-OvCa).[NRAS mutations are much rarer than KRAS, meta-analysis of colorectal cancer patients carrying NRAS mutations had shorter progression-free survival (HR 2.30; 95% CI 1.30–4.07) and shorter overall survival (HR 1.85;95% CI 1.23–2.78).[ In RAS-driven cancers, KRas4B commanded most of the attention. Oncogenic mutants of KRas4B are abundant, particularly in adenocarcinomas, appearing in staggering frequencies, although KRas4A is currently also reappraised.[ RAS-driven cancer cells predominantly proliferate by the combined action of 2 pathways: MAPK and PI3Kα/Akt.[ Several studies indicated that numerous components of the phosphatidylinositol-3-kinase (PI3K)/AKT pathway were targeted by amplification, mutation, and translocation more frequently than any other pathway in cancer patients, leading to pathway activation.

In conclusion, our bioinformatics analysis detected stage-related DEGs and they might play a central role in the development and metastasis of ovarian cancer. In this study, a total of 1294 DEGs and 6 stage-related coexpression network modules were selected, and COL3A1, COL1A1, COL1A2, KRAS, and NRAS might be the prognostic genes of ovarian cancer. To get more accurate correlation results, we need to make a series of further verification experiments later to prove the results of this prediction.

Acknowledgments

The authors appreciate the researchers who provided their data and analyzed tools for this analysis, and thank Xu Wang for correcting the R code.

Author contributions

YY and LY conceived and designed the experiment. YY, LY, JJ, and LS analyzed data. LY wrote this manuscript. All authors reviewed and approved the final manuscript.

Conceptualization: Lili Yang, Ying Yue.

Data curation: Lili Yang, Liqun Sun.

Formal analysis: Lili Yang, Jili Jing, Liqun Sun.

Investigation: Jili Jing.

Methodology: Lili Yang, Jili Jing, Liqun Sun, Ying Yue.

Resources: Lili Yang.

Software: Lili Yang.

Supervision: Ying Yue.

Validation: Lili Yang, Ying Yue.

Visualization: Lili Yang, Jili Jing.

Writing – original draft: Lili Yang.

Writing – review & editing: Lili Yang, Jili Jing, Liqun Sun, Ying Yue.