fix-calculatePvalues.R

## Fix as best as possible the bug in calculatePvalues() regarding the default
## maxRegionGap value.

## Usage:
# Step 1: open R in the main results folder
# Step 2: use fixChrs() with options of your choice
# fixChrs('chr')

## Alternatively use fixRegions() manually and re-run 
## bumphunter::annotateNearest()

## Once all chrs have been fixed, re-run mergeResults()

library('derfinder')


fixRegions <- function(regions, fstats, models, coveragePrep, chr,
    cutoffFstat = 1e-08, cutoffType = 'theoretical',
    significantCut = c(0.05, 0.1), maxRegionGap = 0L, maxClusterGap = 300L,
    lowMemDir = NULL) {
    stopifnot(identical(names(regions),
        c('regions', 'nullStats', 'nullWidths', 'nullPermutation')))
    cutoff <- derfinder:::.calcFstatCutoff(cutoffType, cutoffFstat, fstats,
        models)
    message(paste(Sys.time(), "F-stat cutoff to be used:", cutoff))
    
    new <- calculatePvalues(coveragePrep, models, fstats, nPermute = 0,
        chr = chr, cutoff = cutoff, maxRegionGap = maxRegionGap,
        maxClusterGap = maxClusterGap, lowMemDir = lowMemDir)
    
    new$regions$pvalues <- derfinder:::.calcPval(new$regions$area,
        as.numeric(regions$nullStats * regions$nullWidths))
    new$regions$significant <- factor(new$regions$pvalues < significantCut[1], 
        levels = c(TRUE, FALSE))
    
    qvalues <- qvalue::qvalue(new$regions$pvalues)
    if (is(qvalues, "qvalue")) {
        qvalues <- qvalues$qvalues
        sigQval <- factor(qvalues < significantCut[2], levels = c(TRUE, 
            FALSE))
    } else {
        message(paste(Sys.time(), "skipping q-value calculation."))
        qvalues <- rep(NA, length(regs$pvalues))
        sigQval <- rep(NA, length(regs$pvalues))
    }
    new$regions$qvalues <- qvalues
    new$regions$significantQval <- sigQval
    
    new$nullStats <- regions$nullStats
    new$nullWidths <- regions$nullWidths
    new$nullPermutation <- regions$nullPermutation
    new$note <- "Object fixed for bug in calculatePvalues(). Areas from null regions are approximated."
    return(new)
}

fixChrs <- function(chr.pattern = 'chr', significantCut = c(0.05, 0.1),
    maxRegionGap = 0L, maxClusterGap = 300L, subject = 'hg19',
    runAnnotation = TRUE, lowMemDir = NULL) {
    dirs <- dir(pattern = chr.pattern)
    for(chr in dirs) {
        message(paste(Sys.time(), 'fixing chromosome', chr))
        load(file.path(chr, 'coveragePrep.Rdata'))
        load(file.path(chr, 'regions.Rdata'))
        load(file.path(chr, 'fstats.Rdata'))
        load(file.path(chr, 'optionsStats.Rdata'))
        regions.original <- regions
        
        ## Save original regions for backup purposes
        save(regions.original, file = file.path(chr, 'regions-original.Rdata'))
        
        ## Fix the regions. Note p-values will be approximated given that
        ## there is no way to truly fix the null areas without re-running
        ## derfinder
        regions <- fixRegions(regions = regions.original, fstats = fstats,
            models = optionsStats$models, coveragePrep = prep,
            chr = chr, cutoffFstat = optionsStats$cutoffFstat,
            cutoffType = optionsStats$cutoffType,
            significantCut = significantCut, maxRegionGap = maxRegionGap,
            maxClusterGap = maxClusterGap, lowMemDir = lowMemDir)
        save(regions, file = file.path(chr, 'regions.Rdata'))
        
        if(runAnnotation) {
            load(file.path(chr, 'annotation.Rdata'))
            annotation.original <- annotation
            save(annotation.original, file = file.path(chr,
                'annotation-original.Rdata'))
            annotation <- bumphunter::annotateNearest(regions$regions, subject)
            save(annotation, file = file.path(chr, 'annotation.Rdata'))
        }
    }
    return(invisible(NULL))
}