Identification of candidate biomarkers correlated with the diagnosis and prognosis of cervical cancer via integrated bioinformatics analysis

Introduction

Globally, cervical cancer is the fourth frequently occurring malignancy in women in the world, accounting for almost 12% of all female cancers. There are nearly 5,27,000 new cases and 2,65,000 deaths owing to cervical cancer annually. In fact, nearly 84% of cervical cancer cases occurred in developing countries, including China. Cervical cancer is a high burden for China, in virtue of its large population with geographical and socioeconomic inequities.¹ At present, cervical cancer screening methods commonly used in clinical practice include Pap Smear, TCT, Hybrid Capture 2 (HC2) test, and colposcopy. The therapy of cervical cancer mainly includes surgical treatment, radiotherapy, and chemotherapy. These screening and treatment tests have greatly decreased the incidence and mortality rates for cervical cancer.² However, the 5-year survival rate is still low, especially for advanced patients, 5-year disease-free survival rate (DFS) is only 50%.³ More importantly, there is still no specific biomarkers for the diagnosis and detection of cervical cancer prognosis. Therefore, it is extremely urgent to study the potential molecular mechanisms of malignant biological behavior of cervical cancer, and find more reliable and prospective biomarkers for diagnosing and monitoring recurrence and evaluating prognosis.

As an important component of bioinformatics, gene expression profiling is a promising and powerful tool in medical oncology with great potential clinical applications: from molecular diagnosis to molecular screening of cancers, from new targets discovery to tumor response prediction, from patients classification to prognosis prediction.⁴ Today, it has been widely used to study molecular changes involved in tumor progression to better understand the development processes in human carcinogenesis, and discover new prognostic markers and new therapeutic targets in various types of cancer.⁵ At present, a large number of studies on the expression profile of cervical cancer have been carried out to identify the differentially expressed genes (DEGs) related to the development process or prognosis of cancer. However, owing to different data processing methods and different sample sizes, the results of DEGs are incompatible and even contradictory.^6,7 Therefore, the results obtained by different raw data processing methods still have some limitations and insufficiency. Comprehensive bioinformatics methods have been applied to cancer research, and a large amount of valuable biological information has been discovered, making it possible to find effective and reliable molecular markers.⁸

In this study, we downloaded three original mRNA microarray datasets from the NCBI-Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), GSE7803,⁹ GSE9750,¹⁰ and GSE63514,¹¹ a total of 140 samples containing 58 normal cervical samples and 82 cervical cancer specimens. At the same time, a microRNA microarray, GSE30656,¹² was downloaded from the GEO database, containing 10 normal cervical samples and 19 cervical cancer samples. The differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs) in normal tissues and tumor samples were obtained by R software. DEGs protein interaction was performed on the STRING online website. Subsequently, Cytoscape was applied for further visualization. The HMDD online website was effectively used to search for the miRNAs proved to be associated with cervical cancer. And miRNA target genes were obtained using miRWalk2.0 online database. Finally, we analyzed whether the expression of these target genes was consistent with the expression of cervical tissue from the TCGA and GTEx databases, and these targets were associated with patient prognosis. In the end, we verified that the expression of KLF4 and ESR1 was significantly lower in cervical cancer tissues by immunohistochemistry method. In this article, we identified that KLF4 and ESR1 may be the crucial genes in the cervical carcinogenesis and prognosis.

Materials and methods

Microarray data

