Feedback/santi changes

.Rbuildignore

@@ -0,0 +1,11 @@
+^\.Rproj$
+^\.Rproj\.user$
+^Makefile$
+^CLAUDE\.md$
+^AGENTS\.md$
+^\.claude$
+^\.semquery$
+^\.beads$
+^local_data$
+^\.git$
+^\.gitignore$

.gitignore

@@ -5,3 +5,4 @@
 *.o
 *.so
 *.dll
+tests/*.pdf

DESCRIPTION

@@ -22,6 +22,7 @@ RoxygenNote: 7.3.3
 Depends: R (>= 4.3.0)
 Imports:
     igraph,
+    magrittr,
     plotly,
     visNetwork,
     matrixStats,
@@ -33,7 +34,9 @@ Suggests:
     testthat (>= 3.0.0),
     knitr,
     rmarkdown,
-    BiocStyle
+    BiocStyle,
+    Rtsne,
+    uwot
 VignetteBuilder: knitr
 Config/testthat/edition: 3
 Config/roxygen2/version: 8.0.0

Makefile

@@ -1,4 +1,4 @@
-PKGNAME := WGCNAplus
+PKGNAME := lasagna
 
 .PHONY: doc build install test check clean
 

NAMESPACE

@@ -14,6 +14,7 @@ importFrom(grDevices,adjustcolor)
 importFrom(grDevices,rainbow)
 importFrom(graphics,plot)
 importFrom(graphics,text)
+importFrom(magrittr,"%>%")
 importFrom(matrixStats,rowSds)
 importFrom(plotly,add_markers)
 importFrom(plotly,add_mesh)

R/graph.R

@@ -43,61 +43,127 @@ prune_graph <- function(graph,
                         select = NULL,
                         prune = TRUE) {
 
-  if (is.null(layers)) layers <- graph$layers
-  if (is.null(layers)) layers <- unique(igraph::V(graph)$layer)
-  layers <- setdiff(layers, c("SOURCE", "SINK"))
+  if (is.null(graph)) {
+    message("[lasagna::prune_graph]: graph is NULL. Please check run of create_model() and solve(). Exiting.")
+    return(NULL)
+  }
+
+  layers <- resolve_layers(graph, layers)
   graph <- igraph::subgraph(graph, igraph::V(graph)$layer %in% layers)
 
   if (!"value" %in% names(igraph::vertex_attr(graph))) {
     stop("vertex must have 'value' attribute")
   }
 
-  ## select nodes/modules
-  if (!is.null(select)) {
-    v1 <- (igraph::V(graph)$name %in% select)
-    v2 <- (sub(".*:", "", igraph::V(graph)$name) %in% select)
-    v3 <- (sub(":.*", "", igraph::V(graph)$name) %in% select)
-    graph <- igraph::subgraph(graph, vids = which(v1 | v2 | v3))
-  }
-
   if (!is.null(filter)) {
     if (!is.list(filter)) stop("filter must be a named list")
     if (is.null(names(filter))) stop("filter must be a named list")
-    for (k in names(filter)) {
-      vv <- igraph::V(graph)$name
-      filt <- filter[[k]]
-      vids <- igraph::V(graph)$layer != k | grepl(filt, vv, ignore.case = TRUE)
-      graph <- igraph::subgraph(graph, which(vids))
-    }
   }
 
-  ## select ntop features
-  fc <- igraph::V(graph)$value
-  names(fc) <- igraph::V(graph)$name
-  if (!is.null(ntop) && ntop > 0) {
-    ii <- tapply(
-      seq_along(fc), igraph::V(graph)$layer,
-      function(i) utils::head(i[order(-abs(fc[i]))], ntop)
-    )
-    ii <- unlist(ii[names(ii) %in% layers])
-    fc <- fc[ii]
-    graph <- igraph::subgraph(graph, igraph::V(graph)[ii])
-  }
+  ## reduce the vertex set
+  if (!is.null(select)) graph <- keep_named_nodes(graph, select)
+  if (!is.null(filter)) graph <- keep_matching_nodes(graph, filter)
+  if (!is.null(ntop) && ntop > 0) graph <- keep_top_nodes(graph, ntop, layers)
+
+  ## zero the weight of every unwanted edge, then drop them in one go
+  if (normalize.edges) graph <- normalize_edge_weights(graph)
+  if (min.rho > 0) graph <- zero_weak_edges(graph, min.rho)
+  graph <- zero_edges_by_sign(graph, edge.sign, layers)
+  graph <- zero_edges_by_type(graph, edge.type)
+
+  return(drop_zero_edges(graph, prune))
+
+}
+
+## Layers the graph is to be reduced to: the caller's choice, otherwise the
+## layers attached to the graph object, otherwise those found on its
+## vertices. SOURCE and SINK are helper nodes and are never kept.
+resolve_layers <- function(graph, layers) {
+
+  if (is.null(layers)) layers <- graph$layers
+  if (is.null(layers)) layers <- unique(igraph::V(graph)$layer)
 
-  if (normalize.edges) {
-    for (e in unique(igraph::E(graph)$connection_type)) {
-      ii <- which(igraph::E(graph)$connection_type == e)
-      max.wt <- max(abs(igraph::E(graph)$weight[ii]), na.rm = TRUE) + 1e-3
-      igraph::E(graph)$weight[ii] <- igraph::E(graph)$weight[ii] / max.wt
-    }
+  return(setdiff(layers, c("SOURCE", "SINK")))
+
+}
+
+## Keep the nodes named in 'select', matching either the full node name or
+## the part before or after its layer prefix, so that a module can be
+## selected by bare name as well as by qualified name.
+keep_named_nodes <- function(graph, select) {
+
+  vv <- igraph::V(graph)$name
+  v1 <- (vv %in% select)
+  v2 <- (sub(".*:", "", vv) %in% select)
+  v3 <- (sub(":.*", "", vv) %in% select)
+
+  return(igraph::subgraph(graph, vids = which(v1 | v2 | v3)))
+
+}
+
+## Keep the nodes matching the per-layer regular expressions in 'filter', a
+## list of patterns named by layer. Each pattern only constrains its own
+## layer; nodes of the layers not named in the list are always kept.
+keep_matching_nodes <- function(graph, filter) {
+
+  for (k in names(filter)) {
+    vv <- igraph::V(graph)$name
+    filt <- filter[[k]]
+    vids <- igraph::V(graph)$layer != k | grepl(filt, vv, ignore.case = TRUE)
+    graph <- igraph::subgraph(graph, which(vids))
   }
 
-  if (min.rho > 0) {
-    ii <- which(abs(igraph::E(graph)$weight) < min.rho)
-    if (length(ii)) igraph::E(graph)$weight[ii] <- 0
+  return(graph)
+
+}
+
+## Keep the 'ntop' nodes of largest absolute value within each layer. Layers
+## outside 'layers' are dropped entirely, which matters for the SOURCE/SINK
+## nodes that survive an explicit layer selection.
+keep_top_nodes <- function(graph, ntop, layers) {
+
+  fc <- igraph::V(graph)$value
+  ii <- tapply(
+    seq_along(fc), igraph::V(graph)$layer,
+    function(i) utils::head(i[order(-abs(fc[i]))], ntop)
+  )
+  ii <- unlist(ii[names(ii) %in% layers])
+
+  return(igraph::subgraph(graph, igraph::V(graph)[ii]))
+
+}
+
+## Rescale the edge weights of each connection type onto a common scale, so
+## that weights of different connection types can be compared.
+normalize_edge_weights <- function(graph) {
+
+  for (e in unique(igraph::E(graph)$connection_type)) {
+    ii <- which(igraph::E(graph)$connection_type == e)
+    max.wt <- max(abs(igraph::E(graph)$weight[ii]), na.rm = TRUE) + 1e-3
+    igraph::E(graph)$weight[ii] <- igraph::E(graph)$weight[ii] / max.wt
   }
 
+  return(graph)
+
+}
+
+## Zero the edges weaker than 'min.rho'.
+zero_weak_edges <- function(graph, min.rho) {
+
+  ii <- which(abs(igraph::E(graph)$weight) < min.rho)
+  if (length(ii)) igraph::E(graph)$weight[ii] <- 0
+
+  return(graph)
+
+}
+
+## Zero the edges of unwanted sign. "consensus" keeps the edges whose sign
+## agrees with the expected direction of their source layer, which is
+## inverted for the miRNA layers as those anticorrelate with their targets.
+zero_edges_by_sign <- function(graph, edge.sign, layers) {
+
   ewt <- igraph::E(graph)$weight
+
   if (grepl("pos", edge.sign)) {
     igraph::E(graph)$weight[ewt < 0] <- 0
   } else if (grepl("neg", edge.sign)) {
@@ -111,8 +177,17 @@ prune_graph <- function(graph,
     igraph::E(graph)$weight <- ewt * (sign(ewt) == esign)
   }
 
-  ## delete intra or inter edges
+  return(graph)
+
+}
+
+## Zero the edges of unwanted type, inter-layer edges being those whose
+## connection type names two layers. "both2" keeps both types but drops the
+## negative intra-layer edges.
+zero_edges_by_type <- function(graph, edge.type) {
+
   ic <- grepl("->", igraph::E(graph)$connection_type)
+
   if (edge.type == "inter") {
     igraph::E(graph)$weight[!ic] <- 0
   } else if (edge.type == "intra") {
@@ -121,8 +196,16 @@ prune_graph <- function(graph,
     sel <- (!ic & igraph::E(graph)$weight < 0)
     igraph::E(graph)$weight[sel] <- 0
   }
-  graph <- igraph::delete_edges(graph, which(igraph::E(graph)$weight == 0))
 
+  return(graph)
+
+}
+
+## Drop the edges that the filters zeroed out and, with prune=TRUE, the
+## vertices that are left without any edge.
+drop_zero_edges <- function(graph, prune) {
+
+  graph <- igraph::delete_edges(graph, which(igraph::E(graph)$weight == 0))
   if (prune) graph <- igraph::subgraph_from_edges(graph, igraph::E(graph))
 
   return(graph)

R/lasagna-package.R

@@ -2,6 +2,7 @@
 "_PACKAGE"
 
 #' @import igraph
+#' @importFrom magrittr %>%
 #' @importFrom matrixStats rowSds
 #' @importFrom stats cor median model.matrix sd
 #' @importFrom utils head tail type.convert