R/ChIP.enhancer.R

In robinweide/RHWlib: Simple functions to make life... simpler.

#' ChIP.enhancer
#'
#' Get a GRanges with enhancer-positions
#'
#' @param H3K27ac A BED-df with H3k27ac-peak positions.
#' @param activePeak A BED-df with peak positions of an active enhancer mark (e.g. H3Kme1).
#' @param ignorePeak A BED-df with peak positions of an mark not associated with enhancers (e.g. H3Kme3).
#' @param padding Extend the peaks.
#' @return A genomicRanges-object
#' @export
ChIP.enhancer <- function(H3K27ac, activePeak = NULL, ignorePeak = NULL, padding = 100, verbose = F){
  require(RHWlib)
  require(GenomicRanges)

  AC <- df2gr(H3K27ac, keepMeta = F)
  ME3 <- NULL
  ME1 <- NULL
  activeEnhancersFiltered <- NULL
  activeEnhancers <- NULL
  enhancersFiltered <- NULL

  # If no active and ignore, warn that the fuction will just output H3K27ac!
  if(is.null(activePeak)){
    if(is.null(ignorePeak)){
      warning("No peaks to ignore or to denote active. Outputting the H3K27ac peaks!")
    } else {
      ME3 <- df2gr(ignorePeak, keepMeta = F)
      warning("No peaks given to denote active enhancers.")
    }
  } else { # active is given
    ME1 <- df2gr(activePeak, keepMeta = F)
    if(is.null(ignorePeak)){
      warning("No peaks to ignore. Enhancers could also potentially overlap promoters!")
    } else {
      ME3 <- df2gr(ignorePeak, keepMeta = F)
      if(verbose){message("Both active sites and non-enhancer peaks given.\nLet's get to work!")}
    }
  }

  # # add padding
  # AC <- GenomicRanges::flank(AC,width = padding, both = T)
  # if(!is.null(ME1)){
  #   ME1 <- GenomicRanges::flank(ME1,width = padding, both = T)
  # }
  # if(!is.null(ME3)){
  #   ME3 <- GenomicRanges::flank(ME3,width = padding, both = T)
  # }

  # Get Active enhancers
  if(!is.null(ME1)){ # active is there

    activeEnhancersHits <- suppressWarnings(findOverlaps(AC,ME1, maxgap = padding))
    activeEnhancers <- AC[activeEnhancersHits@from]


    if(!is.null(ME3)){ # ignore is also there

      EPhits <-  suppressWarnings(findOverlaps(activeEnhancers,ME3, maxgap = padding))
      activeEnhancersFiltered <- activeEnhancers[ - EPhits@from]

    }

  } else {


    if(!is.null(ME3)){ # ignore is there

      ACME3hits <-  suppressWarnings(findOverlaps(AC,ME3, maxgap = padding))
      enhancersFiltered <- AC[ - ACME3hits@from]

    }
  }


  # eerst proberen beide
  # dan degene met ignore
  # dan degene met active
  # als laatste ac

  if(!is.null(activeEnhancersFiltered)){
    return(unique(activeEnhancersFiltered))
    } else if(!is.null(enhancersFiltered)){
      return(unique(enhancersFiltered))
    } else if(!is.null(activeEnhancers)){
      return(unique(activeEnhancers))
    } else {
      return(unique(AC))
    }
  # I have to test this output... IGV!
}