diff --git a/NAMESPACE b/NAMESPACE
index 345b48d..de4276e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -27,8 +27,10 @@ export(subgroup_apply)
 export(subgroup_to_tbl)
 export(trade_posts)
 export(weave)
+export(weave_coverage)
 if (getRversion() < "4.3.0") importFrom("S7", "@")
 import(S7)
+importFrom(Matrix,Matrix)
 importFrom(Matrix,colSums)
 importFrom(Matrix,crossprod)
 importFrom(Matrix,rowSums)
@@ -36,13 +38,13 @@ importFrom(Matrix,sparseMatrix)
 importFrom(Matrix,t)
 importFrom(Matrix,which)
 importFrom(S7,"prop<-")
+importFrom(S7,S7_data)
 importFrom(S7,S7_dispatch)
 importFrom(S7,new_class)
 importFrom(S7,new_property)
 importFrom(forcats,lvls_expand)
 importFrom(forcats,lvls_union)
 importFrom(generics,augment)
-importFrom(generics,tidy)
 importFrom(igraph,"edge_attr<-")
 importFrom(igraph,E)
 importFrom(igraph,V)
@@ -61,7 +63,6 @@ importFrom(rlang,dots_list)
 importFrom(rlang,f_lhs)
 importFrom(rlang,f_rhs)
 importFrom(rlang,is_missing)
-importFrom(stats,fisher.test)
 importFrom(stats,reformulate)
 importFrom(utils,count.fields)
 importFrom(utils,data)
diff --git a/R/subgroup_apply.R b/R/subgroup_apply.R
index 9a67743..93cd527 100644
--- a/R/subgroup_apply.R
+++ b/R/subgroup_apply.R
@@ -49,8 +49,9 @@ subgroup_apply <- function( X, LINK, BY, FUN = NULL, ..., INDEX = "row.names" )
 #' @name .index_tbl_by
 #' @rdname index_tbl
 #' @description
-#' `.index_tbl_by()` finds row indices for input table `x`, based on a
-#' user-defined path through a MultiFactor (`link` arg).
+#' `.index_tbl()` is a utility function that flexibly takes a table as input and
+#' returns a special `LinkMap` object that has (1) row names and (2)
+#' corresponding row indices of that table as columns.
 #'
 #' @param type `Character scalar` Specifies name of the feature type, will be
 #'     used as name of the first column in output.
@@ -75,9 +76,8 @@ subgroup_apply <- function( X, LINK, BY, FUN = NULL, ..., INDEX = "row.names" )
 #' @name .index_tbl_by
 #' @rdname index_tbl
 #' @description
-#' `.index_tbl()` is a utility function that flexibly takes a table as input and
-#' returns a special `LinkMap` object that has (1) row names and (2)
-#' corresponding row indices of that table as columns.
+#' `.index_tbl_by()` finds row indices for input table `x`, based on a
+#' user-defined path through a MultiFactor (`link` arg).
 #'
 #' @param link a `MultiFactor` object.
 #' @param .path either a `formula` or a `character vector`` of length 2 with the
@@ -107,5 +107,21 @@ subgroup_apply <- function( X, LINK, BY, FUN = NULL, ..., INDEX = "row.names" )
     as.integer(x[[2L]]), x[[1L]], drop
     )
 
-
+#' Take the coverage .data arg and return a character vector.
+#' @noRd
+#'
+.data_coverage_to_vector <- function(.data, .i = "row.names") {
+    if( !length(.data) ) {
+        res <- NULL
+    } else if( is.character(.data) ) {
+        res <- .data
+    } else {
+        if( .i == "row.names") {
+            res <- row.names(.data)
+        } else {
+            res <- .data[[.i]]
+        }
+    }
+    return(res)
+}
 
diff --git a/R/weave-coverage.R b/R/weave-coverage.R
index a20dde0..6933ed2 100644
--- a/R/weave-coverage.R
+++ b/R/weave-coverage.R
@@ -1,17 +1,19 @@
 #' Perform enrichtment analysis from a weave
