scGNN is a novel graph neural network framework for single-cell RNA-Seq analyses — TEI-tagged view

Juexin WangpersName, Anjun MapersName, Yuzhou ChangpersName, Jianting GongpersName, Yuexu JiangpersName, Ren QipersName, Cankun WangpersName, Hongjun FupersName, Qin MapersName, Dong XupersName (2021date). “scGNN is a novel graph neural network framework for single-cell RNA-Seq analyses.” Nature Communicationstitle·j.

Object of study — single-cell RNA-sequencing data with sparsityterm·dom and complex expression patterns

Research question — whether a graph neural networkterm·mth over cell-cell relationships improves imputation and clustering

Contribution — scGNNterm·mth, a hypothesis-free GNNterm·mth framework with iterative multi-modal autoencodersterm·mth

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

scGNNterm·mth is a novel graph neural networkterm·mth framework for single-cell RNA-Seqterm·key analyses Abstract Single-cell RNA-sequencing (scRNA-Seq) is widely used to reveal the heterogeneity and dynamics of tissues, organisms, and complex diseases, but its analyses still suffer from multiple grand challenges, including the sequencing sparsityterm·dom and complex differential patterns in gene expression. We introduce the scGNNterm·mth (single-cell graph neural networkterm·mth) to provide a hypothesis-free deep learning framework for scRNA-Seq analyses. This framework formulates and aggregates cell–cell relationships with graph neural networksterm·mth and models heterogeneous gene expression patterns using a left-truncated mixture Gaussian model. scGNNterm·mth integrates three iterative multi-modal autoencodersterm·mth and outperforms existing tools for gene imputationterm·mth and cell clusteringterm·key on four benchmark scRNA-Seq datasets. In an Alzheimer’s diseaseterm·key study with 13,214 single nuclei from postmortem brain tissues, scGNNterm·mth successfully illustrated disease-related neural development and the differential mechanism. scGNNterm·mth provides an effective representation of gene expression and cell–cell relationships. It is also a powerful framework that can be applied to general scRNA-Seq analyses.

Introduction Single-cell RNA-sequencing (scRNA-seq) techniques enable transcriptome-wide gene expression measurement in individual cells, which are essential for identifying cell-type clusters, inferring the arrangement of cell populations according to trajectory topologies, and highlighting somatic clonal structures while characterizing cellular heterogeneity in complex diseases 1 , 2 . scRNA-seq analysis for biological inference remains challenging due to its complex and un-determined data distribution, which has a large volume and high rate of dropoutterm·dom events. Some pioneer methodologies, e.g., Phenograph 3 , MAGIC 4 , and Seurat 5 use a k-nearest-neighborterm·mth (KNNterm·mth) graph to model the relationships between cells. However, such a graph representation may over-simplify the complex cell and gene relationships of the global cell population. Recently, the emerging graph neural networkterm·mth (GNNterm·mth) has deconvoluted node relationships in a graph through neighbor information propagation in a deep learning architecture 6 – 8 . Compared with other autoencodersterm·mth used in the scRNA-Seq analysis 9 – 12 for revealing an effective representation of scRNA-Seq data via recreating its own input, the unique feature of graph autoencoderterm·mth is in being able to learn a low-dimensional representation of the graph topology and train node relationships in a global view of the whole graph 13 .

We introduce a multi-modal framework scGNNterm·mth (single-cell graph neural networkterm·mth) for modeling heterogeneous cell–cell relationships and their underlying complex gene expression patterns from scRNA-Seq. scGNNterm·mth trains low-dimensional feature vectors (i.e., embeddingterm·mth) to represent relationships among cells through topological abstraction based on both gene expression and transcriptional regulation information. There are three unique features in scGNNterm·mth: (i) scGNNterm·mth utilizes GNNterm·mth with multi-modal autoencodersterm·mth to formulate and aggregate cell–cell relationships, providing a hypothesis-free framework to derive biologically meaningful relationships. The framework does not need to assume any statistical distribution or relationships for gene expression data or dropoutterm·dom events. (ii) Cell-type-specific regulatory signals are modeled in building a cell graph, equipped with a left-truncated mixture Gaussian (LTMG) model for scRNA-Seq data 14 . This can improve the signal-to-noise ratio in terms of embeddingterm·mth biologically meaningful information. (iii) Bottom-up cell relationships are formulated from a dynamically pruned GNNterm·mth cell graph. The entire graph can be represented by pooling on learned graph embeddingterm·mth of all nodes in the graph. The graph embeddingterm·mth can be used as low-dimensional features with tolerance to noises for the preservation of topological relationships in the cell graph. The derived cell–cell relationships are adopted as regularizers in the autoencoderterm·mth training to recover gene expression values.

scGNNterm·mth has great potential in capturing biological cell–cell relationships in terms of cell-type clustering, cell trajectory inference, cell lineages formation, and cells transitioning between states. In this paper, we mainly focus on discovering its applicative power in two fundamental aspects from scRNA-Seq data, i.e., gene imputationterm·mth and cell clusteringterm·key. Gene imputationterm·mth aims to solve the dropoutterm·dom issue which commonly exists in scRNA-Seq data where the expressions of a large number of active genes are marked as zeros 15 – 17 . The excess of zero values often needs to be recovered or handled to avoid the exaggeration of the dropoutterm·dom events in many downstream biological analyses and interpretations. Existing imputation methods 18 , such as MAGIC 4 and SAVER 19 , have an issue in generating biased estimates of gene expression and tend to induce false-positive and biased gene correlations that could possibly eliminate some meaningful biological variations 20 , 21 . On the other hand, many studies, including Seurat 5 and Phenograph 3 , have explored the cell–cell relationships using raw scRNA-seq data, and built cell graphs with reduced data dimensions and detected cell clusters by applying the Louvain modularity optimization. Accurate cell–cell relationships obey the rule that cells are more homogeneous within a cell type and more heterogeneous among different cell types 22 , The scGNNterm·mth model provides a global perspective in exploring cell relationships by integrating cell neighbors on the whole population.

