index.Rmd

---
title: "derfinder tutorial"
author: "L. Collado-Torres"
date: "December 9th, 2014"
output: ioslides_presentation
---


```{r vignetteSetup, echo=FALSE, message=FALSE, warning = FALSE}
## Track time spent on making the vignette
startTime <- Sys.time()

## Bib setup
library('knitcitations')

## Load knitcitations with a clean bibliography
cleanbib()
cite_options(hyperlink = 'to.doc', citation_format = 'text', style = 'html')
# Note links won't show for now due to the following issue
# https://github.com/cboettig/knitcitations/issues/63

## Write bibliography information
bibs <- c(knitcitations = citation('knitcitations'),
    derfinder = citation('derfinder')[1], 
    knitrBootstrap = citation('knitrBootstrap'), 
    knitr = citation('knitr')[3],
    rmarkdown = citation('rmarkdown'),
    brainspan = RefManageR::BibEntry(bibtype = 'Unpublished', key = 'brainspan', title = 'Atlas of the Developing Human Brain [Internet]. Funded by ARRA Awards 1RC2MH089921-01, 1RC2MH090047-01, and 1RC2MH089929-01.', author = 'BrainSpan', year = 2011, url = 'http://developinghumanbrain.org'),
    originalder = citation('derfinder')[2],
    R = citation(),
    IRanges = citation('IRanges'),
    devtools = citation('devtools'),
    testthat = citation('testthat'),
    GenomeInfoDb = citation('GenomeInfoDb'),
    GenomicRanges = citation('GenomicRanges'),
    ggplot2 = citation('ggplot2'),
    biovizBase = citation('biovizBase'),
    bumphunter = citation('bumphunter'),
    TxDb.Hsapiens.UCSC.hg19.knownGene = citation('TxDb.Hsapiens.UCSC.hg19.knownGene'),
    AnnotationDbi = citation('AnnotationDbi'),
    BiocParallel = citation('BiocParallel'),
    derfinderHelper = citation('derfinderHelper'),
    GenomicAlignments = citation('GenomicAlignments'),
    GenomicFeatures = citation('GenomicFeatures'),
    GenomicFiles = citation('GenomicFiles'),
    Hmisc = citation('Hmisc'),
    qvalue = citation('qvalue'),
    Rsamtools = citation('Rsamtools'),
    rtracklayer = citation('rtracklayer'),
    S4Vectors = citation('S4Vectors'),
    bumphunterPaper = RefManageR::BibEntry(bibtype = 'article', key = 'bumphunterPaper', title = 'Bump hunting to identify differentially methylated regions in epigenetic epidemiology studies', author = 'Jaffe, Andrew E and Murakami, Peter and Lee, Hwajin and Leek, Jeffrey T and Fallin, M Daniele and Feinberg, Andrew P and Irizarry, Rafael A', year = 2012, journal = 'International Journal of Epidemiology'),
    derfinderData = citation('derfinderData'),
    polyester = citation('polyester')
)

write.bibtex(bibs,
    file = 'derTutorRef.bib')
bib <- read.bibtex('derTutorRef.bib')

## Assign short names
names(bib) <- names(bibs)

## Working on Windows?
windowsFlag <- .Platform$OS.type == 'windows'

evalDER <- TRUE
```

##

