In jefferis/elmr: Support for working with light and electron microscopy fly brain data

Intro

One of the core features of elmr is to enable data in the FAFB EM space to be transformed into one of a number of standard light level template brains. This transformation depends on manually selected landmarks. We can inspect these landmarks and to try to measure their consistency.

Setup

library(elmr)
library(knitr)
opts_chunk$set(fig.width=5, fig.height=5)
# set up for 3d plots based on rgl package
rgl::setupKnitr()
# frontal view
view3d(userMatrix=rgl::rotationMatrix(angle = pi, 1,0,0), zoom=0.7)

Landmarks

The standard landmarks used by elmr were specified by Davi Bock on the FAFB and JFRC2013 brains using the elm tool developed by Stefan Saalfeld and John Bogovic. We can plot the landmarks in the context of the JFRC2013 brain as follows.

plot3d(JFRC2013)
# note that the landmarks are in raw voxel coordinates and must be scaled
xyz=scale(elm.landmarks[,c("X","Y","Z")], scale = 1/voxdims(JFRC2013), center = FALSE)
xyz=as.data.frame(xyz)
rownames(xyz)=sub("Pt-","",elm.landmarks$Label)
spheres3d(xyz, col = ifelse(elm.landmarks$Use, "green", 'red'), radius = 4)

Landmarks based transform

The standard elm.landmarks set is used within to the elmr package to register a transformation that can map FAFB->JFRC2013 space in conjunction with the xform_brain functions provided by the nat.templatebrains packages.

Let's first look at the FAFB landmarks

# note that these positions are in nm units and do not need scaling
xyz.fafb=elm.landmarks[,c("X1","Y1","Z1")]
rownames(xyz.fafb)=sub("Pt-","",elm.landmarks$Label)
kable(head(xyz.fafb))

Now we can try to look at the consistency of the transform between each landmark pair as follows:

1. Drop out each landmark pair in turn
2. Calculate position of that landmark in target space using remaining landmark pairs
3. Measure difference between calculated position and manually defined position

# record transformed position with leave one out
xyzt.loo=xyz
for(i in 1:nrow(xyz)) {
  thisreg <- reglist(tpsreg(xyz.fafb[-i, ], xyz[-i, ]))
  xyzt.loo[i,]=xform(xyz.fafb[i, , drop=FALSE], thisreg)
}
deltas=sqrt(rowSums((xyzt.loo-xyz)^2))
hist(deltas)

Let's looks at the deltas (distance between leave one out vs manually defined landmark position) for each point.

kable(cbind(elm.landmarks, delta=deltas))

And then visualise them in 3D, colouring each point by the magnitude of the delta.

clear3d()
view3d(userMatrix=rgl::rotationMatrix(angle = pi, 1,0,0), zoom=0.7)
jet.colors <-
  colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                     "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

plot3d(JFRC2013)
spheres3d(xyz, col = jet.colors(10)[cut(deltas, 10)], radius = 4)
# nb these are numbered from 0 like the elm landmarks
texts3d(xyz, texts = rownames(xyz), adj = c(1,1))

Based on this, I would definitely suggest checking the following landmarks

75/48 43 46 68 23

Compare the leave one out landmark position with manually defined position

clear3d()
view3d(userMatrix=rgl::rotationMatrix(angle = pi, 1,0,0), zoom=0.7)
for(i in 1:nrow(xyz)){
  segments3d(rbind(xyz[i, , drop=FALSE], xyzt.loo[i, , drop=FALSE]), col=jet.colors(10)[cut(deltas, 10)][i], lwd=3)
}
points3d(xyz)
texts3d(xyz, texts = rownames(xyz), adj = c(1,1))
plot3d(JFRC2013)

Now of course this is not just a measure of consistency in each landmark point, but it will also depend strongly on the extent to which the neighbouring points can define the required registration; in areas requiring substantial non-rigid deformation this may be a problem; this will be further emphasised when the neighbouring points are distant. We can check the latter point like so:

library(nabor)
nearest3=rowMeans(nabor::knn(xyz, k=4)$nn.dists[,-1])
plot(deltas~nearest3)
deltabynearest3=lm(deltas~nearest3)
abline(deltabynearest3)
summary(deltabynearest3)

So there is a correlation between neighbour distance and our calculated delta; with an R^2 of about 0.4 it's reasonably strong. This may be explained by Davi placing more landmarks in areas where he is more certain.

Perhaps the ideal thing would be to find a set of landmarks on one side of the brain and then find exactly the same landmarks on the other side of the brain. For the JFRC2013 brain we can exactly mirror points from one side to the other using an image-based registration distributed with the nat.flybrains package.

We can find the mirror image position of all the points that we used like so:

#  nb this is conditioned on being able to find CMTK
if(isTRUE(nzchar(cmtk.bindir()))){
  xyzm=mirror_brain(xyz, brain = JFRC2013)

  clear3d()
  spheres3d(xyz, col='red', rad=2)
  spheres3d(xyzm, col='green', rad=2)
  plot3d(JFRC2013)
  view3d(userMatrix=rgl::rotationMatrix(angle = pi, 1,0,0), zoom=0.7)

}

One could then write out a table of landmarks including these mirrored positions like so:

xyzm.pixels=scale(xyzm, scale = voxdims(JFRC2013), center = FALSE)
colnames(xyzm.pixels)=c("XM","YM","ZM")
write.table(cbind(elm.landmarks, xyzm.pixels), 'elm.landmarks-with-mirror.tsv',
            col.names = T, sep='\t', row.names = F)

If one were then to take the mirror image landmarks in JFRC2013 space and manually pick the corresponding locations in EM space, one should be able to compare two different paths to the same point (one mirroring in FAFB, the other mirroring in JFRC2013) and see whether there is a difference in the calculated position. Differences would presumably reflect the uncertainty in picking positions in the EM volume (and a small contribution from the mirroring uncertainty in JFRC2013, which nevertheless expect to be <1 µm based on analysis in Manton et al 2014).

jefferis/elmr documentation built on Sept. 9, 2023, 1:54 p.m.