An integrative exploratory analysis of –omics data from the ICGC cancer genomes lung adenocarcinoma study

Abstract

It is known that all agents that cause cancer (carcinogens) also cause a change in the DNA sequence. In order to identify such often subtle changes, we attempt to integrate multiple molecular profile data sets released by the International Cancer Genome Consortium (ICGC). The list of data sets includes matched gene and microRNA expression profiles, somatic copy number variation, DNA methylation, and protein expression profiles for lung adenocarcinoma patients receiving treatments. We consider both unsupervised and supervised learning techniques (clustering and penalized regression) to identify interesting molecular markers corresponding to each type of –omics profiles that can differentiate patients. Associations between important markers of 2 types have been studied. An adaptive ensemble binary regression model has been presented that uses the entirety of available –omics profiles leading to a more accurate clinical prognosis for the patients in the given sample. This integrated study provides a more comprehensive picture of lung adenocarcinoma.

Keywords:

• chromosome
• cluster
• gene
• ICGC
• miRNA
• protein

Introduction

According to Stratton et al,¹ identification of changes in the DNA sequence of a genome is central to the study of all cancers. There are many ways of observing biological data related to the DNA sequence, including measuring gene expressions, microRNA (abbreviated miRNA hereafter) expressions, protein expressions, somatic copy number variation, and DNA methylation profiles for individual subjects. Each type of data can give insights into related disease processes; finding ways to integrate statistical analyses of data sets of multiple types has the potential to identify connections between important genes, miRNAs, proteins, chromosomal segments, and methylation patterns associated with the disease process rather than merely identifying those which are associated with the disease process from each individual molecular profile.

In recent years, there have been some attempts to understand the disease process comprehensively by integrating various molecular profiles of patients (Kristensen et al.²). In this paper, we consider several methods of analyzing matched data on genes, miRNAs, proteins, and copy number variation for the CAMDA 2014 lung adenocarcinoma challenge data provided by ICGC. We also explore the methylation patterns for this data.

We begin with separate analyses of data sets of different molecular types, where we cluster subjects with similar profiles (section- Exploratory Cluster Analysis). Subsequently, the clustering results are compared 2 at a time, along with clinical outcomes, using an overlapping proportion measure. In addition, penalized logistic regression models based on the elastic net regularization (Zou et al.³) are used to detect statistically important variables (such as genes or miRNAs) for each molecular profile type (section- Prediction of Clinical Outcomes) for disease prognosis that is characterized by 2 clinical outcomes following treatment: disease progression and disease remission. An optimal model with integrated data is found in terms of accuracy rates (section- Integrated Penalized Regression) followed by ensembling. Finally, a correlation analysis is performed to identify associations between the significant variables detected by the regression analyses leading to a more complete picture of the disease process. Following the format of this journal, the results are presented in the next section. All the above mentioned underlying procedures are described in the Materials and Methods section of the manuscript. The paper ends with a Conclusion section.

Results

As part of an exploratory analysis, we first report the plots of normalized expression values for genes, miRNAs, and proteins (Fig. 1A–C respectively) for each data set (molecular profiles) to obtain overall summaries of the data sets and to identify features that stand out (perhaps indicating strong signals). The two genes with the highest mean expression values are SFTPC and SFTPA1 (Fig. 1A); SFTPC regulation is known to be associated with lung diseases (see refs. ^4-6) and SFTPA1 is associated with newborn respiratory distress syndrome and lung sarcoma (see refs. ^7-13). The miRNA with the highest mean expression value is hsa-mir-21 (Fig. 1B); this is well-known as an oncomir associated with a wide variety of cancers (see refs. ^14-19). The protein with the highest mean expression value is FN1 (Fig. 1C). Its expression is up-regulated in non-small cell lung carcinoma (NSCLC) and it may promote lung tumor growth/survival and resistance to therapy (Han et al.²⁰). A plot of the average segment mean by chromosome for the progression and complete remission groups is shown in Figure 1D. In both groups, chromosomes 1, 5, 14, 17, and 20 have higher average segment means while chromosomes 9, 13 and 18 have lower averages. Oxnard et al.²¹ identified that chromosomes 1 and 17 are the locations of some potential oncogenes for non-small cell lung cancer. More importantly, the average segment means of chromosomes 9, 22 and X are noticeably different in the 2 groups of subjects with 2 different clinical outcomes or 2 different stages of the disease process (disease progression and the remission from the disease). According to Liu et al.,²² chromosome 9 is the location of the LHX6 gene, which acts as a potential tumor-suppressor with epigenetic inactivation in lung cancer. This chromosome has significantly lower average segment mean in the disease progression group compared to the complete remission group. We paid attention to these results for more in depth analysis later. The segment means for Chromosome 22 appear to be positively (negatively) associated with expressions for genes (proteins) that are predictive of the clinical outcome (see, e.g., Fig. 2).

