recoup: An R package for the creation of complex genomic profile plots

Documented in recoup

recoup <- function(
    input,
    design=NULL,
    region=c("genebody","tss","tes","utr3","custom"),
    type=c("chipseq","rnaseq"),
    signal=c("coverage","rpm"),
    genome=c("hg18","hg19","hg38","mm9","mm10","rn5","rn6","dm3","dm6",
        "danrer7","danrer10","pantro4","pantro5","susscr3","susscr11",
        "ecucab2","tair10"),
    version="auto",
    refdb=c("ensembl","ucsc","refseq"),
    flank=c(2000,2000),
    onFlankFail=c("drop","trim"),
    fraction=1,
    orderBy=list(
        what=c("none","suma","sumn","maxa","maxn","avga","avgn","hcn"),
        order=c("descending","ascending"),
        custom=NULL
    ),
    #orderBy=getDefaultListArgs("orderBy"),
    binParams=list(
        flankBinSize=0,
        regionBinSize=0,
        sumStat=c("mean","median"),
        interpolation=c("auto","spline","linear","neighborhood"),
        binType=c("variable","fixed"),
        forceHeatmapBinning=TRUE,
        forcedBinSize=c(50,200),
        chunking=FALSE#,
        #seed=42
    ),
    #binParams=getDefaultListArgs("binParams"),
    selector=NULL,
    #selector=getDefaultListArgs("selector"),
    preprocessParams=list(
        fragLen=NA,
        cleanLevel=c(0,1,2,3),
        normalize=c("none","linear","downsample","sampleto"),
        sampleTo=1e+6,
        spliceAction=c("split","keep","remove"),
        spliceRemoveQ=0.75,
        #bedGenome=ifelse(genome %in% c("hg18","hg19","hg38","mm9","mm10","rn5",
        #    "dm3","danrer7","pantro4","susscr3"),genome,NA),
        bedGenome=NA#,
        #seed=42
    ),
    #preprocessParams=getDefaultListArgs("preprocessParams"),
    plotParams=list(
        plot=TRUE,
        profile=TRUE,
        heatmap=TRUE,
        correlation=TRUE,
        signalScale=c("natural","log2"),
        heatmapScale=c("common","each"),
        heatmapFactor=1,
        corrScale=c("normalized","each"),
        sumStat=c("mean","median"),
        smooth=TRUE,
        corrSmoothPar=ifelse(is.null(design),0.1,0.5),
        singleFacet=c("none","wrap","grid"),
        multiFacet=c("wrap","grid"),
        singleFacetDirection=c("horizontal","vertical"),
        conf=TRUE,
        device=c("x11","png","jpg","tiff","bmp","pdf","ps"),
        outputDir=".",
        outputBase=NULL
    ),
    #plotParams=getDefaultListArgs("plotParams",design),
    saveParams=list(
        ranges=TRUE,
        coverage=TRUE,
        profile=TRUE,
        profilePlot=TRUE,
        heatmapPlot=TRUE,
        correlationPlot=TRUE
    ),
    #saveParams=getDefaultListArgs("saveParams"),
    kmParams=list(
        k=0, # Do not perform kmeans
        nstart=20,
        algorithm=c("Hartigan-Wong","Lloyd","Forgy","MacQueen"),
        iterMax=20,
        reference=NULL#,
        #seed=42
    ),
    #kmParams=getDefaultListArgs("kmParams"),
    strandedParams=list(
        strand=NULL,
        ignoreStrand=TRUE
    ),
    #strandedParams=getDefaultListArgs("strandedParams"),
    ggplotParams=list(
        title=element_text(size=12),
        axis.title.x=element_text(size=10,face="bold"),
        axis.title.y=element_text(size=10,face="bold"),
        axis.text.x=element_text(size=9,face="bold"),
        axis.text.y=element_text(size=10,face="bold"),
        strip.text.x=element_text(size=10,face="bold"),
        strip.text.y=element_text(size=10,face="bold"),
        legend.position="bottom",
        panel.spacing=grid::unit(1,"lines")
    ),
    #ggplotParams=getDefaultListArgs("ggplotParams"),
    complexHeatmapParams=list(
        main=list(
            cluster_rows=ifelse(length(grep("hc",orderBy$what))>0,TRUE,FALSE),
            cluster_columns=FALSE,
            column_title_gp=grid::gpar(fontsize=10,font=2),
            show_row_names=FALSE,
            show_column_names=FALSE,
            heatmap_legend_param=list(
                color_bar="continuous"
            )
        ),
        group=list(
            cluster_rows=ifelse(length(grep("hc",orderBy$what))>0,TRUE,FALSE),
            cluster_columns=FALSE,
            column_title_gp=grid::gpar(fontsize=10,font=2),
            show_row_names=FALSE,
            show_column_names=FALSE,
            row_title_gp=grid::gpar(fontsize=8,font=2),
            gap=unit(5,"mm"),
            heatmap_legend_param=list(
                color_bar="continuous"
            )
        )
    ),
    #complexHeatmapParams=getDefaultListArgs("complexHeatmapParams"),
    bamParams=NULL,
    onTheFly=FALSE, # Directly from BAM w/o Ranges storing, also N/A with BED,
    #localDb=file.path(path.expand("~"),".recoup"),
    localDb=file.path(system.file(package="recoup"),"annotation.sqlite"),
    rc=NULL
) {
    # Simple check to throw error if user us using the old databasing system
    if (any(dir.exists(file.path(dirname(localDb),refdb))))
        stop("Possible old recoup database system detected. Please rebuild ",
            "using the new system or download a pre-built database. See also ",
            "the vignettes.")
    
    ############################################################################
    # Begin of simple parameter checking for a new or a restored object
    if (is.list(input) && !is.null(input$data)) { # Refeeding recoup object
        message("Detected a previous recoup output as input. Any existing ",
            "data requiring\nlengthy calculations (short reads, coverages and ",
            "profile matrices will be\nreused depending on other parameters. ",
            "If you want a complete recalculation,\neither input a fresh list ",
            "of arguments or remove any of the above with the\nremoveData ",
            "function.\n")
        prevCallParams <- input$callopts
        refRanges <- input$refRanges
        input <- input$data
    }
    else
        prevCallParams <- refRanges <- NULL
        
    if (!is.list(input) && file.exists(input))
        input <- readConfig(input)
    checkInput(input)
    if (is.null(names(input)))
        names(input) <- vapply(input,function(x) return(x$id),character(1))
    
    # Check if there are any mispelled or invalid parameters and throw a warning
    checkMainArgs(as.list(match.call()))
    
    # Check rest of arguments
    region <- tolower(region[1])
    refdb <- tolower(refdb[1])
    type <- tolower(type[1])
    onFlankFail <- tolower(onFlankFail[1])
    signal <- tolower(signal[1])
    if (!is.null(design) && !is.data.frame(design))
        checkFileArgs("design",design)
    if (!is.data.frame(genome) && !is.list(genome) && file.exists(genome))
        checkFileArgs("genome",genome)
    else if (is.character(genome)) {
        if (is.character(localDb) && file.exists(localDb)) {
            if (!.userOrg(genome,localDb))
                checkTextArgs("genome",genome,getSupportedOrganisms(),
                    multiarg=FALSE)
        }
    }
    if (is.character(localDb) && file.exists(localDb)) {
        if (!.userRefdb(refdb,localDb))
            checkTextArgs("refdb",refdb,getSupportedRefDbs(),multiarg=FALSE)
    }
    checkTextArgs("type",type,c("chipseq","rnaseq"),multiarg=FALSE)
    checkTextArgs("onFlankFail",onFlankFail,c("drop","trim"),multiarg=FALSE)
    checkTextArgs("signal",signal,c("coverage","rpm"),multiarg=FALSE)
    checkNumArgs("fraction",fraction,"numeric",c(0,1),"botheq")
    if (any(flank<0))
        stop("The minimum flanking allowed is 0 bp")
    if (any(flank>50000))
        stop("The maximum flanking allowed is 50000 bp")
    # Check the version argument
    version <- version[1]
    if (is.character(version)) {
        version <- tolower(version)
        checkTextArgs("version",version,c("auto"),multiarg=FALSE)
    }
    else
        checkNumArgs("version",version,"numeric")
    # If type is rnaseq, only genebody plots are valid
    if (type=="rnaseq" && region!="genebody") {
        warning("When type is \"rnaseq\", plots can be created only on ",
            "genebodies! Switching to genebody regions...",immediate.=TRUE)
        region <- "genebody"
    }
    
    # If type is rnaseq, read extension is illegal because of potential splicing
    if (type=="rnaseq" && !is.null(preprocessParams$fragLen)
        && !is.na(preprocessParams$fragLen)) {
        warning("When type is \"rnaseq\", read extension/reduction should not ",
            "happen because of potential splicing! Ignoring...",immediate.=TRUE)
        preprocessParams$fragLen <- NA
    }
    
    # If type is rnaseq, the only allowed genomes are the ones supported by
    # recoup for the time being. In the future, a custom RNA-Seq genome may be
    # constructed from a GFF or like...
    if (type=="rnaseq" && is.character(localDb) && file.exists(localDb) 
        && !.userOrg(genome,localDb) && !(genome %in% getSupportedOrganisms()))
        stop("When type is \"rnaseq\", only the supported or user imported ",
            "genomes are allowed!")

    # annotation must be a list to be fed to buildCustomAnnotation
    if (is.list(genome) && !is.data.frame(genome)) {
        # members are checked by buildCustomAnnotation if required
        # We only need to check that the gtfFile exists here
        if (!("gtf" %in% names(genome)))
            stop("A gtf field must be provided with an existing GTF file ",
                "when providing a list with custom annotation elements!")
        if ("gtf" %in% names(genome) && is.character(genome$gtf)
            && !file.exists(genome$gtf))
            stop("An existing GTF file must be provided when providing a ",
                "list with custom annotation elements!")
    }
        
    # Check what is going on with genome, first check if file, then localDb
    if (!is.data.frame(genome) && !is.list(genome) && is.character(genome)
        && !file.exists(genome) && is.character(localDb) 
        && file.exists(localDb)) {
        if (!.userOrg(genome,localDb))
            checkTextArgs("genome",genome,c("hg18","hg19","hg38","mm9","mm10",
                "rn5","rn6","dm3","dm6","danrer7","danrer10","pantro4",
                "pantro5","susscr3","susscr11","ecucab2","tair10"),
                multiarg=FALSE)
        if (!.userRefdb(refdb,localDb))
            checkTextArgs("refdb",refdb,c("ensembl","ucsc","refseq"),
                multiarg=FALSE)
    }
    
    # and list arguments
    orderByDefault <- getDefaultListArgs("orderBy")
    binParamsDefault <- getDefaultListArgs("binParams")
    selectorDefault <- getDefaultListArgs("selector")
    preprocessParamsDefault <- getDefaultListArgs("preprocessParams",
        genome=genome)
    plotParamsDefault <- getDefaultListArgs("plotParams",design=design)
    saveParamsDefault <- getDefaultListArgs("saveParams")
    kmParamsDefault <- getDefaultListArgs("kmParams")
    strandedParamsDefault <- getDefaultListArgs("strandedParams")
    
    orderBy <- setArg(orderByDefault,orderBy)
    binParams <- setArg(binParamsDefault,binParams)
    if (!is.null(selector))
        selector <- setArg(selectorDefault,selector)
    preprocessParams <- setArg(preprocessParamsDefault,preprocessParams)
    plotParams <- setArg(plotParamsDefault,plotParams)
    saveParams <- setArg(saveParamsDefault,saveParams)
    kmParams <- setArg(kmParamsDefault,kmParams)
    strandedParams <- setArg(strandedParamsDefault,strandedParams)
    
    orderBy <- validateListArgs("orderBy",orderBy)
    binParams <- validateListArgs("binParams",binParams)
    if (!is.null(selector))
        selector <- validateListArgs("selector",selector)
    preprocessParams <- validateListArgs("preprocessParams",preprocessParams)
    plotParams <- validateListArgs("plotParams",plotParams)
    saveParams <- validateListArgs("saveParams",saveParams)
    kmParams <- validateListArgs("kmParams",kmParams)
    strandedParams <- validateListArgs("strandedParams",strandedParams)
    
    if (is.null(plotParams$outputBase))
        plotParams$outputBase <- paste(vapply(input,function(x) return(x$id),
            character(1)),collapse="-")
    
    bbb <- vapply(input,function(x) return(tolower(x$format)=="bed"),logical(1))
    if (any(bbb) && !(preprocessParams$bedGenome %in% c("hg18","hg19","hg38",
            "mm9","mm10","rn5","dm3","danrer7","pantro4","susscr3")))
        stop("When short read files are in BED format, either the genome ",
            "parameter should be one\nof the supported organisms, or ",
            "preprocessParams$bedGenome must be specified as one of them.")
   
    # End of simple parameter checking for a new or a restored object
    ############################################################################
    
    ############################################################################
    # Begin of complex parameter storage procedure and parameter recall from a
    # previous call
    
    # Store all parameters to an obect for later reference to a next call, after
    # checking if something has changed (e.g. using setr)...
    thisCall <- as.list(match.call())[-1]
    if (!is.null(prevCallParams)) {
        callParams <- list(
            region=if (is.null(thisCall$region)) prevCallParams$region else
                region,
            type=if (is.null(thisCall$type)) prevCallParams$type else type,
            onFlankFail=if (is.null(thisCall$onFlankFail)) 
                prevCallParams$onFlankFail else onFlankFail,
            signal=if (is.null(thisCall$signal)) prevCallParams$signal 
                else signal,
            genome=if (is.null(thisCall$genome)) prevCallParams$genome else
                genome,
            refdb=if (is.null(thisCall$refdb)) prevCallParams$refdb else
                refdb,
            version=if (is.null(thisCall$version)) prevCallParams$version else
                version,
            flank=if (is.null(thisCall$flank)) prevCallParams$flank else flank,
            fraction=if (is.null(thisCall$fraction)) prevCallParams$fraction
                else fraction,
            orderBy=if (is.null(thisCall$orderBy)) prevCallParams$orderBy else
                orderBy,
            binParams=if (is.null(thisCall$binParams)) prevCallParams$binParams
                else binParams,
            selector=selector, # selector is a special case
            preprocessParams=if (is.null(thisCall$preprocessParams))
                prevCallParams$preprocessParams else preprocessParams,
            plotParams=if (is.null(thisCall$plotParams)) 
                prevCallParams$plotParams else plotParams,
            saveParams=if (is.null(thisCall$saveParams))
                prevCallParams$saveParams else saveParams,
            kmParams=if (is.null(thisCall$kmParams))
                prevCallParams$kmParams else kmParams,
            strandedParams=if (is.null(thisCall$strandedParams))
                prevCallParams$strandedParams else strandedParams,
            ggplotParams=if (is.null(thisCall$ggplotParams))
                prevCallParams$ggplotParams else ggplotParams,
            complexHeatmapParams=if (is.null(thisCall$complexHeatmapParams))
                prevCallParams$complexHeatmapParams else complexHeatmapParams,
            #bamParams=if (is.null(thisCall$bamParams),
            #    prevCallParams$bamParams,bamParams), ## Unused
            onTheFly=if (is.null(thisCall$onTheFly)) prevCallParams$onTheFly 
                else onTheFly,
            localDb=if (is.null(thisCall$localDb)) 
                prevCallParams$localDb else localDb,
            rc=if (is.null(thisCall$rc)) prevCallParams$rc else rc,
            customIsBase=NULL # Additional placeholder
        )
    }
    else {
        callParams <- list(
            region=region,
            type=type,
            onFlankFail=onFlankFail,
            signal=signal,
            genome=genome,
            version=version,
            refdb=refdb,
            flank=flank,
            fraction=fraction,
            orderBy=orderBy,
            binParams=binParams,
            selector=selector,
            preprocessParams=preprocessParams,
            plotParams=plotParams,
            saveParams=saveParams,
            kmParams=kmParams,
            strandedParams=strandedParams,
            ggplotParams=ggplotParams,
            complexHeatmapParams=complexHeatmapParams,
            bamParams=bamParams,
            onTheFly=onTheFly,
            localDb=localDb,
            rc=rc,
            customIsBase=NULL # Additional placeholder
        )
    }
    
    # ...and check if there has been a previous call and decide on big
    # recalculations...
    input <- decideChanges(input,callParams,prevCallParams)
    
    # Redefine all final arguments for this call
    region <- callParams$region
    type <- callParams$type
    onFlankFail <- callParams$onFlankFail
    signal <- callParams$signal
    genome <- callParams$genome
    version <- callParams$version
    refdb <- callParams$refdb
    flank <- callParams$flank
    fraction <- callParams$fraction
    orderBy <- callParams$orderBy
    binParams <- callParams$binParams
    selector <- callParams$selector
    preprocessParams <- callParams$preprocessParams
    plotParams <- callParams$plotParams
    saveParams <- callParams$saveParams
    kmParams <- callParams$kmParams
    strandedParams <- callParams$strandedParams
    ggplotParams <- callParams$ggplotParams
    complexHeatmapParams <- callParams$complexHeatmapParams
    bamParams <- callParams$bamParams
    onTheFly <- callParams$onTheFly
    localDb <- callParams$localDb
    rc <- callParams$rc
    # End of complex parameter storage procedure and parameter recall from a
    # previous call
    ############################################################################
    
    # Continue with actual work
    # If annotation already here and not removed from the input object because
    # of change in region, flank or type, no need to reload it and thus, recoup
    # object become even more mobile and portable as localDb is not needed.
    if (!is.null(refRanges)) {
        genomeRanges <- refRanges$genomeRanges
        helperRanges <- refRanges$helperRanges
    }
    else {
        # Annotation case #1: provided as a BED-like data frame with loci
        if (is.data.frame(genome)) {
            if (type=="rnaseq")
                stop("When type=\"rnaseq\", only the usage of a supported or ",
                    "user-imported organism from GTF file is allowed!")
            rownames(genome) <- as.character(genome[,4])
            genomeRanges <- makeGRangesFromDataFrame(
                df=genome,
                keep.extra.columns=TRUE
            )
            helperRanges <- NULL
        }
        # Annotation case #2: provided strictly as a BED-like file with loci
        else if (!is.list(genome) && !is.data.frame(genome) 
            && is.character(genome)) {
            if (file.exists(genome)) {
                if (type=="rnaseq")
                    stop("When type=\"rnaseq\", only the usage of a supported ",
                        "or user-imported organism from GTF file is allowed!")
                genome <- read.delim(genome)
                rownames(genome) <- as.character(genome[,4])
                genomeRanges <- makeGRangesFromDataFrame(
                    df=genome,
                    keep.extra.columns=TRUE
                )
                # Need to assigne Seqinfo...
                sf <- .chromInfoToSeqInfoDf(.chromInfoFromBAM(input[[1]]$file),
                    asSeqinfo=TRUE)
                seqlevels(genomeRanges) <- seqlevels(sf)
                seqinfo(genomeRanges) <- sf
                
                helperRanges <- NULL
            }
            # Annotation case #3: use local database or automatically on-the-fly
            else {
                if (type=="chipseq") {
                    if (region == "utr3")
                        genomeRanges <- tryCatch(loadAnnotation(genome,refdb,
                            type="utr",version=version,summarized=TRUE,
                            db=localDb,rc=rc),
                            error=function(e) {
                                tryCatch({
                                    gtfFile <- genome$gtf
                                    metadata <- genome
                                    metadata$gtf <- NULL
                                    genomeRanges <- importCustomAnnotation(
                                        gtfFile,metadata,"utr")
                                },error=function(e) {
                                    message("Error ",e)
                                    stop("Please provide an existing organism ",
                                        "or a list with annotation metadata ",
                                        "and GTF file!")
                                },finally="")
                            },finally="")
                    else
                        genomeRanges <- tryCatch(loadAnnotation(genome,refdb,
                            type="gene",version=version,db=localDb,rc=rc),
                            error=function(e) {
                                tryCatch({
                                    gtfFile <- genome$gtf
                                    metadata <- genome
                                    metadata$gtf <- NULL
                                    geneData <- importCustomAnnotation(gtfFile,
                                        metadata,"gene")
                                },error=function(e) {
                                    message("Error ",e)
                                    stop("Please provide an existing organism ",
                                        "or a list with annotation metadata ",
                                        "and GTF file!")
                                },finally="")
                            },finally="")

                    helperRanges <- NULL
                }
                else if (type=="rnaseq") {
                    genomeRanges <- tryCatch(loadAnnotation(genome,refdb,
                        type="exon",version=version,summarized=TRUE,db=localDb,
                        rc=rc),
                        error=function(e) {
                            tryCatch({
                                gtfFile <- genome$gtf
                                metadata <- genome
                                metadata$gtf <- NULL
                                genomeRanges <- importCustomAnnotation(gtfFile,
                                    metadata,"exon")
                            },error=function(e) {
                                message("Error ",e)
                                stop("Please provide an existing organism or ",
                                        "a list with annotation metadata and ",
                                        "GTF file!")
                            },finally="")
                        },finally="")
                    helperRanges <- tryCatch(loadAnnotation(genome,refdb,
                        type="gene",version=version,db=localDb,rc=rc),
                        error=function(e) {
                            tryCatch({
                                gtfFile <- genome$gtf
                                metadata <- genome
                                metadata$gtf <- NULL
                                geneData <- importCustomAnnotation(gtfFile,
                                    metadata,"gene")
                            },error=function(e) {
                                message("Error ",e)
                                stop("Please provide an existing organism or ",
                                    "a list with annotation metadata and GTF ",
                                    "file!")
                            },finally="")
                        },finally="")
                    if (!is(genomeRanges,"GRangesList"))
                        genomeRanges <- split(genomeRanges,genomeRanges$gene_id)
                    helperRanges <- helperRanges[names(genomeRanges)]
                }
            }
        }
        
        # For saving...
        refRanges <- list(genomeRanges=genomeRanges,helperRanges=helperRanges)
    }
    
    # After genome read, check if we have a valid custom orderer
    if (!is.null(orderBy$custom)) {
        if (length(orderBy$custom) != length(genomeRanges)) {
            warning("The custom orderer provided with orderBy parameter does ",
                "not have length equal to the number\nof elements in the ",
                "interrogated genomic regions and will be ignored!",
                immediate.=TRUE)
            orderBy$custom <- NULL
        }
    }
    
    # Read and check design compatibilities. Check if k-means is requested and
    # message accordingly. If k-means is requested it will be added to the 
    # design data frame
    #hasDesign <- FALSE
    if (!is.null(design)) {
        if (!is.data.frame(design))
            design <- read.delim(design,row.names=1)
        nfac <- ncol(design)
        if (length(input)>1 && nfac>2)
            stop("When more than one files are provided for coverage ",
                "generation, the maximum number of allowed design factors is 2")
        if (length(input)>1 && nfac>1 && kmParams$k>0)
            stop("When more than one files are provided for coverage ",
                "generation and k-means clustering is also requested, the ",
                "maximum number of allowed design factors is 1")
        if (length(input)==1 && nfac>3)
            stop("The maximum number of allowed design factors is 3")
        if (length(input)==1 && nfac>2 && kmParams$k>0)
            stop("The maximum number of allowed design factors when k-means ",
                "clustering is requested is 2")
        if (length(input)==1 && nfac>2 && plotParams$singleFacet!="none")
            stop("The maximum number of allowed design factors whith one ",
                "sample but wishing a gridded profile layout is 2")
        if (length(input)==1 && nfac>1 && kmParams$k>0 
            && plotParams$singleFacet!="none")
            stop("The maximum number of allowed design factors whith one ",
                "sample but wishing for k-means clustering and gridded ",
                "profile layout is 1")
        # Reduce the genomeRanges according to design or the other way
        if (nrow(design)>length(genomeRanges))
            design <- tryCatch({
                design[names(genomeRanges),,drop=FALSE]
            },error=function(e) {
                stop("Unexpected error occured! Are you sure that element ",
                    "(row) names in the design file are of the same type as ",
                    "the genome file?")
            },finally={})
        else if (nrow(design)<=length(genomeRanges)) {
            genomeRanges <- tryCatch({
                genomeRanges[rownames(design)]
            },error=function(e) {
                stop("Unexpected error occured! Are you sure that element ",
                    "(row) names in the design file are of the same type as ",
                    "the genome file?")
            },finally={})
            if (type=="rnaseq")
                helperRanges <- tryCatch({
                    helperRanges[rownames(design)]
                },error=function(e) {
                    stop("Unexpected error occured! Are you sure that element ",
                        "(row) names in the design file are of the same type ",
                        "as the genome file?")
                },finally={})
        }
        # ...but maybe the names are incompatible
        if (length(genomeRanges)==0)
            stop("No ranges left after using the identifiers provided with ",
                "the design file. Are you sure the identifiers between the ",
                "two files are compatible?")
        if (nrow(design)==0)
            stop("No design elements left after using the identifiers ",
                "provided with the genome file. Are you sure the identifiers ",
                "between the two files are compatible?")
    }
    
    # Apply the rest of the filters if any to reduce later computational burden
    if (!is.null(selector)) {
        if (type=="chipseq")
            genomeRanges <- applySelectors(genomeRanges,selector,rc=rc)
        if (type=="rnaseq") {
            helperRanges <- applySelectors(helperRanges,selector)
            # Filter and align names if we have helperRanges around
            genomeRanges <- genomeRanges[names(helperRanges)]
        }
    }
    
    # Now we must follow two paths according to region type, genebody and custom
    # areas with equal/unequal widths and utr3 or tss, tes and 1-width custom 
    # areas.
    callParams$customIsBase <- FALSE
    if (region=="custom" && all(width(genomeRanges)==1)) {
        if (all(flank==0)) {
            warning("Flanking cannot be zero bp in both directions when the ",
                "reference region is only 1bp! Setting to default ",
                "(2000,2000)...",immediate.=TRUE)
            flank <- c(2000,2000)
            callParams$flank <- flank
        }
        callParams$customIsBase <- TRUE
    }
    if (region %in% c("tss","tes")) {
        if (all(flank==0)) {
            warning("Flanking cannot be zero bp in both directions when the ",
                "reference region is \"tss\" or \"tes\"! Setting to default ",
                "(2000,2000)...", immediate.=TRUE)
            flank <- c(2000,2000)
            callParams$flank <- flank
        }
    }
    
    # Here we must determine the final ranges to pass to preprocessRanges and
    # then work with these later on
    intermRanges <- getMainRanges(genomeRanges,helperRanges=helperRanges,type,
        region,flank,onFlankFail,rc=rc)
    # It is very important that all the Ranges have Seqinfo objects attached as
    # they have to be trimmed for potential extensions (due to flanking regions)
    # beyond chromosome lengths. The built-in annotations do provide this option
    # otherwise when a custom genome/areas is/are provided, the genome of 
    # interest should be provided too. Or if not provided, the user will be
    # responsible for any crash.
    # TODO: When input ranges are custom, genome should be given to make Seqinfo
    
    # Should not be needed anymore with the new onFlankFail option
    #if (type == "chipseq") # Otherwise, throw error with GRangesList
    #    mainRanges <- trim(intermRanges$mainRanges)
    #else if (type == "rnaseq")
    #    mainRanges <- intermRanges$mainRanges
    #bamRanges <- trim(intermRanges$bamRanges)
    
    mainRanges <- intermRanges$mainRanges
    bamRanges <- intermRanges$bamRanges
    hadOffBound <- intermRanges$offBound
    # We must also correct the design in case of drop
    if (!is.null(design) && onFlankFail == "drop" && hadOffBound)
        design <- design[names(mainRanges),,drop=FALSE]
    # Say something if off bound detected
    if (hadOffBound) {
        if (onFlankFail == "drop")
        warning("Off bound regions detected after flanking! Offending ",
            "regions will be dropped from the reference\n  and the design if ",
            "provided...",immediate.=TRUE)
    else if (onFlankFail == "trim")
        warning("Off bound regions detected after flanking! Offending ",
            "reference regions will be trimmed\n  and the flanking regions ",
            "will be merged with the main...",immediate.=TRUE)
    }
    
    # Here we must write code for the reading and normalization of bam files
    # The preprocessRanges function looks if there is a valid (not null) ranges
    # field in input
    if (!onTheFly)
        input <- preprocessRanges(input,preprocessParams,genome,bamRanges,
            bamParams=bamParams,rc=rc)
    # At this point we must apply the fraction parameter if <1. We choose this
    # point in order not to restrict later usage of the read ranges and since it
    # does not take much time to load into memory. In addition, this operation
    # cannot be applied when streaming BAMs. In that case, user must subsample
    # outside recoup.
    if (!onTheFly && fraction<1) {
        newSize <- round(fraction*length(mainRanges))
        #set.seed(preprocessParams$seed)
        refIndex <- sort(sample(length(mainRanges),newSize))
        mainRanges <- mainRanges[refIndex]
        for (i in seq_len(length(input))) {
            if (!is.null(input[[i]]$ranges)) {
                newSize <- round(fraction*length(input[[i]]$ranges))
                #set.seed(preprocessParams$seed)
                fracIndex <- sort(sample(length(input[[i]]$ranges),newSize))
                input[[i]]$ranges <- input[[i]]$ranges[fracIndex]
            }
            if (!is.null(input[[i]]$coverage))
                input[[i]]$coverage <- input[[i]]$coverage[refIndex]
            if (!is.null(input[[i]]$profile))
                input[[i]]$profile <- input[[i]]$profile[refIndex,]
        }
    }

    # Remove unwanted seqnames from reference ranges (if not on the fly) and
    # properly subset ranges if not all chromosomes are provided with BAM
    if (!onTheFly) {
        chrs <- unique(unlist(lapply(input,function(x) {
            if (x$format %in% c("bam","bed"))
                return(as.character(runValue(seqnames(x$ranges))))
            else if (x$format=="bigwig") {
                if (!requireNamespace("rtracklayer"))
                    stop("R package rtracklayer is required to read and ",
                        "import BigWig files!")
                return(as.character(seqnames(seqinfo(BigWigFile(x$file)))))
            }
        })))
        if (type=="chipseq") {
            #keep <- which(as.character(seqnames(mainRanges)) %in% chrs)
            #mainRanges <- mainRanges[keep]
            mainRanges <- .subsetGRangesByChrs(mainRanges,chrs)
        }
        else if (type=="rnaseq") {
            #keeph <- which(as.character(seqnames(helperRanges)) %in% chrs)
            #helperRanges <- helperRanges[keeph]
            #mainRanges <- mainRanges[names(helperRanges)]
            helperRanges <- .subsetGRangesByChrs(helperRanges,chrs)
			mainRanges <- .subsetGRangesByChrs(mainRanges,chrs)
        }
    }
    
    if(signal == "coverage")
        input <- coverageRef(input,mainRanges,strandedParams,rc=rc)
    else if (signal == "rpm")
        input <- rpMatrix(input,mainRanges,flank,binParams,strandedParams,rc=rc)
    
    # If normalization method is linear, we must adjust the coverages
    # TODO: Check for onTheFly in the beginning
    if (preprocessParams$normalize == "linear" && !onTheFly) {
        message("Calculating normalization factors...")
        linFac <- calcLinearFactors(input,preprocessParams)
        for (n in names(input)) {
            message("  normalizing ",n)
            if (linFac[n]==1)
                next
            if (signal == "coverage")
                input[[n]]$coverage <- cmclapply(input[[n]]$coverage,
                function(x,f) {
                    return(x*f)
                },linFac[n])
            else if (signal == "rpm")
                input[[n]]$profile <- apply(input[[n]]$profile,1,function(x,f) {
                    return(x*f)
                },linFac[n])
        }
    }
    
    # Now we must summarize and create the matrices. If genebody, utr3 or 
    # unequal custom lengths, bin is a must, else we leave to user
    mustBin <- FALSE
    if (region=="genebody" || region=="utr3")
        mustBin <- TRUE
    if (region=="custom") {
        w <- width(mainRanges)
        if (any(w!=w[1]))
            mustBin <- TRUE
    }
    if (mustBin) {
        if (binParams$regionBinSize==0) {
            warning("Central region bin size not set for a region that must ",
                "be binned! Setting to 200...",immediate.=TRUE)
            binParams$regionBinSize <- 200
        }
    }
    
    # Chunking?
    if (signal == "coverage") { # Otherwise, matrix calculate by rpm
        if (binParams$chunking) {
            if (type=="chipseq")
                binParams$chunks <- split(seq_len(length(genomeRanges)),
                    as.character(seqnames(genomeRanges)))
            else if (type=="rnaseq")
                binParams$chunks <- split(seq_len(length(helperRanges)),
                    as.character(seqnames(helperRanges)))
        }
        # hadOffBounds must be passed on input at some point
        fe0 <- FALSE
        if (onFlankFail == "trim" && hadOffBound)
            fe0 <- TRUE
        input <- profileMatrix(input,flank,binParams,rc,.feNoSplit=fe0)
    }
    
    # In some strange glimpses, we are getting very few NaNs in profile matrix
    # which I was unable to reproduce on a gene by gene basis. If no NaNs are
    # detected, no action is performed in the input object. Also, these NaNs
    # seem to appear only on zero count gene profiles, so it's safe to impute
    # by zero.
    #input <- imputeProfile(input,method="simple",rc)
    input <- imputeZero(input)
    
    # Perform the k-means clustering if requested and append to design (which
    # has been checked, if we are allowed to do so)
    if (kmParams$k>0)
        design <- kmeansDesign(input,design,kmParams)
    
    # Coverages and profiles calculated... Now depending on plot option, we go 
    # further or return the enriched input object for saving
    if (!plotParams$profile && !plotParams$heatmap && !plotParams$correlation) {
        recoupObj <- toOutput(input,design,saveParams,callParams=callParams,
            refRanges=refRanges)
        return(recoupObj)
    }
    else {
        recoupObj <- toOutput(input,design,list(ranges=TRUE,coverage=TRUE,
            profile=TRUE),callParams=callParams,refRanges=refRanges)
    }
            
    ## Our plot objects
    recoupPlots <- list()
    
    # We must pass the matrices to plotting function
    if (plotParams$profile) {
        message("Constructing genomic coverage profile curve(s)")
        #theProfile <- recoupProfile(recoupObj,rc=rc)
        recoupObj <- recoupProfile(recoupObj,rc=rc)
        theProfile <- getr(recoupObj,"profile")
        recoupPlots$profilePlot <- theProfile
    }
    
    # Some default binning MUST be applied for the heatmap... Otherwise it could
    # take insanely long time and space to draw/store
    if (plotParams$heatmap) {
        # Inform the user about enforced binning (or not)
        if (region %in% c("tss","tes") || callParams$customIsBase) {
            if (binParams$regionBinSize==0 && binParams$forceHeatmapBinning) {
                message("The resolution of the requested profiles will be ",
                    "lowered to avoid\nincreased computation time and/or ",
                    "storage space for heatmap profiles...")
                
            }
            else if (binParams$regionBinSize==0
                && !binParams$forceHeatmapBinning)
                warning("forceHeatmapBinning is turned off for high ",
                    "resolution plotting. Be prepared for\nlong computational ",
                    "times and big figures!",immediate.=TRUE)
        }
        else {
            if ((binParams$regionBinSize==0 || binParams$flankBinSize==0)
                && binParams$forceHeatmapBinning) {
                message("The resolution of the requested profiles will be ",
                    "lowered to avoid\nincreased computation time and/or ",
                    "storage space for heatmap profiles...")
                
            }
            else if ((binParams$regionBinSize==0 || binParams$flankBinSize==0)
                && !binParams$forceHeatmapBinning)
                warning("forceHeatmapBinning is turned off for high ",
                    "resolution plotting. Be prepared for\nlong computational ",
                    "times and big figures!",immediate.=TRUE)
        }
        
        if (binParams$forceHeatmapBinning 
            && (binParams$regionBinSize==0 || binParams$flankBinSize==0)) {
            helpInput <- recoupObj$data
            if (region %in% c("tss","tes") || callParams$customIsBase) {
                for (n in names(helpInput)) {
                    message("Calculating ",region," profile for ",
                        helpInput[[n]]$name)
                    helpInput[[n]]$profile <- 
                        binCoverageMatrix(helpInput[[n]]$coverage,
                            binSize=binParams$forcedBinSize[2],
                            stat=binParams$sumStat,rc=rc)
                    helpInput[[n]]$profile <- helpInput[[n]]$profile
                }
            }
            else {
                for (n in names(helpInput)) {
                    message("Calculating ",region," profile for ",
                        helpInput[[n]]$name)
                    message(" center")
                    #center <- binCoverageMatrix(helpInput[[n]]$coverage$center,
                    #    binSize=binParams$forcedBinSize[2],
                    #    stat=binParams$sumStat,rc=rc)
                    center <- binCoverageMatrix(
                        input[[n]]$coverage,binSize=binParams$forcedBinSize[2],
                        stat=binParams$sumStat,
                        interpolation=binParams$interpolation,
                        flank=flank,where="center",rc=rc
                    )
                    message(" upstream")
                    #left <- binCoverageMatrix(helpInput[[n]]$coverage$upstream,
                    #    binSize=binParams$forcedBinSize[1],
                    #    stat=binParams$sumStat,rc=rc)
                    left <- binCoverageMatrix(
                        input[[n]]$coverage,binSize=binParams$forcedBinSize[1],
                        stat=binParams$sumStat,
                        interpolation=binParams$interpolation,flank=flank,
                        where="upstream",rc=rc
                    )
                    message(" downstream")
                    #right <- binCoverageMatrix(
                    #    helpInput[[n]]$coverage$downstream,
                    #    binSize=binParams$forcedBinSize[1],
                    #    stat=binParams$sumStat,rc=rc)
                    right <- binCoverageMatrix(
                        input[[n]]$coverage,binSize=binParams$forcedBinSize[1],
                        stat=binParams$sumStat,
                        interpolation=binParams$interpolation,flank=flank,
                        where="downstream",rc=rc
                    )
                    helpInput[[n]]$profile <- cbind(left,center,right)
                    rownames(helpInput[[n]]$profile) <-
                        names(input[[n]]$coverage)
                    helpInput[[n]]$profile <- helpInput[[n]]$profile
                }
            }
        }
        else
            helpInput <- recoupObj$data
        
        helpObj <- recoupObj
        helpObj$data <- helpInput
        message("Constructing genomic coverage heatmap(s)")
        #theHeatmap <- recoupHeatmap(helpObj,rc=rc)
        helpObj <- recoupHeatmap(helpObj,rc=rc)
        theHeatmap <- getr(helpObj,"heatmap")
        recoupObj <- setr(recoupObj,"heatmap",theHeatmap)
        recoupPlots$heatmapPlot <- theHeatmap
        
        # Derive the main heatmap in case of hierarchical clustering        
        mainh <- 1
        if (length(grep("hc",orderBy$what))>0) {
            nc <- nchar(orderBy$what)
            mh <- suppressWarnings(as.numeric(substr(orderBy$what,nc,nc)))
            if (is.na(mh))
                warning("Reference profile for hierarchical clustering order ",
                    "not recognized! Using the 1st...",immediate.=TRUE)
            else if (mh > length(input)) {
                warning("Reference profile (",mh,") for hierarchical ",
                    "clustering order does not exist (the input has only ",
                    length(input)," sources! Using the last...",
                    immediate.=TRUE)
                    mainh <- length(input)
            }
            else
                mainh <- mh
        }
    }
    
    # We must pass the matrices to plotting function
    if (plotParams$correlation) {
        message("Constructing coverage correlation profile curve(s)")
        recoupObj <- recoupCorrelation(recoupObj,rc=rc)
        theCorrelation <- getr(recoupObj,"correlation")
        recoupPlots$correlationPlot <- theCorrelation
    }
    
    # Overwrite final object so as to return it
    recoupObj <- toOutput(input,design,saveParams,recoupPlots,callParams,
        refRanges)
    
    # Make any plots asked
    if (plotParams$plot) {
        what <- character(0)
        if (plotParams$profile)
            what <- c(what,"profile")
        if (plotParams$heatmap)
            what <- c(what,"heatmap")
        if (plotParams$correlation)
            what <- c(what,"correlation")
        if (length(what)>0)
            recoupPlot(recoupObj,what,plotParams$device,plotParams$outputDir,
                plotParams$outputBase,mainh)
    }
    
    # Return the enriched input object according to save options
    return(recoupObj)
}

