How to analyze raw USV data

Rodent communication can be measured by recording calls made in the ultrasonic range. Mice pups will start emitting ultrasonic vocalizations (USV) shortly after birth, peaking between PND8-12. These calls elicit pup retrieval by the dam, maternal licking, and crouching behavior. To quantitatively assess social communication in pups we can record and measure the number, duration, frequency, and amplitude of calls. To record the calls, each pup is placed by itself into a plexiglass chamber attached to a microphone. The recordings are later analyzed with an appropriate computer program (UltraVox). From UltraVox, we can get summary statitics (i.e. number of calls, average call duration, etc). We can also extract the duration for every single call made by each pup during the recording period, which can be more informative, and look at patterns of distribution, rathern than summary statistics, as described here: https://github.com/GRousselet/rogme/blob/master/docs/hsf.md.

load your packages:

suppressPackageStartupMessages({
  library(ggplot2)  
  library(lme4)
  library(lmerTest)
  library(graphics)
  library(rogme)
  library(summarytools)
  library(compareGroups)
  library(tidyverse)
})

load your demographics data:

demographics = read.csv('demographics.csv')

load in your usv data. this should be a .csv file with a column for subject id and a column for call duration (i.e. each call has an entry)

usv_data = read.csv("allcalls.csv")

Merge your data frames:

data=merge(demographics, usv_data, by="subject_id",all.y=TRUE)

Do any releveling or setting of reference groups if you have both sexes or a treatment:

data$new_groups= relevel(data$new_groups, ref = "SAL") 
data$Sex= relevel(data$Sex, ref = "Male") 
data$Cohort= as.factor(data$Cohort)

Filter calls below 300 ms in duration, anything above that is typically considered noise (see Scattoni et al., 2008: https://europepmc.org/article/med/18771687)

filtered.data=subset(data, data$duration <300)

Now we can use the standard shift function (https://github.com/GRousselet/rogme/blob/master/docs/hsf.md) to compare distributions:

Turn your data frame into a tibble:

filtered.data$tx=interaction(filtered.data$new_group)
filtered.data.tibble <- as_tibble(filtered.data)

Set the comparisons you want to make. In my case, I had 3 groups, SAL, POL E, and POL L:

sf <- shifthd_pbci(data = filtered.data.tibble, formula = duration ~ tx, todo = list(c("SAL","POL_E"),c("SAL","POL_L"),c("POL_L","POL_E")))
sf

Plot your results based on deciles of distribution:

plot_sf(sf)
p = plot_sf(sf)
p = add_sf_lab(p, sf=sf)

sf2 <- shifthd_pbci(data = filtered.data.tibble, formula = duration ~ new_groups, todo = list(c("POL_E","SAL"),c("POL_L","SAL"),c("POL_L","POL_E")))
sf2
plot_sf(sf2)
p = plot_sf(sf2)
p = add_sf_lab(p, sf=sf2)
p


p <- plot_scat2(data = filtered.data.tibble, formula = duration ~ new_groups,
                xlabel = "",
                ylabel = "Call length (ms)",
                alpha = 1,
                shape = 21,
                colour = "grey20",
                fill = "grey90",
                size = 3) +
  scale_x_discrete(breaks=c("SAL", "POL_E", "POL_L"),
                   labels=c("SAL","POL_E", "POL_L")) +
  theme(axis.text.y = element_text(angle = 90, hjust = .5))
p <- plot_hd_bars(p,
                  col = "black",
                  q_size = 0.5,
                  md_size = 1.5,
                  alpha = 1)
p <- p + coord_flip() #> flip axes
pscat <- p
pscat