Using the keywords “cervical cancer” to search on the GEO database (https://www.ncbi.nlm.nih.gov/geo/), the raw gene expression profiles of GSE7803, GSE9750, and GSE63514 were downloaded. The platform for GSE7803 is GPL96, [HG-U133A] Affymetrix Human Genome U133A Array, which contains 10 normal cervical samples and 21 cervical cancer specimens; The platform for GSE9750 is also GPL96, containing 24 normal cervical epithelial tissues and 33 cervical cancer tissues. The platform for GSE63514 is the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, which consists of 24 normal cervical epithelial tissues and 28 cervical cancer tissues. The miRNA dataset GSE30656 was downloaded from GEO and its platform was GPL6955, Agilent-016436 Human miRNA Microarray 1.0 (Feature Number version), which included 10 normal cervical epithelial tissues and 19 cervical cancer tissues. The dataset information is shown in Table 1.

Table 1 Details for GEO cervical cancer data

Screening for DEGs

The downloaded platform and series of matrix files were converted using the R language software and annotation package. The ID corresponding to the probe name was converted into an international standard name for genes (gene-symbol) and saved in a CSV file. z-SCORE is a commonly used standardization method that can normalize values to a matrix with a standard deviation of 1 and a mean of 0 for subsequent integration of multiple sets of data. After processing the original data through the RMA function and removing the background in the R/Bioconductor software, three gene expression profiles were normalized by the z-SCORE function and processed into one data set, which contained 58 normal cervical tissue samples and 82 tumor tissue samples. The differentially expressed genes (DEGs) and the differentially expressed miRNA (DEMs) were obtained and saved in CSV files by the limma package function. | log2FC| ≥1 and adjust P-value <0.05 were considered to indicate a significant difference.

PPI network

The STRING database (http://www.bork.embl-heidelberg.de/STRING/) is a computerized and powerful global resource for studying the interactions between the predicted and experimental interactions of proteins.¹³ In the present study, the DEGs expression products in cervical cancer were constructed using the STRING database to construct PPI networks, and a combined score of ≥0.4 was used as the cut-off value. After removed the isolated and partially connected nodes, a complex network of DEGs was constructed, the excel of string-interactions was downloaded and saved in a CSV file. String-interactions file and DEGs expression file were imported into cytoscape software for further processing.

Ascertain of miRNA target genes

The HMDD tool (http://www.cuilab.cn/hmdd) was performed to identify the validated miRNAs in cervical cancer, the ‘Calculate and draw custom Venn diagrams’ tool was applied to identify the DEMs within the list of validated miRNAs in cervical cancer. and miRWALK2.0 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) was used to predict the target genes of DEMs.

DEGs basic expression in normal and cancer tissues

The Cancer Genome Atlas (TCGA) is a collaboration between the National Cancer Institute (NCI) and the National Human Genome Research Institute (NHGRI) that has generated comprehensive, multi-dimensional maps of the key genomic changes in 33 types of cancer and corresponding normal tissues.¹⁴ GEPIA (http://gepia.cancer-pku.cn/index.html) is a new web-based tool that uses standardized, integrated processing to analyze tumor and normal gene expression from cancer genome mapping (TCGA) and genotype tissue expression (GTEx). It includes 319 cervical tissue samples (306 tumor samples and 13 normal tissue samples) from the TCGA database.¹⁵ We validated the expression of DEGs in these cervical tissue samples.

Survival analysis

OncoLnc database (www.oncolnc.org) contains survival data for 8,647 patients from 21 cancer studies from TCGA.⁶ The tool was used to conduct an overall survival analysis (OS) for patients with cervical cancer. The 5-year survival of patients with cervical cancer was analyzed by the HPA dataset (https://www.proteinatlas.org/). They allow the production of Kaplan–Meier plots stratified by gene expression levels. In the OncoLnc database, Log-rank P-values in survival analysis were recorded, 90th (upper) percentiles and 10th (lower) percentiles were considered as high and low groups.

Immunohistochemistry

Six paraffin-embedded human cervix normal and tumor tissues were cut into 5 μm thick, 3 mm diameter sections to construct tissue microarrays, and then they were placed in 3% hydrogen peroxide solution to block endogenous peroxidase activity for 10 mins, immunohistochemical staining was carried out using a rabbit polyclonal primary antibody to KLF4 (1:200; Mitaka,11880–1-AP) and ESR1 (1:300; Mitaka, 21244–1-AP) for incubating at 4°C overnight. DAB solution was used for coloration, and hematoxylin was used for counterstaining. The final immunoreactive scores were evaluated independently by two pathologists blinded to the tissue parameters. The scores were determined by the sum of intensity and extent of the staining. The integrated option density (IOD) represents the score of intracellular ESR1 expression. The percentage of nuclear KLF4 staining of cervical tissue was grouped on a scale of <5%, 5–25%, 25–50%, and >50%.

Statistical analysis

GraphPad Prism5 was used for statistical analysis, independent samples t-test was used to analyze the correlation between KLF4, ESR1 expression, and tissue parameters. P≤0.05 was considered as statistical significant.

Results

Identification of DEGs and DEMs in cervical cancer

The cervical cancer raw expression microarray datasets GSE7803, GSE9750, and GSE63514 were processed and normalized, and the results are shown in Figure 1. Three datasets were integrated into one, including 58 normal cervical samples and 82 cervical cancer specimens, which were screened by the limma package, 83 differentially expressed genes (DEGs) were obtained (Data are shown in Table 2). Total 12 upregulated genes and 71 downregulated genes were identified. Among them, CDKN2A, SPP1, SYCP2, MMP12, and CXCL8 were upregulated genes, and ESR1, PPP1R3C, NSG1, TMPRSS11D, GYS2, CDKN2A, ENDOU, and KLF4 were downregulated. The log FC of top 10 up and down-DEGs is shown in Table S1. Volcano plots of gene expression profile data in 140 cervical samples are shown in Figure 1SA. The cluster heatmaps of 83 DEGs in all samples are shown in Figure S2A. The miRNA dataset GSE30656 comprised 10 normal cervix tissue samples, and 19 squamous cell carcinoma or adenocarcinoma tissue samples, 14 differentially expressed miRNAs (DEMs) were obtained using the limma package, including 3 upregulated miRNAs and 11 downregulated miRNAs (Table 3). Volcano plot of miRNA expression profile data and heatmap of DEMs are shown in Figures S1B and S2B. The PPI network complex contains 52 nodes and 100 edges (Figure 2).

Table 2 Screening DEGs in cervical cancer by integrated microarray

Table 3 Screening DEMs in cervical cancer

Figure 1 Normalization of gene expression. (A) The normalization of GSE7803 data, (B) The normalization of GSE9750 data, and (C) The normalization of GSE63514 data. The blue bar represents the data before normalization, and the red bar represents the normalized data.

Figure 2 PPI network. Circles represent genes, lines represent the interaction of proteins between genes. Red indicates that the expression of genes is relatively upregulated, green indicates that the expression of genes is relatively downregulated. The size of the node represents the P-value. (The smaller the P-value, the larger the node). The line color represents the combined score of the interaction between the proteins (Brown represents stronger contact, and yellow means weaker contact). Abbreviation: PPI, protein–protein interaction.

Validated miRNAs in cervical cancer

The HMDD tool was fully utilized for searching the validated cancer-related miRNAs, there are total 47 miRNAs which have been reported and validated in cervical cancer (Table 4). The ‘Calculate and draw custom Venn diagrams’ tool was used to identify the DEMs within the list of validated miRNAs, as shown in Figure 3, there are three identified DEMs within the validated miRNAs, including hsa-miR-21, hsa-miR-16, and hsa-miR-205. In the three DEMs, hsa-miR-21 and hsa-miR-16 are upregulated and hsa-miR-205 is downregulated.

Table 4 miRNAs which have been reported and validated in cervical cancer

Figure 3 The identified DEMs within the list of validated miRNAs.

miRNA‑DEG pairs

MicroRNAs (miRNAs) are non-coding RNAs of approximately 20–25 nucleotides in length that recognize target mRNA by base pairing and up- or downregulate translation of target mRNA. The target genes of the validated DEMs in cervical cancer were performed by online databases miRWALK2.0. It was identified that three validated DEM targets were DEGs. ESR1, PPP1R3C, NSG1,and TMPRSS11D were the gene targets of hsa-miR-21 (Figure 4A); GYS2, CDKN2A, ENDOU, and KLF4 were the targets of hsa-miR-16 (Figure 4B), while hsa-miR-205 targets were DSG1, SPINK5, and IVL (Figure 4C). Among them, CDKN2A was upregulated, while other target genes were downregulated in the DEGs; hsa-miR-21 and hsa-miR-16 were upregulated, while hsa-miR-205 was downregulated in the DEMs. In view of miRNAs negatively regulating the expression of their targets, we analyzed the correlation of upregulated genes with downregulated miRNAs, downregulated genes and upregulated miRNAs. Then, we made the visual networks of hsa-miR-21 (Figure 5A), hsa-miR-16 (Figure 5B), and their targets by cytoscape software.

Figure 4 The target DEGs within the validated miRNAs. (A) hsa-miR-21; (B) hsa-miR-16; (C) hsa-miR-205. Blue represents the 83 identified DEGs, red represents the target genes of miRNA, and brown represents the intersection of them. Abbreviation: DEGs, differentially expressed genes.

Figure 5 The networks of validated miRNAs and their target genes in cervical cancer. (A) hsa-miR-21; (B) hsa-miR-16. Red nodes represent upregulated miRNAs, while green nodes represent downregulated genes. Abbreviation: DEGs, differentially expressed genes.

DEGs expression state in normal and cervical cancer samples

GEPIA, an open-access website of the genome-wide expression in normal and cancer tissues from TCGA and GTEx, was used to identify the basic expression level of differentially expressed genes in 319 cervix uteri tissues. In the present study, we would like to validate whether the expression of miRNA target genes from the GEO database is consistent with their expression in these cervical tissues from TCGA and GTEx. As expected, KLF4, the target genes of hsa-miR-16 was significantly downregulated in cervical cancer samples, as were ESR1 and PPP1R3C. However, the targets ENDOU, GYS2 and TMPRSS11D had no statistical significance, and the gene NSG1 cannot be detected, as shown in Figure 6. Therefore, we chose KLF4, ESR1, and PPP1R3C for subsequent research.

Figure 6 The expression of target genes in cervical tissues from TCGA data and GTEx data. The target genes of hsa-miR-16 (A) The expression of KLF4, (B) The expression of ENDOU, (C) The expression of GYS2; the target genes of hsa-miR-21 (D) The expression of ESR1, (E) The expression of PPPIR3C, (F) The expression of TMPRSS11D. Red represents tumor samples, while gray represents normal samples. Abbreviations: KLF4, Kruppel like factor 4; ENDOU, Endonuclease; GYS2, Glycogen synthase 2; ESR1, Estrogen receptor 1; PPPIR3C, Protein phosphatase 1 regulatory subunit 3C; TMPRSS11D, Transmembrane protease, serine 11D; CESC, Cervical squamous cell carcinoma and endocervical adenocarcinoma.

Survival analysis

In this study, the relationship between the expression of KLF4, ESR1, PPP1R3C, and the overall survival of patients with cervical cancer was analyzed by OncoLnc tool. We found that the low mRNA expression of KLF4 was associated with poor overall survival rate (log-rank P=0.00966; Figure 7A). However, there was no significant in ESR1 and PPP1R3C (log-rank P<0.05; Figure 7C, Figure 7E). At the same time, the 5-year survival of patients with cervical cancer was analyzed by the HPA dataset; it was obvious that the low mRNA expression of KLF4 and ESR1 was associated with poor 5-year survival rate (P<0.05; Figure 7B and D). However, the high mRNA expression of PPP1R3C had a poor 5-year survival rate (Figure 7F). Overall, KLF4 is closely related to overall survival, and ESR1 is closely related to 5-year survival. KLF4 and ESR1 might be critical biomarkers for the prognosis of cervical cancer.

Figure 7 The prognostic value of three differentially expressed genes in cervical cancer patients. The overall survival rate (A) Prognostic value of KLF4, (C) Prognostic value of ESR1, (E) Prognostic value of PPPIR3C; the 5-year survival rate (B) Prognostic value of KLF4, (D) Prognostic value of ESR1, (F) Prognostic value of PPPIR3C. Abbreviations: KLF4, Kruppel like factor 4; ESR1, Estrogen receptor 1; PPPIR3C, Protein phosphatase 1 regulatory subunit 3C.

The expression of KLF4 and ESR1 in cervical tissues by immunohistochemistry

To further determine whether the two key genes obtained above were downregulated in cervical cancer, the protein expression in cervical tissues was detected by immunohistochemistry method. The results were consistent with those obtained by bioinformatics, the expression of KLF4 and ESR1 in cervical cancer tissues was significantly lower than that in normal cervical tissues (P<0.05; Figure 8).

Figure 8 The immunohistochemical analysis of ESR1 and KLF4 in cervical cancer. (A, B) A normal cervical tissue showing high ESR1 and KLF4 expression. (C, D) A cervical cancer tissue showing low ESR1 and KLF4 expression. (E, F) Immunohistochemistry score results indicated that ESR1 and KLF4 staining scores were lower in tumor tissues than non-tumor tissues. Blue represents non-tumor tissues, while red represents tumor tissues. Abbreviations: KLF4, Kruppel like factor 4; ESR1, Estrogen receptor 1.

Discussion

As one of the most common malignant gynecological tumors, cervical cancer is one of the leading causes of cancer death in women. According to the latest statistics released in 2015, China ranks second in new cases of cervical cancer every year in the world.¹⁶ According to the statistics, 80% of patients have developed aggressive cancers at the time of diagnosis, and the age of diagnosis is slowly decreasing. In addition, cervical cancer still carries a high risk of morbidity and mortality on account of the metastasis and recurrence.¹⁷ Meanwhile, reliable and specific biomarkers for diagnosis and prognosis of cervical cancer are deficient and unexplored. Therefore, it is essential and urgent to search for biomarkers of diagnosis and novel treatment targets to predict the survival of cervical cancer. As an indispensable part of bioinformatics, gene expression microarrays have been extensively used to study cancer-related genes and provide broad prospects for drug-based molecular targeting, molecular prediction, and molecular therapy.¹⁸ A series of biomarkers have been presented potentially as the diagnosis and prognosis in cervical cancer and other cancers.^19,20

In the present study, three raw mRNA and one miRNA expression profile datasets were used to identify potential biomarkers in cervical cancer. The expression patterns obtained from the GEO database included 3 human cervical cancer raw mRNA profiles and 1 human cervical cancer miRNA expression dataset, all of which contained a comparison between cervical cancer and normal cervix tissue samples. After using R software and bioinformatics to integrate and deeply analyze these datasets, there are total of 83 DEGs and 14 DEMs that were identified.

MiRNA, a non-coding RNA of approximately 22 nucleotides in length, regulates expression of various oncogenes and tumor suppressor genes by targeting the 3’UTR of the target mRNA, resulting in degradation or inhibition of translation.²¹ In the study, we obtained the validated miRNAs by reference to existing studies using HMDD and identified target genes of the DEMs by miRWALK2.0. As miRNAs negatively regulate the expression of their target genes, we identified that ESR1, PPP1R3C, NSG1, and TMPRSS11D were potentially the targets of miRNA-21, KLF4, ENDOU, and GYS2 were the gene targets of miRNA −16. Among them, miRNA-21 and miRNA-16 were upregulated, while the targets ESR1, PPP1R3C, NSG1, KLF4, ENDOU, and GYS2 were downregulated.

To further study whether there is a difference in the expression of above genes in normal and tumor tissues, we analyzed their expression in 319 cervical tissue samples from the TCGA and GTEx databases, including 306 tumor samples and 13 normal tissue samples. As expected, the target gene KLF4 was significantly downregulated in cervical cancer samples, as were ESR1 and PPP1R3C. However, the expression of ENDOU, GYS2, TMPRSS11D, and NSG1 is upregulated or missing in cancer tissues. Subsequently, further analysis was performed using the survival data from TCGA via OncoLnc and HPA website to identify the relationship between the postoperative survival of patients and the expression of above three genes. It suggested that the low expression of KLF4 was closely associated with poor overall survival and 5-year survival time of cervical cancer patients (P<0.05), the low expression of ESR1 and the high expression of PPP1R3C were closely associated with poor 5-year survival rate (P<0.05). However, ESR1 and PPP1R3C expressions were not significantly associated with OS (P>0.05). Overall, KLF4 and ESR1 are closely related to overall survival and 5-year survival. They might be of great clinical significance for the prognosis of cervical cancer and provide effective targets for the treatment of cervical cancer.

Kruppel like factor 4 (KLF4) located at 9q31 is a critical DNA-binding transcription regulator.²² It is highly expressed in the gastrointestinal tract and other epithelial tissues, and it has been reported to regulate cell proliferation, apoptosis, differentiation and the maintenance of telomerase activity.^23,24 A number of studies have suggested that KLF4 functions as a tumor suppressor in various cancers, such as human pancreatic cancer and colorectal cancer.^25,26 Yang et al indicated that KLF4 was downregulated in the development and progression of cervical carcinoma, and the overexpression of exogenous KLF4 protein inhibited the cell growth and tumor formation by activating the cell cycle suppressor p27, which supports the hypothesis that KLF4 is a tumor suppressor in cervical carcinoma.²⁷ However, there is still no report on the relationship between KLF4 and the prognosis of cervical cancer.

Estrogen receptor 1 (ESR1) is a receptor protein, which is usually enriched expression in breast, endometrium, and normal cervical tissue, which is necessary for sexual development and reproductive function, and performs an important role in cellular development and differentiation. ESR1 protein is closely related to a variety of gynecological oncology diseases. It is well known that loss of ESR in breast cancer patients indicates invasiveness and poor prognosis.²⁸ In ovarian carcinoma, ESR1 methylation influences ovarian cancer development, and its promoter methylation was negatively correlated with survival in the subgroup of low-grade ovarian carcinoma patients.²⁹ ESR1 also plays a critical and diverse role in cervical cancer. Many previous studies indicated that ESR1 was downregulated in cervical malignant cells and human papillomavirus (HPV) DNA positive cervical intraepithelial neoplasia tissues.^30,31 On the other hand, ESR1 plays a major role in mediating cervical cancer invasion and progression, and loss of it enhances cervical cancer invasion.³² However, the expression of ESR1 in advanced cervical cancer and its impact on survival needs further study.

A variety of miRNAs have also been reported to downregulate KLF4 and ESR1 in various types of disease. In gastric cancer, KLF4 was a direct target of miR-103; And in non-small cell lung cancer cells, miR-25 could enhance cell migration and invasion via ERK signaling pathway by inhibiting KLF4.^33,34Similarly, it has been reported that miR-301a-3p suppresses estrogen signaling by directly inhibiting ESR1 in ERα positive breast cancer.³⁵ However, the relationship between miRNAs and KLF4 and ESR1 in cervical cancer is rarely reported.

miRNA-21 and miRNA-16 have been validated to be closely related to cervical cancer. The expression of miR-21 was upregulated by the HPV16 E7 oncoprotein in vivo, and the high expression of miR-21 in cervical cancer may play a possible role of oncogenes.^36,37 miR-16-1 belongs to the miR-16 cluster and has been implicated in various aspects of carcinogenesis including cell proliferation and regulation of apoptosis. In human cervical cancer cells, miR-16-1 post-transcriptionally downregulates Cyclin E1 (CCNE1) gene expression to deregulate the cell cycle .³⁸ The above indicates that miRNA-21 and miRNA-16 might act as oncogenes to promote carcinogenesis in a variety of pathways. In this study, using comprehensive bioinformatics analysis methods, we found that miRNA16 and miRNA-21 might lead to the pathogenesis of cervical cancer by inhibiting the expression of KLF4 and ESR1 in cervical tissues, respectively, and cervical cancer patients with low expression of KLF4 and ESR1 may be associated with poor overall survival and 5-year survival prognosis. At the end, we found the expression of KLF4 and ESR1 protein was significantly lower in cervical cancer tissues via experiment method.

Conclusion

In summary, in this study, we find that KLF4 and ESR1 are downregulated by the upregulated miR21 and miRNA16 in cervical cancer, respectively, and the lower expression of KLF4 and ESR1 is closely related to the poor prognosis. They might be critical biomarker for diagnosis and predicting cervical cancer progression. Our study has important clinical significance for better understanding the development and prognosis of cervical cancer. Of course, whether these genes can be used as the biomarkers for diagnosis and prognosis of cervical cancer needs further clinical study.