Merge pull request #88 from natverse/fix/neuprint_ROI_connectivity

Fix/neuprint roi connectivity
natverse · Mar 2, 2020 · bcf5977 · bcf5977
2 parents 7b42f48 + 34c0a25
commit bcf5977
Show file tree

Hide file tree

Showing 5 changed files with 130 additions and 46 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -37,7 +37,8 @@ Imports:
     memoise,
     dplyr,
     bit64,
-    stringr
+    stringr,
+    tibble
 Remotes: 
     natverse/drvid,
     natverse/nat

diff --git a/R/connectivity.R b/R/connectivity.R
@@ -476,24 +476,20 @@ extract_connectivity_df <- function(rois, json){
   if(is.null(json)){
     return(NULL)
   }
-  a <- unlist(jsonlite::fromJSON(json))
-  values <-  data.frame(row.names = 1)
   rois <- unique(rois) #this takes care if both the input and output ROIs are same..
+  a <- unlist(jsonlite::fromJSON(json))
+  roicols <- c(t(outer(rois, c("pre", "post"), paste, sep=".")))
+  values <- tibble::as_tibble(as.list(structure(rep(0, length(roicols)), .Names=roicols)))
   for(roi in rois){
-    d <-  data.frame(0,0)
-    colnames(d) <- paste0(roi,c(".pre",".post"))
+    thisroicols <- paste0(roi,c(".pre",".post"))
     if (!is.null(a)){
       b <-  a[startsWith(names(a),paste0(roi,"."))]
-      d[names(b)] <-  b
+      values[names(b)] <-  b
     }
-    values <- cbind(values,d)
   }
   values
 }
 
-
-
-
 ##' @title Get a matrix for connectivity between neuron/neuronlist objects
 #'
 #' @description Get an adjacency matrix for the synaptic connectivity between \code{nat::neuron}/\code{nat::neuronlist} objects. This function does not query a neuPrint server.

diff --git a/R/roi.R b/R/roi.R
@@ -72,55 +72,99 @@ neuprint_bodies_in_ROI <- function(roi, dataset = NULL, all_segments = FALSE, co
   df
 }
 
-#' @title Get the connectivity between ROIs in a neuPrint dataset
+#' Get connectivity between ROIs (summary or data frame of connecting neurons)
 #'
+#' @details When requesting summary connectivity data between ROIs, we recommend
+#'   setting \code{cached=FALSE}. We have noticed small differences in the
+#'   connections weights, but computation times can get very long for more than
+#'   a handful of ROIs.
 #' @param rois regions of interest for a dataset
-#' @param cached pull cached results (TRUE) or recalculate the connectivity (FALSE)?
-#' @param full return all neurons involved (TRUE, the default) or give a ROI summary (FALSE, default behavior if `cached` is TRUE)
-#' @param statistic either "weight" or count" (default "weight"). Which number to return (see neuprint explorer for details) for summary results (either `full` is FALSE or `cached` is TRUE)
+#' @param full return all neurons involved (TRUE, the default) or give a numeric
+#'   ROI summary (FALSE)
+#' @param statistic either "weight" or count" (default "weight"). Which number
+#'   to return (see neuprint explorer for details) for summary results (either
+#'   (when \code{full=FALSE})
+#' @param cached pull precomputed results (TRUE) or ask server to recalculate
+#'   the connectivity (FALSE). Only applicable to summary results when
+#'   \code{full=FALSE}.
 #' @param ... methods passed to \code{neuprint_login}
 #' @inheritParams neuprint_fetch_custom
-#' @seealso \code{\link{neuprint_simple_connectivity}}, \code{\link{neuprint_common_connectivity}}
+#' @seealso \code{\link{neuprint_simple_connectivity}},
+#'   \code{\link{neuprint_common_connectivity}}
 #' @export