-#' @param x `MultiFactor` with two `LinkMap` objects, the first one representing
-#'     "seen" features and second the "full" set.
-#' @param observed `list` with one `named character vector`. Name indicated
-#'     feature type, content indicates observed features of said type.
-#' @param set `character` name of the data type to be used as set.
-#' @param alternative `Character scalar` indicates the alternative hypothesis
-#'     and must be one of "two.sided", "greater" or "less".
-#' @param raw `Boolean scalar`. Whether to return the 'untidy' list or call
-#'     `tidy` on the result (default).
-#' @returns a data.frame containing enrichment ratios and p-value following a
-#'     hypergeometric test. see ?phyper
-#' @importFrom stats fisher.test
-#' @importFrom generics tidy
+#' @param x a `MultiFactor`
+#' @param .path Either a `formula` or a `character vector` of length 2 with the
+#'     names of the desired combination of feature types.
+#' @param .data Optional `Character vector`. Lists observed features from the
+#'     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'`.
+#' @param out.format `Character scalar`.
+#'     One of `'LinkMap'`, `'matrix'`.
+#' @param .data_column `Character scalar`. if `.data` is a table, where to find
+#'     feature IDs
+#' @returns a `LinkMap` or `matrix` with the desired coverage information in
+#'     the @metadata slot.
+#' @export
 #' @examples
 #' # Generate random data
 #' x <- randomMultiFactor(n_features = 20)
@@ -22,33 +24,42 @@
 #' )
 #'
 #' # Now enrich test input
