LIDER: cell embedding based deep neural network classifier for supervised cell type identification

Datasets

Table 1 illustrated the distribution of cell types from single-cell RNA sequencing datasets for model evaluation of LIDER. We downloaded the gene expression matrix and cell type annotation of mouse cortical data from the Heisberg Group scRNA-seq dataset (https://hemberg-lab.github.io/scRNA.seq.datasets/) (Zeisel et al., 2015), those of pancreatic islets dataset from ArrayExpress (https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-5061/)  (Segerstolpe et al., 2016), those from GSE71585 (Tasic et al., 2018), those of a PBMC dataset (https://support.10xgenomics.com/single-cell-gene-expression/datasets/2.1.0/pbmc4k) (Zheng et al., 2017), those of a human liver dataset from GSE115469 (MacParland et al., 2018), those from GSE36552 (Yan et al., 2013), those of a mouse bladder cell dataset from the Mouse Cell Atlas project (Han et al., 2018), those of colorectal cancer (CRC) from GSE81861 (Li et al., 2017), those of peripheral bipolar cell from GSE118480 (Peng et al., 2019). The mouse cell atlas datasets Tabula Muris were collected from https://tabula-muris.ds.czbiohub.org. We performed the Z-score transformation $\overline{f}=\frac{f-E\left(f\right)}{\mathrm{Std}\left(f\right)}$ to standardize the gene expression values of each cell for each dataset, where f represented the expression value of each gene, E(f) denoted the mean of each f, and Std(f) represented the standard deviation of each f respectively. This process made the expression levels comparable across genes.

Overview of LIDER development and evaluation workflow

Figure 1 illustrates the overview of the supervised cell type identification for single-cell transcriptomic data. LIDER aims to improve the cell type identification of single cell data by using a combination of stacked denoising autoencoder and deep neural network classifier. The cell embeddings based on stacked denoising autoencoders are used for developing deep neural network classifier (Fig. 1A). The low-dimensional representation by training the first level denoising autoencoders (DAE1) is utilized to train the second level denoising autoencoders (DAE2) and this process is repeated with several times to develop finally stacked denoising autoencoders (SDAE) (Fig. 1B). The process of LIDER for cell type identification could be divided into three main steps: (i) the standardized single-cell transcriptomic data is used as an input of stacked denoising autoencoders; (ii) the stacked denoising autoencoders for identifying cell embedding; (iii) the development of deep neural network classifier based on a parameter optimization algorithm Adam.

Stacked denoising autoencoder for identifying cell embedding

Denoising autoencoder was developed to learn a robust representation for partial corruption of the input pattern (Vincent et al., 2008). In this analysis, we used denoising autoencoder as feature selection method for identifying the low-dimensional representation called cell embedding. For a given single-cell sample, $\overline{X}ϵR^{n\times d}$ represents the gene expression matrix where n is the number of single cells and d is the number of the genes. The output $H^{\left(\mathfrak{l}\right)}$ of encoding layer 𝔩 is calculated as: (1)${\displaystyle H^{\left(\mathfrak{l}\right)}=\varphi \left(W_e^{\left(\mathfrak{l}\right)}{H}^{\left(\mathfrak{l}-1\right)}+b_e^{\left(\mathfrak{l}\right)}\right)}$

where ϕ is an activation function, $W_e^{\left(\mathfrak{l}\right)}$ is the weight matrix, $b_e^{\left(\mathfrak{l}\right)}$ is bias parameters and $H^{\left(\mathfrak{l}-1\right)}$ represents the output of encoding layer (𝔩 − 1) respectively. Specifically, $H^{\left(0\right)}=\overline{X}+\varepsilon$ where ɛ is a given type of noise. We randomly mask some observed input features to zeros to develop the corrupted ones, and thus denoising autoencoders are supposed to learn more robust representations.

The decoded layers could be defined as: (2)${\displaystyle H^{\left(\mathfrak{l}\right)}=\varphi \left(W_d^{\left(\mathfrak{l}\right)}{H}^{\left(\mathfrak{l}-1\right)}+b_d^{\left(\mathfrak{l}\right)}\right)}$

where ϕ is an activation function, $W_d^{\left(\mathfrak{l}\right)}$ is the weight matrix, $b_d^{\left(\mathfrak{l}\right)}$ is bias parameters and $H^{\left(\mathfrak{l}-1\right)}$ represents the output of decoding layer (𝔩 − 1) respectively. The weight matrix and bias parameters are optimized according to the following loss function (3)${\displaystyle L=\frac{1}{n}\sum _{i=1}^n\sum _{j=1}^d{\left({\overline{X}}_{ij}-X_{ij}\right)}^2}$ where X is the reconstructed data as the output of the last layer. As a loss function, we use the mean squared error to implement the denoising autoencoder and optimize the parameters.

Stacked denoising autoencoders are usually developed by multiple stacking layers of denoising autoencoders. After training the first layer of denoising autoencoders, its resulting representation is usually used to train denoising autoencoders of the second level. This multi-layer cycle process develops stack multiple denoising autoencoders for the task of layer-by-layer feature extraction, and eventually makes the features more representative (Vincent et al., 2010). A stacked denoising autoencoder is used to learn encoder–decoder structures for identifying cell embedding. The activation function ReLU = max(0, x) is used for the encoder and decoder networks except for the bottleneck layer and last decoder layer, in which we use $tanh\left(x\right)=\frac{1-e^{-2x}}{1+e^{-2x}}$ as the activation function. The loss function of stacked denoising autoencoders is defined as the reconstruction loss function Eq. (3) between the input data $\overline{X}$ and the output data X. After training each layer by minimizing the reconstruction loss of each layer, all encoder layers are connected with all decoder layers in reverse layer-by-layer training order, in which all parameters are fine-tuned by minimizing the loss. After training the first level denoising autoencoders, the obtained representation is used to train the second level denoising autoencoders. This procedure could be repeated to develop stacked denoising autoencoders for learning and stacking several layers. The embedded features were obtained from the representation H^L of stacked denoising autoencoders (Fig. 1B).

Deep neural network classifier

Given ${\left\{\left(x_i,y_i\right)\right\}}_{i=1}^n$, ${\left\{x_i\right\}}_{i=1}^n$ represents the embedded features for cell i and ${\left\{y_i\right\}}_{i=1}^n$ indicates a cell type label of each cell, where n is the number of single cells. A neural network contains an input layer, two hidden layers and one output layer (Fig. 1A). The number of nodes in the input layer is equal to the number of derived features from stacked denoising autoencoder in a training set. The number of nodes in the output layer is equal to the number of cell types of single-cell transcriptomic data. The neural network could be described as follows (Chen et al., 2020a):

(4)${\displaystyle O_1=\mathrm{sigmoid}\left(W_{1}X+b_1\right)}$ (5)${\displaystyle O_m=\mathrm{sigmoid}\left(W_m{O}_{m-1}+b_m\right),2\le m\le M}$ (6)${\displaystyle \overline{y}=\mathrm{softmax}\left(W_{m+1}{O}_M+b_{m+1}\right)}$

where X represents the output from stacked denoising autoencoders, W_m is the weight matrix of the mth layer, b_m is the bias vector of the mth layer, O_m−1 is the output of the (m −1)th layer, m represents the layer of a neural network and $\overline{y}$ is the output of the classification layer. For the input and hidden layers, we use the sigmoid function as the activation function, which is defined as follows: (7)${\displaystyle \mathrm{sigmoid}\left(x\right)=\frac{1}{1+e^{-x}}}$

For the output layer we use the softmax activation function, which is defined as follows: (8)${\displaystyle \mathrm{softmax}\left(x_j\right)=\frac{exp\left(x_j\right)}{\sum _{j=1}^{J}exp\left(x_j\right)}}$

where J element indicates a total of J cell types in the training set. We define the cross-entropy function as classification loss function measuring the discrepancy between the predicted and true class label: (9)${\displaystyle L=-\sum _{i=1}^n\sum _{j=1}^J{y}_{ji}log{\overline{y}}_{ji}}$

where y_ji is the true label of cell type for x_i and ${\overline{y}}_{ji}$ represents the probability of x_i belongs to the j th cell type. We utilize a multilayer neural network to identify cell types as a supervised-learning problem, that is, to accurately predict the class labels of cell types.

Adam is a promising stochastic optimization algorithm for first-order gradient-based optimization with stochastic objective functions. This approach is invariant to diagonal rescaling of the gradients and computationally efficient with little memory requirements (Kingma & Ba, 2014). In this analysis, we use the Adam algorithm to develop deep neural network classifier. $f\left(\theta \right)$ is a stochastic objective function with respect to parameters $\theta =\left(W,b\right)$ and our aim is to minimize the expected value of the objective function according to the cross-entropy function Eq. (9).

Experimental setting and LIDER implementation

We implemented the stacked denoising autoencoders using Keras (https://github.com/fchollet/keras) and the code was written in Python. For all single-cell datasets, the encoder network of the stacked denoising autoencoders was set to size d-500-500-2000-1000 for the fully connected multilayer perceptron (MLP) where d represents a number of genes from single-cell transcriptomic data, while the decoder network was a MLP of size 1,000-2,000-500-500-d. In this analysis, the neurons hold 500, 500, 2,000 and 1,000 for each layer in the encoder network and 1,000, 2,000, 500 and 500 for each layer in the decoder network respectively. In the process of pretraining stacked denoising autoencoders, the drop ratio for constructing i th stacked denoising autoencoder was set to 0.2 where drop ratio represented the proportion of random disconnection for dropout, and each denoising autoencoder was trained for 200 epochs during the layer-by-layer pre-training period using the stochastic gradient descent (SGD) algorithm. After pre-training stacked denoising autoencoders, we then trained 400 epochs using the SGD algorithm for overall fine-tuning, and the learning rate was set to 0.1 with decayed 10 times every 80 epochs. In this analysis, the selected features of 1,000 dimension from single-cell data as low-dimensional representations were derived to develop a classifier for identifying cell types.

The deep neural network classifier was implemented by Python using TensorFlow in the process of model development (https://www.tensorflow.org). We designed a four-layer neural network with the number of nodes in the input layer (1000), two hidden layers (528 and 256 respectively), and the number of nodes in the output layer according to the number of categories of the single-dell data. We then set the learning rate to 1e−3, the training period for supervised model initialization to 1200, and the batch size to 256 respectively. Good default settings for the implementation of Adam algorithm are α = 0.001, β₁ = 0.9, β₂ = 0.999 and ɛ = 10⁻⁸ respectively. In our analysis, we randomly divided the single-cell transcriptomic data into training and test datasets, containing 80% and 20% of the samples respectively. The deep neural network classifier was trained on the training set using the optimal parameters, validated on the test dataset, and evaluated based on accuracy.

Baseline methods

We compared the performance of LIDER with five baseline methods, such as logistic regression multiclassification algorithm (LR), Moana (Wagner & Yanai, 2018), SingleCellNet (Tan & Cahan, 2019) and ACTINN (Ma & Pellegrini, 2020). Logistic regression multiclassification algorithm (LR) is developed by Python for further comparison. L2 regularization is used in logistic regression modeling. An implementation of Moana is made by Python and found at https://github.com/yanailab/moana. All codes about ACTINN are implemented in Python and available at https://github.com/mafeiyang/ACTINN. The SingleCellNet code by R is available from https://github.com/pcahan1/singleCellNet.

LIDER improves prediction performance for cell type identification

To test whether autoencoder based deep learning model could improve prediction performance, we developed LIDER by training classifiers and assessing their performance when applied to test cohorts. LIDER was used to classify the real single-cell transcriptomic data including Zeisel, Segerstolpe, Tasic, PBMC, MacParland, Yan, Mouse and CRC respectively. The prioritized features based on stacked denoising autoencoders was utilized for the training set and deep neural network classifier was developed to identify cell types of single-cell transcriptomic data. Based on the optimized parameters from training data set, LIDER was developed to classify the test data set of single-cell transcriptomic data. As shown in Fig. 2, the proposed LIDER model achieved an accuracy between 0.857 and 1.0 and an average value of accuracy 0.934 for eight types of single-cell datasets.

To further demonstrate the effectiveness of the method, we compared the proposed LIDER with logistic regression multiclassification algorithm (LR), Moana, SingleCellNet and ACTINN respectively. These four baseline models were trained on the training data set and validated on the test data set. LIDER achieved an accuracy in four test datasets (Figs. 3A–3D), which yielded 1.6%, 4.6%, 10% and 5.8% (Fig. 3A), 2.6%, 5.2%, 7.6% and 6.0% (Fig. 3B), 1.1%, 3.4%,11.7% and 5.2% (Fig. 3C), 1.5%, 2.9%, 5.1% and 4.8% (Fig. 3D) increase in accuracy as compared to LR, Moana, SingleCellNet and ACTINN respectively. An accuracy in four test datasets was yielded by LIDER (Figs. 3E–3H), which was 1.6%, 5.4%, 11.3% and 7.8% (Fig. 3E), 2.9%, 7.5%, 9.8% and 8.9% (Fig. 3F), 1.3%, 4.1%, 7.1% and 5.7% (Fig. 3G), 1.6%, 3.2%, 8.4% and 5.7% (Fig. 3H) higher than that obtained by LR, Moana, SingleCellNet and ACTINN respectively. Overall, independent test results clearly indicated that LIDER could perform better than a diverse panel of cell type identification methods. The performance demonstrated great potential in the deep supervised learning. LIDER converts gene-level expressions into an encoded lower-dimensional representation of single-cell RNA-seq data. The autoencoder removes noise and leaves a high-value representation from the input data. The classifier can perform better because the algorithm is able to learn the patterns in the data from a smaller set of a high-value input.

LIDER achieves similar cell type identification with the true identity

To investigate the validity of the proposed method LIDER, we compared it with true cell types for the Zeise, Segerstolpe, Tasic, PBMC, MacParland, Yan, and Mouse datasets respectively. We applied principal component analysis (PCA) and t-distributed stochastic neighbor embedding (tSNE) to reduce the dimension and plotted the data points on a two-dimensional plane. A graph-based clustering approach was used to cluster the cells for cell type identification. Figure 4 shows the 2D visualizations of the true cell types and the predicted labels from LIDER. In the identified cell embedding of the proposed method LIDER, cells with the same type are correctly separated with only few outliers (Fig. 4). The vast majority of the predicted cell type from LIDER is similar with the true identity of each cell from single-cell datasets. These results suggested that deep supervised learning could help improve prediction performance.

LIDER accurately identifies cell types by using stacked denoising autoencoder

In supervised learning methods, feature selection usually plays an important role in improving classification performance. Therefore, we explore whether stacked denoising autoencoders could improve the prediction performance for cell types identification. As a reference, we used principal component analysis (PCA) based neural network classifier where the architecture of the neural network was similar with the classifier development in LIDER. Specifically, PCA is a popular method for feature selection and shows great potential in dimension reduction of high dimensional data. The PCA algorithm was implemented by Python from PCA package in Sklearn. The performance accuracy is defined by $\mathrm{accuracy}=\frac{1}{n}{\sum }_{i=1}^n\delta \left(y_i^p,y_i\right)$, where δ = 1 for $y_i^p=y_i$ and δ = 0 for $y_i^p\ne y_i$ respectively. $y_i^p$ represents the predicted cell type of each single cell and y_i represents the true cell type of each single cell.

Figure 5 illustrated the accuracy between the proposed method LIDER and the PCA-based neural network classifier from the test dataset of eight types of single-cell transcriptomic data. LIDER achieved an accuracy from seven datasets, which was 7.6% (Fig. 5A), 6.3% (Fig. 5B), 1.6% (Fig. 5D), 1.2%(Fig. 5E), 3.8%(Fig. 5F), 5.7% (Fig. 5G), and1.2% (Fig. 5H) higher than that from the PCA-based neural network classifier respectively. Moreover, LIDER achieved an accuracy in the Tasic dataset, which was comparable with that obtained by the PCA-based neural network classifier (Fig. 5C). The results suggested that LIDER was able to identify cell embeddings from the high-dimensional single-cell transcriptomic data by using stacked denoising autoencoder and achieved significantly better performance. Compared with PCA, LIDER developed stack multiple denoising autoencoders which made the selected features more representative of single-cell RNA-seq data.

LIDER suggests comparable robust to batch effects

Although scRNA-seq studies shows great potential with the characteristic of unprecedented high resolution, a challenging task is to analyze scRNA-seq data with different batches from the generation at different times. Failure to remove batch effects may lead to a false explanation of true biological variations in downstream analysis. To evaluate the performance when reducing the impact of batch effects, we used the Macaque dataset containing animal-level batch effects for further analysis. We still chose 80% as the training set and 20% as the test set. LIDER is effective in removing batch effects and achieves a high accuracy 0.972 when validated on the test data set. To test performance of LIDER accounting for removing batch effects which are due to different techniques, we trained it on single-cell data developed by one platform and tested it on single-cell data developed by another platform. LIDER was trained on the 10X cells and then tested on the SS2 cells. It achieved a testing accuracy of 1.0. Among the seven incorrectly predicted cells, three B cells were predicted as hepatocytes, two monocytes were predicted as stromal cells, and there were remaining two mispredictions. We then trained the proposed method on the SS2 dataset and tested it on the 10X dataset. LIDER obtained a testing accuracy of 0.999. Among the thirty-three incorrectly predicted cells, four B cells were predicted as cardiac muscles, four B cells were predicted as epithelial cells, three epidermises were predicted as natural killer (NK) cells, three epidermises were predicted as hepatocytes, two cardiac muscles were predicted as epithelial cells, two cardiac muscles were predicted as T cells, two epidermises were predicted as epithelial cells, two epidermises were predicted as endothelial cells, two T cells were predicted as B cells, and there were also several mispredictions. The results showed that LIDER fairly removed batch effects from single-cell datasets introduced by different techniques.