1 Introduction

netSmooth implements a network-smoothing framework to smooth single-cell gene expression data as well as other omics datasets. The algorithm is a graph based diffusion process on networks. The intuition behind the algorithm is that gene networks encoding coexpression patterns may be used to smooth scRNA-seq expression data, since the gene expression values of connected nodes in the network will be predictive of each other. Protein-protein interaction (PPI) networks and coexpression networks are among the networks that could be used for such procedure.

More precisely, netSmooth works as follows. First, the gene expression values or other quantitative values per gene from each sample is projected on to the provided network. Then, the diffusion process is used to smooth the expression values of adjacent genes in the graph, so that a genes expression value represent an estimate of expression levels based the gene it self, as well as the expression values of the neighbors in the graph. 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 implementation has one free parameter, alpha, which controls if the diffusion will be local or will reach further in the graph. Higher the value, the further the diffusion will reach. The netSmooth package implements strategies to optimize the value of alpha.

Network-smoothing concept

Figure 1: Network-smoothing concept

In summary, netSmooth enables users to smooth quantitative values associated with genes using a gene interaction network such as a protein-protein interaction network. The following sections of this vignette demonstrate functionality of netSmooth package.

2 Smoothing single-cell gene expression data with netSmooth() function

The workhorse of the netSmooth package is the netSmooth() function. This function takes at least two arguments, a network and genes-by-samples matrix as input, and performs smoothing on genes-by-samples matrix. The network should be organized as an adjacency matrix and its row and column names should match the row names of genes-by-samples matrix.

We will demonstrate the usage of the netSmooth() function using a subset of human PPI and a subset of single-cell RNA-seq data from GSE44183-GPL11154. We will first load the example datasets that are available through netSmooth package.

data(smallPPI)
data(smallscRNAseq)

We can now smooth the gene expression network now with netSmooth() function. We will use alpha=0.5.

smallscRNAseq.sm.se <- netSmooth(smallscRNAseq, smallPPI, alpha=0.5)
## Using given alpha: 0.5
smallscRNAseq.sm.sce <- SingleCellExperiment(
    assays=list(counts=assay(smallscRNAseq.sm.se)),
    colData=colData(smallscRNAseq.sm.se)
)

Now, we can look at the smoothed and raw expression values using a heatmap.

anno.df <- data.frame(cell.type=colData(smallscRNAseq)$source_name_ch1)
rownames(anno.df) <- colnames(smallscRNAseq)
pheatmap(log2(assay(smallscRNAseq)+1), annotation_col = anno.df,
         show_rownames = FALSE, show_colnames = FALSE,
         main="before netSmooth")

pheatmap(log2(assay(smallscRNAseq.sm.sce)+1), annotation_col = anno.df,
         show_rownames = FALSE, show_colnames = FALSE,
         main="after netSmooth")

2.1 Optimizing the smoothing parameter alpha

By default, the parameter alpha will be optimized using a robust clustering statistic. Briefly, this approach will try different clustering algorithms and/or parameters and find clusters that can be reproduced with different algorithms. The netSmooth() function will try different alpha values controlled by additional arguments to maximize the number of samples in robust clusters.

Now, we smooth the expression values using automated alpha optimization and plot the heatmaps of raw and smooth versions.

smallscRNAseq.sm.se <- netSmooth(smallscRNAseq, smallPPI, alpha='auto')
smallscRNAseq.sm.sce <- SingleCellExperiment(
    assays=list(counts=assay(smallscRNAseq.sm.se)),
    colData=colData(smallscRNAseq.sm.se)
)

pheatmap(log2(assay(smallscRNAseq.sm.sce)+1), annotation_col = anno.df,
         show_rownames = FALSE, show_colnames = FALSE,
         main="after netSmooth (optimal alpha)")

3 Getting robust clusters from data

There is no standard method especially 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 have implemented a robust clustering routine based on clusterExperiment. The clusterExperiment framework for robust clustering is based on consensus clustering of clustering assignments obtained from different views of the data and different clustering algorithms. 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. We implemented a clustering framework using the components of clusterExperiment and different dimensionality reduction methods.

We can directly use the robust clustering function robustClusters.

