Deep Cell-Type Deconvolution from Bulk Gene Expression Data Using DECODE

2.1.

Algorithmic background

We assume a bulk expression matrix V of n genes by m samples, where n rows represent the average of cell-type specific gene expression profiles in a sample, weighted by their abundances in that sample. Suppose there are k cell types and S is a signature matrix, where k columns are the gene expression profiles of those cell types. Our goal is to infer a matrix F, where m columns are probability vectors denoting the fraction of each cell type in the corresponding sample such that V ∼ S ⋅ F.

Our algorithm, DECODE is based on deep unfolding approach for NMF, DNMF [14]. DNMF is a deep learning algorithm for the decomposition of a non-negative matrix V_n×m into a product of two non-negative matrices S_n×k and F_k×m such that ∥V − SF∥₂ is minimal. In more detail, DNMF contains several (l) layers that are updated based on the multiplicative update rule of Lee and Seung [15]. This rule iteratively updates S and F as follows:

where ⊙, represent entry wise multiplication and division. However, in DNMF this rule is relaxed to allow learning better solutions. Specifically, the algorithm contains a layer for each iteration and aims to optimize F without explicitly representing S. Focusing on some column f of F and the corresponding column v of V , the output f_i of layer i is (up to regularization): where A_i+1 and B_i+1 are learnable matrices that correspond, in the original update rule, to the matrices S^T and S^T S, respectively (ignoring the dependency between the two latter matrices). Starting from an initial matrix F₀, the algorithm rolls it through t layers of the network, imitating the Lee and Seung iterative algorithm, until an output F_t matrix is produced.

DNMF has two model variants: supervised and unsupervised. The supervised variant assumes a known fraction matrix F_real which is either experimentally measured or generated in simulations. It is trained with an L₂ loss w.r.t. this matrix: ∥F_l − F_real∥₂. When no real matrix is available, DNMF uses an unsupervised variant where the loss is the reconstruction error ∥V − SF_l∥₂. Here we harness these two variants and develop a combined supervised-unsupervised DNMF model that maps from bulk expression (V ) to cell fractions (F).

2.2.

Model details and training process

We introduce several novelties into the DNMF architecture and training process. First, to account for the fact that each column f of F should represent a probability vector, we normalize f_i after each iteration to sum to one. Second, in order to make use of the given signature matrix S, we initialize a DNMF model M (S) where we set the weights of the first layer by A₁ = S^T and B₁ = S^T S to reflect one iteration of Lee and Seung’s update rule. In addition, we initialize F₀ with the result of applying NNLS to V and S. All other parameters are initialized to one as suggested in [14].

Last, we use a combined supervised and unsupervised training process, where the former is based on synthetic data while the latter is based on both synthetic and real data. We first train the model using synthetic data in a supervised fashion for one epoch with I = 60,000 iterations (or batches). Let M_k be the resulting model after the k-th iteration (we omit the signature matrix S from the notation for clarity). Next, we take M_I and train it for another I iterations in an unsupervised fashion on synthetic data to yield model M^U. Next, we train M^U for additional I_r = 100 epochs on real data and report the final model. This process is executed for each of a given list of signature matrices and the model with lowest unsupervised error (on the real data) overall is output.

2.3.

Supervised training on synthetic data

To deal with data scarcity, we first trained the algorithm on synthetic data. Specifically, we collected known ranges of cell frequencies from [16] and used values within these ranges as parameters for a Dirichlet distribution. We then drew random fraction vectors from this distribution. The resulting fraction matrix was multiplied by the given signature matrix to produce a synthetic bulk expression matrix. We drew multiple matrices in this fashion and fed them sequentially to the training process, viewing each such matrix as representing a batch. We further added a small normally-distributed noise to each bulk expression matrix with zero mean and small standard deviation: for the i-th batch the standard deviation is i∕10,000.

2.4.

Unsupervised training on synthetic data

In the unsupervised training, we trained the DECODE model (M_I) to minimize the reconstruction error ∥V − SF_model∥₂ for I iterations. Before we trained the model, we computed the NNLS supervised loss on the synthetic data and denote its loss by e_nnls. For each unsupervised iteration of DECODE, if the supervised loss exceeds e_nnls we stopped at this iteration (and not after I).

2.5.

Hyperparameter tuning

DECODE has several hyperparameters that need tuning. Due to the small size and number of data sets we opted for using an independent data for tuning the hyperparameters. For number of layers we used 4 for speed considerations, as the model’s accuracy is robust to the specific number used, and we used learning rate of 0.001 [14]. Since the problem we are trying to tackle is unsupervised in nature, we followed [14] and did not use regularization. In order to set the number of training iterations (I and I_r) we performed a grid search using an independent stromal dataset from [2] (see next section for detailed description). Specifically, we tried values from 40,000 to 80,000 for I and 50 to 200 for I_r, arriving at the best combination of I = 60,000 and I_r = 100. When applied to the stromal dataset, it can be seen that each of the algorithmic steps aids in improving DECODE’s performance (Figure 2).

2.6.

Data and preprocessing

Good training data for the deconvolution problem is scarce. Our main data source is a recent benchmark paper [2] which has three available datasets, all are results of in-silico simulations. Two of the datasets, PBMC1 and PBMC2, represent 200 simulated mixtures of single-cell peripheral blood mononuclear cell (PBMC) expression profiles (100 mixtures in each dataset). For each mixture, individual cells were randomly chosen and their expression profiles summed. Both contain five common PBMC cell types (B, CD4 T, CD8 K, NK and monocytes). A third independent dataset, STROMAL, contains 100 simulated mixtures of stromal cell types (B, CD4 T, CD8 T, macrophage, mast, endothelial and fibroblast cells). We use the first two for testing and the third for hyperparameter tuning.

In addition, we retrieved a real, GSE65133, dataset from [12]. This dataset contains 20 samples of real PBMC cell fractions that were measured by flow cytometry.

To complement the expression datasets, we used known signature matrices from [2]. This study contains a comparative analysis of 9 deconvolution methods with respect to 10 signature matrices. Among the top performing methods were CIBERSORT, NNLS and GEDIT which we use for comparison. We averaged the accuracy values reported in Figure 1 of [2] for each of the signature matrices and selected the four signatures with the highest accuracy value: Lm22, Skin Signatures, HPCA-Blood, and BlueCode, which are focused in our study. In detail, LM22 [12] contains 22 cell types and 547 genes; Human Primary Cell Atlas (HPCA-Blood) [17] contains 7 cells and 19,715 genes; Blue-Code [18] contains 34 cells and 13,299 genes; and Skin Signatures [19] contains 21 cells and 20,307 genes.

We preprocess the data using GEDIT’s approach [13] which removes cell types that are not present in either the input expression or signature matrix, runs quantile normalization on both matrices—such that each column follows the same distribution as every other, removes genes that are missing from either matrix, and selects a subset of 50 genes with lowest entropy for each cell to focus on (for each cell we want the genes that are expressed in a cell type-specific manner. Entropy is minimized when expression is detected only in a single cell type).

2.7.

Algorithm comparison

We used several performance measures to compare DECODE to four existing cell deconvolution algorithms: CIBERSORTX, NNLS, GEDIT and SCADEN. We ran GEDIT with its R source code. We ran CIBERSORTX from its official website (https://cibersortx.stanford.edu/). We ran NNLS with its R function. We ran Scaden with its Python source code and kept its default training datasets (as it does not train with a signature matrix). To compare the performance of the five deconvolution algorithms, we measured both RMSE (root mean squared error) and Pearson correlation coefficient, comparing real and predicted cell fractions estimates.