benchmarks.R

#!/usr/bin/env Rscript

# source('http://bioconductor.org/biocLite.R')
# bioClite('recount')
# bioClite('GenomicRanges')
# bioCLite('LieberInstitute/recount.bwtool')

library('GenomicRanges')
library('recount')
library('recount.bwtool')

gtex_study <- 'SRP000599'  #'SRP012682'
geuv_study <- 'SRP000542'  #'ERP001942'

do_gtex <- T
cpus <- 8

tab_to_granges <- function(tab) {
    GRanges(tab$chr, IRanges(tab$begin, tab$end), strand='+')
}

gc25 <- tab_to_granges(read.table('gencode_25.txt', sep=' ', header=T))
gc <- read.table('gencode.txt', sep=' ', header=F)
colnames(gc) <- c('xscript_id', 'gene_id', 'chr', 'begin', 'end')
gc <- tab_to_granges(gc)

experiment <- function(study, query, query_name, keep_sums=F) {
    sumsdir <- paste(study, query_name, sep='_')
    t <- proc.time()
    mat <- coverage_matrix_bwtool(study, query, verbose=T, sumsdir=sumsdir, bpparam=MulticoreParam(cpus))
    t <- proc.time() - t
    if(!keep_sums) {
        unlink(sumsdir, recursive=T)
    }
    list(time=t, mat=mat)
}

geuv_gc25 <- experiment(geuv_study, gc25, 'gencode25')
geuv_gc   <- experiment(geuv_study,   gc, 'gencode')
gtex_gc25 <- NULL
gtex_gc   <- NULL
if(do_gtex) {
    gtex_gc25 <- experiment(qtex_study, gc25, 'gencode25')
    gtex_gc   <- experiment(qtex_study,   gc, 'gencode')
}

save(cpus, geuv_gc, geuv_gc25, gtex_gc, gtex_gc25, file = "all.RData")