R/GeneUpDownStream.R

In genomaths/MethylIT.utils: MethylIT utility

#' @rdname GeneUpDownStream
#' @title Get Genes plus Up and Down Stream Regions
#' @description Given a genes region or genomic region (GR), this function
#'     yields the GR plus the especified amount of DNA bases upstream and
#'     downstream the GR. 
#' @details Users can select whether to request only upstream, only downstream,
#'      or both, upstream and downstream. Please notice that for a gene on the
#'      negative strand, 'the start of the gene' corresponds to the 'end' of the
#'      gene in the GRanges object and the 'end of the gene' correspond to the
#'      'start' of the gene in the GRanges object.
#' @param GR A \code{\link[GenomicRanges]{GRanges-class}} object containing the
#'     ranges of the genes or genomic regions to be extended upstream/downstream
#' @param upstream Integer (Default: 0). The amount of DNA bases (bps) upstream
#'     of the GR.
#' @param downstream Integer (Default: 0). The amount of DNA bases (bps)
#'     downstream of the GR.
#' @param extend Integer (Default: NULL). If upstream == downstream, then 
#'     simply you may use extend. 
#' @param fix A string with one of the three possible values: "start", "end" or
#'     "center" to denoteg what to use as an anchor for each element in GR.
#' @param onlyUP Logic (Default: FALSE). If TRUE returns the region upstream the
#'     GR.
#' @param onlyDown Logic (Default: FALSE). If TRUE returns the region downstream
#'     the GR.
#' @examples
#' starts = c(65419, 450703, 923928, 944204)
#' ends = c(71585, 451697, 944581, 959309)
#' chrs = c(rep("chr1", 2), rep("chr2", 2))
#' gr = makeGRangesFromDataFrame(
#'         data.frame(seqnames = chrs, start = starts, end = ends,
#'                 strand = c("+", "-", "+", "-"),
#'                 genes = c("A", "B", "C", "D")), keep.extra.columns = TRUE)
#'
#' gr1 = GeneUpDownStream(GR = gr, upstream = 2000, downstream = 1000)

#' @importFrom GenomeInfoDb seqnames
#' @importFrom BiocGenerics start end
#' @importFrom S4Vectors mcols
#' @export

GeneUpDownStream <- function(GR, upstream=0, downstream=0, extend = NULL,
                               fix = NULL, onlyUP=FALSE, onlyDown=FALSE) {
   if (!is.null(extend)) {
       if (is.numeric(extend)) extend <- as.integer(extend)
       else 
           stop("*** 'extend' must be an integer or coercible to an integer")
       upstream = extend
       downstream = extend
   }
   
   strands <- as.character( strand( GR ) )
   if (!is.null(fix)) {
       if (fix == "start") {
          starts <- start( GR )
          ends <- starts
       }
       if (fix == "end") {
          ends <- end( GR )
          starts <- ends
       }
       if (fix == "center") {
          starts <- start( GR )
          ends <- end( GR )
          starts <- round((ends - starts)/2)
          ends <- starts
       } 
   } else {
       starts <- start( GR )
       ends <- end( GR )
   }
   
   chrs <- as.character(seqnames( GR ))
   
   if (upstream > 0 && !onlyUP && !onlyDown) {
      
       ind <- which( strands == "+" )
       starts[ ind ] <- starts[ ind ] - upstream

       ind <- which( strands == "-" )
       ends[ ind ] <- ends[ ind ] + upstream

       GR.up <- GRanges( seqnames = chrs,
                       ranges = IRanges( start = starts, end = ends ),
                       strand = strands )
       mcols(GR.up) <- mcols(GR)
       GR <- GR.up; rm(GR.up); gc()
   }

   if (downstream > 0 && !onlyUP && !onlyDown) {

       ind <- which( strands == "+" )
       ends[ ind ] <- ends[ ind ] + downstream

       ind <- which( strands == "-" )
       starts[ ind ] <- starts[ ind ] - downstream

       GR.down <- GRanges( seqnames = chrs,
                       ranges = IRanges( start = starts, end = ends ),
                       strand = strands )
       mcols(GR.down) <- mcols(GR)
       GR <- GR.down; rm(GR.down); gc()
   }

   if (onlyUP) {
       if (onlyDown) stop("* If onlyUP is TRUE, then onlyDown must be FALSE")

       ind <- which( strands == "+" )

       ends[ ind ] <- starts[ ind ] - 1
       starts[ ind ] <- starts[ ind ] - upstream

       ind <- which( strands == "-" )
       starts[ ind ] <- ends[ ind ] + 1
       ends[ ind ] <- ends[ ind ] + upstream

       GR.up <- GRanges( seqnames = chrs,
                        ranges = IRanges( start = starts, end = ends ),
                        strand = strands )
       mcols(GR.up) <- mcols(GR)
       GR <- GR.up; rm(GR.up); gc()
   }

   if (onlyDown) {
       if (onlyUP) stop("* If onlyDown is TRUE, then onlyUP must be FALSE")

       ind <- which( strands == "+" )
       starts[ ind ] <- ends[ ind ] + 1
       ends[ ind ] <- ends[ ind ] + downstream

       ind <- which( strands == "-" )
       ends[ ind ] <- starts[ ind ] - 1
       starts[ ind ] <- starts[ ind ] - downstream

       GR.down <- GRanges(seqnames = chrs,
                       ranges = IRanges( start = starts, end = ends ),
                       strand = strands )
       mcols(GR.down) <- mcols(GR)
       GR <- GR.down; rm(GR.down); gc()
   }
   return(GR)
}