---
title: "BEARscc: Using spike-ins to assess single cell cluster robustness"
author: "David T. Severson"
output: BiocStyle::pdf_document
vignette: >
    %\VignetteIndexEntry{Vignette Title}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)
```
\newpage


```{r prepare_libaries, eval=TRUE, include=FALSE}
library("data.table")
library("SingleCellExperiment")
```

# Introduction

## Scope
Single-cell transcriptome sequencing data are subject to substantial 
technical variation and batch effects that can confound the classification 
of cellular sub-types. Unfortunately, current clustering algorithms do not 
account for this uncertainty. To address this shortcoming, we have developed 
a noise perturbation algorithm called **BEARscc** that is designed 
to determine the extent to which classifications by existing clustering 
algorithms are robust to observed technical variation.

**BEARscc** makes use of spike-in measurements to model 
technical variance as a function of gene expression and technical 
dropout effects on lowly expressed genes. In our benchmarks, we found 
that BEARscc accurately models read count fluctuations and drop-out effects 
across transcripts with diverse expression levels. Applying our approach to 
publicly available single-cell transcriptome data of mouse brain and 
intestine, we have demonstrated that BEARscc identified cells that cluster 
consistently, irrespective of technical variation. For more details, see the
[manuscript on bioRxiv](http://biorxiv.org/content/early/2017/03/21/118919).

Importantly, **BEARscc** should not be considered another clustering algorithm.
Specifically, this package is designed to supply users with an organic
tool to evaluate and explore the impact of noise on their single cell 
cluster interpretations. The package provides users with a way to clarify
the precision of a single cell experiment with respect to grouping cells
into clusters that are biologically meaningful. In this way, **BEARscc** 
allows users to achieve confidence in clusters and relevant cells that 
consistently cluster together invariant to the noise inherent to a single
cell experiment. Conversely, the algorithm provides a mechanism to identify
cells and clusters which cannot be resolved given the precision of the 
experiment in conjuction with the clustering algorithm of choice. It is 
our hope that **BEARscc** will enable users to proceed with clusters, or 
biological groups, in which they are confident are robust to noise and in 
which they have an intimate understanding of those cells and clusters that
may be less precisely assigned to the putative biological role. 


## Installation
BEARscc is now available on Bioconductor and can be installed using the 
syntax below. 



```{r install_BEARscc, include=TRUE, eval=FALSE}
if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("BEARscc")
```


## Citation
BEARscc and its associated manuscript are currently under review for 
publication at a peer-reviewed journal. For now, please cite the 
bioRxiv pre-print:

    Severson, DT. Owen, RP. White, MJ. Lu, X. Schuster-Boeckler, B. 
    BEARscc determines robustness of single-cell clusters using simulated
    technical replicates. doi: https://doi.org/10.1101/118919

\newpage

# Tutorial

## Overview

BEARscc relies upon spike-in count measurements in single-cell transcriptome 
experiments to estimate experimental noise and produce simulated technical
replicates to provide a quantitative understanding of the robustness of
proposed single cell cluster labels to experimental noise. In principal, 
the algorithm is compatible with any clustering algorithm. The following
should provide users with a comprehensive tutorial of the use and utlity
of BEARscc as a tool for vetting single cell clusters with respect to 
experimental noise. 

Before getting started, we need to load some example single cell data.BEARscc 
is equipped with a set of sample data for the purpose of testing functions, 
examples in help files, and this nifty tutorial. The data may be loaded as 
follows:

```{r load_data, eval=TRUE}
library("BEARscc")
data("BEARscc_examples")
```

The loaded file `BEARscc_examples` is equipped with separate `data.frame` 
objects including ERCC spike-in observations (`ERCC.counts.df`), 
endogenous count observations (`data.counts.df`), and the expected or 
actual spike-in concentrations (`ERCC.meta.df`) as well as a 
`SingleCellExperiment` object that contains all of the above seperate
components (`BEAR_examples.sce`) as shown below:

```{r display_data, eval=TRUE}
head(ERCC.counts.df[,1:2])

head(data.counts.df[,1:2]) 

head(ERCC.meta.df)

BEAR_examples.sce 
```

In the event we were working with a new set of data, the spike-in 
concentrations `data.frame` can be computed from industry reported 
concentrations and the relevant dilution protocol utilized in the 
experiment. Count tables would need to be mapped and counted with 
preferred software, and the spike-in control counts (ERCC or otherwise) 
would need to be identified from the endogenous counts. 

Below is how one would create a `SingleCellExperiment` object from 
spike-in count, endogenous count, and spike-in concentration `data.frame` 
objects. Note how we create the `observed_expression` assay object in the 
following code. This object is essential in that all estimation and simulation
of replicates occurs assuming these are the observed counts (normalized or 
otherwise). Without them BEARscc will throw an error indicating that 
`"observed_expression"` not in `names(assays(<SingleCellExperiment>))`.
Also, `data.counts.df` in this example includes both spike-in genes and 
endogenous genes. In general, we recommend spike-in genes be simulated
with endogenous genes as a control, and this will be done by default
when the `SingleCellExperiment` object is used. 

```{r create_SCEList, eval=TRUE}
BEAR.se <- SummarizedExperiment(list(counts=as.matrix(data.counts.df)))
BEAR_examples.sce<-as(BEAR.se, "SingleCellExperiment")
metadata(BEAR_examples.sce)<-list(spikeConcentrations=ERCC.meta.df)
assay(BEAR_examples.sce, "observed_expression")<-counts(BEAR_examples.sce)
isSpike(BEAR_examples.sce, "ERCC_spikes")<-grepl("^ERCC-", 
    rownames(BEAR_examples.sce))
