R/hemibrain_synapses.R

In natverse/hemibrainr: Code for Working with Data from Janelia FlyEM's Hemibrain Project and the flywire data sets``

############################################################################
################################ Synapses ##################################
############################################################################

#' Extract synapse location and labels, and edgelists from a neuron/neuronlist
#'
#' @description Extract a single data frame describing synapse/connection types,
#'   partners, locations and position on a neuron's axon/dendrite. You can
#'   either get all synapses returned or all unitary connections to a neuron's
#'   partners returned. Broken down by axon/dendrite (\code{label}), and
#'   pre/postsynapses or pre/postsynaptic partners. Note that \code{hemibrain_extract_compartment_edgelist} will
#'   only return connections between neurons given in the argument \code{x}.
#'
#' @inheritParams flow_centrality
#'
#' @param prepost whether to get presynapses, postsynapses or both
#' @param meta a data frame giving neuron metrics (i.e. as generated by \code{hemibrain_compartment_metrics}) and transmitter predictions (column \code{top_nt}) for the neurons supplied to these functions.
#'
#' @return a \code{data.frame}. Depending on which synapse function was called, it can contain the columns:
#'
#' \itemize{
#'
#'   \item{"treenode_id"} { - the position of the node in the SWC-style table found at \code{neuron$d}, where the neuron is the skeleton for \code{bodyid}.}
#'
#'   \item{"connector_id"}{ - the unique ID for a pre/post synapse, as read from neuPrint. If this is not given, you are looking at a connection not a synapse.
#'   In this case \code{count} should be given, which shows the number of synapses in this connection.}
#'
#'   \item{"prepost"}{ - whether the given synapse is a pre-synapse (0, output synapse) or postsynapse (1, input synapse). Alternatively, if a connection is given,
#'   whether this connection is presynaptic to \code{bodyid} (0, \code{bodyid} is target) or postsynaptic (1, \code{bodyid} is source).}
#'
#'   \item{"x"}{ - x coordinate for the root point.}
#'
#'   \item{"y"}{ - y coordinate for the root point.}
#'
#'   \item{"z"}{ - z coordinate for the root point.}
#'
#'   \item{"confidence"}{ - FlyEM's confidence level. The lower the score, the more likely this synapse is an artefact.}
#'
#'   \item{"bodyid"}{ - The neuPrint neuron/body related to the synapse/connection given in each row.}
#'
#'   \item{"partner"}{ - The neuron connecting to \code{bodyid} by the given synapse/connection.}
#'
#'   \item{"pre"}{ - The body ID for the presynaptic (source) neuron.}
#'
#'   \item{"partner"}{ - The body ID for the presynaptic (target) neuron.}
#'
#'   \item{"label"}{ - The compartment of the \code{bodyid} neuron on which the synapse is placed / which receives/makes the given connection.
#'   See \code{?standardise}.}
#'
#'   \item{"partner_label"}{ - The compartment of the \code{partner} neuron on which the synapse is placed / which receives/makes the given connection.}
#'
#'   \item{"count"}{ - The number of synapses that make the given connection. Sometimes referred to as 'weight'.}
#'
#'   \item{"norm"}{ - The normalised synapse weight. \code{count} is divided by the total number of inputs that the
#'   target neuron's (\code{post}) compartment (\code{label}) has. I.e. this normalisation is by total_inputs onto a dendrite or axon, not the whole neuron.}
#'
#'   \item{"connection"}{ - The type of compartment-compartment connection specified by this row. The first compartment is the source (pre), the second, the target (post).}
#'
#'}
#'
#' @export
#' @rdname hemibrain_extract_connections
#' @seealso \code{\link{flow_centrality}}
#' @importFrom nat nlapply is.neuron is.neuronlist
#' @examples
#' \dontrun{
#' # Choose bodyids
#' some.pns  = sample(pn.ids, 20)
#'
#' # Read in these neurons
#' neurons = neuprintr::neuprint_read_neurons(some.pns)
#'
#' # Re-root
#' neurons.flow = flow_centrality(neurons,
#' polypre = TRUE,
#' mode = "centrifugal", split = "distance")
#'
#' # Let's check that this worked
#' syns = hemibrain_extract_synapses(neurons.flow)
#'
#' # Get the edgelist by compartment
#' elist = hemibrain_extract_compartment_edgelist(neurons.flow)
#'
#' # See result
#' library(nat)
#' plot3d_split(neuron.flow)
#' points3d(xyzmatrix(subset(syns,prepost==1)), col = "cyan")
#' points3d(xyzmatrix(subset(syns,prepost==0)), col = "red")
#' }
hemibrain_extract_synapses <- function(x,
                                       prepost = c("BOTH","PRE","POST"),
                                       ...){
  prepost = match.arg(prepost)
  if(nat::is.neuronlist(x)){
    if("bodyid"%in%colnames(x[,])){
      x = add_field_seq(x,x[,"bodyid"],field="bodyid")
    }else if("root_id"%in%colnames(x[,])){
      x = add_field_seq(x,x[,"root_id"],field="root_id")
    }else if("skid"%in%colnames(x[,])){
      x = add_field_seq(x,x[,"skid"],field="skid")
    }else if(!is.null(names(x))){
      x = add_field_seq(x,names(x),field="id")
    }
    syns.list = nat::nlapply(x,extract_synapses, unitary = FALSE, ...)
    syns = do.call(plyr::rbind.fill,syns.list)
  }else if (nat::is.neuron(x)){
    syns = extract_synapses(x, unitary = FALSE)
  }else if(is.data.frame(x) && "prepost"%in%colnames(x)){
    syns = extract_synapses(x, unitary = FALSE)
  }else{
    stop("x must be a neuron or neuronlist object")
  }
  if(prepost=="PRE"){
    syns = syns[syns$prepost==0,]
  }else if (prepost=="POST"){
    syns = syns[syns$prepost==1,]
  }
  if(!is.null(syns)){
    rownames(syns) = 1:nrow(syns)
  }
  syns
}

