step8-coverageToExon.R

## Get exon coverage table for UCSC or Ensembl annotation

## Load libraries
library('getopt')

## Available at http://www.bioconductor.org/packages/release/bioc/html/derfinder.html
library('derfinder')

## Specify parameters
spec <- matrix(c(
    'experiment', 'e', 1, 'character', 'Experiment. Either stem, brainspan, snyder, hippo, or simulation',
    'annotation', 'a', 1, 'character', 'Annotation to use. Either ensembl or ucsc',
    'readlen', 'r', 1, 'integer', 'Read length',
    'mc.cores', 'c', 1, 'integer', 'Number of cores to use',
	'help' , 'h', 0, 'logical', 'Display help'
), byrow=TRUE, ncol=5)
opt <- getopt(spec)


## if help was asked for print a friendly message
## and exit with a non-zero error code
if (!is.null(opt$help)) {
	cat(getopt(spec, usage=TRUE))
	q(status=1)
}

## Check experiment input
stopifnot(opt$experiment %in% c('stem', 'brainspan', 'snyder', 'hippo', 'simulation'))



## Setup
message(paste(Sys.time(), 'loading coverage data'))
load("../CoverageInfo/fullCov.Rdata")

message(paste(Sys.time(), 'loading annotation data'))
if(opt$annotation == 'ensembl') {
    load('/home/epi/ajaffe/Lieber/Projects/RNAseq/derannotator/rdas/GenomicState.Hsapiens.ensembl.GRCh37.p11.rda')
    anno <- GenomicState.Hsapiens.ensembl.GRCh37.p11$fullGenome

} else if (opt$annotation == 'ucsc') {
    load('/home/epi/ajaffe/Lieber/Projects/RNAseq/derannotator/rdas/GenomicState.Hsapiens.UCSC.hg19.knownGene.rda')
    anno <- GenomicState.Hsapiens.UCSC.hg19.knownGene$fullGenome
}


## get table
message(paste(Sys.time(), 'running coverageToExon'))
covToEx <- coverageToExon(fullCov, anno, L=opt$readlen, strandCores = 1, mc.cores = opt$mc.cores)

## Save results
message(paste(Sys.time(), paste0("saving covToEx-", opt$annotation, ".Rdata")))
save(covToEx, file=paste0("covToEx-", opt$annotation, ".Rdata"))

## Done
proc.time()
sessionInfo()