From 783829e1cbf44b3ad92f9abf82bc7700768783fd Mon Sep 17 00:00:00 2001 From: thomaz Date: Sun, 26 Jul 2026 12:46:17 +0200 Subject: [PATCH] add test_enrichment --- .github/workflows/pkgdown.yaml | 2 +- .github/workflows/test.yml | 2 +- NAMESPACE | 5 ++ R/weave_coverage.R | 139 +++++++++++++++++++-------------- man/weave_coverage.Rd | 21 +++-- 5 files changed, 102 insertions(+), 67 deletions(-) diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index c0e24cd..98ea6b9 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -1,5 +1,5 @@ name: build -'on': +on: push: branches: - devel diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index a670a81..a2ccd6b 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -1,5 +1,5 @@ name: R-BiocCheck -'on': +on: pull_request: branches: - devel diff --git a/NAMESPACE b/NAMESPACE index 0084217..c478ace 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -29,6 +29,7 @@ export(randomLinkMap) export(randomMultiFactor) export(read_adjacency_list) export(select_path) +export(test_enrichment) export(trade_posts) export(weave) export(weave_apply) @@ -36,6 +37,7 @@ 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) diff --git a/R/weave_coverage.R b/R/weave_coverage.R index 9a5309d..7554676 100644 --- a/R/weave_coverage.R +++ b/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) } diff --git a/man/weave_coverage.Rd b/man/weave_coverage.Rd index 9af1655..aead8ee 100644 --- a/man/weave_coverage.Rd +++ b/man/weave_coverage.Rd @@ -2,19 +2,22 @@ % Please edit documentation in R/weave_coverage.R \name{weave_coverage} \alias{weave_coverage} +\alias{test_enrichment} \title{Perform enrichtment analysis from a weave} \usage{ weave_coverage( 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" ) + +test_enrichment(x, log_base = 2L) } \arguments{ -\item{x}{a \code{MultiFactor}} +\item{x}{A LinkMap with coverage metadata from \code{weave_coverage()}.} \item{.path}{Either a \code{formula} or a \verb{character vector} of length 2 with the names of the desired combination of feature types.} @@ -24,13 +27,15 @@ found within the first element of \code{.path}. Alternatively, a \code{data.fram with the same information. (Also see \code{.data_column} argument).} \item{metric}{\verb{Character scalar}. -One of \code{'count'}, \code{'coverage'}, \code{'complete'}.} +One or more of \code{'set_count'}, \code{'set_size'}, \code{'coverage'}, \code{'complete'}.} \item{out.format}{\verb{Character scalar}. One of \code{'LinkMap'}, \code{'matrix'}.} \item{.data_column}{\verb{Character scalar}. if \code{.data} is a table, where to find feature IDs} + +\item{log_base}{\code{Integer}. Base of logarithm for log-fold (default: 2).} } \value{ a \code{LinkMap} or \code{matrix} with the desired coverage information in @@ -49,16 +54,18 @@ scores # 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) }