```

Running the above code should yield a `SingleCellExperiment` object identical
to the one that loads with `data("BEARscc_examples")`.


## Building the noise model
We will now estimate the single-cell noise present in the experiment using
spike-in controls. In this tutorial, we rely upon a subsample of artificial
control data found in `BEARscc_examples`; however, users are encouraged to 
work through the tutorial with their own single cell data provided some form 
of spike-ins were included in the experiment. Building the noise models with 
BEARscc is relatively straightforward with `estimate_noiseparameters()`. We
simply provide the function with the now adequately annotated 
`SingleCellExperiment` object, `BEAR_examples.sce`. Here, the parameter 
`alpha_resolution is set to 0.1 to speed things along, but we suggest values 
between 0.001 and 0.01 be used in real applications of BEARscc.

```{r estimate_noise, eval=TRUE}
BEAR_examples.sce <- estimate_noiseparameters(BEAR_examples.sce,
    max_cumprob=0.9999, alpha_resolution = 0.1, bins=10,
    write.noise.model=FALSE, file="BEAR_examples")
```

Several options exist for `estimate_noiseparameters()`. These are 
fully documented in the help page `?estimate_noiseparameters`. 

## Simulating technical replicates

Following estimation of noise, the parameters computed are then used to 
generate a simulate a technical replicate. Here we will simulate replicates 
on our local computer, but frequently users will want to utilize the methods 
described in Section 2.4. Notably the necessary parameters are conveniently 
stored in the `metadata` of our `SingleCellExpression` object,
`BEAR_examples.sce` following noise estimation, and so we simply run: 
    
```{r simulate_replicates, eval=TRUE}
    BEAR_examples.sce <- simulate_replicates(BEAR_examples.sce,
        max_cumprob=0.9999, n=10)
```

Recall that `BEAR_examples.sce` is our `SingleCellExperiment` object that we 
recently annoated with model parameters describing experimental noise using the
function `estimate_noiseparameters()` and note that the variable `n` is the 
desired number of simulated technical replicates (e.g. 10). Finally, the
`maxcum_prob` is identical to its use in the noise estimation. If the user
deviated from the default parameter, it is highly recommended that this value
be identical to the value utilized in `estimate_noiseparmaters()` The resulting
object is a list, where each element is a simulated technical replicate, 
and one element is the original counts matrix. 


## Simulation of replicates for larger datasets
For larger datasets, we set `write.noise.model=TRUE` when running 
`estimate_noiseparameters()`and copy the written bayesian drop-out and 
noise estimate files with the observed counts table to a high performance 
computing environment. The following code provides an example:
    
```{r estimate_noise_HPC, eval=FALSE, include=TRUE}
BEAR_examples.sce <- estimate_noiseparameters(BEAR_examples.sce,
                                    write.noise.model=TRUE,
                                    file="tutorial_example",
                                    model_view=c("Observed","Optimized"))
```

After running the above code, then within the current working directory 
(if unsure use `getwd()`), we should find the two tab-delimited files
that together completely describe the BEARscc noise model. These are the 
parameters describing the mixed model of technical variation 
(`tutorial_example_parameters4randomize.xls`, see Section 3.1.1) and the 
parameters describing the drop-out model 
(`tutorial_example_bayesianestimates.xls`, see Section 3.1.2). In addition, 
we should find our `"observed_expression"`matrix in the form of 
`tutorial_example_counts4clusterperturbation.xls`. Note that the `xls` 
subscript just allows users to quickly open these tab-delimited files in 
Microsoft Excel if desired, but these can be readily viewed on the terminal 
or in simple text editors as well. 

With the original counts file and noise model prepared, we then copy these
files to our high performance compute cluster. The following code provides
a sense of how to proceed once these files have been copied to a high
performance cluster; however, the job submission structure of each user's
environment will dictate the precise syntax for the following procedure.

The script `HPC_generate_noise_matrices` contains an analogous 
`simulate_replicates()` for a high performance computational node.
To utilize these functions for simulating technical replicates on a 
cluster, please install BEARscc on the relevant cluster. 
The user should write an R script to load the BEARscc library and
run the clustering. The following code provides a suggested format for
both calling the R script with a bash job script and the relevant R 
invocation of BEARscc and may also be
found as a stand alone script in `inst/example/HPC_run_example.R`.

Our cluster utilizes a job submission format that interacts seamlessly
with bash code; therefore, the `$SGE_TASKID` represents an array id for  
jobs to conveniently generate 100 simulated technical replicates in a
single job array. In any case, this variable should be treated as the
index for the simulated technical replicate as we recommend from
experience that users generate 50 to 100 such simulated technical replicates 
to reach a stable noise consensus matrix solution. 

The following bash code could be included in one such job script:

```{bash run_sim_rep_HPC, eval=FALSE, include=TRUE}
Rscript --vanilla HPC_run_example.R $SGE_TASK_ID 
```

Noting that the file `HPC_run_example.R` contains the following 
suggested code to run BEARscc:

```{r HPC_example, eval=FALSE, include=TRUE}
library("BEARscc")

#### Load data ####
ITERATION<-commandArgs(trailingOnly=TRUE)[1]
counts.df<-read.delim("tutorial_example_counts4clusterperturbation.xls")
#filter out zero counts to speed up algorithm
counts.df<-counts.df[rowSums(counts.df)>0,]
probs4detection<-fread("tutorial_example_bayesianestimates.xls")
parameters<-fread("tutorial_example_parameters4randomize.xls")

#### Simulate replicates ####
counts.error<-HPC_simulate_replicates(counts_matrix=counts.df,
    dropout_parameters=dropout_parameters,
    spikein_parameters=spikein_parameters)
write.table(counts.error, file=paste("simulated_replicates/",
    paste(ITERATION,"sim_replicate_counts.txt",sep="_"),
    sep=""), quote =FALSE, row.names=TRUE)
#########
```