#' @export
#' @rdname hemibrain_extract_connections
hemibrain_extract_connections <- function(x,
                                       prepost = c("BOTH","PRE","POST"),
                                       meta = NULL,
                                       ...){
  prepost = match.arg(prepost)
  if(nat::is.neuronlist(x)){
    if("bodyid"%in%colnames(x[,])){
      x = add_field_seq(x,x[,"bodyid"],field="bodyid")
    }else if("root_id"%in%colnames(x[,])){
      x = add_field_seq(x,x[,"root_id"],field="root_id")
    }else if("skid"%in%colnames(x[,])){
      x = add_field_seq(x,x[,"skid"],field="skid")
    }else if(!is.null(names(x))){
      x = add_field_seq(x,names(x),field="id")
    }
    syns = nat::nlapply(x, extract_synapses, unitary = TRUE, meta = meta, ...)
    syns = do.call(rbind,syns)
  }else if (nat::is.neuron(x)){
    syns = extract_synapses(x, unitary = TRUE)
  }else if(is.data.frame(x) && "prepost"%in%colnames(x)){
    syns = extract_synapses(x, unitary = TRUE)
  }else{
    stop("x must be a neuron or neuronlist object")
  }
  if(prepost=="PRE"){
    syns = syns[syns$prepost==0,]
  }else if (prepost=="POST"){
    syns = syns[syns$prepost==1,]
  }
  if(length(syns)){
    syns$label = standard_compartments(syns$label)
    rownames(syns) = 1:nrow(syns)
    syns = syns[order(syns$count, decreasing = TRUE),]
  }
  syns
}

