Topological NBLAST

DESCRIPTION

@@ -7,7 +7,9 @@ Authors@R: c(
     person("Gregory", "Jefferis", email= "[email protected]",
       role = c("aut"), comment = c(ORCID = "0000-0002-0587-9355")),
     person("James", "Manton", email="[email protected]",
-      role = c("aut", "cre"), comment = c(ORCID = "0000-0001-9260-3156"))
+      role = c("aut", "cre"), comment = c(ORCID = "0000-0001-9260-3156")),
+    person("Dominik", "Krzeminski", role = c("ctb"),
+      comment = c(ORCID = "0000-0003-4568-0583"))
     )
 Description: Extends package 'nat' (NeuroAnatomy Toolbox) by providing a
     collection of NBLAST-related functions for neuronal morphology comparison (Costa et al. (2016) <doi: 10.1016/j.neuron.2016.06.012>).
@@ -33,6 +35,6 @@ Suggests:
 License: GPL-3
 LazyData: yes
 VignetteBuilder: knitr
-RoxygenNote: 7.1.1
+RoxygenNote: 7.2.1
 Language: en-GB
 Encoding: UTF-8

R/neuriteblast.r

@@ -46,6 +46,8 @@
 #' @param UseAlpha whether to weight the similarity score for each matched
 #'   segment to emphasise long range neurites rather then arbours (default:
 #'   FALSE, see \bold{\code{UseAlpha}} section for details).
+#' @param UseTopo whether to use directed dotprops vectors (pointing towards
+#'    soma) and features in neuron nodes (default: FALSE).
 #' @param OmitFailures Whether to omit neurons for which \code{FUN} gives an
 #'   error. The default value (\code{NA}) will result in \code{nblast} stopping
 #'   with an error message the moment there is an error. For other values, see