The script generates seperate simulated technical replicate files, 
which can be loaded into R as a list for clustering or, in the case of more 
computationally intense clustering algorithms, re-clustered 
individually in a high performance compute environment. 

## Forming a noise consensus
After generating simulated technical replicates, these should be re-clustered 
    using the clustering method applied to the original dataset. 
    For simplicity, here we use hierarchical clustering on a euclidean 
    distance metric to identify two clusters. In our experience, 
    some published clustering algorithms are sensitive to cell order, 
    so we suggest scrambling the order of cells for each noise iteration 
    as we do below in the function, `recluster()`. The function below will
    serve our purposes for this tutorial, but BEARscc may be used with 
    any clustering algorithm (e.g. SC3, RaceID2, or BackSPIN). 

To quickly recluster a list of simulated technical replicates, we define a 
reclustering function:

```{r recluster, include=TRUE, eval=TRUE}
recluster <- function(x) {
    x <- data.frame(x)
    scramble <- sample(colnames(x), size=length(colnames(x)), replace=FALSE)
    x <- x[,scramble]
    clust <- hclust(dist(t(x), method="euclidean"),method="complete")
    clust <- cutree(clust,2)
    data.frame(clust)
}
``` 

First, we need to determine how the observed data clusters under our algorithm
of choice (e.g. `recluster()`, which is a simple hiearchical clustering
for illustrative purposes). These are the clusters that can be evaluated by 
`BEARscc` and compared to `BEARscc` meta-clustering:

```{r original_clusters, include=TRUE, eval=TRUE}
OG_clusters<-recluster(data.frame(assay(BEAR_examples.sce, 
    "observed_expression")))
colnames(OG_clusters)<-"Original"
```

We then apply the function `recluster()` to all simulated technical replicates 
and the original counts matrix and manipulate the list into a `data.frame`. 

```{r recluster_go, include=TRUE, evalue=TRUE}
cluster.list<-lapply(metadata(BEAR_examples.sce)$simulated_replicates, 
    `recluster`)
clusters.df<-do.call("cbind", cluster.list)
colnames(clusters.df)<-names(cluster.list)
```

Note: If running clustering algorithms on a seperate high performance cluster, 
    the user should retrieve labels and format as a `data.frame` of 
    cluster labels, where the last column must be the original cluster labels 
    derived from the observed count data. As an example, examine the file,
    [inst/example/example_clusters.tsv](inst/example/example_clusters.tsv).

Using the cluster labels file generated by clustering simulated technical
replicates on a high performance compute environment or with `recluster()` 
as described above , we can generate a noise consensus matrix as shown below.
An illustrative three rows and columns of an example noise consensus matrix 
are displayed:

```{r compute_consensus, eval=TRUE, include=TRUE}
noise_consensus <- compute_consensus(clusters.df)
head(noise_consensus[,1:3], n=3)
```

Using the `aheatmap()` function in the `NMF` library, the consensus matrix 
result of 10 simulated replicates by BEARscc:

To reproduce the plot run:
```{r plot_noise_consensus, eval=TRUE, include=TRUE}
library("NMF")
aheatmap(noise_consensus, breaks=0.5)
```

Although, 10 simulated replicates is sparse (we recommend 50 to 100), we can
already see that these samples likely belong to a single cluster. Indeed, 
these samples were prepared from a single ground truth of dilute whole tissue
brain RNA-seq data from two batches of experimental data. If desired, we could
annotate the above heatmap with relevant metadata concerning sample batch, 
origin, etc.

## Evaluating the noise consensus
In order to further interpret the noise consensus, we have defined three 
cluster (and analagous cell) metrics. Stability indicates the propensity for a 
putative cluster to contain the same cells across noise-injected counts 
matrices. Promiscuity indicates a tendency for cells in a putative cluster 
to associate with other clusters across noise-injected counts matrices. 
Score represents the promiscuity substracted from the stability. 

We have found it useful to inform the optimal number of clusters in terms 
of resiliance to noise by examining these metrics by cutting hierarchical 
clustering dendograms of the noise consensus and comparing the results to 
the original clustering labels. To do this create a vector containing 
each number of clusters one wishes to examine (the function automatically 
determines the results for the dataset as a single cluster) and then 
cluster the consensus with various cluster numbers using `cluster_consensus()`:

```{r cluster_consensus, include=TRUE, eval=TRUE}
vector <- seq(from=2, to=5, by=1)
BEARscc_clusts.df <- cluster_consensus(noise_consensus,vector)
```

We add the original clustering to the `data.frame` to evaluate its 
robustness as well as the suggested BEARscc clusters:

```{r add_original, include=TRUE, eval=TRUE}
BEARscc_clusts.df <- cbind(BEARscc_clusts.df, 
    Original=OG_clusters)
```

### Understanding robustness at the cluster level

Now we can compute cluster metrics for each of the BEARscc cluster 
number scenarious and the original clustering; indeed, any cluster labels
of the users choosing could be supplied to vet with the information provided
by the noise consensus. We accomplish this by running the command and
displaying the first 5 rows of the resulting melted `data.frame`:

```{r compute_cluster_scores, include=TRUE, eval=TRUE}
cluster_scores.df <- report_cluster_metrics(BEARscc_clusts.df, 
    noise_consensus, plot=FALSE)
head(cluster_scores.df, n=5)
```

Above is a melted `data.frame` that displays the name of each cluster, 
    the size of each cluster, the metric (Score, Promiscuity, Stability), 
    the value of each metric for the respective cluster and clustering, 
    the clustering in question (1,2,...,Original), whether the cluster 
    consists of only one cell, and finally the mean of each metric 
    across all clusters in a clustering.  

In the previous example, we display all metrics for generating 1
clusters from the data given the previously computed noise consensus from
10 simulated technical replicates and the same for the score metric 
generating 2 clusters from the data. Importantly, by definition, 1 cluster 
scenarios have a promiscuity value of 0. Observe that the score for the
cluster in the 1 cluster scenario is much larger than either cluster of the
2 cluster scenario, and that this is reflected in the average clustering 
column. 

We can examine the BEARcc metrics for the original cluster using `ggplot2`
and by subsetting the `data.frame` to the original cluster and the score
metric:

```{r plot_clusterscores_original, include=TRUE, eval=TRUE}
library("ggplot2")
library("cowplot")
original_scores.df<-cluster_scores.df[
    cluster_scores.df$Clustering=="Original",]
ggplot(original_scores.df[original_scores.df$Metric=="Score",], 
    aes(x=`Cluster.identity`, y=Value) )+
    geom_bar(aes(fill=`Cluster.identity`), stat="identity")+
    xlab("Cluster identity")+ylab("Cluster score")+
    ggtitle("Original cluster scores")+guides(fill=FALSE)
```

We can see that this initial clustering is terrible, which makes sense given
the ground truth consists of a single biological entity. The stability and 
promiscuity metrics bear this out:

```{r plot_clusterstability_original, include=TRUE, eval=TRUE}
ggplot(original_scores.df[original_scores.df$Metric=="Stability",], 
    aes(x=`Cluster.identity`, y=Value) )+
    geom_bar(aes(fill=`Cluster.identity`), stat="identity")+
    xlab("Cluster identity")+ylab("Cluster stability")+
    ggtitle("Original cluster stability")+guides(fill=FALSE)
```

The high stability exhibited in the example above is not suprising as samples
in this example should have strong association with one another. However, the
promuscuity below reveals the reason for the low score:

```{r plot_clusterpromiscuity_original, include=TRUE, eval=TRUE}
ggplot(original_scores.df[original_scores.df$Metric=="Promiscuity",], 
    aes(x=`Cluster.identity`, y=Value) )+
    geom_bar(aes(fill=`Cluster.identity`), stat="identity")+
    xlab("Cluster identity")+ylab("Cluster promiscuity")+
    ggtitle("Original cluster promiscuity")+guides(fill=FALSE)
```

Despite the high stability, the samples within each cluster have high 
association with cells in the other cluster, which results in a high
promiscuity reported from the noise consensus. As a net result, the scores
in the original 2 cluster case for each cluster are subpar. Again, this
is consistent with the ground truth of the example data. 

### Understanding robustness at the sample level
Completely analagous to cluster metrics, the extent to which cells belong
within a given cluster may be evaluated with respect to the noise consensus.
Below we demonstrate how to compute the cell metrics and display 4 cells for
illustrative purposes. 

```{r compute_cell_scores, include=TRUE, eval=TRUE}
cell_scores.df <- report_cell_metrics(BEARscc_clusts.df, noise_consensus)
head(cell_scores.df, n=4)
```

The output is a melted `data.frame` that displays the name of each cluster 
    to which the cell belongs, the cell label, the size of each cluster, 
    the metric (Score, Promiscuity, Stability), the value for each metric, 
    and finally the clustering in question (1,2,...,Original).

As with cluster metrics, these results can be plotted to visualize cells 
in the context of the original clusters using `ggplot2`. The score metric
is plotted below to illustrate this:

```{r plot_cellscore_original, include=TRUE, eval=TRUE}
original_cell_scores.df<-cell_scores.df[cell_scores.df$Clustering=="Original",]
ggplot(original_cell_scores.df[original_cell_scores.df$Metric=="Score",],
    aes(x=factor(`Cluster.identity`), y=Value) )+
    geom_jitter(aes(color=factor(`Cluster.identity`)), stat="identity")+
    xlab("Cluster identity")+ylab("Cell scores")+
    ggtitle("Original clustering cell scores")+guides(color=FALSE)
```

