Introduction to netSmooth package
2023-04-25
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.
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.5smallscRNAseq.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, makeConsensusMinSize=2, makeConsensusProportion=.9)$clusters## Picked dimReduceFlavor: umap## 6 parameter combinations, 0 use sequential method, 0 use subsampling method
## Running Clustering on Parameter Combinations...
## done.## Note: Merging will be done on ' makeConsensus ', with clustering index 1yhat.sm <- robustClusters(smallscRNAseq.sm.se, makeConsensusMinSize=2, makeConsensusProportion=.9)$clusters## Picked dimReduceFlavor: pca## 6 parameter combinations, 0 use sequential method, 0 use subsampling method
## Running Clustering on Parameter Combinations...
## done.## Warning in (function (A, nv = 5, nu = nv, maxit = 1000, work = nv + 7, reorth =
## TRUE, : You're computing too large a percentage of total singular values, use a
## standard svd instead.## Note: Merging will be done on ' makeConsensus ', with clustering index 1cell.types <- colData(smallscRNAseq)$source_name_ch1
knitr::kable(
table(cell.types, yhat), caption = 'Cell types and `robustClusters` in the raw data.'
)| -1 | 1 | 2 | |
|---|---|---|---|
| 2-cell blastomere | 1 | 0 | 2 |
| 4-cell blastomere | 2 | 0 | 2 |
| 8-cell blastomere | 6 | 4 | 0 |
| 8-cell embryo | 1 | 0 | 0 |
| morula | 0 | 3 | 0 |
| oocyte | 1 | 0 | 2 |
| pronucleus | 1 | 0 | 2 |
| zygote | 2 | 0 | 0 |
knitr::kable(
table(cell.types, yhat.sm), caption = 'Cell types and `robustClusters` in the smoothed data.'
)| -1 | 1 | 2 | 3 | 4 | |
|---|---|---|---|---|---|
| 2-cell blastomere | 1 | 0 | 2 | 0 | 0 |
| 4-cell blastomere | 3 | 0 | 1 | 0 | 0 |
| 8-cell blastomere | 6 | 0 | 0 | 2 | 2 |
| 8-cell embryo | 1 | 0 | 0 | 0 | 0 |
| morula | 0 | 3 | 0 | 0 | 0 |
| oocyte | 0 | 0 | 3 | 0 | 0 |
| pronucleus | 0 | 0 | 3 | 0 | 0 |
| zygote | 0 | 0 | 2 | 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.
smallscRNAseq <- runPCA(smallscRNAseq, ncomponents=2)
smallscRNAseq <- runTSNE(smallscRNAseq, ncomponents=2)
smallscRNAseq <- runUMAP(smallscRNAseq, ncomponents=2)
plotPCA(smallscRNAseq, colour_by='source_name_ch1') + ggtitle("PCA plot")
plotTSNE(smallscRNAseq, colour_by='source_name_ch1') + ggtitle("tSNE plot")
plotUMAP(smallscRNAseq, colour_by='source_name_ch1') + ggtitle("UMAP plot")
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 4.3.0 RC (2023-04-13 r84269)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB 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
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] pheatmap_1.0.12 netSmooth_1.20.0
## [3] clusterExperiment_2.20.0 bigmemory_4.6.1
## [5] scater_1.28.0 ggplot2_3.4.2
## [7] scuttle_1.10.0 SingleCellExperiment_1.22.0
## [9] SummarizedExperiment_1.30.0 Biobase_2.60.0
## [11] GenomicRanges_1.52.0 GenomeInfoDb_1.36.0
## [13] IRanges_2.34.0 S4Vectors_0.38.0
## [15] BiocGenerics_0.46.0 MatrixGenerics_1.12.0
## [17] matrixStats_0.63.0 BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.3.0 bitops_1.0-7
## [3] tibble_3.2.1 XML_3.99-0.14
## [5] lifecycle_1.0.3 edgeR_3.42.0
## [7] doParallel_1.0.17 lattice_0.21-8
## [9] MASS_7.3-59 MAST_1.26.0
## [11] magrittr_2.0.3 limma_3.56.0
## [13] sass_0.4.5 rmarkdown_2.21
## [15] jquerylib_0.1.4 yaml_2.3.7
## [17] NMF_0.26 cowplot_1.1.1
## [19] zinbwave_1.22.0 DBI_1.1.3
## [21] RColorBrewer_1.1-3 ade4_1.7-22
## [23] abind_1.4-5 zlibbioc_1.46.0
## [25] Rtsne_0.16 purrr_1.0.1
## [27] RCurl_1.98-1.12 GenomeInfoDbData_1.2.10
## [29] ggrepel_0.9.3 howmany_0.3-1
## [31] irlba_2.3.5.1 genefilter_1.82.0
## [33] annotate_1.78.0 DelayedMatrixStats_1.22.0
## [35] codetools_0.2-19 DelayedArray_0.26.0
## [37] xml2_1.3.3 tidyselect_1.2.0
## [39] RNeXML_2.4.11 locfdr_1.1-8
## [41] farver_2.1.1 ScaledMatrix_1.8.0
## [43] viridis_0.6.2 jsonlite_1.8.4
## [45] BiocNeighbors_1.18.0 phylobase_0.8.10
## [47] survival_3.5-5 iterators_1.0.14
## [49] foreach_1.5.2 tools_4.3.0
## [51] progress_1.2.2 Rcpp_1.0.10
## [53] glue_1.6.2 gridExtra_2.3
## [55] xfun_0.39 dplyr_1.1.2
## [57] HDF5Array_1.28.0 withr_2.5.0
## [59] BiocManager_1.30.20 fastmap_1.1.1
## [61] rhdf5filters_1.12.0 fansi_1.0.4
## [63] entropy_1.3.1 digest_0.6.31
## [65] rsvd_1.0.5 R6_2.5.1
## [67] colorspace_2.1-0 RSQLite_2.3.1
## [69] utf8_1.2.3 tidyr_1.3.0
## [71] generics_0.1.3 data.table_1.14.8
## [73] FNN_1.1.3.2 prettyunits_1.1.1
## [75] httr_1.4.5 uwot_0.1.14
## [77] pkgconfig_2.0.3 gtable_0.3.3
## [79] blob_1.2.4 registry_0.5-1
## [81] XVector_0.40.0 htmltools_0.5.5
## [83] bookdown_0.33 scales_1.2.1
## [85] png_0.1-8 bigmemory.sri_0.1.6
## [87] knitr_1.42 reshape2_1.4.4
## [89] rncl_0.8.7 uuid_1.1-0
## [91] nlme_3.1-162 cachem_1.0.7
## [93] rhdf5_2.44.0 stringr_1.5.0
## [95] parallel_4.3.0 vipor_0.4.5
## [97] softImpute_1.4-1 AnnotationDbi_1.62.0
## [99] pillar_1.9.0 grid_4.3.0
## [101] vctrs_0.6.2 BiocSingular_1.16.0
## [103] beachmat_2.16.0 xtable_1.8-4
## [105] cluster_2.1.4 beeswarm_0.4.0
## [107] evaluate_0.20 magick_2.7.4
## [109] cli_3.6.1 locfit_1.5-9.7
## [111] compiler_4.3.0 rlang_1.1.0
## [113] crayon_1.5.2 rngtools_1.5.2
## [115] labeling_0.4.2 plyr_1.8.8
## [117] ggbeeswarm_0.7.1 stringi_1.7.12
## [119] viridisLite_0.4.1 gridBase_0.4-7
## [121] BiocParallel_1.34.0 munsell_0.5.0
## [123] Biostrings_2.68.0 Matrix_1.5-4
## [125] hms_1.1.3 sparseMatrixStats_1.12.0
## [127] bit64_4.0.5 Rhdf5lib_1.22.0
## [129] KEGGREST_1.40.0 highr_0.10
## [131] kernlab_0.9-32 memoise_2.0.1
## [133] bslib_0.4.2 bit_4.0.5
## [135] ape_5.7-1