Figure 1. Plots of various molecular profiles averaged over all available samples used in our analysis. Only some form of normalized molecular profiles were provided in the source data. (A) Plot of average expression values for different genes. (B) Plot of average expression values for different miRNAs. (C) Plot of average expression values for different proteins. (D) Plot of average segment means for various chromosomes.

Figure 2. Diagram showing connections between some predictive genes, miRNAs, proteins, and chromosomes. All connected pairs have statistically significant correlations with P values< 0.05, except for one.

It is now well known that DNA methylation is essential for proper mammalian development and plays a role in many important activities such as X-chromosome inactivation. Hence we attempt to find the average methylation patterns in every chromosome, especially, for chromosomes 9, 22 and X as they seem to have different segment means in the 2 groups of patients. We plot the methylation patterns for several portions of the data (not shown here). Additionally, we find that on chromosomes 9, 22 and X, there is a statistically significant difference between the average methylation patterns in the subjects in complete remission and progression groups, based on a 2-sample t-test (with test statistic values 7.79, 3.80,−18.07 and p-values 6.7e-15, 0.0001 and <2.2e-16, respectively). Tsou et al.²³ suggested that chromosomal instability caused by hyper or hypo methylation could potentially become a causal epigenetic factor in cancer. Results of this analysis suggest that there may be a role of methylation in altering the copy number variation in chromosomes 9, 22 and X and may play an important role in NSCLC cancer.

Clustering and co-clustering results

For each data set containing a specific molecular profile, internal validation, described in Brock et al.²⁴ and Datta et al.,²⁵ is used to select the appropriate clustering algorithm and the optimal number of clusters. Plots are made of each internal validation measure (see Eqns. 1-3; Fig. 3). Although one can visually inspect such plots to identify an optimal clustering algorithm and the number of clusters, such a task becomes increasingly difficult when the number of competing clustering algorithms is large and so is the number of validation measures. Thus we use a stochastic optimization and rank aggregation based method developed earlier by Pihur et al.²⁶ and Pihur et al.²⁷ to find the best clustering algorithm and the optimal size of clusters for this data.

Figure 3. Internal validation measures for 3 different clustering algorithms using proteins. The optimal cluster size is 2 and the k-means clustering optimizes each of the validation measures for the protein data.

Complete linkage is used for gene and miRNA profiles, but for protein and copy number profiles, Ward's linkage is used for clustering the data. This is done to maintain relative uniformity between the sizes of obtained clusters. We find that hierarchical clustering with 2 clusters for each of the gene, miRNA and copy number profiles and k-means clustering with 2 clusters for the protein profile provide the best clustering of the subjects. A representative dendrogram based on hierarchical clustering is shown in Figure 4. Interestingly, the optimal cluster number matches with the number of patient groups as characterized by the 2 possible clinical outcomes.

Figure 4. Dendrogram of hierarchical clustering using miRNA profiles. The divisions for the 2 main clusters are indicated by dashed lines; the two main clusters are indicated by 2 different colors.

In order to examine the consistency of the subject clustering using pairs of molecular profiles, we calculate the overlap proportions of the cluster profiles. The overlap proportions of subjects between each pair of cluster profiles (Eqn. 4(4) ${\displaystyle R_{j,k}=\frac{1}{n_{j,k}}{\displaystyle \sum _{i=1}^{n_{j,k}}\frac{\left|C_{(i)}^j\cap C_{(i)}^k\right|}{\left|(C_{(i)}^j\cup C_{(i)}^k)\cap I_{j,k}\right|}}}$(4) ) are shown in Table 1. Empirically (as well as mathematically), it can be seen that the overlap proportion is expected to be roughly 0.33 if the group assignments are made randomly. As all the proportions in Table 1 are larger than 0.33, there is some amount of consistency in the subject clustering results based on various –omics profiles data. Moreover, the overlap proportions between clustering based on gene and clinical data, as well as miRNA and clinical data, are quite substantial. These indicate that there are noticeable differences in the –omics profiles between the 2 groups of patients who would end up in remission versus those for whom the cancer would progress.

