R/flywire_synapses.R

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

#' @rdname flywire_googledrive_data
#' @export
flywire_elist <- function(local = FALSE, folder = "flywire_neurons/", sql = TRUE, ...){
  savedir = good_savedir(local = local)
  if(sql){
    gcsv = find_gsql(savedir = savedir, tab = "flywire_edgelist", sql.db = "flywire_data.sqlite", folder = folder, ...)
  }else{
    gfile = find_gfile(savedir = savedir, file = "flywire_edgelist", folder = folder)
    gcsv = suppressWarnings(readr::read_csv(gfile, col_types = sql_col_types))
  }
  gcsv
}

#' @rdname flywire_googledrive_data
#' @export
flywire_connections <- function(local = FALSE, folder = "flywire_neurons/", sql = TRUE, ...){
  savedir = good_savedir(local = local)
  if(sql){
    gcsv = find_gsql(savedir = savedir, tab = "flywire_connections", sql.db = "flywire_data.sqlite", folder = folder, ...)
  }else{
    gfile = find_gfile(savedir = savedir, file = "flywire_connections", folder = folder)
    gcsv = suppressWarnings(readr::read_csv(gfile, col_types = sql_col_types))
  }
}

#' @rdname flywire_googledrive_data
#' @export
flywire_synapses <- function(local = FALSE, folder = "flywire_neurons/", simplified = FALSE, sql = TRUE, ...){
  savedir = good_savedir(local = local)
  if(sql){
    if(simplified){
      gcsv = find_gsql(savedir = savedir, tab = "flywire_simplified_synapses", sql.db = "flywire_data.sqlite", folder = folder, ...)
    }else{
      gcsv = find_gsql(savedir = savedir, tab = "flywire_synapses", sql.db = "flywire_data.sqlite", folder = folder, ...)
    }
  }else{
    if(simplified){
      gfile = find_gfile(savedir = savedir, file = "flywire_simplified_synapses", folder = folder)
    }else{
      gfile = find_gfile(savedir = savedir, file = "flywire_synapses", folder = folder)
    }
    gcsv = suppressWarnings(readr::read_csv(gfile, col_types = sql_col_types))
  }
  gcsv
}

#' Collapse Buhmann predicted flywire synapses into smaller units
#'
#' @description Many predicted Buhmann synapses may be located at one real synapse. Collapse Buhmann predicted flywire synapses into smaller units using a hierarchical clustering
#' based on the Euclidean distances between predicted 'presynapses'.
#'
#' @param x a  \code{data.frame} of synapses for FlyWire neurons, as generated by, for example, \code{fafbseg::flywire_partners}.
#' @param method when synapses are collapsed, they are either represented by the predicted synapse amongst them with the best score, cleft_score or a mean of
#' the meat data for each of these predicted synapses.
#' @param collapse logical, whether or not to collapse synapses in the same cluster to a single synapse. If \code{FALSE} then cluster groupings are given in the returned \code{data.frame}.
#' @param fast logical, whether or not to use \code{fastcluster::hclust} rather than \code{stats::hclust} to cluster synapses.
#'
#' @return a \code{data.frame} of collapsed synapses.
#'
#'
#' @examples
#' \donttest{
#' \dontrun{
#' syns=fafbseg::flywire_partners("720575940621039145")
#' syns.simp = flywire_synapse_simplify(syns)
#' }}
#' @seealso \code{\link{flywire_googledrive_data}}
#' @export
#' @importFrom stats dist
flywire_synapse_simplify <- function(x, method = c("cleft_scores","scores","mean"), collapse = TRUE, fast = TRUE){
  method = match.arg(method)
  h_clust = if(fast) {
    check_package_available('fastcluster')
    fastcluster::hclust
  } else stats::hclust

  # for each neuron
  if("pre_id"%in%colnames(x)){
    if(length(unique(x$pre_id))>1){
      pb <- progress::progress_bar$new(
        format = "  collapsing synapses [:bar] :current/:total eta: :eta",
        total = length(unique(x$pre_id)), clear = FALSE, show_after = 1
      )
      res = list()
      for(preid in unique(x$pre_id)){
        pb$tick()
        z = x[x$pre_id==preid,]
        y = flywire_synapse_simplify(x=z, method=method)
        preid = as.character(preid)
        res[[preid]]=y
      }
      return(do.call(plyr::rbind.fill,res))
    }
  }

  # Get 3D points
  x = x[order(x$scores,decreasing=TRUE),]
  x = x[order(x$cleft_scores,decreasing=TRUE),]
  if("prepost"%in%colnames(x)){
    pre = x[x$prepost==0,]
    post = x[x$prepost==1,]
  }else{
    pre = x
    post = data.frame()
  }
  xyz = tryCatch(nat::xyzmatrix(pre), error = function(e){
    pre[,c("x","y","z")] = pre[,c("pre_x","pre_y","pre_z")]
    nat::xyzmatrix(pre)
  })

  # Filter synapses
  near = try(nabor::knn(query = xyz, data = xyz, k = min(50,nrow(xyz)), radius = 5000),silent = TRUE)
  if('try-error'%in%class(near)){
    syns = plyr::rbind.fill(pre,post)
    syns$cluster = NA
    return(syns)
  }
  near$nn.dists[is.infinite(near$nn.dists)] = 5000
  near$nn.dists[is.na(near$nn.dists)] = 5000
  hc = try(h_clust(dist(near$nn.dists)), silent = TRUE)
  if(inherits(hc, "try-error")){
    ct = 1
    pre$cluster = 1
  }else{
    ct = dendextend::cutree(tree=hc, h = 2500)
    pre$cluster = ct
  }

  # For synapse in cluster, combine
  if(collapse){
    pre.new = data.frame()
    for(t in 1:max(ct)){
      np = pre[pre$cluster==t,]
      if(method=="mean"){
        nl = unlist(lapply(np,class))%in%c("numeric","integer")
        collap = as.data.frame(t(colMeans(np[,nl], na.rm = TRUE)))
        collap$offset = np$offset[1]
        collap$Label = np$Label[1]
        poss.nts=c("gaba", "acetylcholine", "glutamate", "octopamine", "serotonin","dopamine")
        tops = colSums(collap[,poss.nts])
        collap$top_p = max(tops)
        collap$top_nt = names(which.max(tops))
      }else{
        np = np[order(np[[method]],decreasing=TRUE),]
        collap = np[1,]
      }
      if(nrow(post)){
        post[post$pre_id%in%np$pre_id,"pre_id"] = collap$pre_id
        post[post$pre_svid%in%np$pre_svid,"pre_id"] = collap$pre_svid
      }
      pre.new = rbind(pre.new,collap)
    }
    post[!(duplicated(post$post_svid)&duplicated(post$pre_svid)),]
    plyr::rbind.fill(pre.new,post)
  }else{
    plyr::rbind.fill(pre,post)
  }
}