@@ -164,7 +166,7 @@
 #' }
 nblast <- function(query, target=getOption("nat.default.neuronlist"),
                    smat=NULL, sd=3, version=c(2, 1), normalised=FALSE,
-                   UseAlpha=FALSE, OmitFailures=NA, ...) {
+                   UseAlpha=FALSE, UseTopo = FALSE, OmitFailures=NA, ...) {
   version <- as.character(version)
   version <- match.arg(version, c('2', '1'))
 
@@ -187,10 +189,11 @@ nblast <- function(query, target=getOption("nat.default.neuronlist"),
     }
     if(is.character(smat)) smat=get(smat)
     NeuriteBlast(query=query, target=target, NNDistFun=lodsby2dhist, smat=smat,
-                 UseAlpha=UseAlpha, normalised=normalised, OmitFailures=OmitFailures, ...)
+                 UseAlpha=UseAlpha, UseTopo=UseTopo, normalised=normalised,
+                 OmitFailures=OmitFailures, ...)
   } else if(version == '1') {
     NeuriteBlast(query=query, target=target, NNDistFun=WeightedNNBasedLinesetDistFun,
-                 UseAlpha=UseAlpha, sd=sd, normalised=normalised,
+                 UseAlpha=UseAlpha, UseTopo=UseTopo, sd=sd, normalised=normalised,
                  OmitFailures=OmitFailures, ...)
   } else {
     stop("Only NBLAST versions 1 and 2 are currently implemented. For more advanced control, see NeuriteBlast.")
@@ -370,6 +373,8 @@ WeightedNNBasedLinesetMatching <- function(target, query, ...) {
 #'   WeightedNNBasedLinesetMatching. These will be used to scale the dot
 #'   products of the direction vectors for nearest neighbour pairs.
 #' @param UseAlpha Whether to scale dot product of tangent vectors (default=F)
+#' @param UseTopo Whether to use topological information to scale dot products
+#'  of tangent vectors (default=F)
 #' @param ... extra arguments to pass to the distance function.
 #' @export
 #' @seealso \code{\link[nat]{dotprops}}
@@ -380,12 +385,18 @@ WeightedNNBasedLinesetMatching <- function(target, query, ...) {
 #' segvals=WeightedNNBasedLinesetMatching(kcs20[[1]], kcs20[[2]], NNDistFun=list)
 #' names(segvals)=c("dist", "adotprod")
 #' pairs(segvals)
-WeightedNNBasedLinesetMatching.dotprops<-function(target, query, UseAlpha=FALSE, ...) {
+WeightedNNBasedLinesetMatching.dotprops<-function(target, query, UseAlpha=FALSE,
+                                                  UseTopo=FALSE,
+                                                  ...) {
   if(!"dotprops" %in% class(query)) query <- dotprops(query)
   if(UseAlpha)
     WeightedNNBasedLinesetMatching(target$points,query$points,dvs1=target$vect,dvs2=query$vect,
                                    alphas1=target$alpha,alphas2=query$alpha,...)
-  else
+  else if(UseTopo) {
+    if(!"topo.dotprops" %in% class(query)) stop("query must be `topo.dotprops`")
+    WeightedNNBasedLinesetMatching(target$points,query$points,dvs1=target$vect,dvs2=query$vect,
+                                   alphas1=target$topo,alphas2=query$topo,...)
+  } else
     WeightedNNBasedLinesetMatching(target$points,query$points,dvs1=target$vect,dvs2=query$vect,...)
 }
 
@@ -396,10 +407,13 @@ WeightedNNBasedLinesetMatching.dotprops<-function(target, query, UseAlpha=FALSE,
 #' @rdname WeightedNNBasedLinesetMatching
 #' @importFrom nat dotprops
 WeightedNNBasedLinesetMatching.neuron<-function(target, query, UseAlpha=FALSE,
+                                                UseTopo=FALSE,
                                                 OnlyClosestPoints=FALSE,...) {
   if(!"neuron" %in% class(query)) {
     target <- dotprops(target)
-    return(WeightedNNBasedLinesetMatching(target=target, query=query, UseAlpha=UseAlpha, OnlyClosestPoints=OnlyClosestPoints, ...))
+    return(WeightedNNBasedLinesetMatching(target=target, query=query,
+                                          UseAlpha=UseAlpha, UseTopo=FALSE,
+                                          OnlyClosestPoints=OnlyClosestPoints, ...))
   }
   if(UseAlpha)
     stop("UseAlpha is not yet implemented for neurons!")
@@ -467,21 +481,44 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL
       idxArray=idxArray[!targetdupes,,drop=FALSE]
       nntarget$nn.dists=nntarget$nn.dists[!targetdupes]
     }
-    dps=abs(dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],]))
+    if (!is.list(alphas1)) {
+      dps = abs(dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],]))
+    } else {
+      dps = dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],])
+      dps[dps < 0] = 0
+    }
     if(!is.null(alphas1) && !is.null(alphas2)){
       # for perfectly aligned points, alpha = 1, at worst alpha = 0
       # sqrt seems reasonable since if alpha1=alpha2=0.5 then
       # the scalefac will be 0.5
       # zapsmall to make sure there are no tiny negative numbers etc.
-      scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]]))
+      #scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]]))
+      if (is.numeric(alphas1) && is.numeric(alphas2)) {
+        scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]]))
+      } else if (is.list(alphas1) && is.list(alphas2)) {
+        # use local neuron topology
+        if (any(names(alphas1) != names(alphas2)))
+          stop("Topo features not consistent between neurons.")
+        compare_feature <- function(feat) {
+          #' compares the features between points according to metric
+          #' 1 - ( |f_1 - f_2| / max(features))^2
+          maxlen = max(max(alphas1[[feat]]), max(alphas2[[feat]])) + .Machine$double.eps
+          zapsmall(1 - (abs(alphas1[[feat]][idxArray[,1]] - alphas2[[feat]][idxArray[,2]])/maxlen)^2)
+        }
+        sim_scores <- sapply(names(alphas1), compare_feature,
+                             USE.NAMES = TRUE, simplify = TRUE)
+        # Aggregate feature scores
+        scalefac = apply(sim_scores, 1, mean)
+      } else{
+        warning("Unknown alpha format")
+        scalefac = rep(1, length(dps))
+      }
       dps=dps*scalefac
     }
   }
-
   NNDistFun(as.vector(nntarget$nn.dists),dps,...)
 }
 