-#' weave(x, a ~ b ~ c)
-#'     #test_set_enrichment()
-#' @noRd
-test_set_enrichment <- function(
-        x, observed = NULL, set = NULL, alternative = "greater", raw = FALSE
-        ) {
+#' weave_coverage(x, a ~ b ~ c)
+#'
+weave_coverage <- function(
+        x, .path, .data = NULL, metric = c("count", "coverage", "complete"),
+        out.format = c("LinkMap", "matrix"), .data_column = "row.names"
+) {
+    metric <- match.arg(metric, c("count", "coverage", "complete"))
+    out.format <- match.arg(out.format, c("LinkMap", "matrix"))
+
+    .p_check <- .check_path(.path)
+
+    if(.p_check["vars"] == "complex") {
+        stop("weave_coverage() '.path' cannot contain '+'.\n",
+             "Use stack() to prepare input.")
+    }
+    full_path <- .path_ordinary_to_full(x, .path)[[1L]]
+
     stopifnot(
-        "x must be a MultiFactor" = inherits(x, "MultiFactor::MultiFactor")
-    )
-    if(is.null(observed)) {
-        stopifnot(
-            "If 'observed' is missing 'x' must have two LinkMaps" =
-                length(x) == 2L
+        "weave_coverage() '.path' must be 2 or 3 steps long." =
+            length(full_path) %in% c(2L, 3L)
         )
-        param_df <- .weave_contingency_params(x)
-    } else {
+    x <- subsetByPath(x, full_path)
 
-        param_df <- .char_contingency_params(x, observed, set)
+    .data <- .data_coverage_to_vector(.data, .data_column)
+    if( length(full_path) == 2L ) {
+        res <- .weave_coverage_two(x, full_path, .data, metric, out.format)
+    } else if( length(full_path) == 3L ) {
+        res <- .weave_coverage_three(x, full_path, .data, metric, out.format)
     }
 
-    cont_mats <- apply(param_df, 1L, .as_cont_matrix, simplify = FALSE)
-    res <- lapply(cont_mats, fisher.test, alternative = alternative)
-    if(!raw) res <- dplyr::bind_rows(lapply(res, tidy), .id = "term")
-
     return(res)
+
+
 }
 
+
 .weave_contingency_params <- function(x) {
     shared <- do.call(intersect, unname(lapply(x, names)))
     stopifnot( "LinkMaps must share exactly one column" = length(shared) == 1L )
@@ -110,45 +121,111 @@ test_set_enrichment <- 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.cpt <- .weave_coverage_three(x, all_terms = at, metric = "complete")
 #'
-.weave_summarize_links <- function(
-        x, all_terms, metric = c("count", "coverage", "complete")
+#' x.cov@metadata
+#' x.cnt@metadata
+#' x.cpt@metadata
+#'
+.weave_coverage_three <- function(
+        x, all_terms, .data, metric, out.format
         ) {
-    metric <- match.arg(metric, c("count", "coverage", "complete"))
-    seen <- x[[1L]]
-    full <- x[[2L]]
     # All steps in order
     from <- all_terms[[1L]]
     shared <- all_terms[[2L]]
     to <- all_terms[[3L]]
 
+    seen <- x[[rowsWithCol(x@map, c(shared, from))]]
+    full <- x[[rowsWithCol(x@map, c(shared, to))  ]]
+
     # Ensure LinkMap order
-    shared2from <- as.matrix(seen[, c(shared, from)])
-    shared2to <- as.matrix(full[, c(shared, to)])
+    shared2from <- as.matrix(seen, terms = c(shared, from))
+    shared2to <- as.matrix(full, terms = c(shared, to))
+
+    if( length(.data) ) { shared2from <- shared2from[, .data] }
     # Link observed features to sets
-    res <- Matrix::crossprod(
-        shared2to, shared2from != 0L
-    )
-    if(metric != "count") tot_set <- pmax(Matrix::colSums(shared2to), 1L)
-    if(metric == "coverage") res <- res/tot_set
-    if(metric == "complete") res <- res == tot_set
+    res <- Matrix::crossprod( shared2to, shared2from != 0L )
+    res <- .calc_coverage(res, bg = shared2to, metric)
+
+    if( out.format == "LinkMap" ) {
+        val <- data.frame(res@x)
+        colnames(val) <- metric
+        res <- .res_weave_matrix_to_LinkMap(res, levels(x)[c(from, to)])
+        res <- LinkMap(res, metadata = val)
+    }
+    return(res)
+}
+
+.weave_coverage_two <- function(
+        x, all_terms, .data, metric, out.format
+) {
+    bg <- as.matrix(x[[1L]], terms = all_terms)
+    if( length(.data) ) {
+        obs <- bg[.data, ]
+    } else {
+        obs <- bg
+    }
+    res <- .calc_coverage(obs, bg, metric)
+    if( out.format == "LinkMap" ) {
+        val <- data.frame(res@x)
+        colnames(val) <- metric
+        res <- .res_weave_matrix_to_LinkMap(res, levels(x)[all_terms])
+        res <- LinkMap(res, metadata = val)
+    }
+    return(res)
+}
+
+#' res (observed) and bg (background) are two sparse matrices.
+#' @noRd
+#' @importFrom Matrix colSums t Matrix
+#'
+.calc_coverage <- function(res, bg, metric) {
+    if( metric != "count" ) tot_set <- pmax(Matrix::colSums(bg), 1L)
+    if( metric == "coverage" ) res <- res/tot_set
+    if( metric == "complete" ) res <- res == tot_set
 
     res <- Matrix::t( Matrix::Matrix( res, sparse = TRUE ) )
     return(res)
 }
 
 
-# .weave_calc_coverage <- function(seen, full, to, from, shared) {
-#     # Ensure LinkMap order
-#     shared2from <- as.matrix(seen[, c(shared, from)])
-#     shared2to <- as.matrix(full[, c(shared, to)])
-#     # Link observed features to sets
-#     Matrix::Matrix(
-#         Matrix::crossprod(shared2to, shared2from!=0L) /
-#             Matrix::colSums(shared2to),
-#         sparse = TRUE
+# Legacy
+# path_coverage <- function(x, path, out.format = "matrix") {
+#     # rename to all_terms for internal consistency with .weave_*
+#     all_terms <- path
+#     # tolerate single path result in list
+#     if(length(all_terms) == 1L && is.list(all_terms)) all_terms <- all_terms[[1L]]
+#     stopifnot(
+#         "'path' must be a character vector of steps to take, in order." =
+#             is.character(all_terms)
 #     )
-# }
+#     stopifnot(
+#         "All entries in 'path' must be found in colnames(x)." =
+#             all( all_terms %in% colnames(x) )
+#     )
+#     stopifnot("length( path ) must be 3." = length( all_terms ) == 3L )
 #
-
-
+#     x <- subsetByPath(x, all_terms)
+#     terms <- all_terms[c(1L, length(all_terms))]
+#     # Compute coverage
+#     res <- .weave_coverage_three(x, all_terms, .data, "coverage")
+#
+#     # Check if we're done
+#     if(out.format == "matrix") {
+#         dimnames(res) <- levels(x)[terms]
+#         return(res)
+#     }
+#     # Otherwise, make a LinkMap
+#     res <- .res_weave_matrix_to_LinkMap(res, levels(x)[terms])
+#
+#     return(res)
+#
+# }
diff --git a/R/weave.R b/R/weave.R
index c57e046..4478d2f 100644
--- a/R/weave.R
+++ b/R/weave.R
@@ -67,37 +67,6 @@ weave_along_path <- function(x, path, out.format = "LinkMap") {
   res
 }
 
-path_coverage <- function(x, path, out.format = "matrix") {
-  # rename to all_terms for internal consistency with .weave_*
-  all_terms <- path
-  # tolerate single path result in list
-  if(length(all_terms) == 1L && is.list(all_terms)) all_terms <- all_terms[[1L]]
-  stopifnot(
-    "'path' must be a character vector of steps to take, in order." =
-      is.character(all_terms)
-  )
-  stopifnot(
-    "All entries in 'path' must be found in colnames(x)." =
-      all( all_terms %in% colnames(x) )
-  )
-  stopifnot("length( path ) must be 3." = length( all_terms ) == 3L )
-
-  x <- subsetByPath(x, all_terms)
-  terms <- all_terms[c(1L, length(all_terms))]
-  # Compute coverage
-  res <- .weave_summarize_links(x, all_terms, "coverage")
-
-  # Check if we're done
-  if(out.format == "matrix") {
-    dimnames(res) <- levels(x)[terms]
-    return(res)
-  }
-  # Otherwise, make a LinkMap
-  res <- .res_weave_matrix_to_LinkMap(res, levels(x)[terms])
-
-  return(res)
-
-}
 
 .weave_complex_formula <- function(x, .path, out.format) {
   .path_list <- .path_prep_complex(.path)
@@ -160,7 +129,8 @@ path_coverage <- function(x, path, out.format = "matrix") {
 #' @importFrom Matrix which
 #'
 .res_weave_matrix_to_LinkMap <- function(res, lvs) {
-  res <- as.data.frame.matrix(Matrix::which(res, arr.ind = TRUE))
+
+  res <- as.data.frame.matrix(Matrix::which(res > 0L, arr.ind = TRUE))
   res[] <- mapply(FUN = function(x, y) {
     attr(x, "levels") <- y
     `class<-`(x, "factor")
diff --git a/man/weave_coverage.Rd b/man/weave_coverage.Rd
new file mode 100644
index 0000000..b036a22
--- /dev/null
+++ b/man/weave_coverage.Rd
@@ -0,0 +1,54 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/weave-coverage.R
+\name{weave_coverage}
+\alias{weave_coverage}
+\title{Perform enrichtment analysis from a weave}
+\usage{
+weave_coverage(
+  x,
+  .path,
+  .data = NULL,
+  metric = c("count", "coverage", "complete"),
+  out.format = c("LinkMap", "matrix"),
+  .data_column = "row.names"
+)
+}
+\arguments{
+\item{x}{a \code{MultiFactor}}
+
+\item{.path}{Either a \code{formula} or a \verb{character vector} of length 2 with the
+names of the desired combination of feature types.}
+
+\item{.data}{Optional \verb{Character vector}. Lists observed features from the
+found within the first element of \code{.path}. Alternatively, a \code{data.frame}
+with the same information. (Also see \code{.data_column} argument).}
+
+\item{metric}{\verb{Character scalar}.
+One of \code{'count'}, \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}
+}
+\value{
+a \code{LinkMap} or \code{matrix} with the desired coverage information in
+the @metadata slot.
+}
+\description{
+Perform enrichtment analysis from a weave
+}
+\examples{
+# Generate random data
+x <- randomMultiFactor(n_features = 20)
+
+# Spike in a lower number of a observations
+x <- MultiFactor(
+    list(x[[1]][sample(size = 5, 1:NROW(x[[1]])),], x[[2]])
+)
+
+# Now enrich test input
+weave_coverage(x, a ~ b ~ c)
+
+}
diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml
index f8442f2..f31e270 100644
--- a/pkgdown/_pkgdown.yml
+++ b/pkgdown/_pkgdown.yml
@@ -19,6 +19,7 @@ reference:
 - contents:
   - weave
   - weave-methods
+  - weave_coverage
   - stack
   - select_path
   - subgroup_apply
diff --git a/vignettes/MultiFactor.Rmd b/vignettes/MultiFactor.Rmd
index 1c04443..2f9c378 100644
--- a/vignettes/MultiFactor.Rmd
+++ b/vignettes/MultiFactor.Rmd
@@ -68,7 +68,7 @@ e <- tapply(trade_posts(raw.data = TRUE), ~group, `[[<-`, "group", value = NULL)
     e[[l]][match(x[[l]], levels(x)[[l]]), ],
     e[[r]][match(x[[r]], levels(x)[[r]]), ]
     ), 
-  y = rep(seq_len(NROW(x)), 2), 
+  y = rep(rev(seq_len(NROW(x))), 2), 
   x = rep(c(-1, -1/2), each = NROW(x))
   ), 
   data.frame(
@@ -153,8 +153,8 @@ ggplot() +
     rev(plot_dfs[-1]), function(x) list(
       geom_tile(data = x, aes(fill = header), colour = "black", linewidth = 1/2), 
       geom_text(data = filter(x, header),  aes(label = name), size = 4),
-      geom_text(data = filter(x, !header), aes(label = emoji), size = 5)
-      )
+      geom_text(data = filter(x, !header), aes(label = emoji), size = 5, family = "Noto Color Emoji")
+      ) 
 ) +
   geom_tile( data = first_linkmap, aes(fill = header), colour = "black", linewidth = 1/2 ) +
   geom_text( data = filter(first_linkmap, header), aes(label = name), size = 4 ) +