---
title: Detecting rat ultrasonic vocalizations
subtitle: Curating count data
author: <a href="https://marce10.github.io/">Marcelo Araya-Salas, PhD</a>
date: "`r Sys.Date()`"
toc: true
toc-depth: 3
toc-location: left
number-sections: true
highlight-style: pygments
format:
html:
df-print: kable
code-fold: true
code-tools: true
editor_options:
chunk_output_type: console
---
::: {.alert .alert-success}
# Objetive {.unnumbered .unlisted}
- Curate and double check OH data
:::
```{r add link to github repo}
#| eval: true
#| output: asis
#| echo: false
# print link to github repo if any
if (file.exists("./.git/config")){
config <- readLines("./.git/config")
url <- grep("url", config, value = TRUE)
url <- gsub("\\turl = |.git$", "", url)
if (nchar(url) > 1)
cat("\nSource code and acoustic data found at [", url, "](", url, ")", sep = "")
}
```
# Load packages {.unnumbered .unlisted}
```{r functions and global parameters}
#| eval: true
#| echo: true
#| message: false
#| warning: false
## add 'developer/' to packages to install from github
x <- c(
"parallel",
"kableExtra",
"knitr",
"DT",
"viridis",
"ggplot2",
"warbleR"
)
sketchy::load_packages(x)
opts_knit$set(root.dir = "..")
options(knitr.kable.NA = '')
print <- function(x) {
kb <- kable(x, row.names = FALSE, digits = 4, "html")
kb <- kable_styling(kb,
bootstrap_options = c("striped", "hover", "condensed", "responsive"))
scroll_box(kb, width = "100%")
}
theme_set(theme_classic())
options(
DT.options = list(
pageLength = 200,
scrollX = TRUE,
scrollY = "800px",
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel')
)
)
datatable2 <- function(x, ...) datatable(data = x, extensions = 'Buttons', height = 200, ...)
source("~/Dropbox/R_package_testing/warbleR/R/acoustic_activity.R")
source("~/Dropbox/R_package_testing/warbleR/R/internal_functions.R")
```
# Read data
```{r}
#| eval: true
#| echo: true
# read bout 1
bout1_dat <- readRDS(
file.path(
"/home/m/Dropbox/Projects/rat_vocalization_alcohol/data/processed/",
"bout1_USV_counts_per_minute_and_detections.RDS"
)
)
# read bout 2
bout2_dat <- readRDS(
file.path(
"/home/m/Dropbox/Projects/rat_vocalization_alcohol/data/processed/",
"bout2_USV_counts_per_minute_and_detections.RDS"
)
)
# table(bout2_dat$wide_count_min$experiment)
# table(bout1_dat$wide_count_min$experiment)
obesidad_counts <- bout2_dat$wide_count_min[bout2_dat$wide_count_min$experiment == "Audios Obesidad 3", ]
oh_counts <- rbind(bout1_dat$wide_count_min, bout2_dat$wide_count_min[bout2_dat$wide_count_min$experiment != "Audios Obesidad 3", ])
oh_metadata <- rbind(bout1_dat$metadata, bout2_dat$metadata[grep("Obesidad 3", bout2_dat$metadata$original_dir, invert = TRUE), ])
# remove those shorter than 1 min
oh_counts <- oh_counts[oh_counts$sound.file.duration >= 60, ]
# table(oh_counts$experiment)
```
# Add metadata
Observations:
* bout 172 only has 2 cages
```{r}
oh_counts$sound.file.duration.min <- round(oh_counts$sound.file.duration/ 60, 1)
oh_counts$date <- sapply(oh_counts$sound.files, function(x) {
# extract date
date <- strsplit(split = "202",
x = x,
fixed = TRUE)[[1]]
date <- ifelse(grepl("EXP2_OH", date[length(date)]), date[2], date[length(date)])
date <- paste0("202", substr(date, 1, 7))
date
})
# unique(oh_counts$date)
oh_counts$date <- as.Date(oh_counts$date)
oh_counts$time <- sapply(oh_counts$sound.files, function(x) {
# extract date
date <- strsplit(split = "202",
x = x,
fixed = TRUE)[[1]]
date <- ifelse(grepl("EXP2_OH", date[length(date)]), date[2], date[length(date)])
date <- substr(date, 9, 16)
date
})
# unique(oh_counts$time)
oh_counts <- oh_counts[order(oh_counts$date, oh_counts$time), ]
oh_counts$directory <- sapply(oh_counts$sound.files, function(x){
oh_metadata$original_dir[oh_metadata$new_name == x]
}
)
oh_counts$cage <- substr(oh_counts$sound.files, 1, 1)
oh_counts$cage[oh_counts$cage == "E"] <- NA
oh_counts$cage <- ifelse(oh_counts$cage == "2",substr(basename(oh_counts$directory),1, 1), oh_counts$cage)
oh_counts$cage <- ifelse(grepl("CT_PRUEBA",oh_counts$sound.files),substr(oh_counts$sound.files, 11, 11), oh_counts$cage)
# table(oh_counts$cage)
num_time <- as.numeric(gsub("-", "", oh_counts$time))
oh_counts$bout <- NA
oh_counts$bout[1] <- 1
for (i in 2:nrow(oh_counts)) {
if (!is.na(num_time[i])){
if (num_time[i] - num_time[i - 1] <= 41 & oh_counts$date[i] == oh_counts$date[i - 1]){
oh_counts$bout[i] <- oh_counts$bout[i - 1]
} else {
oh_counts$bout[i] <- oh_counts$bout[i - 1] + 1
}
} else{ oh_counts$bout[i] <- NA}
}
tab <- table(oh_counts$bout)
# all(tab == 4)
#
# which(tab < 4)
num_time <- as.numeric(gsub("-", "", oh_counts$time))
oh_counts$day.bout[1] <- NA
oh_counts$day.bout[1] <- 1
for (i in 2:nrow(oh_counts)) {
if (!is.na(oh_counts$date[i])){
if (oh_counts$date[i] != oh_counts$date[i - 1]) {
oh_counts$day.bout[i] <- 1
} else {
if (num_time[i] - num_time[i - 1] <= 41 & oh_counts$date[i] == oh_counts$date[i - 1]){
oh_counts$day.bout[i] <- oh_counts$day.bout[i - 1]
} else {
oh_counts$day.bout[i] <- oh_counts$day.bout[i - 1] + 1
}
}
} else{ oh_counts$day.bout[i] <- NA}
}
tab <- table(oh_counts$day.bout, oh_counts$date)
# all(tab <= 4)
oh_counts$period <- ifelse(oh_counts$sound.file.duration < 660, "pre", "post")
oh_counts$period <- ifelse(grepl("Día 0", oh_counts$day), "day 0", oh_counts$period)
# table(oh_counts$period)
min_cols <- paste("min", 1:30)
# order columns
oh_counts <- oh_counts[,c(
setdiff(names(oh_counts), min_cols),
min_cols
)]
```
# Graphs summarizing data
```{r}
agg <- aggregate(total ~ period + day + experiment, data = oh_counts, length)
agg$day_num <- sapply(strsplit(x = gsub("Día ", "", agg$day), split = " ", fixed = TRUE), "[[", 1)
agg$day_num <- factor(agg$day_num, levels = sort(as.numeric(unique(agg$day_num))))
agg$period <- factor(agg$period, levels = c("day 0", "pre", "post"))
# bar graph colored by period
ggplot(agg, aes(x = day_num, y = total, fill = period)) +
geom_bar(stat = "identity", position = "dodge") +
labs(
title = "Total recordings per day",
x = "Day",
y = "# of recordings"
) +
scale_fill_viridis_d(option = "G", alpha = 0.9) +
facet_wrap(~ experiment, ncol = 2)
```
```{r}
agg <- aggregate(total ~ period + day + experiment, data = oh_counts, sum)
agg$day_num <- sapply(strsplit(x = gsub("Día ", "", agg$day), split = " ", fixed = TRUE), "[[", 1)
agg$day_num <- factor(agg$day_num, levels = sort(as.numeric(unique(agg$day_num))))
agg$period <- factor(agg$period, levels = c("day 0", "pre", "post"))
# bar graph colored by period
ggplot(agg, aes(x = day_num, y = total, fill = period)) +
geom_bar(stat = "identity", position = "dodge") +
labs(
title = "Total USVs counts per day all minutes",
x = "Day",
y = "# of USVs"
) +
scale_fill_viridis_d(option = "G", alpha = 0.9) +
facet_wrap(~ experiment, ncol = 2)
```
```{r}
post_last15 <- oh_counts[oh_counts$period == "post", c("experiment", "day", paste("min", 15:30))]
post_last15$total <- rowSums(post_last15[, paste("min", 15:30)], na.rm = TRUE)
agg <- aggregate(total ~ day + experiment, data = post_last15, sum)
agg$day_num <- sapply(strsplit(x = gsub("Día ", "", agg$day), split = " ", fixed = TRUE), "[[", 1)
agg$day_num <- factor(agg$day_num, levels = sort(as.numeric(unique(agg$day_num))))
# bar graph colored by period
ggplot(agg, aes(x = day_num, y = total)) +
geom_bar(stat = "identity", position = "dodge", fill = "#54C9ADB3") +
labs(
title = "Total USV counts per day",
x = "Day",
y = "Total USV counts"
) +
facet_wrap(~ experiment, ncol = 2)
```
# Activity by minute
```{r}
detections <- rbind(bout1_dat$detections, bout2_dat$detections)
file_durations <- rbind(bout1_dat$wide_count_min[, c("sound.files", "sound.file.duration")], bout2_dat$wide_count_min[, c("sound.files", "sound.file.duration")])
names(file_durations)[2] <- "duration"
# counting
count_high_min <- acoustic_activity(
X = detections,
time.window = 60,
hop.size = 1,
Y = file_durations
)
# rate per minute
count_high_min$rate <- count_high_min$rate * 60
## add metadata
all(count_high_min$sound.files %in% oh_counts$sound.files)
count_high_min <- count_high_min[count_high_min$sound.files %in% oh_counts$sound.files,]
setdiff(unique(count_high_min$sound.files), unique(oh_counts$sound.files))
count_high_min$experiment <- sapply(count_high_min$sound.files, function(x) {
oh_counts$experiment[oh_counts$sound.files == x][1]
})
count_high_min$day <- sapply(count_high_min$sound.files, function(x) {
oh_counts$day[oh_counts$sound.files == x][1]
})
count_high_min$period <- sapply(count_high_min$sound.files, function(x) {
oh_counts$period[oh_counts$sound.files == x][1]
})
# get highest per sound file
sub_count_high_min <- count_high_min[count_high_min$duration == 60, ]
sub_count_high_min_list <- lapply(unique(sub_count_high_min$sound.files), function(x){
Y <- sub_count_high_min[sub_count_high_min$sound.files == x,]
Y$norm.rate <- Y$rate / max(Y$rate, na.rm = TRUE)
return(Y)
})
sub_count_high_min <- do.call(rbind, sub_count_high_min_list)
sub_count_high_min$period <- factor(sub_count_high_min$period, levels = c("day 0", "pre", "post"))
# aggregate by start, day and period
agg_count <- aggregate(norm.rate ~ start + period + experiment, data = sub_count_high_min, mean)
tab <- table(agg_count$start)
agg_count <- agg_count[agg_count$start %in% names(tab)[tab > 3], ]
sub_count_high_min <- sub_count_high_min[sub_count_high_min$start %in% names(tab)[tab > 3],]
sub_count_high_min$experiment <- factor(sub_count_high_min$experiment, levels = c("2020 OH3 M-SENS - USVs", "2022 OH5 H-SENS - USVs", "2022 OH4 M-INCUB - USVs", "2023 OH6 H-INCUB - USVs"))
# ggplot of rate per start faceted by sound file
ggplot(sub_count_high_min, aes(x = start / 60, y = norm.rate, group = sound.files, color = period)) +
geom_line() +
facet_wrap( ~ experiment+ period, scales = "free_y", ncol = 3) +
labs(
title = "Rate of USV detections across recordings",
x = "Time (min)",
y = "Normalized rate (USVs/min)"
) +
geom_line(data = agg_count, aes(x = start / 60, y = norm.rate), size = 1, inherit.aes = FALSE) +
scale_color_viridis_d(option = "G", alpha = 0.5, end = 0.8, guide = NULL) +
theme_classic()
```
Highest min at the first, second and third segments of the file:
```{r}
count_sections_list <- lapply(unique(count_high_min$sound.files), function(x){
count_high_min_file <- count_high_min[
count_high_min$sound.files == x,
]
max_dur <- max(count_high_min_file$start)
cut_points <- c(0, max_dur / 3, 2 * max_dur / 3, max_dur)
out_sections <- lapply(1:(length(cut_points) - 1), function(i) {
count_high_min_file_segment <- count_high_min_file[
count_high_min_file$start >= cut_points[i] & count_high_min_file$start < cut_points[i + 1],
]
Y <- count_high_min_file_segment[which.max(count_high_min_file_segment$rate),]
Y$section <- ifelse(i == 1, "begin", ifelse(i == 2, "middle", "end"))
return(Y)
})
out_sections <- do.call(rbind, out_sections)
return(out_sections)
})
count_sections <- do.call(rbind, count_sections_list)
```
```{r}
#| eval: true
#| echo: false
datatable2(count_sections)
```
# Session information {.unnumbered .unlisted}
<details>
<summary>Click to see</summary>
```{r session info}
#| echo: false
# if devtools is installed use devtools::session_info()
if (requireNamespace("devtools", quietly = TRUE)) {
devtools::session_info()
} else {
sessionInfo()
}
```
</details>
