The intuition behind the netSmooth algorithm is that gene networks encoding co-expression patterns can be used to smooth scRNA-seq data, pushing its coexpression patterns in a biologically meaningful direction. We demonstrate this using protein-protein interaction networks, which are predictive of coexpression ⁹ . We produced a PPI graph of high-confidence interactions based on the PPI database STRING ¹⁰ .

There are 2 inputs to the method: (1) a gene expression matrix, N genes by M cells, and (2) a graph where genes are nodes, and edges indicate genes which are expected to be co-expressed. The edges may be weighed, indicating the strength or direction of a relationship; an edge weight of 2 indicates stronger expected co-expression than an edge weight of 1, and an edge weight of −1 indicates negative expected co-expression, such as one gene being a repressor for another. The expression profile of each cell is then projected onto the graph, and a diffusion process is used to smooth the expression values, within each sample, of adjacent genes in the graph ( Figure 1). In this way, post-smoothing values of genes represent an estimate of activity levels based on reads aligned to that gene, as well as those aligned to its neighbors in the graph. Thus, a gene with a low read count (possible technical drop-out), whose neighbors in the graph are highly expressed, will get a higher value post smoothing. The rate at which expression values of genes diffuse to their neighbors is degree-normalized, so that genes with many edges will affect their neighbors less than genes with more specific interactions. The diffusion is done using a "random walks with restarts" (RWR) process ¹³ , where a conceptual random walker starts in some node in the graph, and at each iteration moves to a neighboring node with a probability determined by the edge weight between the nodes, or, with some probability, restarts the walk from the original node. The network-smoothed value is the stationary distribution of this process. The RWR process has one free parameter, the restart rate. A low value for the restart rate allows diffusion to reach further in the graph; a high restart rate will lead to more local diffusions. For more details see the Methods section.

We first assess netSmooth on a dataset of 1645 mouse hematopoietic stem/progenitor cells (HSPCs) assayed using flow cytometry as well as scRNA-seq ¹⁵ . The cells are FACS-sorted into 12 common HSPC phenotypes. This presents an atlas of the hematopoiesis process at a single cell resolution, showing the differentiation paths taken by E-SLAM HSCs as they differentiate to E, GM, and L progenitors. The authors of this study demonstrate that upon clustering the data, some clusters corresponds to cell types. However, the clusters are not noise free and do not fully recapitulate cell type identity. We obtained clusterings of the cells from the normalized counts, as well as after application of netSmooth, MAGIC ⁷ , and scImpute ⁶ , using a robust clustering procedure based on the clusterExperiment R package ¹⁶ (See Methods). After clustering, we used the edgeR-QLF test ¹⁷ to identify genes that are differentially expressed in any of the discovered clusters. Figure 2a,b shows that after network-smoothing, we are able to identify clusters with a more pronounced differential expression profile. Further, many more of the genes identified as differentially expressed between the clusters (without smoothing) seem to have low and uninformative expression values overall. MAGIC and scImpute also improve this pattern ( Figure 2c,d). MAGIC seems to do the strongest transformation to the data, as seen in lower dimension embeddings ( Figure S2, Figure S3).

As this dataset has cells with labels independent of the RNAseq (FACS-sorted phenotypes), it presents us with an opportunity to compare the gene expression levels (as measured by RNAseq), to a meaningful phenotypic variable, i.e. the cell type. The cell type discrimination of a clustering result is compared using a cluster purity metric and and the adjusted mutual information (AMI). The cluster purity measures how cell-type specific clusters are by comparing homogeneity of the external labels (FACS-defined cell types), within clusters provided by scRNA-seq data. AMI is a chance-adjusted information-theoretic measure of agreement between two labellings. This method accounts for artificially high mutual information between external labels and clusters when there are high number of clusters (See Methods for details on metrics). We also measured number of cells in robust clusters as quantitative metric. The robust clustering procedure allows cells to be omitted (not be assigned to a cluster) if they cannot be placed in a cluster across multiple clustering methods and/or parameters (See Methods). Only MAGIC is able to increase the proportion of cells in this dataset which fall into robust clusters ( Figure 3a), but only netSmooth leads to more biologically meaningful clusters, in terms of purity and AMI ( Figures 3b,c), demonstrating that netSmooth can assist in cell type identification, and outperformed both MAGIC and scImpute in this task. The higher clusterability following application of MAGIC than netSmooth, might indicate that MAGIC was overzealous in its transformation, squeezing more cells into the same space. This might lead to more robust clusters, but less reliable cell type identification.

