Screening and prognostic value of potential biomarkers for ovarian cancer.

BACKGROUND: Ovarian cancer is a common gynecological malignant tumor that greatly threatens women's health, so we screened potential biomarkers of ovarian cancer and analyzed their prognostic value.
METHODS: The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) datasets were used to analyze the ovarian cancer-related genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to analyze the function of ovarian cancer-related genes. The survival-related genes were screened out through the least absolute shrinkage and selection operator (LASSO) method. Multivariate Cox regression model and stepwise regression analysis were performed to construct the risk model. The receiver operating characteristic (ROC) and the area under the ROC curve (AUC) were used to evaluate the prediction accuracy of risk score model. Finally, gene set enrichment analysis (GSEA) and immune cell infiltration analysis were performed to investigate the biological function and immune cell infiltration.
RESULTS: A total of 111 genes were found to have common effects on survival. These genes were mainly involved in metabolism, protein phosphorylation and immune-related signaling pathways. Seven risk genes (AP3D1, DCAF10, FBXO16, LRFN4, PTPN2, SAYSD1, ZNF426) were screened out. Among these genes, AP3D1 and LRFN4 are risk genes and DCAF10, FBXO16, PTPN2, SAYSD1, and ZNF426 are protective genes. These findings suggest that risk status may be an independent prognostic factor. The risk score had a high predictive value for the prognosis of ovarian cancer. In addition, GSEA revealed that the biological function of genes expressed in patients at a high risk was mostly related to immune-related function. The contents of CD4+ T cells, macrophages, myeloid dendritic cells (mDC) and neutrophils were high in samples at a high risk for ovarian cancer.
CONCLUSIONS: The abnormal expression of AP3D1, DCAF10, FBXO16, LRFN4, PTPN2, SAYSD1 and ZNF426 is highly related to the progression of ovarian cancer. These seven genes can be used as independent prognostic markers of ovarian cancer. This study not only adds evidence to the pathogenesis of ovarian cancer but also provides scientific basis for judging the prognosis of ovarian cancer. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Ovarian cancer is the most common cause of death among the gynecologic malignancies, and mainly occurs in women aged over 50 years. This disease is frequently diagnosed at the late stage because patients are often asymptomatic, which seriously affects prognosis (1). Approximately 90% of ovarian cancers are epithelial tumors, which are classified into five histotypes: serous, endometrioid, mucinous, clear cells and undifferentiated. Advanced serous cancer is the most common and fatal of the ovarian cancers. In North America, advanced serous cancer accounts for 70% of all ovarian cancers and 90% of deaths from ovarian cancer (2,3). Endometrial carcinoma accounts for ~10% of ovarian cancers, and is usually low stage and low grade. Undifferentiated serous ovarian cancer has a high mortality rate, and is often associated with metastasis and drug resistance (4). Mucinous ovarian cancer is a rare tumor, accounting for ~3% of all epithelial ovarian cancers (5,6). Each type of ovarian cancer has its own molecular phenotype and unique treatment, so it is challenging to diagnose and treat this disease. To control the progression of ovarian cancer, its pathogenesis must be understood, which will help to better screen the populations at high risk for ovarian cancer, and the identification of potentially changeable pathological factors provides an opportunity to reduce morbidity (7). At present, the molecular mechanism of ovarian cancer remains unclear, which hinders prognosis and the development of treatment strategies. Therefore, it is an urgent task to provide new biomarkers or biological targets for ovarian cancer.