scGNNterm·mth achieves promising performance in gene imputationterm·mth and cell cluster prediction on four scRNA-Seq data sets with gold-standard cell labels 23 – 26 , compared to nine existing imputation and four clustering tools (Supplementary Table 1 ). We believe that the superior performance in gene imputationterm·mth and cell cluster prediction benefits from (i) our integrative autoencoderterm·mth framework, which synergistically determines cell clusters based on a bottom-up integration of detailed pairwise cell–cell relationships and the convergence of predicted clusters, and (ii) the integration of both gene regulatory signals and cell network representations in hidden layers as regularizers of our autoencodersterm·mth. To further demonstrate the power of scGNNterm·mth in complex disease studies, we applied it to an Alzheimer’s diseaseterm·key (AD) data set containing 13,214 single nuclei, which elucidated its application power on cell-type identification and recovering gene expression values 27 . We claim that such a GNN-based framework is powerful and flexible enough to have great potential in integrating scMulti-Omics data.

Results The architecture of scGNNterm·mth comprises stacked autoencodersterm·mth The main architecture of scGNNterm·mth is used to seek effective representations of cells and genes that are useful for performing different tasks in scRNA-Seq data analyses (Fig. 1 and Supplementary Fig. 1 ). It has three comprehensive computational components in an iteration process, including gene regulation integration in a feature autoencoderterm·mth, cell graph representation in a graph autoencoderterm·mth, gene expression updating in a set of parallel cell-type-specific cluster autoencodersterm·mth, as well as the final gene expression recovery in an imputation autoencoderterm·mth (Fig. 1 ). Fig. 1 The architecture of scGNNterm·mth. It takes the gene expression matrix generated from scRNA-Seq as the input. LTMG can translate the input gene expression data into a discretized regulatory signal as the regularizer for the feature autoencoderterm·mth. The feature autoencoderterm·mth learns a dimensional representation of the input as embeddingterm·mth, upon which a cell graph is constructed and pruned. The graph autoencoderterm·mth learns a topological graph embeddingterm·mth of the cell graph, which is used for cell-type clustering. The cells in each cell type have an individual cluster autoencoderterm·mth to reconstruct gene expression values. The framework treats the reconstructed expression as a new input iteratively until converging. Finally, the imputed gene expression values are obtained by the feature autoencoderterm·mth regularized by the cell–cell relationships in the learned cell graph on the original pre-processed raw expression matrix through the imputation autoencoderterm·mth. LTMG is abbreviated for the left-truncated mixed Gaussian model.

The feature autoencoderterm·mth intakes the pre-processed gene expression matrix after the removal of low-quality cells and genes, normalization, and variable gene ranking (Fig. 2a ). First, the LTMG model 14 , 28 is adopted to the top 2,000 variable genes to quantify gene regulatory signals encoded among diverse cell states in scRNA-Seq data (see “Methods” section and Supplementary Fig. 2 ). This model was built based on the kinetic relationships between the transcriptional regulatory inputs and mRNA metabolism and abundance, which can infer the expression of multi-modalities across single cells. The captured signals have a better signal-to-noise ratio to be used as a high-order restraint to regularize the feature autoencoderterm·mth. The aim of this regularization is to treat each gene differently based on their individual regulation status through a penalty in the loss function. The feature autoencoderterm·mth learns a low-dimensional embeddingterm·mth by the gene expression reconstruction together with the regularization. A cell–cell graphterm·mth is generated from the learned embeddingterm·mth via the KNNterm·mth graph, where nodes represent individual cells and the edges represent neighborhood relations among these cells 29 , 30 . Then, the cell graph is pruned from selecting an adaptive number of neighbors for each node on the KNNterm·mth graph by removing the noisy edges 3 . Fig. 2 The architecture of scGNNterm·mth Autoencodersterm·mth. a The feature autoencoderterm·mth takes the expression matrix as the input, regularized by LTMG signals. The dimensions of the encoder and decoder layers are 512 × 128 and 128 × 512, respectively. The feature autoencoderterm·mth is trained by minimizing the difference between the input matrix and the output matrix. b The graph autoencoderterm·mth takes the adjacency matrix of the pruned graph as the input. The encoder consists of two layers of GNNsterm·mth. In each layer, each node of the graph aggregates information from its neighbors. The encoder learns a low-dimensional presentation (i.e., graph embeddingterm·mth) of the pruned cell graph. The decoder reconstructs the adjacency matrix of the graph by dot products of the learned graph embeddingterm·mth followed by a sigmoid activation function. The graph autoencoderterm·mth is trained by minimizing the cross-entropy loss between the input and the reconstructed graph. Cell clusters are obtained by applying k-means and Louvain on the graph embeddingterm·mth. c The cluster autoencoderterm·mth takes a reconstructed expression matrix from the feature autoencoderterm·mth as the input. An individual encoder is built on the cells in each of the identified clusters, and each autoencoderterm·mth is trained individually. The concatenation of the results from all clusters is treated as the reconstructed matrix.

