diff --git a/R/calcNhoodDistance.R b/R/calcNhoodDistance.R
index 0bb5e29..5db7782 100644
--- a/R/calcNhoodDistance.R
+++ b/R/calcNhoodDistance.R
@@ -42,7 +42,7 @@
 #' @importFrom irlba prcomp_irlba
 #' @importFrom SummarizedExperiment assay
 #' @importFrom Matrix which
-calcNhoodDistance <- function(x, d, reduced.dim=NULL, use.assay="logcounts"){
+calcNhoodDistance <- function(x, d, reduced.dim=NULL, use.assay="logcounts", ncores=1L){
     if(is(x, "Milo")){
         # check for reducedDims
         if((length(reducedDimNames(x)) == 0) & is.null(reduced.dim)){
@@ -58,33 +58,43 @@ calcNhoodDistance <- function(x, d, reduced.dim=NULL, use.assay="logcounts"){
         stop("Input is not a valid Milo object")
     }
 
-    non.zero.nhoods <- which(nhoods(x)!=0, arr.ind = TRUE)
+    nhood_matrix <- nhoods(x)
+    non.zero.nhoods <- which(nhood_matrix != 0, arr.ind = TRUE)
 
+    if(ncores > 1){
+        if (Sys.info()['sysname'] == "Windows") 
+          bpp <- BiocParallel::SnowParam(workers=ncores)
+        else
+          bpp <- BiocParallel::MulticoreParam(workers=ncores)
+    }else{
+         bpp <- BiocParallel::SerialParam()
+    }
+    
     if(is.character(reduced.dim)){
         # check if it exists in the slot
         if(!any(names(reducedDims(x)) %in% reduced.dim)){
             stop(reduced.dim, " not found in the reducedDim slot")
         }
-        nhood.dists <- sapply(seq_len(ncol(nhoods(x))),
-                              function(X) .calc_distance(reducedDim(x, reduced.dim)[non.zero.nhoods[non.zero.nhoods[,'col']==X,'row'],
-                                                                                    seq_len(d),drop=FALSE]))
+        nhood.dists <- BiocParallel::bplapply(seq_len(ncol(nhood_matrix)),
+                              function(X) as(proxyC::dist(reducedDim(x, reduced.dim)[non.zero.nhoods[non.zero.nhoods[,'col']==X,'row'],
+                                                                                    seq_len(d),drop=FALSE]), "dsCMatrix"), BPPARAM = bpp)
         names(nhood.dists) <- nhoodIndex(x)
     } else if(is(reduced.dim, "matrix")){
-        nhood.dists <- sapply(seq_len(ncol(nhoods(x))),
-                              function(X) .calc_distance(reduced.dim[non.zero.nhoods[non.zero.nhoods[,'col']==X,'row'],
-                                                                     seq_len(d),drop=FALSE]))
+        nhood.dists <- BiocParallel::bplapply(seq_len(ncol(nhood_matrix)),
+                              function(X) as(proxyC::dist(reduced.dim[non.zero.nhoods[non.zero.nhoods[,'col']==X,'row'],
+                                                                     seq_len(d),drop=FALSE]), "dsCMatrix"), BPPARAM = bpp)
     } else if(is.null(reduced.dim)){
         if(any(names(reducedDims(x)) %in% c("PCA"))){
-            nhood.dists <- sapply(seq_len(ncol(nhoods(x))),
-                                  function(X) .calc_distance(reducedDim(x, "PCA")[non.zero.nhoods[non.zero.nhoods[,'col']==X,'row'],
-                                                                                  seq_len(d),drop=FALSE]))
+            nhood.dists <- BiocParallel::bplapply(seq_len(ncol(nhood_matrix)),
+                               function(X) as(proxyC::dist(reducedDim(x, "PCA")[non.zero.nhoods[non.zero.nhoods[,'col']==X,'row'],
+                                                                                  seq_len(d),drop=FALSE]), "dsCMatrix"), BPPARAM = bpp)
             names(nhood.dists) <- nhoodIndex(x)
 
         } else{
             stop("No reduced.dim slot specified")
         }
     }
-
+                                                  
     nhoodDistances(x) <- nhood.dists
 
     return(x)
diff --git a/R/findNhoodGroupMarkers.R b/R/findNhoodGroupMarkers.R
index 4d207d3..2871b85 100644
--- a/R/findNhoodGroupMarkers.R
+++ b/R/findNhoodGroupMarkers.R
@@ -159,28 +159,15 @@ findNhoodGroupMarkers <- function(x, da.res, assay="logcounts",
     fake.meta[,"sample_id"] <- colData(x)[fake.meta$CellID,sample_col]
     fake.meta[,'sample_group'] <- paste(fake.meta[,"sample_id"], fake.meta[,"Nhood.Group"], sep="_")
 
-    sample_gr_mat <- matrix(0, nrow=nrow(fake.meta), ncol=length(unique(fake.meta$sample_group)))
-    colnames(sample_gr_mat) <- unique(fake.meta$sample_group)
+    sample_gr_mat <- model.matrix(~ 0 + sample_group, data=fake.meta)
+    colnames(sample_gr_mat) <- gsub("sample_group", "", colnames(sample_gr_mat))
     rownames(sample_gr_mat) <- rownames(fake.meta)
 
-    for (s in colnames(sample_gr_mat)) {
-      sample_gr_mat[which(fake.meta$sample_group == s),s] <- 1
-    }
-
     ## Summarise expression by sample
-    exprs_smp <- matrix(0, nrow=nrow(exprs), ncol=ncol(sample_gr_mat))
-    if (assay=='counts') {
-      summFunc <- rowSums
-    } else {
-      summFunc <- rowMeans
-    }
-
-    for (i in seq_len(ncol(sample_gr_mat))){
-      if (sum(sample_gr_mat[,i]) > 1) {
-        exprs_smp[,i] <- summFunc(exprs[,which(sample_gr_mat[,i] > 0)])
-      } else {
-        exprs_smp[,i] <- exprs[,which(sample_gr_mat[,i] > 0)]
-      }
+    exprs_smp <- exprs %*% sample_gr_mat
+    
+    if (assay != "counts"){
+      exprs_smp <- exprs_smp / rep(colSums(sample_gr_mat), each = nrow(exprs))
     }
     rownames(exprs_smp) <- rownames(exprs)
     colnames(exprs_smp) <- colnames(sample_gr_mat)
@@ -188,7 +175,7 @@ findNhoodGroupMarkers <- function(x, da.res, assay="logcounts",
     smp_meta <- unique(fake.meta[,c("sample_group","Nhood.Group")])
     rownames(smp_meta) <- smp_meta[,"sample_group"]
 
-    fake.meta <- smp_meta
+    fake.meta <- smp_meta[colnames(exprs_smp),]
     exprs <- exprs_smp
   }
 
@@ -198,7 +185,7 @@ findNhoodGroupMarkers <- function(x, da.res, assay="logcounts",
     nhood.gr <- subset.groups
   }
 
-  for(i in seq_along(nhood.gr)){
+  for(i in seq_len(min(length(unique(fake.meta$Nhood.Group)), length(unique(nhood.gr))))){
     i.meta <- fake.meta
     i.meta$Test <- "Ref"
     i.meta$Test[fake.meta$Nhood.Group == nhood.gr[i]] <- "Test"
diff --git a/R/plotNhoods.R b/R/plotNhoods.R
index 4d8e273..410eaee 100644
--- a/R/plotNhoods.R
+++ b/R/plotNhoods.R
@@ -213,7 +213,8 @@ plotNhoodGraph <- function(x, layout="UMAP", colour_by=NA, subset.nhoods=NULL, s
   }
 
   if (is.numeric(V(nh_graph)$colour_by)) {
-    pl <- pl + scale_fill_gradient2(name=colour_by)
+    pl <- pl + scale_fill_gradient2(name=colour_by, low = scales::muted("blue"),
+                     mid = "white", high = scales::muted("red"))
   } else {
       if(is.factor(V(nh_graph)$colour_by)){
           mycolors <- colorRampPalette(brewer.pal(11, "Spectral"))(length(levels(V(nh_graph)$colour_by)))
@@ -734,39 +735,20 @@ plotDAbeeswarm <- function(da.res, group.by=NULL, alpha=0.1, subset.nhoods=NULL)
     da.res <- da.res[subset.nhoods,]
   }
 
-  # Get position with ggbeeswarm
-  beeswarm_pos <- ggplot_build(
-    da.res %>%
-      mutate(is_signif = ifelse(SpatialFDR < alpha, 1, 0)) %>%
-      arrange(group_by) %>%
-      ggplot(aes(group_by, logFC)) +
-      geom_quasirandom()
-  )
-
-  pos_x <- beeswarm_pos$data[[1]]$x
-  pos_y <- beeswarm_pos$data[[1]]$y
-
-  n_groups <- unique(da.res$group_by) %>% length()
-
-  da.res %>%
-    mutate(is_signif = ifelse(SpatialFDR < alpha, 1, 0)) %>%
-    mutate(logFC_color = ifelse(is_signif==1, logFC, NA)) %>%
-    arrange(group_by) %>%
-    mutate(Nhood=factor(Nhood, levels=unique(Nhood))) %>%
-    mutate(pos_x = pos_x, pos_y=pos_y) %>%
-    ggplot(aes(pos_x, pos_y, color=logFC_color)) +
-    scale_color_gradient2() +
+  da.res %>% dplyr::mutate(logFC_color = ifelse(SpatialFDR < alpha, logFC, NA)) %>% 
+    dplyr::arrange(group_by, !is.na(logFC_color)) %>% 
+    ggplot(aes(group_by, logFC, color=logFC_color)) + 
+    geom_hline(yintercept=0, linewidth = 0.2, color = "black") +
+    ggbeeswarm::geom_quasirandom() +
+    scale_color_gradient2(low  = scales::muted("blue"), mid  = "grey80", 
+                          high = scales::muted("red"), na.value = "grey80") +
     guides(color="none") +
-    xlab(group.by) + ylab("Log Fold Change") +
-    scale_x_continuous(
-      breaks = seq(1,n_groups),
-      labels = setNames(levels(da.res$group_by), seq(1,n_groups))
-      ) +
-    geom_point() +
-    coord_flip() +
-    theme_bw(base_size=22) +
-    theme(strip.text.y =  element_text(angle=0))
-
+    xlab(NULL) +
+    ylab("Log Fold Change") +
+    geom_boxplot(alpha=0,outlier.shape=NA,width=0.25) +
+    theme_bw(base_size = 12) +
+    theme(axis.text = element_text(color = "black"),
+          axis.text.x = element_text(angle=90, hjust = 1, vjust = 0.5))
 }
 
 
diff --git a/R/testNhoods.R b/R/testNhoods.R
index c943fcd..6ebc1f5 100644
--- a/R/testNhoods.R
+++ b/R/testNhoods.R
@@ -572,7 +572,8 @@ testNhoods <- function(x, design, design.df, kinship=NULL,
                                 "SE"= unlist(lapply(lapply(fit, `[[`, "SE"), function(x) x[ret.beta])),
                                 "tvalue" = unlist(lapply(lapply(fit, `[[`, "t"), function(x) x[ret.beta])),
                                 "PValue" = unlist(lapply(lapply(fit, `[[`, "PVALS"), function(x) x[ret.beta])),
-                                matrix(unlist(lapply(fit, `[[`, "Sigma")), ncol=length(rand.levels), byrow=TRUE),
+                                as.matrix(t(do.call(cbind, lapply(fit, `[[`, "Sigma")))),
+                                # matrix(unlist(lapply(fit, `[[`, "Sigma")), ncol=length(rand.levels), byrow=TRUE),
                                 "Converged"=unlist(lapply(fit, `[[`, "converged")), "Dispersion" = unlist(lapply(fit, `[[`, "Dispersion")),
                                 "Logliklihood"=unlist(lapply(fit, `[[`, "LOGLIHOOD")))
 
diff --git a/R/utils.R b/R/utils.R
index 6a02bf1..a13ad9e 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -134,8 +134,12 @@
 #' @importFrom Matrix crossprod
 .build_nhood_adjacency <- function(nhoods, overlap=1){
     nh_intersect_mat <- Matrix::crossprod(nhoods)
-    nh_intersect_mat[nh_intersect_mat < overlap] <- 0
-
+    if(inherits(nh_intersect_mat, c("sparseMatrix"))){
+        nh_intersect_mat@x[nh_intersect_mat@x < overlap] <- 0
+        nh_intersect_mat <- Matrix::drop0(nh_intersect_mat)
+    }else{
+        nh_intersect_mat[nh_intersect_mat < overlap] <- 0    
+    }
     rownames(nh_intersect_mat) <- colnames(nhoods)
     colnames(nh_intersect_mat) <- colnames(nhoods)
     return(nh_intersect_mat)