### Using BEARscc to inform cluster number, $k$, choice
While BEARscc certainly does not claim to provide a definitive solution to the 
question concerning the number of clusters by any means, we provide what we 
believe to be a useful perspective on the matter. Specifically, we have found 
that by examining the cluster metrics across various hiearchical clusterings 
of the BEARscc noise consensus, the cluster number $k$ with the highest score 
tends to provide a number of clusters that resembles ground truth more closely 
than simple gene sampling or relevant algorithms along as evidenced by 
experiments  in control samples, c. elegans , and murine brain data (see our 
manuscript on [bioRxiv](http://biorxiv.org/content/early/2017/03/21/118919). 
Importantly, this only provides a heuristic for determining cluster number $k$ 
that takes into count the inherent noise of the single cell experiment, and the
heuristic is dependent upon the clustering algorithm of choosing (the better
the utilized clustering algorithm, the better the BEARscc $k$ heuristic) and
represents a form of "meta-clustering" rather than a new way to determine the
number of clusters, $k$, per se. As an illustration of utilizing this 
heuristic, we plot the average cluster score values for various cluster number
scenarios below:

```{r choosing_k, include=TRUE, eval=TRUE}
ggplot(cluster_scores.df[cluster_scores.df$Metric=="Score",], 
    aes(x=`Clustering`, y=Value) )+ geom_boxplot(aes(fill=Clustering,  
    color=`Clustering`))+ xlab("Clustering")+
    ylab("Cluster score distribution")+
    ggtitle("Original cluster scores for each clustering")+
    guides(fill=FALSE, color=FALSE)
```

As we can see, the 1 cluster scenario provides the best cluster score
as determined by the noise consensus. As mentioned previously,
the single cluster solution resembles the biological ground truth. 

\newpage
# Algorithm and theory
In order to simulate technical replicates, BEARscc first builds a statistical 
model of expression variance and drop-out frequency, which is dependent only 
on observed gene expression. The parameters of this model are estimated from 
spike-in counts.  Expression-dependent variance is approximated by fitting 
read counts of each spike-in transcript across cells to a mixture model 
comprised of a Poisson and negative binomial distribution (Section 3.1.1). 
The drop-out model (Section 3.1.2) in BEARscc has two distinct parts: 
the \emph{drop-out injection distribution} models the likelihood that a 
given transcript concentration will result in a drop-out, and the 
\emph{drop-out recovery distribution} models the likelihood that an observed 
drop-out resulted from a given transcript concentration. The drop-out 
injection distribution is taken to be the observed drop-out rate in spike-in 
controls as a function of actual spike-in transcript concentration. 
This distribution is then used to estimate the drop-out recovery distribution 
density via Bayes'  theorem and a an empirically informed set of priors and 
assumptions. Briefly, BEARscc utilizes the drop-out injection distribution 
and the number of observed zeroes for each endogenous gene to infer a 
gene-specific probability distribution describing the likelihood that an 
observed drop-out should in fact have been some non-zero value, given 
the drop-out rate of the endogenous gene. This entire process is facilitated 
by the function, `estimate_noiseparameters()`.

In the second step, BEARscc generates simulated technical replicates by 
applying the models described in the first step (Section 3.2). For every 
observed count in the range of values where drop-outs occurred amongst the 
spike-in transcripts, BEARscc uses the drop-out injection distribution from 
Step 1 to determine whether to convert the count to zero.  For observations 
where the count is zero, the drop-out recovery distribution is used to 
estimate a new value, given the overall drop-out frequency for the gene 
(Section 3.2). Next, BEARscc substitutes all values larger than zero with a 
value generated from the derived model of expression variability, 
parameterized to the observed count for that gene. This procedure can then be 
repeated any number of times to generate a collection of simulated technical 
replicates. This step is carried out by `create_noiseinjected_counts()`.

In the third step, the simulated technical replicates are then re-clustered, 
using exactly the same method as for the observed data; this re-clustering 
for each simulated technical replicate is described as an 
\emph{association matrix} where each element indicates whether two cells 
share a cluster identity (1) or cluster apart from each other (0). 
The association matrices for each simulated technical replicate are averaged 
to form a noise consensus matrix that can be easily interpreted (Section 3.3).
This is accomplished with the function `compute_consensus()`. Each element of 
the noise consensus matrix represents the fraction of simulated technical 
replicates that, upon applying the clustering method of choice, resulted in 
two cells clustering together (the \emph{association frequency}). Then, the 
functions `report_cell_metrics()` and `report_cluster_metrics()` may be used 
to explore and quantitate the noise consensus matrix at the cell sample 
and cluster levels, respectively. 

## Noise estimation
As mentioned previously, BEARscc uses spike-ins to estimate the noise of the 
experiment for the purpose of producing simulated technical replicates. 
BEARscc models overall technical variation with a mixture-model 
(Section 3.1.1) and inferred drop-out effects (Section 3.1.2) independently 
using the spike-in observations. However, a single function in 
BEARscc `estimate_noiseparameters()` accomplishes this task.  

### Estimating transcript variation
Technical variance is modeled in BEARsccby fitting a single parameter 
mixture model, $Z(c)$, to the spike-ins' observed count distributions. 
The noise model is fit independently for each spike-in transcript and 
subsequently regressed onto spike-in mean expression to define a 
generalized noise model. This is accomplished in three steps: 

\begin{enumerate}
    \item Define a mixture model composed of \emph{poisson} and 
\emph{negative binomial} random variables: 
        \begin{equation}
            Z \sim (1-\alpha)*Pois(\mu)+\alpha*NBin(\mu,\sigma)
        \label{eq:mixture_model}\end{equation}
    \item Empirically fit the parameter, $\alpha_{i}$, in a spike-in specific 
mixture-model, $Z_{i}$, to the observed distribution of counts for each ERCC 
spike-in transcript, $i$, where $\mu_i$ and $\sigma_i$ are the observed mean 
and variance of the given spike-in. The parameter, $\alpha_{i}$, is chosen 
such that the error between the observed and mixture-model is minimized.
    \item  Generalize the mixture-model by regressing $\alpha_{i}$ parameters
and the observed variance, $\sigma_i$, onto the observed spike-in mean 
expression, $\mu_i$. Thus the mixture model describing the noise observed 
in ERCC transcripts is defined solely by $\mu$, which is treated as the 
count transformation parameter, $c$, in the generation of simulated 
technical replicates. 
\end{enumerate}

In step 2, a mixture model distribution is defined for each spike-in, $i$:
\begin{equation}
    Z_{i}(\alpha_{i},\mu_{i},\sigma_{i})\sim (1-\alpha_{i})*Pois(\mu_i )+
\alpha_{i}*NBin(\mu_i,\sigma_i ).
\end{equation}
The distribution, $Z_i$, is fit to the observed counts of the respective 
spike-in, where $\alpha_i$ is an empirically fitted parameter, such that 
the $\alpha_i$ minimizes the difference between the observed count 
distribution of the spike-in and the respective fitted model, $Z_i$. 
Specifically, for each spike-in transcript, $\mu_i$ and $\sigma_i$ 
are taken to be the mean and standard deviation, respectively, of the 
observed counts for spike-in transcript, $i$. Then, $\alpha_i$ is computed 
by empirical parameter optimization; $\alpha_i$ is taken to be the 
$\alpha_{i,j}$ in the mixture-model, 
\begin{equation}
    Z_{i,j}(\alpha_{i,j},\mu_i,\sigma_i)\sim (1-\alpha_{i,j})*
    Pois(\mu_i )+\alpha_{i,j}*NBin(\mu_i,\sigma_i ),
\end{equation}
found to have the least absolute total difference between the observed count 
density and the density of the fitted model, $Z_i$. In the case of ties, the 
minimum $\alpha_{i,j}$ is chosen.

In step 3, $\alpha(c)$ is then defined with a linear fit, 
$\alpha_i=\alpha*log2(\mu_i)+b$. $\sigma(c)$  was similarly defined, 
$log2(\sigma_i) = \alpha*log2(\mu_i)+b$. In this way, the observed 
distribution of counts in spike-in transcripts defines the single 
parameter mixture-model, $Z(c)$, used to transform counts during 
generation of simulated technical replicates:

\begin{equation}
    Z(c) \sim (1-\alpha(c))*Pois(c)+\alpha(c)*NBin(c,\sigma(c))
\label{eq:noisefunction}\end{equation}


During technical replicate simulation, the parameter $c$ is set to the 
observed count value, $a$, and the transformed count in the simulated 
replicate was determined by sampling a single value from $Z(c=a)$. 


### Defining the drop-out models
A model of the drop-outs is developed by BEARscc in order to inform the 
permutation of zeros during noise injection. The observed zeros in spike-in 
transcripts as a function of actual transcript concentration and Bayes' 
theorem are used to define two models: the 
\emph{drop-out injection distribution} and the 
\emph{drop-out recovery distribution}. 

The drop-out injection distribution is described by  $Prob(X=0 | Y=y)$, 
where $X$ is the distribution of observed counts and $Y$ is the distribution 
of actual transcript counts; the density is computed by regressing the 
fraction of zeros observed in each sample, $D_i$, for a given spike-in, $i$, 
onto the expected number spike-in molecules in the sample, $y_i$, 
e.g. $D=a*y+b$. Then, $D$ describes the density of zero-observations 
conditioned on actual transcript number, $y$, or $Prob(X=0 | Y=y)$. 
Notably, each gene was treated with an identical density distribution 
for drop-out injection.

In contrast, the density of the drop-out recovery distribution, 
$Prob(Y_j=y | X_j=0)$, is specific to each gene, $j$, where $X_j$ is the 
distribution of the observed counts and $Y_j$ is the distribution of actual 
transcript counts for a given gene. The gene-specific drop-out recovery 
distribution is inferred from drop-out injection distribution using Bayes' 
theorem and a prior.  This is accomplished in 3 steps:

\begin{enumerate}
    \item For the purpose of applying Bayes' theorem, the gene-specific 
distribution, $Prob(X_j=0 | Y_j=y)$, is taken to be the the drop-out 
injection density for all genes, $j$. 
    \item The probability that a specific transcript count is present 
in the sample, $Prob(Y_j=y)$, is a necessary, but empirically unknowable
prior. Therefore, the prior was defined using the law of total probability, 
an assumption of uniformity, and the probability that a zero was observed 
in a given gene, $Prob(X_j=0)$. The probability, $Prob(X_j=0)$, is taken to 
be the fraction of observations that are zero for a given gene. BEARscc does 
this in order to better inform the density estimation of the gene-specific 
drop-out recovery distribution.
    \item The drop-out recovery distribution density is then computed 
by applying Bayes' theorem:

\begin{equation}
    Prob(Y_j=y | X_j=0)=\frac{Prob(X_j=0 | Y_j=y)*Prob(Y_j=y)}{Prob(X_j=0)},
\label{eq:bayesianeq}\end{equation}

\end{enumerate}

In the second step, the law of total probability, an assumption of uniformity, 
and the fraction of zero observations in a given gene are leveraged to define 
the prior, $Prob(Y_j=y)$. First, a threshold of expected number of 
transcripts, $k$ in $Y$, is chosen such that $k$ was the maximum value for 
which the drop-out injection density was non-zero. Next, uniformity is 
assumed for all expected number of transcript values, $y$ greater than zero 
and less than or equal to $k$; that is $Prob(Y_j=y)$ is defined to be some 
constant probability, $n$. Furthermore, $Prob(Y_j=y)$ is defined to be 0 for 
all $y>k$. In order to inform $Prob(Y_j=y)$ empirically, $Prob(Y_j=0)$ and 
$n$ are derived by imposing the law of total probability 
(Equation \ref{eq:totalprob}) and unity (Equation \ref{eq:unity}) 
yielding a system of equations:

\begin{equation}
    Prob(X_j=0)=\Sigma_{y=0}^{k-1}{ (Prob(X_j=0 | Y_j=y)*Prob(Y_j=y))} 
\label{eq:totalprob}\end{equation}

\begin{equation}
    \Sigma_{y=0}^{k-1}{ Prob(Y_j=y)}=Prob(Y_j=0)+(k-1)*{n}=1  
\label{eq:unity}\end{equation}

The probability that a zero is observed given there are no transcripts in 
the sample, $Prob(X_j=0 | Y_j=0)$, is assumed to be 1. With the preceding 
assumption, solving for $Prob(Y_j=0)$ and $n$ give:

\begin{equation}
n=\frac{1-Prob(Y_j=0)}{k-1}
\label{eq:nsolved}\end{equation}

\begin{equation}
Prob(Y_j=0)=\frac{Prob(X_j=0)-\frac{1}{k-1}*
\Sigma_{y=1}^{k-1}( Prob(X_j=0 | Y_j=y)) }{(1-\frac{1}{k-1}*
\Sigma_{y=1}^{k-1} (Prob(X_j=0 | Y_j=y)) }
\label{eq:py0solved}\end{equation} \\

In this way, $Prob(Y_j=0)$ is defined by (Equation \ref{eq:nsolved}) 
for $y$ in $Y_j$ less than or equal to $k$ and greater than zero, 
and defined by (Equation \ref{eq:py0solved}) for $y$ in $Y_j$ equal to zero. 
For $y$ in $Y_j$ greater than $k$, the prior $Prob(Y_j=y)$ is defined to 
be equal to zero.

In the third step, the previously computed prior, $Prob(Y_j=y)$, the fraction 
of zero observations in a given gene, $Prob(X_j=0)$, and the drop-out 
injection distribution, $Prob(X_j=0 | Y_j=y)$, are utilized to estimate, with 
Bayes's theorem, the density of the drop-out recovery distribution, 
$Prob(Y_j=y | X_j=0)$. During the generation of simulated technical replicates 
for zero observations and count observations less than or equal to $k$, 
values are sampled from the drop-out recovery and injection distributions as 
described in the pseudocode of the BEARscc algorithm for 
simulating technical replicates.


## Simulating technical replicates
Simulated technical replicates were generated from the noise mixture-model and 
two drop-out models. For each gene, the count value of each sample is 
systematically transformed using the mixture-model, $Z(c)$, and the 
drop-out injection, $Prob(X=0 | Y=y)$, and recovery, $Prob(Y_j=y | X_j=0)$, 
distributions in order to generate simulated technical replicates as indicated 
by the following pseudocode:


```
FOR EACH gene, $j$
    FOR EACH count, $c$ 
        IF $c=0$ 
            $n \leftarrow SAMPLE$ one count, y, from $Prob(Y_j=y | X_j=0)$
            IF $n=0$ 
                $c \leftarrow 0$ 
            ELSE
                $c \leftarrow SAMPLE$ one count from $Z(n)$
            ENDIF
        ELSE
            IF $c\leq k$ 
                $dropout \leftarrow TRUE$ with probability, $Prob(X=0 | Y=k)$ 
                IF $dropout=TRUE$ 
                    $c \leftarrow 0$ 
                ELSE 
                    $c \leftarrow SAMPLE$ one count from $Z(c)$ 
                ENDIF 
            ELSE 
                $c \leftarrow SAMPLE$ one count from $Z(c)$ 
            ENDIF 
        ENDIF 
        RETURN $c$
    DONE
DONE 
```

\newpage
# List of functions

## `estimate_noiseparameters()`
### Description
Estimates the drop-out model and technical variance 
    from spike-ins present in the sample.

For greater detail, please see help file `?estimate_noiseparameters()`.

### Usage
```R
data(BEARscc_examples)
BEAR_examples.sce <- estimate_noiseparameters(BEAR_examples.sce,
alpha_resolution=0.25, write.noise.model=FALSE)

BEAR_examples.sce
```

### Output
The resulting output of `estimate_noiseparameters()` is a long list,
    which is enumerated in the function's package help page.

### Note
The above usage is for execution of `simulate_replicates` 
    on a local machine. To save results as files for us of
    `prepare_probabilities` on high performance computing
    environment, then use:
```R 
estimate_noiseparameters(BEAR_examples.sce,
    write.noise.model=TRUE, alpha_resolution=0.25,
    file="noise_estimation", model_view=c("Observed","Optimized"))
```


## `simulate_replicates()`
### Description
Computes BEARscc simulated technical replicates from the 
    previously estimated noise parameters computed with the 
    function `estimate_noise_parameters()`.
    
For greater detail, please see help file `?simulate_replicates()`.

### Usage
```R
data(analysis_examples)
BEAR_examples.sce<-simulate_replicates(BEAR_examples.sce, n=3)

BEAR_examples.sce
```
### Output
The resulting object is a list of counts data, where each element of the list 
    is a `data.frame` of the counts representing a BEARscc simulated technical 
    replicate. For further details refer to the function help page.
    
### Note
This function is the in-package analog of the high-performance 
    computing function,`prepare_probabilities`.

## `HPC_simulate_replicates()`
### Description
The high-performance computing function analog to 
    `simulate_replicates()`.
    
### Usage
Please refere to section 2.4. 

### Output
The resulting objects would normally be output to a tab-delimited file, where
    each file results from a `data.frame` of the counts representing a BEARscc
    simulated technical replicate. 

### Note
This function has no help file, but is referred to in the section 2.4
    of this document on simulating for larger datasets.


## `compute_consensus()`
### Description
Computes the consensus matrix using a `data.frame` of cluster labels across 
    different BEARscc simulated technical replicates.The consensus matrix is
    is a visual and quantitaive representation of the clustering variation on 
    a cell-by-cell level created by using cluster labels to compute 
    the number of times any given pair of cells associates in the same 
    cluster; this forms the 'noise consensus matrix'. Each element of 
    this matrix represents the fraction of simulated technical replicates 
    in which two cells cluster together (the 'association frequency'), 
    after using a clustering method of the user's choice to generate a 
    data.frame of clustering labels. This consensus matrix may be used 
    to compute BEARscc metrics at both the cluster and cell level.
    
For greater detail, please see help file `?compute_consensus()`.

### Usage
```R
data("analysis_examples")
noise_consensus <- compute_consensus(clusters.df)
noise_consensus
```

### Output
When the number of samples are \emph{n}, then the noise consensus 
    resulting from this function is an \emph{n} x \emph{n} matrix 
    describing the fraction of simulated technical replicates in 
    which each cell of the experiment associates with another cell.

## `cluster_consensus()`
### Description
This function will perform hierarchical clustering on the noise 
    consensus matrix allowing the user to investigate the appropriate 
    number of clusters, \emph{k}, considering the noise within the 
    experiment. Frequently one will want to assess multiple possible 
    cluster number situations at once. In this case it is recommended 
    that one use a `lapply` in conjunction with a vector of all 
    biologically reasonable cluster numbers to fulfill the task of 
    attempting to identify the optimal cluster number.
    
For greater detail, please see help file `?cluster_consensus()`.

### Usage
```R
data(analysis_examples)
vector <- seq(from=2, to=5, by=1)
BEARscc_clusts.df <- cluster_consensus(consensus_matrix=noise_consensus, 
    vector)
BEARscc_clusts.df
```

### Output
The output is a vector of cluster labels based on hierarchical clustering
    of the noise consensus. In the event that a vector is supplied for
    number of clusters in conjunction with `lapply`, then the output
    is a data.frame of the cluster labels for each of the various number
    of clusters deemed biologically reasonable by the user.
    
## `report_cell_metrics()`
### Description
To quantitatively evaluate the results, three metrics are calculated 
    from the noise consensus matrix: 'stability' is the average 
    frequency with which a cell within a cluster associates with other 
    cells within the same cluster across simulated replicates; 
    'promiscuity' measures the average association frequency of a cell 
    within a cluster with the \emph{n} cells outside of the cluster 
    with the strongest association with the cell in question; 
    and 'score' is the difference between 'stability' and 'promiscuity'. 
    Importantly, 'score' reflects the overall "robustness" of a given cell's 
    assignment to a user-provided cluster label with respect to technical 
    variance. Together these metrics provide a quantitative measure of 
    the extent to which cluster labels provided by the user are invariant 
    across simulated technical replicates.
    
For greater detail, please see help file `?report_cell_metrics()`. 

### Usage
```R
data(analysis_examples)
cell_scores.dt <- report_cell_metrics(BEARscc_clusts.df, noise_consensus)
cell_scores.dt
```

### Output
A melted `data.frame` describing the BEARscc metrics for each cell, 
    where the columns are enumerated in the help file.

## `report_cluster_metrics()`
### Description
To quantitatively evaluate the results, three metrics are calculated 
    from the noise consensus matrix: 'stability' is the average 
    frequency with which cells within a cluster associate with each 
    other across simulated replicates; 'promiscuity' measures the 
    association frequency between cells within a cluster and those 
    outside of it; and 'score' is the difference between 'stability' 
    and 'promiscuity'. Importantly, 'score' reflects the overall "robustness" 
    of a cluster to technical variance. Together these metrics provide
    a quantitative measure of the extent to which cluster labels provided
    by the user are invariant across simulated technical replicates.
    
For greater detail, please see help file `?report_cluster_metrics()`. 

### Usage
```R
data(analysis_examples)
cluster_scores.dt <- report_cluster_metrics(BEARscc_clusts.df,noise_consensus,
    plot=TRUE, file="example")
cluster_scores.dt
```

### Output
A melted `data.frame` describing the BEARscc metrics for each cluster, 
    where the columns are enumerated in the help file.

\newpage
#Example data
Within the package there are data subsampled from single cell sequencing 
    protocol applied to water samples containing ERCC spike-ins (blanks) 
    and dilute RNA from brain whole tissue (brain) discussed at length 
    in in a manuscript on 
    [bioRxiv](http://biorxiv.org/content/early/2017/03/21/118919)

## `BEARscc_examples`
### Description
A toy dataset for applying BEARscc functions as described in the README on 
    https://bitbucket.org/bsblabludwig/bearscc.git.These data are a subset of 
    observations made by Drs. Michael White and Richard Owen 
    in the Xin Lu Lab. Samples were sequenced by the Wellcome Trust Center 
    for Genomics, Oxford, UK. These data are available in full with GEO 
    accession number, GSE95155.

For greater detail, please see help file `?BEARscc_examples`.
    
### Usage
```R
data("BEARscc_examples")
```

## `analysis_examples`
### Description
BEARscc downstream example objects: The `analysis_examples` Rdata object 
    contains downstream data objects for use in various help pages for 
    dynamic execution resulting from running tutorial in README and 
    vignette on `BEARscc_examples`. The objects are a result of applying 
    BEARscc functions as described in the README found at
    https://bitbucket.org/bsblabludwig/bearscc.git.

For greater detail, please see help file `?analysis_examples`.

### Usage
```R
data("analysis_examples")
```

\newpage
# License

This software is made available under the terms of the 
[GNU General Public License v3](http://www.gnu.org/licenses/gpl-3.0.html)

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 
    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND 
    NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR 
    ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR 
    OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING 
    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 
    DEALINGS IN THE SOFTWARE.


