---
output:
  knitrBootstrap::bootstrap_document:
    theme.chooser: TRUE
    highlight.chooser: TRUE
---

Supplementary Website
===================================

```{r citationsSetup, echo=FALSE, message=FALSE, warning=FALSE, bootstrap.show.code=FALSE}
## Track time spent on making the report
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

bibs <- c("knitcitations" = citation("knitcitations"),
    "derfinder" = citation("derfinder")[1],
    "regionReport" = citation("regionReport")[1],
    "GenomicRanges" = citation("GenomicRanges"),
    "DESeq" = citation("DESeq"),
    "knitrBootstrap" = citation("knitrBootstrap"),
    'rmarkdown' = citation('rmarkdown'),
    'knitr' = citation('knitr')[3],
    'eff' = RefManageR::BibEntry('manual', key = 'eff', title = 'Efficiency analysis of Sun Grid Engine batch jobs', author = 'Alyssa Frazee', year = 2014, url = 'http://dx.doi.org/10.6084/m9.figshare.878000'),
    'zhou2011' = RefManageR::BibEntry('article', key = 'zhou2011', author = "Zhou, Zhifeng and Yuan, Qiaoping and Mash, Deborah C and Goldman, David", title = "Substance-specific and shared transcription and epigenetic changes in the human hippocampus chronically exposed to cocaine and alcohol", journal = "Proceedings of the National Academy of Sciences of the United States of America", year = 2011, volume = "108", number = "16", pages = "6626-6631"))
    
write.bibtex(bibs, file = 'index.bib')
bib <- read.bibtex('index.bib')

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

This page describes the supplementary material for the `derfinder` software paper which includes several HTML reports as well as code files for reproducing the results.

# Code to reproduce analyses


All the `bash`, `R` and `R Markdown` source files used to analyze the data for this project as well as generate the HTML reports are available in this website. However, it is easier to view them at [github.com/leekgroup/derSoftware](https://github.com/leekgroup/derSoftware).

There are 9 main `bash` scripts named _step1-*_ through _step9-*_. 

1. _fullCoverage_ loads the data from the raw files. See [step1-fullCoverage.sh](step1-fullCoverage.sh) and [step1-fullCoverage.R](step1-fullCoverage.R).
1. _makeModels_ creates the models used for the single-level base analysis.  See [step2-makeModels.sh](step2-makeModels.sh) and [step2-makeModels.R](step2-makeModels.R).
1. _analyzeChr_ runs the single-base level analysis by chromosome.  See [step3-analyzeChr.sh](step3-analyzeChr.sh) and [step3-analyzeChr.R](step3-analyzeChr.R).
1. _mergeResults_ merges the single-base level analysis results for all the chromosomes. See [step4-mergeResults.sh](step4-mergeResults.sh).
1. _derfinderReport_ generates a HTML report for the single-base level DERs. See [step5-derfinderReport.sh](step5-derfinderReport.sh).
1. _regionMatrix_ identifies the expressed regions for the expressed-regions level approach. See [step6-regionMatrix.sh](step6-regionMatrix.sh).
1. _regMatVsDERs_ creates a simple HTML report comparing the single-base and expressed-regions approaches. See [step7-regMatVsDERs.sh](step7-regMatVsDERs.sh) and [step7-regMatVsDERs.Rmd](step7-regMatVsDERs.Rmd).
1. _coverageToExon_ creates an exon count table using known annotation information. See [step8-coverageToExon.sh](step8-coverageToExon.sh) and [step8-coverageToExon.R](step8-coverageToExon.R).
1. _summaryInfo_ creates a HTML report with brief summary information for the given experiment. See [step9-summaryInfo.sh](step9-summaryInfo.sh), [step9-summaryInfo.R](step9-summaryInfo.R), and [step9-summaryInfo.Rmd](step9-summaryInfo.Rmd).

There are also 3 optional `bash` scripts used when BAM files are available.

1. _sortSam_ creates sorted by sequence name SAM files. See [optional1-sortSam.sh](optional1-sortSam.sh).
1. _HTSeq_ creates the exon count tables using `HTSeq`. See [optional2-HTSeq.sh](optional2-HTSeq.sh).
1. _summOv_ uses `GenomicRanges` to create the exon count tables. See [optional3-summOv.sh](optional3-summOv.sh) and [optional3-summOv.R](optional3-summOv.R).

A final `bash` script, [run-all.sh](run-all.sh), can be used to run the main 9 steps (or a subset of them).

All 13 `bash` scripts show at the top the way they were used. Some of them generate intermediate small `bash` scripts, for example one script per chromosome for the _analyzeChr_ step. For some steps, there is a companion `R` or `R Markdown` code file when the code is more involved or an HTML file is generated in the particular step.



The [check-analysis-time.R](check-analysis-time.R) script was useful for checking the progress of the _step3-analyzeChr_ jobs and detect whenever a node in the cluster was presenting problems.



We expect that these scripts will be useful to `derfinder` users who want to automate the single-base level and/or expressed-regions level analyses for several data sets and/or have the jobs run automatically without having to check if each step has finished running.


Note that all `bash` scripts are tailored for the cluster we have access to which administer job queues with Sun Grid Engine (SGE).


# Single-base level

## Quick overview HTML reports

These HTML reports contain  basic information on the `derfinder` `r citep(bib[["derfinder"]])` results from the public data sets used (_BrainSpan_, _Simulation_, and _Hippo_). These reports answer basic questions such as:

* What is the number of filtered bases?
* What is the number of candidate regions?
* How many candidate regions are significant?

They also illustrate three clusters of candidate differentially expressed regions (DERs) from the single-base level analysis. You can view the reports by following these links:

* [BrainSpan](brainspan/summaryInfo/run4-v1.0.10/summaryInfo.html)
* [Simulation](simulation/summaryInfo/run2-v1.0.10/summaryInfo.html)
* [Hippo](hippo/summaryInfo/run3-v1.0.10/summaryInfo.html)

## CSV files and annotation comparison

This HTML report has the code for loading the R data files and generating the CSV files. The report also has Venn diagrams showing the number of candidate DERs from the single-base level analysis that overlap known exons, introns and intergenic regions using the UCSC hg19 annotation. It also includes a detailed description of the columns in the CSV files.

View the [venn](venn/venn.html) report or its `R Markdown` source file [venn.Rmd](venn/venn.Rmd). 

# Expressed-regions level


For each experiment, we made a simple comparison between the single-base level DERs and the expressed regions identified via `regionMatrix()`. Note that you would need to use `limma`, `edgeR`, `DESeq`, or another count-level differential expression package to determine which expressed regions are differentially expressed. In the simulation case we did so using `limma` as is described in the corresponding section. For some data sets we used more than one mean cutoff (shown in parenthesis below) for determining the expressed regions.

These reports show how many bases are picked up in each of the approaches and different overlap comparisons. You can view the reports by following these links:


* [BrainSpan (cutoff 0.25)](brainspan/regionMatrix-vs-DERs/cut0.25-vs-run4-v1.0.10/step7-regMatVsDERs.html)
* [BrainSpan (cutoff 0.1)](brainspan/regionMatrix-vs-DERs/cut0.1-vs-run4-v1.0.10/step7-regMatVsDERs.html)
* [Simulation (cutoff 0)](simulation/regionMatrix-vs-DERs/cut0-vs-run2-v1.0.10/step7-regMatVsDERs.html)
* [Simulation (cutoff 5)](simulation/regionMatrix-vs-DERs/cut5-vs-run2-v1.0.10/step7-regMatVsDERs.html)
* [Hippo (cutoff 3)](hippo/regionMatrix-vs-DERs/cut3-vs-run3-v1.0.10/step7-regMatVsDERs.html)

# Simulation

## Generating reads

The code for generating the simulated RNA-seq reads and the chosen setup is described in the [generateReads](simulation/generateReads.html) report. This report is generated by the `R Markdown` [generateReads.Rmd](simulation/generateReads.Rmd) file.

The code for aligning the reads to the genome with TopHat is in the [run-paired-tophat.sh](simulation/thout/run-paired-tophat.sh) and [makeBai.sh](simulation/thout/makeBai.sh) scripts.

There is also code for exporting the coverage data to BigWig files, which was necessary for a tutorial on how to use `derfinder`: [derTutor](http://lcolladotor.github.io/derTutor/). The code is available in the [bigwig.sh](simulation/bw/bigwig.sh) script.

## Single-base level evaluation

A through evaluation of the simulation results from the single-base level analysis is described in the [evaluate](simulation/evaluate/evaluate.html) report. Several comparisons are made at the gene, transcript and exonic segment levels. 

The `R Markdown` source file for this report is [evaluate.Rmd](simulation/evaluate/evaluate.Rmd)

## Expressed-regions level evaluation

Similarly to the previous report, the expressed-regions level analysis evaluation is described in the [evaluate-regionMatrix](simulation/evaluate-regionMat/evaluate-regionMatrix.html) report. The code and language were slightly modified from the previous report. 

The `R Markdown` source file for this report is [evaluate-regionMatrix.Rmd](simulation/evaluate-regionMat/evaluate-regionMatrix.Rmd)

## Comparing versus other methods

We have several scripts and reports for the comparison between `derfinder` and  `DESeq2`, `edgeR`-robust as well as the original implementation of _derfinder_. 

This first set of scripts were used to run the original implementation of _derfinder_.

* [dercount.sh](simulation/dercount/dercount.sh) is the shell script for counting that is required for the original implementation of _derfinder_ analysis.
* [dertable.sh](simulation/dertable/dertable.sh) and [makeTableDB.R](simulation/dertable/makeTableDB.R) are the bash and R scripts needed for generating the SQL databases required for the original implementation of _derfinder_.
* [derfinder-analysis.R](simulation/derfinder-original/derfinder-analysis.R) and [deranalysis.sh](simulation/derfinder-original/deranalysis.sh) are the R and bash scripts for running the original implementation of _derfinder_ on the simulated data set.

The following scripts and reports show the comparison between these methods.

* [counts-based](simulation/counts-based/counts-based.html) is an initial comparison of the disjoint exon counts-based methods, the original implementation of _derfinder_, and `derfinder` itself.  The `R Markdown` source is [counts-based.Rmd](simulation/counts-based/counts-based.Rmd).
* [all-exons](simulation/all-exons/all-exons.html) determines differentially expressed disjoint exons for five exon sets. This was necessary to check the performance of different exon-level pipelines. The `R Markdown` source is [all-exons.Rmd](simulation/all-exons/all-exons.Rmd).
* [counts-gene](simulation/counts-gene/counts-gene.html) is where the different gene sets are defined: a complete one, an incomplete one with only one transcript per gene, and one where 20 percent of transcripts are missing.  The `R Markdown` source is [counts-gene.Rmd](simulation/counts-gene/counts-gene.Rmd).
* [counts-gene-eval](simulation/counts-gene-eval/counts-gene-eval.html) compares `derfinder` against the gene-level counts-based methods.  The `R Markdown` source is [counts-gene-eval.Rmd](simulation/counts-gene-eval/counts-gene-eval.Rmd).




# Timing information


This HTML report has code for reading and processing the time and memory information for each job extracted with [efficiency_analytics](https://github.com/alyssafrazee/efficiency_analytics) `r citep(bib[["eff"]])`. Several plots exploring the relationship between wall time and memory used by the cluster jobs are included. Some of the plots make explicit the number of cores used by each job. The report contains a detailed description of the analysis steps shown on the plots. It also contains tables summarizing the maximum memory and time for each analysis step if all the jobs for that particular step were running simultaneously. Finally, there is an interactive table with the timing results.

View the [timing](timing/timing.html) report or check the `R Markdown` file [timing.Rmd](timing/timing.Rmd).

 

# Miscellaneous

## Figure 1

The code used for generating the panels using in Figure 1 of the paper is available in the [figure1.R](figure1/figure1.R) file.


## Figure 2

The code used for generating the panels using in Figure 2 of the paper is available in the [figure2.R](figure2/figure2.R) file.

## Additional analyses

The following `R` source files have the code for reproducing additional analyses described in the paper

* [brainspan_regionLevel.R](additional-analyses/brainspan_regionLevel.R) Analysis of BrainSpan expressed-regions level DERs.
* [characterize_brainspan_DERs.R](additional-analyses/characterize_brainspan_DERs.R) Analysis of BrainSpan single-base level DERs.
* [feature_counts.R](additional-analyses/feature_counts.R) Feature counts analysis of Hippo and Snyder data sets.

These scripts also include other exploratory code.

## Hippo vs previous results

[compareVsPNAS](hippo/pnas/compareVsPNAS.html) is an HTML report comparing 29 regions that were previously found to be differentially expressed `r citep(bib[["zhou2011"]])` versus the `derfinder` single-base level results. It also has code for identified differentially expressed disjoint exons. The additional script [counts-gene.R](hippo/counts-gene/counts-gene.R) has the code for gene counting with `summarizeOverlaps()`. [compareVsPNAS-gene](hippo/pnas/compareVsPNAS-gene.html) compares the results between `DESeq2` and `edgeR`-robust against `derfinder` at the gene level with 40 total plots: 10 for each case of agreement/disagreement. 

View the [compareVsPNAS](hippo/pnas/compareVsPNAS.html) report or check the `R Markdown` file [compareVsPNAS.Rmd](hippo/pnas/compareVsPNAS.Rmd) run by the [runComparison.sh](hippo/pnas/runComparison.sh) script. Also view the [compareVsPNAS-gene](hippo/pnas/compareVsPNAS-gene.html) report and its linked `R Markdown` file [compareVsPNAS-gene.Rmd](hippo/pnas/compareVsPNAS-gene.Rmd).


# Reproducibility

Date this page was generated.

```{r reproducibility1, echo=FALSE, bootstrap.show.code=FALSE}
## Date the report was generated
Sys.time()
```

Wallclock time spent generating the report.

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

`R` session information.

```{r "reproducibility3", echo=FALSE, bootstrap.show.code=FALSE, bootstrap.show.message=FALSE}
## Session info
options(width=120)
devtools::session_info()
```

You can view the source `R Markdown` file for this page at [index.Rmd](index.Rmd).

# Bibliography

This report was generated using `knitrBootstrap` `r citep(bib[['knitrBootstrap']])`
with `knitr` `r citep(bib[['knitr']])` and `rmarkdown` `r citep(bib[['rmarkdown']])` running behind the scenes.  

Citations were made with  `knitcitations` `r citep(bib[['knitcitations']])`.

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


```{r clean, echo=FALSE, results="hide"}
## Clean up
file.remove("index.bib")
```