Table 1. Overlap proportions between the subject cluster profiles combining the results of clustering of 2 molcular profiles at a time

	Gene	miRNA	Protein	Chromosome
miRNA	0.52
Protein	0.44	0.39
Chromosome Clinical data	0.370.47	0.380.44	0.360.39	0.36

Penalized regression results

Next we investigate the predictive performances of the molecular profiles. The elastic net regularization for the penalized logistic regression model (see Eqn. 5) is used to predict the clinical outcomes/disease status (progression or remission) of the subjects. Each data set contains the molecular profiles of the subjects, along with their age and gender. Overall, 19 subjects have missing age and 104 subjects have missing disease status in the clinical data set; these subjects are excluded from regression analysis which adjusts for age and gender.

The optimal values of the elastic net tuning parameters (α, λ) (see equation Eqn. 6(6) $\begin{array}{l}-k^{-1}\text{log likelihood}+\lambda \times \text{penalty, where penalty}\\ =0.5\left(1-\alpha \right)\Vert \beta {\Vert }_2^2+\alpha \Vert \beta {\Vert }_1\end{array}$(6) ) are determined to be (1, .087), (.7, .102), (.5, .145), and (.3, .145) for individual regression models based on the gene expression, miRNA expression, protein expression, and copy number variation data sets, respectively. These models were used for variable selection for a correlation analysis (see Section 2.3). Next, we ran additional elastic net regressions where we integrate all possible combinations of the individual molecular profiles as covariates for predicting the clinical outcomes and then report the accuracy of prediction of clinical status of the subjects (Table 2). The accuracy rates for all models range from 64.3% to 74.3%. Finally, we use an adaptive ensemble binary regression model, described in Datta et al.²⁸ and Shah et al.,²⁹ where all individual penalized logistic regressions based on each molecular profile are considered as its components (see the methods section for details) leading to an overall accuracy rate of 75.7%.

Table 2. Accuracy rates for models with different combinations –omics profiles. Here G, M, P and C stand for gene, miRNA, protein and copy number variation, respectively

	Covariates included
	G	M	P	C
Accuracy Rate	0.643	0.686	0.686	0.743
	Covariates included
	G, C	M, P	M, C	P, C	G, M	G, P
Accuracy Rate	0.643	0.657	0.700	0.714	0.657	0.643
	Covariates included
	G, P, C	G, M, C	G, M, P	M, P, C
Accuracy Rate	0.643	0.657	0.657	0.657
	Covariates included
	All	Ensemble Model
Accuracy Rate	0.657	0.757

Correlation analysis results

The Spearman's rank correlation coefficients (see Eqn. 8) are computed between each pair among the genes, miRNAs, and proteins deemed most important by the penalized logistic regression analyses, and with the segment means for each chromosome. A diagram illustrating all interesting connections corresponding to strong (positive or negative) correlation is given in Figure 2. Each pair of connected items is statistically significant (P-value <0.05) except for the CDKN1B –chromosome 22 pair, which is borderline significant (P-value = 0.06).

