{'fraction_non_alphanumeric': 0.03634164299989481, 'fraction_numerical': 0.018775639002840013, 'mean_word_length': 4.79357099329677, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Imputing single-cell RNA-seq data by combining graph convolution and autoencoder neural networks

Summary Single-cell RNA sequencing technology promotes the profiling of single-cell transcriptomes at an unprecedented throughput and resolution. However, in scRNA-seq studies, only a low amount of sequenced mRNA in each cell leads to missing detection for a portion of mRNA molecules, i.e. the dropout problem which hinders various downstream analyses. Therefore, it is necessary to develop robust and effective imputation methods for the increasing scRNA-seq data. In this study, we have developed an imputation method (GraphSCI) to impute the dropout events in scRNA-seq data based on the graph convolution networks. Extensive experiments demonstrated that GraphSCI outperforms other state-of-the-art methods for imputation on both simulated and real scRNA-seq data. Meanwhile, GraphSCI is able to accurately infer gene-to-gene relationships and the inferred gene-to-gene relationships could also provide powerful assistance for imputation dynamically during the training process, which is a key promotion of GraphSCI compared with other imputation algorithms.


INTRODUCTION
Compared to bulk cell RNA sequencing (Wang et al., 2009) (RNA-seq), single-cell RNA sequencing technology (Kolodziejczyk et al., 2015) (scRNA-seq) has greatly promoted the profiling of transcriptomes at single-cell level and helped researchers to improve understanding of complex biological questions. It allows people to study cell-to-cell variability at a much higher throughput and resolution, such as studies of cell heterogeneity, differentiation and developmental trajectories (Saliba et al., 2014).
Despite its improvements, various technical deviations occurred due to the upgrade of sequencing techniques from bulk samples to single cells. Typically, the low RNA capture rate and sequencing efficiency lead to a large proportion of expressed genes with false zero counts in some cells, defined as 'dropout' event (Svensson et al., 2017;Kharchenko et al., 2014). For example, protocols based on droplet microfluidics (Zilionis et al., 2017) and Fluidigm C1 platform usually have a high dropout rate in the scRNA-seq data due to their technical limitations. And new droplet-based protocols, such as inDrop (Klein et al., 2015) and 10X Genomics (Zheng et al., 2017), have improved the detection rates but still have relatively low sensitivity, leading to the dropout events. On the other hand, although many of the zero counts represent the true absence of gene expression in specific cells, a considerable fraction is due to the dropout phenomenon where a truly expressed gene is undetected in some cells, resulting in zero or low read counts. Therefore, it is important to note the distinction between the truly expressed zeros and the false zeros in statistical analysis. Not all zeros can be considered as the missing values to be imputed. Imputation methods should impute the non-zero space but preserve the true-zero expression.
As a result, methods such as MAGIC (Van Dijk et al., 2018), SAVER (Huang et al., 2018), scImpute (Li and Li, 2018), scVI (Lopez et al., 2018), DCA (Eraslan et al., 2019), and DeepImpute (Arisdakessian et al., 2019) have been developed to correct the false zero read counts in order to recover true expression levels in scRNAseq data. These approaches estimate ''corrected'' gene expressions by borrowing information across similar genes or cells. For example, MAGIC imputes gene expression data for each gene across similar cells based on Markov transition matrix, while SAVER takes advantage of gene-to-gene relationships by using Bayesian approach to infer the denoised expression. Both MAGIC and SAVER would recover the expression level of each gene in each cell including those unaffected by dropout events. ScImpute, on the other hand, determines the dropout entries based on a mixture model and imputes only the likely dropout entries across similar cells. However, MAGIC and SAVER fail to learn the non-linear relationships and the counting structures in the scRNA-seq data. Thus, with the development of deep learning, neural network-based imputation methods have been proposed such as SAVER-X (Wang et al., 2019), DCA and DeepImpute. By combining a deep autoencoder with a Bayesian model, SAVER-X extracts transferable geneÀgene relationships to impute scRNA-seq data sets. DCA proposes an imputation method based on an Autoencoder, a kind of deep neural networks used to reconstruct data in an unsupervised manner. DeepImpute, in another way, constructs multiple sub-neural networks to impute genes in a divide-andconquer approach, which utilizes dropout layers and loss functions to learn patterns in the data. In the imputation of each cell, DCA minimizes the zero-inflated negative binomial (ZINB) (Risso et al., 2018) model-based loss function to learn gene-specific distribution in scRNA-seq data.
However, these existing imputation methods for scRNA-seq aim at learning the similarity of cells or genes but not considering gene-to-gene relationships and cell-to-cell correlations simultaneously, resulting in the fact that they cannot retain biological variation across cells or genes. And decades of molecular biology research have taught us much about the principles of gene interaction and their influence on gene expression (Bhardwaj and Lu, 2005;Fraser et al., 2004). For example, the gene is truly not expressed due to gene regulation, but imputed by similar cells, which makes it difficult to study cell-to-cell variation and downstream analysis. This means that our imputation method not only needs to take advantage of the information between similar cells but also gene-to-gene relationships. More importantly, as imputation proceeds, the imputed gene expression matrix could infer more accurate gene-to-gene relationships while the inferred gene-to-gene relationship helps improve the accuracy of imputation. Therefore, our imputation method needs to be able to dynamically integrate the imputation of gene expressions and inference of the gene-to-gene relationships during the training process. Accordingly, in this paper, we developed a Single-Cell Imputation method that combines Graph convolution network (GCN) and Autoencoder neural networks, called GraphSCI, to impute the dropout events in scRNA-seq by systematically integrating the gene expression with gene-to-gene relationships. We will use gene-to-gene relationships as prior knowledge to recover gene expression in a single cell because geneto-gene interactions are likely to affect gene expression sensitively. And the combination of GCN and autoencoder neural networks makes it possible for us to dynamically utilize the increasingly accurate gene-to-gene relationships to impute gene expressions. By stacking the GCN and autoencoder network, GraphSCI is capable of exploring the gene-to-gene relationships in an explicit way, so as to impute the dropout events effectively. Furthermore, the deep generative model with gene-specific distribution such as ZINB and NB distribution could learn the true data distribution of scRNA-seq data and then impute the dropout events and avoid overfitting.
The gene-to-gene relationships can be regarded as a gene graph, in which the gene is the node and the edge is the relationship. As a consequence, the imputation task of gene expression can be converted into the node recovering problem on graphs. GCN (Kipf and Welling, 2017) is a very powerful neural network architecture for machine learning on graphs. It was designed to learn hidden layer representations that encode both local graph structure and features of nodes and edges. A number of recent studies describe applications of GCN such as node recovering problem (Meng et al., 2019;Gong et al., 2014;Chakrabarti et al., 2014;Yang et al., 2017). Inspired by the co-embedding attributed network (Meng et al., 2019), we combine GCN and autoencoder neural network to systematically learn the low-dimensional embedded representations of genes and cells. GCN exploits the spatial feature of gene-to-gene relationships effectively while Autoencoder neural network learns the non-linear relationships of cells and counting structures of scRNA-seq data, and thus the deep learning framework reconstructs gene expressions by integrating gene expressions and gene-to-gene relationships dynamically in the backward propagation of neural networks.
Our proposed method was shown to outperform competing methods over both simulated and real data sets by diverse downstream analyses. To assess the performance of the imputation methods, we evaluate their improvement on several downstream analyses. Firstly, we perform cell clustering and use clustering metrics to demonstrate their effectiveness to impute the dropout events. And then we also perform the differential expression analysis to evaluate their improvement of the identification of differentially expressed genes (DEGs). The evaluation performance illustrates the rationality and effectiveness of our iScience 24, 102393, May 21, 2021 iScience Article proposed method. Furthermore, our method takes advantage of the gene-to-gene relationships in the framework that infers new reliable relationships simultaneously. Altogether, we demonstrate that our proposed method is highly scalable and parallelizable via graphical processing units (GPUs).

Overview of the GraphSCI algorithm
GraphSCI is a deep neural network model that combines the GCN and autoencoder neural network to impute gene expression levels in scRNA-seq data. The overview of our method is shown in Figure 1 and the detailed model architecture is shown in Figure 2.
Usually, the expressions of genes are correlated by their related genes or interacting genes because the coexpressed genes are controlled by the same transcriptional regulatory program, functionally related, or members of the same pathway or protein complex (Weirauch, 2011). Therefore, given the log-normalized expression data X, we first construct gene co-expression networks, called Gene Graph, from the raw expression data through the Pearson correlation coefficient (PCC). When the PCC between two genes is greater than 0.3 or less than À0.3, we assume that the two genes are co-expressed and there is an initial edge between them in the gene network. Obviously, our initial gene co-expression networks have a high rate of false positives because of dropout events in scRNA-seq data. We therefore combine the GCN and autoencoder neural network (AE) to dynamically integrate the imputation of gene expression and the inference of the gene co-expressed network where GCN encodes the gene co-expressed network with expression levels to the latent vector Z and then reconstructs the edges in gene co-expression network. AE encodes the gene expression matrix with gene co-expression network and finally sample Z from ZINB or NB distributions to reconstruct gene expression matrix.
This model enables us to utilize gene-to-gene relationships to impute the dropout events and further refine the gene co-expression network. The gene-to-gene network is an undirected graph, where each node corresponds to a gene and each edge between two genes indicates there is a significant co-expression relationship between them (Stuart et al., 2003). And the excellent characteristics of GCN allow us to regard the gene expression levels in different samples as node (gene) features in gene co-expression network and utilize them in the learning of gene network.

GraphSCI identifies cell types in simulated data
In order to assess our method, we followed the same way as the previous study (Eraslan et al., 2019) to construct two simulated data sets by Splatter (Zappia et al., 2017) package: (1) 2000 cells belonging to two types clustered by expression data of 3000 genes (namely SIM-T2) and (2) 3000 cells belonging to six types of cells clustered by expression data of 5000 genes (namely SIM-T6). On the SIM-T2 data with a simpler case, GraphSCI achieved imputed expression with a mean absolute error of 0.226, which is 21.2% lower than 0.274 by DCA. We further projected the imputed gene expression by t-SNE and clustered  Table S1. By comparison, the clustering over original expression data without any imputation achieved 0.716, 0.508, 0.342 for ARI, CA and SC. Figure 3B shows 2000 cells by using the first two principle components obtained from t-SNE (Maaten and Hinton, 2008). Obviously, GraphSCI clearly separates two types of cells, while both DeepImpute and DCA have a small number of cells mixed together. The original data can't separate the cells at all.
When tested on the SIM-T6 data set with six cell groups, similar results were obtained. As shown in Figure S1A, our method achieved 0.818, 0.859, and 0.34 for ARI, CA, and SC values, respectively. These are 5.1, 3.2, and 16.4% higher than those by DeepImpute, and 6.6%, 2.0%, and 19.7% higher than DCA. Table  S1 details results by the raw data, SAVER-X and scImpute. The visualization consistently indicated that our method separates the six types of cells better than other methods ( Figure S1B). Figure S2 shows the image of gene expression matrix (X) before and after imputation ( b X ) in our simulated experiments. This comparison again demonstrates that GraphSCI could recover the original cell types effectively both in the Sim-T2 and the Sim-T6 data sets.

GraphSCI recovers transcriptome dynamics in real single-cell data
Another key criterion to evaluate the imputation methods is their ability to recover transcriptome dynamics in real single-cell data set. Therefore, we applied our method to three real scRNA-seq data sets and made comparisons with other methods. The first data set was obtained from mouse ES cells (Klein et al., 2015), which were measured to analyze the heterogeneity of mouse embryonic stem cells in different stages after leukemia inhibitory factor (LIF) withdrawal. We selected four different LIF withdrawal intervals (0, 2, 4, and 7 days) and put all cells together as the input of imputation. The imputed data were clustered by t-SNE. As shown by Figure 4A, GraphSCI separated the four stages of mouse ES cells clearly except that a few blue samples were mixed with the yellow. In comparison, the clustering obtained from the scImpute and DCA methods seriously mixed the blue samples with the yellow ones. As indicated by Figure 4B, ARI, CA, and SC of GraphSCI were significantly higher than DeepImpute and DCA. Results of ARI, CA, and SC with all methods are detailed in Table S2.
GraphSCI was further applied to two large data sets generated by the 10X scRNA-seq platform (Zheng et al., 2017), one of which is involved by the transcriptome of peripheral blood mononuclear cells (PBMCs) from a healthy donor. The data set contains 5247 PBMCs of 11 cell types. Because the same type of cells has similar expression profiles, we randomly selected 80% of PBMCs to train the model and used the remained for the independent test set. The imputed data on the independent set was conducted with dimension reduction results by t-SNE for visualizations. Figure 5A shows that the imputations by GraphSCI could iScience Article separate brown and dark samples well in the low-dimension representation. The orange samples had a diving line with other samples. In comparison, the results obtained by DeepImpute and DCA didn't show obvious differences among the black, red, and green samples. The results of ARI, CA, and SC on the independent test set also showed that GraphSCI outperformed other methods ( Figure 5B). In details, our method achieved 0.472, 0.552, and 0.177 for ARI, CA, and SC values, respectively. These are 1.7%, 0.7%, 4.7% and 14.0, 9.7, 34.1% greater than those by DeepImpute and DCA respectively. More detailed results are shown in Table S2.
In the E18 Mouse data set, $12,000 brain cells of 16 cell types were profiled from the 10X scRNA-seq platform. We applied GraphSCI to the large-scale scRNA-seq data set to demonstrate its robustness and scalability. As shown in Figure S3A, GraphSCI is able to separate cells of 16 cell types effectively in the low-dimension representation, while DeepImpute and DCA have mixed many subcellular types together. GraphSCI again achieved ARI, CA, and SC of 0.316, 0.422, and 0.030, respectively, consistently the greatest among all methods ( Figure S3B and Table S2).

GraphSCI recovers gene expression levels in bulk RNA-seq data set
The efficacies of GraphSCI in recovering gene expression levels were further evaluated by a real RNA sequencing data set. The RNA sequencing data was obtained from C. elegans development experiments by Francesconi et al. (Francesconi and Lehner, 2014), which was used to simulate single-cell RNA-seq data with dropout rates ranging from 50% to 70%. The three data sets were generated by adding the single-cell specific noises through gene-wise subtracting values drawn from the exponential distribution. Since bulk iScience Article RNA-seq data contains less noise than scRNA-seq, PCC was used to evaluate the effectiveness of imputation on real RNA-seq data set. As shown by Figure 6, GraphSCI outperformed DCA in recovering the gene expression levels in real RNA-seq data set. In details, the median of PCCs reached by GraphSCI in three data sets are 0.858, 0.807, and 0.788 respectively, consistently greater than 0.821, 0.765, and 0.742 achieved by DCA, and 0.787, 0.718, and 0.604 achieved by SAVER.

GraphSCI improves differential expression analysis on scRNA-seq data
An effective imputation method should lead to an improvement in differential expression analysis because scRNA-seq data provide insight into gene expression in a single cell. To evaluate whether the identification of DEGs are more accurate after imputation, we utilized a scRNA-seq data set with corresponding bulk RNA-seq data to compared differential expression analysis results using DESeq2. This data set, generated by Chu et al. from H1 human embryonic stem cells (H1) differentiated into definitive endoderm cells (DEC), has six samples of bulk RNA-seq and 350 samples of scRNA-seq (212 for H1 ESC and 138 for DEC). We applied GraphSCI and DeepImpute to impute the gene expression on scRNA-seq data and performed DE analysis on the raw data and the imputed data respectively. The percentages of zero gene expression are 49.1% in raw scRNA-seq data that results in the lowest DEGs identification results. In contrast, GraphSCI and DeepImpute have improved the identification of DEGs and share more DEGs with bulk samples. In Figure 7A, we defined more quantitative evaluation metrics such as the area under the receiver operating characteristic curve (AUC), the accuracy (ACC), and F-scores for DEGs detection. In detail, our method achieved 0.913, 0.782, and 0.608 for AUC, ACC, and F-score values, respectively. Moreover, Figure 7B and Figure 7C show that the expression profiles of DEC and ESC marker genes (SOX2 and LEFTY1) after GraphSCI imputation could better reflect the gene expression signatures on recovering the expression patterns of signature genes. The performance of other signature genes including NANOG, DNMT3B, GATA6, and CXCR4 et al. has been shown in Figure S4. iScience Article GraphSCI infers gene-to-gene relationships from scRNA-seq data GraphSCI can not only impute gene expression data of scRNA-seq effectively, but also infer gene-to-gene relationships from the data. Due to the dropout events in raw scRNA-seq data (Zheng et al., 2017;Iacono et al., 2019), it is challenging to obtain accurate gene interactions directly from correlation coefficients between gene expression (Aibar et al., 2017). Here, we applied our method to raw scRNA-seq data sets and reconstructed gene relations during imputation. By compared with the known interactions from the STRINGdb (Szklarczyk et al., 2017), the gene interactions constructed by GraphSCI had a precision of 0.713 with the threshold of 0.5. Specifically, the true positive (TP) was 232647 and the false positive (FP) was 93,812. Figure 8 shows the imputed gene-to-gene relationships obtained by Cytoscape (Smoot et al., 2010). The true positive (TP) is an outcome where the model correctly predicts the gene-gene relationships and the false positive is an outcome where the model incorrectly predicts the gene-gene relationships. They showed the accuracy of our model to infer the gene-gene relationships from the raw scRNA-seq data. Similar results were also observed in previous experiments on the mouse ES cells data set. The constructed gene relations had a precision of 0.682 with 199291 true positives and 92,924 false positives. As a comparison, we utilized the PCC to infer gene-gene relationships from the raw scRNA-seq data, which obtains the precision of 0.598 and 0.492 respectively. It again verifies the effectiveness of our method, empirically showing that it facilitates the inference of the gene-to-gene relationships during the training process.

DISCUSSION
In this study, we presented an imputation method, GraphSCI, based on GCNs, which are particularly suitable for single-cell RNA-seq data. Our method focused on imputing gene expression levels by integrating the gene expression with gene-to-gene relationships. By using gene-to-gene relationships as prior To our best knowledge, this is the first study to integrate gene-to-gene relationships into deep learning framework for imputations on scRNA-seq. It is also the first attempt to employ GCN for learning the representation of gene-to-gene relationships in imputation study. Most importantly, extensive experiments were conducted on different kind of scRNA-seq data sets to demonstrate the superiority of our method.
GraphSCI was evaluated on both simulated and real data, which was identified with the best performances on diverse downstream analyses in comparison with other methods. In simulated data sets, GraphSCI was found to outperform other methods on the data with both small and large numbers of cells and cell types. The better performance of GraphSCI was further observed when it was applied in real data sets like bulk RNA-seq data and scRNA-seq data. In addition, another advantage of GraphSCI is its ability to infer the new gene-to-gene relationships, which is an absence of existing methods.
Applications of GCNs and exploiting gene-to-gene relationships for imputation, however, may also bring uncontrollable errors. For instance, the reliability of gene-to-gene relationships may influence the results of imputation. To solve this problem, we tried a variety of methods to build the gene-to-gene relationships, such as setting different thresholds to build edges or selecting original co-expressed samples to calculate PCC. We found that better performance could be achieved with the adjacency matrix obtained by selecting original co-expressed samples and PCC of >0.3 to determine edges. Figures S5 and S6 show the influence of the input gene-to-gene relationships on the overall results.
Another challenge for real data is that the evaluation of imputation may be difficult due to lack of ground truth. Therefore, we performed many clustering metrics, such as ARI, CA, and SC, to describe the effectiveness and robustness of competing methods, while we also utilized visualization to make the results clearer and more convincing. As shown in Figure S7, we could find that t-SNE showed better display results and GraphSCI consistently yields better performance with different clustering approaches.
The current GraphSCI was tested on data sets including simulated data and real data. The imputation power could be further improved with the increasing number of cells in the training set. Additionally, the deep learning networks by GraphSCI enable parallelization using GPUs to speed up training on large scRNA-seq data sets ( Figure S8).

Resource availability Lead contact
Yuedong Yang (yangyd25@mail.sysu.edu.cn) is the lead contact for this work. iScience Article

Materials availability
This study does not generate any new materials.

Data and code availability
The scRNA-seq data sets used in this manuscript are publicly available and their details are summarized in Table S3. The C. elegans time course experimental data was provided by the supplementary material of Francesconi. et al. The mouse embryonic stem cells data was downloaded from GSE65525. The 5k PBMC from a healthy donor and 10K brain cells from an E18 Mouse were provided by the 10X scRNA-seq platform and the website of the data are https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.1.0/5k_pbmc_protein_v3 and https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.0.0/neuron_10k_v3. The human embryos cells scRNA-seq data was downloaded from GSE44183. The Human ESC scRNA-seq data set for differential expression analysis was downloaded from GSE75748. The code generated during this study is available at https://github.com/biomed-AI/GraphSCI. We tuned model hyper-parameters based on the experimental results on simulated data sets and used them across all data sets ( Figures S9 and S10).

METHODS
All methods can be found in the accompanying transparent methods supplemental file.

ACKNOWLEDGMENTS
This study has been supported by the National Key R&D Program of China (2020YFB0204803), National Natural Science Foundation of China (61772566), Guangdong Key Field R&D Plan (2019B020228001 and       From the visualization and the clustering performance, we could find that the gene regulatory network inference tools such as PIDC could facilitates the imputation of scRNA-seq data using GraphSCI. We attribute the remarkable improvement to the accuracy of the input gene-to-gene relationships. We could find that t-SNE showed better display results with closer inner-group distance and larger between-group distances. (B) We compared different clustering approaches (PCAreduce, SC3 and KMeans) through the clustering performance (ARI). We observed that GraphSCI consistently yields better performance with different clustering approaches, showing that our algorithm could achieve stable and better results under the same conditions. It again illustrates the rationality and effectiveness of our algorithm.

Figure S8. The runtimes for imputation with different numbers of cells down-sampled from 1.3 million mouse brain cells, Related to Figure 1-2.
We analyzed the largest scRNA-seq data set in our experiments, which consists of 1.3 million mouse brain cells from 10X Genomics. The 1.3 million cell data matrix was down-sampled to 1,000, 2,000, 5,000, 10,000 and 100,000 cells and each subsampled matrix was imputed, and the runtime measured. We could find that the runtime of DCA and GraphSCI scaled linearly with the number of cells and the other methods took hours to impute 100,000 cells. It makes sense since DCA and GraphSCI are the neural network-based method that could be accelerated by GPU and the other methods failed to run due to the memory limitations on the large dataset.

Transparent Methods
The proposed model GraphSCI imputes gene expression levels in scRNA-seq data based on a combination of the graph convolution network and Autoencoder neural network, with the input of gene expression matrix and gene-to-gene relationships . In our framework, GCN encodes the gene-togene network with expression matrix to the latent vector and then reconstructs the edges in geneto-gene network. AE encodes the gene expression matrix with gene-to-gene network and finally sample from ZINB or NB distributions to reconstruct gene expression matrix.
By using single cells RNA-seq data with genes, an undirected gene graph with gene expressions and gene-to-gene relationships can be constructed. Let and ℳ be a set of genes and samples respectively, an undirected gene graph can be denoted as = ( , ℳ, ℰ), where ℰ is the set of gene-to-gene relationships. Thus, we introduce an adjacency matrix ∈ ℝ W×W and a gene expression matrix ∈ ℝ W×Z for , with _` representing the edge of the -th gene and the -th gene and _` being the expression value with rows representing genes and columns representing cells. Table S4 summarizes our main notations for scRNA-seq data.
Data processing and normalization. There are two inputs to our proposed model: (1) a gene expression matrix ∈ ℝ W×Z , (2) an adjacency matrix ∈ ℝ W×W , and our final goal is to construct an imputed gene expression matrix $ with the same dimensions. First, in raw scRNA-seq read count matrix Y , genes with no reads in any cell would be filtered out. Then, the library size of cell is denoted as _ and is calculated as the total number of read counts of cell . The size factor _ of cell is _ divide by the median of total counts per cell. Therefore, we make a normalized matrix by taking the log transformation with a pseudo count and scale of the read counts: where = 1,2, … , representing each gene and = 1,2, … , representing each sample.
Secondly, we attempt to obtain the adjacency matrix ∈ ℝ W×W from a graph where genes are nodes and edges indicate genes which are likely to be co-expressed. For the simulated datasets generated from Splatter (Zappia et al., 2017), we introduce the adjacency matrix ∈ ℝ W×W by Pearson correlation coefficient (PCC) as: where ( , ) and g is the covariance between and and the standard deviation of respectively.
Imputation based on graph convolution network. The preprocessed gene expression matrix and adjacency matrix are treated as the input for GraphSCI. Two neural network models, i.e., the inference model … and the generative model ‡ were used to constructed the model for the probabilistic encoder … and probabilistic decoder ‡ respectively, to preform gradient descent for learning all trainable parameters.
To infer the embeddings of cells and genes, we apply a two-layer graph convolution network and a two-layer fully connected neural network mapping the adjacency matrix A and the gene expression matrix X to the low-dimensional representations of the posterior distribution (i.e. Gaussian distributions and ZINB distributions) respectively. In particular, the two-layer GCN is defined as: In particularly, the ZINB distribution is applied for count data that exhibit over-dispersion and excess zeros, which is parameterized with the mean ( ) and dispersion ( ) of the negative binomial distribution as well as the dropout probability ( ) representing the probability of zeros (dropout events). But dropletbased scRNA-seq (such as 10X) are supposed to follow a NB distribution. A count matrix X that is ZINB-distributed with ( , , ) or NB-distributed with ( , ) are denoted as: where Γ(x) and ' ( ) is the Gamma function and Dirac function respectively. Therefore, we could estimate the parameters , , of ZINB distribution from the hidden layer in Eq. (6): where sigmoid(⋅) is the sigmoid function.
Optimization. The optimization was performed to obtain accurate embeddings of both genes and cells in an unsupervised way. For this purpose, and ℳ were optimized by the variational lower bound ℒ: where ¿ À is the cross entropy function with the probabilistic distribution … and ‡ and is the Kullback-Leibler (KL) divergence between q(⋅) and p(⋅). In the above equation, … ( ℳ | , ) and … ( | , ) is defined as the probabilistic encoder with the input of , and , respectively, aiming at producing the representations , ℳ . Similarly, ‡ o r $ _`· _ ,`ℳs and ‡ oe r ] _`· _ ,`s are the probabilistic decoders for construct the imputed gene expression matrix $ and gene-to-gene relationships ] . Furthermore, the KL divergence in optimization could be interpreted as the regularization to make the predicted posterior distributions closer to the prior distributions ( ℳ ), ( ).
With the help of reparameterization trick, we could represent the distributions with deterministic variables: Therefore, with the optimization, the gradient-based optimization techniques can be used to train the end-to-end model.