-
 dotprod=function(a,b){
   # expects 2 matrices with n cols each
   c=a*b

tests/testthat/test-NBLAST1.r

@@ -25,3 +25,10 @@ test_that("we can calculate normalised nblast v1 scores", {
   scores.norm=scale(scores, scale=c(scores[1,1],scores[2,2]), center = FALSE)
   expect_equivalent(nblast(testneurons[1:2], testneurons, version=1, normalised=TRUE), scores.norm)
 })
+
+testtopodps <- readRDS('testdata/testtopodps.rds')
+
+test_that("topoNBLAST works correct", {
+  scores <- nblast_allbyall(testtopodps, version=1, UseTopo = TRUE)
+  expect_true(scores["18820","20262"] < 0.1)
+})

tests/testthat/test-NBLAST2.r

@@ -100,3 +100,14 @@ test_that("we can handle all combinations of dotprops and neurons, both as neuro
   expect_is(nblast(Cell07PNs[[1]], testneurons[[1]]), 'numeric')
   expect_is(nblast(Cell07PNs[1:3], testneurons[[1]]), 'numeric')
 })
+
+
+testtopodps <- readRDS('testdata/testtopodps.rds')
+
+test_that("topoNBLAST works correct", {
+  scores <- nblast_allbyall(testtopodps, version=2, UseTopo = TRUE)
+  expect_true(scores["18820","20262"] > 0)
+  scores <- nblast_allbyall(testtopodps, version=2,
+                            normalisation = "normalised", UseTopo = TRUE)
+  expect_true(scores["18820","20262"] > 0.9)
+})

utils.R

@@ -0,0 +1,27 @@
+# utils
+
+get_dist_to_soma <- function(nrn) {
+  gw <- as.ngraph(nrn, weights=TRUE)
+  dst <- igraph::distances(gw, v = rootpoints(nrn))
+  as.numeric(dst)
+}
+
+make_topo_dotprops <- function(nrn, resample = 1, k = 5) {
+  tdps <- nat::dotprops(nrn, resample = resample, k = k, .parallel=TRUE)
+  tdps$alpha <- get_dist_to_soma(nrn)
+  tdps
+}
+
+make_sotopo_dotprops <- function(nrn, resample = 1, k = 5) {
+  tdps <- nat::dotprops(nrn, resample = resample, k = k, .parallel=TRUE)
+  tdps$alpha <- list()
+  tdps$alpha$distance <- get_dist_to_soma(nrn)
+  so <- strahler_order(nrn)
+  tdps$alpha$so <- abs(so$points-max(so$points)) # normalizing so the main branch is always 0
+  tdps
+}
+
+encode_ordinal <- function(x, order = unique(x)) {
+  x <- as.numeric(factor(x, levels = order, exclude = NULL))
+  x
+}

workbook.R

@@ -0,0 +1,72 @@
+library(nat)
+library(catmaid)
+library(dendroextras)
+
+devtools::load_all()
+
+catmaid_login()
+
+neurons <- read_catmaid_selection("../tnblast/data/DA1s.json", readNeurons = TRUE)
+neurons <- resample(prune_twigs(neurons/1e3, twig_length = 2), stepsize=1)
+
+plot3d(neurons, soma = TRUE)
+
+#somaid(neurons) # resampled neurons don't have soma
+#igraph::topo_sort(as.ngraph(nrn4))
+#max(igraph::diameter(gw1), igraph::diameter(gw2))
+
+dps_list <- nlapply(neurons, make_topo_dotprops)
+
+dps_aba1 <- nblast_allbyall(dps_list, normalisation = "raw")
+dps_aba2 <- nblast_allbyall(dps_list, UseAlpha = T, normalisation = "raw")
+
+par(mfrow = c(1,2))
+image(dps_aba1 / diag(dps_aba1), zlim= c(0,1), yaxt='n', xaxt='n')
+title("NBLAST")
+image(dps_aba2 / diag(dps_aba2), zlim= c(0,1), yaxt='n', xaxt='n')
+title("TNBLAST")
+
+dps_list2 <- nlapply(neurons, make_sotopo_dotprops)
+
+dps_aba21 <- nblast_allbyall(dps_list2, normalisation = "raw")
+dps_aba22 <- nblast_allbyall(dps_list2, UseAlpha = T, normalisation = "raw")
+
+par(mfrow = c(1,2))
+image(dps_aba21 / diag(dps_aba21), zlim= c(0,1), yaxt='n', xaxt='n')
+title("NBLAST")
+image(dps_aba22 / diag(dps_aba22), zlim= c(0,1), yaxt='n', xaxt='n')
+title("TNBLAST + SO")
+
+# ---------------------- artificial neurons
+nrn3x <- neurons[[3]]
+nrn3x$StartPoint <- endpoints(neurons[[3]])[length(endpoints(neurons[[3]]))]
+
+nrn3x <- neurons[[3]]
+nrn3x$StartPoint <- endpoints(neurons[[3]])[106]
+nrn4x <- neurons[[4]]
+nrn4x$StartPoint <- endpoints(neurons[[4]])[length(endpoints(neurons[[4]]))]
+
+newnrns <- neuronlist(neurons[[1]], neurons[[2]], nrn3x, nrn4x)
+
+plot(newnrns,soma = T)
+
+dps_listT1 <- nlapply(newnrns, make_topo_dotprops)
+
+dps_aba1 <- nblast_allbyall(dps_listT1, normalisation = "raw")
+dps_aba2 <- nblast_allbyall(dps_listT1, UseAlpha = T, normalisation = "raw")
+
+par(mfrow = c(1,2))
+image(dps_aba1 / diag(dps_aba1), zlim= c(0,1), yaxt='n', xaxt='n')
+title("NBLAST")
+image(dps_aba2 / diag(dps_aba2), zlim= c(0,1), yaxt='n', xaxt='n')
+title("TNBLAST")
+
+dps_listT2 <- nlapply(newnrns, make_sotopo_dotprops)
+
+dps_aba3 <- nblast_allbyall(dps_listT2, UseAlpha = T, normalisation = "raw")
+
+par(mfrow = c(1,2))
+image(dps_aba1 / diag(dps_aba1), zlim= c(0,1), yaxt='n', xaxt='n')
+title("NBLAST")
+image(dps_aba3 / diag(dps_aba2), zlim= c(0,1), yaxt='n', xaxt='n')
+title("TNBLAST + SO")