Taking the pruned cell graph as input, the encoder of the graph autoencoderterm·mth uses GNNterm·mth to learn a low-dimensional embeddingterm·mth of each node and then regenerates the whole graph structure through the decoder of the graph autoencoderterm·mth (Fig. 2b ). Based on the topological properties of the cell graph, the graph autoencoderterm·mth abstracts intrinsic high-order cell–cell relationships propagated on the global graph. The low-dimensional graph embeddingterm·mth integrates the essential pairwise cell–cell relationships and the global cell–cell graphterm·mth topology using a graph formulation by regenerating the topological structure of the input cell graph. Then the k-means clustering method is used to cluster cells on the learned graph embeddingterm·mth 31 , where the number of clusters is determined by the Louvain algorithm 31 on the cell graph.

The expression matrix in each cell cluster from the feature autoencoderterm·mth is reconstructed through the cluster autoencoderterm·mth. Using the inferred cell-type information from the graph autoencoderterm·mth, the cluster autoencoderterm·mth treats different cell types specifically and regenerates expression in the same cell cluster (Fig. 2c ). The cluster autoencoderterm·mth helps discover cell-type-specific information for each cell type in its individualized learning. Accompanied by the feature autoencoderterm·mth, the cluster autoencoderterm·mth leverages the inferences between global and cell-type-specific representation learning. Iteratively, the reconstructed matrix is fed back into the feature autoencoderterm·mth. The iteration process stops until it converges with no change in cell clusteringterm·key and this cell clusteringterm·key result is recognized as the final results of cell-type prediction.

After the iteration stops, this imputation autoencoderterm·mth takes the original gene expression matrix as input and is trained with the additional L1 regularizer of the inferred cell–cell relationships. The regularizers (see “Methods” section) are generated based on edges in the learned cell graph in the last iteration and their co-occurrences in the same predicted cell type. Besides, the L1 penalty term is applied to increase the model generalization by squeezing more zeroes into the autoencoderterm·mth model weights. The sparsityterm·dom brought by the L1 termterm·mth benefits the expression imputationterm·mth in dropoutterm·dom effects. Finally, the reconstructed gene expression values are used as the final imputation output.

scGNNterm·mth can effectively impute scRNA-Seq data and accurately predict cell clusters To assess the imputation and cell clusteringterm·key performance of scGNNterm·mth, four scRNA data sets (i.e., Chung 26 , Kolodziejczy 23 , Klein 24 , and Zeisel 25 ) with gold-standard cell-type labels are chosen as the benchmarks (more performance evaluation on other data sets can be found in Supplementary Data 1 – 2 ). We simulated the dropoutterm·dom effects by randomly flipping a number of the non-zero entries to zeros. The synthetic dropoutterm·dom simulation was based on the same leave-one-out strategy used in scVI 32 (Supplementary Fig. 3 ). Median L1 distance, cosine similarity, and root-mean-square-deviation (RMSE) scores between the original data set and the imputed values for these synthetic entries were calculated to compare scGNNterm·mth with MAGIC 4 , SAUCIE 10 , SAVER 19 , scImpute 33 , scVI 32 , DCA 11 , DeepImpute 34 , scIGANs 35 , and netNMF-sc 36 (see “Methods” section). scGNNterm·mth achieves the best results in recovering gene expressions in terms of median L1 distance, and RMSE at the 10 and 30% synthetic dropoutterm·dom rate, respectively. While the cosine similarity score of scGNNterm·mth ranks at the top place for 10% rate and the third place for 30% rate. (Fig. 3a and Supplementary Data 1 ). Furthermore, scGNNterm·mth can recover the underlying gene–gene relationships missed in the raw expression data due to the sparsityterm·dom of scRNA-Seq. For example, two pluripotency epiblast gene pairs, Ccnd3 versus Pou5f1 and Nanog versus Trim28 , are lowly correlated in the original raw data but show strong correlations relations, which are differentiated by time points after scGNNterm·mth imputation and, therefore, perform with a consistency leading to the desired results sought in the original paper 24 (Fig. 3b ). The recovered relations of four more gene pairs are also showcased in Supplementary Figure 4 . scGNNterm·mth amplifies differentially expressed genes (DEGs) signals with a higher fold change than the original, using an imputed matrix to confidently depict the cluster heterogeneity (Fig. 3c ). We also compared the DEG signal changespersName before and after imputation using other imputation tools. As an example, 744 DEGs (logFC > 0.25) identified in Microglia (benchmark cell label) of Zeisel data were compared logFC value change before and after imputation (Supplementary Fig. 5 ). The result turned out that scGNNterm·mth is the only tool that increases all most all DEG signals in Microglia with the strongest Pearson’s correlation coefficient to the original data. Other tools showed weaker coefficients and signals in some of the genes were decreased, indicating imputation bias in these tools. Our results indicate that scGNNterm·mth can accurately restore expression values, capture true gene–gene relations, and increase DEG signals, without inducing additional noises. Fig. 3 Comparison of the imputation performance. a Comparison of the cosine similarity, median L1 distance, and RMSE scores between scGNNterm·mth and other nine imputation tools under 10 and 30% synthetic dropoutterm·dom rate. Darker color indicates better performances. The highest score in each column is highlighted with the yellow box. RMSE scores were scaled by multiplying by 100. b Co-expression patterns can be addressed more explicitly after applying scGNNterm·mth on the Klein data. No clear gene pair relationship of Ccnd3 versus Pou5f1 (upper panel) and Nanog versus Trim28 (lower panel) is observed in the raw data (left) compared to the observation of unambiguous correlations within each cell type after scGNNterm·mth imputation (right). c Comparison of DEG logFC scores using the original expression value ( x axis) and the scGNNterm·mth imputed expression values ( y axis) identified in day 1 cells of the Klein data (up) and microglial cells of the Zeisel data (bottom). The differentiation signals are amplified after imputation.