#' @importFrom magrittr %>%
#' @export
magrittr::`%>%`
# hidden
#' @importFrom dplyr filter mutate group_by distinct select n case_when
#' @importFrom rlang .data
extract_synapses <-function(x, pre_id = "pre_id", unitary = FALSE, meta = NULL){ # top_nt may be taken on a per synapse, rather than per neuron, basis
  if(nat::is.neuron(x)){
    syn = x$connectors
  }else{
    syn = x
  }
  if(!nrow(syn)){
    warning("neuron ", x$id," has no synapses")
    return(NULL)
  }
  if(!is.null(x$root_id)){
    id = "root_id"
  }else if(!is.null(x$skid)){
    id = "skid"
  } else if(!is.null(x$bodyid)){
      id = "bodyid"
  }else{
    id = "id"
  }
  syn[[id]] = nullToNA(as.character(x[[id]]))
  syn[[id]] = gsub(" ","",syn[[id]])
  colnames(syn) = snakecase::to_snake_case(colnames(syn))
  if(is.null(syn$label)){
    syn$label = nullToNA(x$d$Label[match(syn$treenode_id,x$d$PointNo)])
  }
  poss.nts=c("gaba", "acetylcholine", "glutamate", "octopamine", "serotonin", "dopamine", "neither")
  if(!all(c(poss.nts,"top_nt")%in%colnames(syn))){
    for(pnt in setdiff(poss.nts,colnames(syn))){
      syn[[pnt]] = NA
    }
  }
  poss.nts = intersect(poss.nts, colnames(syn))
  if(!is.null(syn$top_nt)){ # top_nt on a per synapse level
    syn$syn_top_nt = syn$top_nt
    syn$top_nt = NULL
    syn$syn_top_p = syn$top_p
    syn$top_p = NULL
  }else{
    syn$syn_top_nt = NA
    syn$syn_top_p = NA
  }
  # Now work out top_nt at the neuron level
  if(!is.null(meta)){
    if(all(c(id,"top_nt")%in%colnames(meta))){
      meta[[pre_id]]  = meta[[id]]
      if(bit64::is.integer64(syn[[pre_id]])){
        meta[[pre_id]] = bit64::as.integer64(meta[[pre_id]])
      }else if(is.character(syn[[pre_id]])){
        meta[[pre_id]] = as.character(meta[[pre_id]])
      }else if(is.numeric(syn[[pre_id]])){
        syn[[pre_id]] = as.character(syn[[pre_id]])
      }
      syn = dplyr::inner_join(syn, meta[,c(pre_id,"top_nt")],
                              by = pre_id,
                              copy = TRUE,
                              auto_index = TRUE)
    }else{
      warning("top_nt must be in colnames(meta)")
    }
  }else if(length(poss.nts) && is.null(syn$top_nt)){
    bad = tryCatch(sum(syn[,colnames(syn)%in%poss.nts],na.rm = TRUE)==0, error = function(e) 1)
    if(bad){
      syn$top_nt = "unknown"
      syn$top_p = 0
    }else{
      if(!is.null(x$ntpred)){
        topnt = top_nt(x, id = id)
      }else if (id == "bodyid"){
        topnt = get_top_nt(syn)
      }
      syn$top_nt = topnt$top_nt
      syn$top_p = topnt$top_p
    }
  }else{
    syn$top_nt = "unknown"
    syn$top_p = 0
  }
  if(unitary){ # connections, rather than synapses
    if(id == "root_id"){
      syn %>%
        dplyr::mutate(partner = dplyr::case_when(
          .data$prepost==0 ~ .data$post_id,
          .data$prepost==1  ~ .data[[pre_id]]
        )) %>%
        dplyr::mutate(prepost = dplyr::case_when(
          .data$prepost==0 ~ 1,
          .data$prepost==1  ~ 0
        )) %>% # i.e. switch perspective, presynapses connect to postsynaptic partners
        dplyr::group_by(.data[[id]], .data$partner, .data$prepost, .data$label) %>%
        dplyr::mutate(count = dplyr::n()) %>%
        dplyr::distinct(.data[[id]], .data$partner, .data$prepost, .data$label, .data$count, .data$top_nt, .keep_all = FALSE) %>%
        as.data.frame(stringsAsFactors = FALSE) ->
        syn
    }else{
      syn %>%
        dplyr::mutate(prepost = dplyr::case_when(
          .data$prepost==0 ~ 1,
          .data$prepost==1  ~ 0
        )) %>% # i.e. switch perspective, presynapses connect to postsynaptic partners
        dplyr::group_by(.data[[id]], .data$partner, .data$prepost, .data$label) %>%
        dplyr::mutate(count = dplyr::n()) %>%
        dplyr::distinct(.data[[id]], .data$partner, .data$prepost, .data$label, .data$count, .data$top_nt, .keep_all = FALSE) %>%
        as.data.frame(stringsAsFactors = FALSE) ->
        syn
    }
  }
  syn$label = standard_compartments(syn$label)
  syn
}

