Imputation of Gene Expression Data in Blood Cancer and Its Significance in Inferring Biological Pathways

Purpose: Gene expression data generated from microarray technology is often analyzed for disease diagnostics and treatment. However, this data suffers with missing values that may lead to inaccurate findings. Since data capture is expensive, time consuming, and is required to be collected from subjects, it is worthwhile to recover missing values instead of re-collecting the data. In this paper, a novel but simple method, namely, DSNN (Doubly Sparse DCT domain with Nuclear Norm minimization) has been proposed for imputing missing values in microarray data. Extensive experiments including pathway enrichment have been carried out on four blood cancer dataset to validate the method as well as to establish the significance of imputation.

Methods: A new method, namely, DSNN, was proposed for missing value imputation on gene expression data. Method was validated on four dataset, CLL, AML, MM (Spanish data), and MM (Indian data). All the dataset were downloaded from GEO repository. Missing values were introduced in the original data from 10 to 90% in steps of 10% because method validation requires ground truth. Quantitative results on normalized mean square error (NMSE) between the ground truth and imputed data were computed. To further validate and establish the significance of the proposed imputation method, two experiments were carried out on the data imputed with the proposed method, data imputed with the state-of-art methods, and data with missing values. In the first experiment, classification of normal vs. cancer subjects was carried out. In the second experiment, biological significance of imputation was ascertained by identifying top candidate tumor drivers using the existing state-of-the-art SPARROW algorithm, followed by gene list enrichment analysis on top candidate drivers.

Results: Quantitative NMSE results of the DSNN method were compared with three state-of-the-art imputation methods. DSNN method was observed to perform better compared to these other methods both at high as well as low observable data. Experiment-1 demonstrated superior results on classification with imputation compared to that performed on missing data matrix as well as compared to classification on imputed data with existing methods. In experiment-2, cancer affected pathways were discovered with higher significance in the data imputed with the proposed method compared to those discovered with the missing data matrix.

Conclusion: Missing value problem in microarray data is a serious problem and can adversely influence downstream analysis. A novel method, namely, DSNN is proposed for missing value imputation. The method is validated quantitatively on the application of classification and biologically by performing pathway enrichment analysis.

1. Introduction

1.1. Motivation

High dimensional gene expression data helps in determining gene-to-gene interactions in different biological pathways. Molecular techniques such as “Microarrays” are used to measure expression levels of genes. Different machine learning algorithms and statistical methods are applied to gene expression data to extract relevant information for applications such as disease diagnosis and classification of clinical subtypes. These analyses assist in developing effective drugs for specific diseases because treatment procedures differ from disease to disease and case to case. For example, in blood cancer, the drugs are required to be targeted toward the specific type of cancer. Similarly, in personalized therapy, the response of a subject to a particular drug can be captured and its correlation with the mutation profiles of the subject can be examined to design targeted medicine. Thus, the therapy can be customized according to the genetic built of the subject.

A persistent problem associated with microarray dataset is the presence of varying number of missing values in the data that may arise owing to poor slide quality (dusty or scratchy), poor image quality, or insufficient resolution (1). Subsequent downstream analysis on incomplete gene expression matrices may be highly inaccurate. One of the ways of dealing with the problem of missing values is to capture microarray data again but it does not guarantee complete data matrix. Moreover, the entire process is expensive and time consuming. An alternate solution to this problem is to remove the genes containing missing values from the analysis. However, this can result in loss of information and may lead to inaccurate findings on driver genes and/or altered biological pathways. Therefore, it is worthwhile to apply advanced computational methods for the imputation of missing values in microarray data prior to any analysis.

1.2. Background

Numerous methods have been developed in the recent times for imputation of gene expression data. These can be broadly categorized into four classes: hybrid methods, local methods, global methods, and knowledge assisted methods (Table 1). Some of the early methods developed to account for the missing values are ZEROimpute, ROWimpute and COLimpute (18). In ZEROimpute, missing values are replaced with zeros. In ROWimpute and COLimpute, missing values are replaced with the averaged values of the observed entries of the corresponding rows or columns. These methods do not take into consideration the correlation present among genes and therefore, do not perform optimally. Gene expression matrix is highly correlated. Therefore, it is important to consider correlation among genes. Several methods exist in literature based on correlation among genes. These are categorized into local and global approaches based on the type of correlation utilized by them. As shown in Table 1, local approaches impute missing values by considering the group of genes that show high correlation with the gene containing missing values. Such methods perform optimally when the data is heterogeneous. k nearest-neighbor imputation (KNNimpute) (2) is one of the earliest local approach method to impute missing value. It first estimates k nearest group of genes that are similar to the missing target gene, followed by averaging of these genes to impute the missing value of the target gene. SKNNimpute (Sequential KNNimpute) (3) and IKNNimpute (iterative KNNimpute) (4) are variations of KNNimpute. Gaussian mixture clustering imputation (GMCimpute) (5), least square imputation (LSimpute) (6), and variations to LLSimpute, sequential LLSimpute (SLLSimpute) (7), iterative LLSimpute (ILLSImpute) (8), robust least square estimation with principal components (RLSP) (9), Bayesian gene selection BGSregress (10), collateral missing value imputation (CMVE) (11), and auto-regressive least square imputation (ARLS) (12) are all examples of local approaches. On the other hand, SVDimpute (Singular Value Decomposition) (2), Bayesian Principal Component Analysis (BPCA) (13) are the examples of global approach and utilize the global correlation present in the entire gene expression matrix. Hybrid approaches include methods like LinCmb (14), HPM-MI (15), and tri-imputation (16). GOimpute (17), HAimpute (1), and (iMISS) (19) are knowledge-assisted methods that combine the already existing domain knowledge to imputation techniques for imputing missing values in gene expression data, thereby, increasing their imputation accuracy. Gene Ontology based similarity measure has been recently used for missing value imputation in miRNA microarray data (20). A brief review of all the existing methods is shown in Table 1.

