Identification of potential biomarkers in ovarian carcinoma and an evaluation of their prognostic value.

BACKGROUND: Ovarian cancer is one of the most common malignant tumors in female genital organs, and its incidence rate is high. However, the pathogenesis and prognostic markers of ovarian cancer are unclear. This study sought to screen potential markers of ovarian cancer and to explore their prognostic value.
METHODS: The Cancer Genome Atlas, Gene Expression Omnibus, Gene Ontology and Kyoto Encyclopedia of Genes and Genomes databases were used in this study. The least absolute shrinkage and selection operator (LASSO), multivariate Cox regression and stepwise regression analysis were chosen to screen genes and construct risk model. Gene Set Enrichment Analysis (GSEA) and an immune-infiltration analysis were performed.
RESULTS: One hundred thirty two co-expressed genes were found. They involved in metabolism, protein phosphorylation, mitochondria, and immune signaling pathways. Twelve genes significantly related to the survival of ovarian cancer were identified. Eight risk genes (i.e., CACNB1, FAM120B, HOXB2, MED19, PTPN2, SMU1, WAC.AS1, and BCL2L11) were further screened and used to construct the risk model. The risk status might be an independent prognostic factor of ovarian cancer, and most of the biological functions of genes expressed in high-risk ovarian cancer were related to synapse, adhesion, and immune-related functions. The clusters of CD4+ T cells and M2 macrophages were high in high-risk status samples.
CONCLUSIONS: In ovarian cancer, the abnormal expression of 8 genes, including CACNB1, FAM120B, HOXB2, MED19, PTPN2, SMU1, WAC.AS1, and BCL2L11, is closely related to ovarian cancer progression, and these genes can serve as independent prognosis markers of ovarian cancer. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Ovarian cancer is a common cancer in women worldwide. The 5-year survival rate of ovarian cancer is less than 45% (1). Ovarian cancer is the most common cause of death in gynecological malignancies (2). Because of its asymptomatic development, the disease is often diagnosed at an advanced, incurable stage (3). To treat ovarian cancer effectively, we need to analyze the cause of the disease. This will help to better predict high-risk groups that may be suitable for screening, and identify potential and changing etiologies, thereby providing an opportunity to intervene and thus reduce the incidence rate (1). To date, the underlying molecular mechanisms of ovarian cancer are unclear, which impedes the development of its diagnosis and treatment. Thus, finding new biomarkers or biological targets of ovarian cancer is an urgent task.

Currently, the mechanisms of ovarian cancer-related genes in the development and progression of ovarian cancer have not been widely explored. However, previous research has shown that many ovarian cancer-related genes are involved in the pathogenesis and development of other tumors. For example, in erythroblastic oncogene B Receptor Tyrosine Kinase 2 positive (ERBB2+) gastric cancer (GC), CACNB1 (a drug sensitive biomarker approved by the Food and Drug Administration) can be used to target ERBB2. All of these genes work through the myelocytomatosis oncogene (MYC) signaling pathway. The expression of CACNB1 is negatively correlated with Myc activation (4). In glioma, homeobox B2 (HOXB2) is an independent prognostic marker. HOXB2 is related to the invasion and proliferation of glioma cells in vitro, and has a certain effect on the invasion of glioma cells (5). We examined the regulatory role of the SP1/LINC00339/miR-378a-3p/MED19 axis in colorectal cancer, and provided new insights into the molecular mechanism of colorectal cancer (6). An in vivo clustered regularly interspaced short palindromic repeats (CRISPR) screening showed that protein tyrosine phosphatase, non-receptor type 2 (PTPN2) deletion in tumor cells enhances the effect of immunotherapy by enhancing the antigen presentation and growth inhibition mediated by interferon-γ. PTPN2 can be used as a target of tumor immunotherapy (7). Micro-ribonucleic acid-29b (miR-29b) can be used as a potential therapeutic molecule in prostate cancer. In prostate cancer cell line 3 (PC3) cells with miR-29b overexpression, the pro-apoptotic gene BCL2L11 (BIM) is significantly overexpressed. Additionally, in miR-29b treated xenograft tumors, BIM significantly induces apoptosis (8). In glioblastoma, anoctamin 6 (ANO6) promotes the proliferation and invasion of glioblastoma cells by promoting the nuclear translocation of extracellular signal-regulated kinase (ERK) and activating the ERK signaling pathway (9). In glioma cells, the downregulation of F-box protein 16 (FBXO16) induces β-catenin nuclear accumulation, which activates the wingless and int-1 (WNT) signal, which in turn induces a highly proliferative and malignant state for glioma cells (10). In GC cells, FK866 selectively kills GC cells with EMT gene expression characteristics by inhibiting nicotinamide phosphoribosyl transferase in nicotinic acid phosphoribosyl transferase (NAPRT) deficient cells. In many EMT subtypes of gastric tumors, the loss of NAPRT expression can usually be achieved by promoter methylation (11,12).

