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

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

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.

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

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 potential-gene expression cervical cancer signatures.