TABLE 1

Table 1. Review of existing methods for missing value imputation in gene expression data.

Most of the methods perform missing value imputation in gene expression data at comparatively higher observability, say, when 70% or more data is available (that is equivalent to 30% or less data is missing). Recent developments have made it possible to predict expression data values when the observed data is as low as 10%. Gene expression data is a highly correlated data because of the high level of interdependence between the genes. This interdependence is due to functional relationship between the genes as the group of genes interact together in any biological process. Therefore, it is evident that gene expression matrix is very similar to a low rank matrix that can be embedded into a lower dimensional subspace. Hence, imputation of missing values in data matrix has been projected as the matrix completion problem.

Matrix completion is a popular and challenging area of research in various domains. Many matrix completion methods exist in the literature and out of these methods, LMaFit (Low Rank Matrix Fitting) (21), LogDet (Logarithm determinant) (22), and Robust PCA (RPCA) (23) are three different state-of-the-art matrix completion methods. LMaFit is based on matrix factorization, while LogDet implements nuclear norm minimization. RPCA performs feature reduction and is quite robust to outliers. However, these methods have some limitations. LogDet becomes computationally expensive as the size of the matrix increases. LMaFit and RPCA-GD provide good performance, but their parameters need to be tuned properly for better recovery of missing values. Recently Kapur et al. (24) has used low rank constrained matrix completion method for imputing missing values in genomics.

In this paper, a novel 2-stage method, DSNN (Doubly Sparse DCT domain with Nuclear Norm minimization), has been proposed for predicting missing values in gene expression data using Compressive Sensing (CS) based formulation. In the first stage, missing values were recovered in gene expression data by formulating it as the CS-based reconstruction with double sparsity in the Discrete Cosine Transform (DCT). It has been shown in Gupta et al. (25) that DCT acts as approximate Karhunen-Loève type transform for a large class of signals, particularly, for slowly varying signals. Inspired by this, researchers have used DCT-based CS recovery to impute gene expression data in recent TV-DCT method (26) and CT-NNBI method (27), although only column sparsity was utilized in both these methods, while DSNN utilizes double sparsity.

Matrix obtained in first stage is considered a noisy version of the original matrix. Therefore, in Stage-2, denoising of the matrix recovered from Stage-1 is done by utilizing nuclear norm minimization. It exploits the low rank property of the data matrix. Missing value imputation was performed on four blood cancer dataset at different observability of data (10–90%) using NMSE as evaluation metric. Significance of imputation was validated by two experiments. In the first experiment, classification of normal vs. cancer subjects was carried out. In the second experiment, biological significance of imputation was ascertained by first identifying top 500 genes using SPARROW algorithm (28), followed by KEGG and GO analysis on these top 500 genes. SPARROW (SPARse selected expRessiOn regulators identified With penalized regression) algorithm finds candidate tumor drivers from the “selected expression regulators” (SERs). It defines SERs as the genes that drive dysregulated transcription leading to carcinogenesis. This algorithm regresses the gene expression values on the candidate SERs and provides a rank to each SERs based on the genes expression values of the samples. The method has been described briefly in section 3. Once the ranking was done by SPARROW, top 500 ranked genes from the list were further studied by KEGG (29–31) and GO pathway (32, 33) analysis using a web based application, Enrichr, developed and maintained by Chen et al. (34) and Kuleshov et al. (35).

2. Materials and Methods

2.1. Dataset Description

Four publicly available microarray gene expression dataset of different cancer types and different population have been used. Dataset-1 is Chronic lymphocytic leukemia (CLL) cancer dataset (GSE50006) submitted by Dana-Farber Cancer Institute, USA. CLL dataset contains expression values of 220 subjects across 54675 probe-ids and consists of two classes depending on whether the subject has CLL or not. There are 188 tumor samples and rest 32 are normal samples. Dataset-2 is Acute myeloid leukemia (AML) cancer dataset (GSE9476) (36) from Fred Hutchinson Cancer Research Centre, USA. It contains gene expression values of 64 subjects across 22283 probe-ids. Two classes are present in the data. Label “1” corresponds to person suffering from AML and label “2” corresponds to healthy subject. There are 26 tumor subjects and 38 healthy subjects. Dataset-3 is Multiple Myeloma (MM) cancer dataset (GSE47552) (37) from Centro de Investigacion del Cancer de Salamanca, Spain. It contains gene expression data of 99 subjects across 33297 probe-ids. It has data from 20 subjects with MGUS, 33 with high-risk SMM, 41 with MM, and rest 5 were healthy subjects. Dataset-4 is Multiple Myeloma (MM) cancer dataset (GSE125361) belonging to Indian population. It contains gene expression data of 48 MM subjects across 58341 probe-ids.

Data was pre-processed to convert probe-ids to gene symbols because gene vs. sample information is required for SPARROW analysis. It was observed that several probe-ids showed same gene names. To overcome this problem, gene expression levels of the probe-ids corresponding to the same gene name were averaged and gene vs. sample matrix was created. After pre-processing, CLL dataset had 220 samples with expression values of 23,348 genes. AML dataset had 64 samples with expression values of 13,650 genes. MM-Spanish dataset had 99 samples with expression values of 23,307 genes. MM-Indian dataset had 48 samples with gene expression values of 33,973 genes. Since the range of gene expression values was very high (of the order of 10⁶) for the CLL dataset, data was log transformed to reduce its dynamic range and to ensure that the smaller values were not shadowed by the higher values during the missing data recovery method.

$\begin{array}{ll}{\text{X}}_{\text{log-transformed}}(i,j)=lo{g}_{10}({\text{X}}_{\text{original}}(i,j)+1)& (1)\end{array}$

Matrix imputation was carried out on the sample vs. gene matrices. After matrix imputation, only tumor samples of both the dataset were used for SPARROW analysis.

