diff --git a/R/neuriteblast.r b/R/neuriteblast.r
index 7e990c1..ab83657 100644
--- a/R/neuriteblast.r
+++ b/R/neuriteblast.r
@@ -289,8 +289,13 @@ NeuriteBlast <- function(query, target, targetBinds=NULL, normalised=FALSE,
   if(nat::is.neuronlist(query)) {
     res=plyr::llply(query, NeuriteBlast, target=target, targetBinds=targetBinds,
                   normalised=normalised, OmitFailures=OmitFailures, simplify=simplify, ...=...)
-    if (!identical(simplify, FALSE) && length(res))
+    if (!identical(simplify, FALSE) && length(res)){
+      if(normalised)
+        selfscores=sapply(res, "attr", "scaled:scale")
       res=simplify2array(res, higher = (simplify == "array"))
+      if(normalised)
+        attr(res,'scaled:scale')=selfscores
+    }
     return(res)
   } else {
     if(is.null(targetBinds))
@@ -315,8 +320,12 @@ NeuriteBlast <- function(query, target, targetBinds=NULL, normalised=FALSE,
   }
 
   if(normalised){
-    if(is.list(scores)) stop("Cannot normalise results when they are not a single number")
-    scores=scores/NeuriteBlast(query, neuronlist(query), normalised=FALSE, ...)
+    if(is.list(scores))
+      stop("Cannot normalise results when they are not a single number")
+    # nb roll our own scale, since scale turns vector into matrix
+    selfscore=NeuriteBlast(query, neuronlist(query), normalised=FALSE, ...)
+    scores=scores/selfscore
+    attr(scores,'scaled:scale')=selfscore
   }
   scores
 }
diff --git a/tests/testthat/test-NBLAST2.r b/tests/testthat/test-NBLAST2.r
index e7054cf..0b5b849 100644
--- a/tests/testthat/test-NBLAST2.r
+++ b/tests/testthat/test-NBLAST2.r
@@ -21,15 +21,49 @@ test_that("nblast v2 handles a neuronlist as query", {
   expect_equal(scores, scores.expected)
 })
 
+test_that("we can calculate scores for non square matrices", {
+  expect_equal(scoresaba <- nblast_allbyall(testneurons, version=2),
+               nblast(testneurons, testneurons, version=2))
+  expect_equal(nblast(testneurons[1:2], testneurons[3:5], version=2),
+               scoresaba[3:5, 1:2])
+})
+
 test_that("we can calculate normalised nblast v2 scores", {
-  scores <- nblast(testneurons[1:2], testneurons, version=2)
-  scores.norm=scale(scores, scale=c(scores[1,1],scores[2,2]), center = FALSE)
+  expect_is(scores <- nblast(testneurons[1:2], testneurons, version=2), 'matrix')
+
+  expect_is(scores.norm <-
+              scale(scores, scale=c(scores[1,1],scores[2,2]), center = FALSE),
+            'matrix')
   expect_equivalent(nblast(testneurons[1:2], testneurons, normalised=TRUE), scores.norm)
+
+  expect_equal(scoresaba <- nblast_allbyall(testneurons,
+                                            version=2,
+                                            normalisation = 'normalised'),
+               nblast(testneurons, testneurons, version=2, normalised = TRUE))
+
+  scores.rect <- nblast(testneurons[1:2], testneurons[3:5], version = 2)
+  scores.rect.norm <-
+    nblast(testneurons[1:2],
+           testneurons[3:5],
+           version = 2,
+           normalised = T)
+  expect_equal(scores.rect.norm,
+                    scale(
+                      scores.rect,
+                      center = FALSE,
+                      scale = attr(scores.rect.norm, 'scaled:scale')
+                    ))
+
+  # the simple subset will drop the scaled attribute
+  baseline=scores.rect.norm[, 1, drop = FALSE]
+  attr(baseline, 'scaled:scale')=attr(scores.rect.norm, 'scaled:scale')[1]
+  expect_equal(nblast(testneurons[1], testneurons[3:5], version = 2, normalised = T),
+               baseline)
+
 })
 
 test_that("we can calculate scores for regular neurons",{
-  library(nat)
-  nblast_allbyall(Cell07PNs[1:4])
+  expect_is(nblast_allbyall(Cell07PNs[1:4]), 'matrix')
 })
 
 test_that("we can calculate scores using getOption('nat.default.neuronlist')", {
@@ -57,13 +91,12 @@ test_that("we can handle OmitFailures", {
 })
 
 test_that("we can handle all combinations of dotprops and neurons, both as neuronlists and singly", {
-  nblast(testneurons[[1]], Cell07PNs[1:3])
-  nblast(testneurons[1:3], Cell07PNs[1:3])
-  nblast(testneurons[[1]], Cell07PNs[[1]])
-  nblast(testneurons[1:3], Cell07PNs[[1]])
-
-  nblast(Cell07PNs[[1]], testneurons[1:3])
-  nblast(Cell07PNs[1:3], testneurons[1:3])
-  nblast(Cell07PNs[[1]], testneurons[[1]])
-  nblast(Cell07PNs[1:3], testneurons[[1]])
+  expect_is(nblast(testneurons[[1]], Cell07PNs[1:3]), 'numeric')
+  expect_is(nblast(testneurons[1:3], Cell07PNs[1:3]), 'matrix')
+  expect_is(nblast(testneurons[[1]], Cell07PNs[[1]]), 'numeric')
+  expect_is(nblast(testneurons[1:3], Cell07PNs[[1]]), 'numeric')
+  expect_is(nblast(Cell07PNs[[1]], testneurons[1:3]), 'numeric')
+  expect_is(nblast(Cell07PNs[1:3], testneurons[1:3]), 'matrix')
+  expect_is(nblast(Cell07PNs[[1]], testneurons[[1]]), 'numeric')
+  expect_is(nblast(Cell07PNs[1:3], testneurons[[1]]), 'numeric')
 })