Integration and comparison of different genomic data for outcome prediction in cancer

Background In cancer, large-scale technologies such as next-generation sequencing and microarrays have produced a wide number of genomic features such as DNA copy number alterations (CNA), mRNA expression (EXPR), microRNA expression (MIRNA), and DNA somatic mutations (MUT), among others. Several analyses of a specific type of these genomic data have generated many prognostic biomarkers in cancer. However, it is uncertain which of these data is more powerful and whether the best data-type is cancer-type dependent. Therefore, our purpose is to characterize the prognostic power of models obtained from different genomic data types, cancer types, and algorithms. For this, we compared the prognostic power using the concordance and prognostic index of models obtained from EXPR, MIRNA, CNA, MUT data and their integration for ovarian serous cystadenocarcinoma (OV), multiform glioblastoma (GBM), lung adenocarcinoma (LUAD), and breast cancer (BRCA) datasets from The Cancer Genome Atlas repository. We used three different algorithms for prognostic model selection based on constrained particle swarm optimization (CPSO), network feature selection (NFS), and least absolute shrinkage and selection operator (LASSO). Results The integration of the four genomic data produced models having slightly higher performance than any single genomic data. From the genomic data types, we observed better prediction using EXPR closely followed by MIRNA and CNA depending on the cancer type and method. We observed higher concordance index in BRCA, followed by LUAD, OV, and GBM. We observed very similar results between LASSO and CPSO but smaller values in NFS. Importantly, we observed that model predictions highly concur between algorithms but are highly discordant between data types, which seems to be dependent on the censoring rate of the dataset. Conclusions Gene expression (mRNA) generated higher performances, which is marginally improved when other type of genomic data is considered. The level of concordance in prognosis generated from different genomic data types seems to be dependent on censoring rate. Electronic supplementary material The online version of this article (doi:10.1186/s13040-015-0065-1) contains supplementary material, which is available to authorized users.


Background
Cancer is a public health problem worldwide due to its high prevalence and mortality rates [1]. In the year 2012 alone, there were 14.1 million new cases of cancer, from which 8.2 million resulted in death [2]. Moreover, projections estimate a 20 % and 40 % increase of cancer cases for the years 2020 and 2030, respectively relative to 2010. The cancers of breast and lung cancers are expected to remain within the top cancer diagnoses and leading causes of cancer-related death [3].
Identifying which source of genomic data, or combination, generates the most powerful prognostic biomarker could help to describe cancer etiology [16,19,20]. However, some studies have generated inconsistent results across cancers when evaluating distinct sources of genomic data for prognosis [19,20], probably because of the use of different algorithms. Thus, it is not clear which type of data is the best at predicting cancer prognosis or whether combinations of data types provide some improvement. For example, it has been shown that no significant improvement is obtained adding any genomic measurement once EXPR data and clinical covariates were included in the model [19] using principal components, partial least squares, and a penalization algorithm. On the other hand, a similar study showed that all clinical outcomes were better predicted when integrating multi-layers of genomic data [20] using a graph-based algorithm while others suggest that the clinical improvement of genomic data is limited in magnitude and on cancer types [21] using diverse classification algorithms.
Given the lack of concordance on methods and genomic data provides the best prognostic results and its utility, our purpose is to characterize the prognostic power of models obtained from different genomic data types, cancer types, and algorithms. For this, we tested the prognostic and concordance index of models obtained by three different algorithms from EXPR, MIRNA, CNA, MUT data and their integration for ovarian serous cystadenocarcinoma (OV), multiform glioblastoma (GBM), lung adenocarcinoma (LUAD), and breast cancer (BRCA) datasets from the TCGA repository. The algorithms used are based on very different properties to search for diverse solutions attempting to derive conclusions at certain independency of the algorithms. We used constrained particle swarm optimization (CPSO) [22], which explores combinations of features irrespectively of its biological connections, network feature selection (NFS) [23] that explores combination of features integrating protein-protein interaction information, and the least absolute shrinkage and selection operator (LASSO) [24] that explores penalized models.

Methods
The methodology is summarized in Fig. 1. Briefly, samples from the four cancer types that fulfill a specific inclusion criterion were selected. Features from each database and source were filtered. The resulting four databases were then merged into a metabase (MERGE) for comparisons with single-sourced databases. Predictive models were obtained for each database using three feature selection algorithms that generate a unique model. Finally, the performance of the models was evaluated using the concordance index (c-index) [25].

