Analyze connectivity of FlyWire objects using natverse

[anhd...@gmail.com]

Dear Dr. Jefferis,

I'm curious if there is a way using your natverse to find connections (neuron-neuron synapse or glia-neuron contact) between 2 objects identified from FlyWire.

As you have mentioned in some threads, I can fetch meshes of FlyWire neurons into R using the fafbseg package. 

Once I read in 2 cells (either neuron-neuron or neuron-glia), is there a way for me to find/analyze where they make contacts?

Thank you so much, as always!

Regards,

Mike Le   

[Greg Jefferis]

Dear Mike,

Please find some ideas below for how to find mesh contacts. I have not had to answer exactly your question before, so this is just a quick example rather than a well tested optimal solution. The core operation (finding mesh to mesh distances) can be achieved in just a few lines of code. Deciding when one contact turns into two separate but adjacent contacts strikes me as the complicated part. Note that none of this has anything directly to do with synaptic connectivity.

Best wishes,

Greg.

Mesh to mesh distances

Note

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

Find the distance between meshes

There are two ways that I can immediately think of addressing this

1. Find the distance of every vertex on mesh 1 to its nearest neighbour on mesh 2
2. Find the (signed) distance from every vertex of mesh 1 to the surface of mesh 2

Obviously the ideal situation would be to find the mesh to mesh distances for the two objects (or even do a Euclidean distance transform on the voxel representation of the objects.) but neither of those are so accessible.

Option 1 is probably still a perfectly good option for most purposes since the meshes are quite densely sampled (i.e. the vertices are close together).

Setup and Mesh volume

Let’s use the same example that I mentioned on FlyWire slack.

Hide

library(fafbseg)
va6pn=read_cloudvolume_meshes("720575940633169983")

There we were asked to find the volume of a mesh. The first attempt fails

Hide

Rvcg::vcgVolume(va6pn[[1]])

So we need to clean up the mesh:

Hide

va6pn.clean=Rvcg::vcgClean(va6pn[[1]], sel=0:6, iterate = T)
vol=Rvcg::vcgVolume(va6pn.clean)

Hide

vol
# cubic microns 
vol/1e9 

Compare with neuron’s skeleton available from VFB catmaid

Hide

library(elmr)
va6pn.skel=read.neurons.catmaid("name:VA6.*PN")
summary(va6pn.skel)

Hide

cable.length=summary(va6pn.skel)[1,"cable.length"]
vol/cable.length # in nm
vol/cable.length/1e6 # in µm

Hide

A=vol/cable.length/1e6 # mean area µm2
sqrt(A/(2*pi)) # implied radius µm, seems reasonable

Option 1: Vertex to vertex distance

So let’s try and find the vertex to vertex distance for a pair of neurons. We happen to know that this lateral horn local interneuron is a target of the VA6 PN.

Hide

pv4b4=read_cloudvolume_meshes('720575940619630430')

Now we can use a k-d tree to find the nearest neighbour distances very efficiently. Note also the use of the xyzmatrix function to extract vertex locations - this works for all objects that the natverse knows about.

Hide

# better to make the tree from the larger object 
# and query using points from the smaller object
nvertices(pv4b4)
nvertices(va6pn)
kd=nabor::knn(query = xyzmatrix(pv4b4), data=xyzmatrix(va6pn), k=1)
str(kd)

This gives us:

1. a set of vertex indices for the closest match on the VA6 PN for every vertex on the LHLN
2. the distances to those matching vertices

We can take a look at the distribution of the distances like so:

Hide

hist(kd$nn.dists, br=100)

You can see that there are actually quite a lot of close distances.

Hide

hist(kd$nn.dists, br=500, xlim=c(0,2000))

and I assume that the excess of distances <1000 nm where there is a little dip in the distribution is related to this being a pair of neurons that are actually synaptically connected.

Locations of closest approach

OK. So where do those meshes actually approach each other. Well we could look at the locations of the very closest distances <100 nm. How many are there?

Hide

sum(kd$nn.dists<100)

Still a lot! This is because at a location of close approach between the neurons, there will be multiple vertices that are very close together. Separating out these distinct contact locations is a more complex problem than simply calculating all of these vertex distances. The first thing we could do is throw out cases where the same vertex on the target neuron is a match.

Hide

length(unique(kd$nn.idx[kd$nn.dists<100]))

Looks like that would reduce by a third.

Visualising contacts

Let’s just take a quick detour and see what the close approaches look like:

Hide

rgl::setupKnitr()

Hide

nopen3d()
wire3d(pv4b4)
wire3d(va6pn, col='red')
spheres3d(xyzmatrix(pv4b4)[kd$nn.dists<100,], col='green', r=500)

Show close nodes on the query mesh

Hide

nclear3d()
wire3d(pv4b4, col=ifelse(kd$nn.dists<100, 'red', 'black'))
wire3d(va6pn, col='grey')

LHLN + PN

Just the query mesh

Hide

nclear3d()
wire3d(pv4b4, col=ifelse(kd$nn.dists<100, 'red', 'black'))

LHLN with close contact areas in red

Identifying discrete contact zones

There are a range of ways this could be done, essentially what we have is a number of local minima in the distance field to find. It might also be the case that one “patch” on one neuron actually matches several patches on the other neuron. However we’ll just try to pick a simpler approach as an example.

Let’s take the points below threshold distance on the query neuron and find the geodesic distance (i.e. the distance along the mesh surface) between them.

We need to do a bit of prep work for this, to make a graph representation of the mesh.