In this study, The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases were used to analyze eight ovarian cancer-related genes. The results of a bioinformatics analysis revealed that 8 ovarian cancer-related genes played crucial roles in the tumorigenesis, invasion, and metastasis of ovarian cancer. We also conducted a univariate analysis, multivariate regression analysis, and receiver operating characteristic (ROC) curve to explore clinical prognostic value. Gene Set Enrichment Analysis (GSEA) were used to analyze the corresponding biological functions of high- and low-risk patients, and the prognostic values of these markers in ovarian cancer were further analyzed. We present the following article in accordance with the REMARK reporting checklist (available at https://dx.doi.org/10.21037/atm-21-4606).

Methods

Data acquisition

The clinical information and gene expression information related to ovarian cancer were downloaded from TCGA. All samples were tumor samples. The GEO (GSE63885) data set only contained serous ovarian cancer samples. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Data processing

Data normalization processing: Tumor, node, metastasis data from TCGA was used. The normalization method was employed, and the GEO data were normalized by robust multi-array average (RMA) (13). Removal of low expression samples: sequencing data will be executed according to RPKM ≥1 in at least 1/3 of the samples; the normalized expression value of the microarray data in at least 1/3 samples was ≥ the median value of the overall expression.

Survival analysis

To further screen the genes that were significantly related to the survival of ovarian cancer, a univariate Cox proportional hazards regression analysis was conducted to analyze TCGA and GEO tumor data sets. The threshold value was set as P<0.05. 1,422 and 1,448 genes were screened in each data set, respectively. We took 132 genes co-expressed by TCGA and GEO tumor data sets and made 300 calculations using the least absolute shrinkage and selection operator (LASSO) method (14). Ultimately, 12 genes were identified, and the P values of the Cox regression analysis results for these 12 genes from the 2 data sets were compared. These 12 genes were found to be significantly correlated with survival.

Construction of the risk model

A multivariable Cox regression analysis was performed to analyze the 12 selected genes, and the Akaike information criterion (AIC) information statistics were used as the criterion for the stepwise regression analysis (15). By selecting the smallest AIC information statistic to delete genes, 8 genes were ultimately selected to construct the risk model and risk formula. The following risk formula was used:

GSEA

Next, the log2 fold changes of the genes detected in the ovarian cancer database were analyzed. The genes were sorted according to the multiple of log2 fold change (from large to small). To test a set of related genes in coordinated way, GSEA of KEGG was enforced by “gseKEGG” function of “clusterProfiler” R package. The pathway with P adj <0.05 and absolute value of Normalized Enrichment Scores (NES) ≥1 was considered as the enriched pathway (16).

Immune-infiltration analysis

To explore the immune microenvironment in tumor sites between the high- and low-risk status samples, TCGA data sets were analyzed by Timer (Tumor Immune Estimation Resource) for tumor-infiltrating immune cell signatures and Cibersort (Cell type Identification By Estimating Relative Subsets Of RNA Transcripts) for relative fraction of 22 types of tumor-infiltrating immune cells (17-19).

Statistical analysis

To assess the accuracy of the risk assessment model, a receiver operating characteristic (ROC) curve was drawn and the area under the curve (AUC) wad calculated. Both univariate Cox regression and multivariate Cox regression analysis were performed to examine the potential prognostic factors, such as age, grade, stage and risk status in TCGA and GEO data sets. A P value less than 0.05 was considered statistically significant.

Results

Functional analysis of co-expressed genes in the 2 data sets

There were 1,422 and 1,448 genes screened in TCGA and GEO tumor data sets, respectively, by a univariate Cox proportional hazards region analysis. The results showed that 132 genes were co-expressed (see ). Then Gene Ontology (GO) and KEGG functions of these 132 co-expressed genes were analyzed. A biological process (BP) analysis showed that these co-expressed genes were mainly involved in metabolism, protein phosphorylation, mitochondria, and other functions. The KEGG analysis showed that these genes were mainly involved in the immune pathway (see ).

Figure 1

Functional analysis of co-expressed genes. (A) A univariate Cox proportional hazards region analysis was performed on TCGA and GEO tumor data sets; 132 co-expressed genes were identified by setting a threshold value of P<0.05. (B) The BP results showed that these genes were mainly involved in metabolism, protein phosphorylation, mitochondria, and other functions. (C) The KEGG analysis showed that these genes are involved in immune-related pathways. TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Screening of genes related to the survival of ovarian cancer

To analyze the relationship between the ovarian cancer-related genes and the survival of ovarian cancer patients, 132 co-expressed genes in the 2 data sets were further screened by a LASSO analysis. 39 and 24 genes were selected as potential survival-related genes in TCGA and GEO data sets, respectively (see ). Interestingly, we found 12 common genes (see ), and compared the p-values from the Cox regression analysis of these 12 genes in the 2 data sets (see ). PTPN2 (protein tyrosine phosphatase nonreceptor 2), SMU1 (suppressor of Mec-8 and Unc-52 protein homolog), WAC.AS1 (WAC Antisense RNA 1), CACNB1 (beta 1 subunit of voltage-dependent calcium channel), FAM120B (Family With Sequence Similarity 120B), ANO6 (Anoctamin-6), FBXO16 (F-Box Protein 16), MED19 (Mediator Complex Subunit 19), NAPRT (Nicotinate Phosphoribosyltransferase), BCL2L11 (BCL2 Like 11), HOX (homeobox gene), and FBXO27 (F-Box Protein 27) were found to be significantly correlated with the survival of patients with ovarian cancer.

Figure 2

Screening of survival-related genes. In TCGA (A) and GEO (B) data sets, 132 genes were further screened using the LASSO method. Genes with frequency >150 were selected as potential survival-related genes. (C) 12 common genes were screened by the LASSO analysis. (D) The Cox regression P values of these 12 genes in the two groups were compared. The y-axis represents the –log10 (P value) value, and the dotted line indicates P=0.05. TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; LASSO, least absolute shrinkage and selection operator.

Prognostic value of markers in ovarian cancer

To study the prognosis value of the markers in ovarian cancer, we first performed a multivariate Cox regression model analysis on the above-mentioned 12 genes, and then used a stepwise regression analysis to screen out 8 risk genes (see ). Next, the patients were grouped into the high- and low-risk group according to the risk median of the TCGA and GEO data sets, and survival curves were drawn (see ). As shows, patients with a low risk had a good living condition. Next, a univariate analysis and a multivariate regression analysis were conducted to assess the age, stage, grade, and risk status in TCGA and GEO data sets (see ). The results for both the single-factor and multi-factor analyses showed that risk status was related to ovarian cancer survival. Thus, risk status can serve as an independent prognostic factor for ovarian cancer. To assess the accuracy of the risk-score model, a time-dependent ROC curve based on the risk score of the 8 genes was drawn and the AUC was calculated (see ). In TCGA data set, we found that the AUCs of 1-, 3-, and 5-year ROC curves were 0.681, 0.667, and 0.695, respectively. In the GEO data set, we found that the AUCs of the 1-, 3-, and 5-year ROC curves of the prognostic model were 0.587, 0.625, and 0.66, respectively. In TCGA and GEO data sets, the predictive power of the 1-, 3- and 5-year risk-score models was very high.

Table 1

8 genes were selected to construct a risk model by a multivariate Cox regression analysis and a stepwise regression analysis

Id	Coef	HR	HR.95L	HR.95H	P value
CACNB1	1.4	3.9	1.2	13	0.023
FAM120B	2.8	16	2.2	118	0.0060
HOXB2	0.84	2.3	1.0	5.1	0.040
MED19	–2.4	0.094	0.010	0.89	0.034
PTPN2	–4.2	0.016	0.0013	0.19	0.0012
SMU1	–3.0	0.049	0.0032	0.75	0.030
WAC.AS1	–2.5	0.080	0.011	0.58	0.012
BCL2L11	2.4	11	1.1	103	0.040

Figure 3

Patients were divided into the high- and low-risk groups according to the median risk, and survival curves were drawn. (A) A survival analysis of the high- and low-risk groups in TCGA data set. (B) A survival analysis of the high- and low-risk groups in the GEO data set. TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus.

Table 2

Survival analysis of age, stage, and risk status in TCGA data set by a univariate Cox regression analysis

Clinicopathological characteristics	Beta	HR	Lower .95	Upper .95	P value
Age (>60)	0.31	1.4	1.1	1.8	0.019
Stage (3+4)	0.76	2.1	0.95	4.8	0.067
Risk_status (low)	–0.68	0.5	0.39	0.66	5.30E-07

TCGA, The Cancer Genome Atlas.

Table 3

Survival analysis of age, stage, and risk status in TCGA data set by a multivariate Cox regression analysis

Clinicopathological characteristics	P value	HR	Lower .95	Upper .95
Age (>60)	0.018	1.4	1.0	1.8
Stage (3+4)	0.058	2.2	0.98	5.0
Risk status (low)	5.92E-07	0.50	0.39	0.66

TCGA, The Cancer Genome Atlas.

Table 4

Survival analysis of grade, stage, and risk status in the GEO data set by a univariate Cox regression analysis

Clinicopathological characteristics	Beta	HR	CI.low	CI.up	P value
Grade	0.86	2.4	1.1	5.3	0.036
Stage (3+4)	1.4	3.9	0.54	28	0.18
Risk status (low)	–0.71	0.49	0.29	0.82	0.0069

GEO, Gene Expression Omnibus.

Table 5

Survival analysis of grade, stage, and risk status in the GEO data set by a multivariate Cox regression analysis

Clinicopathological characteristics	P value	HR	CI.low	CI.up
Grade	0.068	3.18	0.924	11
Stage (3+4)	0.35	2.6	0.35	19
Risk status (low)	0.028	0.56	0.33	0.94

GEO, Gene Expression Omnibus.

Figure 4

ROC curve analysis of risk score in predicting the prognosis of lung cancer. TCGA (A) and GEO (B) data sets were used for risk score based on 8 genes, and a time-dependent ROC analysis was performed. ROC curves of 1, 3, and 5 years. ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus.

GSEA was conducted to explore the signaling pathways

In relation to the biological functions of genes expressed in ovarian cancer in high- and low-risk patients, the results of the GESA showed that synapses, adhesion junctions, and immune-related pathways were upregulated in high-risk samples (see ). Additionally, the first 6 pathways (sorted by P value from small to large) were chosen to draw the distribution curves of the enriched scores (see ). As shows, the genes involved in phosphatidylinositol-3-kinase and protein kinase B (PI3K/Akt), mitogen-activated protein kinase (MAPK), cyclic adenosine monophosphate (cAMP), calcium signaling pathway, focal adhesion, and neuroactive ligand-receptor interaction pathways tend to gather in the high expression region. Thus, the biological functions of genes expressed in ovarian cancer from high-risk patients are mostly involved in tumor proliferation and tumor immune.

Figure 5

A GSEA was conducted to analyze the biological functions of genes expressed in high- and low-risk patients. The enriched pathways were screened and sorted from small to large according to the P value. The top 15 pathways are shown in a histogram. POS means that genes in this pathway tend to be upregulated in high-risk samples. NES, Normalized Enrichment Scores; GSEA, Gene Set Enrichment Analysis.

Figure 6

Distribution curves of enriched scores. The first 6 (sorted by P value from small to large) were selected to make the distribution curves of the enriched scores.

Immunocyte infiltration was analyzed by Timer and Cibersort

To assess the proportion of immune cells between ovarian cancer from high- and low-risk patients, Timer and Cibersort were used to analyze the infiltration of immune cells (see ). Combining the results of the Timer and Cibersort analysis, it was clear that the high-risk ovarian cancer patients tended to have a high abundance of immune cells. The expression of B cells was low in high-risk ovarian cancer patient, while the expression of CD4+ T cells and M2 macrophages was high.

Figure 7

The proportion of immune cell infiltration in high- and low-risk status samples analyzed by Timer. The horizontal axis represents high- and low-risk status. The vertical axis represents the log2 cell fraction score value. The higher the value, the higher the abundance of cells.

Figure 8

The proportion of immune cell infiltration in high- and low-risk status samples analyzed by Cibersort. The horizontal axis represents high- and low-risk status. The vertical axis represents the log2 cell fraction score value. The higher the value, the greater the abundance of cells.

Discussion

Ovarian cancer is one of the most common malignant tumors and has a poor prognosis (20). Studies have shown that an increase of tumor infiltrating lymphocytes in ovarian cancer is associated with improved survival rates, while an increase of immunosuppressive regulatory T cells is associated with a poor prognosis (21). In addition, the upregulation or downregulation of genes is related to the oncogenesis and development of ovarian cancer. For example, RAD51C and RAD51D are genes that are moderately susceptible to ovarian cancer. They have a certain level of risk for invasive epithelial ovarian cancer, and may be used for routine clinical gene detection together with BRCA1 and BRCA2 (22). circPLEKHM3 inhibits ovarian cancer by targeting the miR-9/BRCA1/DNAJB6/KLF4/AKT1axis, and thus might be developed into a therapeutic target for ovarian cancer. The new strategy of Paclitaxel combined with Akt inhibitor MK-2206 in the treatment of ovarian cancer warrants further study, especially for ovarian cancer patients with a loss of circPLEKHM3 expression (23). YTHDF1 can enhance the translation of eIF3c in an m6A dependent manner, and promote the oncogenesis, invasion, and metastasis of ovarian cancer. The upregulation of YTHDF1 is related to the poor prognosis of ovarian cancer patients (24). However, due to the variety of pathological types of ovarian cancer, the pathogenesis and etiology remain unclear. There is no ideal tumor marker with high specificity and sensitivity for the early clinical diagnosis of ovarian cancer (25).

In this study, we investigated the relationship between 8 ovarian cancer-related genes and the oncogenesis and development of ovarian cancer. The results showed that the 8 ovarian cancer-related genes are closely associated with the malignant degree and prognosis of ovarian cancer. Our study further examined the mechanism of ovarian cancer development. To investigate the relationship between the ovarian cancer-related genes and the survival rate of ovarian cancer patients, 12 genes associated with ovarian cancer survival were identified. Based on these 12 genes, a follow-up study was conducted. First, we selected 8 risk genes to construct the risk model by conducting a multivariate Cox regression analysis and a stepwise regression analysis. The results showed that risk status may be an independent prognostic factor. In addition, ROC curves were used to analyze the accuracy of the risk score for the prognosis of ovarian cancer. The AUCs were calculated, and we found that the predictive value was still very high. We believe that CACNB1, FAM120B, HOXB2, MED19, PTPN2, SMU1, WAC.AS1, and BCL2L11 can serve as independent markers for ovarian cancer. Finally, GSEA and immune invasion were chosen to analyze the biological functions and explore these biological functions’ significance in patients with ovarian cancer.

CACNB1 is a calcium channel protein. In colon cancer, CACNB1 can be used to divide patients with colon cancer into high- and low-risk groups. Its high expression indicates a poor prognosis for colon cancer (26). HOXB2 is the top gene regulated by expression quantitative trait loci, which has many polymorphisms related to its messenger ribonucleic acid (mRNA) expression level in lung tissue. The HOXB gene mutation may increase a patient’s susceptibility to ovarian cancer by regulating the expression of HOXB in the lung (27). In bladder cancer, the overexpression of HOXB2 in RT112 cells induces the growth and adhesion of bladder cancer cells (28). In colorectal cancer, the overexpression of miR-4324 inhibits the proliferation, migration, and invasion of colorectal cancer cells. Conversely, the overexpression of HOXB2 promotes the malignant behavior of colorectal cancer cells. miR-4324 regulates the proliferation, migration, and invasion of colorectal cancer cells by targeting HOXB2 (29). In addition, in patients with Wilms’ tumors, high levels of HOXB2 and FOXC1 indicate a higher risk of advanced stage and lymph node metastasis, and a worse prognosis (30).

MED19 is a member of the mediator complex, which forms a bridge between the activator of transcription and RNA polymerase II (31). In breast cancer, MED19 was significantly increased, and was significantly correlated with larger tumors, higher malignancy, and a poor prognosis. MED19 promotes the proliferation, epithelial-mesenchymal transition, invasion, and migration of breast cancer cells (32). In cervical cancer, miR-4778-3p reduces the viability and the ability of proliferation and migration of cervical cancer cells, and may regulate the radio sensitivity of cervical cancer by targeting NR2C2 and MED19 (33). In GC, the downregulation of MED19 significantly reduces the ability of cell proliferation and colony formation, inducing cell cycle arrest in the G1 phase. In conclusion, MED19 has the function of promoting cell growth and may be an effective target for the treatment of malignant GC (31,34). In bladder cancer, the expression of MED19 is upregulated in tissues, and lentivirus-mediated MED19 inhibition via short hairpin RNA may be an effective method for the treatment of bladder cancer (35).

The variation of the PTPN2 gene locus is associated with inflammatory diseases, such as inflammatory bowel disease, rheumatoid arthritis and type 1 diabetes (36). An increased number of cytotoxic Tim-3+CD8+ T cells effectively improve anti-tumor immunity; thus, PTPN2 is an attractive target of tumor immunotherapy (37). In breast cancer, PTPN2 is not expressed in half of breast cancer patients. These findings show that PTPN2 plays a crucial role in the different subtypes of breast cancer, and affects the prognosis and treatment response of breast cancer (38). In addition, in thyroid cancer, the inflammatory response or oxidative stress induces the upregulation of PTPN2, leading to the progression of thyroid cancer (39).

BCL2L11 belongs to the BCL-2 family, and is the central regulatory factor of the endogenous apoptosis cascade response, which mediates apoptosis. In GC, miR-24 regulates the expression of BCL2L11 by directly binding to the 3’untranslated region (3’UTR) of mRNA, thereby inhibiting apoptosis and promoting cell growth and migration (40). In ovarian cancer, the overexpression of HAND2-AS1 or BCL2L11 or the downregulation of miR-340-5p decreases the invasion and migration of ovarian cancer cells, decreases cell proliferation, and increases apoptosis (41). In childhood acute lymphoblastic leukemia, BCL2L11 gene polymorphism encoding proapoptotic Bim protein affects the overall survival of childhood acute lymphoblastic leukemia patients treated with corticosteroids in a dose-dependent manner (42).

In this study, we analyzed the relationship between the expression of 8 ovarian cancer-related genes and the tumorigenesis, development and prognosis of ovarian cancer. The abnormal expression of CACNB1, FAM120B, HOXB2, MED19, PTPN2, SMU1, WAC.AS1, and BCL2L11 was significantly associated with the progression of ovarian cancer. These 8 genes can serve as independent markers to evaluate the prognosis of ovarian cancer, and indicate that there is severe inflammatory infiltration in the tumors of high-risk patients. This study not only extends our knowledge of the pathogenesis and prognostic indicators of ovarian cancer, but also provides potential targets for targeted therapy in ovarian cancer. However, the study still has some limitations, that is, these conclusions have not been further verified by experiments, which must be paid attention to in our later research.

The article’s supplementary files as