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

View article


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) ( 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.

Table 1:
Details of the four microarray datasets used in this research.
Study Platform Samples Genes Variables
Noordhuis et al. (2011) (GSE26511) HG-U133_Plus_2 39 54,675 Lymph node status: negative, positive{19}
Bachtiary et al. (2006) (GSE5787) HG-U133_Plus_2 33 54,675 Biopsy: 1{5}, 2{6}, 3{7}, 4{5}, 5{5}, 6{3}, 7{2}
Scotto et al. (2008) (GSE9750) HG-U133A 61 22,283 tissue: normal{24}, cell line{5}, cancer{28}
Den Boon et al. (2015) (GSE63514) HG-U133_Plus_2 130 54,675 tissue: normal{24}, CIN1 lesions{14}, CIN2 lesions{22}, CIN3 lesions{40}, cancer{28}
DOI: 10.7717/peerj.5285/table-1

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.

Framework for the proposed integrative approach of meta-analysis and machine learning in gene expression profiling of cervical cancer.

Figure 1: Framework for the proposed integrative approach of meta-analysis and machine learning in gene expression profiling of cervical cancer.

*DE- differential expression value

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 [HG-U133_Plus_2]. 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 ( (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.


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.

Next, the analysis of the 117 genes by RankProd identified 32 upregulated genes (FC > 1, pfp < 0.01) and 33 down regulated genes (FC < 1, pfp < 0.01), with p-value < 0.05 and FDR of 0.05 as the cutoff points.

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.

Table 2:
Results of RankProd analysis-32 upregulated genes.
ID_REF Identifier gene.index RP/Rsum FC: (class1/class2) pfp P.value
206731_at CNKSR2 102 2.697 1.13 7.77E−155 6.64E−157
209243_s_at PEG3 107 2.846 1.13 1.38E−149 2.36E−151
210549_s_at CCL23 110 4.156 1.19 3.15E−114 8.08E−116
216108_at LOC105373738 22 4.363 1.22 3.41E−110 1.17E−111
220298_s_at SPATA6 37 6.35 1.29 6.39E−81 2.73E−82
212797_at SORT1 20 6.371 1.33 9.09E−81 4.66E−82
207996_s_at LDLRAD4 71 10.32 1.66 1.72E−49 1.03E−50
208762_at SUMO1 59 11.02 1.66 8.37E−46 5.72E−47
207257_at EPO 84 12.25 1.48 4.26E−40 3.28E−41
202035_s_at SFRP1 73 12.53 1.50 5.99E−39 5.12E−40
204672_s_at ANKRD6 86 13.23 1.57 3.28E−36 3.09E−37
220994_s_at STXBP6 62 13.62 1.57 9.11E−35 9.34E−36
213652_at PCSK5 88 14.12 1.67 4.76E−33 5.29E−34
213994_s_at SPON1 50 14.35 1.61 2.72E−32 3.25E−33
221606_s_at HMGN5 35 15.39 1.68 5.13E−29 6.57E−30
203440_at CDH2 38 15.74 1,65 5.24E−28 7.17E−29
207046_at HIST2H4B 33 16.79 1.66 3.64E−25 5.29E−26
211865_s_at FZR1 99 19.24 1.69 1.47E−19 2.26E−20
215401_at AU147698 17 19.79 1.74 1.72E−18 2.79E−19
213068_at DPT 4 20.42 1.75 2.48E−17 4.24E−18
214117_s_at BTD 34 22.69 1.80 1.28E−13 2.30E−14
34187_at AK026407 100 23.65 1.83 2.78E−12 5.24E−13
213611_at AQP5 82 25.09 1.86 1.82E−10 3.58E−11
209465_x_at PTN 53 26.18 2.00 3.00E−09 6.16E−10
208606_s_at WNT4 13 27.42 1.92 5.31E−08 1.14E−08
211814_s_at CCNE2 52 33.28 2.00 0.001144 0.0002542
208790_s_at PTRF 81 33.49 2.00 0.00144 0.0003322
200908_s_at RPLP2 40 33.97 2.04 0.002479 0.0005933
219795_at SLC6A14 57 34.09 2.22 0.002761 0.0006842
205730_s_at ABLIM3 64 35.48 2.05 0.01187 0.003044
202648_at TCF3 32 37.13 2.10 0.05047 0.01337
216267_s_at TMEM115 95 37.72 2.07 0.07745 0.02118
DOI: 10.7717/peerj.5285/table-2
Table 3:
Results of RankProd analysis-33 downregulated genes.
ID_REF Identifier gene.index RP/Rsum FC:(class1/class2) pfp P.value
221475_s_at RPL15 104 3.889 0.67 7.27E−120 1.24E−121
203535_at S100A9 103 4.437 0.71 1.20E−108 3.09E−110
209719_x_at SERPINB3 3 6.49 0.46 2.94E−79 1.00E−80
200074_s_at GUK1 55 7.095 0.45 4.17E−73 1.78E−74
209351_at KRT14 26 7.505 0.393 2.33E−69 1.20E−70
201097_s_at ARF4 11 8.908 0.32 2.43E−58 1.46E−59
217845_x_at HIGD1A 2 10.75 0.247 3.25E−47 2.22E−48
210413_x_at SERPINB4 7 10.93 0.27 2.59E−46 1.99E−47
210835_s_at CTBP2 77 12.68 0.17 2.39E−38 2.05E−39
209720_s_at SERPINB3 1 13.56 0.13 5.96E−35 5.60E−36
200761_s_at ARL6IP5 65 13.75 0.13 2.67E−34 2.74E−35
201619_at PRDX3 69 14.18 0.10 7.82E−33 8.69E−34
206276_at LY6D 47 14.37 0.15 3.27E−32 3.91E−33
201098_at COPB2 42 15.65 0.07 3.17E−28 4.07E−29
201653_at CNIH1 87 16.31 0.067 2.05E−26 2.80E−27
211906_s_at SERPINB4 8 16.42 0.003 3.80E−26 5.52E−27
202753_at PSMD6 56 17.6 0.01 3.54E−23 5.45E−24
202209_at LSM3 15 17.92 0.004 1.97E−22 3.20E−23
211023_at PDHB 9 17.95 0.018 2.14E−22 3.66E−23
221896_s_at HIGD1A 5 19.47 954.6 3.69E−19 6.62E−20
213164_at SLC5A3 63 22.18 787.9 2.11E−14 3.97E−15
201863_at FAM32A 106 22.74 767.2 1.42E−13 2.78E−14
209694_at PTS 66 24.92 688.4 1.09E−10 2.23E−11
218845_at DUSP22 10 26.22 675.4 3.21E−09 6.86E−10
203315_at NCK2 80 27.84 584.5 1.28E−07 2.84E−08
213357_at GTF2H5 114 28.09 612.8 2.09E−07 4.83E−08
217850_at SNORD19B 29 28.27 592.6 2.95E−07 7.06E−08
218283_at SS18L2 12 28.72 578.6 6.99E−07 1.73E−07
203282_at GBE1 14 31.25 526.7 5.77E−05 1.48E−05
212488_at COL5A1 68 32.86 426.3 0.0005592 0.000148
219389_at SUSD4 48 34.36 439.2 0.003389 0.000927
218115_at ASF1B 36 35.81 405 0.00481 0.004177
201692_at SIGMAR1 117 36.53 393 0.00775 0.008063
DOI: 10.7717/peerj.5285/table-3


Description of Tables 2 and 3: gene.index-index of genes; RP/Rsum-computed Rank Product statistics; FC: (class1/class2)-average expression levels’ computed fold change under two conditions (upregulated class and downregulated class); pfp-estimated false positive predictions value of the genes; and P-value-estimated p-value of each gene.

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.

Hierarchical clustering of the selected genes.

Figure 2: Hierarchical clustering of the selected genes.

(A) Hierarchical clustering of the 32 upregulated genes. (B) Hierarchical clustering of the 33 downregulated genes. In the figures, each gene is represented by the rows and the samples are represented by the columns. The dendogram at the side indicates the relation between the pattern of the gene expression while the top-dendogram indicates the relation between the samples used. The level of expression of the genes, relative to the mean of the gene across all samples, is indicated by the color key, with the green representing the higher expression of the genes. The color bar represents the cluster of the genes after they have been cut at 1.5 of the height of the tree so that the clustering of the genes is seen more clearly. The color of the bar is indicated by the color key similarly.
Table 4:
Ranking of each gene (upregulated) using proposed FS methods.
ID_REF Identifier HC PCC Relief-F SFS SVM-RFE CFS RF IG Score
206731_at CNKSR2 1 1 1 1 1 1 1 1 8
209243_s_at PEG3 1 0 1 1 1 1 1 1 7
210549_s_at CCL23 1 1 1 1 0 0 0 1 5
216108_at LOC105373738 1 1 1 0 0 1 0 1 5
220298_s_at SPATA6 1 1 1 1 1 0 1 1 7
212797_at SORT1 0 1 0 1 0 0 1 1 4
207996_s_at LDLRAD4 0 1 0 1 0 0 0 0 2
208762_at SUMO1 0 1 0 0 0 0 1 1 3
207257_at EPO 1 1 0 1 0 1 0 1 5
202035_s_at SFRP1 1 0 0 0 1 0 1 1 4
204672_s_at ANKRD6 0 0 1 1 1 0 1 1 5
220994_s_at STXBP6 0 0 1 1 1 0 1 1 5
213652_at PCSK5 1 0 0 0 1 1 0 0 3
213994_s_at SPON1 1 0 1 1 1 1 1 0 6
221606_s_at HMGN5 0 0 1 1 1 1 1 0 5
203440_at CDH2 1 0 1 1 0 1 0 1 5
207046_at HIST2H4B 1 0 0 0 0 0 0 0 1
211865_s_at FZR1 0 1 1 1 0 1 0 1 5
215401_at AU147698 1 0 1 0 1 1 0 1 5
213068_at DPT 1 0 0 0 1 1 1 5
214117_s_at BTD 1 1 0 1 1 1 1 1 7
34187_at AK026407 1 1 0 1 1 0 1 1 6
213611_at AQP5 1 0 0 1 1 1 1 0 5
209465_x_at PTN 0 0 1 0 1 1 0 0 3
208606_s_at WNT4 0 0 1 1 1 1 1 0 5
211814_s_at CCNE2 0 0 1 1 1 1 1 0 5
208790_s_at PTRF 0 0 0 0 0 1 0 0 1
200908_s_at RPLP2 1 1 1 1 1 0 1 1 7
219795_at SLC6A14 0 0 0 0 1 1 0 0 2
205730_s_at ABLIM3 1 0 1 0 1 1 1 0 5
202648_at TCF3 1 1 1 1 0 1 1 1 7
216267_s_at TMEM115 1 1 1 0 1 0 1 1 6
DOI: 10.7717/peerj.5285/table-4


The nine upregulated genes selected by most (≥6) of the methods are CNKSR2, PEG3, SPATA6, SPON1, BTD, AK026407 RPLP2, TCF3 and TMEM115.

Table 5:
Ranking of each gene (downregulated) using proposed FS methods.
ID_REF Identifier HC PCC Relief-F SFS SVM-RFE CFS RF IG Score
221475_s_at RPL15 1 1 1 0 1 1 1 1 7
203535_at S100A9 0 0 0 1 1 1 1 0 4
209719_x_at SERPINB3 0 0 0 0 1 1 0 0 2
200074_s_at GUK1 1 0 1 0 0 0 1 0 3
209351_at KRT14 0 0 1 0 0 1 0 0 2
201097_s_at ARF4 1 1 0 1 0 0 0 1 4
217845_x_at HIGD1A 1 1 1 1 1 0 1 1 7
210413_x_at SERPINB4 0 1 0 0 1 1 1 0 4
210835_s_at CTBP2 1 0 0 0 0 1 1 1 4
209720_s_at SERPINB3 0 1 0 1 0 0 1 0 3
200761_s_at ARL6IP5 1 1 1 1 1 0 1 0 6
201619_at PRDX3 1 1 1 1 1 1 0 1 7
206276_at LY6D 0 0 1 0 1 1 0 0 3
201098_at COPB2 1 1 1 1 1 1 1 7
201653_at CNIH1 1 1 1 1 0 1 1 1 7
211906_s_at SERPINB4 0 1 0 1 0 1 0 0 3
202753_at PSMD6 1 1 0 1 1 0 1 1 6
202209_at LSM3 1 1 1 1 1 1 0 0 6
211023_at PDHB 1 1 1 1 1 1 1 1 8
221896_s_at HIGD1A 1 1 1 1 1 0 1 0 6
213164_at SLC5A3 1 0 0 1 1 1 1 1 6
201863_at FAM32A 1 1 1 0 0 0 0 4
209694_at PTS 1 1 1 1 1 1 0 0 6
218845_at DUSP22 1 0 1 0 1 0 1 1 5
203315_at NCK2 1 0 0 1 1 0 1 1 5
213357_at GTF2H5 1 1 1 0 0 0 0 1 4
217850_at SNORD19B 1 1 0 1 0 1 0 1 5
218283_at SS18L2 1 1 0 0 0 0 0 1 3
203282_at GBE1 1 1 1 0 0 1 0 1 3
212488_at COL5A1 1 0 1 0 0 0 0 3
219389_at SUSD4 1 0 0 1 1 0 1 0 4
218115_at ASF1B 1 0 1 1 1 1 1 1 7
201692_at SIGMAR1 1 0 1 0 0 1 1 1 5
DOI: 10.7717/peerj.5285/table-5


The 12 downregulated genes selected by most (≥6) of the methods are RPL15, HIGD1A, ARL6IP5, PRDX3, COPB2, CNIH1, PSMD6, LSM3, PDHB, SLC5A3, PTS and ASF1B. As gene HIGD1A is selected twice with different ID REF, it was only considered once here.

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 gene sets were indicated in the colour nodes while the line connecting the nodes are the pathways that the enriched genes found in.

Table 6:
Results of GSEA analysis using the nine upregulated genes.
Identifier Highest normalized enrichment score (NES) FDR q-value Status
CNKSR2 1.00 0.859 Rejected
PEG3 1.33 0.000 Pass
SPATA6 1.00 0.823 Rejected
SPON1 1.60 0.053 Pass
BTD 1.50 0.026 Pass
RPLP2 1.40 0.000 Pass
DOI: 10.7717/peerj.5285/table-6


There was no matching of gene AK026407, TCF3 and TMEM115 with the MSigDB.

Table 7:
Results of GSEA analysis using the 12 downregulated genes.
Identifier Highest normalized enrichment score (ES) FDR q-value Status
RPL15 1.21 0.779 Rejected
HIGD1A 1.25 1.00 Rejected
PRDX3 1.41 0.020 Pass
COPB2 1.41 0.020 Pass
PSMD6 1.44 0.560 Rejected
LSM3 1.41 0.020 Pass
PDHB 1.33 0.839 Rejected
SLC5A3 1.50 0.000 Pass
ASF1B 1.43 0.000 Pass
DOI: 10.7717/peerj.5285/table-7


There was no matching of gene ARL6IP5, CNIH1 and PTS with the MSigDB.

Pathway network of the upregulated genes.

Figure 3: Pathway network of the upregulated genes.

Pathway network of the downregulated genes.

Figure 4: Pathway network of the downregulated genes.

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 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.

Functional pathway analysis- GSEA

According to the annotation of the 21 genes (nine upregulated genes and 12 downregulated genes) to the C4 cancerous gene set collection in MSigDB GSEA, it showed that six upregulated genes (CNKSR2, PEG3, SPATA6, SPON1, BTD and TCF3) and nine downregulated genes (RPL15, HIGD1A, PRDX3, COPB2, PSMD6, LSM3, PDHB, SLC5A3 and ASF1B) were enriched with the computational cancerous gene set. Particularly, four upregulated genes (SPON1, PEG3, RPLP2 and BTD) and five downregulate genes (PRDX3, COPB2, LSM3, SLCA3 and ASF1B) showed significant enrichment to the cancerous modules and cancerous genes neighbourhood gene sets with the FDR <0.25. This indicates 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 8:
Function, related diseases and description of the most significant (upregulated).
Genes Full name Function Related diseases and description
PEG3 Paternally Expressed Gene 3 • protein coding gene. • Hypermethylation of PEG3 is linked to epigenetics mechanism and HPV infection in cervical intraepithelial neoplasia (CIN) and Invasive Cervical Cancer. (Nye et al., 2013)
• Cell proliferation and p53-mediated apoptosis. (Yamaguchi et al., 2002)
• Involvement in the glioma and ovarian cancers (Jiang et al., 2010)
SPON1 F-spondin 1 • extracellular matrix organization regulation, interaction between cells and axon guidance (Burstyn-Cohen et al., 1999) • Showed extreme expression activities in the identification of colorectal biomarkers (Pagnotta et al., 2013)
• Down regulation activity in the cervical cell carcinoma (Zhai et al., 2010)
BTD bitonidase • gene expression, proliferation and differentiation of cells, gene signaling (Faith & Abraham, 2013) • gene expression signature that marks pelvic lymph node metastasis (PLNM) in cervical carcinoma (Huang et al., 2011)
• Involvement in thyroid and lung cancers (So et al., 2012; Scheerger & Zempleni, 2003)
• Potential serological marker in breast cancer plasma (Kang et al., 2010)
RPLP2 Ribosomal Protein Lateral Stalk Subunit P2 • encodes a 60s subunit ribosomal protein. • 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)
DOI: 10.7717/peerj.5285/table-8

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.

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).
• Showed immunopositivity to be significantly higher in cervical cancer cells (Hu, Gao & Li, 2013)
COPB2 Coatomer Protein Complex Subunit Beta 2 • protein coding gene.
• 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).
DOI: 10.7717/peerj.5285/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 potential-gene expression cervical cancer signatures.


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.

Supplemental Information

Result of LIMMA analysis

DOI: 10.7717/peerj.5285/supp-1

Gene expression dataset of combined dataset

DOI: 10.7717/peerj.5285/supp-2

GSEA analysis result

DOI: 10.7717/peerj.5285/supp-3
25 Citations   Views   Downloads