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 shulha or epimap',
	'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)
}

## For testing
if(FALSE) {
    # qrsh -l mem_free=4G,h_vmem=6G,h_fsize=10G -pe local 2
    # module load R/3.3
    # cd /dcl01/lieber/ajaffe/derRuns/derChIP/epimap/derAnalysis/run1-v1.5.38-H3K27ac
    opt <- list(experiment = 'epimap', chr = 'chrY', mcores = 2,
    CovFile = '/dcl01/lieber/ajaffe/derRuns/derChIP/epimap/CoverageInfo/chrYCovInfo-filtered.Rdata'
    )
}

## Check experiment input
stopifnot(opt$experiment %in% c('shulha', 'epimap'))

## Format chromosome name appropriately
if(opt$chr == 'M') opt$chr <- 'chrM' else 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 == 'shulha') {
    run <- gsub('/dcl01/lieber/ajaffe/derRuns/derChIP/shulha/derAnalysis/', '', getwd())
    if(run == 'run1-v1.0.10') {
        analyzeChr(chr = opt$chr, coverageInfo = covData, models = models,
            cutoffFstat = 1e-04, colsubset = colsubset, cutoffPre = 5,
            nPermute = 100, seeds = seq_len(100) + 150422, maxClusterGap = 3000,
            groupInfo = groupInfo, mc.cores = opt$mcores,
            lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'))
    } else if (run == 'run2-v1.0.10') {
        analyzeChr(chr = opt$chr, coverageInfo = covData, models = models,
            cutoffFstat = 0.05, colsubset = colsubset, cutoffPre = 2,
            nPermute = 100, seeds = seq_len(100) + 150422, maxClusterGap = 3000,
            groupInfo = groupInfo, mc.cores = opt$mcores,
            lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'))
    } else if (run == 'run3-v1.0.10') {
        analyzeChr(chr = opt$chr, coverageInfo = covData, models = models,
            cutoffFstat = 0.01, colsubset = colsubset, cutoffPre = 2,
            nPermute = 100, seeds = seq_len(100) + 150422, maxClusterGap = 3000,
            groupInfo = groupInfo, mc.cores = opt$mcores,
            lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'))
    } else if (run == 'run4-v1.0.10') {
        analyzeChr(chr = opt$chr, coverageInfo = covData, models = models,
            cutoffFstat = 0.01, colsubset = colsubset, cutoffPre = 2,
            nPermute = 100, seeds = seq_len(100) + 150422, maxClusterGap = 3000,
            groupInfo = groupInfo, mc.cores = opt$mcores,
            lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'), scalefac=1)
    } else if(run == 'run5-v1.5.34') {
        analyzeChr(chr = opt$chr, coverageInfo = covData, models = models,
            cutoffFstat = 0.05, colsubset = colsubset, cutoffPre = 2,
            nPermute = 100, seeds = seq_len(100) + 150422, maxClusterGap = 3000,
            groupInfo = groupInfo, mc.cores = opt$mcores,
            lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'),
            smooth = TRUE, minNum = 36, bpSpan = 1000, minInSpan = 36)
    } else if (run == 'run6-v1.5.34') {
        analyzeChr(chr = opt$chr, coverageInfo = covData, models = models,
            cutoffFstat = 0.01, colsubset = colsubset, cutoffPre = 2,
            nPermute = 100, seeds = seq_len(100) + 150422, maxClusterGap = 3000,
            groupInfo = groupInfo, mc.cores = opt$mcores,
            lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'),
            smooth = TRUE, minNum = 36, bpSpan = 1000, minInSpan = 36)
    } else if (run == 'run7-v1.5.35') {
        analyzeChr(chr = opt$chr, coverageInfo = covData, models = models,
            cutoffFstat = 0.05, colsubset = colsubset, cutoffPre = 2,
            nPermute = 100, seeds = seq_len(100) + 150422, maxClusterGap = 3000,
            groupInfo = groupInfo, mc.cores = opt$mcores,
            lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'),
            smooth = TRUE, minNum = 36, bpSpan = 300, minInSpan = 36)
    } else if (run == 'run8-v1.5.35') {
        analyzeChr(chr = opt$chr, coverageInfo = covData, models = models,
            cutoffFstat = 0.01, colsubset = colsubset, cutoffPre = 2,
            nPermute = 100, seeds = seq_len(100) + 150422, maxClusterGap = 3000,
            groupInfo = groupInfo, mc.cores = opt$mcores,
            lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'),
            smooth = TRUE, minNum = 36, bpSpan = 300, minInSpan = 36)
    }
} else if (opt$experiment == 'epimap') {
    analyzeChr(chr = opt$chr, coverageInfo = covData, models = models,
        cutoffFstat = 0.01, colsubset = colsubset, cutoffPre = 10,
        nPermute = 100, seeds = seq_len(100) + 160427, maxClusterGap = 3000,
        groupInfo = groupInfo, mc.cores = opt$mcores,
        lowMemDir = file.path(tempdir(), opt$chr, 'chunksDir'), chunksize = 5e5,
        smooth = TRUE, minNum = 100, bpSpan = 300, minInSpan = 100)
}


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