Integrative machine learning analysis of multiple gene expression profiles in cervical cancer

Although most of the cervical cancer cases are reported to be closely related to the Human Papillomavirus (HPV) infection, there is a need to study genes that stand up differentially in the final actualization of cervical cancers following HPV infection. In this study, we proposed an integrative machine learning approach to analyse multiple gene expression profiles in cervical cancer in order to identify a set of genetic markers that are associated with and may eventually aid in the diagnosis or prognosis of cervical cancers. The proposed integrative analysis is composed of three steps: namely, (i) gene expression analysis of individual dataset; (ii) meta-analysis of multiple datasets; and (iii) feature selection and machine learning analysis. As a result, 21 gene expressions were identified through the integrative machine learning analysis which including seven supervised and one unsupervised methods. A functional analysis with GSEA (Gene Set Enrichment Analysis) was performed on the selected 21-gene expression set and showed significant enrichment in a nine-potential gene expression signature, namely PEG3, SPON1, BTD and RPLP2 (upregulated genes) and PRDX3, COPB2, LSM3, SLC5A3 and AS1B (downregulated genes).


INTRODUCTION
Cervical cancer, ranked 4th most common cancer (World Cancer Research Fund International, 2017) among women worldwide, caused 279,311 mortality cases in the year 2015 (World Health Organization, 2016). It is the 3rd most common cancer (7.7%) among Malaysian females (AbM, Nor Saleh & Noor Hashima, 2016). Most of the development of cervical cancer is linked to the infection of the Human Papillomavirus (HPV) that passes to the human host during sexual activity. The other mentioned risk factors include smoking, imbalanced diet and family history (Plummer et al., 2003;Ames & Gold, 1997;Lee et al., 2003). However, as in many cancers, it is conceivable that genetic pathways play a role in the development of cervical cancer.
Microarray technology has been used to simultaneously identify a large number of gene expressions that may be associated with a particular cancer. Reports of prognostically predictive genes have been identified by whole genome microarray in cancer. Hu et al. (2010) showed that miR-200a and miR-9 were associated with cervical cancer outcome while Aziz et al. (2016), reported on 19 common genes that were predictive in colorectal cancers. Meta-analysis is touted to solve the problems of low statistical power of studies, to be able to analyze intersected genes between datasets, extract the optimum values across multiple datasets and identify the most commonly expressed genes (Warnat, Eils & Brors, 2005;Brazma et al., 2001). Based on our literature review, we believed that this is one of the first few studies which involves meta-analysis using multiple cervical cancer microarray datasets.
Despite the many approaches which were introduced in the previous meta-analysis study, a biased issue that is caused by merging datasets from different platforms and experiments is a hot topic in the field. Alles et al. (2009), Chow, Alias & Jamal (2017, Dong et al. (2010) and Li et al. (2016) carried out meta-analysis in various type of cancers to find the significant genes in certain diseases such as breast cancer, paediatric B-acute lymphoblastic leukaemia (B-ALL), ovarian cancer, etc. Trendily, most of these studies were done by using the original datasets that varied from different platforms. As the original datasets contained noisy data such as the extreme values, thus the bias effects from each of the respective datasets might be strong.
Furthermore, until the present, most of the previous studies either used statistical analysis methods (Alles et al., 2009;Chow, Alias & Jamal, 2017;Dong et al., 2010;Li et al., 2016) or machine learning-based methods (Chow, Alias & Jamal, 2017;Dong et al., 2010;Li et al., 2016;Grützmann et al., 2005;Liu et al., 2004;Obayashi & Kinoshita, 2009;Wang & Makedon, 2004;Díaz-Uriarte & De Andres, 2006) in carrying out the meta-analysis. Due to different measurements being used in analysing the expression activities in the genes, there is no way in comparing the feasibility of both methods in a meta-analysis. Hence, by combining both the statistical and machine learning methods, meaningful insights could be discovered. Favourable features from both of the methods are integrated and fully utilized in carrying out the meta-analysis in this study.
To address the shortcomings of meta-analysis we attempted to build a system which integrates machine learning with meta-analysis to improve the quality of the data acquired. The whole process encompasses (i) gene expression analysis of individual dataset, (ii) meta-analysis of multiple datasets, and (iii) feature selection (FS) & machine learning analysis. We report our findings using our proposed integrative analysis in identifying genes which could be integral in the development of cervical cancers, a cancer in which understanding has moved significantly forward and awaits unravelling of possible genetic intermediaries in the carcinogenetic process.

MATERIALS AND METHODS
In brief, our approach commences with the gene expression analysis of individual cervical cancer microarray datasets obtained from publicly available sources and to detect the genes that perform most significantly in each of the dataset. The most significant gene sets from each dataset are combined and intersected to form a larger dataset in the first stage. Microarray gene expression meta-analysis is performed using the combined dataset in the second stage. For the final stage, several FS and machine learning methods as well as a functional analysis by gene set enrichment analysis (GSEA) will be carried out to identify a more precise set of potential gene markers in cervical cancer biology.

Microarray datasets
All the data used are from the public website National Center Biotechnology Information Gene Expression Omnibus (NCBI GEO) (http://www.ncbi.nlm.nih.gov/geo). A search using the keywords ''cervical cancer'' and ''microarray'' was carried out. To ensure the accurateness and purity of the resultant 104 (before 5th January 2017) datasets, filtration was performed by manually reviewing each dataset using the following criteria. To gain entry into this study, the dataset must (1) contain expression profiling by array types, (2) use CEL raw files, and (3) be studies that involve cervical cancer as the main subject. Finally, only four cervical cancer microarray gene expression datasets (see Table 1) could be entered for further analysis. All analyses of the data were performed in R Core Team (2016).

Integrative approach of meta-analysis and machine learning
The aim of the proposed integrative approach is to fully utilize the benefits from both of the statistical analysis and the machine learning approaches in carrying out gene expression meta-analysis, in order to find out the most significant gene markers in cervical cancer. Generally, the proposed integrative analysis is composed of three stages, namely: (i) gene expression analysis of individual dataset; (ii) meta-analysis of multiple datasets and (iii) feature selection and machine learning analysis. The proposed integrative method is illustrated in Fig. 1. At the first stage, gene expression analysis of individual dataset was carried out where 10,000 differentially expressed genes from each dataset were selected by using Linear Models for Microarray Analysis (LIMMA). This process was repeated on four datasets. In stage (ii), meta-analysis of four datasets was carried out. Intersection process was done in order to obtain the intersected genes among the four datasets. A set of upregulated and downregulated genes were selected and ranked accordingly by using Rank Prod at the end of stage (ii). In stage (iii), eight FS and machine learning approaches were applied in order to identify the common optimum genes. Stages of the integrative approach are as detailed below:

Stage (i): gene expression analysis of individual dataset
The data were first filtered to select out cancer samples (removal of normal samples) to ensure that the activities of the genes presented in the normal samples would not affect the analysis. Among the four sets of microarray data used in our study, only microarray data of GSE9750 (Scotto et al., 2008) is based on the Affymetrix Human Genome U133A Array [HG-U133A] platform. The other three datasets are based on Affymetrix Human Genome U133 Plus 2.0 Array . Therefore, it is necessary to pre-process and normalize all of the raw level probe data across the platforms to reduce bias caused by the different platforms. A quantile algorithm was applied to normalize the data and the data were log2-transformed according to all the samples median.
Then, gene expression analysis was performed on each of the microarray datasets. RMA normalization from affy function in Bioconductor package in R (Gentleman et al., 2004) was carried out for every sample in the dataset to create expression values that are normalized and background corrected.
After that, differential gene expression analysis was carried out to select the top 10,000 differentially expressed genes from individual dataset using Linear Models for Microarray Analysis (LIMMA) (Smyth, 2004). In this stage, the rare and over-extremely expressed genes in every dataset were eliminated by using moderated t -statistics of LIMMA so that bias errors from each dataset could be reduced when combining the datasets to become a larger dataset.
The R package limma (Ritchie et al., 2015) was used to rank and select out the top differential expressed genes. The differentially-expressed genes were filtered using adjusted p-value of <0.05 and ranked based according to the p-values. In this study, the same gene expression analysis was repeated on all the four datasets.
Hence, this stage is treated as quality control of the dataset which is going to use in stage (ii) analysis as genes produced by this stage are consistently different between the group according to the LIMMA statistics.

Stage (ii): meta-analysis for multiple datasets
The aim of this stage is to select the genes that show upregulated and downregulated activities in cervical cancer. Four microarray gene expression datasets that contained top 10,000 genes each (selected in stage(i)) were combined in the intersection process by identifying the intersected genes found among the four datasets and merging all samples to form a single, larger sample size dataset (meta-dataset).
Next, meta-analysis was conducted by using RankProd (Breitling et al., 2004) package in R to analyse the intersected genes. With the default value of p < 0.05, FDR < 0.05 as the cut-off point, RankProd ranked the genes into two group: upregulated (FC < 1, php < 0.01) and downregulated genes (FC < 1, php < 0.01).

Stage (iii): feature selection (FS) and machine learning methods
Stage (iii) is to select the optimum genes to ensure the final set of genes selected is with higher standard and high quality in term of genetic evaluation. This is also to fully utilize the measurements of different machine learning methods (supervised and unsupervised) by collaborating them in selecting the significant genes. The machine learning methods chosen in this study are based on the literature review done which obtained good result in the related studies (Urbanowicz et al., 2017;Wang et al., 2005;Martín-Valdivia et al., 2008;Guyon et al., 2002;Breiman, 2001;Díaz-Uriarte & De Andres, 2006;Sneath & Sokal, 1973;Eisen et al., 1998).
Seven supervised machine learning FS methods and an unsupervised one were used in this study. Among the seven supervised machine learning FS methods, four filter approaches (Pearson Correlation Coefficient(PCC), Relief-F Feature Selection (Relief-F), Correlation Feature Selection (CFS) and Information Gain (IG)) were used in this study for their simplicity and high in good succession in performing FS. The features (genes) that selected by the filter approach are based on the relevance of the genes with the surrounding data (Chandrashekar & Sahin, 2014;Saeys, Inza & Larrañaga, 2007). In addition, one wrapper approach, Sequential Forward Selection (SFS), was carried out in this study as it acts as a searching algorithm that could find a subset of genes that is maximum according to the objective function (classification performance) (Chandrashekar & Sahin, 2014;Saeys, Inza & Larrañaga, 2007). Furthermore, this study also applied two embedded approaches, Support Vector Machine-based Recursive Feature Elimination (SVM-RFE) and Random Forest (RF), which embed variable selection as part of their learning process and evaluate them according to the classification result of specific machine learning method (Chandrashekar & Sahin, 2014;Saeys, Inza & Larrañaga, 2007). Hierarchical Clustering (HC) was applied as the unsupervised model in this study for its ability to select the genes without the knowledge of the class labels (Chandrashekar & Sahin, 2014).
As the supervised machine learning methods involves the classification according to the class labels, the normal samples were added back at stage (ii) so that the ''learning procedures'' of the machine learning can be completed.
The machine learning approaches applied in this study are more on exploratory models which the performances of the machine learning methods are not being emphasize here. The main objective of this stage is to obtain the optimum genes that commonly selected using different machine learning methods, regardless of the differences in the algorithm used.
The results (selected genes) of different FS and machine learning will be combined and converted in the scoring matrix. Scores will be determined and counted by the frequency of the genes that was being selected by various FS and machine learning methods.

Functional Pathway Analysis
Gene Set Enrichment Analysis (GSEA) GSEA (http://software.broadinstitute.org/gsea/index.jsp) (Subramanian et al., 2005;Mootha et al., 2003) computationally evaluates the statistical significance of a priori defined set of genes with their biological states (i.e., phenotypes). All the genes in the datasets are ranked according to their differential expression between phenotype groups. In this study, normal tissue samples are also added into the dataset to enable the GSEA to calculate the Enrichment Scores (ES) between the phenotypes. A computational annotated gene set (C4) that consists of large collections of cancer-orientated microarray data is provided by Molecular Signature Database (MSigDB) and linked to GSEA in order to calculate the ES on our selected set of genes. ES is computed for every set of genes by going through the rank list from the top. A gene set will have a higher ES if the encountered gene in the gene input list is a member of the annotated gene set. A gene set is said to be significant if the normalized enrichment scores (NES) is high and its corresponding false discovery rate (FDR) is low. In this study, FDR value of <0.25 was set as the threshold to consider the significance of the gene enriched. Furthermore, pathway network could be visualized through Cytoscape (Shannon et al., 2003) software based on the GSEA results.

RESULT Stage (i): gene expression analysis of individual dataset
LIMMA was used to perform the analysis in order to select out the top 10,000 differentially expressed genes from each of the individual dataset (Results shown in Appendix S1).
LIMMA tests the average difference in log expression levels of two groups per gene by using moderated t -test. Therefore, genes selected by LIMMA are moderately different with average differential values, rather than extremely different.

Stage (ii): meta-analysis of multiple dataset
In the intersection process, 117 intersected genes among four cervical cancer datasets were identified. All cervical cancer samples from four datasets (total of 142 samples) were combined to form a meta-dataset with these 117 intersected genes. The meta-dataset is attached in Appendix S2.
RankProd estimated the FDR for each gene and ranked the gene according to the value that has equal or smaller p-value in the whole list of genes (estimated percentage of false predictions). As the product of RankProd is affected more by extreme expression values, genes that are extremely expressed in at least one sample in the dataset tend to be selected. Tables 2 and 3 show the results of RankProd analysis.

Stage (iii): feature selection and machine learning methods
To find the optimum set of genes, eight FS and machine learning methods namely PCC, Relief-F, CFS, IG, SFS, SVM-RFE, RF and HC, were carried out using the set of genes resulting from stage (ii) of meta-analysis.
For HC analysis, the expression patterns of the 32 upregulated genes are shown in the heatmap and were further sub-clustered by cutting tree at the maximum heights of 1.5 of the row cluster. As a result, a cluster group that contained 20 genes was found to be most correlated to each other ( Fig. 2A). The same rule was applied to the 33 downregulated genes. The most correlated group that contained 26 genes was also identified (Fig. 2B). The genes that were selected in the cluster that had the highest scores are presented in Tables 4 and 5 respectively.
As for the scoring measures, the final score obtained by each of the genes is the frequency it was selected by the listed FS and machine learning methods (1 score for selected, 0 score for not selected). The higher the score of a gene, the more frequently the gene was selected by the different FS methods, and the higher the potential for the gene to be a gene marker of cervical cancer. Tables 4 and 5 show the set of genes selected by the FS methods. The up-and downregulated genes with score ≥ 6 (full score = 8, threshold = 75%) were identified and selected. In total, 21 gene expression set consisted of 9 upregulated genes and 12 downregulated genes were selected.

Functional Pathway analysis with GSEA
GSEA was carried out on the selected 21-gene expression set for functional analysis to measure the statistical significance of their biological phenotypes. In our study, we are assessing the enrichment of the selected 9 upregulated and 12 downregulated genes on their natural phenotypes, which are tumour and non-tumour. Therefore, samples of normal cells that were filtered at the beginning of the study were added back to complete the study. A computational gene set (C4) consisted of large collections of cancer-orientated microarray data was downloaded from MSigDB to calculate the enrichment scores of the nine upregulated genes and 12 downregulated genes among the annotated gene sets of C4. We evaluated the normalized enrichment scores (NES) and the corresponding false discovery rate (FDR) of each gene on all of the samples used in this study. In this study, a gene will be rejected when the FDR q-value is more than 0.25 (not significant). As NES values are dependent on the permutation method and size of the expression dataset, there is no standard threshold for this value in analysing the gene. A gene could match multiple sets of gene sets in the database, hence we selected the highest NES values for each respective gene as shown in Tables 6 and 7. The details of GSEA report on the annotation of the 21-gene expression set with the description of each gene sets are attached in Appendix S3. Among the nine upregulated genes, six genes were matched with the cancerous gene set collections, as shown in Table 6. Notably, PEG3, SPON1, BTD and RPLP2 passed the defined threshold of FDR <0.25, indicating the potential of these genes in showing their significances in cervical cancer. As there was very limited information regarding genes AK026407, TCF3 and TMEM115, there are no matching genes with the MSigDB. Besides that, nine downregulated genes were matched with the C4 gene sets in MSigDB (Table 7). This showed that the nine genes were enriched and highly related to cancers. Among enriched genes, five downregulated genes, PRDX3, COPB2, LSM3, SLC5A3 and AS1B, passed the default threshold with FDR <0.25. These also provided some attraction points for the selected downregulated genes in proving their significances in cervical cancer. However, there are no matching of the genes ARL6IP5, CNIH1 and PTS with the MSigDB due to the low information found regarding these genes. Figures 3 and 4 show the pathway network constructed according to the GSEA report generated in Appendix S3. The threshold that used in the pathway analysis was set to FDR <0.25 in order to standardize with GSEA. Each of the nodes in the networks represented the gene sets that the genes were involved in. The genes that showed enrichment with the Table 5 Ranking of each gene (downregulated) using proposed FS methods. gene sets were indicated in the colour nodes while the line connecting the nodes are the pathways that the enriched genes found in. In Fig. 3, the upregulated genes PEG3, SPON1, BTD and RPLP2 were shown to be related with the pathway of gene expression, disease or cancer, and protein metabolism. Similarly, as shown in Fig. 4, downregulated genes PRDX3, COPB2, LSM3, SLC5A3 and

Notes.
There was no matching of gene AK026407, TCF3 and TMEM115 with the MSigDB. ASF1B were involving in the gene expression, disease or cancer, protein metabolism and cell cycle pathways.

Integrative approach of meta-analysis and machine learning
The aim of this study was to identify a set of more precise, accurate and reliable potential cervical cancer gene markers. Generally, there are three stages of analysis in this study: (i) gene expression analysis of individual dataset; (ii) meta-analysis of multiple datasets; and (iii) FS and machine learning analysis. At the first stage, individual analysis of each cervical cancer microarray gene expression dataset was carried out using the LIMMA package in R. In LIMMA, the average difference of log expression levels between two groups of information provided for each gene was calculated using moderated t-statistic, to ensure that the final list of genes produced consists of genes that are consistently different between the groups. In other words, LIMMA avoids choosing the extremely different genes in the datasets. In this study, Benjamini and Hochberg's method in controlling the FDR value was used as the default adjustment method of the adjusted p-value. The adjusted p-value is set to be <0.05 in expecting the proportion of the false discoveries in the selected gene group to be less than 5%.
In the meta-analysis stage, RankProd package was chosen for use for two reasons: to analyse the combination of datasets from different origins and to select genes with higher biological significances. RankProd package is able to avoid errors that are caused from batch and platform differences in combining expression values of various studies, by performing ranking of the genes by their differential expression (Hong et al., 2006). RankProd estimated the gene ranking by randomly permutating the list of genes. False discovery rates (FDR) for the whole list of genes which are equal or smaller than the p-value (known as false positive prediction, pfp) were calculated by this package and the genes were ranked accordingly.
The selection of genes in the RankProd analysis in stage (ii) was based on the criteria of FC > 1 (upregulated genes) or FC < 1 (downregulated genes). In log expression level, a gene is said to be significantly upregulated when the FC expression is >1 while it is significantly downregulated when the FC <1. An FC with value of two simply brings the meaning that the particular genes expressed twice more than the genes in the original control. In addition the php value is set to <0.01 in order to select the genes that sited on the top 1% of the truly significant genes among all the genes to be studied.
Following meta-analysis, 32 most significantly up-regulated and 33 down-regulated genes were selected.
To further narrow down the range of the genes selected and also to fully adopt the benefits from the statistical and also the machine-based analysis method, four filter approaches (PCC, Relief-F, CFS and IG), a wrapper approach (SFS), two embedded approach (SVM-RFE and RF) and an unsupervised approach of HC were carried out. As the algorithms in selecting genes by each method were different, the resultant genes selected by each method varies. In this study, most of the parameters set are based on the default values or the minimum values required as the focus is on the most significant genes that selected by the machine learning method and not the predictive performances of the machine learning methods).
We ensured the accuracy of our results by creating a scoring system to calculate the score of each gene. The threshold for the scoring system was set to 75%, where only the genes with score ≥6 (out of 8) would be selected. In total, 21-gene expression set which consisted of nine upregulated genes and 12 downregulated genes were selected. The upregulated genes that were most commonly selected by the stated FS methods were PEG3, SPON1, BTD, RPLP2, CNKSR2, SPATA6, AK026407, TCF3 and TMEM115 while the FS-commonly-selected-downregulated-genes were PRDX3, COPB2, LSM3, SLC5A3, ASF1B, RPL15, HIGD1A, ARL6IP5, CNIH1, PSMD6, PDHB, and PTS.
• prognostic marker for gastric cancer (Zhang et al., 2011) • involved in the carcinogenesis and progression of various cancers (Sharp et al., 1990) • associated with other ribosomal protein genes with colorectal carcinomas (CRC) (Mao-De & Jing, 2007) that these nine genes were highly correlated with the formation of cancer and were selected as the nine potential gene expression signature in this study. The pathway network (Figs. 3 and 4) were constructed using the GSEA report that attached as Appendix S3. The networks shown were based on the correlation between the gene sets. Similar with GSEA, the threshold was set to FDR <0.25 in this analysis. By comparing the map constructed and the GSEA result, it can be clearly seen that both results are supporting each other. The genes that showed significant enrichment with the C4 gene sets were shown and constructed in the network too. Moreover, the genes involved in the pathways of gene expression, disease or cancer, protein metabolism and cell cycle (Figs. 3  and 4). This further strengthened the role of the selected genes to be the nine potential signature in relation with cervical cancer development.
The nine potential signature will be discussed and supported by the previous studies done in the following section.
Upregulated genes are the genes that are known to be overexpressed during the gene regulation stage, causing the development of the tumour or potential cancerous cell. The most significant upregulated genes set that formed after the GSEA analysis are PEG3, SPON1, BTD and RPLP2. Table 8 shows the function of the upregulated genes selected, together with the related diseases and description. Table 9 Function, related diseases and description of the most significant genes (downregulated).

Genes
Full Name Function related diseases and description PRDX3 Peroxiredoxin3 • encodes antioxidant function protein, provides protection to mitochondria from oxidative stress (Safaeian et al., 2012) • regulates the cellular Reactive Oxidative Species (ROS) level of the cell (Jiang, Sang & Qiu, 2017) • showed correlation with the severity of the cervical carcinoma (Kim et al., 2009).
• involves in the transduction of signal for G protein (Hu & Gao, 2012) • One of the differentially expressed genes in LNCaP prostate cancer cell lines (Coutinho-Camillo et al., 2006) LSM3

Hypothetical protein LOC285378
• Involves in the process of MRNA splicing • Involved in the rapid proliferation, invasiveness, oxidative phosphorylation and tumor size of cervical cancers (Lyng et al., 2006) SLC5A3 Solute Carrier Family 5 Member 3 • Involves in the osmoregulation of cells (Wright & Turk, 2004) • prevents the impairment of cell cellular function (Database GCHG, 2017) • differentially expressed between Parental SiHa Cells and SiHa/R Cells of cervical cancer (Chung et al., 2005) • involved in the formation of fusion transcripts gene in some prostate cancer cases (Narod, Seth & Nam, 2008) ASF1B Anti-Silencing Function 1B Histone Chaperone • codes for the substrate protein of the cycle regulated-kinase cell.
• derepress the overexpression of transcriptional silencing (Le et al., 1997) • associated with the aggressiveness of breast tumor (Corpet et al., 2011) • leaded to the poor disease outcome in the study of the cervical cancer proliferation cluster that corresponded to the 163 transcripts (Rosty et al., 2005).
Despite the focus of gene expression research in cancers are towards the upregulated genes, the role of downregulated genes should not be ignored. Lu et al. (2005) found that in general, tumours contained more downregulated miRNAs compared with the normal tissues, leading them to hypothesize that down-regulated expression of some miRNAs may play a key role in the development of tumours and that gene downregulation is as important an aspect to be studied as up-regulation.
The most significant identified downregulated genes in cervical cancer in this study are PRDX3, COPB2, LSM3, SLC5A3 and ASF1B. The biological functions, related diseases and description are shown in Table 9.
Our integrated approach of meta-analysis and machine-learning method successfully identified the potential genetic markers in cervical cancer, along with the evidence provided from the GSEA results and the previous study done. The identified genes can be considered as potential gene signatures in cervical cancer. However, further evaluation using other experimental methods are needed in order to validate the significance of the nine potentialgene expression cervical cancer signatures.

CONCLUSION
In this study, we proposed an integrative machine learning analysis of multiple gene expression profiles in cervical cancer. This integrative analysis is divided into three stages: (i) gene expression analysis of individual dataset; (ii) meta-analysis of multiple datasets; and (iii) FS & machine learning analysis with PCC, Relief-F, CFS, IG, SFS, SVM-RFE, RF and HC. As a result, nine upregulated genes and 12 downregulated genes, were selected through the three-step analysis. The result of the functional pathway analysis in GSEA showed the significant enrichment in nine potential gene expression signature namely PEG3, SPON1, BTD, RPLP2 (as upregulated genes) and PRDX3, COPB2, LSM3, SLC5A3, ASF1B (as downregulated genes). More experiments are needed in order to confirm the significance of the selected nine potential gene expression signature in cervical cancer.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This study was supported by the University of Malaya research grants with the project number of RP038C-15AET & BK041-2014. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: University of Malaya: RP038C-15AET, BK041-2014.