---
title: "scDDboost"
author: 
- name: Xiuyu Ma
  email:  watsonforfun@gmail.com
output: 
  BiocStyle::html_document:
  toc_float: true
vignette: >
  %\VignetteIndexEntry{scDDboost Tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r options, include=FALSE, echo=FALSE}
library(BiocStyle)
knitr::opts_chunk$set(warning=FALSE, error=FALSE, message=FALSE)
```

# Installation

The package can be installed from bioconductor
```{r,eval=FALSE}
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("scDDboost")
```
Issue can be reported at "https://github.com/wiscstatman/scDDboost/issues"


# Introduction
scDDboost scores evidence of a gene being differentially distributed(DD) across two conditions for single cell RNA-seq data. Higher resolution brings several chanllenges for analyzing the data, specifically, the distribution of gene expression tends to have high prevalence of zero and multi-modes.
To account for those characteristics and utilizing some biological intuition,
we view the expression values sampled from a pool of cells mixed by distinct cellular subtypes blind to condition label. Consequently, the distributional change can be fully determined by the
the change of subtype proportions. One tricky part is that not any change of proportions will lead to a distributional change. Given that some genes could be equivalent expressed across several subtypes, even the individual subytpe proportion may differ between conditions but as long as the aggregated proportions over those subtypes remain the same between conditions, it will not introduce different distribution.
For example


```{r,echo=FALSE}
library(ggplot2)
df <- data.frame(
  group = as.factor(rep(c("subtype1", "subtype2", "subtype3"),2)),
proportion = c(1/3,1/3,1/3,1/6,1/2,1/3),
  label = as.factor(rep(c("condition1","condition2"),each = 3))
  )
bp<- ggplot(df, aes(x=proportion,y = "", fill=group))+
geom_bar(width = 2, stat = "identity") + facet_grid(cols = vars(label)) 
pie <- bp + coord_polar("x", start=0) + xlab("subtype proportions") + ylab("")
pie +
  theme(panel.background = element_rect(
    fill = 'white', colour = 'black'),
    panel.grid.minor.x = element_line(size = 0.5),
    panel.grid.minor.y = element_line(size = 0.5),
    panel.grid.major.x = element_line(size = 0.5),
    panel.grid.major.y = element_line(size = 0.5),
    panel.grid.major = element_line(colour = "grey"))
```

Proportions of subtypes 1 and 2 changed between the 2 conditions. 
The gene is not DD if subtype 1 and 2 have the same expression level

```{r,echo = F}
x1 <- rnorm(33,10,2)
x2 <- rnorm(33,10,2)
x3 <- rnorm(33,20,2)

x4 <- rnorm(17,10,2)
x5 <- rnorm(50,10,2)
x6 <- rnorm(33,20,2)

df <- data.frame(
  group = c(rep(c("subtype1", "subtype2", "subtype3"),each = 33), c(rep("subtype1",17),rep("subtype2",50),rep("subtype3",33))),
  conditions = factor(c(rep(1,99),rep(2,100))),
  values = c(x1,x2,x3,x4,x5,x6)
  )
p <- ggplot(df,aes(x = conditions,y = values)) + geom_violin() + geom_jitter(aes(colour = group)) + xlab("conditions") + 
  ylab("gene expressions")
  

p + theme(panel.background = element_rect(
    fill = 'white', colour = 'black'),
    panel.grid.minor.x = element_line(size = 0.5),
    panel.grid.minor.y = element_line(size = 0.5),
    panel.grid.major.x = element_line(size = 0.5),
    panel.grid.major.y = element_line(size = 0.5),
    panel.grid.major = element_line(colour = "grey"))
```

For subtype 1 and 2 have different expression level, there is different distribution
```{r, echo = F}
x1 <- rnorm(33,10,2)
x2 <- rnorm(33,20,2)
x3 <-  rnorm(33,30,2)

x4 <- rnorm(17,10,2)
x5 <- rnorm(50,20,2)
x6 <- rnorm(33,30,2)

df <- data.frame(
  group = c(rep(c("subtype1", "subtype2", "subtype3"),each = 33), c(rep("subtype1",17),rep("subtype2",50),rep("subtype3",33))),
  conditions = factor(c(rep(1,99),rep(2,100))),
  values = c(x1,x2,x3,x4,x5,x6)
  )
p <- ggplot(df,aes(x = conditions,y = values)) + geom_violin() + geom_jitter(aes(colour = group)) + xlab("conditions") + 
  ylab("gene expressions")
  

p + theme(panel.background = element_rect(
    fill = 'white', colour = 'black'),
    panel.grid.minor.x = element_line(size = 0.5),
    panel.grid.minor.y = element_line(size = 0.5),
    panel.grid.major.x = element_line(size = 0.5),
    panel.grid.major.y = element_line(size = 0.5),
    panel.grid.major = element_line(colour = "grey"))





```




# Posterior probability of a gene being DD 


`pdd` is the core function developed to quantify the
posterior probabilities of DD for input genes. 
Let's look at an example,

```{r}
suppressMessages(library(scDDboost))
```
Next, we load the toy simulated example a \textit{SingleCellExperiment} object that we will use for
identifying and classifying DD genes.
```{r}
data(sim_dat)
```
Verify that this object is a member of the SingleCellExperiment class and that it contains 200
cells and 1000 genes. The colData slot (which contains a dataframe of metadata for the
cells) should have a column that contains the biological condition or grouping of interest. In
this example data, that variable is the `condition` variable. Note that the input gene set needs
to be a matrix of normalized counts.
We run the function `pdd` 
```{r}
data_counts <- SummarizedExperiment::assays(sim_dat)$counts
conditions <- SummarizedExperiment::colData(sim_dat)$conditions
rownames(data_counts) <- seq_len(1000)

##here we use 2 cores to compute the distance matrix
bp <- BiocParallel::MulticoreParam(2)
D_c <- calD(data_counts,bp)

ProbDD <- pdd(data = data_counts,cd = conditions, bp = bp, D = D_c)
```
There are 4 input parameters needed to be specified by user, the dataset, the condition label, number of cpu cores used for computation and a distance matrix
of cells. Other input parameters have default settings.

## clustering of cells

We provide a default method of getting the distance matrix, archived by `calD`, in general `pdd` accept all valid distance matrix. User can also input a cluster label rather than distance matrix for the argument `D`, but the random distancing mechanism which relies on distance matrix will be disabled and `random` should be set to false.  

For the number of sutypes, we provide a default function `detK`, which consider the smallest number of sutypes such that the ratio of difference within cluster between difference between clusters become smaller than a threshold (default setting is 1). 

If user have other ways to determine $K$, $K$ should be specified in `pdd`. 
```{r}
## determine the number of subtypes
K <- detK(D_c)
```


If we set threshold to be 5% then we have estimated DD genes 
```{r}
EDD <- which(ProbDD > 0.95)
```

Notice that, pdd is actually local false discovery rate, this is a conservative estimation of DD genes. 
We could gain further power, 
let index gene by $g = 1,2,...,G$ and let $p_g = P(DD_g | \text{data})$, $p_{(1)},...,p_{(G)}$ be ranked local false discovery rate from small to large. 
To control the false discovery rate at 5\%, our positive set is those genes with the $s^*$ smallest lFDR, where
$$s^* = \text{argmax}_s\{s,\frac{\Sigma_{i = 1}^s p_{(i)}}{s} \leq 0.05\}$$
```{r}
EDD <- getDD(ProbDD,0.05)
```

Function `getDD` extracts the estimated DD genes using the above transformation. 



# Session Information {-}

```{r}
sessionInfo()
```