* Presentation: [lcolladotor.github.io/derTutor/](http://lcolladotor.github.io/derTutor/)
* R code: [lcolladotor.github.io/derTutor/index.R](http://lcolladotor.github.io/derTutor/index.R)

## Pre-reqs

This is a short tutorial on how to use `derfinder` `r citep(bib[['derfinder']])`.

In order to run this tutorial, you will need R 3.1.2 or newer (available from [CRAN](http://cran.r-project.org/)) and `derfinder` which is available from [Bioconductor](http://www.bioconductor.org/). It can be installed by:

```{r 'install', eval = FALSE}
## Install derfinder from BioC
source('http://bioconductor.org/biocLite.R')
biocLite('derfinder')
```

You will also need [TxDb.Hsapiens.UCSC.hg19.knownGene](http://bioconductor.org/packages/release/data/annotation/html/TxDb.Hsapiens.UCSC.hg19.knownGene.html)

```{r 'install2', eval = FALSE}
biocLite('TxDb.Hsapiens.UCSC.hg19.knownGene')
```



## Docs

You can find the latest documentation at

* [Intro to derfinder](http://www.bioconductor.org/packages/devel/bioc/vignettes/derfinder/inst/doc/derfinder.html) vignette
* [Advanced details and usage](http://www.bioconductor.org/packages/devel/bioc/vignettes/derfinder/inst/doc/derfinderAdvanced.html) vignette
* [derfinderPlot](http://www.bioconductor.org/packages/devel/bioc/html/derfinderPlot.html)
* [regionReport](http://www.bioconductor.org/packages/devel/bioc/html/regionReport.html)
* help pages

```{r 'rhelp', eval = FALSE}
help(package = 'derfinder')
```


Coming soon:

* `derfinder` paper


## Why should I use `derfinder`?

> * You want to perform differential expression analysis at base-pair resolution
> * Because it's annotation-agnostic
> * Can handle your data set (has been used with ~500 samples)
> * Others have gotten grants partially thanks to it

## {.flexbox .vcenter}

<img alt = "tim's grant" height='560' src='images/tim-grant.PNG' />


## What will I need?

* Most likely a high performance computing environment. It'll depend on the size of your data set
* Familiarity with other packages such as [GenomicRanges](http://www.bioconductor.org/packages/release/bioc/html/GenomicRanges.html)
* Experience writing `bash` scripts can help, although its not required


## What does it do?

* It calculates F-statistics at every base: uses nested models
* Let $y_{ij} = \log_2\left( \text{coverage}_{ij} + \text{scaling factor} \right)$ for position $i$ and sample $j$
* Alternative model:
$$ y_{ij} =  \alpha_i + \sum_{p=1}^n \beta_{ip} X_{jp}  + \sum_{q=1}^m \gamma_{iq} Z_{jq} + \epsilon_{ij} $$
* Null model doesn't have $\beta$ terms.
* F-statistics are given by:
$$ F_i = \frac{(\text{RSS}0_i - \text{RSS}1_i) / (\text{df}_1 - \text{df}_0) }{ \text{offset} + (\text{RSS}1_i / (n - \text{df}_1)) } $$


## {.flexbox .vcenter}

<img alt = 'fstats' height='560' src='images/fig1.png' />

## Main functions {.flexbox .vcenter}

<img alt = 'flow' height='530' src='http://lcolladotor.github.io/derfinder/fig/DERpathway.png' />


## Finding DERs {.flexbox .vcenter}

<img alt = 'analyzechr chart' height='530' src='http://lcolladotor.github.io/derfinder/fig/analyzeChr.png' />

## Data

Simulated data using [polyester](http://www.bioconductor.org/packages/release/bioc/html/polyester.html) `r citep(bib[['polyester']])` for 3 groups, 10 samples per group, fold changes 2x and $\frac{1}{2}$x.

* 24 single transcript genes
    * 12 set to be DE: 4 per group (2 high, 2 low, other groups normal)
    * 12 not DE
* 36 two transcript genes
    * 12 not DE
    * 12 one transcript DE
    * 12 both transcripts DE

Aligned with TopHat 2.0.13. Saved coverage in BigWig files.

## Load data


```{r 'loadDer', message = FALSE, eval = evalDER}
library('derfinder')
library('TxDb.Hsapiens.UCSC.hg19.knownGene')

## Files
files <- paste0('http://lcolladotor.github.io/derTutor/data/sample', 
    1:30, '.bw')
names(files) <- paste0('sample', 1:30)

## Load data
system.time( fullCov <- fullCoverage(files, 'chr22', verbose = FALSE) )

## You are using Windows?
# Use http://lcolladotor.github.io/derTutor/data/fullCov.Rdata
```

## Lib size {.smaller}


```{r 'adjustments', eval = evalDER}
## Calculate library size adjustment
system.time( collapsedFull <- collapseFullCoverage(fullCov) )
lapply(collapsedFull[[1]], head)
sampleDepths <- sampleDepth(collapsedFull, probs = 1)

```


## Build models {.smaller}


```{r 'models', eval = evalDER}
## Create models
groupInfo <- factor(rep(c('A', 'B', 'C'), each = 10))
models <- makeModels(sampleDepths = sampleDepths, testvars = groupInfo)
```

## Run analysis {.smaller}


```{r 'analysis', eval = evalDER}
## Filter Data
covData <- filterData(fullCov$chr22, cutoff = 0)

## Run analysis for chr22
dir.create('chr22', showWarnings = FALSE)
system.time(
    res <- analyzeChr(chr = 'chr22', coverageInfo = covData, models = models, 
        cutoffFstat = 1e-03, cutoffPre = 0,
        nPermute = 100, seeds = seq_len(100) + 20141202, maxClusterGap = 3000,
        groupInfo = groupInfo, mc.cores = 1,
        lowMemDir = file.path(tempdir(), 'chr22', 'chunksDir'),
        writeOutput = TRUE, returnOutput = TRUE)
)
```

## Explore results {.smaller}

```{r 'exploreRes', eval = evalDER}
names(res)
dir('chr22')
names(res$coveragePrep)
names(res$regions)
```


## Merge {.smaller}

```{r 'merge', eval = evalDER}
## Genomic state
system.time(gs <- makeGenomicState(txdb=TxDb.Hsapiens.UCSC.hg19.knownGene, chrs='chr22'))

## Merge results from different chrs
mergeResults('chr22', genomicState = gs$fullGenome, optionsStats = res$optionsStats)
```

## Main results {.smaller}

```{r 'final', eval = evalDER}
dir(pattern = 'Rdata')
load('fullRegions.Rdata')
class(fullRegions)
length(fullRegions)
```

##  {.smaller}

```{r 'final2', eval = evalDER}
colnames(mcols(fullRegions))
table(fullRegions$significantFWER)
```


## Need help?

* Post in the [Bioconductor support site](https://support.bioconductor.org/) -- remember to use the `derfinder` tag!
* If you think you found a bug, please report it on the [derfinder GitHub repo](https://github.com/lcolladotor/derfinder/issues)


## Citing `derfinder`

```{r 'citation'}
## Citation info
citation('derfinder')
```

## Reproducibility {.smaller}

Code for creating this page

```{r createVignette, eval=FALSE}
## Create this page
library('rmarkdown')
render('index.Rmd')

## Clean up
file.remove('derTutorRef.bib')

## Extract the R code
library('knitr')
knit('index.Rmd', tangle = TRUE)
```

Date this tutorial was generated.

```{r reproducibility1, echo=FALSE}
## Date the vignette was generated
Sys.time()
```

Wallclock time spent running this tutorial.

```{r reproducibility2, echo=FALSE}
## Processing time in seconds
totalTime <- diff(c(startTime, Sys.time()))
round(totalTime, digits=3)
```

## {.smaller}

`R` session information.

```{r reproducibility3, echo=FALSE}
## Session info
library('devtools')
options(width = 120)
session_info()$platform
```

## {.smaller}

```{r reproducibility4, echo=FALSE}
## Session info packages
subset(session_info()$packages, package %in% c('derfinder', unique(names(getNamespaceImports('derfinder'))),  'TxDb.Hsapiens.UCSC.hg19.knownGene', 'rmarkdown', 'knitcitations', 'Matrix'))
```

## Bibliography {.smaller}

This tutorial was generated using `rmarkdown` `r citep(bib[['rmarkdown']])` and `knitcitations` `r citep(bib[['knitcitations']])`.

```{r vignetteBiblio, results = 'asis', echo = FALSE, warning = FALSE}
## Print bibliography
bibliography()
```