The miRNA hsa-mir-223 has significant positive correlations with 3 genes TMOD3, TBC1D1, and PLEC. Each of these genes has a significant negative correlation with the protein CDKN1B while the genes TMOD3 and TBC1D1 have significant negative correlations with the protein ESR1 (part of estrogen signaling pathway). p27 is an enzyme inhibitor which is encoded by the CDKN1B gene. Past research suggests that p27 is regulated by miRNAs and is associated with the cell cycle pathway. This is associated with multiple carcinomas. According to Blain et al.,³⁰ p27 is highly correlated with the prognosis of patients and hence has been explored as a potential target for cancer therapy. For example, low levels of p27 may indicate that a cancer is not that aggressive. This also has a significant negative association with chromosome 22. Baik et al.³¹ identified that ESR1 is a part of estrogen signaling pathway which is being considered as a novel therapeutic target for NSCLC in women. The gene TMOD3 has been associated with estrogen excess and transitional cell carcinoma as noted by Pawlak et al.³² It is connected with the protein ESR1 which has been linked to estrogen resistance syndrome and uterine disease (ref. ³³) and NSCLC as stated before. TMOD3 is also connected with the protein CDKN1B which has been found to be connected with breast and prostate cancer as noted by Canbay et al.³⁴ and Chang et al.³⁵ From Wiesner et al.,³⁶ we know that loss of CDKN1B has also been frequently detected in tumors for blastic plasmacytoid dendritic cell neoplasms. The PLEC gene is also highly correlated with the protein CDKN1B. According to Charlesworth et al.,³⁷ it can have mutations which cause autosomal recessive forms of epidermolysis bullosa simplex with muscular dystrophy or pyloric atresia. Additionally, PLEC is known to be positively correlated with the Epidermal growth factor receptor (EGFR) signaling pathway. Rosell et al.³⁸ found that EGFR mutation related lung cancers are associated with lung cancers in women and never-smokers. Genetic variation in the TBC1D1 gene causes obesity in women. Obesity is typically associated with decreased insulin sensitivity and elevated circulating concentrations of glucose and insulin. Estrogen signaling plays a role in this process.

Materials and Methods

We consider pre-processed and normalized data from the International Cancer Genomic Consortium (ICGC) on gene expression, miRNA expression, protein expression, somatic copy number variation, and DNA methylation profiles for subjects with lung adenocarcinoma. As the number of genes is too large, a variance based filtration of the gene profile is used. Through this technique, we retain the top 75% of the genes based on the variances of their expression values. Two subjects (DO23996 and DO25020) with different sets of genes are dropped leading to a data set with 15916 genes on 131 subjects. Also, there are expression values for 709 miRNAs on 379 subjects (after removing the miRNAs which had all zero values for normalized expression levels) and expression values for 139 proteins on 237 subjects. In addition, there are a few subjects with missing chromosomal segment means in the copy number variation data, and after removing these subjects, segment means for each of the 24 chromosomes (taking X and Y as different chromosomes) are available for 383 subjects. Finally, methylation data is available on all chromosomes, and clinical data for disease status (progression or complete remission), age, and gender are also available. In summary, there are 2 groups of patients, namely, progression and complete remission, with multiple molecular profiles.

Exploratory cluster analysis

In order to find natural groupings in the data we use cluster analysis. Cluster analysis groups subjects in terms of the similarity of their molecular profiles. As different clustering algorithms group the data differently, we use 3 different clustering algorithms, hierarchical (with Ward and/or complete linkage), k-means, and PAM, and obtained the best clustering algorithm for each of these molecular profiles in the following manner. We evaluate the performances of these clustering algorithms using 3 different measures for internal validation, namely, connectivity (see refs. ^24,39), silhouette width (see refs. ^24,40), and the Dunn index (see refs. ^24,41). Then a stochastic optimization technique named rank aggregation^26,27 is used to combine the results of above internal validation measures to determine the optimal clustering algorithm as well as an appropriate cluster size. The definitions of the internal validation measures considered by us are given below.

Connectivity: For the i^th observation in the data set, let its ${\displaystyle j^{th}}$ nearest neighbor be denoted by ${\displaystyle n{n}_{i(j)}}$ and let ${\displaystyle x_{in{n}_{i(j)}}}$ be defined as 0 if ${\displaystyle i}$ and ${\displaystyle n{n}_{i(j)}}$ are in the same cluster and ${\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$j$}\right.}$ otherwise. Here, the nearest neighbor to the ${\displaystyle i^{th}}$ observation is that observation (other than the ${\displaystyle i^{th}}$ observation itself) which has the smallest distance from the ${\displaystyle i^{th}}$ observation. This is the interpretation for ${\displaystyle j\hspace{0.17em}=\hspace{0.17em}1}$. Similarly, the ${\displaystyle j^{th}}$ nearest neighbor to an observation, for ${\displaystyle j\hspace{0.17em}>1}$, can also be defined. For a particular clustering partition of N observations into K disjoint clusters, say, ${\displaystyle C=\left\{C_1,C_2,\mathrm{\dots .},C_K\right\}}$, the connectivity is defined as(1) ${\displaystyle Conn(C)={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^L{x}_{in{n}_{i(j)}}}}}$(1) Here, L determines the number of neighbors that contribute to the connectivity measure. Note that, a neighbor to an observation can belong to a cluster different from that of the observation. Connectivity ranges from 0 to ${\displaystyle \infty }$; it should be minimized for obtaining the ‘optimal’ clustering partition.

Silhouette width: For any observation ${\displaystyle i}$, the silhouette value is defined as(2) ${\displaystyle S(i)=\frac{b_i-a_i}{\max (b_i,a_i)}}$(2) Here, $a_i=\frac{1}{n(C(i))}\sum _{j\in C(i)}dist(i,j)$ is the average distance between ${\displaystyle i}$ and all other observations in the same cluster, $b_i=\underset{C_K\in C\backslash C(i)}{\min }\sum _{j\in C_K}\frac{dist(i,j)}{n(C_K)}$ is the average distance between ${\displaystyle i}$ and the observations in the “nearest neighboring cluster,” where ${\displaystyle C(i)}$ is the cluster which contains the ${\displaystyle i^{th}}$ observation, ${\displaystyle dist(i,j)}$ is a distance (e.g., Euclidean, Manhattan) between the ${\displaystyle i^{th}}$ and the ${\displaystyle j^{th}}$ observations and n(C) is the number of observations in cluster C. The silhouette width is the average of the silhouette values for all observations and ranges from −1 to 1. The silhouette width should be maximized for obtaining the best clustering partition.

Dunn index: The Dunn index is computed as the ratio of the smallest distance between observations not in the same cluster to the largest intra-cluster distance which is given by(3) ${\displaystyle D(C)=\frac{\underset{C_K,C_l\in C,C_K\ne C_l}{\min }\left(\underset{i\in C_K,j\in C_l}{\min }dist(i,j)\right)}{\underset{C_m\in C}{\max }diam(C_m)}}$(3) Here, ${\displaystyle diam(C_m)}$ is the maximum distance between observations in cluster ${\displaystyle C_m}$. The Dunn index ranges from 0 to ${\displaystyle \infty }$ and needs be maximized for obtaining the ‘best’ clustering partition.

In order to maintain brevity of this paper, we do not include the details of the rank aggregation technique which are available in Pihur et al.²⁶ and Pihur et al.²⁷ For each molecular profile (e.g., gene, miRNA, protein, copy number variation), the rank aggregation technique is used to determine the optimal clustering algorithm among hierarchical, k-means and PAM based on the 3 internal validation measures described above. This technique also determines the optimal cluster size for each profile.

To measure the similarity between the optimal clustering partitions of data sets for different molecular profiles, the overlapping proportion(4) ${\displaystyle R_{j,k}=\frac{1}{n_{j,k}}{\displaystyle \sum _{i=1}^{n_{j,k}}\frac{\left|C_{(i)}^j\cap C_{(i)}^k\right|}{\left|(C_{(i)}^j\cup C_{(i)}^k)\cap I_{j,k}\right|}}}$(4) is computed for each pair of profiles. Here, ${\displaystyle C_{\left(i\right)}^j}$ is the set of subjects in the cluster containing the ith subject based on the jth profile, ${\displaystyle n_{j,k}=\left|I_{j,k}\right|}$ is the number of subjects common to the jth and kth data sets, and ${\displaystyle I_{j,k}}$ is the set of common subjects.

Prediction of clinical outcomes

For each data set, a penalized logistic regression model for predicting clinical outcomes for the disease status (progression or disease remission) is considered based on each molecular profile data such as expression values or chromosomal segment means of the subjects. Additionally, the ages and gender of the subjects are also used in the model. The model is given by(5) ${\displaystyle \text{logit}\left(p_{j,i}\right)={\beta }_{j,0}+{\beta }_{j,A}Ag{e}_{j,i}+{\beta }_{j,G}Gende{r}_{j,i}+{\beta }_{j,1}{X}_{j,i,1}+\dots +{\beta }_{j,m}{X}_{j,i,m}}$(5) for ${\displaystyle i=1, \dots , n_j}$ where, for the ${\displaystyle i}$th subject based on the ${\displaystyle j}$th profile data set, ${\displaystyle p_{j,i}}$ is the probability for progression of the disease, ${\displaystyle X_{j,i,k}}$ is the ${\displaystyle k}$th expression value or chromosomal segment mean, and the β's are the corresponding regression coefficients in the model. We fit the elastic net regression model, which estimates the parameter vector by minimizing(6) $\begin{array}{l}-k^{-1}\text{log likelihood}+\lambda \times \text{penalty, where penalty}\\ =0.5\left(1-\alpha \right)\Vert \beta {\Vert }_2^2+\alpha \Vert \beta {\Vert }_1\end{array}$(6) .

