add test_enrichment

.github/workflows/pkgdown.yaml

@@ -1,5 +1,5 @@
 name: build
-'on':
+on:
   push:
     branches:
     - devel

.github/workflows/test.yml

@@ -1,5 +1,5 @@
 name: R-BiocCheck
-'on':
+on:
   pull_request:
     branches:
     - devel

NAMESPACE

@@ -29,13 +29,15 @@ export(randomLinkMap)
 export(randomMultiFactor)
 export(read_adjacency_list)
 export(select_path)
+export(test_enrichment)
 export(trade_posts)
 export(weave)
 export(weave_apply)
 export(weave_coverage)
 export(weave_to_tbl)
 if (getRversion() < "4.3.0") importFrom("S7", "@")
 import(S7)
+importFrom(Matrix,Matrix)
 importFrom(Matrix,colSums)
 importFrom(Matrix,crossprod)
 importFrom(Matrix,rowSums)
@@ -46,6 +48,7 @@ importFrom(S7,"S7_data<-")
 importFrom(S7,"prop<-")
 importFrom(S7,S7_data)
 importFrom(S7,S7_dispatch)
+importFrom(S7,S7_inherits)
 importFrom(S7,class_character)
 importFrom(S7,new_class)
 importFrom(S7,new_generic)
@@ -73,6 +76,8 @@ importFrom(rlang,dots_list)
 importFrom(rlang,f_lhs)
 importFrom(rlang,f_rhs)
 importFrom(rlang,is_missing)
+importFrom(stats,p.adjust)
+importFrom(stats,phyper)
 importFrom(stats,reformulate)
 importFrom(stats,terms)
 importFrom(utils,count.fields)

R/weave_coverage.R

@@ -6,7 +6,7 @@
 #'     found within the first element of `.path`. Alternatively, a `data.frame`
 #'     with the same information. (Also see `.data_column` argument).
 #' @param metric `Character scalar`.
-#'     One of `'count'`, `'coverage'`, `'complete'`.
+#'     One or more of `'set_count'`, `'set_size'`, `'coverage'`, `'complete'`.
 #' @param out.format `Character scalar`.
 #'     One of `'LinkMap'`, `'matrix'`.
 #' @param .data_column `Character scalar`. if `.data` is a table, where to find
@@ -24,25 +24,27 @@
 #'
 #' # Now enrich input
 #' result <- weave_coverage(
-#'     x = scores, .path = card ~ suit, .data = drawn, metric = "count"
+#' x = scores, .path = card ~ suit, .data = drawn, metric = "set_count"
 #' )
 #'
 #' result
 #'
-#' # Now let's spike a hand
+#' # Now let's spike our hand with a royal straight flush
 #' cheat <- draw_cards()[c(10, 11, 12, 13, 1)]
 #' cheat
 #'