Besides the artificial dropoutterm·dom benchmarks, we continued to evaluate the clustering performance of scGNNterm·mth and the nine imputation tools on the same two data sets. The predicted cell labels were systematically evaluated using 10 criteria including an adjusted Rand index (ARI) 37 , Silhouette 38 , and eight other criteria (Fig. 4a and Supplementary Data 2 ). By visualizing cell clusteringterm·key results on UMAPs 39 , one can observe more apparent closeness of cells within the same cluster and separation among different clusters when using scGNNterm·mth embeddingsterm·mth compared to the other nine imputation tools (Fig. 4b ). We also observed that compared to the tSNE 40 and PHATE 41 visualization methods, UMAP showed better display results with closer inner-group distance and larger between-group distances (Supplementary Fig. 6 ). The expression patterns show heterogeneity along with embryonic stem cell development. In the case of Klein’s time-series data, scGNNterm·mth recovered a complex structure that was not well represented by the raw data, showing a well-aligned trajectory path of cell development from Day 1 to Day 7 (Fig. 4c ). Moreover, scGNNterm·mth showed significant enhancement in cell clusteringterm·key compared to the existing scRNA-Seq analytical framework (e.g., Seurat using the Louvain method) when using the raw data (Supplementary Fig. 7 ). We hypothesized that the cell–cell graphterm·mth constructed from scGNNterm·mth can reflect cell–cell communications based on ligand–receptor pairs. Using CellChat 42 and curated receptor–ligand pairs, we proved that aggregated interaction probability of cell pairs defined in an scGNNterm·mth cell–cell graphterm·mth is significantly higher than randomly selected cell pairs, which strongly indicates the capability of scGNNterm·mth in capturing the real cell–cell communications and interactions (Supplementary Fig. 8 ). Fig. 4 Cell clusteringterm·key and trajectory evaluations. a Comparison of ARI and Silhouette scores among scGNNterm·mth and nine tools using Klein and Zeisel data sets. b Comparison of UMAP visualizations on the same two data sets, indicating that when scGNNterm·mth embeddingsterm·mth are utilized, cells are more closely grouped within the same cluster but when other tools are used, cells are more separated between clusters. Cells were clustered via the Louvain method and visualized using UMAP. c Pseudotime analysis using the raw expression matrix and scGNNterm·mth imputed matrix of the Klein data set via Monocle. d Justification of using the graph autoencoderterm·mth, the cluster autoencoderterm·mth, and the top 2000 variable genes on the Klein data set in the scGNNterm·mth framework, in terms of ARI. scGNNterm·mth CA- shows the results of the graph autoencoderterm·mth’s ablation, CA- shows the results of the cluster autoencoderterm·mth’s ablation, and AG shows the results after using all genes in the framework.

On top of that, to address the significance of using the graph autoencoderterm·mth and cluster autoencoderterm·mth in scGNNterm·mth, we performed ablation tests to bypass each autoencoderterm·mth and compare the ARI results on the Klein data set (Fig. 4d and Supplementary Fig. 9 ). The results showed that removing either of these two autoencodersterm·mth dramatically decreased the performance of scGNNterm·mth in terms of cell clusteringterm·key accuracy. Another test using all genes rather than the top 2,000 variable genes also showed poor performance in the results and doubled the runtime of scGNNterm·mth, indicating that those low variable genes may reduce the signal-to-noise ratio and negatively affect the accuracy of scGNNterm·mth. The design and comprehensive results of the ablation studies on both clustering and imputation are detailed in Supplementary Methods 1 – 4 , Supplementary Table 2 , and Supplementary Data 3 – 8 . We also extensively studied the parameter selection in Supplementary Data 9 – 12 .