applySelectors <- function(ranges,selector,rc=NULL) {
    if (!is.null(selector$id)) {
        ranges <- ranges[selector$ids]
        if (length(ranges)==0)
            stop("No ranges left after using the identifiers provided with ",
                "the selector field. Are you sure the identifiers between the ",
                "two files are compatible?")
    }
    if (!is.null(selector$biotype) && !is.null(ranges$biotype)) {
        good <- which(ranges$biotype %in% selector$biotype)
        ranges <- ranges[good]
        if (length(ranges)==0)
            stop("No ranges left after using the biotypes provided with the ",
                "selector field. Are you sure the identifiers between the two ",
                "files are compatible?")
    }
    if (!is.null(selector$exonType) && !is.null(ranges$exon_type)) {
        good <- which(ranges$exonType %in% selector$exonType)
        ranges <- ranges[good]
        if (length(ranges)==0)
            stop("No ranges left after using the exon types provided with ",
                "the selector field. Are you sure the identifiers between the ",
                "two files are compatible?")
    }
    return(ranges)
}

toOutput <- function(input,design=NULL,saveParams,plotObjs=NULL,
    callParams=NULL,refRanges=NULL) {
    if (!saveParams$ranges)
        input <- removeData(input,"ranges")
    if (!saveParams$coverage)
        input <- removeData(input,"coverage")
    if (!saveParams$profile)
        input <- removeData(input,"profile")
    plots <- list(profile=NULL,heatmap=NULL,correlation=NULL)
    if (!is.null(plotObjs) && saveParams$profilePlot)
        plots$profile <- plotObjs$profilePlot
    if (!is.null(plotObjs) && saveParams$heatmapPlot)
        plots$heatmap <- plotObjs$heatmapPlot
    if (!is.null(plotObjs) && saveParams$correlationPlot)
        plots$correlation <- plotObjs$correlationPlot
    return(list(
        data=input,
        design=design,
        plots=plots,
        callopts=callParams,
        refRanges=refRanges
    ))
}