Evaluation of FGFR1 as a diagnostic biomarker for ovarian cancer using TCGA and GEO datasets

TCGA data description

To obtain a dataset based on studies listed in The Cancer Gene Atlas (TCGA) database, patients’ gene expression and clinical information were downloaded from the publicly available TCGA Data Portal at https://tcga-data.nci.nih.gov/tcga/. For the TCGA cohort, gene expression profiles were studied in 303 female patients with histologically confirmed ovarian cancer who had undergone prior surgical resection and received no pretreatment. Gene expression profiles were measured experimentally using an Illumina HiSeq2000 RNA Sequencing instrument at the University of North Carolina TCGA genome characterization center. Level 3 data was downloaded from the TCGA data coordination center. This dataset shows gene-level transcription estimates, as in log2(x+1) transformed RSEM normalized counts. Genes were mapped onto the human genome coordinates using HUGO probeMap. Reference to method description from University of North Carolina TCGA genome characterization center. All analyses were performed using FDR (False Discovery Rate) corrected q-values adjusted to 0.05.

GEO dataset selection

We obtained OC microarray profiles from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). The following keywords were used to query the GEO database: (ovary OR ovarium OR oophoron OR ootheca OR germarium) AND (cancer OR carcinoma OR tumor OR tumour OR neoplas* OR malignan*). All the functional genomics data of FGFR1 were requested and assembled from the GEO Database (http://www.ncbi.nlm.nih.gov/geo/) with the closing date of 10 February 2019. The inclusion and exclusion criteria of this study were as follows: (1) only human tissue studies were included, and samples from cell lines and animal models were excluded; (2) the expression data of the experimental and control groups must be provided or calculated; (3) dual-channel microarray studies were excluded; (4) DNA methylation arrays studies were excluded; (5) studies without a control group, as well as studies on other cancers were excluded; (6) studies with less than 20 cases were excluded. Differentially expressed genes (DEGs) identified by comparing expression in ovarian cancer and normal samples were analyzed using the R language software. Samples with both corrected P-values <0.05 and log fold changes (FC) >1 were deemed to be DEGs.

Gene set enrichment analysis (GSEA)

GSEA was performed using GSEA software 3.0 from the Broad Institute as previously described (Subramanian et al., 2005). GSEA evaluates gene expression profiles from samples belonging to normal and tumorous samples, and analyzes data in terms of gene sets. Normalized enrichment scores (NES) were acquired by analyzing genes with permutations 1,000 times. A gene set was considered to be significantly enriched if it had a normal p-value <0.05. However, if the number of samples per group was fewer than 7, gene set was selected as permutation type, and FDR <0.05 was set as the criterion of statistical significance.

Statistical analysis

Overall survival was calculated from the data of surgery to the date of death or last follow-up. Recurrence-free survival was defined as the time from the date of surgery to the date of first recurrence or last follow-up. Patients without events or death were recorded as censored at the time of last follow-up. Pathological stage and grade were considered to be distinct categorical variables. Age, longest dimension and residual tumor were included in the model as continuous variables. Gene expression levels from the TCGA database were included as continuous variables.

All data are displayed as mean ± standard deviation (SD) for each group. Student’s t-tests were performed to test for differences in mean values of vaviables of interest between two groups, whereas one-way analyses of variance (ANOVAs) were used to test for differences in means among three or more groups. Standardized mean difference (SMD) was applied to evaluate the association between FGFR1 levels and OC using RevMan 5.3.0. We pooled SMDs across GEO datasets using the Mantel-Haenszel formula (fixed-effect model) or the DerSimonian-Laird formula (random-effect model). A fixed-effect model was adopted when the Q statistic was considered significant (p > 0.1, or I²<50%); otherwise, a random-effect model was used.

Study Characteristics

The present study consists of several processes sequentially (Fig. 1), that is, TCGA-based RNA-seq data aggregation and clinical values, GEO-based data verification, meta-analyses based on GEO and TCGA, and multiple bioinformatics analyses. A total of twelve GEO datasets (GSE105437, GSE66957, GSE66387, GSE40595, GSE29450, GSE29156, GSE27651, GSE26712, GSE18521, GSE18520, GSE17308, and GSE12470) were collected for use in our study. Two datasets (GSE105437 and GSE12470) contained detailed information on Asian populations, and these were used in this study. The platform for GSE105437 was GPL570, the [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, includes 5 normal ovarian surface epithelial samples and 10 high-grade stage III invasive serous ovarian cancer samples. The platform for GSE12470 was GPL887, the Agilent-012097 Human 1A Microarray (V2) G4110B; this platform, included 10 normal peritoneum samples, 35 advanced stage serous ovarian cancer samples and 8 early stage serous ovarian cancer samples. In addition, 308 cases with EGFR or FGFR expression data were screened in TCGA database excluding 5 cases with recurrence, and 303 patients form this database were included in this study. A flowchart and other detailed information regarding the studies included in the meta-analysis are shown in Fig. 1.

FGFR expression in OC based on GEO

FGFR or EGFR expression was initially assessed in a series of ovarian tumor and normal tissues based on GEO dataset (Fig. 2). The datasets were analyzed with GeneSpringGX and mapped with Graphpad 6.0 software. The expression levels of FGFR1 in ovarian cancer tissues were significantly lower than in non-cancer control tissues in both GSE105437 and GSE12470 with analyses performing Fold change cut-off adjusted to 2.0 and corrected p-values cut-off adjusted to 0.05(p < 0.0001 and p = 0.0032, respectively; Figs. 2A and 2B). The expression of FGFR1 was significantly different in advanced ovarian cancer, but not in early ovarian cancer in GSE12470 (p = 0.0124; Fig. 2C). The expression levels of EGFR in ovarian cancer tissues were significantly lower than in non-cancer control tissues in GSE105437 (p = 0.0008; Fig. 2D), whereas no difference was found in GSE12470. The expression levels of FGFR2 and FGFR3 were significantly different only in GSE12470 (p = 0.0046 and p = 0.0049, respectively; Figs. 2E and 2F). Therefore, FGFR1 was chosen for subsequent analysis.

To have a more comprehensive understanding of the biomedical predictive value of FGFR1, ROC Curves were provided to investigate the diagnostic value of GSE105437 and GSE12470 in distinguishing OC tissues from normal controls. As shown in Fig. 3, the Area Under the Curve (AUC) of GSE105437 outperformed GSE12470 (AUC=1, p = 0.0009 and AUC = 0.8256, p = 0.0015 respectively). Therefore, we considered FGFR1 might play an important role in diagnosing ovarian cancer. In brief, no significant differences were found between OC and nontumor groups based on the GSE105437 and GSE12470 data considered here (SMD = −6.22; 95% CI, −7.48 to −4.96; p = 0.93, Fig. 4A). Moreover, we found no significant heterogeneity by fixed-effected model (p < 0.00001, I² = 0%). Finally, the funnel plot shown in Fig. 4B did not imply significant publication bias.

Identification of DEGs in ovarian cancer using integrated bioinformatics

The two ovarian cancer gene expression microarray datasets GSE105437 and GSE12470 were analyzed using the R limma package and were sorted according to logfold-change values with corrected p-values <0.05. The results of DEGs from the two databases showed in Fig. 5, and FGFR was identified both in the down-regulated genes (Green points).

GSEA enrichment analysis of DEGs

To find enriched pathways related to FGFR1 and to identify its potential function, we performed a GSEA (Fig. 6). First, common pathways were found by comparing the GSE105437 (Fig. 6A) and GSE12470 (Fig. 6B) datasets. Although, many of these common pathways were not significantly associated with FGFR1, we found that the KEGG adherent junction signaling pathway was significantly correlated with FGFR1 in the GSE105437 dataset (Fig. 6C).

Clinical characteristics of patients from the TCGA

The TCGA cohort consisted of 303 female patients with histologically confirmed ovarian cancer who had undergone prior surgical resection and received no pretreatment. Summary data of clinical indicators such as age, stage, grade, the longest dimension of the tumor, ANS, and residual tumor were shown in Table 1. We tested whether these variables were correlated with the relative expression of FGFR1 mRNA in OC tissues. We found no significant differences in the clinical features for all the tested parameters. The median follow-up time was 26.11 months and 182 patients died during follow-up. The results of the overall survival (OS) analyses, as calculated using the Kaplan–Meier method are shown in Table 2 and Fig. 7; we found significant differences in survival between all groups. We combined stage I and II into a single group, as well as grades G1 and G2 into another group for statistical analysis, since only 1 patient was included in the dataset for both stage I and grade G1. Residual tumors were found to be independently associated with increased OS according to univariate analysis (p = 0.001). In contrast, we found no statistical evidence that other characteristics, including age, stage, grade, tumor longest dimension, and ANS were independent prognostic factors. In addition, our results showed that patients with larger residual tumors, high FGFR expression, and advanced cancer stages were associated with shorter OS (p = 0.001, p = 0.233, p = 0.148, respectively; Fig. 7).