The R package glmnet, developed by Friedman et al.,⁴² is used to select the tuning parameters ${\displaystyle \alpha }$ and ${\displaystyle \lambda }$ by minimizing the cross validation error. In order to reduce the sensitivity of the model fitting and variable selection procedure, a bootstrap step is added (Efron et al.⁴³). Covariates with non-zero coefficients for each bootstrap sample are recorded. The bootstrap sampling process is repeated 1000 times, and the covariates are ranked based on how often they were selected by the elastic net model for the bootstrap samples. This provides us with a ranked list of variables corresponding to each molecular data type. The variables selected by this algorithm may represent the genes, miRNAs, and proteins that are most useful for prediction of the clinical outcome. We use the top 20 genes, miRNAs, and proteins for a correlation analysis to derive a more complete picture.

Integrated penalized regression

Next, the 4 data sets (gene expression, miRNA expression, protein expression and somatic copy number variation) are merged together and a penalized logistic regression model is fitted for predicting the clinical outcome of disease status based on age, gender, expression values and chromosomal segment means. The accuracy rate (see Eqn. 7(7) $\text{Accuracy Rate}=1-{\displaystyle \sum }_i\frac{\left|observe{d}_i-predicte{d}_i\right|}{total}$(7) ) for this model is computed and compared to that of the models with individual molecular profiles as covariates as described in the Penalized Regression Results section based on each of the 4 individual data sets.

Next, to understand which combinations of integrated molecular profiles provide more accurate prediction of the clinical outcomes of the subjects, we compute the accuracy rates for all the models involving all possible combinations of the individual data sets.(7) $\text{Accuracy Rate}=1-{\displaystyle \sum }_i\frac{\left|observe{d}_i-predicte{d}_i\right|}{total}$(7) where ‘observed’ and ‘predicted’ in the above display denotes the observed and predicted clinical outcomes, respectively, for every subject, and ‘i’ denotes a typical subject.

Finally, we build a super binary regression model to predict the clinical outcome where we integrate all data types in the following way. Bootstrap resamples of the original data are used as training data to fit penalized regression models based on individual –omics profiles. These models are then evaluated on the “out of bag” sample (this is a term used in the bagging literature to indicate part of the original sample that was not selected in a bootstrap sample) using multiple performance measures. An overall best predictive model is identified for each bootstrap replicate via rank aggregation of the performances of the models on various measures and stored for future prediction. Finally, the clinical outcome for a patient is predicted by running these individual best models and combining the predictions by majority voting.

In our case, 4 penalized logistic regression models with a common set of observations (70 subjects) based on each of the 4 data sets (gene expression, miRNA expression, protein expression and somatic copy number variation) are considered. Also, we consider 4 performance measures which are to be used (sensitivity, specificity, false discovery rate and false non-discovery rate). Data of one subject is removed from each of the 4 data sets and is kept aside for outcome prediction by the ensemble model. A bootstrap sample of size 69 is drawn (from the remaining observations) using simple random sampling with replacement. For each of the 4 data sets, a penalized logistic regression model is fitted using the bootstrap sample which is used to predict the binary outcomes of the out of bag (OOB) samples. Since the true responses (clinical outcome) for the OOB samples are known to us, we can compute the 4 performance measures (sensitivity, specificity, false discovery rate and false non-discovery rate) from the predicted outcomes. The four models are ranked according to their performances under each of the 4 performance measures and 4 ordered lists of models ${\displaystyle (L_1,L_2,L_3,L_4)}$, each of size 4, are obtained. These four ordered lists are aggregated using the weighted rank aggregation procedure^26,27 which determines the best performing model ${\displaystyle M_{(1)}^1}$. We repeat the procedure 100 times with new bootstrap samples and to obtain 100 best models. Using these 100 best models, we predict the clinical outcome of the left out observation. Each model cast a vote (either 0 or 1) regarding the clinical outcome and the final prediction is done on the basis of majority voting. The whole algorithm is repeated for each of the 70 observations. A schematic diagram describing the steps of computing the ensemble binary regression model is shown in Figure 5. The accuracy rate is calculated for the ensemble binary regression using the formula given in (7)(7) $\text{Accuracy Rate}=1-{\displaystyle \sum }_i\frac{\left|observe{d}_i-predicte{d}_i\right|}{total}$(7) .