yhat <- robustClusters(smallscRNAseq, combineManyMinSize=2, combineManyProportion=.9)$clusters
## Picked dimReduceFlavor: tsne
## Warning: 'combineMany' is deprecated.
## Use 'makeConsensus' instead.
## See help("Deprecated")
## Note: Merging will be done on ' combineMany ', with clustering index 1
yhat.sm <- robustClusters(smallscRNAseq.sm.se, combineManyMinSize=2, combineManyProportion=.9)$clusters
## Picked dimReduceFlavor: pca
## Warning: 'combineMany' is deprecated.
## Use 'makeConsensus' instead.
## See help("Deprecated")
## Note: Merging will be done on ' combineMany ', with clustering index 1
cell.types <- colData(smallscRNAseq)$source_name_ch1
knitr::kable(
  table(cell.types, yhat), caption = 'Cell types and `robustClusters` in the raw data.'
)
Table 1: Cell types and robustClusters in the raw data
-1 1
2-cell blastomere 0 3
4-cell blastomere 1 3
8-cell blastomere 6 4
8-cell embryo 1 0
morula 0 3
oocyte 0 3
pronucleus 1 2
zygote 0 2 
knitr::kable(
  table(cell.types, yhat.sm), caption = 'Cell types and `robustClusters` in the smoothed data.'
)
Table 1: Cell types and robustClusters in the smoothed data
-1 1 2 3 4 5
2-cell blastomere 1 2 0 0 0 0
4-cell blastomere 1 0 0 0 3 0
8-cell blastomere 8 0 0 0 0 2
8-cell embryo 1 0 0 0 0 0
morula 0 0 0 3 0 0
oocyte 0 1 2 0 0 0
pronucleus 0 2 1 0 0 0
zygote 0 1 1 0 0 0

A cluster assignment of -1 indicates that the cell could not be placed in a robust cluster, and has consequently been omitted. We see that the clusters are completely uninformative in the raw data, while the smoothed data at least permitted the robustClusters procedure to identify a subset of the 8-cell blastomeres as a separate cluster.

4 Deciding for the best dimension reduction method for visualization and clustering

The robustClusters() function works by clustering samples in a lower dimension embedding using either PCA or t-SNE. Different single cell datasets might respond better to different dimensionality reduction techniques. 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. The entropy in the 2D embedding is a measure for the information captured by it. We pick the embedding with the highest information content. pickDimReduction() function implements this strategy and returns the best embedding according to this strategy.

Below, we pick the best embedding for our example dataset and plot scatter plots for different 2D embedding methods.

plotPCA(smallscRNAseq, colour_by='source_name_ch1') + ggtitle("PCA plot")
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.

plotTSNE(smallscRNAseq, colour_by='source_name_ch1') + ggtitle("tSNE plot")
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.

The pickDimReduction method picks the dimensionality reduction method which produces the highest entropy embedding:

pickDimReduction(smallscRNAseq)
## [1] "tsne"

5 Frequently asked questions

5.0.1 How can I make smoothing faster ?

Make sure you compile R with openBLAS or variants that are faster.

5.0.2 What happens if all the genes are not in my network ?

The smoothing will only be done using the genes in the network then unsmoothed genes will be attached to the gene expression matrix.