scGNNterm·mth illustrates AD-related neural development and the underlying regulatory mechanism To further demonstrate the applicative power of scGNNterm·mth, we applied it to a scRNA-Seq data set (GEO accession number GSE138852 ) containing 13,214 single nuclei collected from six AD and six control brains 27 . scGNNterm·mth identifies 10 cell clusters, including microglia, neurons, oligodendrocyte progenitor cells (OPCs), astrocytes, and six sub-clusters of oligodendrocytes (Fig. 5a ). Specifically, the proportions of these six oligodendrocyte sub-clusters differ between AD patients (Oligos 2, 3, and 4) and healthy controls (Oligos 1, 5, and 6) (Fig. 5b ). Moreover, the difference between AD and the control in the proportion of astrocyte and OPCs is observed, indicating the change of cell population in AD patients compared to healthy controls (Fig. 5b ). We then combined these six oligodendrocyte sub-clusters into one to discover DEGs. Since scGNNterm·mth can significantly increase true signals in the raw data set, DEG patterns are more explicit (Supplementary Fig. 10 ). Among all DEGs, we confirmed 22 genes as cell-type-specific markers for astrocytes, OPCs, oligodendrocytes, and neurons, in that order 43 (Fig. 5c ). A biological pathway enrichment analysis shows several highly positive enrichments in AD cells compared to control cells among all five cell types. These enrichments include oxidative phosphorylation and pathways associated with AD, Parkinson’s disease, and Huntington disease 44 (Fig. 5d and Supplementary Fig. 11 ). Interestingly, we observed a strong negative enrichment of the MAPK (mitogen-activated protein kinase) signaling pathway in the microglia cells, suggesting a relatively low MAPK regulation in microglia than other cells. Fig. 5 Alzheimer’s diseaseterm·key data set ( GSE138852 ) analysis based on scGNNterm·mth. a Cell clusteringterm·key UMAP. Labeled with scGNNterm·mth clusters (left) and AD/control samples (right). b Comparison of cell proportions in AD/control samples (left) and each cluster (right). c Heatmap of DEGs (logFC > 0.25) in each cluster. Six oligodendrocyte sub-clusters are merged as one to compare with other cell types. Marker genes identified in DEGs are listed on the right. d Selected AD-related enrichment pathways in each cell type in the comparison between AD and control cells. e Underlying TFs are responsible for the cell-type-specific gene regulations identified by IRIS3.

In order to investigate the regulatory mechanisms underlying the AD-related neural development, we applied the imputed matrix of scGNNterm·mth to IRIS3 (an integrated cell-type-specific regulon inference server from single-cell RNA-Seqterm·key) and identified 21 cell-type-specific regulons (CTSR) in five cell types 45 (Fig. 5e and Supplementary Data 13 ; IRIS3 job ID: 20200626160833). Not surprisingly, we identified several AD-related transcription factors (TFs) and target genes that have been reported to be involved in the development of AD. SP2 is a common TF identified in both oligodendrocytes and astrocytes. It has been shown to regulate the ABCA7 gene, which is an IGAP (International Genomics of Alzheimer’s Project) gene that is highly associated with late-onset AD 46 . We also observed an SP2 CTSR in astrocytes that regulate APOE , AQP4 , SLC1A2 , GJA1 , and FGFR3 . All of these five targeted genes are marker genes of astrocytes, which have been reported to be associated with AD 47 , 48 . In addition, the SP3 TF, which can regulate the synaptic function in neurons is identified in all cell clusters, and it is highly activated in AD 49 , 50 . We identified CTSRs regulated by SP3 in OPCs, astrocytes, and neurons suggesting significant SP3-related regulation shifts in these three clusters. We observed 26, 60, and 22 genes that were uniquely regulated in OPCs, astrocytes, and neurons, as well as 60 genes shared among the three clusters (Supplementary Data 14 ). Such findings provide a direction for the discovery of SP3 function in AD studies.

Discussion It is still a fundamental challenge to explore cellular heterogeneity in high-volume, high-sparsity, and noisy scRNA-Seq data, where the high-order topological relationships of the whole-cell graph are still not well explored and formulated. The key innovations of scGNNterm·mth are incorporating global propagated topological features of the cells through GNNsterm·mth, together with integrating gene regulatory signals in an iterative process for scRNA-Seq data analysis. The benefits of GNNterm·mth are its intrinsic learnable properties of propagating and aggregating attributes to capture relationships across the whole cell–cell graphterm·mth. Hence, the learned graph embeddingterm·mth can be treated as the high-order representations of cell–cell relationships in scRNA-Seq data in the context of graph topology. Unlike the previous autoencoderterm·mth applications in scRNA-Seq data analysis, which only captures the top-down distributions of the overall cells, scGNNterm·mth can effectively aggregate detailed relationships between similar cells using a bottom-up approach. We also observed that the imputation of scGNNterm·mth can decrease batch effectsterm·dom introduced by different sequencing technologies (Supplementary Fig. 12 ), which makes scGNNterm·mth a good choice for data imputation prior to multiple scRNA-Seq data integration 51 . Furthermore, scGNNterm·mth integrates gene regulatory signals efficiently by representing them discretely in LTMG in the feature autoencoderterm·mth regularization. These gene regulatory signals can help identify biologically meaningful gene–gene relationships as they apply to our framework and eventually, they are proven capable of enhancing performance. Technically, scGNNterm·mth adopts multi-modal autoencodersterm·mth in an iterative manner to recover gene expression values and cell-type prediction simultaneously. Notably, scGNNterm·mth is a hypothesis-free deep learning framework on a data-driven cell graph model, and it is flexible to incorporate different statistical models (e.g., LTMG) to analyze complex scRNA-Seq data sets.

Some limitations can still be found in scGNNterm·mth. (i) It is prone to achieve better results with large data sets, compared to relatively small data sets (e.g., <1000 cells), as it is designed to learn better representations with many cells from scRNA-Seq data, as shown in the benchmark results, and (ii) Compared with statistical model-based methods, the iterative autoencoderterm·mth framework needs more computational resources, which is more time-consuming (Supplementary Data 15 ). In the future, we will investigate creating a more efficient scGNNterm·mth model with a lighter and more compressed architecture.

In the future, we will continue to enhance scGNNterm·mth by implementing heterogeneous graphs to support the integration of single-cell multi-omics data (e.g., the intra-modality of Smart-Seq2 and Droplet scRNA-Seq data; and the inter-modality integration of scRNA-Seq and scATAC-Seq data). We will also incorporate attention mechanisms and graph transformer models 52 to make the analyses more explainable. Specifically, by allowing the integration of scRNA-Seq and scATAC-Seq data, scGNNterm·mth has the potential to elucidate cell-type-specific gene regulatory mechanisms 53 . On the other hand, T cell receptor repertoires are considered as unique identifiers of T cell ancestries that can improve both the accuracy and robustness of predictions regarding cell–cell interactions 54 . scGNNterm·mth can also facilitate batch effectsterm·dom and build connections across diverse sequencing technologies, experiments, and modalities. Moreover, scGNNterm·mth can be applied to analyze spatial transcription data sets regarding spatial coordinates as additional regularizers to infer the cell neighborhood representation and better prune the cell graph. We plan to develop a more user-friendly software system from our scGNNterm·mth model, together with modularized analytical functions in support of standardizing the data format, quality control, data integration, multi-functional scMulti-seq analyses, performance evaluations, and interactive visualizations.

Methods Data set preprocessing scGNNterm·mth takes the scRNA-Seq gene expression profile as the input. Data filtering and quality control are the first steps of data preprocessing. Due to the high dropoutterm·dom rate of scRNA-seq expression data, only genes expressed as non-zero in more than 1% of cells, and cells expressed as non-zero in more than 1% of genes are kept. Then, genes are ranked by standard deviation, i.e., the top 2000 genes in variances are used for the study. All the data are log-transformed.

Left-truncated mixed Gaussian (LTMG) modeling A mixed Gaussian model with left truncation assumption is used to explore the regulatory signals from gene expression 14 . The normalized expression values of gene X over N cells are denoted as X = { x 1 ,… x N }, where x j ∈ X is assumed to follow a mixture of k Gaussian distributions, corresponding to k possible gene regulatory signals (TRSs). The density function of X is: 1 p X ; Θ = ∏ j = 1 N p ( x j ; Θ ) = ∏ j = 1 N ∑ i = 1 k α i p x j ; θ i = ∏ j = 1 N ∑ i = 1 k α i 1 2 π σ i e − x j − μ i 2 2 σ i 2 = L Θ ; X where α i is the mixing weight, μ i and σ i are the mean and standard deviation of the i th Gaussian distribution, which can be estimated by: Θ * = arg max L ( Θ ; X ) Θ to model the errors at zero and the low expression values. With the left truncation assumption, the gene expression profile is split into M , which is a truly measured expression of values, and N − M representing left-censored gene expressions for N conditions. The parameter Θ maximizes the likelihood function and can be estimated by an expectation-maximization algorithm. The number of Gaussian components is selected by the Bayesian Information Criterion; then, the original gene expression values are labeled to the most likely distribution under each cell. In detail, the probability that x j belongs to distribution i is formulated by: 2 p x j ∈ TRS i ∣ K , Θ * ∝ α i 2 π σ j 2 e − x j − μ i 2 2 σ i 2 where x j is labeled by TRS i if p x j ∈ TRS i ∣ K , Θ * = max i = 1 , … , K ( p x j ∈ TRS i ∣ K , Θ * ) . Thus, the discrete values (1,2, …, K ) for each gene are generated.

Feature autoencoderterm·mth The feature autoencoderterm·mth is proposed to learn the representative embeddingterm·mth of the scRNA expression through stacked two layers of dense networks in both the encoder and decoder. The encoder constructs the low-dimensional embeddingterm·mth of X ′ from the input gene expression X , and the encoder reconstructs the expression X ^ from the embeddingterm·mth; thus, X , X ^ ∈ R N × M and X ′ ∈ R N × M ′ , where M is the number of input genes, M ′ is the dimension of the learned embeddingterm·mth, and M ′ < M . The objective of training the feature autoencoderterm·mth is to achieve a maximum similarity between the original and reconstructed through minimizing the loss function, in which ∑ X − X ^ 2 is the main term serving as the mean squared error (MSE) between the original and the reconstructed expressions.

Regularization Regularization is adopted to integrate gene regulation information during the feature autoencoderterm·mth training process. The aim of this regularization is to treat each gene differently based on their individual gene regulation role through penalizing it in the loss function. The MSE is defined as: 3 α ∑ X − X ^ 2 ∘ TRS where TRS ∈ R N × M ; α is a parameter used to control the strength of gene regulation regularization; α ∈ [0,1]. ∘ denotes element-wise multiplication. Thus, the loss function of the feature autoencoderterm·mth is shown as Eq.( 4 ). 4 Loss = ( 1 − α ) ∑ X − X ^ 2 + α ∑ X − X ^ 2 ∘ TRS In the encoder, the output dimensions of the first and second layers are set as 512 and 128, respectively. Each layer is followed by the ReLU activation function. In the decoder, the output dimensions of the first and second layers are 128 and 512, respectively. Each layer is followed by a sigmoid activation function. The learning rate is set as 0.001. The cluster autoencoderterm·mth has the same architecture as the feature autoencoderterm·mth, but without gene regulation regularization in the loss function.