-#' @rdname neuprint_ROI_connectivity
-neuprint_ROI_connectivity <- function(rois, cached = FALSE, full=TRUE, statistic = c("weight","count"),dataset = NULL, conn = NULL, ...){
+#' @examples
+#' \donttest{
+#' aba <- neuprint_ROI_connectivity(neuprint_ROIs(superLevel = TRUE),
+#'   full=FALSE)
+#' heatmap(aba)
+#' }
+neuprint_ROI_connectivity <- function(rois, full=TRUE,
+                                      statistic = c("weight","count"),
+                                      cached = !full,
+                                      dataset = NULL, conn = NULL, ...) {
   statistic <- match.arg(statistic)
+  if(isTRUE(full) && isTRUE(cached))
+    stop("It is not possible to return a full list of connecting neurons when ",
+         "`cached=TRUE`!\nPlease leave `cached` with its default value (FALSE).")
   roicheck <- neuprint_check_roi(rois=rois, dataset = dataset, conn = conn, ...)
-  if (cached){
-    results <-matrix(nrow=length(rois),ncol=length(rois),dimnames = list(inputs=rois,outputs=rois))
+  if (cached) {
+    results <-matrix(ifelse(statistic == 'count', 0L, 0),
+                     nrow=length(rois), ncol=length(rois),
+                     dimnames = list(inputs=rois,outputs=rois))
     roi.conn = neuprint_fetch(path = 'api/cached/roiconnectivity', conn = conn, ...)
-    for (inp in rois){
-      for (out in rois){
-         results[inp,out] <-  roi.conn$weights[[paste(inp,out,sep="=>")]][[statistic]]
+    missing=setdiff(rois, unlist(roi.conn$roi_names))
+    if(length(missing))
+      warning("Dropping missing rois:", paste(missing, collapse = " "))
+    allpairs = names(roi.conn$weights)
+    for (inp in rois) {
+      for (out in rois) {
+        edgename=paste(inp, out, sep="=>")
+        if(edgename %in% allpairs)
+          results[inp,out] <- roi.conn$weights[[edgename]][[statistic]]
       }
     }
-  }else{
+  } else {
     Payload = noquote(sprintf('{"dataset":"%s","rois":%s}',
                             dataset,
                             ifelse(is.null(rois),jsonlite::toJSON(list()),jsonlite::toJSON(rois))))
     class(Payload) = "json"
     roi.conn <- neuprint_fetch(path = 'api/npexplorer/roiconnectivity', body = Payload, conn = conn, ...)
-    connData <- roi.conn$data[sapply(roi.conn$data,function(d) any(sapply(rois,function(r) grepl(paste0("\"",r,"\""), d[[2]]))))]
-    connections <-lapply(connData, function(rc) extract_connectivity_df(rois=rois,json=rc[[2]]))
-    resultsD <- dplyr::bind_rows(connections)
-    resultsD$bodyid <- as.character(sapply(connData, function(d) d[[1]]))
-    if (!full){
-      results <- matrix(nrow=length(rois),ncol=length(rois),dimnames = list(inputs=rois,outputs=rois))
-      if (statistic == "count"){
-        for (inp in rois){
-          for (out in rois){
-            results[inp,out] <- length(which(resultsD[[paste0(inp,".post")]]>0 & resultsD[[paste0(out,".pre")]]>0))
+    ll <- neuprint_list2df(roi.conn)
+    # running fromJSON on many separate strings is slow, so start by
+    # selecting strings that actually contain the selected ROIs
+    hasroi=sapply(rois, function(roi)
+      stringr::str_detect(ll$roiInfo, stringr::fixed(paste0('"',roi,'"'))))
+    if(is.matrix(hasroi)) hasroi=rowSums(hasroi)>0
+
+    connections <-lapply(ll$roiInfo[hasroi],
+                         function(x) extract_connectivity_df(rois=rois,json=x))
+    resultsD <- cbind(ll[hasroi, 1, drop=FALSE], dplyr::bind_rows(connections))
+    if (!full) {
+      results <-
+        matrix(
+          nrow = length(rois),
+          ncol = length(rois),
+          dimnames = list(inputs = rois, outputs = rois)
+        )
+      if (statistic == "count") {
+        for (inp in rois) {
+          for (out in rois) {
+            results[inp, out] <-
+              length(which(resultsD[[paste0(inp, ".post")]] > 0 &
+                             resultsD[[paste0(out, ".pre")]] > 0))
           }
         }
-      }else{
+      } else{
         totalInputs <- neuprint_get_meta(resultsD$bodyid)$post
-        for (inp in rois){
-          for (out in rois){
-            results[inp,out] <- sum((resultsD[[paste0(out,".pre")]]*resultsD[[paste0(inp,".post")]]/totalInputs)[totalInputs>0])
+        for (inp in rois) {
+          for (out in rois) {
+            results[inp, out] <-
+              sum((resultsD[[paste0(out, ".pre")]] * resultsD[[paste0(inp, ".post")]] /
+                     totalInputs)[totalInputs > 0])
           }
         }
       }
-    }else{
+    } else {
       results <- resultsD
     }
   }

diff --git a/man/neuprint_ROI_connectivity.Rd b/man/neuprint_ROI_connectivity.Rd
diff --git a/tests/testthat/test-roi.R b/tests/testthat/test-roi.R
@@ -1,3 +1,13 @@
+test_that("", {
+  json='{"SNP(R)": {"pre": 71, "post": 155}, "SLP(R)": {"pre": 67, "post": 153}, "SIP(R)": {"pre": 4, "post": 2}, "LH(R)": {"pre": 20, "post": 25}, "VLNP(R)": {"pre": 1}, "PLP(R)": {"pre": 1}, "AL(R)": {"pre": 1, "post": 162}}'
+  rois=c("AL(R)", "LH(R)")
+
+  baseline <- structure(list(`AL(R).pre` = 1L, `AL(R).post` = 162L,
+                             `LH(R).pre` = 20L, `LH(R).post` = 25L),
+                        class = "data.frame", row.names = "1")
+  expect_equal(xdf <- extract_connectivity_df(rois = rois, json=json), baseline)
+})
+
 skip_if(as.logical(Sys.getenv("SKIP_NP_SERVER_TESTS")))
 
 test_that("neuprint_bodies_in_ROI works", {
@@ -9,8 +19,22 @@ test_that("neuprint_bodies_in_ROI works", {
     'data.frame')
 })
 
-
 test_that("neuprint_ROI_mesh works", {
   expect_is(m <- neuprint_ROI_mesh('ATL(L)'),
     'mesh3d')
 })
+
+test_that("neuprint_ROI_connectivity works", {
+  rois <- neuprint_ROIs(superLevel = TRUE)[1:2]
+  expect_is(m <- neuprint_ROI_connectivity(rois, full=F), 'matrix')
+  expect_equal(dimnames(m),
+               list(inputs = rois, outputs = rois))
+  expect_error(neuprint_ROI_connectivity(rois, full = T, cached = T))
+
+  expect_is(m2 <- neuprint_ROI_connectivity(rois, full=F, cached=F), 'matrix')
+  # note low tolerance as the cached and recomputed results are not identical
+  expect_equal(m2, m, tolerance = 1e-2)
+
+  expect_is(df <- neuprint_ROI_connectivity(rois[1], full=T), 'data.frame')
+  expect_true(ncol(df)==3)
+})