workbook_fw.R

@@ -0,0 +1,51 @@
+library(nat)
+library(fafbseg)
+library(dendroextras)
+
+devtools::load_all()
+
+OUT_DIR <- "meshes"
+
+fw_pair <- flywire_fetch("https://globalv1.flywire-daf.com/nglstate/5478428614590464", return="text")
+fw_scene <- flywire_fetch("https://globalv1.flywire-daf.com/nglstate/6012751403024384", return="text")
+
+fw_nids1 <- ngl_segments(ngl_decode_scene(fw_pair), as_character = TRUE)
+fw_nids2 <- ngl_segments(ngl_decode_scene(fw_scene), as_character = TRUE)
+#save_cloudvolume_meshes(fw_nids, savedir = OUT_DIR, format = 'obj')
+
+nids <- c(fw_nids1, fw_nids2)
+
+fw_nrns <- skeletor(nids, method = "wavefront")
+
+saveRDS(fw_nrns, "fw_neurons.rds")
+fw_nrns <- readRDS("fw_neurons.rds")
+
+fw_nrns_so <- nlapply(fw_nrns, function(x) reroot_hairball(x, k.soma.search = 50, radius.soma.search = 2500))
+
+fw_nrns_cl <- nlapply(fw_nrns_so/1e3, stitch_neurons_mst, OmitFailures=T) %>%
+  prune_twigs(twig_length=2, OmitFailures=T)
+
+fw_nrns_sp <- nlapply(fw_nrns_cl, function(x) simplify_neuron(x, n=2))
+
+plot3d(fw_nrns_sp, soma = T)
+
+dps_list <- nlapply(fw_nrns_sp, make_topo_dotprops)
+
+dps_aba1 <- nblast_allbyall(dps_list, normalisation = "normalised")
+dps_aba2 <- nblast_allbyall(dps_list, UseAlpha = T, normalisation = "normalised")
+
+par(mfrow = c(1,2))
+image(dps_aba1, zlim= c(-1,1), yaxt='n', xaxt='n')
+title("NBLAST")
+image(dps_aba2, zlim= c(-1,1), yaxt='n', xaxt='n')
+title("TNBLAST")
+
+hckcs1 <- nhclust(scoremat=dps_aba1)
+dev.off()
+dkcs1 <- colour_clusters(hckcs1, k=3)
+plot(dkcs1)
+
+hckcs2 <- nhclust(scoremat=dps_aba2)
+dkcs2 <- colour_clusters(hckcs2, k=3)
+plot(dkcs2)
+