-#' .path = c("card", "rank", "straight")
-#' weave_coverage(scores, c("card", "rank", "straight"), .data = cheat)
+#' .path = card ~ straight
+#' result <- weave_coverage(scores, card ~ straight, .data = cheat)
+#' result
+#' test_enrichment(result)
 #'
 weave_coverage <- function(
         x, .path, .data = NULL,
-        metric = c("count", "size", "coverage", "complete"),
+        metric = c("set_count", "set_size", "coverage", "complete"),
         out.format = c("LinkMap", "matrix"), .data_column = "row.names"
 ) {
     metric <- match.arg(
-        metric, c("count", "size", "coverage", "complete"),
+        metric, c("set_count", "set_size", "coverage", "complete"),
         several.ok = TRUE
     )
     out.format <- match.arg(out.format, c("LinkMap", "matrix"))
@@ -123,28 +125,54 @@ weave_coverage <- function(
     data.frame(q, m, n, k)
 }
 
-#' @importFrom stats reformulate
+#' @importFrom S7 S7_inherits
+#' @importFrom stats p.adjust
+#' @param x A LinkMap with coverage metadata from `weave_coverage()`.
+#' @param log_base `Integer`. Base of logarithm for log-fold (default: 2).
+#' @export
+#' @rdname weave_coverage
+#' @name test_enrichment
 #'
-.char_contingency_params <- function(x, observed, set) {
-    stopifnot("'observed' must be a named list" = length(names(observed)) == 1L)
-
-    shared <- names(observed)
-    to <- set
+test_enrichment <- function(x, log_base = 2L) {
     stopifnot(
-        "Both set and observed must be found in 'x'" =
-            all(c(shared, set) %in% colnames(x))
+        "'x' must be a LinkMap." = S7::S7_inherits(x, LinkMap),
+        "No coverage data found in metadata. Run `weave_coverage()` first." =
+            all(c("observed", "set_count", "set_size") %in% colnames(x@metadata))
     )
-    full <- weave(x, stats::reformulate(to, shared))
-    hits <- full[full[[shared]] %in% observed[[1L]] ,]
+    i <- stats::reformulate(colnames(x)[[2L]])
+    mm <- as.data.frame(x)
 
-    q <- tapply(X = hits, INDEX = stats::reformulate(to), FUN = NROW)
-    m <- tapply(X = full, INDEX = stats::reformulate(to), FUN = NROW)
-    n <- nlevels(full)[[to]] - m
-    k <- length( unique(hits[[1L]]) )
+    # q = actual hits, m = possible hits, n = possible misses, k = n_draws
+    q <- tapply(mm, i, function(xx) unique(xx[["set_count"]]))
+    m <- tapply(mm, i, function(xx) unique(xx[["set_size"]]))
+    n <- nlevels(mm[[1L]]) - unlist(m)
+    k <- length(unique(mm[mm[["observed"]],1L]))
 
-    data.frame(q, m, n, k)
+    cont_df <- data.frame(q, m, n, k)
+
+    expected  <- k / (m + n)
+    observed  <- q / m
+    fold      <- observed / expected
+    log_fold  <- log(fold, log_base)
+    p.value   <- apply(cont_df, 1L, .apply_phyper)
+    p.adj     <- stats::p.adjust(p.value, method = "BH")
+
+    res <- data.frame(expected, observed, fold, log_fold, p.value, p.adj)
+
+    # Some residual code
+    # cont_list <- apply(cont_df, 1L, .as_cont_matrix, simplify = FALSE)
+    #
+    # lapply(cont_list, fisher.test)
+    # lapply(cont_list, chisq.test)
+
+    return(res)
 }
 
+#' @importFrom stats phyper
+.apply_phyper <- function(x) phyper(
+    x[[1L]] -1L, x[[2L]], x[[3L]], x[[4L]], lower.tail = FALSE
+    )
+
 .as_cont_matrix <- function(x) {
     q <- x[1L]
     m <- x[2L]
@@ -157,25 +185,21 @@ weave_coverage <- function(
 
 
 
-
-#' @importFrom Matrix crossprod colSums t
-#' @importFrom S7 S7_data
-#'
 #' @noRd
 #' @examples
 #' x <- randomMultiFactor()[seq_len(2L)]
 #' at <- c("a", "b", "c")
 #' x.cov <- .weave_coverage_three(x, all_terms = at, metric = "coverage")
-#' x.cnt <- .weave_coverage_three(x, all_terms = at, metric = "count")
+#' x.cnt <- .weave_coverage_three(x, all_terms = at, metric = "set_count")
 #' x.cpt <- .weave_coverage_three(x, all_terms = at, metric = "complete")
 #'
 #' x.cov@metadata
 #' x.cnt@metadata
 #' x.cpt@metadata
 #'
-.weave_coverage_three <- function(
-        x, all_terms, .data, metric
-) {
+#' @importFrom Matrix Matrix crossprod colSums t
+#'
+.weave_coverage_three <- function(x, all_terms, .data, metric ) {
     # All steps in order
     from <- all_terms[[1L]]
     shared <- all_terms[[2L]]
@@ -188,63 +212,62 @@ weave_coverage <- function(
     shared2from <- as.matrix(seen, terms = c(shared, from))
     shared2to   <- as.matrix(full, terms = c(shared, to))
 
-    if( length(.data) ) {
-        shared2from[, !colnames(shared2from) %in% .data] <- FALSE
-        }
     # Link observed features to sets
-    res <- Matrix::crossprod( shared2to, shared2from != 0L )
-    res <- Matrix::t(Matrix::Matrix( res, sparse = TRUE ))
-
-    obs_set <- Matrix::colSums(res)
-
-    tot_set <- pmax(Matrix::colSums(shared2to), 1L)
-
-    metadata <- data.frame(
-        count    = obs_set,
-        size     = tot_set,
-        coverage = obs_set / tot_set,
-        complete = obs_set == tot_set
+    bg_mat <- Matrix::t(
+        Matrix::Matrix(
+            Matrix::crossprod( shared2to, shared2from != 0L ),
+            sparse = TRUE )
     )
 
-    res <- .res_weave_matrix_to_LinkMap(res, levels(x)[c(from, to)])
-    metadata <- metadata[as.integer(res[[to]]), metric]
-    res <- cbind.data.frame(res, metadata)
-
-    row.names(res) <- NULL
+    bg <- .res_weave_matrix_to_LinkMap(bg_mat, levels(x)[c(from, to)])
 
+    if( length(.data) ) {
+        obs <- bg[bg[[from]] %in% .data , ]
+    } else {
+        obs <- bg
+    }
+    res <- .weave_coverage_cont_table(bg, obs, from, to, metric)
     return(res)
+
 }
 
 #' @importFrom Matrix colSums
 #' @importFrom stats reformulate
+#' @importFrom S7 S7_data
 #'
 .weave_coverage_two <- function(
         x, all_terms, .data, metric
 ) {
     set_unit <- all_terms[[1L]]
     set_full <- all_terms[[2L]]
-    x.lm <-x[[1L]]
+    x.lm <- x[[1L]]
     bg <- `class<-`(S7::S7_data(x.lm), "data.frame")
     if( length(.data) ) {
         obs <- bg[bg[[set_unit]] %in% .data , ]
     } else {
         obs <- bg
     }
+    res <- .weave_coverage_cont_table(bg, obs, set_unit, set_full, metric)
+
+    return(res)
+}
+
+.weave_coverage_cont_table <- function(bg, obs, set_unit, set_full, metric) {
     tot_set <- pmax.int(
         c(tapply(bg, INDEX = stats::reformulate(set_full), FUN = NROW)), 1L
-        )
+    )
     obs_set <- c(tapply(obs, INDEX = stats::reformulate(set_full), FUN = NROW))
     val <- data.frame(
-        count    = obs_set,
-        size     = tot_set,
-        coverage = obs_set / tot_set,
-        complete =  obs_set == tot_set
+        set_count = obs_set,
+        set_size  = tot_set,
+        coverage  = obs_set / tot_set,
+        complete  = obs_set == tot_set
     )[metric]
 
     metadata <- val[match(bg[[set_full]], row.names(val)), , drop = FALSE]
-    res <- cbind.data.frame(x.lm, metadata)
+    observed <- bg[[1L]] %in% obs[[1L]]
+    res <- cbind.data.frame(bg, observed, metadata)
 
     row.names(res) <- NULL
-
     return(res)
 }