---
title: "Introduction to cinaR"
author: "[E Onur Karakaslar](https://eonurk.github.io)"
date: "`r format(Sys.Date(), '%m/%d/%Y')`"
output:
  rmarkdown::html_document:
    highlight: pygments
    toc: true 
    fig_width: 5 
    
vignette: >
  %\VignetteIndexEntry{cinaR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, echo=FALSE, results="hide"}
knitr::opts_chunk$set(tidy = FALSE,
                      cache = FALSE,
                      dev = "png",
                      message = FALSE, error = FALSE, warning = TRUE,
                      fig.dpi = 96)
```

## Quick Start

```{r}
library(cinaR)
data("atac_seq_consensus_bm")
```

Bed formatted consensus matrix (chr, start, end and samples)
```{r}
dim(bed)
```


```{r}
# bed formatted file
head(bed[,1:4])
```

Create the contrasts you want to compare, here we create contrasts for 
22 mice samples from different strains.
```{r}
# create contrast vector which will be compared.
contrasts<- c("B6", "B6", "B6", "B6", "B6", "NZO", "NZO", "NZO", "NZO", "NZO", "NZO", 
              "B6", "B6", "B6", "B6", "B6", "NZO", "NZO", "NZO", "NZO", "NZO", "NZO")
```

`cinaR` function directly computes the differentially accessible peaks.
```{r}
# If reference genome is not set hg38 will be used!
results <- cinaR(bed, contrasts, reference.genome = "mm10")
```

Now, you can access differential accessibility (DA) and enrichment results.
```{r}
names(results)
```

Inside `DA.results`, you have the consensus peaks (cp) and differentially accessible (DA) peaks.
If batch correction was run, then `cp` will be a batch-corrected consensus matrix, otherwise it is the filtered and normalized version of original consensus peaks you provided.

```{r}
names(results$DA.results)
```
There are many information `cinaR` provides such as adjusted p value, 
log fold-changes, gene names etc for each peak:
```{r}
colnames(results$DA.results$DA.peaks$B6_NZO)
```


Here is an overview of those DA peaks:
```{r }
head(results$DA.results$DA.peaks$B6_NZO[,1:5])
```
> Since the comparison is `B6_NZO`, if fold-changes are positive it means they are more accesible in B6 compared to NZO and vice versa for negative values!

and here is a little overview for enrichment analyses results:
```{r}
head(results$Enrichment.Results$B6_NZO[,c("module.name","overlapping.genes", "adj.p")])
```

### PCA Plots
You can easily get the PCA plots of the samples:
```{r}
pca_plot(results, contrasts, show.names = F)
```

You can overlay different information onto PCA plots as well!
```{r}
# Overlaid information
overlaid.info <- c("B6-18mo", "B6-18mo", "B6-18mo", "B6-18mo", "B6-18mo", 
                   "NZO-18mo", "NZO-18mo", "NZO-18mo", "NZO-18mo", "NZO-18mo", "NZO-18mo", 
                   "B6-3mo", "B6-3mo", "B6-3mo", "B6-3mo", "B6-3mo", 
                   "NZO-3mo", "NZO-3mo", "NZO-3mo", "NZO-3mo", "NZO-3mo", "NZO-3mo")
# Sample IDs
sample.names <- c("S01783", "S01780", "S01781", "S01778", "S01779", 
"S03804", "S03805", "S03806", "S03807", "S03808", 
"S03809", "S04678", "S04679", "S04680", "S04681", 
"S04682", "S10918", "S10916", "S10919", "S10921", 
"S10917", "S10920")
```

```{r }
pca_plot(results, overlaid.info, sample.names)
```

### Heatmaps


#### Differential peaks

You can see the available comparisons using:
```{r}
show_comparisons(results)
```

Then, plot the differential peaks for a selected contrast using:
```{r}
heatmap_differential(results, comparison = "B6_NZO")
```
  
Also, you can configure your heatmaps using the additional arguments of `pheatmap` function. For more information check out `?pheatmap`.
```{r}
heatmap_differential(results, comparison = "B6_NZO", show_colnames = FALSE)
```

#### Most variable peaks
You can also plot most variable 100 peaks for all samples:
```{r}
heatmap_var_peaks(results)
```

Plus, you can set the number of peaks to be used in these plots, and again you can change the additional arguments of `pheatmap` function. For more information check out `?pheatmap`.
```{r}
heatmap_var_peaks(results, heatmap.peak.count = 200, cluster_cols = F)
```

### Enrichment Plots
You can plot your enrichment results using:
```{r }
dot_plot(results)
```

if it gets too crowded, you can get rid of the irrelevant pathways as well:
```{r }
dot_plot(results, filter.pathways = T)
```

## Creating different contrasts

Note that you can further divide the resolution of contrasts, for instance this
is also a valid vector
```{r}
contrasts <- sapply(strsplit(colnames(bed), split = "-", fixed = T), 
                    function(x){paste(x[1:4], collapse = ".")})[4:25]