Hide

#' Make an igraph representation of the edge network of a mesh object
#'
#' @param x A \code{mesh3d} object or an object for which an \code{as.mesh3d}
#'   method is defined. As a convenience the first mesh object is extracted when
#'   \code{x} is a \code{shapelist3d} object.
#' @param ... additional arguments passed to \code{as.mesh3d}
#'
#' @return a \code{igraph} object whose vertices correspond to the vertices of the
#'   original mesh object and whose edges have weights matching the edge lengths
#'   in the mesh.
#' @export
mesh2graph <- function(x, ...) {
  if(inherits(x, 'shapelist3d')) x=x[[1]]
  if(!inherits(x, 'mesh3d')) x=as.mesh3d(x, ...)
  if(nrow(x$it)!=3) stop("I only work for triangulat meshes right now")
  
  # make the edge list
  v1=x$it[1,]
  v2=x$it[2,]
  v3=x$it[3,]
  el=rbind(cbind(v1,v2), cbind(v1,v3), cbind(v2,v3))
  g=igraph::graph_from_edgelist(el, directed = FALSE)
  
  # calculate lengths of each edge
  xyz=xyzmatrix(x)
  edgelengths = function(v1, v2) {
    deltas=xyz[v1,,drop=F]-xyz[v2,,drop=F]
    sqrt(rowSums(deltas^2))
  }
  weights=c(edgelengths(v1,v2), edgelengths(v1,v3), edgelengths(v2,v3))
  igraph::E(g)$weight=weights
  g
}

Convert the LHLN to this mesh representation

Hide

g=mesh2graph(pv4b4)

Now we can find the geodesic distance between all selected nodes.

Hide

dmat=igraph::distances(g, v=which(kd$nn.dists<100), to = which(kd$nn.dists<100))

Let’s try clustering that distance matrix

Hide

dd=hclust(as.dist(dmat))
plot(dd, labels = F)

Hide

# cut using a height relating to a separation of 3000nm
groups=cutree(dd, h=3000)
# sample randomises colours so that neighbouring clusters have different colours
set.seed(42)
groupcols=sample(rainbow(max(groups)))[groups]

nclear3d()
spheres3d(xyzmatrix(pv4b4)[kd$nn.dists<100,], col=groupcols, rad=200)
wire3d(pv4b4, col=ifelse(kd$nn.dists<100, 'red', 'black'))

Close contacts in different colours

Ps the same file as an attachment and the markdown source.

[Greg Jefferis]

[andy cole]

Hi Greg, 

If I want to find the shortest distances between my 3d objects from an Amira .surf file, how would you suggest I modify this rmd file? Is it even possible? Sorry if this question doesn't make sense, the data structure of surfaces from Amira is still a mystery for me after all these years. 

Best, 

Andy

[Gregory Jefferis]

Hi Andy, I only have time to give a couple of pointers, but read in the surfaces using nat,

Use as.mesh3d to convert to the more widely supported mesh3d object type. 

Then use the Rvcg package (which should already installed). The functions 

vcgClost

Or

vcgClostKD

Should help you find closest contacts. 

All the best, Greg. 

Sent from my iPhone

On 7 Dec 2023, at 21:14, andy cole <cole....@gmail.com> wrote:

Hi Greg, 

--
You received this message because you are subscribed to the Google Groups "nat-user" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nat-user+u...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/nat-user/b2d53b43-400c-49aa-8320-3501dd8f1067n%40googlegroups.com.

[andy cole]

Thats exactly what I was looking for. Thanks for this Greg, and all the work you do! It's extremely valuable!

Best,

Andy

[andy cole]

Unfortunately, I this method didn't work. I assume it's because we're exporting as binary surfaces (perhaps their not considered closed surfaces?). However, I was able to workout a workaround, in case someone else has issues with getting Amira files to work in a normal way. Below is the backbone of the method, but essentially, I make a grid of points from vertices coords. Then one doesn't need to deal with the meshes anymore, just math. It's an estimate, but it's accurate to within our resolution limit, IMO.

# Convert all surfaces to meshes

library(nat)
library(Rvcg)

# Initialize an empty list to store all mesh3d objects
all_meshes <- list()

# Loop through each hxsurf binary object
for (surf_name in names(df_surfs)) {
  surf <- df_surfs[[surf_name]]
  all_meshes[[surf_name]] <- list()

  # Check if the hxsurf object has regions and create mesh for each region
  if (!is.null(surf$Regions)) {
    for (region_name in names(surf$Regions)) {
      mesh <- as.mesh3d(surf, Regions = region_name)
      all_meshes[[surf_name]][[region_name]] <- mesh
    }
  }
}
rm(surf, mesh)

# Extract the vertices
vertices <- all_meshes$boo$hoo$vb[1:3,]  # Extract only XYZ coordinates
# Select a subset of vertices for the grid
# Here, we choose every nth vertex for simplicity
n <- 100  # Adjust this based on the density of your mesh, lower is higher resolution
grid_points <- vertices[, seq(1, ncol(vertices), by = n)]


# Visualize the actual object (slow)
open3d()
plot3d(all_meshes$boo_membrane.surf$`Post-synaptic_membrane`, col = "blue")
# Visualize the mesh and the grid points
points3d(grid_points[1,], grid_points[2,], grid_points[3,], col = "red", size = 1)
aspect3d(1, 1, zsize)
# Extract coordinates of the grid points
grid_coords <- t(grid_points)

rm(grid_points,vertices)