---
title: "DEGage Tutorial"
output: 
    rmarkdown::html_vignette:
        number_sections: true
        toc: true
        toc_depth: 2
       
        
vignette: >
  %\VignetteIndexEntry{my-vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

This package allows for the differential expression analysis of scRNA-seq and other NGS count data. It employs a novel family of discrete distributions for describing the difference of two NB distributions (named DOTNB). DEGage takes the raw counts of scRNA-seq as inputs, and thus avoids introducing artificially bias in normalization steps in current methods. A workflow is shown as follows.

# Setup Instructions

## Installation
To install DEGage on your local machine, ensure devtools is installed and run the following commands
```{r eval=FALSE}
library(devtools)
install_github("chenyongrowan/DEGage")
```

## Dependencies 

DEGage require the following packages to be installed and up to date. The installation will be performed automatically if these packages are not already installed. 
```{r eval=FALSE}
MASS >= 7.3-58.2
pscl >= 1.5.5
hypergeo >= 1.2-13
parallel >= 4.2.2
doParallel >= 1.0.17
stats >= 4.2.2
SingleR >= 2.0.0
Seurat >= 4.3.0
parallel >= 4.3.0
doParallel >= 1.0.17
foreach >=1.5.2
```

# DEGage Family Functions

Currently, the DEGage package has five primary functions: DEGage, DEGage_preprocess, DEGage_multitest, DEGage_complete, and DEGage_simulate.

## DEGage
DEGage performs pairwise differential expression analysis on scRNA-seq data. The input is typically a dataframe where columns contain samples and rows contain genes. Counts do not need to be normalized prior to use with DEGage. 
```{r eval=FALSE}
DEGage(counts, group, perm.preprocess = TRUE,
       gene.filter.threshold = 1, nperms = 2000,
       nsubsample = NA, perm.pval = 0.1, ncores = 4)
```

### Parameters
- `counts`: An object or path to directory containing read counts. The following inputs are accepted:   
    + A dataframe where columns are samples and rows are genes
    + A SeuratObject with an assay titled "RNA" containing counts
    + A SingleCellExperiment containing an assay titled "counts" 
    + A path to an mtx directory
- `group`: A factor of a numeric vector that indicates which samples belong to which pairwise condition. An example is shown below under the Example Usage section  
- `perm.preprocess`: A boolean indicating whether a permutation test is used to pre-filter genes  
- `gene.filter.threshold`: A value between 0-1 indicating the maximum proportion of zero counts a gene can have before being filtered out. A value of 1 means only genes with all zero counts are filtered out, where as a value of 0.5 would filter out genes where half the counts are zeros  
- `nsubsample`: The number of cells to sub-sample for each group. If left as NA, the default is 125+0.05n, where n is the number of samples in the data set.   
- `nperms`: The number of permutations performed during the permutation test  
- `perm.pval`: A p value used by the permutation test to pre-filter genes
- `ncores`: The number of cores to use for parallel computing

### Output 
A dataframe containing the following information about each gene is output:  
- `r1`: The r parameter for the NB distribution fit to the first condition of counts  
- `p1`: The p parameter for the NB distribution fit to the first condition of counts   
- `mu1`: The mean of the counts for the first condition  
- `r2`: The r parameter for the NB distribution fit to the second condition of counts  
- `p2`: The p parameter for the NB distribution fit to the second condition of counts  
- `mu2`: The mean of the counts for the second condition  
- `fit.method`: Indicates whether a gene was fit according to a NB or ZINB distribution  
- `permPvals`: p-value generated from a pre-filtering permutation test. It is recommended to filter genes according to a permutation p-value of 0.1 to minimize false positive rates.   
- `pval`: p-values generated according to DOTNB  
- `FDR`: FDR adjusted p-values  

## DEGage_multitest
Similar to DEGage, except it can perform an indefinite number of pairwise comparisons. 
```{r eval=FALSE}
DEGage_multitest(counts, group, perm.preprocess = FALSE,
                 gene.filter.threshold = 1, nperms = 2000,
                 nsubsample = NA, perm.pval = 0.1, ncores = 4,
                 writing.dir = NULL)
```

### Parameters
- `counts`: An object or path to directory containing read counts. The following inputs are accepted:   
    + A dataframe where columns are samples and rows are genes
    + A SeuratObject with an assay titled "RNA" containing counts
    + A SingleCellExperiment containing an assay titled "counts" 
    + A path to an mtx directory 
- `group`: A factor of a numeric vector that indicates which samples belong to which pairwise condition. An example is shown below under the Example Usage section  
- `perm.preprocess`: A boolean indicating whether a permutation test is used to pre-filter genes  
- `gene.filter.threshold`: A value between 0-1 indicating the maximum proportion of zero counts a gene can have before being filtered out. A value of 1 means only genes with all zero counts are filtered out, where as a value of 0.5 would filter out genes where half the counts are zeros  
- `nsubsample`: The number of cells to sub-sample for each group. If left as NA, the default is 125+0.05n, where n is the number of samples in the data set.   
- `nperms`: The number of permutations performed during the permutation test  
- `perm.pval`: A p value used by the permutation test to pre-filter genes
- `ncores`: The number of cores to use for parallel computing
- `writing.dir`: A directory where the results of each comparison are to be written to. It is highly recommended to provide a directory for large numbers of comparisons.   

### Output
If writing.dir is NULL, a list of Dataframes structured in the same way as the output of DEGage() that each contain a single comparison are returned. If writing.dir is not null and is a valid directory, nothing is output, and the results can be retrieved with list.files("directory")   

## DEGage_preprocess
Performs a simple scRNA-seq count pre-processing workflow with Seurat and generates automatic cell type annotations with SingleR. Is primarily used in DEGage_complete, however, it is available for use individually.  
```{r eval=FALSE}
DEGage_preprocess(input, dir.type = 'mtx', min.nFeatureRNA = 200,
                  max.nFeatureRNA = 8000, mt.percent = .2,
                  cell.annotations = TRUE,
                  celltype.min = 20, output = "Seurat")
```
### Parameters
- `input`: An object or path to directory containing read counts. The following inputs are accepted:   
    + A dataframe where columns are samples and rows are genes
    + A SeuratObject with an assay titled "RNA" containing counts
    + A SingleCellExperiment containing an assay titled "counts" 
    + A path to an mtx directory
- `dir.type`: If input is a directory, this is used to specify the directory type. If input is not a directory, 
then this variable will be ignored. Directory types that are currently supported include:
    + "mtx"
- `min.nFeatureRNA`: The minimum number of transcripts that must be present in a cell before it is filtered   
- `max.nFeatureRNA`: The maximum number of transcripts that must be present in a cell before it is filtered    
- `mt.percent`: A value between 0 and 1 that indicates the maxiumum proprotion of mitochondrial genes permitted before a cell is filtered    
- `cell.annotations`: A Boolean indicating whether or not cell type annotations are performed. For datasets with very few cell, it is recommended to set this to false, because singleR cannot annotate small datasets.    
- `celltypemin`: After cell type annotation are generated, cells of celltypes less than "celltypemin" will be filtered out.     
- `output`:  Determines the format of the output. The following outputs are currently supported:    
    + "Seurat"
    + "df"
    
### Output
Depending on the value of 'output', either a data frame or SeuratObject containing processed counts will be output  

## DEGage_complete
Takes unprocessed and non-annotated scRNA-seq counts, automatically processes and generates cell type annotations for them, then performs differential expression analysis on the groups identified through clustering.    
```{r eval=FALSE}
DEGage_complete(input, dir.type = 'mtx',
                min.nFeatureRNA = 200, max.nFeatureRNA = 8000, 
                mt.percent = .2, cell.annotations = TRUE,
                celltype.min = 20, perm.preprocess = FALSE,
                gene.filter.threshold = 1, nperms = 2000, 
                nsubsample = NA, perm.pval = 0.1,
                ncores = 4, writing.dir = NULL)
```
### Input
- `input`: An object or path to directory containing read counts. The following inputs are accepted:   
    + A dataframe where columns are samples and rows are genes
    + A SeuratObject with an assay titled "RNA" containing counts
    + A SingleCellExperiment containing an assay titled "counts" 
    + A path to an mtx directory
- `dir.type`: If input is a directory, this is used to specify the directory type. If input is not a directory, 
then this variable will be ignored. Directory types that are currently supported include:
    + "mtx"
- `min.nFeatureRNA`: The minimum number of transcripts that must be present in a cell before it is filtered   
- `max.nFeatureRNA`: The maximum number of transcripts that must be present in a cell before it is filtered    
- `mt.percent`: A value between 0 and 1 that indicates the maxiumum proprotion of mitochondrial genes permitted before a cell is filtered    
- `cell.annotations`: A Boolean indicating whether or not cell type annotations are performed. For datasets with very few cell, it is recommended to set this to false, because singleR cannot annotate small datasets.    
- `celltypemin`: After cell type annotation are generated, cells of celltypes less than "celltypemin" will be filtered out.  
- `perm.preprocess`: A boolean indicating whether a permutation test is used to pre-filter genes  
- `gene.filter.threshold`: A value between 0-1 indicating the maximum proportion of zero counts a gene can have before being filtered out. A value of 1 means only genes with all zero counts are filtered out, where as a value of 0.5 would filter out genes where half the counts are zeros  
- `nsubsample`: The number of cells to sub-sample for each group. If left as NA, the default is 125+0.05n, where n is the number of samples in the data set.   
- `nperms`: The number of permutations performed during the permutation test  
- `perm.pval`: A p value used by the permutation test to pre-filter genes
- `ncores`: The number of cores to use for parallel computing
- `writing.dir`: A directory where the results of each comparison are to be written to. For DEGage_complete() specfically, there are likely to be many comparisons so it is *highly* recommended to provide a writing directory  

### Output
If writing.dir is nul, a list of dataframes containing the DEGage() output for each comparison is output. If a writing directory is provided, the results are written to files and NULL is returned. 

## DEGage_simulation
Generates very simplified simulated scRNA-Seq counts following an NB distribution with pre-defined proportions of dropouts.   
```{r eval=FALSE}
DEGage_Simulation(ngenes, ndegs, cellgroups, lfc = 1, 
                  prop.zeros = .3, seed = NULL, ncores = 4)

```

### Input
- `ngenes`: The number of genes to be simulated in total   
- `ndegs`: The number of genes to be differentially expressed   
- `cellgroups`: A factor of a vector that is the same length as the number of cells that are to be simulated. The vector must contain two different integers to correspond to two sample conditions. See an example in the "Example Usage" section   
- `lfc`: Either an integer, or a vector that is the same length as 'ndegs' indicating the log2 fold change that differtially expressed genes are to be up/down regulated by   
- `prop.zeros`: The proportion of counts that are to be replaced with zeros to simulate dropouts   
- `seed`: Sets seed for random number generation for reproducible results   
- `ncores`: The number of corse to use for parallel computing to increase speed   

### Output

A dataframe where cells are represented by columns and genes are represented by rows. Differentially expressed genes have the prefix "DE", while equivalently expressed genes have the prefix "EE"   
  

# Example Usage

In this section, we will detail how to use DEGage functions

First, we will simulate a small data frame to pass through DEGage() using DEGage simulate: 
```{r eval=FALSE}
library(DEGage)
cellgroups <- factor( c( rep(1,5), rep(2,5) ) )
df <- DEGage_Simulation(ngenes = 5, ndegs = 2, cellgroups = cellgroups)
```

Next, we will pass this through DEGage:
```{r eval=FALSE}
results <- DEGage(counts = df, group = cellgroups)
```

To test DEGage_multitest, we will simulate a second dataframe of counts, merge them together, and pass them through DEGage_multitest(): 
```{r eval=FALSE}
cellgroups2 <- factor( c( rep(3,5), rep(4,5) ) )
df2 <- DEGage_Simulation(ngenes = 5, ndegs = 2, cellgroups = cellgroups2)

df <- cbind(df, df2)
cellgroups <- factor(c(cellgroups, cellgroups2))

multitest.results <- DEGage_multitest(df, cellgroups)
```