Figure 5. A Schematic Diagram of the Ensemble Binary Regression (adapted from Datta et al., 2010, BMC Bioinformatics, 11, 427).

Correlation analysis

Spearman's rank correlation coefficients are computed for each pair of profile data sets using the formula(8) $${\rho }_{j,k}=1-\left(6\underset{i=1}{\overset{n_{j,k}}{{{{\displaystyle \sum }}^{\text{}}}^{\text{}}}}{d}_{j,k,i}^2\right)/\left(n_{j,k}^3-n_{j,k}\right),$$(8) where ${\displaystyle d_{j,k,i}}$ is the difference between the ranks of the ith subject that is common to the jth and kth profile data sets and ${\displaystyle n_{j,k}}$ denotes the number of such subjects (patients). Using the asymptotic t-distribution, approximate marginal p-values are computed for these correlation coefficients to identify statistically significant associations between important genes, miRNAs, proteins, and chromosomes that are predictive of the clinical outcomes.

Conclusion

Although there has been a vast number of research articles which examine how to connect individual molecular profiles with clinical outcomes, there are relatively few results on the integration of molecular profiles. However, in recent years there have been some attempts to integrate different molecular profiles to get a comprehensive knowledge of the disease process for a complex disease like cancer. In this paper, we use multiple methods of data integration to provide an integrated analysis of matched gene expression, miRNA expression, protein expression, and somatic copy number variation profiles and their effects on the clinical outcomes using NSCLC lung adenocarcinoma, CAMDA 2014 challenge data. It is to be noted that the data available were already normalized and we did not have access to the raw data. In general, due to ICGC policies regarding public disclosure of clinical information, very limited information was provided about the data. Additionally, the data lacks a control group and thus all these subjects have had cancer. So, we expect a lot of similarity in their molecular profiles. In spite of all these limitations, we are able to show that there are subtle but important and biologically meaningful differences between the –omics profiles of patients who would enter remission vs. those for whom the cancer would progress even after the treatment. The results are quite illuminating; however, further investigation is needed to claim the causation of the clinical outcome from such retrospective investigations. We are nevertheless optimistic that this data analysis pipeline can be used for integrating multiple molecular profiles data toward disease status prediction.

We find that the cluster analysis along with the rank aggregation technique is able to identify the number of groups in the data consistently across the molecular profiles which also matched with the number of clinical groups in the sample. Also, the co-clustering results show that not only individually but also collectively, the molecular profiles are somewhat consistent in terms of their ability to group the subjects into their respective clinical outcome groups. Secondly, the comparative results of penalized regression models using different individual molecular profiles, the combinations of multiple molecular profiles, and the ensemble binary regression (Table 2) show that optimally ensembling different individual regression models may provide the best predictive model for the clinical outcome; in our case it led to an observed accuracy rate of 75.7%. Thirdly, from our correlation analysis of top molecular features that are predictive of the clinical separation, we note that the estrogen signaling pathway and the epidermal growth factor receptor (EGFR) signaling pathways are 2 of the most differentiating pathways between the 2 groups of patients. Estrogen signaling is known to play a significant role in lung cancer, especially for women with estrogen therapy. EGFR pathway has long been established as therapeutic targets of lung cancer. It is also known that estrogen signaling interacts with EGFR in lung cancer. The existing results and the present analysis may mean that both of these pathways may play a role not only in terms of disease causation but also in the treatment outcome for NSCLC.

Overall, it is our hope that this paper makes a convincing case in favor of more integrated studies of –omics molecular profiles in the future since such studies yield a bigger biological picture than a single assay type which may in turn help us make better prognoses and provide better treatment regimes.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Acknowledgments

We thank an anonymous reviewer for many useful comments.

Funding

This research work is partially supported by NIH grant CA 170091−01A1 (Su. Datta).