Movement Classification and Visualization

[zzaidi148]

Does anyone have any guidance as to how I can go about classifying and visualizing fly movements using RStudio as either 'walking', 'micromovement' or 'immobile' as shown in the original ethoscope publication. I'd like to replicate the box plot and frequency plot similar to this figure:

Thank you! @qgeissmann

[qgeissmann]

Hello, I put together an example for you. It actually just uses custom ggplot. Feel free to adapt :smile: :

rm(list=ls())
library(ggetho)
library(sleepr)

# mock data
metadata <- data.frame(id = paste0("toy_experiment|",1:20),
                       region_id = 1:20,
                       condition = c("A", "B", "C", "D")
                       )
dt <- toy_ethoscope_data(metadata, duration=days(1))
dt <- dt[, sleep_annotation(.SD), by=id]

# we can express x relatively to the food/...
dt[,x_rel:=ifelse(xmv(region_id) > 10, 1-x, x)]

# custom behaviour definition
dt_ord_2 <- dt[,
               {
                 dd <- copy(.SD[, .(t, x_rel, moving)])
                 dd[, t := (floor(t/60))*60]
                 dd[,
                    .(
                      x_rel=mean(x_rel),
                      walked_dist = sum(abs(diff(x_rel))),
                      behaviour = ifelse(all(!moving),1,2)
                    ),
                    by="t"]
               },
               by="id"]

dt_ord_2[, behaviour := ifelse(behaviour != 1, (walked_dist > .25) + 2  , 1)]
dt_ord_2[, behaviour:= ordered(c("q","m","w")[behaviour], levels = c("q","m","w"))]
dt_ord_2[, micro.mov. := (behaviour == "m")]
dt_ord_2[, walking := (behaviour == "w")]
dt_ord_2[, quiescent := (behaviour == "q")]
dt_ord_2[, label:= as.factor(xmv(region_id))]

# this could also be a named list, which is better
BEHAVIOUR_STATES_PALETTE <- c("#999999ff", "#4e9a06ff", "#0070b0ff")

png(h=1920/2,w=1080/2)
pl_exple_behaviour  <- 
  ggplot(dt_ord_2, aes(t ,x_rel)) + 
  geom_rect(mapping=aes(xmin=(t) , xmax=(60+t),fill=behaviour), ymin=-1, ymax=1, alpha=.90)+
  geom_line(size = .5) +
  scale_fill_manual(values=BEHAVIOUR_STATES_PALETTE) +
  facet_grid(label ~ .)  + 
  scale_x_hours() + scale_y_continuous(breaks=c(0,  0.5), name="Position (rel.)") + 
  theme(panel.spacing = unit(0, "cm"), strip.text.y = element_text(angle = -0)) +
  coord_cartesian(ylim = c(0,1))
pl_exple_behaviour
dev.off()

That gives that:

[zzaidi148]

Thank you @qgeissmann! This worked beautifully! What I am now hoping to do is make a similar box plot to how you did in your PLOS paper:

Also, I would greatly appreciate it if you knew of some way we could develop a heatmap for individual flies in R based on the velocity data used in the previous plot. Here's an example of what I mean and the R code that goes along with it adapted from AnimApp(another drosophila tracking software). https://github.com/sraorao/AnimApp/blob/master/fly_plot.R https://www.nature.com/articles/s41598-019-48841-7#Sec2

[zzaidi148]

Also @qgeissmann I think more importantly for my lab would be to understand how to visualize the following plots rather than the ones I asked about earlier (Figures A, D, and E are of particular interest):

I'm terribly sorry for inundating you with questions!!!

[antortjim]

Adding to @qgeissmann beautiful MWE, here is the code that would produce the boxplot @zzaidi is asking for.