We present the following article in accordance with the REMARK reporting checklist (available at https://dx.doi.org/10.21037/atm-21-2627).

Methods

Data collection

Ovarian cancer-related clinical information and gene expression information were downloaded from the TCGA. All samples were ovarian cancer. Only samples of severe ovarian cancer are kept in the GEO dataset. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Data processing

Trimmed mean of M values (TMM) was used to normalize the TCGA sequencing data (8), and Robust multi-array average (RMA) was used to process the GEO data (9). After removing low-expression samples, the following criteria were used: for sequencing data, reads per kb per million (rpkm) ≥1 in at least one-third of the samples; for chip data, normalized expression value ≥ median value of overall expression in at least one-third of the samples.

Survival analysis

To further screen the genes greatly associated with ovarian cancer survival, we first analyzed the TCGA and GEO datasets using univariable Cox proportional hazards regression analysis. The threshold value was set to P<0.05. The key genes in the two datasets were screened, and key genes co-expressed by the two datasets were selected. The genes with common effect on tumor survival in the two datasets [log2 (HR) values >0 or <0] were further selected. The LASSO method was used to set 300 calculations (10). The genes with frequency ≥150 were defined as potentially highly related to survival. Finally, the P values of the Cox regression model of these potential genes highly related to survival were analyzed and compared between the TCGA and GEO datasets.

Risk model construction

A multivariate Cox regression model was used to analyze the genes that were potentially highly related to survival. Stepwise regression analysis was used to select the smallest akaike information criteria (AIC) information statistics for gene deletion. The remaining genes were used for risk model construction. The risk formula was generated according to the Cox regression results. The risk coefficient of each patient was calculated according to the formula. The patients were classified according to the median risk value: low-risk group (low median risk value) and high-risk group (high median risk value).

GSEA analysis

The log2 difference folds (high vs. low) of the genes detected in the datasets were analyzed. The genes were sequenced according to log2 difference folds (from high to low). The database used for GSEA analysis was the KEGG dataset (data from https://www.gsea-msigdb.org/gsea/msigdb). The pathways with Padj <0.05 and absolute value of normalized enrichment score (NES) ≥1 were considered as enriched pathways (11). The first six of them (arranged by P value from low to high) were selected to create the distribution curve of the enrichment score.

Immune cell infiltration analysis

To analyze the differences in the proportion of immune cells between samples of high and low risk for ovarian cancer, immune cell infiltration analysis of the TCGA datasets was performed by the Timer (tumor immune estimation resource; https://cistrome.shinyapps.io/timer) (12) and Cibersort (https://cibersort.stanford.edu/) (13), respectively. The Wilcoxon test was used to analyze the differences in cell proportions between samples of high- and low-risk ovarian cancer (14).

Statistical analysis

To evaluate the predictive accuracy of the risk score model, the receiver operating characteristic (ROC) curve was plotted and the area under the ROC (AUC) was calculated. The effects of age, grade, stage, risk status and other potential prognostic factors in the TCGA and GEO datasets were analyzed using univariate and multivariate Cox regression analyses.

Results

Principal component analysis (PCA)

To investigate the feasibility of subsequent study, PCA was applied to the TCGA and GEO datasets (). There were no outliers in the PCA results, indicating that sample homogeneity was good.

Figure 1

Principal component analysis (PCA) on ovarian cancer samples. (A) Ovarian cancer samples in The Cancer Genome Atlas (TCGA), (B) ovarian cancer samples in Gene Expression Omnibus (GEO) datasets.

Functional analysis of co-expressed genes

Univariate Cox proportional hazards regression was used to analyze the ovarian cancer data in the TCGA and GEO datasets and screen out co-expressed genes. According to the log2(HR) value, 111 co-expressed genes with common effects on ovarian cancer survival were screened (). The function of these co-expressed genes was analyzed based on the GO and KEGG datasets. Based on the results of biological process (BP), a type of GO analysis, these 111 co-expressed genes mainly participated in metabolism, protein phosphorylation, and other functions. According to the KEGG results, these genes were involved in immune-related signaling pathways ().

Table 1

Construction of risk models of seven genes screened out using multivariate Cox regression model and stepwise regression analysis

Gene	coeff	HR	HR.95L	HR.95H	P value
AP3D1	2.7	15	1.4	156	0.023
DCAF10	−2.7	0.067	0.0050	0.91	0.042
FBXO16	−1.0	0.36	0.14	0.95	0.039
LRFN4	1.3	3.5	0.61	20	0.16
PTPN2	−2.9	0.053	0.0050	0.57	0.015
SAYSD1	−2.11	0.12	0.020	0.75	0.023
ZNF426	−1.8	0.17	0.029	1.0	0.052

Figure 2

Functional analysis of co-expressed genes. (A) According to biological process results, 111 co-expressed genes mainly participated in metabolism, protein phosphorylation, and other functions. (B) Based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) results, the 111 co-expressed genes were involved in immune-related signaling pathways.

Screening of ovarian cancer survival-related genes

To investigate the relationship between ovarian cancer-related genes and survival, the LASSO method was used to further analyze the 111 co-expressed genes from the TCGA and GEO datasets. A total of 42 and 18 genes with a frequency ≥150 that were highly related to ovarian cancer survival were selected from the TCGA and CEO datasets, respectively (). Interestingly, nine common genes were found (). Finally, the P value of the Cox proportional hazards regression of these nine genes was compared between the TCGA and GEO datasets (). The results showed that AP3D1, DCAF10, FBXO16, LRFN4, PTPN2, SAYSD1, SMU1, WAC.AS1 and ZNF426 were highly related to ovarian cancer survival. Subsequent study of ovarian cancer biomarkers should be performed based on these nine genes.

Figure 3

Screening of ovarian cancer survival-related genes. The 111 genes screened out in the The Cancer Genome Atlas (TCGA; A) and Gene Expression Omnibus (GEO; B) datasets were further analyzed using the LASSO method: 300 calculations were designated. The genes with frequency >150 were selected as potentially highly related to ovarian cancer survival. (C) Nine common genes were screened out using the LASSO method. (D) Comparison of P values of Cox proportional hazards regression of these nine genes between the TCGA and GEO datasets. Y-axis represents –log10(P value) value. The dotted line is the position where P=0.05. It can be seen that AP3D1, DCAF10, FBXO16, LRFN4, PTPN2, SAYSD1, SMU1, WAC.AS1 and ZNF426 are highly related to ovarian cancer survival.

Prognostic value of ovarian cancer markers

To investigate the prognostic value of ovarian cancer markers, the nine genes were analyzed using a multivariate Cox regression model, and seven risk genes for ovarian cancer were finally screened out using stepwise regression analysis (). According to the Cox regression results, the risk formula was constructed as below: 2.7090 * AP3D1-2.6960 * DCAF10-1.0232 * FBXO16 +1.2613 * LRFN4-2.9331 * PTPN2-2.1094 * SAYSD1-1.7618 * ZNF426. According to the median risk value, patients in the TCGA and GEO datasets were assigned to low-risk and high-risk groups, and the ovarian cancer survival curve was plotted (), which showed that patients at a low risk of ovarian cancer had good survival status. The heat map () shows the expression level of risk genes in the low-risk and high-risk groups from the TCGA dataset and the distribution of clinical pathological features. Accordingly, AP3D1 and LRFN4 are highly expressed in patients with low survival and are considered risk genes for ovarian cancer. DCAF10, FBXO16, PTPN2, SAYSD1, SMU1, WAC.AS1 and ZNF426 are considered protective genes. Univariate and multivariate Cox regression analyses were performed to analyze the effects of age, stage, grade, and risk status on ovarian cancer survival based on the TCGA and CEO datasets (). Both univariate and multivariate Cox regression analysis results revealed that risk status was related to survival. These findings demonstrate that risk status may be an independent prognostic factor. To analyze the predictive value of a risk score for ovarian cancer prognosis, a time-dependent ROC curve was plotted according to the risk score of seven genes screened out of the TCGA and CEO datasets (). In the TCGA dataset, the AUC for 1-, 3- and 5-year prognosis models was 0.6, 0.628 and 0.667, respectively, and 0.741, 0.9 and 0.955 in the GEO dataset. Therefore, in the two datasets, the predictive power of 1-, 3- and 5-year prognosis models was very high. Generally speaking, the AUC values were mostly >0.6, and therefore, the 1-, 3- and 5-year prognosis models had high predictive value.

Figure 4

Patients were assigned to a high-risk group or low-risk group according to median risk value, and time-dependent survival curves were plotted. (A) Analysis of survival of patients at high and low risk for ovarian cancer in The Cancer Genome Atlas dataset. (B) Analysis on survival of patients at high and low risk for ovarian cancer in the Gene Expression Omnibus dataset.

Figure 5

Heat map displaying the expression level of risk genes in the low-risk and high-risk groups from The Cancer Genome Atlas (TCGA) dataset and the distribution of clinical pathological features.

Table 2

Analysis of the effects of age, stage, and risk status on ovarian cancer survival based on The Cancer Genome Atlas (TCGA) dataset using univariate and multivariate Cox regression analyses

	HR	CI.low	CI.up	P value
Univariate Cox regression analysis
Age (>60 vs. ≤60 years)	1.4	1.1	1.8	0.019
Stage (III/IV vs. I/II)	2.1	0.95	4.8	0.067
Risk (low vs. high)	0.51	0.39	0.67	1.00E-06
Multivariable Cox regression analysis
Age (>60 vs. ≤60 years)	1.5	1.2	2.0	0.0015
Stage (III/IV vs. I/II)	2.3	1.0	5.1	0.049
Risk (low vs. high)	0.48	0.37	0.63	1.39E-07

Table 3

Analysis of the effects of age, stage, and risk status on ovarian cancer survival based on the Gene Expression Omnibus (GEO) dataset using univariate and multivariate Cox regression analyses

	HR	CI.low	CI.up	P value
Univariate Cox regression analysis
Stage (III/IV vs. II)	3.9	0.54	28	0.18
Grade (3/4 vs. 2)	2.4	1.1	5.3	0.036
Risk (low vs. high)	0.11	0.053	0.21	2.80E-10
Multivariable Cox regression analysis
Stage (III/IV vs. II)	2.4	0.32	18	0.40
Grade (3/4 vs. 2)	1.6	0.54	4.6	0.41
Risk (low vs. high)	0.12	0.059	0.25	5.84E-09

Figure 6

Analysis of the predictive value of a risk score for ovarian cancer prognosis. The time-dependent ROC curves of risk scores of seven genes screened out from The Cancer Genome Atlas (A) and Gene Expression Omnibus (B) datasets were plotted using 1-, 3-, and 5-year prognosis models.

Investigation of the signaling pathways involved in ovarian cancer using GSEA

The GSEA results revealed that in the samples from patients at high risk of ovarian cancer the immune-related pathways were upregulated and ribosomal pathway was downregulated (). The first six pathways (arranged by P values from low to high) were used to plot the distribution curves of the enrichment score (), which showed that the genes in “complement and coagulation cascades”, “ECM-receptor interaction”, “Hematopoietic cell lineage”, “Th1 and Th2 cell differentiation” and “Toll-like receptor signaling pathway” were clustered in the high-expression region, and the genes in the ribosomal pathway were clustered in the low-expression region. Generally speaking, the biological functions of the signaling pathways involved in ovarian cancer in high-risk patients were mostly related to the patient’s immune function.

Figure 7

Gene set enrichment analysis (GSEA) of the biological functions of the signaling pathways involved in ovarian cancer in high- and low-risk patients. The first 15 signaling pathways arranged by P value from low to high were selected to create histograms. Pos refers to that the genes in the signaling pathway that were upregulated in the samples from patients at a high risk for ovarian cancer.

Figure 8

Distribution curves of enrichment scores. The first six signaling pathways (arranged by P value from low to high) were selected to plot the distribution curves of enrichment scores. ECM, extracellular matrix.

Immune cell infiltration analysis using Timer and Cibersort

Immune cell infiltration analysis using Timer and Cibersort () revealed that high-risk status samples tended to have high abundance of most immune cells. B cells, mDC activated cells and NK(natural killer)-activated cells were less expressed in high-risk status samples, and CD4+ T cells, macrophages, mDC cells and neutrophils were highly expressed in high-risk status samples.

Figure 9

Analysis of the proportions of infiltrating immune cells in high- and low-risk samples using Timer. x-axis indicates high and low. y-axis represents log2 cell fraction score. Higher log2 cell fraction score indicates greater cell abundance. mDC, myeloid dendritic cells.

Figure 10

Analysis of the proportions of infiltrating immune cells in high- and low-risk samples using Cibersort. x-axis indicates high and low. y-axis represents log2 cell fraction score. Higher log2 cell fraction score indicates higher cell abundance. mDC, myeloid dendritic cells; NK, natural killer; Tregs, regulatory T cells.

Discussion

Ovarian cancer is one of the most common gynecological malignancies and carries a poor prognosis (15). The molecular mechanism underlying the occurrence and development of ovarian cancer remains unclear. In this study, we investigated the relationship between seven ovarian cancer-related genes and the occurrence and development of ovarian cancer. More importantly, we found that these ovarian cancer-related genes were closely related to the degree of malignancy and prognosis of ovarian cancer. Our findings add more evidence to existing knowledge of the occurrence and development of ovarian cancer.

To investigate the relationship between ovarian cancer-related genes and survival, we screened nine genes related to survival and further investigated them. Firstly, we used a multivariate Cox regression model and stepwise regression analysis to screen seven risk genes to construct the risk model, then univariate and multivariate Cox regression analyses determined that risk status may be an independent prognostic factor. In addition, we used ROC curves to analyze the predictive value of a risk score for ovarian cancer prognosis. We calculated the AUC value and found that the risk score had a high predictive value. We found that AP3D1, DCAF10, FBXO16, LRFN4, PTPN2, SAYSD1 and ZNF426 can be used as independent markers of ovarian cancer. Finally, we performed GSEA and immune cell infiltration analyses to analyze the biological functions of these genes and confirmed the important significance of these genes as new biomarkers of ovarian cancer. We found that AP3D1 and LRFN4 are risk genes, and their high expression inhibits the survival of ovarian cancer patients. AP3D1 and LRFN4 promote the occurrence and development of many cancers or diseases. In colorectal cancer, the loss of the optic nerve element leads to immune evasion and intrinsic immunotherapy resistance. Unexpectedly, the optic nerve element interacts with AP3D1 to prevent the sorting and degradation of palmitoylated IFNGR1 lysosomes, thereby maintaining the integrity of interferon gamma and MHC (major histocompatibility complex)-I signaling (16). Hermansky-Pudlak syndrome represents a group of immune dysfunction syndromes, and mutation of AP3D1 encoding the main subunit AP-3 (δ) is the cause of Hermansky-Pudlak syndrome type 10 (17). In addition, leucine-rich repeat and fibronectin type III domain-containing (LRFN) family proteins are considered to be neuron-specific proteins. LRFN4 is expressed in cells in various cancers and leukemia. LRFN4 signaling plays an important role in monocyte/macrophage migration (18).

However, it is reported that AP3D1 and LRFN4 also play a protective role in many diseases. For example, envelope glycoprotein 51 (gp51) is necessary for bovine leukosis virus to enter bovine B lymphocytes. Bovine adaptor protein 3 complex subunit δ-1 is considered to be a potential receptor. There is evidence that the N-terminal part of gp51 interacts with the AP3D1 receptor in vitro, and provides a reliable silicon interaction model (19). In colorectal cancer, LRFN4 expression is upregulated, which strongly correlates with clinical pathological features and prognosis. High expression of LRFN4 reduces the risk of death in colorectal cancer patients (20).

We found that DCAF10, FBXO16, PTPN2, SAYSD1 and ZNF426 are protective genes in ovarian cancer patients, and their high expression promotes survival. These genes inhibit the occurrence and development of many cancers and diseases. For example, FBXO16 functions as a tumor suppressor, and is a component of the skp1-cullin1-f-box protein complex that targets nuclear β-catenin and promotes proteasome degradation through 26S proteasome. Depletion of FBXO16 results in increased levels of β-catenin, which in turn promotes cancer cell invasion, tumor growth and epithelial-mesenchymal transition (21). In breast cancer, FBXO16 is considered to be a potential clinical target and prognostic biomarker for patients with different molecular types of breast cancer (22). In addition, it has been reported that increasing the number of cytotoxic Tim-3+/CD8+ T cells can promote effective anti-tumor immunity, and PTPN2 is an attractive tumor immunotherapy target in immune cells (23). In breast cancer, PTPN2 protein is lost in half of breast cancer patients. The role of PTPN2 is related to the subtype of breast cancer, and PTPN2 affects the prognosis and treatment response of breast cancer (24).

However, these genes can also promote the occurrence and development of some diseases. For example, mutations in the gene locus encoding protein tyrosine phosphatase non receptor type 2 (i.e., PTPN2) are associated with inflammatory diseases, including inflammatory bowel disease, rheumatoid arthritis, and type 1 diabetes (25). In thyroid cancer, inflammation or oxidative stress induces the upregulation of PTPN2, which promotes the progression of thyroid cancer (26).

Conclusions

We identified that among the ovarian cancer-related genes, AP3D1, DCAF10, FBXO16, LRFN4, PTPN2, SAYSD1, and ZNF426 were greatly associated with the progression of ovarian cancer. AP3D1 and LRFN4 were risk genes, and DCAF10, FBXO16, PTPN2, SAYSD1 and ZNF426 were protective genes, which could be used as independent prognostic markers of ovarian cancer. Moreover, there was severe inflammatory infiltration in the tumors of patients at a high risk of ovarian cancer. This study adds more evidence to confirming the pathogenesis of ovarian cancer, determining prognostic indicators of ovarian cancer, and providing potential targets for targeted therapy of ovarian cancer.

The article’s supplementary files as