diff --git a/DESCRIPTION b/DESCRIPTION index 3e9be63..7ce476a 100644 --- a/DESCRIPTION +++ b/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) ). @@ -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 diff --git a/R/neuriteblast.r b/R/neuriteblast.r index 775b4c2..1a91137 100644 --- a/R/neuriteblast.r +++ b/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 diff --git a/man/WeightedNNBasedLinesetMatching.Rd b/man/WeightedNNBasedLinesetMatching.Rd index 55e0b78..21c5483 100644 --- a/man/WeightedNNBasedLinesetMatching.Rd +++ b/man/WeightedNNBasedLinesetMatching.Rd @@ -8,12 +8,19 @@ \usage{ WeightedNNBasedLinesetMatching(target, query, ...) -\method{WeightedNNBasedLinesetMatching}{dotprops}(target, query, UseAlpha = FALSE, ...) +\method{WeightedNNBasedLinesetMatching}{dotprops}( + target, + query, + UseAlpha = FALSE, + UseTopo = FALSE, + ... +) \method{WeightedNNBasedLinesetMatching}{neuron}( target, query, UseAlpha = FALSE, + UseTopo = FALSE, OnlyClosestPoints = FALSE, ... ) @@ -26,6 +33,9 @@ compare (must be of the same class)} \item{UseAlpha}{Whether to scale dot product of tangent vectors (default=F)} +\item{UseTopo}{Whether to use topological information to scale dot products +of tangent vectors (default=F)} + \item{OnlyClosestPoints}{Whether to restrict searches to the closest points in the target (default FALSE, only implemented for \code{dotprops}).} } diff --git a/man/nblast.Rd b/man/nblast.Rd index 74a7aba..0fcfd46 100644 --- a/man/nblast.Rd +++ b/man/nblast.Rd @@ -12,6 +12,7 @@ nblast( version = c(2, 1), normalised = FALSE, UseAlpha = FALSE, + UseTopo = FALSE, OmitFailures = NA, ... ) @@ -38,6 +39,9 @@ query} segment to emphasise long range neurites rather then arbours (default: FALSE, see \bold{\code{UseAlpha}} section for details).} +\item{UseTopo}{whether to use directed dotprops vectors (pointing towards +soma) and features in neuron nodes (default: FALSE).} + \item{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 diff --git a/tests/testthat/test-NBLAST1.r b/tests/testthat/test-NBLAST1.r index 3be6293..0e25d72 100644 --- a/tests/testthat/test-NBLAST1.r +++ 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) +}) diff --git a/tests/testthat/test-NBLAST2.r b/tests/testthat/test-NBLAST2.r index 0b5b849..bdea92b 100644 --- a/tests/testthat/test-NBLAST2.r +++ b/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) +}) diff --git a/tests/testthat/testdata/testtopodps.rds b/tests/testthat/testdata/testtopodps.rds new file mode 100644 index 0000000..d09bf61 Binary files /dev/null and b/tests/testthat/testdata/testtopodps.rds differ diff --git a/utils.R b/utils.R new file mode 100644 index 0000000..d24d136 --- /dev/null +++ b/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 +} diff --git a/workbook.R b/workbook.R new file mode 100644 index 0000000..2438947 --- /dev/null +++ b/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") diff --git a/workbook_fw.R b/workbook_fw.R new file mode 100644 index 0000000..dc57151 --- /dev/null +++ b/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) + diff --git a/workbook_hemi.R b/workbook_hemi.R new file mode 100644 index 0000000..1bfd1ca --- /dev/null +++ b/workbook_hemi.R @@ -0,0 +1,126 @@ +library(nat) +library(fafbseg) +library(hemibrainr) +library(dendroextras) + +SKEL_DIR <- "/Users/dominik/projects/nblast_scoring_data/datasets/hemibrain_mismatch" +hemi_nrns <- read.neurons(SKEL_DIR) +plot3d(hemi_nrns, soma=T) + +sc_pairs <- read.csv(paste(SKEL_DIR, "../hemibrain_mismatch_pairs.csv", sep="/")) + +hemi_types <- unlist(lapply(hemi_nrns, function(nn) strsplit(nn$NeuronName,"_")[[1]][[2]])) +names(hemi_nrns) <- nlapply(hemi_nrns, function(nn) strsplit(nn$NeuronName,"_")[[1]][[3]]) + +hemi_nrns_raw <- hemi_nrns +hemi_nrns <- hemi_nrns/125 +hemi_nrns <- nlapply(hemi_nrns,stitch_neurons_mst) + +htst <- make_topo_dotprops(hemi_nrns[[1]], resample = 1, k=5) +plot3d(htst) + +library(doMC) +registerDoMC(cores=4) +hemi_dps <- nlapply(hemi_nrns, make_topo_dotprops) +hemi_dps_so <- nlapply(hemi_nrns, make_sotopo_dotprops) + +saveRDS(hemi_dps, "hemi_dotprops.rds") +hemi_dps <- readRDS("hemi_dotprops.rds") + +dps_aba1 <- nblast_allbyall(hemi_dps, normalisation = "normalised") +dps_aba2 <- nblast_allbyall(hemi_dps, UseAlpha = T, normalisation = "normalised") +dps_aba3 <- nblast_allbyall(hemi_dps_so, UseAlpha = T, normalisation = "normalised") + +par(mfrow = c(1,3)) +image(dps_aba1, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("NBLAST") +image(dps_aba2, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("TNBLAST") +image(dps_aba3, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("TNBLAST+SO") +dev.off() + +hcpn1 <- nhclust(scoremat=dps_aba1) +hcpn2 <- nhclust(scoremat=dps_aba2) +hcpn3 <- nhclust(scoremat=dps_aba3) + +plot3d(hemi_nrns["5813027132"], soma=T) +plot3d(hemi_nrns["2009184596"], soma=T, add=T, color="blue") + +val1 <- c() +val2 <- c() +val3 <- c() +for (i in 1:dim(sc_pairs)[[1]]) { + sc1 <- dps_aba1[as.character(sc_pairs[i,]$query_id), as.character(sc_pairs[i,]$target_id)] + sc2 <- dps_aba2[as.character(sc_pairs[i,]$query_id), as.character(sc_pairs[i,]$target_id)] + sc3 <- dps_aba3[as.character(sc_pairs[i,]$query_id), as.character(sc_pairs[i,]$target_id)] + val1 <- c(val1, sc1) + val2 <- c(val2, sc2) + val3 <- c(val3,sc3) +} + +ii<-5 +plot3d(hemi_nrns[as.character(sc_pairs[ii,]$query_id)], soma=T) +plot3d(hemi_nrns[as.character(sc_pairs[ii,]$target_id)], soma=T, add=T, color="blue") + + +tanglegram(dkcs1,dkcs2) + +library('fossil') +library('clValid') + +k = length(unique(hemi_types)) + +cutk1 <- cutree(hcpn1, k = k) +rand.index(cutk1, encode_ordinal(hemi_types)) + +cutk2 <- cutree(hcpn2, k = k) +rand.index(cutk2, encode_ordinal(hemi_types)) + +cutk3 <- cutree(hcpn3, k = k) +rand.index(cutk3, encode_ordinal(hemi_types)) + +# test on neurons represented in multiple types + +nrns_types <- data.frame(type=hemi_types) +multi_types <- nrns_types %>% group_by(type) %>% count() %>% filter(n>1) + +m_hemi_dps <- hemi_dps[nrns_types$type %in% multi_types$type] +m_hemi_dps_so <- hemi_dps_so[nrns_types$type %in% multi_types$type] + +nrns_types_m <- hemi_types[nrns_types$type %in% multi_types$type] + +dps_aba1 <- nblast_allbyall(m_hemi_dps, normalisation = "normalised") +dps_aba2 <- nblast_allbyall(m_hemi_dps, UseAlpha = T, normalisation = "normalised") +dps_aba3 <- nblast_allbyall(m_hemi_dps_so, UseAlpha = T, normalisation = "normalised") + +hc1 <- nhclust(scoremat=dps_aba1) +hc2 <- nhclust(scoremat=dps_aba2) +hc3 <- nhclust(scoremat=dps_aba3) + +k = length(multi_types$type) + +dkcs1 = colour_clusters(hc1, k=k) +labels(dkcs1) <- nrns_types_m +plot(dkcs1) +dkcs2 <- colour_clusters(hc2, k=k) +labels(dkcs2) <- nrns_types_m +plot(dkcs2) +dkcs3 <- colour_clusters(hc3, k=k) +labels(dkcs3) <- nrns_types_m +plot(dkcs3) + +cutk1 <- cutree(hc1, k = k) +rand.index(cutk1, encode_ordinal(nrns_types_m)) + +cutk2 <- cutree(hc2, k = k) +rand.index(cutk2, encode_ordinal(nrns_types_m)) + +cutk3 <- cutree(hc3, k = k) +rand.index(cutk3, encode_ordinal(nrns_types_m)) + +dunn((dps_aba1-1)*-1, cutk1) + +dunn((dps_aba2-1)*-1, cutk2) + +dunn((dps_aba3-1)*-1, cutk3) diff --git a/workbook_swc.R b/workbook_swc.R new file mode 100644 index 0000000..d616396 --- /dev/null +++ b/workbook_swc.R @@ -0,0 +1,81 @@ +library(nat) +library(fafbseg) +library(dendroextras) +library(dendextend) +library(elmr) +library('fossil') +library('clValid') + +devtools::load_all() + +SKEL_DIR <- "/Users/dominik/projects/nblast_scoring_data/datasets/FAFB_RHS_uPNs" +nrns <- read.neurons(SKEL_DIR) +plot3d(nrns, soma=T) + +lin_type <- unlist(nlapply(nrns, function(n) strsplit(n$NeuronName, "PN_")[[1]][[1]])) + +nrns <- nrns/1e3 +nrns_s <- nlapply(nrns, function(x) elmr::unspike(x, threshold=4)) +nrns_s <- nlapply(nrns_s, function(x) smooth_neuron(x, sigma=0.8)) +nrns_cl <- prune_twigs(nrns_s, twig_length=5, OmitFailures=T) + +plot3d(nrns_cl, soma=T) + +dps_list <- nlapply(nrns_cl, make_topo_dotprops) + +plot3d(dps_list) + +dps_aba1 <- nblast_allbyall(dps_list, normalisation = "normalised") +dps_aba2 <- nblast_allbyall(dps_list, UseAlpha = T, normalisation = "normalised") + +dps_list_so <- nlapply(nrns_cl, make_sotopo_dotprops) + +dps_aba3 <- nblast_allbyall(dps_list_so, UseAlpha = T, normalisation = "normalised") + +par(mfrow = c(1,2)) +image(dps_aba1, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("NBLAST") +image(dps_aba2, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("TNBLAST") + +par(mfrow = c(1,2)) +image(dps_aba1, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("NBLAST") +image(dps_aba3, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("TNBLAST+SO") +dev.off() + +hcpn1 <- nhclust(scoremat=dps_aba1) +hcpn2 <- nhclust(scoremat=dps_aba2) +hcpn3 <- nhclust(scoremat=dps_aba3) + +k = length(unique(lin_type)) + +dkcs1 = colour_clusters(hcpn1, k=k) +labels(dkcs1) <- lin_type +plot(dkcs1) + +plot3d(hcpn1, k=k, db=dps_list) + +dkcs2 = colour_clusters(hcpn2, k=k) +labels(dkcs2) <- lin_type +plot(dkcs2) + +dkcs3 = colour_clusters(hcpn3, k=k) +labels(dkcs3) <- lin_type +plot(dkcs3) + +cutk1 <- cutree(hcpn1, k = k) +rand.index(cutk1, encode_ordinal(lin_type)) + +cutk2 <- cutree(hcpn2, k = k) +rand.index(cutk2, encode_ordinal(lin_type)) + +cutk3 <- cutree(hcpn3, k = k) +rand.index(cutk3, encode_ordinal(lin_type)) + +dunn((dps_aba1-1)*-1, cutk1) + +dunn((dps_aba2-1)*-1, cutk2) + +dunn((dps_aba3-1)*-1, cutk3)