#' @export
#' @rdname hemibrain_extract_connections
hemibrain_extract_compartment_edgelist <- function(x, meta = NULL, ...){
  if(nat::is.neuronlist(x)){
    y = x[[1]]
    if(!is.null(y$root_id)){
      id = "root_id"
      partner = "pre_id" #"post_id" ????
    } else if(!is.null(y$skid)){
      id = "skid"
      partner = "partner"
    } else if(!is.null(y$bodyid)){
      id = "bodyid"
      partner = "partner"
    }else{
      id = "id"
      partner = "partner"
    }
    x = add_field_seq(x,names(x),field=id)
    syns.list = nat::nlapply(x, extract_synapses, pre_id = partner, unitary = FALSE, meta = meta, ...)
  }else{
    syns.list = x
    warning("x should be a neuronlist object")
  }
  if(!is.null(meta)){
    lookup.nt = meta[,c(id,"top_nt","top_p")]
    lookup.nt = as.data.frame(lookup.nt)
    rownames(lookup.nt) = lookup.nt[[id]]
  }else{
    lookup.nt = nat::nlapply(syns.list, get_top_nt)
    lookup.nt = lapply(names(lookup.nt), function(n){
      v = lookup.nt[[n]]
      rownames(v) = n
      v[[id]] = n
      v
    })
    lookup.nt = do.call(plyr::rbind.fill,lookup.nt)
    rownames(lookup.nt) = lookup.nt[[id]]
  }
  names(syns.list) = NULL
  lookup = nat::nlapply(syns.list, extract_lookup,...)
  lookup = unlist(lookup)
  elists = nat::nlapply(syns.list, extract_elist, lookup = lookup, lookup.nt = lookup.nt, id = id, partner = partner, ...)
  elist = do.call(rbind, elists)
  if(length(elist)){
    rownames(elist) = 1:nrow(elist)
    elist = elist[order(elist$norm, decreasing = TRUE),]
    elist = elist[order(elist$count, decreasing = TRUE),]
  }
  elist
}

# hidden, for one pre neuron
extract_elist <- function(syns, lookup, lookup.nt = NULL, id = "bodyid", partner = "partner"){
  colnames(syns) = snakecase::to_snake_case(colnames(syns))
  ids = unique(nullToNA(as.character(syns[[id]])))
  compartment.inputs = c()
  for(l in unique(syns$label)){
    d.post = nrow(subset(syns, syns$prepost==1 & syns$label==l))
    names(d.post) = l
    compartment.inputs = c(compartment.inputs, d.post)
  }
  syns %>%
    # Re-name for clarity
    dplyr::filter(.data$prepost==1) %>%
    dplyr::rename(post = .data[[id]]) %>%
    dplyr::rename(pre = .data[[partner]]) %>%
    dplyr::rename(post_label = .data$label) %>%
    # Compartment labels
    dplyr::mutate(pre_label = lookup[as.character(.data$connector_id)]) %>%
    dplyr::mutate(pre_label = ifelse(is.na(.data$pre_label),"unknown",.data$pre_label)) %>%
    # Transmitter
    { if(is.null(lookup.nt)){
      dplyr::mutate(.,post_top_nt = "unknown") %>%
      dplyr::mutate(.,pre_top_nt = "unknown") %>%
      dplyr::mutate(.,post_top_p = 0) %>%
      dplyr::mutate(.,pre_top_p = 0)
    }else{
      dplyr::mutate(.,pre_top_nt = lookup.nt[as.character(.data$pre),"top_nt"]) %>%
      dplyr::mutate(.,post_top_nt = lookup.nt[as.character(.data$post),"top_nt"]) %>%
      dplyr::mutate(.,pre_top_p = lookup.nt[as.character(.data$pre),"top_p"]) %>%
      dplyr::mutate(.,post_top_p = lookup.nt[as.character(.data$post),"top_p"])
    } }  %>%
    # Synapse counts
    dplyr::group_by(.data$post, .data$pre, .data$post_label, .data$pre_label) %>%
    dplyr::mutate(count = dplyr::n()) %>%
    # Normalised synapses, by compartment
    dplyr::ungroup() %>%
    dplyr::group_by(.data$post,.data$post_label) %>%
    dplyr::mutate(norm = .data$count/compartment.inputs[as.character(.data$post_label)]) %>%
    # Clean up
    dplyr::distinct(.data$post, .data$pre,.data$post_label, .data$pre_label, .data$count, .data$norm,
                    .data$pre_top_nt, .data$pre_top_p, .data$post_top_nt, .data$post_top_p, .keep_all = FALSE) %>%
    dplyr::filter(!is.na(.data$pre_label) & .data$count > 0) %>%
    as.data.frame(stringsAsFactors = FALSE) ->
    elist
  elist$post_top_nt[is.na(elist$post_top_nt)] = "unknown"
  elist$pre_top_nt[is.na(elist$pre_top_nt)] = "unknown"
  elist$post_top_p[is.na(elist$post_top_p)] = 0
  elist$pre_top_p[is.na(elist$pre_top_p)] = 0
  if(nrow(elist)){
    rownames(elist) = 1:nrow(elist)
    elist$post_label = standard_compartments(elist$post_label)
    elist$pre_label = standard_compartments(elist$pre_label)
    elist$connection = paste(elist$pre_label,elist$post_label,sep="-")
    elist
  }else{
    NULL
  }
}

# hidden
extract_lookup <- function(syns, type = c("label","top_nt")){
  type = match.arg(type)
  colnames(syns) = snakecase::to_snake_case(colnames(syns))
  syns %>%
    dplyr::filter(.data$prepost==0) %>%
    dplyr::distinct(.data$connector_id, .data[[type]], .keep_all = FALSE) %>%
    as.data.frame(stringsAsFactors = FALSE) ->
    conn.lookup
  lookup = conn.lookup[[type]]
  names(lookup) = as.character(conn.lookup$connector_id)
  lookup = lookup[!names(lookup)%in%c("0","NA")]
  lookup
}

# read synapes
hemibrain_syns_nt <- function(nt.file = "/Volumes/GoogleDrive/Shared drives/hemibrain/fafbsynapses/synister_hemi_whole_volume_v0_t8.feather"){
  if(grepl("parquet$",nt.file)){
    arf = arrow::read_parquet(nt.file)
    colnames(arf) = gsub("nts_8\\.","",colnames(arf))
    arf %>%
      dplyr::rename(connector_id = id)
  }else{
    arf = arrow::read_feather(nt.file)
    colnames(arf) = gsub("nts_3\\.","",colnames(arf))
    arf %>%
      dplyr::rename(connector_id = id)
  }
}

# read synapes
flywire_syns_nt <- function(nt.file = "/Volumes/GoogleDrive/Shared drives/hemibrain/fafbsynapses/synister_fafb_whole_volume_v3_t11.parquet"){
  if(grepl("feather$",nt.file)){
    arf = arrow::read_feather(nt.file)
    colnames(arf) = gsub("nts_3\\.","",colnames(arf))
    arf %>%
      dplyr::rename(connector_id = id)
  }else{
    arf = arrow::read_parquet(nt.file)
    colnames(arf) = gsub("nts_8\\.","",colnames(arf))
    arf %>%
      dplyr::rename(connector_id = id)
  }
}

# hidden
hemibrain_add_nt <- function(x,
                   nts = hemibrain_syns_nt(nt.file),
                   nt.file = "/Volumes/GoogleDrive/Shared drives/hemibrain/fafbsynapses/synister_hemi_whole_volume_v0_t8.feather",
                   classic = FALSE,
                   ...) UseMethod("hemibrain_add_nt")

# hidden
hemibrain_add_nt.neuronlist <- function(x, nts = NULL, nt.file = "/Volumes/GoogleDrive/Shared drives/hemibrain/fafbsynapses/synister_hemi_whole_volume_v0_t8.feather", classic = FALSE, ...){
  if(is.null(nts)){
    nts = hemibrain_syns_nt(nt.file = nt.file)
  }
  x = nat::nlapply(x, hemibrain_add_nt.neuron, nts = nts, nt.file=nt.file, classic=classic)
  df.nt = do.call(plyr::rbind.fill, lapply(x, function(y) y$ntpred))
  df.nt[is.na(df.nt)] = 0
  x[,] = cbind(x[,setdiff(colnames(x[,]),colnames(df.nt))],df.nt)
  x
}

# hidden
hemibrain_add_nt.neuron <- function(x,
                                    nts = NULL,
                                    nt.file = "/Volumes/GoogleDrive/Shared drives/hemibrain/fafbsynapses/synister_hemi_whole_volume_v0_t8.feather",
                                    classic = FALSE,
                                    ...){
  if(is.null(nts)){
    nts = hemibrain_syns_nt(nt.file=nt.file)
  }
  syns = x$connectors
  poss.nts = c("gaba", "acetylcholine", "glutamate", "serotonin", "octopamine", "dopamine", "neither")
  syns[,colnames(syns)%in%poss.nts] = NULL
  syns.nt = syns %>%
    dplyr::left_join(nts[,intersect(colnames(nts),c("x", "y", "z", poss.nts))], by = c("x","y","z"))
  tops = syns.nt[,c(poss.nts)]
  tops[,'top_p']=do.call(pmax, as.list(tops[poss.nts]))
  top.col=max.col(tops[poss.nts], ties.method = "first")
  tops[,'top_nt']=poss.nts[top.col]
  syns.nt = cbind(syns.nt,tops[,c(poss.nts,"top_p","top_nt")])
  x$connectors = syns.nt
  x$ntpred = get_top_nt(syns.nt, classic=classic, ...)
  x
}

# find top_nt
get_top_nt <- function(syns.nt,
                       poss.nts=c("gaba", "acetylcholine", "glutamate", "octopamine", "serotonin","dopamine", "neither"),
                       classic = FALSE,
                       confidence.thresh = 0.5,
                       cleft.threshold = 75){
  if("confidence.thresh"%in%colnames(syns.nt)){
    syns.nt = subset(syns.nt, syns.nt$confidence >= confidence.thresh)
    #   syns.nt = syns.nt[!duplicated(syns.nt$confidence),]
  }
  if("cleft.threshold"%in%colnames(syns.nt)){
    syns.nt = subset(syns.nt, syns.nt$confidence >= cleft.threshold)
  }
  if(is.null(syns.nt)){
    return( data.frame(top_nt = "unknown", top_p = "unknown") )
  }else if(!length(syns.nt)){
    return( data.frame(top_nt = "unknown", top_p = "unknown") )
  }else if(!nrow(syns.nt)){
    return( data.frame(top_nt = "unknown", top_p = "unknown") )
  }
  colnames(syns.nt) = snakecase::to_snake_case(colnames(syns.nt))
  if("prepost" %in% colnames(syns.nt)){
    syns.nt = subset(syns.nt, syns.nt$prepost == 0)
  }
  if("label" %in% colnames(syns.nt)){
    syns.nt$label = hemibrainr:::standard_compartments(syns.nt$label)
    syns.nt = subset(syns.nt, `label` %in% c("axon","dendrite"))
  }
  if(classic){
    poss.nts = c("gaba", "acetylcholine", "glutamate", "neither")
  }
  nts = syns.nt[,c(poss.nts)]
  nts[,'top_p']=do.call(pmax, as.list(nts[poss.nts]))
  top.col=max.col(nts[poss.nts], ties.method = "first")
  nts[,'top_nt']=poss.nts[top.col]
  tx=table(nts$top_nt)
  tx=tx/sum(tx)*100
  tops = colSums(nts[,poss.nts])/nrow(nts)
  top_p = max(tops)
  top_nt = names(which.max(tx))
  trans = t(c(tx))
  trans[is.na(trans)] = 0
  trans = as.data.frame(trans)
  trans$top_nt = ifelse(is.null(top_nt),"unknown",top_nt)
  trans$top_p = ifelse(is.null(top_nt),"unknown",top_p)
  trans %>%
    dplyr::mutate_if(is.numeric, round, 3)
}