Cell graph and pruning The cell graph formulates the cell–cell relationships using embeddingterm·mth learned from the feature autoencoderterm·mth. As done in the previous works 4 , 55 , the cell graph is built from a KNNterm·mth graph, where nodes are individual single cells, and the edges are relationships between cells. K is the predefined parameter used to control the scale of the captured interaction between cells. Each node finds its neighbors within the K shortest distances and creates edges between them and itself. Euclidian distance is calculated as the weights of the edges on the learned embeddingterm·mth vectors. The pruning process selects an adaptive number of neighbors for each node on the original KNNterm·mth graph and keeps a more biologically meaningful cell graph. Here, Isolation Forest is applied to prune the graph to detect the outliner in the K -neighbors of each node 56 . Isolation Forest builds individual random forest to check distances from the node to all K -neighbors and only disconnects the outliners.

Graph autoencoderterm·mth The graph autoencoderterm·mth learns to embed and represent the topological information from the pruned cell graph. For the input pruned cell graph, G = ( V , E ) with N = | V | nodes denoting the cells and E representing the edges. A is its adjacency matrix and D is its degree matrix. The node feature matrix of the graph autoencoderterm·mth is the learned embeddingterm·mth X ′ from the feature autoencoderterm·mth.

The graph convolution network (GCN) is defined as GCN X ′ , A = ReLU ( A ~ X ′ W ) , and W is a weight matrix learned from the training. A ~ = D − 1 / 2 A D − 1 / 2 is the symmetrically normalized adjacency matrix and activation function ReLU(∙) = max(0, ∙). The encoder of the graph autoencoderterm·mth is composed of two layers of GCN, and Z is the graph embeddingterm·mth learned through the encoder in Eq.( 5 ). W 1 and W 2 are learned weight matrices in the first and second layers, and the output dimensions of the first and second layers are set at 32 and 16, respectively. The learning rate is set at 0.001. 5 Z = ReLU ( A ~ R e L U A ~ X ′ W 1 W 2 ) The decoder of the graph autoencoderterm·mth is defined as an inner product between the embeddingterm·mth: 6 A ^ = sigmoid ( Z Z T ) where A ^ is the reconstructed adjacency matrix of A . sigmoid(∙) = 1/(1 + e −∙ ) is the sigmoid activation function.

The goal of learning the graph autoencoderterm·mth is to minimize the cross-entropy L between the input adjacency matrix A and the reconstructed matrix A ^ : 7 L A , A ^ = − 1 N × N ∑ i = 1 N ∑ j = 1 N ( a i j * log a ^ i j + 1 − a i j * log ( 1 − a ^ i j ) ) where a ij and a ^ i j are the elements of the adjacency matrix A and A ^ in the i th row and the j th column. As there are N nodes as the cell number in the graph, N × N is the total number of elements in the adjacency matrix.

Iterative process The iterative process aims to build the single-cell graph iteratively until converging. The iterative process of the cell graph can be defined as: 8 A ~ = λ L 0 + 1 − λ A i j ∑ j A i j where L 0 is the normalized adjacency matrix of the initial pruned graph, and L 0 = D 0 − 1 / 2 A 0 D 0 − 1 / 2 , where D 0 is the degree matrix. λ is the parameter to control the converging speed, λ ∈ [0,1]. Each time in iteration t , two criteria are checked to determine whether to stop the iteration: (1) that is, to determine whether the adjacency matrix converges, i.e., A ~ t − A ~ t − 1 < γ 1 A ~ 0 , or (2) whether the inferred cell types are similar enough, i.e., ARI < γ 2 . ARI is the similarity measurement, which is detailed in the next section. In our setting, λ = 0.5 and γ 1 , γ 2 = 0.99. The cell-type clustering results obtained in the last iteration are chosen as the final cell-type results.

Imputation autoencoderterm·mth After the iterative process stops, the imputation autoencoderterm·mth imputes and denoises the raw expression matrix within the inferred cell–cell relationship. The imputation autoencoderterm·mth shares the same architecture as the feature autoencoderterm·mth, but it also uses three additional regularizers from the cell graph in Eq. ( 9 ), cell types in Eq. ( 10 ), and the L1 regularizer in Eq. ( 11 ): 9 γ 1 ∑ ( A ⋅ ( X − X ^ ) 2 ) where A ∈ R N × N is the adjacency matrix from the pruned cell graph in the last iteration. ∙ denotes dot product. Cells within an edge in the pruned graph will be penalized in the training: 10 γ 2 ∑ ( B ⋅ ( X − X ^ ) 2 ) B i j = 1 where i and j in same cell type 0 else where B ∈ R N × N is the relationship matrix between cells, and two cells in the same cell type have a B ij value of 1. Cells within the same inferred cell type will be penalized in the training. γ 1 , γ 2 are the intensities of the regularizers and γ 1 , γ 2 ∈ [0,1]. The L1 regularizer is defined as 11 β ∑ w which brings sparsityterm·dom and increases the generalization performance of the autoencoderterm·mth by reducing the number of non-zero w terms in ∑| w |, where β is a hyper-parameter controlling the intensity of the L1 termterm·mth ( β ∈ [0,1]). Therefore, the loss function of the imputation autoencoderterm·mth is 12 Loss = 1 − α ∑ X − X ^ 2 + α ∑ X − X ^ 2 ∘ TRS + β ∑ w + γ 1 ∑ A ⋅ X − X ^ 2 + γ 2 ∑ ( B ⋅ ( X − X ^ ) 2 ) Benchmark evaluation compared to existing tools Imputation evaluation For benchmarking imputation performance, we performed synthetic dropoutterm·dom simulation to randomly flip 10% of the non-zero entries to zeros. These synthetic dropoutsterm·dom still follow the zero-inflated negative binomial (ZINB) distribution with details shown in Supplementary Method 5 and Data 16 . We evaluated median L1 distance, cosine similarity, and root-mean-squared error (RMSE) between the original data set and the imputed values for these corrupted entries. For all the flipped entries, x is the row vector of the original expression, and y is its corresponding row vector of the imputed expression. The L1 distance is the absolute deviation between the value of the original and imputed expression. A lower L1 distance means a higher similarity. 13 L 1 distance = x − y , L 1 distance ∈ [ 0 , + ∞ ) The cosine similarity computes the dot products between original and imputed expression. 14 Cosine similarity x , y = x y T x y , Cosine similarity ∈ [ 0 , 1 ] The RMSE computes the squared root of the quadratic mean of differences between original and imputed expression. 15 RMSE x , y = ∑ i = 1 N x i − y i 2 N , RMSE ∈ 0 , + ∞ The process is repeated three times, and the mean and standard deviation were selected as a comparison. The scores are compared between scGNNterm·mth and nine imputation tools (i.e., MAGIC 4 , SAUCIE 10 , SAVER 19 , scImpute 33 , scVI 32 , DCA 11 , DeepImpute 34 , scIGANs 35 , and netNMF-sc 36 ), using the default parameters.