Next, we test netSmooth on 269 isolated cells from mouse embryos at different stages of pre-implantation development between oocyte and blastocyst, as well as 5 liver cells and 10 fibroblast cells ¹⁸ . The authors of this study demonstrated that lower dimension embeddings capture much of the developmental trajectory ( Figure 4a, Figure S4a, Figure S5a). We then applied netSmooth, MAGIC, and scImpute. Figure 4b shows the principal component analysis of netSmooth-processed data, and Figures 4c and 4d show the PCA plot following application of MAGIC and scImpute, respectively. netSmooth and scImpute preserve most of the variance structure of the data, while MAGIC seems to push the data onto a completely different manifold ( Figure 4, Figure S5). We used the robust clustering procedure to obtain clusters, and computed the cluster purity and AMI metrics. netSmooth enabled the clustering procedure to place more of the samples into robust clusters ( Figure 5a), and as in the hematopoiesis case, netSmooth is able to assist in identifying the developmental stage or tissue that cells belong to better than the other methods, as evidenced by the higher cluster purities ( Figure 5b) and AMI ( Figure 5c). Although MAGIC and scImpute reduce the 0-count genes further than netSmooth ( Figure S1), they do not add as much clarity to the developmental stage signal inherent in the data. This shows that imputing missing counts based on data from the same experiment is not as powerful as including priors in the quasi-imputation process netSmooth does.

Finally, we demonstrate applicability of netSmooth to cancer research. Patel et al. generated scRNA-seq data of 800 cells from 5 glioblastoma tumors and 2 cell lines ¹⁹ . Lower dimension embedding plots show that cells from different tumors or cell lines generally group together, but some are not wholly distinguishable from other tumors ( Figure 6a, Figure S4a, Figure S5a). Further, the two cell lines group closer to each other than the other patient samples. After applying netSmooth to the data, tumors become easier to distinguish in a lower dimensional embedding ( Figure 6b), indicating that netSmooth improves assignment of each cell to its tumor, cell line, or clone of origin. Again, scImpute also leads to similar reduced dimension embedding ( Figure 6d), while MAGIC distorted the data more than the other methods ( Figure 6c). We used the robust clustering procedure before and after netSmooth, MAGIC, and scImpute. Only MAGIC increase the clusterabitliy of the data ( Figure 7a), but netSmooth leads to the most pure clusters, in terms of tumor or cell line of origin ( Figure 7b, Figure 7c).

Tumor or cell line of origin is an imperfect proxy for phenotypical variation in cancer cells, because some cells cluster by cell type rather than tumor of origin, demonstrating the heterogeneity in these glioblastoma tumors and similarities across origins ¹⁹ . Nevertheless, we chose to compute cluster purity based on the cell origin rather than other labels which might be assigned to them, as it is the only ground truth variable that is independent of the RNAseq experiment. Further, cells do group by origin ( Figure 6, Figure S6), and identification of origin is an interesting question in its own right in the field of cancer genomics, particularly for heterogeneous tumors such as these.

Next, we set out to ensure that the results are not an artifact of the network structure, i.e. that the actual links between genes that we used in the network are important. We expect netSmooth not to perform well when using networks with similar characteristics, but where edges do not represent real interactions. To that effect, we constructed 20 random networks by keeping the same graph structure of the real PPI graph, but shuffling the gene names.

Thus, these random networks share all the characteristics of the real network (degree distribution, community structure), except for the true identity of the nodes. We then used those networks as inputs to netSmooth and ran the benchmarks as before on the hematopoiesis dataset. Using random networks as an input to netSmooth gives cluster purities distributed around a mode given by the cluster purities of the raw data, while the cluster purities given from using the real PPI network lie at the extreme edge of the distribution ( Figure 8a). Further, most random networks result in fewer samples belonging to robust clusters ( Figure 8b). These results demonstrate that it is indeed the information contained in the PPI graph enables netSmooth to transform the gene expression matrix in a more biologically coherent direction, and that the transformation we see can not be explained simply by the network structure.

In addition to using an unweighed (where all edge weights are 1), undirected (where all edge weights are positive) network from string-db, we constructed other gene networks and used them as inputs to netSmooth. We created a directed gene network from only those edges in string-db which are marked as activating or inhibiting ⁱ . We set the edge weights of the activating interactions to +1, and −1 for the inhibiting interactions, allowing gene expression values to be adjusted downwards for genes whose known antagonists are highly expressed. After smoothing, we set all negative smoothed expression values to 0. We also constructed a gene network from string-db using only genes that are known to demonstrate cell-type specific expression. In order to obtain a list of genes with such cell-type specific expression patterns from the Expression Atlas ²⁰ , we used only the genes which show a cell-type specific expression with a mean TPM of at least 1 in some cell type, and used the subset of string-db network containing those genes as an input to netSmooth. Both of those modified graphs perform similarly to the undirected graph from string-db ( Figure 9, Figure S8a, S8b), demonstrating that netSmooth is able to use priors from different types of experiments in order to improve clustering of scRNA-seq.

We also considered other sources for the gene network. We constructed a gene network from HumanNet ²¹ , a functional gene network where edges denote interactions between two genes. We constructed a smoothing graph by taking all edges from HumanNet, and producing a graph where all edge weights are set to 1. We then used this graph as an input to netSmooth on the glioblastoma dataset. It performs similarly to the network from string-db ( Figure 10, Figure S8c), demonstrating that other sources for gene interactions may also be used by netSmooth to improve clustering results of scRNA-seq.

The netSmooth algorithm, given a gene network, has one free parameter - the restart rate of the random walker, (1 − α). Alternatively, α is the complement of the restart rate. An α = 0 indicates a perfect restart rate and consequently no smoothing; an α = 1 corresponds to a random walk without restarts. Intermediate values for α result in increasing levels of smoothing; the value of α determines how far random walks will go on the graph before restarting, or how far along the network a gene’s influence is allowed to reach (See Methods). It is tempting to optimize α with respect to the variable the experiment sets out to measure, e.g. cluster purity. For instance, in the embryonic development dataset, we would choose α = 0.7 as the value that produces the highest cluster purity ( Figure 11b). However, in many experiments the identity of the samples is not known a-priori. Therefore, we propose a data driven workflow to pick a sensible value for α.

One such data-driven statistic is the proportion of samples assigned to robust clusters; following application of netSmooth, the robust clustering procedure is able to assign more samples to statistically robust clusters. For all three datasets, picking the α that gives the highest proportion of cells in robust clusters, also gives the clusters with the highest purity index ( Figure 12). Importantly, this metric is entirely data-driven and does not require external labels, making it feasible for any scRNA-seq study. The results in the previous sections all use the value of α picked to optimize proportion in robust clusters.

The netSmooth algorithm takes a graph G = { V, E} where V = { gene _i } is the set of genes, and E = {( i → j)} is the set of edges between genes. The edge weights are degree-normalized, so that each gene’s outgoing edges’ weights sum to 1. We then define a process of random walk with restarts as in 13, on the PPI graph, where a conceptual random walker starts on a node in the graph (a gene/protein) and at each step walks to an adjacent node with the probability determined by the α times the edge weight. Further, at each step, there is a probability of (1 − α) that the walker restarts to its original node.

Mathematically, given a graph defined by an adjacency matrix A _{[ MxM]}, where A _{i j} is the edge weight between gene i and gene j (and 0 for unconnected genes), and a vector f _{[ Mx1]}, where f i t is the probability that the walker is at node i at step t, the process is defined by

f t + 1 = α A f t + ( 1 − α ) f 0 .

This process is convergent, and the stationary distribution is given by

f ∞ = ( 1 − α ) ( I − α A ) − 1 f 0 .

Hence, the random walk with restarts process is a diffusion process defined on the PPI graph, or through the diffusion kernel (smoothing kernel)

K A α = ( 1 − α ) ( I − α A ) − 1

where (1 − α) is the restart probability, and A is the (column normalized) adjacency matrix of the PPI graph. Consequently, we define the network-smoothed expression profile

E s m = K A α E ,

where E _{[ MxN]} is the normalized count values of the M genes in the N cells.

Clustering analysis features prominently in scRNA-seq analyses; whether recapitulating known results or discovering new cell types, clustering cells by their gene expression profiles is commonly used to identify distinct populations. While some approaches directly take into account the zero-inflation of scRNA-seq data ⁵ , other studies use traditional methods ¹⁸ . There is no standard method for clustering single cell RNAseq data, as different studies produce data with different topologies, which respond differently to the various clustering algorithms.

In order to avoid optimizing different clustering routines for the different datasets we benchmark on, we have implemented a robust clustering routine based on clusterExperiment ^{ii 16} , a framework for robust clustering based on consensus clustering of clustering assignments obtained from different clustering algorithms, different parameters for these algorithms, and different views of the data. The different views are different reduced dimensionality projections of the data based on different techniques. Thus, no single clustering result will dominate the data, and only cluster structures which are robust to different analyses will prevail. The procedure we implemented using the framework is as follows:

Using only the 500 most variable genes insures the biological variation will dominate the technical variation, and enhances the reproducibility of t-SNE ²² .

Importantly, samples that at step (4) don’t have a high enough affinity to any emerging cluster, will not be assigned to any cluster. The clustering is performed using the clusterExperiment::clusterSingle and clusterExperiment::clusterMany functions, the consensus clustering is obtained using the clusterExperiment::combineMany function, and the cluster merging (steps 5 and 6) using the clusterExperiment::makeDendrogram and clusterExperiment::mergeClusters functions. For more details, see 16.

In step (1) above, we cluster cells in a lower dimension embedding using either PCA ²³ or t-SNE ²⁴ , in a dataset-dependent manner. Different single cell datasets respond better to different dimensionality reduction techniques which are better able to tease out the biological cluster structure of the data. In order to pick the right technique algorithmically, we compute the entropy in a 2D embedding. We obtained 2D embeddings from the 500 most variable genes using either PCA or t-SNE, binned them in a 20x20 grid, and computed the entropy using the discretize and entropy functions in the entropy R package ^{iii 25} . The entropy in the 2D embedding is a measure for the information captured by it. For the clustering procedure, we pick the embedding with the highest information content. For the hematopoiesis and glioblastoma datasets, this is t-SNE, while for the embryonic development dataset it is PCA ( Table 2). This method may be used to pick any dimensionality reduction technique other than the ones mentioned here, which might be more suitable for other analyses.

The cluster purity metric displayed above refers to the proportion of the samples in a cluster which are of the dominant cell type in that cluster. The purity for cluster i is given by

P u r i t y i = ∑ j ∈ C i { 1 , if l a b e l j = dom i 0 , otherwise n i

where C _i = { j|cell _j ∈ cluster _i }, label _j is the cell type of cell _j , n _i = | C _i | is the number of cells in cluster i, and

dom i = arg ⁡ max ⁡ l ∑ j ∈ C i { 1 , if l a b e l j = l 0 , otherwise

is the dominant cell type in cluster C _i .

In addition to the cluster purity metric, we computed the Adjusted Mutual Information (AMI) ²⁶ , an information theoretic measure of clustering accuracy which accounts for true positives (two cells of the same type in the same cluster) being caused by chance. The AMI between a clustering C and the true labels L is given by

A M I ( L , C ) = M I ( L , C ) − E [ M I ( L , C ) ] max ⁡ ( H ( L ) , H ( C ) ) − E [ M I ( L , C ) ] ,

where MI( a, b) is the mutual information between labellings a and b, H( a) is entropy of clustering a, and E[ ·] denotes the expectation.

We do not compare the clusterings using the Rand index, as that measure penalizes for so-called false negatives (two cells of the same cell type but in different clusters), which is undesirable as cells from the same cell type might be rightly split into several clusters when a novel cell type is identified.

The PPI graph from which the diffusion kernel was derived was constructed using data from string-db ¹⁰ . For each pair of proteins, string-db provides a combined interaction score, which is a score indicating how confident we can be in the interaction between the proteins, given the different kinds of evidence string-db collates. We subset the links to only those above the 90th percentile of combined interaction scores, only keeping the 10% most confident interactions. For mouse that is 1,020,816 interactions among 17013 genes. For human, 852,722 interactions among 17467 genes.

For all the results presented in this paper, scImpute was run using the default parameters ( drop_thre = 0.5). For MAGIC, we used values for the diffusion time parameter ( T = { 1, 2, 4, 8, 16}). Unlike netSmooth, for MAGIC the proportion of samples in robust clusters and the cluster purities were anti-correlated; thus we picked the one that gave the best cluster purities as the best MAGIC parameter. The chosen T values are given in Table 3. We used MAGIC version 0.1 ^iv and scImpute version 0.0.2 ^v .