sessionInfo()
## R version 3.5.1 Patched (2018-07-12 r74967)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.8-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.8-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] pheatmap_1.0.10             netSmooth_1.2.0            
##  [3] clusterExperiment_2.2.0     bigmemory_4.5.33           
##  [5] scater_1.10.0               ggplot2_3.1.0              
##  [7] SingleCellExperiment_1.4.0  SummarizedExperiment_1.12.0
##  [9] DelayedArray_0.8.0          BiocParallel_1.16.0        
## [11] matrixStats_0.54.0          Biobase_2.42.0             
## [13] GenomicRanges_1.34.0        GenomeInfoDb_1.18.0        
## [15] IRanges_2.16.0              S4Vectors_0.20.0           
## [17] BiocGenerics_0.28.0         BiocStyle_2.10.0           
## 
## loaded via a namespace (and not attached):
##   [1] uuid_0.1-2               backports_1.1.2         
##   [3] copula_0.999-18          NMF_0.21.0              
##   [5] plyr_1.8.4               lazyeval_0.2.1          
##   [7] splines_3.5.1            entropy_1.2.1           
##   [9] pspline_1.0-18           rncl_0.8.3              
##  [11] gridBase_0.4-7           digest_0.6.18           
##  [13] foreach_1.4.4            htmltools_0.3.6         
##  [15] viridis_0.5.1            magrittr_1.5            
##  [17] memoise_1.1.0            cluster_2.0.7-1         
##  [19] doParallel_1.0.14        limma_3.38.0            
##  [21] annotate_1.60.0          stabledist_0.7-1        
##  [23] prettyunits_1.0.2        colorspace_1.3-2        
##  [25] blob_1.1.1               xfun_0.4                
##  [27] dplyr_0.7.7              crayon_1.3.4            
##  [29] RCurl_1.95-4.11          bigmemory.sri_0.1.3     
##  [31] genefilter_1.64.0        bindr_0.1.1             
##  [33] phylobase_0.8.4          survival_2.43-1         
##  [35] iterators_1.0.10         ape_5.2                 
##  [37] glue_1.3.0               registry_0.5            
##  [39] gtable_0.2.0             zlibbioc_1.28.0         
##  [41] XVector_0.22.0           kernlab_0.9-27          
##  [43] Rhdf5lib_1.4.0           prabclus_2.2-6          
##  [45] DEoptimR_1.0-8           HDF5Array_1.10.0        
##  [47] abind_1.4-5              scales_1.0.0            
##  [49] mvtnorm_1.0-8            DBI_1.0.0               
##  [51] edgeR_3.24.0             rngtools_1.3.1          
##  [53] bibtex_0.4.2             Rcpp_0.12.19            
##  [55] viridisLite_0.3.0        xtable_1.8-3            
##  [57] progress_1.2.0           bit_1.1-14              
##  [59] mclust_5.4.1             glmnet_2.0-16           
##  [61] httr_1.3.1               RColorBrewer_1.1-2      
##  [63] fpc_2.1-11.1             modeltools_0.2-22       
##  [65] pkgconfig_2.0.2          XML_3.98-1.16           
##  [67] flexmix_2.3-14           nnet_7.3-12             
##  [69] locfit_1.5-9.1           labeling_0.3            
##  [71] softImpute_1.4           howmany_0.3-1           
##  [73] tidyselect_0.2.5         rlang_0.3.0.1           
##  [75] reshape2_1.4.3           AnnotationDbi_1.44.0    
##  [77] munsell_0.5.0            tools_3.5.1             
##  [79] RSQLite_2.1.1            ade4_1.7-13             
##  [81] evaluate_0.12            stringr_1.3.1           
##  [83] yaml_2.2.0               knitr_1.20              
##  [85] bit64_0.9-7              robustbase_0.93-3       
##  [87] purrr_0.2.5              dendextend_1.9.0        
##  [89] bindrcpp_0.2.2           nlme_3.1-137            
##  [91] whisker_0.3-2            xml2_1.2.0              
##  [93] compiler_3.5.1           beeswarm_0.2.3          
##  [95] tibble_1.4.2             RNeXML_2.1.2            
##  [97] pcaPP_1.9-73             stringi_1.2.4           
##  [99] gsl_1.9-10.3             highr_0.7               
## [101] RSpectra_0.13-1          lattice_0.20-35         
## [103] trimcluster_0.1-2.1      Matrix_1.2-14           
## [105] pillar_1.3.0             ADGofTest_0.3           
## [107] BiocManager_1.30.3       zinbwave_1.4.0          
## [109] cowplot_0.9.3            data.table_1.11.8       
## [111] bitops_1.0-6             R6_2.3.0                
## [113] bookdown_0.7             gridExtra_2.3           
## [115] vipor_0.4.5              codetools_0.2-15        
## [117] MASS_7.3-51              assertthat_0.2.0        
## [119] MAST_1.8.0               rhdf5_2.26.0            
## [121] pkgmaker_0.27            rprojroot_1.3-2         
## [123] withr_2.1.2              GenomeInfoDbData_1.2.0  
## [125] locfdr_1.1-8             diptest_0.75-7          
## [127] hms_0.4.2                grid_3.5.1              
## [129] tidyr_0.8.2              class_7.3-14            
## [131] rmarkdown_1.10           DelayedMatrixStats_1.4.0
## [133] Rtsne_0.13               numDeriv_2016.8-1       
## [135] ggbeeswarm_0.6.0