Clustering evaluation We compared the cell clusteringterm·key results of scGNNterm·mth, the same nine imputation tools, and four scRNA-Seq analytical frameworks, in terms of ten clustering evaluation scores. Noted that, we considered the default cell clusteringterm·key method (i.e., Louvain method 31 in Seurat 5 , Ward.D2 57 method in CIDR 58 , Louvain method in Monocle 59 , and k -means 60 method in RaceID 61 ) in each of the analytical frameworks to compare the cell clusteringterm·key performance with scGNNterm·mth. The default parameters are applied in all test tools. ARI 37 is used to compute similarities by considering all pairs of the samples that are assigned in clusters in the current and previous clustering adjusted by random permutation: 16 ARI = RI − E RI max ( RI ) − E [ RI ] where the unadjusted rand index (RI) is defined as: 17 RI = a + b C n 2 where a is the number of pairs correctly labeled in the same sets, and b is the number of pairs correctly labeled as not in the same data set. C n 2 is the total number of possible pairs. E[RI] is the expected RI of random labeling.

Different from ARI which requires known ground truth labels, the Silhouette coefficient score 38 defines how similar an object is to its own cluster compared to other clusters. It is defined as: 18 Silhouette = b − a max a , b where a is the mean distance between a sample and all other points in the same class, b is the mean distance between a sample and all other points in the next nearest cluster. Silhouette ∈ [−1,1], where 1 indicates the best clustering results and −1 indicates the worst. We calculated the average Silhouette score of all cells in each data set to compare the cell clusteringterm·key results. More quantitative measurements are also used in Supplementary Method 4 .

Statistical validation of cell–cell graphterm·mth topology based on LRPs We used CellChat 42 to predict potential interaction probability scores (ranging from 0 to 1; a higher score indicates the two cells are more likely to interact with each other) of ligand–receptor pairs (LRP) between any two cells. We built a fully connected cell–cell background graph (using all the cells) based on Pearson’s correlation of the raw expression matrix and compared it with the cell–cell graphterm·mth generated from scGNNterm·mth. CellChat calculates an aggregated interaction probability for each linked cell pair based on the expression level of LRPs. For all linked cell pairs in the background graph and scGNNterm·mth cell–cell graphterm·mth, we performed a Wilcoxon test to evaluate the statistical significance between the corresponding aggregated interaction probability. Five scRNA-Seq data sets (i.e., Klein, Zeisel, Kolo, Chung, and AD) were used in this analysis.

Case study of the AD database We applied scGNNterm·mth on public Alzheimer’s diseaseterm·key (AD) scRNA-Seq data with 13,214 cells 27 . The resolution of scGNNterm·mth was set to 1.0, KI was set to 20, and the remaining parameters were kept as default. The AD patient and control labels were provided by the original paper and used to color the cells on the same UMAP coordinates generated from scGNNterm·mth. We simply combined cells in six oligodendrocyte subpopulations into one cluster, referred to as merged oligo. The DEGs were identified in each cell cluster via the Wilcoxon rank-sum test implemented in the Seurat package along with adjusted p -values using the Benjamini-Hochberg procedure with a nominal level of 0.05. DEGs with logFC > 0.25 or <−0.25 were finally selected. We further identified the DEGs between AD and control cells in each cluster using the same strategy and applied GSEA for pathway enrichment analysis 62 . The imputed matrix, which resulted from scGNNterm·mth was then sent to IRIS3 for CTSR prediction, using the predicted cell clusteringterm·key labels with merged oligodendrocytes 45 . The default parameters were served in regulatory analysis in IRIS3.

Software implementation Tools and packages used in this paper include: Python version 3.7.6, numpy version 1.18.1, torch version 1.4.0, networkx version 2.4, pandas version 0.25.3, rpy2 version 3.2.4, matplotlib version 3.1.2, seaborn version 0.9.0, umap-learn version 0.3.10, munkres version 1.1.2, R version 3.6.1, and igraph version 1.2.5. The IRIS3 website is at https://bmbl.bmi.osumc.edu/iris3/index.php .

Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Supplementary information Supplementary Information Description of Additional Supplementary Files Supplementary Data 1–16 Reporting Summary