Minor addition and fixes to code related to transmitters, skeletonisation and CAVE tokens

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: fafbseg
 Title: Support Functions for Analysis of FAFB EM Segmentation
-Version: 0.15.1.9000
+Version: 0.15.2
 Authors@R: 
     c(person(given = "Gregory",
              family = "Jefferis",

NAMESPACE

@@ -7,6 +7,7 @@ S3method(as.character,ngscene)
 S3method(as.data.frame,ng_raw_list)
 S3method(as.mesh3d,ng_raw)
 S3method(as.mesh3d,ng_raw_list)
+S3method(as.neuron,navis.core.dotprop.Dotprops)
 S3method(as.neuron,navis.core.skeleton.TreeNeuron)
 S3method(as.neuronlist,navis.core.base.BaseNeuron)
 S3method(as.neuronlist,navis.core.neuronlist.NeuronList)

R/autosyn.R

@@ -698,8 +698,10 @@ flywire_ntpred <- function(x,
     ntpredictions=ntpredictions_tbl(local=local)
     if(is.null(ntpredictions))
       stop("I cannot find the neurotransmitter predictions sqlite database!")
-
     x = ntpredictions %>%
+      dplyr::select(id,
+                    gaba, acetylcholine, glutamate,
+                    dopamine, serotonin, octopamine) %>%
       dplyr::inner_join(x, copy = TRUE, by=c("id"="offset")) %>%
       dplyr::rename(offset="id")
   }
@@ -867,7 +869,7 @@ flywire_ntplot3d <- function(x, nts=c("gaba", "acetylcholine", "glutamate",
 #'   \code{\link{flywire_partners}}.
 #' @param remove_autapses  whether to remove autapses (defaults to \code{TRUE}).
 #' @param transmitters if \code{TRUE} also attempt to retrieve neurotransmitter
-#'   predictions from Eckstein et al. 2020, for the flywire neuron in question.
+#'   predictions from Eckstein and Bates et al. 2024, for the flywire neuron in question.
 #' @param cleft.threshold select only synaptic connections exceeding this
 #'   confidence threshold (default of 0 uses all synapses; values in the range
 #'   30-100 seem to make sense).
@@ -978,7 +980,7 @@ flywire_neurons_add_synapses.neuron <- function(x,
     # Add synapses
     wanted = c(intersect(c("offset", "prepost","scores", "cleft_scores",
     "segmentid_pre", "segmentid_post", "pre_svid", "post_svid",
-    "pre_id", "post_id"),colnames(synapses)), "x", "y", "z")
+    "pre_id", "post_id", "pre_x","pre_y","pre_z"),colnames(synapses)), "x", "y", "z")
     for(pos in c("pre_x","pre_y","pre_z","post_x","post_y","post_z")){
       if(!pos%in%colnames(synapses)){
         synapses[[pos]] = NA
@@ -999,25 +1001,44 @@ flywire_neurons_add_synapses.neuron <- function(x,
   }
   attr(synapses.xyz, "rootid") = rootid
   # If transmitters
+  pref.order = c("offset", "x", "y", "z", "scores", "cleft_scores",
+                 "top_p", "top_nt",
+                 "gaba", "acetylcholine", "glutamate", "octopamine", "serotonin", "dopamine",
+                 "prepost", "segmentid_pre", "segmentid_post",
+                 "pre_svid", "post_svid", "pre_id", "post_id")
   if(transmitters & nrow(synapses.xyz)){
     if(Verbose){
-      message("adding transmitter prediction information (Eckstein et al. 2020)")
+      message("Adding transmitter prediction information (Eckstein and Bates et al. 2024)")
+    }
+    npred = tryCatch(flywire_ntpred(x=synapses.xyz, local = local, cloudvolume.url = cloudvolume.url),
+                     error = function(e){
+                       warning(as.character(e))
+                       NULL
+                     })
+    if(is.null(npred)){
+      if(all(c("gaba", "acetylcholine","glutamate", "octopamine", "serotonin", "dopamine")%in%colnames(synapses))){
+        ntpred = synapses
+      }else{
+        warning('no transmitter data found')
+      }
     }
-    npred = flywire_ntpred(x=synapses.xyz, local = local, cloudvolume.url = cloudvolume.url)
-    pref.order = c("offset", "x", "y", "z", "scores", "cleft_scores", "top_p", "top_nt", "gaba", "acetylcholine",
-                   "glutamate", "octopamine", "serotonin", "dopamine", "prepost",
-                   "segmentid_pre", "segmentid_post",
-                   "pre_svid", "post_svid", "pre_id", "post_id")
-    pref.order = intersect(pref.order,colnames(npred))
-    if(nrow(npred)){
-      synapses.xyz = npred[,pref.order] %>%
-        dplyr::rename(syn_top_p = top_p,
-                      syn_top_nt = top_nt)
+    if(length(npred)){
+      npred %>%
+        dplyr::filter(.data$cleft_scores >= cleft.threshold) %>%
+        dplyr::arrange(.data$offset) %>%
+        as.data.frame() ->
+        synapses.xyz
     }
   }else if(nrow(synapses.xyz)){
-    synapses.xyz$syn_top_nt = "unknown"
-    synapses.xyz$syn_top_p = 0
+    synapses.xyz$top_nt = "unknown"
+    synapses.xyz$top_p = 0
+  }
+  pref.order.missing = setdiff(pref.order,colnames(synapses.xyz))
+  if(length(pref.order.missing)){
+    warning('missing synapse columns: ', paste(pref.order.missing,collapse=" "))
   }
+  pref.order = intersect(pref.order,colnames(synapses.xyz))
+  synapses.xyz = synapses.xyz[,pref.order]
   # Attach synapses to skeleton
   if(nrow(synapses.xyz)){
     nat::xyzmatrix(synapses.xyz) = fafb2flywire(nat::xyzmatrix(synapses.xyz))
@@ -1029,7 +1050,7 @@ flywire_neurons_add_synapses.neuron <- function(x,
       x$connectors[,colnames(x$connectors)%in%poss.nts] = round(x$connectors[,colnames(x$connectors)%in%poss.nts],digits=2)
     }
     # Get top transmitter result
-    tx=table(subset(synapses.xyz, synapses.xyz$prepost == 0)$syn_top_nt)
+    tx=table(subset(synapses.xyz, synapses.xyz$prepost == 0)$top_nt)
     tx=sort(tx, decreasing = TRUE)/sum(tx)*100
     if(length(tx)){
       x$ntpred = tx
@@ -1084,9 +1105,8 @@ flywire_neurons_add_synapses.neuronlist <- function(x,
     neurons.syn
   }
 }
-# neurons.syns = flywire_neurons_add_synapses(neurons, transmitters = TRUE, local =  "/Volumes/nnautilus/projects/JanFunke")
 
-# extract predictions neurons
+# extract predictions for neurons
 extract_ntpredictions.neuronlist <- function(x,
                                              poss.nts=c("gaba", "acetylcholine", "glutamate", "octopamine", "serotonin","dopamine")){
   nmeta = lapply(x,extract_ntpredictions.neuron,poss.nts=poss.nts)

R/cloudvolume.R

@@ -53,7 +53,6 @@ chunkedgraph_token <- function(url=NULL, cached=TRUE) {
   .chunkedgraph_token.memo(domain=domain)
 }
 
-
 #' Record (and if necessary create) a FlyWire chunkedgraph token
 #'
 #' Writes a token to a standard location on disk so that it will be found by the

R/fafbseg-py.R

@@ -30,6 +30,7 @@ as.neuronlist.navis.core.neuronlist.NeuronList <- function(l, df=NULL, ...) {
   nl
 }
 
+#' @export
 as.neuron.navis.core.dotprop.Dotprops <- function(x, ...) {
   y=x$convert_units('micron')
   dp=nat::as.dotprops(list(points=y$points, alpha=NULL, vect=y$vect))

R/flywire-fetch.R

@@ -20,6 +20,7 @@
 #'   token.
 #' @param config (optional) curl options, see \code{httr::\link[httr]{config}}
 #'   for details.
+#' @param domain the domain name for which your CAVE token is valid, defaults to `url`.
 #'
 #' @return Either an R object based on parsing returned JSON, a character vector
 #'   containing the raw JSON or a \code{httr::\link[httr]{response}} object,
@@ -42,6 +43,7 @@ flywire_fetch <- function(url,
                           retry = 0L,
                           include_headers = FALSE,
                           simplifyVector = TRUE,
+                          domain = url,
                           ...) {
 
   #Step 1: Identify the return type to be sent back..
@@ -54,7 +56,8 @@ flywire_fetch <- function(url,
   if (is.null(config))
     config = httr::config()
   if(is.null(token))
-    token = chunkedgraph_token()
+    domain <- sub("^middleauth\\+", "", domain)
+    token <- chunkedgraph_token(url=domain)
   if(!isTRUE(is.na(token)))
     config = c(config, add_headers(Authorization = paste("Bearer", token)))
 

R/flywire-skeleton.R

@@ -116,6 +116,9 @@
 #'   fellow leaf nodes for the rerooting process. Will be inaccurate at values
 #'   that are too high or too low. Should be about the size of the expected
 #'   soma.
+#' @param reroot_method whether to try to reroot the neuron based on
+#'  mixed direction of vectors in the neuron at nearby point('direction') or
+#'  proximity of points alone ('density').
 #' @param x a \code{nat::neuron} object.
 #' @param brain a \code{mesh3d} or \code{hxsurf} object within which a soma
 #'   cannot occur. For the re-rooting process. (Insect somata tend to lie
@@ -158,6 +161,8 @@
 #' @param inv_dist numeric.For \code{method = "teasar"}. Distance along the mesh
 #'   used for invalidation of vertices. This controls how detailed (or noisy)
 #'   the skeleton will be.
+#' @param resample stepsize by which to resample a neuron once skeletonised, but
+#' before healing or root finding. If \code{NULL}, no resampling attempt is made.
 #' @param cpu double (of length one). Set a limit on the total cpu time in
 #'   seconds.
 #' @param elapsed double (of length one). Set a limit on the total elapsed cpu
@@ -281,6 +286,7 @@ skeletor <- function(segments = NULL,
                      reroot = TRUE,
                      k.soma.search = 10,
                      radius.soma.search = 2500,
+                     reroot_method = c("direction","density"),
                      brain = NULL,
                      n = 5,
                      n_rays = 20,
@@ -295,6 +301,7 @@ skeletor <- function(segments = NULL,
                      inv_dist = 100,
                      cpu = Inf,
                      elapsed = Inf,
+                     resample = NULL,
                     ...){
   if(is.null(segments)&&is.null(obj)){
     stop("Either the argument segments or obj must be given.")
@@ -317,9 +324,10 @@ skeletor <- function(segments = NULL,
   method = match.arg(method)
   projection = match.arg(projection)
   cluster_pos = match.arg(cluster_pos)
+  reroot_method = match.arg(reroot_method)
   segments = unique(segments)
   py_skel_imports()
-  py_cloudvolume(cloudvolume.url, ...)
+  py_cloudvolume(cloudvolume.url) # ...
   neurons = nat::neuronlist()
   pb <- progress::progress_bar$new(
     format = sprintf("  %s [:bar] :current/:total eta: :eta", msg),
@@ -351,6 +359,7 @@ skeletor <- function(segments = NULL,
                                          reroot = reroot,
                                          k.soma.search = k.soma.search,
                                          radius.soma.search = radius.soma.search,
+                                         reroot_method = reroot_method,
                                          brain = brain,
                                          n = n,
                                          n_rays = n_rays,
@@ -363,6 +372,7 @@ skeletor <- function(segments = NULL,
                                          shape_weight = shape_weight,
                                          sample_weight = sample_weight,
                                          inv_dist = inv_dist,
+                                         resample = resample,
                                          ...))),
                           cpu = cpu,
                           elapsed = elapsed)
@@ -440,6 +450,7 @@ py_skeletor <- function(id,
                         reroot = TRUE,
                         k.soma.search = 10,
                         radius.soma.search = 2500,
+                        reroot_method = c('direction','density'),
                         brain = NULL,
                         n = 5,
                         n_rays = 20,
@@ -452,13 +463,15 @@ py_skeletor <- function(id,
                         shape_weight = 1,
                         sample_weight = 0.1,
                         inv_dist = 100,
+                        resample = NULL,
                         ...){
   stopifnot(length(id)==1)
   operator = match.arg(operator)
   method.radii = match.arg(method.radii)
   method = match.arg(method)
   projection = match.arg(projection)
   cluster_pos = match.arg(cluster_pos)
+  reroot_method = match.arg(reroot_method)
   if(grepl("obj$",id)){
     obj.file = TRUE
     reticulate::py_run_string(sprintf("m=tm.load_mesh('%s',process=False)",id), ...)
@@ -473,7 +486,7 @@ py_skeletor <- function(id,
     }))){
       py_cloudvolume(cloudvolume.url=cloudvolume.url)
     }
-    reticulate::py_run_string(sprintf("id=%s",id), ...)
+    reticulate::py_run_string(sprintf("id=%s",id)) # ..
     # this can error out with a bad connection to the server
     counter = 3 # we'll try 3 times
     while(counter>0){
@@ -539,13 +552,24 @@ py_skeletor <- function(id,
   swc = skel$swc
   colnames(swc) = c("PointNo","Parent","X","Y","Z","W")
   neuron = nat::as.neuron(swc)
-  if(heal){
-    neuron = suppressMessages(nat::stitch_neurons_mst(x = neuron, threshold = heal.threshold, k = heal.k))
-  }else{
-    neuron = subtree(neuron)
+  if(!is.null(resample)){
+    neuron <- nat::resample(neuron, stepsize = resample)
+  }
+  if(neuron$nTree>1){
+    if(heal){
+      neuron = suppressMessages(nat::stitch_neurons_mst(x = neuron,
+                                                        threshold = heal.threshold,
+                                                        k = heal.k))
+    }else{
+      neuron = subtree(neuron)
+    }
   }
   if(reroot){
-    neuron = tryCatch(reroot_hairball(neuron, k.soma.search = k.soma.search, radius.soma.search = radius.soma.search, brain = brain),
+    neuron = tryCatch(reroot_hairball(neuron,
+                                      k.soma.search = k.soma.search,
+                                      radius.soma.search = radius.soma.search,
+                                      reroot_method = reroot_method,
+                                      brain = brain),
                       error = function(e){
                         warning(e)
                         neuron
@@ -583,13 +607,15 @@ py_skeletor <- function(id,
 reroot_hairball <- function(x,
                            k.soma.search = 10,
                            radius.soma.search = 2500,
-                           brain = NULL){
+                           brain = NULL,
+                           reroot_method = c("direction","density")){
 
   # Get end and branch points, as vectors
+  reroot_method = match.arg(reroot_method)
   e = nat::endpoints(x)
   if(!is.null(brain)){
     pin = !nat::pointsinside(x = x$d, surf = brain)
-    pin[is.na(pin)||is.infinite(pin)||is.nan(pin)] = FALSE
+    pin[is.na(pin)|is.infinite(pin)|is.nan(pin)] = FALSE
     ins = c(1:nrow(x$d))[pin]
     ee = intersect(e, ins)
     if(length(ee)>=1){
@@ -614,9 +640,15 @@ reroot_hairball <- function(x,
     idx = idx[apply(idx,1,function(r) sum(r>0)>3),]
 
     # Asses vector direction
-    l = lapply(rownames(idx), function(r) sum(abs(apply(v[idx[r,],],1,function(vr) crossprod3D(vr, v[r,],i=3) ) )))
-    u = unlist(l)
-    root = rownames(idx)[which.max(u)]
+    if(reroot_method=='direction'){
+      l = lapply(rownames(idx), function(r) sum(abs(apply(v[idx[r,],],1,function(vr) crossprod3D(vr, v[r,],i=3) ) )))
+      u = unlist(l)
+      root = rownames(idx)[which.max(u)]
+    }else{
+      # Point with the most near points
+      a=apply(dists,1,sum)
+      root=names(which.min(a))
+    }
   }
 
   # Re-root neuron