2.2. Method

Workflow pipeline of the proposed analysis is shown in Figure 1. First of all, pre-processing of raw data was done as described in the previous section. Next, missing value imputation was carried out on four blood cancer dataset at different observability of data using Normalized Mean Square error (NMSE) as evaluation metric. Significance of imputation was validated by two experiments. In the first experiment, classification of normal vs. cancer subjects was carried out. In the second experiment, biological significance of imputation was ascertained using SPARROW algorithm (28) followed by KEGG and GO analysis on the top 500 genes identified by SPARROW.

FIGURE 1

Figure 1. Workflow of the proposed analysis.

2.2.1. Proposed DSNN Method of Matrix Imputation

The proposed “Doubly Sparse DCT domain with Nuclear Norm minimization” (DSNN) method consists of two stages. Stage-1 imputes missing values using a CS-based framework and DCT-based sparsity, while Stage-2 removes noise from the matrix obtained from Stage-1 by using a simple denoising framework.

2.2.1.1. Stage-1: Compressive sensing based matrix completion

In this stage, missing value problem was projected as compressive sensing based reconstruction problem. To understand it better, consider an incomplete matrix Y of size r × s, where r represents the number of subjects and s denotes the number of genes. Since the expression value of any gene will not vary much across subjects, data within a column would be sparse in some transform domain. Similarly, for a sample, gene expression levels of the gene will also be sparse in some transform domain. Columns and rows of the gene expression matrix were studied in the DCT domain and were observed to be highly sparse as shown in Figure 2. Based on this observation, Discrete Cosine Transform was chosen as the sparsifying transform in DSNN method because DCT acts as a KL-type basis for slow-varying signals (25) and data is sparse in the DCT domain.

FIGURE 2

Figure 2. Each curve represents DCT coefficients of a few randomly chosen columns and rows of gene expression matrices of CLL dataset.

Thus, the missing data recovery problem was formulated in a compressive sensing framework, where the sensing matrix Φ was of size r × s and had “0” entries for missing values in data matrix Y, while rest of the entries were “1.” Corresponding to each observed entry (that is not missing) of the ith column, there is a row in Φ_i with an entry “1” for the corresponding position and zeros in the rest of the positions. For example, assume x_missing=${\left[\begin{array}{cccccc}{x}_1& .& x_3& .& .& x_6\end{array}\right]}^T$ is the observed vector where only x₁, x₃, and x₆ are available and, x₂, x₄, and x₅ are missing (denoted as “.” in the vector). Then, the vector x_missing can be re-written as y:

$\begin{array}{ll}\text{y}=\Phi \text{x}& (2)\end{array}$

$\begin{array}{ll}\text{y}=\left[\begin{array}{c}{x}_1\\ x_3\\ x_6\end{array}\right]=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\end{array}\right],& (3)\end{array}$

where the sensing matrix is written as $\Phi =\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]$ and x is the desired vector to be recovered. This is the standard formulation in compressive sensing literature, where it is assumed that only few values of data are sensed. In the above example, these values are x₂, x₄, and x₅. Thus, we have recast the problem of missing values in vector x_missing as the compressively sensed vector y. Now, the task is to recover full data x from compressively sensed data y that will lead to missing value recovery.

Gene expression data was interpreted as a matrix with few observed samples, where the goal was to reconstruct the original matrix from the observed entries using DCT-based sparsity of gene expression data.

The following optimization problem was solved to recover the missing values in Y

$\begin{array}{ll}\underset{\tilde{\text{X}}}{min}({\Vert \text{Y}-\Phi \tilde{\text{X}}\Vert }_2^2+{\lambda }_1{\Vert {\text{D}}_c\tilde{\text{X}}{\text{D}}_r^T\Vert }_1),& (4)\end{array}$

where D_c is columnwise DCT matrix applied on columns of the $\widetilde{\mathbf{\text{X}}}$ and D_r is the rowwise DCT matrix applied on rows of the $\widetilde{\mathbf{\text{X}}}$. $\widetilde{\mathbf{\text{X}}}$ is the matrix to be recovered. The above formulation is also known as analysis-prior and presence of DCT matrices in the formulation makes it non-separable. Using the orthogonal property of DCT transform, analysis prior was transformed to synthesis-prior formulation as

$\begin{array}{ll}\underset{\text{Z}}{min}({\Vert \text{Y}-\Phi {\text{D}}_c^T\text{Z}{\text{D}}_r\Vert }_2^2+{\lambda }_1{\Vert \text{Z}\Vert }_1),& (5)\end{array}$

where ${\mathbf{\text{D}}}_c\widetilde{\mathbf{\text{X}}}{\mathbf{\text{D}}}_r^T=\mathbf{\text{Z}}$. The above optimization problem was solved using the function handle and “SPGL1” solver (38, 39), where the regularization parameter λ₁ was chosen automatically by the “SPGL1” solver.

2.2.1.2. Stage-2: denoising framework

It was assumed that the recovered $\widetilde{\mathbf{\text{X}}}$ from Stage-1 is the noisy version of the original matrix X and hence, the recovered matrix was denoised in Stage-2. Before denoising, $\widetilde{\mathbf{\text{X}}}$ was re-organized into ${\widetilde{\mathbf{\text{X}}}}_{\mathrm{\text{rec}}}$ as