# hidden
hemibrain_ntplot <- function(x,
                             poss.nts=c("gaba", "acetylcholine", "glutamate","octopamine", "serotonin", "dopamine", 'neither'),
                             classic = FALSE,
                             confidence.thresh = 0.5,
                             collapse = FALSE,
                             y.plot = c("confidence", "top_p"),
                             ...) UseMethod("hemibrain_ntplot")

# plot neurotransmitter distributions
hemibrain_ntplot.neuron <- function(x,
                                    poss.nts=c("gaba", "acetylcholine", "glutamate","octopamine", "serotonin", "dopamine", 'neither'),
                                    classic = FALSE,
                                    confidence.thresh = 0.5,
                                    collapse = FALSE,
                                    y.plot = c("confidence", "top_p"),
                                    ...) {
  check_package_available('ggplot2')
  poss.nts=match.arg(poss.nts, several.ok = T)
  syns = x$connectors
  if("prepost" %in% colnames(syns)){
    syns = subset(syns, syns$prepost == 0)
  }
  syns.nt=dplyr::filter(syns, .data$confidence>=confidence.thresh &
                    .data$top_nt %in% poss.nts)
  syns.nt$id = x$bodyid
  hemibrain_ntplot.data.frame(x=syns.nt, syns.nt=syns.nt, poss.nts=poss.nts, classic=classic, confidence.thresh=confidence.thresh,y.plot=y.plot, ...)
}

# Plot from a neuronlist
hemibrain_ntplot.neuronlist <- function(x,
                                    poss.nts=c("gaba", "acetylcholine", "glutamate","octopamine", "serotonin", "dopamine", 'neither'),
                                    classic = FALSE,
                                    confidence.thresh = 0.5,
                                    collapse = FALSE,
                                    y.plot = c("confidence", "top_p"),
                                    ...) {
  check_package_available('gridExtra')
  if(!collapse){
    plist = lapply(x, hemibrain_ntplot.neuron, poss.nts=poss.nts, classic=classic, confidence.thresh=confidence.thresh, y.plot=y.plot,...)
    n <- length(plist)
    nCol <- floor(sqrt(n))
    do.call(gridExtra::grid.arrange, c(plist, ncol=nCol))
  }else{
    syns.nt = lapply(x, function(y) y$connectors)
    syns.nt = do.call(plyr::rbind.fill, syns.nt)
    syns.nt$id = "collapse"
    hemibrain_ntplot.data.frame(syns.nt)
  }
}

