step3-analyzeChr.R

## Run derfinder's analysis steps with timing info

## Load libraries
library('getopt')

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

## Specify parameters
spec <- matrix(c(
    'experiment', 'e', 1, 'character', 'Experiment. Either snyder or hippo',
	'CovFile', 'd', 1, 'character', 'path to the .Rdata file with the results from loadCoverage() with cutoff >= 0',
	'chr', 'c', 1, 'character', 'Chromosome under analysis',
	'mcores', 'm', 1, 'integer', 'Number of cores',
	'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('snyder', 'hippo'))

## Format chromosome name appropriately
opt$chr <- mapSeqlevels(opt$chr, 'UCSC')


message('Loading Rdata file with the output from loadCoverage()')
load(opt$CovFile)

## Make it easy to use the name later. Here I'm assuming the names were generated using output='auto' in loadCoverage()
eval(parse(text=paste0('covData <- ', opt$chr, 'CovInfo')))
eval(parse(text=paste0('rm(', opt$chr, 'CovInfo)')))


## Load the models
load('models.Rdata')

## Load group information
load('groupInfo.Rdata')

if(file.exists('colsubset.Rdata')) {
    load('colsubset.Rdata')
} else {
    colsubset <- NULL
}


## Run the analysis
if (opt$experiment == 'snyder') {
    analyzeChr(chr = opt$chr, coverageInfo = covData, models = models, 
        cutoffFstat = 1e-05, colsubset = colsubset,
        nPermute = 100, seeds = seq_len(100) + 20131212, maxClusterGap = 3000,
        groupInfo = groupInfo, mc.cores = opt$mcores,
        lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'))
} else if (opt$experiment == 'hippo') {
    analyzeChr(chr = opt$chr, coverageInfo = covData, models = models, 
        cutoffFstat = 1e-04, colsubset = colsubset, cutoffPre = 3,
        nPermute = 100, seeds = seq_len(100) + 20131212, maxClusterGap = 3000,
        groupInfo = groupInfo, mc.cores = opt$mcores,
        lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'))
}


## Done
proc.time()
options(width = 120)
session_info()