R/morifEnrichment.R

In LouisKwok-PICB/motifscanr: Tools for Scanning TF motif on DNA regions

motif_enrichment_helper <- function(motif_sites, motif_sites_control){
  motif_sites <- as.matrix(motif_sites)
  motif_sites_control <- as.matrix(motif_sites_control)
  if(any(colnames(motif_sites) != colnames(motif_sites_control))){
    stop('the TF motif is not consistent between position weight matrixs list and input matrix!')
  }
  n_motifs <- ncol(motif_sites)
  enrichresult <- t(as.matrix(as.data.frame(lapply(seq_len(n_motifs), function(n){
    sites <- motif_sites[,n]
    sites_control <- motif_sites_control[,n]
    n_input_total <- as.numeric(length(sites))
    n_control_total <- as.numeric(length(sites_control))
    n_input <- as.numeric(sum(sites > 0))
    n_control <- as.numeric(sum(sites_control > 0))
    if((n_input_total > 0) & (n_control > 0)){
      fold_change <- n_input * n_control_total / n_control / n_input_total
    }else{
      fold_change <- NA
    }
    table <- matrix(c(n_input, n_input_total - n_input, n_control, n_control_total - n_control),
                    nrow = 2, ncol = 2)
    enriched_out <- fisher.test(table, alternative = "greater")
    depleted_out <- fisher.test(table, alternative = "less")
    p_enriched <- enriched_out[['p.value']]
    p_depleted <- depleted_out[['p.value']]
    p_corrected = min(min(p_enriched, p_depleted) * n_motifs, 1)
    return(c(colnames(motif_sites)[n], n_input, n_control, fold_change,
             p_enriched, p_depleted, p_corrected))
  }))))
  colnames(enrichresult) <- c('Motif', 'Num_input_regions', 'Num_control_regions', 'Fold_change',
                              'Enriched_P_value', 'Depleted_P_value', 'Corrected_P_value')
  rownames(enrichresult) <- 1:n_motifs
  enrichresult_out <- data.frame(Motif = enrichresult[,'Motif'],
                                 Num_input_regions = as.integer(enrichresult[,'Num_input_regions']),
                                 Num_control_regions = as.integer(enrichresult[,'Num_control_regions']),
                                 Fold_change = as.numeric(enrichresult[,'Fold_change']),
                                 Enriched_P_value = as.numeric(enrichresult[,'Enriched_P_value']),
                                 Depleted_P_value = as.numeric(enrichresult[,'Depleted_P_value']),
                                 Corrected_P_value = as.numeric(enrichresult[,'Corrected_P_value']))
  return(enrichresult_out)
}

#' motifEnrichment
#'
#' @description function to find enrichment motif between two motifscan result
#'
#' @param InputMotif should be generated by
#'  \code{\link[motifscanR]{motifScan}} with input peaks
#' @param ControlMotif should be generated by
#'  \code{\link[motifscanR]{motifScan}}, with control peaks and
#'  the same parameter of input peaks
#'
#' @return return Motif Enrichment data.frame of InputMotif by ControlMotif
#'
#' @export
#'
#' @examples
#' example_motifs <- getJasparMotifs(species = "Homo sapiens",
#'                                   collection = "CORE")
#' # Make a set of input regions
#' Input <- GenomicRanges::GRanges(seqnames = c("chr1","chr2","chr2"),
#'                                 ranges = IRanges::IRanges(start = c(76585873,42772928,
#'                                                                     100183786),
#'                                                           width = 500))
#' # Make a set of control regions
#' Control <- GenomicRanges::GRanges(seqnames = c("chr1","chr3","chr5"),
#'                                   ranges = IRanges::IRanges(start = c(453123,6524593,
#'                                                                       100184233),
#'                                                             width = 500))
#' # Scan motif for example motifs
#' motif_ix_input <- motifScan(example_motifs, Input, genome = "BSgenome.Hsapiens.UCSC.hg19")
#' motif_ix_control <- motifScan(example_motifs, Control, genome = "BSgenome.Hsapiens.UCSC.hg19")
#'
#' # Find Enrichment motif of input by control
#' Enrichment_result <- motifEnrichment(motif_ix_input, motif_ix_control)
#'
setGeneric("motifEnrichment",
           function(InputMotif, ControlMotif, ...) standardGeneric("motifEnrichment"))


#' @describeIn motifEnrichment RangedSummarizedExperiment/RangedSummarizedExperiment
#' @export
setMethod("motifEnrichment", signature(InputMotif = "RangedSummarizedExperiment",
                                       ControlMotif = "RangedSummarizedExperiment"),
          function(InputMotif,
                   ControlMotif) {
            motif_enrichment_helper(InputMotif, ControlMotif)
          })