unique(contrasts)
```
in this case, each of them will be compared to each other which will result in
28 different comparisons.

## Running for bulk RNA-seq data

To run `cinaR` with RNA-seq experiments, just provide the count matrix, and specify the experiment type:
```{r eval=FALSE}
cinaR(matrix = count.matrix, ..., experiment.type = "RNA-Seq")
```

Note that, `count.matrix` should be in the form of $g \times (1+n)$ where $g$ is the number of genes and $n$ is the number of samples, and plus one is for gene names.

> Note that currently `cinaR` can only handle gene symbols (e.g. FOSL2, FOXA) and ensembl stable IDs (e.g. ENSG00000010404) for both human and mice!


## Running enrichment with different dataset
You can run the enrichment analyses with a custom gene set:
```{r eval=FALSE}
cinaR(..., geneset = new_geneset)
```

#### `cinaRgenesets`
Easiest way to do this is to use [cinaRgenesets](https://github.com/eonurk/cinaR-genesets) package. 
You can select your gene set of interest and just plug it into your pipeline.

#### MSigDB
You can also download different gene sets from [MSigDB site](https://www.gsea-msigdb.org/gsea/downloads.jsp). 
Note that you should use the human gene symbol versions. 

> You can use `read.gmt` function from `qusage` package to read these genesets into your current environment.

#### Custom gene sets

A `geneset` must be a `.gmt` formatted symbol file.  
You can familiarize yourself with the format by checking out :
```{r eval=FALSE}
# default geneset to be used
data("VP2008")
```

> If you have gene and pathway names in a `data.frame`, you can use `split` function to create your own `.gmt` formatted gene sets e.g. `split(df$genes, df$pathways)`.

## Selecting different reference genomes
For now, `cinaR` supports 3 genomes for human and mice models:

- `hg38`
- `hg19`
- `mm10`

You can set your it using `reference.genome` argument.

## Batch Effect Correction
If you suspect your data have unknown batch effects, you can use:
```{r eval=FALSE}
cinaR(..., batch.correction = T)
```
This option will run [Surrogate Variable Analysis](http://bioconductor.org/packages/release/bioc/html/sva.html) (SVA) and try to
adjust your data for unknown batch effects. If however, you already know the batches of the samples, you can simply set the `batch.information` argument as well. This will not run the SVA but just add the batches to design matrix.

```{r eval=FALSE}
# runs SVA
cinaR(..., batch.correction = T)

# runs SVA with 2 surrogate variables
cinaR(..., batch.correction = T, sv.number = 2)

# adds only batch information to the design matrix! (does not run SVA)
# batch.information should be number a vector where
# the length of it equals to the number of samples.
cinaR(..., batch.correction = T, batch.information = c(rep(0, 11), rep(1,11)))
```

> Reminder - In our example data we have 22 samples
  
  
## Adding extra covariates
Sometimes, one might want to add additional covariates to adjust the design matrix further, which affects the down-stream analyses. For instance, ages or sexes of the samples could be additional covariates. 
To account for those:

```{r eval=FALSE}
# Ages of the samples could be not in biological interests but should be accounted for!
cinaR(..., additional.covariates = c((18, 11), (3, 11)))

# More than one covariate for instance, sex and age
sex.info <- c("M", "F", "M", "F", "F", "F", "F", "F", "M", "M", "M", 
              "F", "F", "M", "M", "M", "F", "F", "M", "M", "F", "M")

age.info <- c((18, 11), (3, 11)
covs <- data.frame(Sex = sex.info, Age = age.info)

cinaR(..., additional.covariates = covs)
```

## Saving DA peaks to excel

Setting `save.DA.peaks = TRUE` in `cinaR` function will create a `DApeaks.xlsx` file
in the current directory. This file includes all the comparisons in different tabs.
Additionally, you can set the path/name of the file using 
`DA.peaks.path` argument after setting `save.DA.peaks = TRUE`.

For instance, 

```{r eval=FALSE}
results <- cinaR(..., save.DA.peaks = T, DA.peaks.path = "./Peaks_mice.xlsx")
```

will create an excel file with name `Peaks_mice.xlsx` in the current directory.

## Using different GLM algorithms
Currently, `cinaR` supports 4 different algorithms, namely; 

- edgeR
- limma-voom
- limma-trend
- DESeq2

If not set, it uses `edgeR` for differential analyses. 
You can change the used algorithm by simply setting `DA.choice` argument.
For more information, `?cinaR`

## Some useful arguments
```{r eval=FALSE}
# new FDR threshold for DA peaks
results <- cinaR(..., DA.fdr.threshold = 0.1)

# filters out pathways
results <- cinaR(..., enrichment.FDR.cutoff = 0.1)

# does not run enrichment pipeline
results <- cinaR(..., run.enrichment = FALSE)

# creates the piechart from chIpSeeker package
results <- cinaR(..., show.annotation.pie = TRUE)

# change cut-off value for dot plots
dot_plot(..., fdr.cutoff = 0.05)
```


## References

- Robinson MD, McCarthy DJ, Smyth GK (2010). “edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.” Bioinformatics, 26(1), 139-140. doi: 10.1093/bioinformatics/btp616.

- Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK (2015). “limma powers differential expression analyses for RNA-sequencing and microarray studies.” Nucleic Acids Research, 43(7), e47.

- Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. 10.1186/s13059-014-0550-8


## Session Info
```{r session Info}
sessionInfo()
```

