---
title: "Timing information"
output:
html_document:
toc: true
toc_float: true
code_folding: hide
---
```{r citationsSetup, echo=FALSE, message=FALSE, warning=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"),
"GenomicRanges" = citation("GenomicRanges"),
"DESeq" = citation("DESeq"),
"DT" = citation("DT"),
"ggplot2" = citation("ggplot2"),
'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'),
'rail' = RefManageR::BibEntry('article', key = 'rail', author = 'Abhinav Nellore and Leonardo Collado-Torres and Andrew E. Jaffe and José Alquicira-Hernández and Jacob Pritt and James Morton and Jeffrey T. Leek and Ben Langmead', journal = 'bioRxiv', year = '2015', title = 'Rail-RNA: {Scalable} analysis of {RNA}-seq splicing and coverage'))
write.bibtex(bibs, file = 'timing.bib')
bib <- read.bibtex('timing.bib')
## Assign short names
names(bib) <- names(bibs)
```
This report shows the time and memory used to run `derfinder` `r citep(bib[["derfinder"]])` for single base resolution differential expression analysis. It also shows the same information for going from BAM files to getting ready to run `DESeq` `r citep(bib[["DESeq"]])` by using `samtools` `r citep("http://samtools.sourceforge.net/")` to convert to SAM format and `HTSeq` `r citep("http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html")` to make the count tables. Furthermore, this process was compared to using the `summarizeOverlaps()` function from the `GenomicRanges` `r citep(bib[["GenomicRanges"]])` package as well as using the `coverageToExon()` function included in the `derfinder` package [requires the output from the _fullCov_ step].
# Results
```{r 'effanalytics', eval = FALSE, bootstrap.show.code = FALSE, boostrap.show.output = FALSE}
## Extract information from Gmail
system('cp ../../efficiency_analytics/client_secrets .')
system('python ../../efficiency_analytics/analyze_efficiency.py --email fellgernon@gmail.com --folder "Cluster/derSoftware" --outfile timing-derSoftware.txt')
```
```{r loadLibs, warning = FALSE}
## Load libraries
library("ggplot2")
library("knitr")
```
```{r process}
## Setup
## Define number of cores used
exps <- c('brainspan', 'simulation', 'hippo', 'snyder', 'stem')
## Read data and process it
all <- read.table('timing-derSoftware.txt', header = TRUE, stringsAsFactors = FALSE)
all <- all[!grepl('brainspan.*run3', all$jobid), ] # remove older info
all$step <- gsub('.*th', 'TopHat', sapply(strsplit(all$jobid, "-"), function(x) x[1]))
all$memG <- all$memory
all$memG[all$memunit == "M"] <- all$memG[all$memunit == "M"] / 1024
all$chr <- gsub('.*chr', 'chr', all$jobid)
all$chr[ !grepl('chr', all$chr) ] <- NA
## Experiment info
all$experiment <- NA
for(exp in exps) {
all$experiment[ grepl(exp, tolower(all$jobid)) ] <- exp
}
all$experiment[ all$step %in% c('TopHat', 'bigwig') ] <- 'simulation'
all$experiment[ all$jobid == 'makeBai-Sim' ] <- 'simulation'
## Cores info
all$cores <- mapply(function(chr, exp, step) {
if(step == 'fullCov') {
return(10L)
} else if(step == 'derA') {
if(exp == 'brainspan') {
return(ifelse(chr == 'chrY', 2L, ifelse(chr == 'chr1', 40L, ifelse(chr == 'chr2', 32L, ifelse(chr == 'chr3', 27L, ifelse(chr == 'chr19', 29L, 20L))))))
} else if (exp == 'simulation'){
return(1L)
} else if (exp == 'hippo'){
return(2L)
} else if (exp == 'snyder'){
return(4L)
} else if (exp == 'stem'){
return(8L)
}
} else if(step == 'regMat') {
return(5L)
} else if(step == 'TopHat') {
return(4L)
} else if(step == 'summOv') {
return(ifelse(exp == 'hippo', 24L, 10L))
} else {
return(1L)
}
}, all$chr, all$experiment, all$step)
all$timeByCore <- all$walltime * all$cores
all$memByCore <- all$memG / all$cores
## Add software labels
all$software <- factor(ifelse(all$step %in% c('toSam', 'htseq'), 'HTSeq', ifelse(all$step == 'summOv', 'GenomicRanges', ifelse(all$step == 'TopHat', 'TopHat', ifelse(all$step %in% c('makeBai', 'regVsDERs', 'PNAS', 'summInfo'), 'misc', ifelse(all$step == 'derR', 'regionReport', 'derfinder'))))))
## Experiment and cores groups info
all$experiment <- factor(all$experiment, levels = exps)
all$coresGroups <- all$cores
all$coresGroups[ all$cores >= 20] <- '20+'
all$coresGroups <- factor(all$coresGroups, levels = c(1, 2, 4, 5, 8, 10, '20+'))
## Types of analysis
all$analysis <- factor(ifelse(all$step %in% c('derMod', 'derA', 'derM'), 'Single-base DER', ifelse(all$step %in% c('toSam', 'htseq', 'summOv', 'covToEx'), 'Exon count', ifelse(all$step == 'regMat', 'Expressed-region DER', ifelse(all$step == 'fullCov', 'Load data', ifelse(all$step == 'derR', 'HTML report', 'misc'))))))
## Show only information for the data sets described in this website
all <- subset(all, experiment %in% c('hippo', 'snyder'))
```
## Adjusting by number of cores
The following plots show the wall time and memory used by each job while taking into account the number of cores used by each job. Note that doing so is a crude approximation of how much time and memory each job would have needed had it ran on a single node.
Points are colored by which analysis type they belong to. Note that the loading data step is required for the single-level and expressed-regions DER approaches as well as exon counting (with derfinder).
```{r edaAnalysis, fig.width=10, fig.height=7}
## Walltime and memory adjusted by number of cores (it's an approximation)
ggplot(all, aes(x=timeByCore, y=memByCore, colour=analysis, shape=software)) + geom_point(size = 3) + facet_grid(~ experiment) + xlab("Wall time (hrs) multiplied by the number of cores") + ylab("Memory (GB) divided by the number of cores") + scale_colour_brewer(palette="Dark2") + theme_bw(base_size = 18) + theme(legend.position=c(.5, .75), legend.box = 'horizontal')
ggplot(all, aes(x=log2(timeByCore), y=memByCore, colour=analysis, shape=software)) + geom_point(size = 3) + facet_grid(~ experiment) + xlab("Wall time (hrs) multiplied by the number of cores (log2)") + ylab("Memory (GB) divided by the number of cores") + scale_colour_brewer(palette="Dark2") + theme_bw(base_size = 18) + theme(legend.position=c(.5, .75), legend.box = 'horizontal')
## For supp text
time <- ggplot(subset(all, !software %in% c('TopHat', 'regionReport') & analysis != 'misc'), aes(x=log2(timeByCore), y=log2(memByCore), colour=analysis, shape=software)) + geom_point(size = 3) + facet_grid(~ experiment) + xlab("Wall time (hrs) multiplied by the number of cores (log2)") + ylab("GB memory divided by number of cores (log2)") + scale_colour_brewer(palette="Set1") + theme_bw(base_size = 18) + theme(legend.position=c(.55, .15), legend.box = 'horizontal')
time
pdf(file = 'time.pdf', width = 10)
time
dev.off()
#system('open time.pdf')
```
## Resources by step for each analysis
```{r 'analysisSummary'}
getInfo <- function(df, sumTime = FALSE, peakCores = FALSE) {
memByCore <- max(df$memByCore)
walltime <- ifelse(sumTime, sum(df$walltime), max(df$walltime))
memG <- max(df$memG)
peakCores <- ifelse(peakCores, max(df$peakCores), sum(df$cores))
res <- c(memByCore = memByCore, walltime = walltime, memG = memG, peakCores = peakCores)
return(res)
}
analysisInfo <- list('Single-base DER' = c('Load data', 'Single-base DER'),
'Expressed-region DER' = c('Load data', 'Expressed-region DER'),
'HTML report' = 'HTML report',
'Exon count - derfinder' = 'Load data'
)
analysisInfo <- lapply(analysisInfo, function(x) { which(all$analysis %in% x)})
analysisInfo[[4]] <- c(analysisInfo[[4]], which(all$step == 'covToEx'))
analysisInfo$"Exon count - HTSeq" <- which(all$step %in% c('toSam', 'htseq'))
analysisInfo$"Exon count - GenomicRanges" <- which(all$step == 'summOv')
## Summarize the information for each step of each analysis
analysisSummary <- lapply(names(analysisInfo), function(analysis) {
current <- all[analysisInfo[[analysis]], ]
res_analysis <- lapply(exps, function(exp) {
use <- subset(current, experiment == exp)
if(nrow(use) == 0) return(NULL)
res_exp <- lapply(unique(use$step), function(step) {
res_step <- as.data.frame(t(getInfo(use[use$step == step, ])))
res_step$step <- step
res_step$experiment <- exp
res_step$analysis <- analysis
return(res_step)
})
res_exp <- do.call(rbind, res_exp)
return(res_exp)
})
res_analysis <- do.call(rbind, res_analysis)
return(res_analysis)
})
analysisSummary <- do.call(rbind, analysisSummary)
```
The table shown below shows per analysis the maximum memory used by a job and maximum wall time for that step. This is assuming that all jobs for a given step ran simultaneously. For example, that all jobs running `derfinder::analyzeChr()` were running at the same time. Note that for some analyses relied on the same steps, like loading the data (_fullCov_). This table can be useful to find the peak number of cores (the sum of cores for all jobs running simultaneously) for a given analysis step.
```{r 'analysisSumTab', results = 'asis'}
kable(analysisSummary, format = 'markdown', digits = c(2, 4, 2))
```
## Resources for each analysis
```{r 'peakSummary'}
## Summary the information for each analysis
peaks <- lapply(names(analysisInfo), function(analysis) {
res_analysis <- lapply(exps, function(exp) {
current <- analysisSummary[analysisSummary$analysis == analysis & analysisSummary$experiment == exp, ]
if(nrow(current) == 0) return(NULL)
res_exp <- as.data.frame(t(getInfo(current, sumTime = TRUE, peakCores = TRUE)))
res_exp$experiment <- exp
res_exp$analysis <- analysis
return(res_exp)
})
res_analysis <- do.call(rbind, res_analysis)
return(res_analysis)
})
peaks <- do.call(rbind, peaks)
save(peaks, file = 'peaks.Rdata')
```
We can further summarize the resources used by each analysis by identified the maximum memory used in the steps required for a particular analysis and the total wall time for running all the steps when all the jobs of a particular step are running simultaneously. Thus giving us the total actual wall time to run a specific analysis and the maximum memory required.
The table below shows the final summary. Note that in some analyses, the peak memory is from the _fullCov_ step. We did not focus on reducing the memory load of this step as we sacrificed memory for speed. We know that much lower memory limits can be achieved using 1 core instead of the 10 cores used.
```{r 'peakSumTab', results = 'asis'}
kable(peaks, format = 'markdown', digits = c(2, 3, 2))
```
Regarding the high memory load for the HTML report, this could be significantly lowered by only loading the required coverage data used for the plots instead of the full output from the _fullCov_ step. That is, using the _which_ argument from `fullCoverage()` to create a much smaller _fullCov_ object, which would also reduce the memory used when plotting.
__Note__: since these analyses were done, we have found other ways to run `derfinder::regionMatrix()` that require less memory. In particular, if you have BigWig files (as those generated by `Rail-RNA` `r citep(bib[['rail']])`), we recommend using `railMatrix()`.
## Comparing methods for gene count table generation
The previous table can also be used to compare the sum of the time and peak memory used by the different steps to obtain the exon count table with the following software options.
* `derfinder`: includes resources used for reading coverage data in `R` and then running creating a feature count matrix. We did so for
* UCSC hg19 knownGene annotation
* Ensembl GRCh37 p11 annotation.
* `HTSeq`: includes resources used for generating sorted SAM files and then running HTSeq.
* `summOv`: resources used for running `GenomicRanges::summarizeOverlaps()` directly on the BAM files.
# Details
The following table shows the details of the resources used by the different jobs. It shows the experiment (_experiment_), the analysis step (_step_), wall time used (shown in hours, _walltime_), number of cores used (_cores_), memory in GB used (_memG_), software used (_software_), analysis for which the step is used (_analysis_), and the job name (_jobib_). Furthermore, it shows two simple approximations:
* _timeByCore_ is the wall time (in hours) multiplied by the number of cores used. It is a very simple approximation for the wall time used had the job been ran on a single node. This approximation is known to be false, but it gives a basic idea.
* _memByCore_ is the memory (in GB) divided by the number of cores used. It is an approximation for the memory used had the job been ran on a single node.
These are the following analysis steps:
1. __fullCov__ Extract coverage information from raw files (BAM or BigWig) by chromosome, then filter it, and save it in Rdata files.
1. __derMod__ Calculate the sample depth adjustments and build models appropriate for the experiment.
1. __derA__ Run single base-level analysis by chromosome.
1. __derM__ Merge derfinder analysis results from the different chromosomes, calculate p-values and q-values.
1. __derR__ Generate HTML report with `regionReport`.
1. __regMat__ Run expressed regions-level analysis with `regionMatrix()`.
1. __regsVsDers__ Compare expressed regions-level vs single base-level approaches.
1. __toSam__ Transform BAM files to sorted (by name) SAM files for running HTSeq.
1. __htseq__ Run HTSeq to generate exon count table.
1. __summOv__ Run `GenomicRanges::summarizeOverlaps()` to generate exon count table.
1. __covToExon__ Generate exon table using `derfinder::coverageToExon()` for UCSC hg19 knownGene or GRCh37 p11 Ensembl annotation table.
1. __PNAS__ (Only for _Hippo_) Generate an HTML report comparing the derfinder results vs previously published results (PNAS paper).
1. __summInfo__ Summarize results to then use then in the derfinder software paper.
```{r tables, results="asis"}
library("DT")
## Print whole table
d <- all[, c("experiment", "step", "walltime", "cores", "memG", "timeByCore", "memByCore", "software", "analysis", "jobid")]
datatable(d, options = list(pagingType='full_numbers', pageLength=50, scrollX='100%')) %>% formatRound(columns = c(3, 5:7), digits = 3)
```
Table made using `DT` `r citep(bib[["DT"]])`.
# Reproducibility
Date the report was generated.
```{r reproducibility1, echo=FALSE}
## Date the report was generated
Sys.time()
```
Wallclock time spent generating the report.
```{r "reproducibility2", echo=FALSE}
## Processing time in seconds
totalTime <- diff(c(startTime, Sys.time()))
round(totalTime, digits=3)
```
`R` session information.
```{r "reproducibility3", echo=FALSE}
## Session info
options(width=120)
devtools::session_info()
```
# Bibliography
This report was generated using `rmarkdown` `r citep(bib[['rmarkdown']])` with `knitr` `r citep(bib[['knitr']])` running behind the scenes. Timing information extracted from the SGE reports using `efficiency analytics` `r citep(bib[["eff"]])`. Figures and citations were made using `ggplot2` `r citep(bib[["ggplot2"]])` and `knitcitations` `r citep(bib[['knitcitations']])` respectively.
Citation file: [timing.bib](timing.bib)
```{r vignetteBiblio, results = 'asis', echo = FALSE, warning = FALSE}
## Print bibliography
bibliography()
```