library(netSmooth) library(pheatmap) library(SingleCellExperiment)

*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`

.

# All defaults knitr::include_graphics("bckgrnd.png")

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.

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) 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")

`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)")

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 yhat.sm <- robustClusters(smallscRNAseq.sm.se, makeConsensusMinSize=2, makeConsensusProportion=.9)$clusters cell.types <- colData(smallscRNAseq)$source_name_ch1 knitr::kable( table(cell.types, yhat), caption = 'Cell types and `robustClusters` in the raw data.' ) knitr::kable( table(cell.types, yhat.sm), caption = 'Cell types and `robustClusters` in the smoothed data.' )

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.

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)
```

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

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

```
sessionInfo()
```

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.