topoNBLAST first implementation added

DESCRIPTION

@@ -7,7 +7,9 @@ Authors@R: c(
     person("Gregory", "Jefferis", email= "jefferis@gmail.com",
       role = c("aut"), comment = c(ORCID = "0000-0002-0587-9355")),
     person("James", "Manton", email="ajd.manton@googlemail.com",
-      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>).

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 weight the similarity score for each matched
+#'   segment to relative topology of the neurons (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, sd=sd, normalised=normalised, UseTopo=UseTopo,
                  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,17 @@ 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,
+                                   topo1=target$topo,topo2=query$topo,...)
+  } else
     WeightedNNBasedLinesetMatching(target$points,query$points,dvs1=target$vect,dvs2=query$vect,...)
 }
 
@@ -396,10 +406,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,
+                                          UseTopo=UseTopo, UseAlpha=UseAlpha,
+                                          OnlyClosestPoints=OnlyClosestPoints, ...))
   }
   if(UseAlpha)
     stop("UseAlpha is not yet implemented for neurons!")
@@ -413,7 +426,8 @@ WeightedNNBasedLinesetMatching.neuron<-function(target, query, UseAlpha=FALSE,
 
 WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NULL,alphas1=NULL,
                                         alphas2=NULL,NNDistFun=WeightedNNBasedLinesetDistFun,Verbose=FALSE,
-                                        BothDirections=FALSE,BothDirectionsFun=list,OnlyClosestPoints=FALSE,...){
+                                        BothDirections=FALSE,BothDirectionsFun=list,OnlyClosestPoints=FALSE,
+                                        topo1=NULL,topo2=NULL,...){
   # return a score based on the similarity of nearest neighbour location
   # and the dot product of the direction vectors
 
@@ -476,6 +490,19 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL
       scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]]))
       dps=dps*scalefac
     }
+    if(!is.null(topo1) && !is.null(topo2)){
+      # use of local neuron topology
+      ## distance from body
+      maxlen = max(max(topo1$distance), max(topo2$distance))
+      scalefact_dist = zapsmall(1 - (abs(topo1$distance[idxArray[,1]] - topo2$distance[idxArray[,2]])/maxlen)^2)
+      ## strengthen by reversed Strahler's Order
+      max_rso = max(max(topo1$rso), max(topo2$rso))
+      if (max_rso != 0) {
+        scalefac_rso = zapsmall(1 - (abs(topo1$rso[idxArray[,1]] - topo2$rso[idxArray[,2]])/max_rso)^0.5)
+      } else scalefac_rso = 1
+      scalefac = scalefact_dist * scalefac_rso
+      dps=dps*scalefac
+    }
   }
 
   NNDistFun(as.vector(nntarget$nn.dists),dps,...)

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

@@ -2,6 +2,7 @@ context("NBLAST v2")
 
 testneurons <- readRDS('testdata/testneurons.rds')
 
+
 test_that("nblast v2 produces expected scores", {
   scores <- nblast(testneurons[[1]], testneurons, version=2)
   scores.expected <- structure(c(43518.2468824115, -34509.2778341796, -36608.1290149345, -18366.0246697251, -36782.7167494325), .Names = c("5HT1bMARCM-F000001_seg001", "5HT1bMARCM-F000002_seg001", "5HT1bMARCM-F000003_seg001", "5HT1bMARCM-F000004_seg001", "5HT1bMARCM-F000005_seg001"))
@@ -100,3 +101,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)
+})