Database selection
The data used in this study was downloaded in April 2013 from the TCGA data portal (https://tcga-data.nci.nih.gov/tcga/) including level 2 (MUT) and level 3 (EXPR, MIRNA, and CNA) data. CNA was segmented by regions per sample using the GISTIC algorithm [26]. EXPR and MIRNA data were quantile-normalized before use. Using the TCGA-ID, a tag unique to each subject, only those subjects with available EXPR, MIRNA, CNA, and MUT data were used. The results published here are in whole or

Feature filtering
We filtered features to reduce complexity, to avoid the use of invariant information, and to balance the number of features from each source avoiding to remove predictive information. MIRNA and EXPR databases were filtered using a correlation and quantization strategy. First, features without absolute Spearman correlation coefficient larger than 0.6 were excluded. Second, to remove invariant genes, we split the data into five uniform segments and only those genes having counts in more than two segments were used. For CNA data, we used the 10 % probes having the most significant p-values using the univariate log-rank test from a Cox proportional hazard model splitting the linear predictor at the median. For MUT data, we used the 11.4 % of OV, 12.2 % of LUAD, 9.9 % of BRCA and 30.4 % of GBM, of the most frequently mutated genes. Using LUAD as an overall validation of the filtering procedure, we observed that using LASSO and all features in CNA, EXPR, MIRNA and MERGE, the results were 61, 77, 76, and 78 of concordance index respectively, which are very close to those observed after the filtering.

Metabase generation
A fifth dataset (MERGE) was constructed per cancer type by merging their corresponding EXPR, MIRNA, CNA, and MUT filtered databases. This allowed a direct comparison on which data source is best selected in the presence of other sources. Furthermore, the metabases permitted the identification of predictive models with features from different sources, and compare such compound models with singlesource models.

Feature selection algorithms
We used a multivariate Cox proportional hazard model for the three feature selection algorithms. Beta coefficients were calculated by optimizing either the log-likelihood (NFS and CPSO) or a penalized maximum likelihood function (LASSO) through several iterations using bootstrap (NFS and CPSO) or a 10-fold cross-validation (LASSO) scheme [27,28]. In the case of NFS and CPSO, only two-thirds of the population was used for training while the remaining was used to perform a blind test. Bootstrap consists in randomly sampling the population using a similar fraction per strata in the resampled sets [29]. For MUT databases, we relied on resubstitution because mutation data is sparse where only small number mutations are observed per gene, which may generate sets of training samples with no mutations at all.

Constrained particle swarm optimization (CPSO)
Particle swarm optimization algorithms are based on the biological behavior of swarms. Concisely, these algorithms create a swarm of particles with random positions and velocities. The positions represent parameters of the problem to solve. The particles will update their velocity and position depending on their performance, iteratively. The performance is a function that evaluates the particle position relatively to the swarm [22].
We have customized PSO (CPSO) to handle feature selection problems from large genomic datasets [30]. This algorithm uses a user-defined number of features, k, to generate efficiently a subset of features that is used as the survival model. We used k = 5 and 500 iterations. We ran the algorithm 1,000 times. Models generated contained between 8 and 10 genes. We used the model with the highest c-index estimated by bootstrapping.

Network feature selection (NFS)
Network Feature Selection (NFS) is based on the exploration of protein-protein interaction networks to select features resulting in more biologically coherent models [23]. NFS has recently been used to generate multi-cancer biomarkers [23]. Briefly, each feature is evaluated individually by the p-value of an univariate Cox proportional hazards model. Each gene is then considered as a survival model. Each model grows by considering all possible neighbors according to the interactions provided by a network. The top 5 % of these grown models having higher performance are selected to grow in the next iteration. This procedure is carried on until no model can be further grown, or until 10 iterations. The protein-protein interaction network used was downloaded from the human protein reference database (HPRD, http://www.hprd.org/). Genes having more than 1,000 connections are not allowed to grow (for example, the UBC gene). For MIRNA data, the interactions between miRNA and mRNA were considered as surrogate interactions for the network, where the mRNA was replaced by the miRNA that regulates it. In order to identify the targets of each miRNA and create the miRNA/protein-protein interaction network, the miRTarDatabase (http://mirtarbase.mbc.nctu.edu.tw/) was used. In the MERGE dataset, the gene/protein connections were used irrespective of the data type.

Least absolute shrinkage and selection operator (LASSO)
LASSO is a well-known widely used feature selection algorithm, particularly when the number of samples is considerably smaller than the number of features. This algorithm performs a coefficient penalization in which only well-associated features emerge [28]. The best model containing around 10 features was used.

Performance evaluation
Models were evaluated and compared using the concordance index (c-index) and the p-value of the log-rank. The c-index was used to assess the prediction power of the survival model [25,31]. The log-rank test was used to determine whether low-and highrisk groups were significantly different from each other [25,31]. These statistics were estimated using the blind test subset for the models generated with CPSO and NFS, or using re-substitution for the models generated with LASSO. To compare the agreement of prognostic prediction of two models, we used the Cohen's kappa statistic in R implemented within the package fmsb [32]. For this, we split the prognostic index by the median. The prognostic index is the linear predictor of the exponential function in the Cox model [27].

Results
We used OV, LUAD, BRCA, and GBM datasets that had at least 100 subjects with EXPR, MIRNA, CNA, and MUT data in the TCGA repository at the time of accession.
A brief description of the technologies and clinical and demographic information is included in Additional file 1: Table S1 and Additional file 2: Table S2. The number of features of each dataset before and after filtering is detailed in Table 1.
The results of the c-index and the log-rank test of all cancer types, data types, and algorithms are shown in Table 2. From the genomic data types, we observed better prediction in EXPR closely followed by MIRNA and CNA depending on the cancer type and method (Figs. 2 and 3). In our tests, mutation data generated poor predictions. In average, the results of the MERGE dataset were marginally more predictive that any of the other data types (Figs. 2 and 3). Within the MERGE dataset, we explored which of the dataset was more important. The Table 3 shows the number of features per data type used by the best model in the MERGE database. The results further support that EXPR is the preferred data (54 % of the features) when all other data is present. Surprisingly, EXPR was followed by CNA (27 %) and then MUT (14 %) while MIRNA data was almost not used (6 %).
We observed higher predictions in BRCA, followed by LUAD, OV, and GBM having an average c-index around 0.82, 0.71, 0.63, and 0.60, respectively. These comparisons agree with recent results of multi-cancer gene expression biomarkers [23]. In BRCA and OV, CNA data were more predictive than MIRNA. In LUAD and GBM, the c-index of MIRNA was higher than CNA and comparable with EXPR data. Although the results of MUT were poor, in BRCA and LUAD the predictions were higher than in OV and GBM even though we used more genes in those cancer types.
We observed similar c-index values between LASSO and CPSO but smaller c-index values in NFS (Fig. 3). The MERGE data was more predictive in LASSO and NFS but not in CPSO where EXPR was the best. CNA was clearly more predictive in CPSO than in LASSO and NFS (Fig. 3).
We also compared whether the predictions made by models concur. We used the Kappa statistic that measures the level of concordance of two predictors. Values of Kappa close to 0 correspond to random agreements whereas values close to 1 represent perfect agreement. The results show that MIRNA, CNA, and EXPR models have acceptable agreement in LUAD, OV, and GBM irrespective of the method (Fig. 4). In BRCA, we found agreement in CNA models and partially in MIRNA. In addition, MIRNA slightly agrees with CNA in LUAD and with EXPR in GBM. In general, however, the predictions made by different types of data disagree.

Discussion
Our objective was to compare and characterize the prognostic level of different genomic data sources in cancer. For this, we analyzed four important cancer types (BRCA,   Indicate models whose Kaplan-Meier curves were statistically different at 0.05, 0.01, and 0.001 level respectively using the log-rank test. For this, the population was split by the median using the prognostic index (linear predictor of the Cox model). "-" indicates that no models were generated Overall, MERGE data was the most predictive across the four cancer types (Fig. 2) and the three algorithms (Fig. 3). This result is sensible because MERGE contained all other data types. Nevertheless, in some cases MERGE was not the best. This was the case in CPSO whose performance could be influenced by the increased number of features.
From the genomic data types (EXPR, MIRNA, CNA, MUT), the best performance was obtained with EXPR ( Figs. 2 and 3). The gene expression is the result of complex   Table 2). CNA represents changes in DNA, which are presumably less dynamic than EXPR. Nevertheless, the CNA performance was surprisingly comparable to EXPR suggesting that a considerable component of the survival is dictated by CNA. The performance of MUT data was poor when compared to EXPR, MIRNA, and CNA. Some issues are known relative to this lack of prediction. First, mutation frequencies per gene are generally low suggesting that mutation data is highly disperse [21]. Second, the combination of sparseness and binary data (mutated or not mutated) may generate difficulties in the Cox model fitting. Third, the reports of mutation frequencies do not commonly find associations with survival [33][34][35][36]. We did not observe big differences in performance relative to the algorithm used. CPSO seems to show consistent and highly competitive results, but LASSO seems to report slightly higher results while NFS seems to produce lower performance.
The prognostic values provided by different methods for the same data type were remarkably high suggesting that the algorithm used is a minor source of differences (Fig. 4). However, the observation of the lack of similitude in risk prediction between different data types was surprising and an important result of our study (Fig. 4). It is known that the precision of the prognostic values is highly influenced by the proportion of censoring [37]. We observed higher similitudes between the prognostic values generated by different data types in GBM where the proportion of censoring is the lowest (21.3 %) and lower similitudes in BRCA where the proportion of censoring is the highest (84.4 %). We also observed that the prediction in BRCA is high (around 0.8 of c-index) while in GBM is low (around 0.6 of c-index). The c-index measures how well the model fit the censoring data while kappa measures the consistency of two predictions. We showed that these properties seem to be highly influenced by the proportion of events. More research is needed to determine the lack of consistency.
We used four cancer types, three algorithms, and five data types. There may be some level of interaction between these three components. For instance, the performance of MIRNA data was higher in LUAD and the performance of NFS was generally lower using CNA data. We did not study thoroughly the possible parameters combinations within each algorithm, nor many potential schemes of data type filtering and processing. However, our results suggest some tendencies and the results should be similar to other cancer types and algorithms in similar circumstances to those tested here.

Conclusions
The integration of genomic data produced survival models were marginally higher in performance than those from single genomic data, specially those of mRNA. From the genomic data, the mRNA gene expression generated the highest predictive models and were preferred in models that integrate the four types of genomic data. CNA and miRNA data followed mRNA in performance while mutation data poorly predicted survival. The risk prediction of survival models of different types of data disagrees and the level of agreement seems to be related to the censoring rate.