$\begin{array}{ll}{\tilde{\text{X}}}_{\text{rec}}(j,i)=\{\begin{array}{ll}0,\hfill & \text{if}(|\tilde{\text{x}}(j,i)-\text{mean}({\text{y}}_{\text{i}})|\ge {\lambda }_{\text{2}}\text{std}({\text{y}}_{\text{i}})\hfill \\ \tilde{\text{x}}(j,i),\hfill & \text{otherwise}\hfill \end{array}& (6)\end{array}$

where j ranges from 1 to m (number of rows/ subjects), |.| denotes the absolute value and, mean(y_i) and std(y_i) denote the mean and the standard deviation of the ith column of the initial observed (but incomplete) matrix Y. Parameter λ₂ was determined empirically and was set to value 0.2 for experiments on CLL dataset, MM-Spanish dataset, and MM-Indian dataset. It was set to 0.1 for experiments on AML dataset. Denoising was formulated in the Split-Bregman type optimization as

$\begin{array}{ll}\underset{\text{W}}{min}({\Vert \text{W}\Vert }_{*}+{\lambda }_3{\Vert \text{W}-\widehat{\text{X}}-\text{B}\Vert }_F^2)s.t.\widehat{\text{X}}=\text{W},& (7)\end{array}$

where $\widehat{\mathbf{\text{X}}}$ was initialized as:

$\begin{array}{ll}\widehat{\text{X}}={\tilde{\text{X}}}_{\text{rec}}+{\tilde{\text{X}}}_{\text{inv-rec}}◦rand(m,n),& (8)\end{array}$

where “°” represents the Hadamard product of two matrices with the elements of ${\widetilde{\mathbf{\text{X}}}}_{\mathrm{\text{inv-rec}}}$ defined as

$\begin{array}{ll}{\tilde{\text{X}}}_{\text{inv-rec}}(j,i)=\{\begin{array}{ll}1,\hfill & \text{if }{\tilde{\text{X}}}_{\text{rec}}(j,i)=0,\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}& (9)\end{array}$

Equation (5) was solved in Split Bregman type iterations as

$\begin{array}{ll}{\text{W}}^{k+1}=SV{T}_{{\lambda }_3}({\widehat{\text{X}}}^k+{\text{B}}^k),& (10)\end{array}$

$\begin{array}{ll}{\text{B}}^{k+1}={\widehat{\text{X}}}^k+{\text{B}}^k-{\text{W}}^k,& (11)\end{array}$

$\begin{array}{ll}{\widehat{\text{X}}}^{k+1}={\tilde{\text{X}}}_{\text{rec}}+{\tilde{\text{X}}}_{\text{obs}}◦{\text{W}}^{k+1},& (12)\end{array}$

where “SVT” denotes the soft singular value thresholding method (40) and ${\widetilde{\mathbf{\text{X}}}}_{obs}$ is the observed incomplete matrix. Optimal value of parameter λ₃ was determined using grid search and was set to 100 in all experiments. The complete algorithm for the proposed DSNN method is presented below.

ALGORITHM 1

Algorithm 1. Proposed DSNN Method

3. Results

3.1. Evaluation

For assessing the performance of the proposed DSNN method, some data were dropped randomly to create incomplete matrices with available data ranging from 10 to 90%. Next, incomplete matrices were imputed using the DSNN method. Results were simultaneously generated using three state-of-the-art matrix completion methods for comparative analysis. Normalized mean squared error (NMSE) was used as the evaluation metric and was calculated between the original and the recovered matrix. NMSE is defined as:

$\begin{array}{ll}\text{NMSE}=\frac{{\Vert \text{X}(original)-\widehat{\text{X}}(recovered)\Vert }_F^2}{{\Vert \text{X}(original)\Vert }_F^2}.& (13)\end{array}$

Semi-log plots of NMSE at different stages are shown in Figure 3. Stage-1 results were obtained when missing values in data matrix were imputed using compressive sensing based matrix completion, where double sparsity in DCT domain was exploited. Stage-2 results were obtained when only nuclear norm minimization was used for matrix imputation. DSNN method combined both these stages. Results clearly indicated that the performance of imputation has improved with the two successive stages of DSNN. DSNN method also worked better than the existing methods even at high missing rates of 10% as shown in Figure 4. NMSE reported in the figures is averaged over 30 iterations. For CLL dataset, highest NMSE reported was 0.09 at 10% observed data and lowest NMSE was 0.004 at 90% observed data. For AML dataset, highest NMSE was 0.013 at 10% observed data and lowest NMSE was 0.00056 at 90% observed data. For MM-Spanish dataset, highest NMSE reported was 0.005 at 10% observed data and lowest NMSE was 0.00039 at 90% observed data. For MM-Indian dataset, highest NMSE was 0.0122 at 10% observed data and lowest was 6.25E-04 at 90% observed data.

FIGURE 3

Figure 3. Semi-log plots show NMSE after imputation on (A) CLL, (B) AML, (C) MM-Spanish, and (D) MM-Indian dataset using Stage-1 only, Stage-2 only, and Proposed DSNN method (Stage-1 + Stage-2).

FIGURE 4

Figure 4. Small Semi-log plots showing comparison of the proposed DSNN method with the three state-of-the-art methods in terms of NMSE for (A) CLL, (B) AML, (C) MM-Spanish, and (D) MM-Indian dataset.

3.2. Validation

In order to determine the significance of the DSNN method, two separate experiments were carried out on the original data, incomplete data, and imputed data matrices. In experiment-1, classification of normal vs. cancer subjects was carried out. In experiment-2, biological significance of imputation was ascertained by first identifying top candidate tumor drivers from SPARROW algorithm, followed by gene enrichment analysis on the top-ranked genes using web based application Enrichr.

3.2.1. Experiment 1: Classification

Simulation results on missing value recovery were validated by performing classification on original matrices, matrices with random missing values, and imputed matrices of the CLL and AML dataset. Classification can be either supervised or unsupervised depending on the availability of ground truth labels. In these dataset, ground truth labels were available. Hence, supervised classification was performed to distinguish between two classes, normal and cancer using two different classifiers: linear Support Vector Machine (SVM) and k nearest neighbor (KNN) method with k = 3. Both the dataset had large number of features, therefore, feature reduction was performed to extract important features from the data. Three different methods of feature reduction were used, Mutual Information criterion, Principal Component Analysis (PCA) and Chi-square method. Optimal number of features in each method were estimated by grid search. Further, five-fold cross validation was performed and average accuracy over 20 iterations was reported. Experiments were performed in Python 3 environment with Sklearn 0.20 library. Classification code was written in Python programming language. Scikit-learn is a Python module for machine learning and contains various algorithms related to regression, classification, and clustering. Examples of these algorithms are support vector machines (SVM), random forest (RF), k-means. Classification accuracy and F₁ score were calculated at different sampling ratios from 10 to 90%. The accuracy and F₁ score are defined as

$\begin{array}{ll}\text{Accuracy}=\frac{1}{M}{\displaystyle \sum _{i=1}^{M}1(x_i={\tilde{x}}_i)},& (14)\end{array}$

$\begin{array}{ll}\text{and }{F}_1=\frac{2\times \text{precision }\times \text{ recall}}{\text{precision }+\text{ recall}},& (15)\end{array}$

where M is the total number of samples in the dataset, x_i is the ground truth class label of the ith sample, and ${\tilde{x}}_i$ is the predicted class label. Weighted F₁ score was used in order to account for label imbalance arising out due to unequal number of tumor and normal samples. CLL dataset had 188 tumor and 32 normal samples and AML dataset had 26 tumor samples and 38 normal. Tables 2, 3 clearly indicate that values of classification accuracy and F₁ scores for incomplete matrices are low as compared to the values obtained on imputed matrices. Classification accuracy and F₁ scores were also computed on imputed matrices obtained from the three existing methods on both the dataset and compared with the results of DSNN method as shown in Figures 5, 6. Classification was also performed on MM-Spanish dataset (Results are shown in Table T1). Classification could not be performed on MM-Indian data because it was a single class data, i.e., of tumor samples only.

TABLE 2

Table 2. Classification accuracy and F₁ score for CLL dataset at varying sampling ratios (FR, feature reduction; SR, sampling ratio; Obs., observed; Rec., recovered using DSNN method.

TABLE 3

Table 3. Classification Accuracy and F₁ score for AML dataset at varying sampling ratios (FR, feature reduction; SR, sampling ratio; Obs., observed; Rec., recovered using DSNN method).

FIGURE 5

Figure 5. Comparison of different methods in terms of classification accuracy and F₁ score at varying sampling ratios on CLL dataset.

FIGURE 6

Figure 6. Comparison of different methods in terms of classification accuracy and F₁ score at varying sampling ratios on AML dataset.

3.2.2. Experiment 2: Biological Validation

For biological validation of the results, SPARROW was applied on the original matrix, incomplete matrices, and imputed matrices to identify top candidate tumor driver genes. SPARROW (SPARse selected expRessiOn regulators identified With penalized regression) was proposed by Logsdon et al. (28) and aims to find out candidate tumor drivers from the “selected expression regulators” (SERs). It defines SERs as the genes that drive dysregulated transcription leading to carcinogenesis. In this method, variational Bayesian spike regression model has been used to fit the following model,

$\begin{array}{ll}{y}_{m,n}={\displaystyle \sum x_{m,k}}{\beta }_{k,n}+e_{m,n},& (16)\end{array}$

where y_m,k is the value of expression of the nth gene for the mth subject, e_m,n is a normally distributed error, x_m,k is the value of expression of the kth SER for the mth subject and β_k,n is the additive effect of the expression of the kth SER on the expression of the nth gene. m ranges from 1….M, where M is the total number of subjects and n ranges from 1…….N, where N is the total number of genes. Total SERs used in the analysis were around 3,400 and they were downloaded from the link provided in the original paper. This algorithm provides a rank to each SER based on the gene expression values of the samples. The top-ranked genes from the list can be further studied by gene enrichment analysis.

For finding top 500 candidate driver genes, only the tumor samples from the data matrices were considered for SPARROW analysis. Algorithm was applied on original complete data matrices of all the dataset to identify the top-ranked candidate tumor drivers. This served as the ground truth for our analysis. Further, SPARROW was applied on incomplete and imputed data matrices of both the dataset at sampling ratios of 50 and 70%. Top-ranked candidate drivers from the incomplete and imputed data matrices were obtained. Gene enrichment analysis was performed on top 500 genes. KEGG and GO pathways were studied using web based application, Enrichr, developed, and maintained by Chen et al. (34) and Kuleshov et al. (35). KEGG pathways obtained from gene lists of original dataset were the ground truth. It was observed that when KEGG pathway analysis was done for incomplete matrices, these were not able to predict cancer pathways with a higher significance (low p-value) whereas for imputed matrices, cancer pathways were predicted with a higher significance due to decrease in p-value. Results from KEGG analysis on all dataset are presented in tabular form showing the p-values, combined score for original data and the incomplete and complete matrices in tabular form in Tables T2–T5. p-value was computed from the Fisher exact test. Fisher test was run on random gene sets and ranks were derived at each run. Mean rank was calculated from the different runs and standard deviation of the rank obtained from the expected rank was also calculated for each term in the gene-set library. Finally, a z-score was calculated to estimate the deviation from the expected rank. z-score and p-value were used to compute combined score which is obtained by multiplying z-score with the logarithm of p-value. A detailed analysis for CLL dataset consisting of z-score and combined score has also been shown in Tables T6, T7. Similarly, GO pathways were obtained from the Enrichr and bubble plots were constructed from GO tables as shown in Figures S5–S12.

4. Discussion

4.1. Importance of the Proposed DSNN Method

DSNN, a two stage method proposed for matrix recovery in the paper, was based on Compressive Sensing Framework. In Stage-1, it utilized column and row sparsity of the gene expression matrix in DCT domain for missing value imputation, while in Stage-2, it exploited low rank nature of the matrix for denoising. Expression values of any particular gene would vary slowly across subjects, thereby, exhibiting sparsity in columns in some transformed domain. Similarly, expression values of a subject for most of the genes will also be slowly varying, thereby, exhibiting sparsity in the rows. Since there is a high inter-dependence between the expression levels of the genes, one may consider gene expression matrix as a low rank matrix. Thus, as discussed earlier, both the assumptions used in Stage-1 (of sparsity in DCT domain) and Stage-2 (low rank of matrix) hold true for the given gene expression data. The concept of DCT-based sparsity was recently applied on biological data in methods, TV-DCT (26) and CT-NNBI (27), although only column sparsity in the DCT domain was used. On the other hand, this work utilizes double sparsity, i.e., sparsity on both the columns and the rows. Most of the imputation algorithms developed for missing value imputation such as KNN, LSimpute, LLSimpute, BPCA etc. work at high observability of data, while the proposed DSNN method worked well even when data had very high missing rates of 10–40%. The proposed DSNN method performed better than the other matrix completion methods at all sampling ratios. The state-of-the-art matrix imputation methods that have been used for performance comparison in this work required a lot of parameter tuning for optimal performance, while DSNN method did not require parameter tuning to such a great extent.

4.2. Improvement in Classification Accuracy

It was evident from the results shown in Tables 2, 3 that the classification accuracy and F₁ scores reduced as the number of missing values increased. There were 220 samples in CLL dataset and 64 samples in AML dataset. For smaller dataset like AML, missing values affected the classification accuracy and F₁ scores greatly. Thus, it is necessary to impute missing values in gene expression data to prevent incorrect downstream analysis of the data. When the classification was performed on the imputed data, there was considerable improvement in the classification accuracy, thereby, validating our hypothesis. Classification accuracy and F₁ scores calculated on original complete data matrices (100% sampling ratio) were considered as ground truth values. For CLL dataset, ground truth values of classification accuracy and weighted F₁ score were 0.99 and 0.99, respectively, as shown in Figure 5. For KNN classifier and Chi-square feature selection approach, classification accuracy and F₁ score obtained for 50% observed data was 0.86 and 0.80, respectively as shown in Table 2. After imputation, values improved significantly to 0.99 and 0.99. For AML dataset, ground truth values of classification accuracy and F₁ score were 1.0 and 1.0, respectively as shown in Figure 6. Similarly for Linear SVM classifier and PCA feature selection approach, classification accuracy and F₁ score for 50% observed data was 0.56 and 0.57, respectively, as shown in Table 3. After matrix imputation, classification accuracy and F₁ score improved considerably to 0.99 and 0.99, respectively. For every sampling ratio, consistent results were obtained that validates our method.

4.3. Improvement in Functional Enrichment Analysis for KEGG Pathways

KEGG and GO enrichment analysis was performed on the top 500 ranked genes obtained from SPARROW algorithm to biologically validate our results. As mentioned earlier, KEGG pathways obtained by the top-ranked genes of original matrices were considered the ground truth values. Pathways with p-value < 0.05 were only considered. When KEGG analysis was done on top-ranked genes from incomplete matrices, there was significant decrease in the p-value of the most significant pathways. “Wnt signaling pathway” (41, 42) and “Notch signaling pathway” (43, 44) are important pathways in CLL cancer. An important observation was that p-value for “Notch signaling pathway” was 2.00E-01 at ground truth and it was 5.76E-02 at 70% observed data for CLL dataset. Values were insignificant in both the cases. However, after imputation, p-value became significant with value 1.56E-02 which was <0.05 as shown in Figure 7.

FIGURE 7

Figure 7. Few important KEGG pathways at 70% observed and imputed data for CLL data. Adjusted p-values are shown in brackets.

Similarly, p-value for “Wnt signaling pathway” was 8.33E-05 on original dataset, as shown in Figure S1. At 50% observed data p-value for “Wnt signaling pathway” was 3.10E-02 which was less significant than the ground truth value at 50% observed data. After matrix imputation, p-value became significant with value 2.13E-03. Similarly, p-value became 1.90E-05 after matrix imputation on 70% observed data which was more significant than the p-value 6.66E-5, observed at 70% data. “Fc epsilon RI signaling pathway” is an important pathway in AML cancer (45). This pathway was insignificant for original data with p-value-2.12E-01. At 70% observed data, p-value was 9.40E-02 which was again greater than 0.05. After matrix imputation, the value became significant at 2.75E-02, which was less than 0.05 as shown in Figure 8. Similarly, ‘Ras signaling pathway” is activated in Multiple Myeloma cancer (46). For MM-Spanish data, “Ras signaling pathway” was significant with p-value-0.0052 for original data but became insignificant with p-value-0.23 when 70% data was observed as shown in Figure 9. After matrix imputation, significance of the pathway was restored with p-value 0.04. For MM-Indian dataset,“Transcriptional misregulation in cancer” was found to be insignificant with p-value 0.47 as shown in Figure 10. After imputation, p-value decreased to 1.37E-03 and became more significant than ground truth p-value, 7.8E-03. Additional KEGG analysis results on the dataset CLL, AML and MM Spanish data are provided in the Figures S2–S4.

FIGURE 8

Figure 8. Few important KEGG pathways at 70% observed and imputed data for AML data. Adjusted p-values are shown in brackets.

FIGURE 9

Figure 9. Few important KEGG pathways at 70% observed and imputed data for MM-Spanish data. Adjusted p-values are shown in brackets.

FIGURE 10

Figure 10. Few important KEGG pathways at 70% observed and imputed data for MM-Indian data. Adjusted p-values are shown in brackets.

Thus, DSNN method not only imputed missing entries but also performed some denoising to improve the results. It is quite evident from the analysis that gene enrichment analysis results were partially inaccurate due to incomplete matrices. This was because the genes identified as top-ranked genes by performing SPARROW analysis on complete data matrix were not identified in the top-ranked list obtained from incomplete data matrix. However, when the incomplete matrix was imputed using the proposed DSNN method, top-ranked list of genes obtained from SPARROW analysis was quite similar to the ground truth. Our observations demonstrate the importance of imputing missing values in gene expression data.

5. Conclusion

Microarray data generally has a lot of missing values that can adversely influence the downstream analysis. In this paper, a new method, namely DSNN, is proposed that imputes missing values in the gene expression data using discrete cosine transform based double sparsity and nuclear norm minimization. Method was also validated quantitatively based on the application of classification approach as well as biologically by performing pathway enrichment analysis and showed consistent findings.

Data Availability Statement

Data was downloaded from GEO database with accession numbers, GSE9476, GSE50006, and GSE47552 and GSE125361. Code of the DSNN method proposed in this manuscript is available publicly at github:https://github.com/AkankshaFarswan/DSNN_MatrixImputation.

Author Contributions

AF and AG conceptualized the method of data imputation. RG and GK contributed to the ideas on method application and validation in bioinformatics. AF prepared all codes, generated results, and prepared the first draft of the paper. All authors edited the manuscript for final submission.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

AF would like to thank University Grants Commission, Govt. of India for UGC-Senior Research Fellowship. AG and RG would like to thank the Department of Science and Technology, Govt. of India for providing the research grant support [Grant: DST/ICPS/CPS-Individual/2018/279(G)]. Authors also acknowledge the Infosys Center for Artificial Intelligence, IIIT-Delhi for their research work.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2019.01442/full#supplementary-material

References

1. Xiang Q, Dai X, Deng Y, He C, Wang J, Feng J, et al. Missing value imputation for microarray gene expression data using histone acetylation information. BMC Bioinformatics. (2008) 9:252. doi: 10.1186/1471-2105-9-252

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Troyanskaya O, Cantor M, Sherlock G, Brown P, Hastie T, Tibshirani R, et al. Missing value estimation methods for DNA microarrays. Bioinformatics. (2001) 17:520–5. doi: 10.1093/bioinformatics/17.6.520

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Kim KY, Kim BJ, Yi GS. Reuse of imputed data in microarray analysis increases imputation efficiency. BMC Bioinformatics. (2004) 5:160. doi: 10.1186/1471-2105-5-160

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Brás LP, Menezes JC. Improving cluster-based missing value estimation of DNA microarray data. Biomol Eng. (2007) 24:273–82. doi: 10.1016/j.bioeng.2007.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Ouyang M, Welsh WJ, Georgopoulos P. Gaussian mixture clustering and imputation of microarray data. Bioinformatics. (2004) 20:917–23. doi: 10.1093/bioinformatics/bth007

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Bø TH, Dysvik B, Jonassen I. LSimpute: accurate estimation of missing values in microarray data with least squares methods. Nucleic Acids Res. (2004) 32:e34. doi: 10.1093/nar/gnh026

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Zhang X, Song X, Wang H, Zhang H. Sequential local least squares imputation estimating missing value of microarray data. Comput Biol Med. (2008) 38:1112–20. doi: 10.1016/j.compbiomed.2008.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Cai Z, Heydari M, Lin G. Iterated local least squares microarray missing value imputation. J Bioinformatics Comput Biol. (2006) 4:935–57. doi: 10.1142/s0219720006002302

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Yoon D, Lee EK, Park T. Robust imputation method for missing values in microarray data. BMC Bioinformatics. (2007) 8:S6. doi: 10.1186/1471-2105-8-S2-S6

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Zhou X, Wang X, Dougherty ER. Missing-value estimation using linear and non-linear regression with Bayesian gene selection. Bioinformatics. (2003) 19:2302–7. doi: 10.1093/bioinformatics/btg323

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Sehgal MSB, Gondal I, Dooley LS. Collateral missing value imputation: a new robust missing value estimation algorithm for microarray data. Bioinformatics. (2005) 21:2417–23. doi: 10.1093/bioinformatics/bti345

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Choong MK, Charbit M, Yan H. Autoregressive-model-based missing value estimation for DNA microarray time series data. IEEE Trans Inform Technol Biomed. (2009) 13:131–7. doi: 10.1109/TITB.2008.2007421

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Oba S, Sato Ma, Takemasa I, Monden M, Matsubara Ki, Ishii S. A Bayesian missing value estimation method for gene expression profile data. Bioinformatics. (2003) 19:2088–96. doi: 10.1093/bioinformatics/btg287

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Jörnsten R, Wang HY, Welsh WJ, Ouyang M. DNA microarray data imputation and significance analysis of differential expression. Bioinformatics. (2005) 21:4155–61. doi: 10.1093/bioinformatics/bti638

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Purwar A, Singh SK. Hybrid prediction model with missing value imputation for medical data. Expert Syst Appl. (2015) 42:5621–31. doi: 10.1016/j.eswa.2015.02.050

CrossRef Full Text | Google Scholar

16. He C, Zhao C, Li GZ, Zhu W, Yang W, Yang MQ. A hybrid iterative approach for microarray missing value estimation. In: 2016 IEEE International Conference on Bioinformatics and Biomedicine (BIBM). Shenzhen: IEEE (2016). p. 1350–2.

Google Scholar

17. Tuikkala J, Elo L, Nevalainen OS, Aittokallio T. Improving missing value estimation in microarray data with Gene Ontology. Bioinformatics. (2005) 22:566–72. doi: 10.1093/bioinformatics/btk019

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Alizadeh AA, Eisen MB, Davis RE, Ma C, Lossos IS, Rosenwald A, et al. Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature. (2000) 403:503–11. doi: 10.1038/35000501

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Hu J, Li H, Waterman MS, Zhou XJ. Integrative missing value estimation for microarray data. BMC Bioinformatics. (2006) 7:449. doi: 10.1186/1471-2105-7-449

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Yang Y, Xu Z, Song D. Missing value imputation for microRNA expression data by using a GO-based similarity measure. BMC Bioinformatics. (2016) 17(Suppl. 1):S10. doi: 10.1186/s12859-015-0853-0

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wen Z, Yin W, Zhang Y. Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm. Math Programm Comput. (2012) 4:333–61. doi: 10.1007/s12532-012-0044-1

CrossRef Full Text | Google Scholar

22. Kang Z, Peng C, Cheng Q. Top-N recommender system via matrix completion. In: AAAI. Phoenix, AZ (2016). p. 179–185.

Google Scholar

23. Yi X, Park D, Chen Y, Caramanis C. Fast algorithms for robust PCA via gradient descent. In: Advances in Neural Information Processing Systems. Barcelona (2016). p. 4152–60.

Google Scholar

24. Kapur A, Marwah K, Alterovitz G. Gene expression prediction using low-rank matrix completion. BMC Bioinformatics. (2016) 17:243. doi: 10.1186/s12859-016-1106-6

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Gupta A, Joshi SD, Singh P. On the approximate discrete KLT of fractional Brownian motion and applications. J Franklin Inst. (2018) 355:8989–9016. doi: 10.1016/j.jfranklin.2018.09.023

CrossRef Full Text | Google Scholar

26. Farswan A, Gupta A. TV-DCT: Method to impute gene expression data using DCT based sparsity and total variation denoising. In: ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). Brighton (2019). p. 1244–8.

Google Scholar

27. Gehlot S, Farswan A, Gupta A, Gupta R. CT-NNBI: Method to impute gene expression data using DCT based sparsity and nuclear norm constraint with split bregman iteration. In: 2019 IEEE 16th International Symposium on Biomedical Imaging (ISBI 2019). Venice: IEEE (2019). p. 1315–8.

Google Scholar

28. Logsdon BA, Gentles AJ, Miller CP, Blau CA, Becker PS, Lee SI. Sparse expression bases in cancer reveal tumor drivers. Nucleic Acids Res. (2015) 43:1332–44. doi: 10.1093/nar/gku1290

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. (2000) 28:27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Kanehisa M, Sato Y, Furumichi M, Morishima K, Tanabe M. New approach for understanding genome variations in KEGG. Nucleic Acids Res. (2018) 47:D590–5. doi: 10.1093/nar/gky962

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. (2019) 28:1947–51. doi: 10.1002/pro.3715

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. (2000) 25:25–9. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

33. The Gene Ontology Consortium. The Gene Ontology resource: 20 years and still GOing strong. Nucleic Acids Res. (2018) 47:D330–8. doi: 10.1093/nar/gky1055

CrossRef Full Text | Google Scholar

34. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. (2013) 14:128. doi: 10.1186/1471-2105-14-128

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. (2016) 44:W90–7. doi: 10.1093/nar/gkw377

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Stirewalt DL, Meshinchi S, Kopecky KJ, Fan W, Pogosova-Agadjanyan EL, Engel JH, et al. Identification of genes with abnormal expression changes in acute myeloid leukemia. Genes Chromos Cancer. (2008) 47:8–20. doi: 10.1002/gcc.20500

PubMed Abstract | CrossRef Full Text | Google Scholar

37. López-Corral L, Corchete LA, Sarasquete ME, Mateos MV, García-Sanz R, Fermiñán E, et al. Transcriptome analysis reveals molecular profiles associated with evolving steps of monoclonal gammopathies. Haematologica. (2014) 99:1365–72. doi: 10.3324/haematol.2013.087809

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Van den Berg E, Friedlander MP. Probing the Pareto frontier for basis pursuit solutions. SIAM J Sci Comput. (2008) 31:890–912. doi: 10.1137/080714488

CrossRef Full Text | Google Scholar

39. Van den Berg E, Friedlander MP. SPGL1: A Solver for Large-Scale Sparse Reconstruction (2007). Available online at: https://www.cs.ubc.ca/~mpf/spgl1/download.html

Google Scholar

40. Cai J, Candès EJ, Shen Z. A singular value thresholding algorithm for matrix completion. SIAM J Optim. (2010) 20:1956–82. doi: 10.1137/080738970

CrossRef Full Text | Google Scholar

41. Wang L, Shalek AK, Lawrence M, Ding R, Gaublomme JT, Pochet N, et al. Somatic mutation as a mechanism of Wnt/β-catenin pathway activation in CLL. Blood. (2014) 124:1089–98. doi: 10.1182/blood-2014-01-552067

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Zhan T, Rindtorff N, Boutros M. Wnt signaling in cancer. Oncogene. (2017) 36:1461–73. doi: 10.1038/onc.2016.304

PubMed Abstract | CrossRef Full Text | Google Scholar

43. De Falco F, Sabatini R, Del Papa B, Falzetti F, Di Ianni M, Sportoletti P, et al. Notch signaling sustains the expression of Mcl-1 and the activity of eIF4E to promote cell survival in CLL. Oncotarget. (2015) 6:16559–72. doi: 10.18632/oncotarget.4116

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Wickremasinghe R, Prentice A, Steele A. p53 and Notch signaling in chronic lymphocytic leukemia: clues to identifying novel therapeutic strategies. Leukemia. (2011) 25:1400–7. doi: 10.1038/leu.2011.103

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Gu Tl, Nardone J, Wang Y, Loriaux M, Villén J, Beausoleil S, et al. Survey of activated FLT3 signaling in leukemia. PLoS ONE. (2011) 6:e19169. doi: 10.1371/journal.pone.0019169

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Xu J, Pfarr N, Endris V, Mai EK, Hanafiah NM, Lehners N, et al. Molecular signaling in multiple myeloma: association of RAS/RAF mutations and MEK/ERK pathway activation. Oncogenesis. (2017) 6:e337. doi: 10.1038/oncsis.2017.36

PubMed Abstract | CrossRef Full Text | Google Scholar