Evaluation metrics.
To evaluate the accuracy of imputation, we examine the reconstruction accuracy and clustering performance to the scRNA-seq datasets. The reconstruction accuracy on the simulated dataset can be measured by mean absolute error (MAE), which is the reconstruction error between the true expression matrix and imputed matrix. Clustering performance can be measured by the clustering metrics: adjusted Rand index (ARI) (Hubert and Arabie, 1985), clustering accuracy (CA) and Silhouette Coefficient (Rousseeuw, 1987) (SC). To fairly quantitate the performance of differentially expressed genes (DEGs) detection using scRNA-seq data, we calculated the accuracy (ACC), F-score and AUC for each DEG detection.
The adjusted Rand index (ARI) is the corrected-for-chance version of the Rand index. The Rand index is a measure of the similarity between two data clustering and the ARI is adjusted for the chance grouping of elements. Given a set of n samples, the two clusters of these samples are = { ‹ , • , … , Ï } and = { ‹ , • , … , C } and _` is defined as _`= | _ ∩`|.
The CA is defined as the accuracy of the clustering assignments. Given a sample , let _ be the ground truth label and _ be the assignments of clustering, then the CA is where is the number of samples, is the set of one-by-one mapping between clustering assignments and true labels and (x, y) = 1 if x = y otherwise 0.
The SC measured the similarity between a single cell and its cluster. The silhouette ranges from −1 to +1, where a high value indicates that the object is well matched to its own cluster. It could be defined as In the experiments of differential expression analysis, we took the DEG detection as the problem of predicting a gene is DEG or not, and the gold standard are obtained from bulk RNA-seq. Therefore, the accuracy (ACC), F-score and AUC could be calculated by: The F-score is calculated from the precision and recall of the DEG predictions, where the precision is the number of correctly detected genes divided by the number of all DEGs and the recall is the number of correctly detected genes divided by the number of all DEGs that should have been detected. It could be defined as: where TP is the true positive meaning that the correct DEG have been detected, FP is the false positives and FN is the false negatives.
AUC -ROC curve is a performance measurement for the classification problems at various threshold settings, which could be applied to evaluate the detection of DEGs. The AUC is calculated by the area under the ROC-curve, which represents the degree or measure of separability.
Simulated datasets. Our simulated data are generated by Splatter (Zappia et al., 2017) R package, a widely used package for simulating the scRNA-seq count data. First, we simulated a dataset with two cell groups, 2000 cells of 3000 genes by setting 27% of data values to zero mimicking dropout events.
During the simulation, we set the parameter . ℎ = −1 , . = 0 and . = 0.3 for simulating the dropout events and the other parameters are set to default values.
Hence, we could obtain the true counts before dropout and the raw counts after dropout, which are the simulated scRNA-seq data. Furthermore, we simulated a complex dataset of 3000 cells by 5000 genes to evaluate the robustness of our model, The 3000 cells are divided into six groups and the parameter were elegans development (Francesconi and Lehner, 2014). We analyzed the dataset after simulating singlecell transcriptomics dropout noises and the bulk transcriptomics data can be the ground truth for evaluation. Hence, we compared our method with the existing method DCA (Eraslan et al., 2019) by Pearson correlation coefficient.
Mouse embryonic stem cells data. Klein. et al. profiled the single-cell transcriptomics by dropletmicrofluidic approach and applied it on embryonic stem cells (Klein et al., 2015). They analyzed the heterogeneity of mouse embryonic stem cells differentiation after leukemia inhibitory factor (LIF) withdrawal. Here, we selected the four different LIF withdrawal intervals (0, 2, 4, 7 days) and construct a scRNA-seq dataset with 2717 cells of 24175 detected genes. And the cell types are determined by the intervals of LIF withdrawal.
Human ESC scRNA-seq dataset for differential expression analysis. Chu et al generated bulk and scRNA-seq data from H1 human embryonic stem cells (H1) differentiated into definitive endoderm cells (DEC). This dataset contains six samples of bulk RNA-seq (four for H1 ESC and two for DEC) and scRNA-seq of 350 single cells (212 for H1 ESC and 138 for DEC). The percentage of zero expression is 14.8% for the bulk RNA-Seq dataset and 49.1% for the scRNA-Seq dataset.
5k peripheral blood mononuclear cells (PBMC) from a healthy donor. The dataset was provided by 10X scRNA-seq platform (Zheng et al., 2017), profiling the transcriptome of the peripheral blood mononuclear cells (PBMCs) from a healthy donor. The total number of cells was 5247 after filtering process and the cell types were identified by graph-based clustering on the platform. 10K Brain Cells from an E18 Mouse dataset. The dataset was also provided by 10X scRNA-seq platform, profiling the brain cells from a combined cortex, hippocampus and sub ventricular zone of an E18 mouse. We could obtain the dataset containing 11843 mouse brain cells of 31053 detected genes and the cell types were identified by graph-based clustering on the platform.
Human Embryos cells scRNA-seq data. Xue et al. performed a comprehensive analysis of transcriptome dynamics by weighted gene co-expression network analysis (Xue et al., 2013). Therefore, we could obtain the dataset containing 30 samples from oocyte to morula in human embryos samples from their experiments. Here, we utilized the dataset to demonstrate the effectiveness of our method on inferring the gene-to-gene relationships.
In the training process, we utilized the Adam (Kingma and Ba, 2014) optimizer with an initial learning rate of 0.01 and allowed it to decay exponentially with decay_rate = 0.9 and decay_steps = 50