# hidden
hemibrain_ntplot.data.frame <- function(x,
                                        poss.nts=c("gaba", "acetylcholine", "glutamate","octopamine", "serotonin", "dopamine", 'neither'),
                                        classic = FALSE,
                                        confidence.thresh = 0.5,
                                        collapse = FALSE,
                                        y.plot = c("confidence", "top_p"),
                                        ...) {
  check_package_available('ggplot2')
  poss.nts=match.arg(poss.nts, several.ok = TRUE)
  y.plot=match.arg(y.plot)
  plot.p = y.plot=="top_p"
  if("id" %in% colnames(x)){
    tit = x$id[1]
  }else{
    tit = ""
  }
  if("prepost" %in% colnames(x)){
    x = subset(x, `prepost` == 0)
  }
  if(classic){
    poss.nts = c("gaba", "acetylcholine", "glutamate",'neither')
  }
  syns.nt=dplyr::filter(x, .data$confidence>=confidence.thresh & .data$top_nt %in% poss.nts)
  colnames(syns.nt) = snakecase::to_snake_case(syns.nt)
  ntcols = c(
    gaba = "#E6A749",
    acetylcholine = "#4B506B",
    glutamate = "#70B657",
    octopamine = "#7A4F98",
    serotonin = "#93A3CF",
    dopamine = "#CF6F6C",
    neither = "grey70")[poss.nts]
  if("label" %in% colnames(syns.nt) & !collapse){
    syns.nt$label = hemibrainr:::standard_compartments(syns.nt$label)
    syns.nt = subset(syns.nt, `label` %in% c("axon","dendrite"))
    ggplot2::ggplot(data=syns.nt, ggplot2::aes(x = .data[[y.plot]], fill = `top_nt`)) +
      ggplot2::geom_histogram(binwidth = 0.01) +
      ggplot2::scale_fill_manual(values = ntcols) +
      ggplot2::theme_minimal() +
      ggplot2::theme(legend.position="none") +
      ggplot2::labs(x = tit) +
      ggplot2::facet_wrap(~`label`, ncol=1)
  }else{
    ggplot2::ggplot(data=syns.nt, ggplot2::aes(x= .data[[y.plot]], fill = `top_nt`)) +
      ggplot2::geom_histogram(binwidth = 0.01) +
      ggplot2::scale_fill_manual(values = ntcols) +
      ggplot2::theme_minimal() +
      ggplot2::theme(legend.position="none") +
      ggplot2::labs(x = tit)
  }
}

# Search
hemibrain_ntplot.character <- function(x,
                                       poss.nts=c("gaba", "acetylcholine", "glutamate","octopamine", "serotonin", "dopamine", 'neither'),
                                       classic = FALSE,
                                       confidence.thresh = 0.5,
                                       collapse = FALSE,
                                       y.plot = c("confidence", "top_p"),
                                       ...) {
  n = neuprint_search(x, ...)
  bids = hemibrain_neuron_bodyids()
  bids = intersect(n$bodyid, bids)
  neurons = hemibrain_read_neurons(bids, ...)
  neurons.nt = hemibrain_add_nt.neuronlist(neurons, poss.nts=poss.nts, classic=classic, confidence.thresh=confidence.thresh)
  hemibrain_ntplot.neuronlist(neurons.nt, poss.nts=poss.nts, classic=classic, confidence.thresh=confidence.thresh, collapse=collapse, y.plot=y.plot)
  ntcols = c(
    gaba = "#E6A749",
    acetylcholine = "#4B506B",
    glutamate = "#70B657",
    octopamine = "#7A4F98",
    serotonin = "#93A3CF",
    dopamine = "#CF6F6C",
    neither = "grey70")
  plot3d(neurons.nt, lwd = 2, col = ntcols[neurons.nt[,"top_nt"]])
  print(neurons.nt[,])
  return(neurons.nt[,])
}

# Get top_nt
top_nt <- function(x,  id = "root_id"){
  y = x$ntpred
  y = y[sapply(y,is.numeric)]
  y =  y[!names(y)%in%c("top_p","top_nt")]
  top_nt = names(y)[which.max(y)]
  top_p = max(y)
  if(is.null(x[[id]])){
    df = data.frame(top_nt=top_nt,top_p=top_p)
  }else{
    df = data.frame(id=x[[id]], top_nt=top_nt,top_p=top_p)
    colnames(df) = c(id, "top_nt", "top_p")
  }
  df
}