• The order of the behaviors is controlled by the order of the letters in the vector passed to levels =.
• The velocity can be interpreted as scaled distance because the sampling points are all equidistant in time.
• I like specifying the namespace of the rethomics functions (with ::) so it's clear where they are coming from, but it's not needed as long as the library is loaded with library().
• On this toy data there seems to be no separation between m and w in terms of maximum velocity, but maybe your data does have a separation.
• Finally, I would like to add the following clarification, @qgeissmann correct me if I am wrong Micromovement here is defined as a scenario where a fly may have moved in some 10 second window (total 6 windows in 1 minute), but where the total distance walked in the minute is less than 0.25. If it's more than 0.25 the behaviour is scored as wake for that minute, and if the fly did not move on any window, it is considered as sleep for the minute. Here we are skipping the 5 minute non movement definition of sleep, even though it is applied in the asleep column in the output of sleep_annotation, every behavior in the plot is derived from the moving column, which does not take this into account. I guess 0.25 is an empirical number. If I am correct, it means 0.25 pixels. Keep in mind that the fly can "move" 0.25 pixels because the position of the fly in the tracking algorithm is represented by the center of mass of the segmented foreground, and the center of mass has subpixel resolution. i.e. if you have a blob in pixels 1, 2, 3, 4, the center of mass (mean) is in 2.5, and same applies for 2D. If you know how much real world distance 1 or 100 or whatever px is, then you can extrapolate how many mm is 0.25 px ;).

rm(list=ls())
library(ggetho)
library(sleepr)
library(behavr)

# mock data
metadata <- data.frame(id = paste0("toy_experiment|",1:20),
                       region_id = 1:20,
                       condition = c("A", "B", "C", "D")
)
dt <- behavr::toy_ethoscope_data(metadata, duration=days(1))
dt <- dt[, sleepr::sleep_annotation(.SD), by=id]

# we can express x relatively to the food/...
dt[,x_rel:=ifelse(xmv(region_id) > 10, 1-x, x)]

# custom behaviour definition
dt_ord_2 <- dt[,
               {
                 dd <- copy(.SD[, .(t, x_rel, moving, max_velocity)])
                 dd[, t := (floor(t/60))*60]
                 dd[,
                    .(
                      x_rel=mean(x_rel),
                      walked_dist = sum(abs(diff(x_rel))),
                      behaviour = ifelse(all(!moving),1,2),
                      max_velocity = max(max_velocity) # take the 60 maximum out of the 10 second maximums
                    ),
                    by="t"]
               },
               by="id"]

dt_ord_2[, behaviour := ifelse(behaviour != 1, (walked_dist > .25) + 2  , 1)]
dt_ord_2[, behaviour:= ordered(c("q","m","w")[behaviour], levels = c("q","m","w"))]
dt_ord_2[, micro.mov. := (behaviour == "m")]
dt_ord_2[, walking := (behaviour == "w")]
dt_ord_2[, quiescent := (behaviour == "q")]
dt_ord_2[, label:= as.factor(xmv(region_id))]

# this could also be a named list, which is better
BEHAVIOUR_STATES_PALETTE <- c("#999999ff", "#4e9a06ff", "#0070b0ff")

png(h=1080/2,w=1080/2)
ggplot(data = dt_ord_2, aes(x = behaviour, y = max_velocity, col=behaviour)) +
  geom_boxplot() +
  # comment if you dont want log10 transformation
  scale_y_continuous(trans = "log10") +
  geom_hline(yintercept = 1, linetype="dashed", size=1) +
  # geom_hline(yintercept = ?, linetype="dashed", size =2) +
  geom_jitter(alpha=0.2) +
  labs(y = "Maximal velocity (rel units)", x = "Behaviour",
       caption = "1 rel unit = velocity correction coef fraction (0.003 or 0.3%) of length of tube")
dev.off()

[qgeissmann]

@antortjim yes. thanks for elaborating, you are right about your last point, but the 0.25 (and the position in general), is expressed as a fraction of the tube (ROI) width. It is an empirical number (the distribution of the distance moved given movement is bimodal, and this threshold splits the modes). I described this in my PhD thesis actually, see figure 4.7: https://spiral.imperial.ac.uk/bitstream/10044/1/69514/3/Geissmann-Q-2018-PhD-Thesis.pdf

[zzaidi148]

Thank you @antortjim and @qgeissmann! I am wondering though why my data seems incredibly close together as opposed to what you have presented. These are my WT flies under normal conditions, tracked for 96 hours: