@@ -148,10 +148,12 @@ recoup <- function(

             "of arguments or remove any of the above with the\nremoveData ",

             "function.\n")

         prevCallParams <- input$callopts

+        refRanges <- input$refRanges

         input <- input$data

     }

     else

-        prevCallParams <- NULL

+        prevCallParams <- refRanges <- NULL

+        

     if (!is.list(input) && file.exists(input))

         input <- readConfig(input)

     checkInput(input)

@@ -409,42 +411,91 @@ recoup <- function(

     ############################################################################

     # Continue with actual work

-    # Annotation case #1: provided strictly 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

-        )

+    # 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

     }

-    # 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)) {

+    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!")

-            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")

+        # 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="utr",version=version,summarized=TRUE,db=localDb,

+                        type="exon",version=version,summarized=TRUE,db=localDb,

                         rc=rc),

                         error=function(e) {

                             tryCatch({

@@ -452,16 +503,15 @@ recoup <- function(

                                 metadata <- genome

                                 metadata$gtf <- NULL

                                 genomeRanges <- importCustomAnnotation(gtfFile,

-                                    metadata,"utr")

+                                    metadata,"exon")

                             },error=function(e) {

                                 message("Error ",e)

                                 stop("Please provide an existing organism or ",

-                                    "a list with annotation metadata and GTF ",

-                                    "file!")

+                                        "a list with annotation metadata and ",

+                                        "GTF file!")

                             },finally="")

                         },finally="")

-                else

-                    genomeRanges <- tryCatch(loadAnnotation(genome,refdb,

+                    helperRanges <- tryCatch(loadAnnotation(genome,refdb,

                         type="gene",version=version,db=localDb,rc=rc),

                         error=function(e) {

                             tryCatch({

@@ -477,54 +527,21 @@ recoup <- function(

                                     "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)]

+                    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 provide with orderBy parameter does ",

+            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)

@@ -666,7 +683,7 @@ recoup <- function(

     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 thre main...",immediate.=TRUE)

+            "will be merged with the main...",immediate.=TRUE)

     }

     # Here we must write code for the reading and normalization of bam files

@@ -675,7 +692,6 @@ recoup <- function(

     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

@@ -723,17 +739,7 @@ recoup <- function(

             #helperRanges <- helperRanges[keeph]

             #mainRanges <- mainRanges[names(helperRanges)]

             helperRanges <- .subsetGRangesByChrs(helperRanges,chrs)

-            mainRanges <- .subsetGRangesByChrs(mainRanges,chrs)

-            

-            ####################################################################

-            ## There must be an R bug with `lengths` here as although it runs in

-            ## Rcmd, it does not pass package building or vignette kniting...  

-            ## But for the time being it seems that it is not needed as the name

-            ## filtering works

-            #lens <- which(lengths(genomeRanges)==0)

-            #if (length(lens)>0)

-            #    genomeRanges[lens] <- NULL

-            ####################################################################

+			mainRanges <- .subsetGRangesByChrs(mainRanges,chrs)

         }

     }

@@ -745,8 +751,10 @@ recoup <- function(

     # 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")

@@ -812,12 +820,13 @@ recoup <- function(

     # 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)

+        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)

+            profile=TRUE),callParams=callParams,refRanges=refRanges)

     }

     ## Our plot objects

@@ -961,7 +970,8 @@ recoup <- function(

     }

     # Overwrite final object so as to return it

-    recoupObj <- toOutput(input,design,saveParams,recoupPlots,callParams)

+    recoupObj <- toOutput(input,design,saveParams,recoupPlots,callParams,

+        refRanges)

     # Make any plots asked

     if (plotParams$plot) {

@@ -1009,7 +1019,7 @@ applySelectors <- function(ranges,selector,rc=NULL) {

 }

 toOutput <- function(input,design=NULL,saveParams,plotObjs=NULL,

-    callParams=NULL) {

+    callParams=NULL,refRanges=NULL) {

     if (!saveParams$ranges)

         input <- removeData(input,"ranges")

     if (!saveParams$coverage)

@@ -1027,6 +1037,7 @@ toOutput <- function(input,design=NULL,saveParams,plotObjs=NULL,

         data=input,

         design=design,

         plots=plots,

-        callopts=callParams

+        callopts=callParams,

+        refRanges=refRanges

     ))

 }

@@ -216,10 +216,11 @@ recoup <- function(

     # 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" && !.userOrg(genome,localDb)

-        && !(genome %in% getSupportedOrganisms()))

-        stop("When type is \"rnaseq\", only the supported genomes are allowed!")

-    

+    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

@@ -10,6 +10,7 @@ recoup <- function(

     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"),

@@ -58,6 +59,7 @@ recoup <- function(

         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=".",

@@ -163,15 +165,25 @@ recoup <- function(

     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))

-        checkTextArgs("genome",genome,getSupportedOrganisms(),multiarg=FALSE)

-    checkTextArgs("refdb",refdb,getSupportedRefDbs(),multiarg=FALSE)

+    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))

@@ -204,7 +216,8 @@ recoup <- function(

     # 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" && !(genome %in% getSupportedOrganisms()))

+    if (type=="rnaseq" && !.userOrg(genome,localDb)

+        && !(genome %in% getSupportedOrganisms()))

         stop("When type is \"rnaseq\", only the supported genomes are allowed!")

     # annotation must be a list to be fed to buildCustomAnnotation

@@ -291,6 +304,8 @@ recoup <- function(

             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

@@ -335,6 +350,7 @@ recoup <- function(

         callParams <- list(

             region=region,

             type=type,

+            onFlankFail=onFlankFail,

             signal=signal,

             genome=genome,

             version=version,

@@ -366,7 +382,8 @@ recoup <- function(

     # Redefine all final arguments for this call

     region <- callParams$region

     type <- callParams$type

-    singnal <- callParams$signal

+    onFlankFail <- callParams$onFlankFail

+    signal <- callParams$signal

     genome <- callParams$genome

     version <- callParams$version

     refdb <- callParams$refdb

@@ -517,7 +534,7 @@ recoup <- function(

     # 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

+    #hasDesign <- FALSE

     if (!is.null(design)) {

         if (!is.data.frame(design))

             design <- read.delim(design,row.names=1)

@@ -617,7 +634,7 @@ recoup <- function(

     # 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,rc=rc)

+        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

@@ -625,12 +642,32 @@ recoup <- function(

     # 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

-    if (type == "chipseq") # Otherwise, throw error with GRangesList

-        mainRanges <- trim(intermRanges$mainRanges)

-    else if (type == "rnaseq")

-        mainRanges <- intermRanges$mainRanges

-    bamRanges <- trim(intermRanges$bamRanges)

-

+    

+    # 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 thre 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

@@ -648,7 +685,7 @@ recoup <- function(

         #set.seed(preprocessParams$seed)

         refIndex <- sort(sample(length(mainRanges),newSize))

         mainRanges <- mainRanges[refIndex]

-        for (i in 1:length(input)) {

+        for (i in seq_len(length(input))) {

             if (!is.null(input[[i]]$ranges)) {

                 newSize <- round(fraction*length(input[[i]]$ranges))

                 #set.seed(preprocessParams$seed)

@@ -745,13 +782,17 @@ recoup <- function(

     if (signal == "coverage") { # Otherwise, matrix calculate by rpm

         if (binParams$chunking) {

             if (type=="chipseq")

-                binParams$chunks <- split(1:length(genomeRanges),

+                binParams$chunks <- split(seq_len(length(genomeRanges)),

                     as.character(seqnames(genomeRanges)))

             else if (type=="rnaseq")

-                binParams$chunks <- split(1:length(helperRanges),

+                binParams$chunks <- split(seq_len(length(helperRanges)),

                     as.character(seqnames(helperRanges)))

         }

-        input <- profileMatrix(input,flank,binParams,rc)

+        # 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

@@ -756,8 +756,11 @@ recoup <- function(

     # 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.

-    input <- imputeProfile(input,method="simple",rc)

+    # 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)

@@ -415,6 +415,11 @@ recoup <- function(

                 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

         }

         # Annotation case #3: use local database or automatically on-the-fly

         else {

@@ -752,7 +757,7 @@ recoup <- function(

     # 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.

-    input <- imputeProfile(input,rc)

+    input <- imputeProfile(input,method="simple",rc)

     # Perform the k-means clustering if requested and append to design (which

     # has been checked, if we are allowed to do so)

@@ -212,7 +212,7 @@ recoup <- function(

         # members are checked by buildCustomAnnotation if required

         # We only need to check that the gtfFile exists here

         if (!("gtf" %in% names(genome)))

-            stopp("A gtf field must be provided with an existing GTF file ",

+            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))

@@ -479,8 +479,8 @@ recoup <- function(

                     type="gene",version=version,db=localDb,rc=rc),

                     error=function(e) {

                         tryCatch({

-                            gtfFile <- annotation$gtf

-                            metadata <- annotation

+                            gtfFile <- genome$gtf

+                            metadata <- genome

                             metadata$gtf <- NULL

                             geneData <- importCustomAnnotation(gtfFile,

                                 metadata,"gene")

@@ -749,6 +749,11 @@ recoup <- function(

         input <- profileMatrix(input,flank,binParams,rc)

     }

+    # 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.

+    input <- imputeProfile(input,rc)

+    

     # 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)

@@ -1,7 +1,7 @@

 recoup <- function(

     input,

     design=NULL,

-    region=c("genebody","tss","tes","custom"),

+    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",

@@ -126,9 +126,16 @@ recoup <- function(

     #complexHeatmapParams=getDefaultListArgs("complexHeatmapParams"),

     bamParams=NULL,

     onTheFly=FALSE, # Directly from BAM w/o Ranges storing, also N/A with BED,

-    localDbHome=file.path(path.expand("~"),".recoup"),

+    #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

@@ -147,7 +154,7 @@ recoup <- function(

         input <- readConfig(input)

     checkInput(input)

     if (is.null(names(input)))

-        names(input) <- sapply(input,function(x) return(x$id))

+        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()))

@@ -159,7 +166,7 @@ recoup <- function(

     signal <- tolower(signal[1])

     if (!is.null(design) && !is.data.frame(design))

         checkFileArgs("design",design)

-    if (!is.data.frame(genome) && file.exists(genome))

+    if (!is.data.frame(genome) && !is.list(genome) && file.exists(genome))

         checkFileArgs("genome",genome)

     else if (is.character(genome))

         checkTextArgs("genome",genome,getSupportedOrganisms(),multiarg=FALSE)

@@ -200,6 +207,33 @@ recoup <- function(

     if (type=="rnaseq" && !(genome %in% getSupportedOrganisms()))

         stop("When type is \"rnaseq\", only the supported 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)))

+            stopp("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")

@@ -232,12 +266,12 @@ recoup <- function(

     strandedParams <- validateListArgs("strandedParams",strandedParams)

     if (is.null(plotParams$outputBase))

-        plotParams$outputBase <- paste(sapply(input,function(x) return(x$id)),

-            collapse="-")

+        plotParams$outputBase <- paste(vapply(input,function(x) return(x$id),

+            character(1)),collapse="-")

-    if (any(sapply(input,function(x) return(tolower(x$format)=="bed")))

-        && !(preprocessParams$bedGenome %in% c("hg18","hg19","hg38","mm9",

-            "mm10","rn5","dm3","danrer7","pantro4","susscr3")))

+    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.")

@@ -291,8 +325,8 @@ recoup <- function(

             #    prevCallParams$bamParams,bamParams), ## Unused

             onTheFly=if (is.null(thisCall$onTheFly)) prevCallParams$onTheFly 

                 else onTheFly,

-            localDbHome=if (is.null(thisCall$localDbHome)) 

-                prevCallParams$localDbHome else localDbHome,

+            localDb=if (is.null(thisCall$localDb)) 

+                prevCallParams$localDb else localDb,

             rc=if (is.null(thisCall$rc)) prevCallParams$rc else rc,

             customIsBase=NULL # Additional placeholder

         )

@@ -319,7 +353,7 @@ recoup <- function(

             complexHeatmapParams=complexHeatmapParams,

             bamParams=bamParams,

             onTheFly=onTheFly,

-            localDbHome=localDbHome,

+            localDb=localDb,

             rc=rc,

             customIsBase=NULL # Additional placeholder

         )

@@ -350,152 +384,127 @@ recoup <- function(

     complexHeatmapParams <- callParams$complexHeatmapParams

     bamParams <- callParams$bamParams

     onTheFly <- callParams$onTheFly

-    localDbHome <- callParams$localDbHome

+    localDb <- callParams$localDb

     rc <- callParams$rc

     # End of complex parameter storage procedure and parameter recall from a

     # previous call

     ############################################################################

     # Continue with actual work

-    if (!is.data.frame(genome)) {

+    # Annotation case #1: provided strictly 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

+        )

+    }

+    # 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)) {

-            genome <- read.delim(genome) # Must be bed like

+            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

             )

         }

+        # Annotation case #3: use local database or automatically on-the-fly

         else {

-            # Check if local storage has been set

-            storageHome <- file.path(localDbHome,refdb,genome)

-            if (dir.exists(storageHome)) {

-                # Some compatibility with previous annoation stores

-                if ("gene.rda" %in% dir(storageHome)) { # Older version

-                    warning("Older annotation storage detected! Please ",

-                        "rebuild your annotation storage\nusing the ",

-                        "buildAnnotationStore function so as to have the\n",

-                        "ability of versioning genomic coordinate annotations.",

-                        "\nIn the future, older versions may become unusable.",

-                        "\n",immediate.=TRUE)

-                    storageHomeVersion <- storageHome

-                }

-                else { # Newer version

-                    # Decide on version

-                    if (version != "auto") {

-                        storageHomeVersion <- file.path(storageHome,version)

-                        if (!dir.exists(storageHomeVersion)) {

-                            warning("The annotation directory ",

-                                storageHomeVersion,

-                                " does not seem to exist! Have you run ",

-                                "buildAnnotationStorage? Will use newest ",

-                                "existing version.",immediate.=TRUE)

-                            version <- "auto"

-                        }

-                    }

-                    if (version == "auto") {

-                        vers <- suppressWarnings(as.numeric(dir(storageHome)))

-                        if (any(is.na(vers))) {

-                            of <- vers[which(is.na(vers))]

-                            stop("Corrupted annotation storage directory ",

-                                "contents! ->",of,"<- Either remove offending ",

-                                "files/directories or rebuild.")

-                        }

-                        vers <- sort(vers,decreasing=TRUE)

-                        version <- vers[1]

-                        storageHomeVersion <- file.path(storageHome,version)

-                    }

-                }

-                

-                if (type=="chipseq") {

-                    g <- load(file.path(storageHomeVersion,"gene.rda"))

-                    genomeRanges <- gene

-                    helperRanges <- NULL

-                }

-                else if (type=="rnaseq") {

-                    # Load the helper ranges anyway

-                    g <- load(file.path(storageHomeVersion,"gene.rda"))

-                    helperRanges <- gene

-                    if (all(flank==0)) {

-                        g <- load(file.path(storageHomeVersion,

-                            "summarized_exon.rda"))

-                        genomeRanges <- sexon

-                    }

-                    else if (all(flank==500)) {

-                        g <- load(file.path(storageHomeVersion,

-                            "summarized_exon_F500.rda"))

-                        genomeRanges <- flankedSexon

-                    }

-                    else if (all(flank==1000)) {

-                        g <- load(file.path(storageHomeVersion,

-                            "summarized_exon_F1000.rda"))

-                        genomeRanges <- flankedSexon

-                    }

-                    else if (all(flank==2000)) {

-                        g <- load(file.path(storageHomeVersion,

-                            "summarized_exon_F2000.rda"))

-                        genomeRanges <- flankedSexon

-                    }

-                    else if (all(flank==5000)) {

-                        g <- load(file.path(storageHomeVersion,

-                            "summarized_exon_F5000.rda"))

-                        genomeRanges <- flankedSexon

-                    }

-                    else {

-                        warning("When using recoup in RNA-Seq mode, it is ",

-                            "much faster to use a precalculated set of ",

-                            "regions and their flanks (0, 500, 1000, 2000 and ",

-                            "5000 bps supported) and subset this one for a ",

-                            "custom set of genes. Otherwise calculations may ",

-                            "take too long to complete.",immediate.=TRUE)

-                        genomeRanges <- 

-                            getMainRnaRangesOnTheFly(helperRanges,flank,rc=rc)

-                    }

-                }

+            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 { # On the fly

-                message("Getting annotation on the fly for ",genome," from ",

-                    refdb)

-                if (type=="chipseq") {

-                    genome <- getAnnotation(genome,"gene",refdb=refdb,rc=rc)

-                    helperRanges <- NULL

-                }

-                else if (type=="rnaseq") {

-                    helper <- getAnnotation(genome,"gene",refdb=refdb,rc=rc)

-                    helperRanges <- makeGRangesFromDataFrame(

-                        df=helper,

-                        keep.extra.columns=TRUE,

-                        seqnames.field="chromosome"

-                    )

-                    names(helperRanges) <- as.character(helperRanges$gene_id)

-                    genome <- getAnnotation(genome,"exon",refdb=refdb,rc=rc)

-                    annGr <- makeGRangesFromDataFrame(

-                        df=genome,

-                        keep.extra.columns=TRUE,

-                        seqnames.field="chromosome"

-                    )

-                    message("Merging exons")

-                    annGr <- reduceExons(annGr,rc=rc)

-                    names(annGr) <- as.character(annGr$exon_id)

-                    genomeRanges <- split(annGr,annGr$gene_id)

-                }

+            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="")