R/computePairwiseEnrichments.R

In jsalojar/TFbindR: Identification and Counting of the Occurence of Motifs in a Set of Promoter Sequences. Motif Enrichment Analysis

#' Pairwise motif enrichment analysis across a database of motifs
#' 
#' Test for pairwise motif enrichment in a set of genes of interest. The
#' function performs a fisher exact test to calculte the enrichment of each
#' motif combination for the set of selected genes, from the background genes
#' in promoter.DB.
#' 
#' %% ~~ If necessary, more details than the description above ~~
#' 
#' @param Genes list of the genes of interest
#' @param promoterDB promoterDB object
#' @return An R object with 3 elements. Each of the elements is a data frame
#' with motif combination on the rows and the enrichment statistics on the
#' columns. Each row corresponds to one-sided \code{\link[stats]{fisher.test}}
#' for enirchment of a motif pair in the given gene set, based on combined
#' occurence of the motifs in the plus, minus or both strands combined. The
#' columns of the data frame are:
#' 
#' p.value - the p.value of the Fisher test,
#' 
#' odds.ratio - the estimate of the odds ration from the Fisher test,
#' 
#' genes - subset of the genes that are enirched with the motif and
#' 
#' p.adj.BH - the multiple test correction of the p.value with the "BH"
#' adjustment.
#' @note %% ~~further notes~~
#' @author %% ~~who you are~~
#' @seealso %% ~~objects to See Also as \code{\link{help}}, ~~~
#' @references %% ~put references to the literature/web site here ~
#' @examples
#' 
#' # Read an example gene set. 
#' genes=read.table(system.file("extdata","test_geneSet.txt",package="TFbindR"),stringsAsFactors=F)
#' 
#' # Run the function with the promoterDB generated by the function makePromoterDBfromFIMO.
#' pairwise_enrichment=computePairwiseEnrichments(as.character(unlist(genes)), "TAIR10_500bp_upstream.fimo.db.rda") 
#' 
#' @export computePairwiseEnrichments
computePairwiseEnrichments <- function(Genes, promoterDB = "promoter.DB") {
    if (promoterDB!="TAIR10_500bp_upstream.fimo.db.rda")
    	load(promoterDB)
    Allgenes = rownames(Strands[[1]])
    
    # test if all genes exist in promoter DB
    t1 = which(!(Genes %in% Allgenes))
    if (length(t1) > 0) {
        message("Gene(s) ", paste(Genes[t1], collapse = ""), " not found from promoter DB. Dropping from analysis.")
        Genes = Genes[Genes %in% Allgenes]
    }
    
    # test pairwise interactions
    Ori = c("(-)", "(+)", "(+/-)")
    res = list()
    rescnt = 1
    for (k in 1:3) {
        D1 = Strands[[k]]
        Nmotif1 = ncol(D1)
        for (r in 1:3) {
            D2 = Strands[[r]]
            Nmotif2 = ncol(D2)
            ftest = data.frame(motif1 = rep("", Nmotif1 * Nmotif2), motif2 = rep("", Nmotif1 * 
                Nmotif2), p.value = rep(NA, Nmotif1 * Nmotif2), odds.ratio = rep(NA, Nmotif1 * 
                Nmotif2), p.adj.BH = rep(NA, Nmotif1 * Nmotif2), genes = rep("", Nmotif1 * 
                Nmotif2), stringsAsFactors = F)
            cnt = 1
            motifs1 = colnames(D1)
            for (j in motifs1) {
                motifs2 = setdiff(colnames(D2), motifs1[1:which(motifs1 == j)])
                for (i in motifs2) {
                  if (i != j) 
                    paired.motif = (D1[, j] & D2[, i]) else paired.motif = (D1[, j] > 1)
                  if (sum(paired.motif[Genes] > 0) > 0) {
                    MM = fisher.test(Allgenes %in% Genes, paired.motif)
                    ftest$motif1[cnt] = j
                    ftest$motif2[cnt] = i
                    ftest$p.value[cnt] = MM$p.value
                    ftest$odds.ratio[cnt] = MM$estimate
                    ftest$genes[cnt] = paste(names(which(paired.motif[Genes] > 0)), collapse = ";")
                    cnt = cnt + 1
                  }
                }
            }
            ftest = ftest[which(!is.na(ftest$p.value)), ]
            ftest$p.adj.BH = p.adjust(ftest$p.value, method = "BH")
            res[[rescnt]] = ftest
            names(res)[rescnt] = paste("Enrichment", Ori[k], Ori[r], sep = "")
            rescnt = rescnt + 1
        }
    }
    return(res)
}