In danielkick/mRNA24hLC:
knitr::opts_chunk$set(echo = TRUE)
library(readxl)
library(here)
library(tidyverse)
library(janitor)
library(zoo)
library(readABF)
library(broom)
library(minpack.lm) # getting Tau
library(patchwork)
library(ggsci)
library(cowplot)
library(patchwork)
library(minpack.lm) # for use in fit_tau_1term_exp()
theme_set(ggplot2::theme_minimal())
library(devtools)
devtools::load_all()
abf_dir_path <- "D:/0000_Traces_mRNA24hLC/"
abf_dir_path <- "./inst/extdata/ABFs/"
Handling different kinds of files (Proof of principle)
Load in data (ionic currents)
file_list <- list(FI = "190903_0025.abf",
EPSP = "190903_0023.abf",
GJCC = "190903_0021.abf",
GJVC = "190903_0020.abf",
HTK = "190903_0031.abf",
HTKC = "190903_0034.abf",
A = "190903_0030.abf",
AC = "190903_0033.abf")
trace_list <- map(seq_along(file_list), function(i){
temp <- readABF_as_matrix2(path = paste0(abf_dir_path,
# "C:/Users/Daniel/Documents/Trace_Holding/",
file_list[[i]]),
channels = "all",
relative.time = T) %>%
as.data.frame() %>%
janitor::clean_names(case = "upper_camel") %>%
mutate(Sweep = as.factor(Sweep)) %>%
group_by(Sweep) %>%
mutate(MinTime = min(Time, na.rm = T)) %>%
mutate(Time = Time - MinTime) %>%
select(-MinTime)
})
names(trace_list) <- names(file_list)
Show all
example_plots <- map(trace_list, function(i){
downsample_data(i, len = 5000) %>%
ungroup() %>%
gather(Ch, Value, c("In4","In7","In9","In12")) %>%
mutate(Cell = case_when(Ch %in% c("In4","In7") ~ "Cell A",
Ch %in% c("In9","In12") ~ "Cell B",)) %>%
mutate(ChType = case_when(Ch %in% c("In4","In9") ~ "mV",
Ch %in% c("In7","In12") ~ "nA",)) %>%
mutate(Ch = factor(.$Ch, levels = c("In4","In7","In9","In12"))) %>%
ggplot(aes(Time, Value, group = Sweep, color = Cell))+
geom_line()+
facet_wrap(ChType~Cell, scales = "free_y", ncol = 2)+
labs(y = "")+
theme(legend.position = "")+
scale_color_aaas()
})
plot_grid(plotlist = example_plots)
GJVC
## GJVC: get Ig, Ileak====
process_gjvc_trace <- function(trace = trace_list$GJVC){
### Setup
i_leak <- trace %>%
# filter(Time > 0.75) %>%
mutate(In4 = as.numeric(In4),
In7 = as.numeric(In7),
In9 = as.numeric(In9),
In12 = as.numeric(In12))
i_leak <- i_leak %>%
#When are the step start times?
# 1. In9: ~0.27575
# 2. In4: ~1.02532
mutate(Inj_in = ifelse(Time > 0.8, "In4", "In9"),
Segment = case_when(
Time >= (0.07575 - 0.02532) & Time <= (0.27575 - 0.02532) ~ "PreStep", # Inj In9
Time >= (0.27575+ 0.17468) & Time <= (0.27575+ 0.22468) ~ "Step",
Time >= 0.8 & Time <= 1.0 ~ "PreStep", # Inj In4
Time >= 1.2 & Time <= 1.25 ~ "Step"
))
plt_tagged_regions <- downsample_data(i_leak, len = 10000) %>%
gather("key", "value", c("In4", "In7", "In9", "In12")) %>%
mutate(Cell = case_when(key %in% c("In4", "In7") ~ "Cell1",
key %in% c("In9", "In12") ~ "Cell2"),
Aspect = case_when(key %in% c("In4", "In9") ~ "mV",
key %in% c("In7", "In12") ~ "nA")
) %>%
ggplot(aes(Time, value, color = Segment, group = Sweep))+
geom_path()+
# facet_wrap(key~., scales = "free")
facet_grid(Aspect~Cell, scales = "free")
i_leak_for_In4_to_In9 <- i_leak %>%
filter(Inj_in == "In4") %>%
group_by(Sweep, Segment, Inj_in) %>%
mutate(In4Mean = mean(In4, na.rm = T),
In7Mean = mean(In7, na.rm = T),
In9Mean = mean(In9, na.rm = T),
In12Mean =mean(In12, na.rm = T)) %>%
ungroup()
i_leak_for_In9_to_In4 <- i_leak %>%
filter(Inj_in == "In9") %>%
group_by(Sweep, Segment, Inj_in) %>%
mutate(In4Mean = mean(In4, na.rm = T),
In7Mean = mean(In7, na.rm = T),
In9Mean = mean(In9, na.rm = T),
In12Mean =mean(In12, na.rm = T)) %>%
ungroup()
### Ig ####
i_leak_for_In4_to_In9 <- i_leak_for_In4_to_In9 %>%
select(Sweep, Segment, In4Mean, In7Mean, In9Mean, In12Mean) %>%
filter(Segment != "NA") %>%
distinct() %>%
pivot_wider(names_from = "Segment",
values_from = c("In4Mean", "In12Mean",
"In9Mean", "In7Mean" )) %>%
mutate(In4Mean_Step = In4Mean_Step - In4Mean_PreStep,
In12Mean_Step = In12Mean_Step - In12Mean_PreStep,
In9Mean_Step = In9Mean_Step - In9Mean_PreStep,
In7Mean_Step = In7Mean_Step - In7Mean_PreStep
) %>%
select(Sweep, In4Mean_Step, In12Mean_Step,
In9Mean_Step, In7Mean_Step) %>%
rename(In4Mean = In4Mean_Step, In12Mean = In12Mean_Step,
In9Mean = In9Mean_Step, In7Mean = In7Mean_Step) %>%
filter(abs(In4Mean) >= 4 | abs(In9Mean) >= 4 )# exclude no step condition
i_leak_for_In9_to_In4 <- i_leak_for_In9_to_In4 %>%
select(Sweep, Segment, In4Mean, In7Mean, In9Mean, In12Mean) %>%
filter(Segment != "NA") %>%
distinct() %>%
pivot_wider(names_from = "Segment",
values_from = c("In4Mean", "In12Mean",
"In9Mean", "In7Mean" )) %>%
mutate(In4Mean_Step = In4Mean_Step - In4Mean_PreStep,
In12Mean_Step = In12Mean_Step - In12Mean_PreStep,
In9Mean_Step = In9Mean_Step - In9Mean_PreStep,
In7Mean_Step = In7Mean_Step - In7Mean_PreStep
) %>%
select(Sweep, In4Mean_Step, In12Mean_Step,
In9Mean_Step, In7Mean_Step) %>%
rename(In4Mean = In4Mean_Step, In12Mean = In12Mean_Step,
In9Mean = In9Mean_Step, In7Mean = In7Mean_Step) %>%
filter(abs(In4Mean) >= 4 | abs(In9Mean) >= 4 )# exclude no step condition
plt_ig_slopes <- ggplot()+
geom_smooth(data = i_leak_for_In4_to_In9, aes(In4Mean, In12Mean), method = lm, color = "red")+
geom_point( data = i_leak_for_In4_to_In9, aes(In4Mean, In12Mean), color = "red")+
geom_smooth(data = i_leak_for_In9_to_In4, aes(In9Mean, In7Mean), method = lm)+
geom_point( data = i_leak_for_In9_to_In4, aes(In9Mean, In7Mean))+
labs(x = "Difference in mV", y = "nA")
In4_to_In9_Ig <- as.numeric(broom::tidy(lm(In12Mean ~ In4Mean, data = i_leak_for_In4_to_In9))[2, "estimate"])
In9_to_In4_Ig <- as.numeric(broom::tidy(lm(In7Mean ~ In9Mean, data = i_leak_for_In9_to_In4))[2, "estimate"])
### Ileak ####
In4_R <- as.numeric(broom::tidy(lm(In7Mean ~ In4Mean, data = i_leak_for_In4_to_In9))[2, "estimate"])
In9_R <- as.numeric(broom::tidy(lm(In12Mean ~ In9Mean, data = i_leak_for_In9_to_In4))[2, "estimate"])
In4_R <- In4_R^-1
In9_R <- In9_R^-1
return(
list(diagnostic_plots = list(
trace = plt_tagged_regions,
fit = plt_ig_slopes),
df = data.frame(
In4_to_In9 = In4_to_In9_Ig,
In9_to_In4 = In9_to_In4_Ig,
In4_R1c = In4_R,
In9_R1c = In9_R
)
)
)
}
out <- process_gjvc_trace(trace = trace_list$GJVC)
out$diagnostic_plots$trace
out$diagnostic_plots$fit
out$df
# # Save for review:
# input_file <- "190903_0020.abf"
#
# write.csv(out$df,
# here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjvc", "df.csv", sep = "-")) )
#
# save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjvc", "trace.pdf", sep = "-")),
# out$diagnostic_plots$trace, base_height = 7, base_width = 10)
#
# save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjvc", "fit.pdf", sep = "-")),
# out$diagnostic_plots$fit, base_height = 7, base_width = 10)
GJCC
For reference, we use the following equations to determine the resistances of interest.
$$r_{11}=\frac{v_1}{i_1}$$
$$r_{12}=\frac{v_2}{i_1}$$
$$r_1= \frac{r_{11}*r_{22} - r_{12}^2}{r_{22} - r_{12}}$$
$$r_c = \frac{r_{11}*r_{22} - r_{12}^2}{r_{12}}$$
## GJCC: get CC, R_C, C, Rin ====
process_gjcc_trace <- function(trace = trace_list$GJCC){
## prep ====
temp <- trace
temp_plt <- temp %>%
# filter(Sweep == 4) %>%
mutate(In4 = as.numeric(In4)) %>%
mutate(In7 = as.numeric(In7))
# Show what the protocol looks like
plt1 <- downsample_data(temp_plt, len = 1000) %>%
gather(key, value, c("In4", "In9")) %>%
ggplot(aes(Time, value, color = key, group = interaction(Sweep, key)))+
geom_path()+
scale_color_manual(values = RColorBrewer::brewer.pal(7, name = "PuOr")[c(1, 7)] )
plt2 <- downsample_data(temp_plt, len = 1000) %>%
gather(key, value, c("In7", "In12")) %>%
# mutate(key = factor(.$key, levels = c("In7", "In12"))) %>%
ggplot(aes(Time, value, color = key, group = interaction(Sweep, key)))+
geom_path()+
scale_color_manual(values = RColorBrewer::brewer.pal(7, name = "PuOr")[c(7, 1)] )# Without flipping the factor level, get the right colors.
plt_trace <- plt1/plt2
# annotate
temp <- temp %>%
mutate(Segment = case_when(Time >= 0 & Time <= 0.32820 ~ "A", #"PreStepIn12",
Time >= 1.32836 & Time <= 1.82836 ~ "B", #"StepIn12",
Time >= 2.57744 & Time <= 3.07744 ~ "C", #"PreStepIn7",
Time >= 4.07760 & Time <= 4.57760 ~ "D", #"StepIn7"
))
## resistances/ coupling ====
# $$r_{11}=\frac{v_1}{i_1}$$
# $$r_{12}=\frac{v_2}{i_1}$$
# $$r_1= \frac{r_{11}*r_{22} - r_{12}^2}{r_{22} - r_{12}}$$
# $$r_c = \frac{r_{11}*r_{22} - r_{12}^2}{r_{12}}$$
temp_resist <- temp %>%
group_by(Sweep, Segment) %>%
filter(!is.na(Segment)) %>%
summarise(MeanIn4 = mean(In4, na.rm = T),
MeanIn7 = mean(In7, na.rm = T),
MeanIn9 = mean(In9, na.rm = T),
MeanIn12 = mean(In12, na.rm = T)) %>%
ungroup() %>%
pivot_wider(names_from = Segment, values_from = c("MeanIn4", "MeanIn7", "MeanIn9", "MeanIn12")) %>%
mutate(In4_CC = MeanIn9_D/MeanIn4_D, #inj in In4
In4_R11 = MeanIn4_D/MeanIn7_D,
In4_R12 = MeanIn9_D/MeanIn7_D,
In9_CC = MeanIn4_B/MeanIn9_B, #inj in In9
In9_R11 = MeanIn9_B/MeanIn12_B,
In9_R12 = MeanIn4_B/MeanIn12_B
) %>%
mutate(mean_R12 = (In4_R12+In9_R12)/2) %>%
mutate(In4_R1 = ((In4_R11*In9_R11) - (mean_R12^2)) / (In9_R11-mean_R12),
In4_Rc = ((In4_R11*In9_R11) - (mean_R12^2)) / (mean_R12),
In9_R1 = ((In9_R11*In4_R11) - (mean_R12^2)) / (In4_R11-mean_R12),
In9_Rc = ((In9_R11*In4_R11) - (mean_R12^2)) / (mean_R12)
) %>% select(starts_with("In"), "Sweep")
## V_rest ====
temp_rmp <- temp %>%
filter(Segment == "C") %>%
group_by(Sweep) %>%
summarise(In4_Vrest = min(In4, na.rm = T),
In9_Vrest = min(In9, na.rm = T))
## Tau ====
temp <- temp %>%
mutate(Segment = case_when(Time >= 3.07744 & Time <= 4.47760 ~ "In4Step", #"PreStepIn7",
Time >= 0.32820 & Time <= 1.72836 ~ "In9Step", #"StepIn7"
))
plt1 <-
downsample_data(temp, len = 1000) %>%
gather(Key, Value, c("In4", "In9")) %>%
ggplot(aes(Time, Value, color = Segment, group = interaction(Key, Sweep)))+
geom_path()
plt2 <-
downsample_data(temp, len = 1000) %>%
ungroup() %>%
gather(key, value, c("In7", "In12")) %>%
mutate(key = factor(.$key, levels = c("In7", "In12"))) %>%
ggplot(aes(Time, value, color = Segment, group = interaction(key, Sweep)))+
geom_path()
plt_segments <- plt1/plt2
### in4
temp_in4 <- temp %>%
filter(Segment == "In4Step") %>%
mutate(Time = Time - min(Time, na.rm = T))
test_fits <- map(unique(temp_in4$Sweep), function(i){
fit_tau_1term_exp(df = filter(temp_in4, Sweep == i),
IV = "Time",
DV = "In4")
})
fits_in4 <- map(seq_along(transpose(test_fits)[[1]]), function(i){
df <- transpose(test_fits)[[1]][i][[1]]
df <- df[df$term == "b", c("estimate", "std.error")]
df$Sweep <- i
return(df)
}) %>%
do.call(rbind, .) %>%
rename(In4_Tau_est = estimate, In4_Tau_std.err = std.error) %>%
mutate(Sweep = as.factor(Sweep))
test_fits_plts <- map(seq_along(test_fits), function(i){
ggplot()+
geom_line(data = as.data.frame(test_fits[[i]]$check_fit), aes(x = Time, y = Fit), color = "firebrick")+
geom_line(data = as.data.frame(test_fits[[i]]$check_fit), aes(x = Time, y = In4), alpha = 0.8)+
labs(subtitle = paste("tau =", as.character(
round(test_fits[[i]]$fit[2, 2], digits = 3)
)))
})
plt_tau_in4 <- cowplot::plot_grid(plotlist = test_fits_plts)
### in9
temp_in9 <- temp %>%
filter(Segment == "In9Step") %>%
mutate(Time = Time - min(Time, na.rm = T))
test_fits <- map(unique(temp_in9$Sweep), function(i){
fit_tau_1term_exp(df = filter(temp_in9, Sweep == i),
IV = "Time",
DV = "In9")
})
fits_in9 <- map(seq_along(transpose(test_fits)[[1]]), function(i){
df <- transpose(test_fits)[[1]][i][[1]]
df <- df[df$term == "b", c("estimate", "std.error")]
df$Sweep <- i
return(df)
}) %>%
do.call(rbind, .) %>%
rename(In9_Tau_est = estimate, In9_Tau_std.err = std.error) %>%
mutate(Sweep = as.factor(Sweep))
test_fits_plts <- map(seq_along(test_fits), function(i){
ggplot()+
geom_line(data = as.data.frame(test_fits[[i]]$check_fit), aes(x = Time, y = Fit), color = "firebrick")+
geom_line(data = as.data.frame(test_fits[[i]]$check_fit), aes(x = Time, y = In9), alpha = 0.8)+
labs(subtitle = paste("tau =", as.character(
round(test_fits[[i]]$fit[2, 2], digits = 3)
)))
})
plt_tau_in9 <- cowplot::plot_grid(plotlist = test_fits_plts)
## Condense ====
summary_df <- full_join(temp_resist, temp_rmp) %>% full_join(., fits_in4) %>% full_join(., fits_in9)
return(list(
diagnostic_plots = list(
trace = plt_trace,
segment = plt_segments,
fit_In4 = plt_tau_in4,
fit_In9 = plt_tau_in9
),
df = summary_df
))
}
out <- process_gjcc_trace(trace = trace_list$GJCC)
out$diagnostic_plots$trace
out$diagnostic_plots$segment
out$diagnostic_plots$fit_In4
out$diagnostic_plots$fit_In9
out$df
# # Save for review:
# input_file <- "190903_0021.abf"
#
# write.csv(out$df,
# here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjcc", "df.csv", sep = "-")))
#
# save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjcc", "trace.pdf", sep = "-")),
# out$diagnostic_plots$trace, base_height = 7, base_width = 10)
#
# save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjcc", "segment.pdf", sep = "-")),
# out$diagnostic_plots$segment, base_height = 7, base_width = 10)
#
# save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjcc", "fit_in4.pdf", sep = "-")),
# out$diagnostic_plots$fit_In4, base_height = 7, base_width = 10)
#
# save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjcc", "fit_in9.pdf", sep = "-")),
# out$diagnostic_plots$fit_In9, base_height = 7, base_width = 10)
FI
temp <- trace_list$FI
downsample_data(df = temp, len = 5000) %>%
gather(key, value, c("In9", "In12")) %>%
ungroup() %>%
mutate(key = factor(.$key, levels = c("In9", "In12"))) %>%
ggplot(aes(Time, value, group = Sweep))+
geom_line()+
facet_wrap(key~., scales = "free_y", ncol = 1)
#r11 == ~15.3
#r1 == ~17.9
place_params_for_predict_voltage_trace(tau = 19.9,
Rin = 15.3)
# place_data_for_predict_voltage_trace(input.df = temp[temp$Sweep == 10, ],
# Time.ch = "Time",
# Inj.ch = "In12")
#
# run_predict_voltage_trace_py()
#
# predicted_response <- retrieve_predicted_voltage()
predicted_response_list <- map(seq(1, length(unique(temp$Sweep))), function(i){
place_data_for_predict_voltage_trace(input.df = temp[temp$Sweep == i, ],
Time.ch = "Time",
Inj.ch = "In12")
run_predict_voltage_trace_py()
predicted_response <- retrieve_predicted_voltage()
return(predicted_response)
})
my_plts <- map(
seq_along(predicted_response_list),
function(i){
ggplot()+
geom_line(data = predicted_response_list[[i]], aes(Time, predicted-56.5))+
geom_line(data = temp[temp$Sweep == i, ], aes(Time, In9), color = "red")+
labs(x = "Seconds", y = "mV")+
coord_cartesian(y = c(-150, +120))
})
# cowplot::plot_grid(plotlist = my_plts)
multi_plt <- my_plts[[11]] / my_plts[[11 - 2]] / my_plts[[1 + 2]]
# # Save for review:
# input_file <- "190903_0025.abf"
# saveRDS(predicted_response_list,
# here("data", "predicted_voltage", paste(str_remove(input_file, ".abf"), "fi", "predict_v.rds", sep = "-")))
#
# save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "fi", "multi_plt.pdf", sep = "-")),
# multi_plt, base_height = 7, base_width = 10)
# test_load <- readRDS(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "fi", "predict_v.rds", sep = "-")))
predict_fi_responses <- function(trace = trace_list$FI,
tau = 19.8, #abs(mean(unlist(out$df[, "In9_Tau_est"]))),
rin = 20.0, #abs(mean(unlist(out$df[, "In9_R1"]))),
v_ch = "In9",
i_ch = "In12"){
# trace = trace_list$FI
# tau = 19.8
# rin = 20.0
# v_ch = "In9"
# i_ch = "In12"
temp <- trace
# downsample_data(df = temp, len = 10000) %>%
# gather(key, value, c(v_ch, i_ch)) %>%
# ungroup() %>%
# mutate(key = factor(.$key, levels = c(v_ch, i_ch))) %>%
# ggplot(aes(Time, value, group = Sweep))+
# geom_line()+
# facet_wrap(key~., scales = "free_y", ncol = 1)
place_params_for_predict_voltage_trace(tau = tau,
Rin = rin)
predicted_response_list <- map(seq(1, length(unique(temp$Sweep))), function(i){
place_data_for_predict_voltage_trace(input.df = temp[temp$Sweep == i, ],
Time.ch = "Time",
Inj.ch = i_ch)
run_predict_voltage_trace_py()
predicted_response <- retrieve_predicted_voltage()
pred_v <- predicted_response[predicted_response$Time >= 4, "predicted"]
meas_v <- temp[temp$Sweep == i & temp$Time >= 4, v_ch]
offset <- median(unlist(pred_v - meas_v), na.rm = T)
predicted_response$offset <- offset
return(predicted_response)
})
my_plts <- map(
seq_along(predicted_response_list),
function(i){
ggplot()+
geom_line(data = predicted_response_list[[i]], aes(Time, predicted-offset), color = "red")+
geom_line(data = temp[temp$Sweep == i, ], aes(Time, In9), color = "black")+
labs(x = "Seconds", y = "mV")+
coord_cartesian(y = c(-150, +120))
})
return(list(diagnostic_plots = list(traces = my_plts),
predictions = predicted_response_list))
}
# my_plts[[8]]
#
# cowplot::plot_grid(plotlist = my_plts)
#
# multi_plt <- my_plts[[11]] / my_plts[[11 - 2]] / my_plts[[1 + 2]]
EPSP
(epsp) How could we slice up the trace?
## Here's the stimulus used in the protocol ====
# sweep duration should be 19.687 seconds
epsp_stim <- readABF_as_matrix(here("inst", "extdata", "epsp_stim_files", "08 current injection.abf"),
channels = "Axo1I2")
epsp_stim <- as_tibble(epsp_stim) %>%
mutate(Time = Time - min(Time, na.rm = T)) %>%
rename(Stim = Axo1I2)
shading_annotations <- data.frame(
starts = c(0.40,
4.87,
10.74,
16.74),
next_start = c(4.87,
10.74,
16.74,
19.685),
equal_len = c(0.40,
4.87,
10.74,
16.74) + 2.945,
on_end = c(2,
6.17,
12.54,
18.5),
equal_on = c(0.40,
4.87,
10.74,
16.74) + 1.30
)
annotation_df <- data.frame(x = c(0.5, 0.5, 2.1, 1.8),
y = seq(0.25, -1.25, length.out = 4),
text = c("Period",
"Min. Period",
"Burst",
"Min Burst"))
segmentation_demo_plt <- ggplot()+
geom_vline(xintercept = c(shading_annotations$starts, shading_annotations$equal_len, shading_annotations$next_start, shading_annotations$equal_on), color = "gray")+
geom_path(data = epsp_stim, aes(Time, Stim+1.5))+
geom_rect(data = shading_annotations, aes(xmin = starts, xmax = next_start, ymin = 0, ymax = 0.5), alpha = 0.3, fill = RColorBrewer::brewer.pal(3, "Set1")[1], color = RColorBrewer::brewer.pal(3, "Set1")[1])+
geom_rect(data = shading_annotations, aes(xmin = starts, xmax = equal_len, ymin = -0.5, ymax = 0), alpha = 0.3, fill = RColorBrewer::brewer.pal(3, "Set1")[1], color = RColorBrewer::brewer.pal(3, "Set1")[1])+
geom_rect(data = shading_annotations, aes(xmin = starts, xmax = on_end, ymin = -1, ymax = -0.5), alpha = 0.3, fill = RColorBrewer::brewer.pal(3, "Set1")[2], color = RColorBrewer::brewer.pal(3, "Set1")[2])+
geom_rect(data = shading_annotations, aes(xmin = starts, xmax = equal_on, ymin = -1.5, ymax = -1), alpha = 0.3, fill = RColorBrewer::brewer.pal(3, "Set1")[2], color = RColorBrewer::brewer.pal(3, "Set1")[2])+
theme_classic()+
# geom_vline(xintercept = 12.54)+
geom_text(data = annotation_df[1:2, ], aes(x=x, y=y, label = text), hjust = "inward", color = RColorBrewer::brewer.pal(3, "Set1")[1])+
geom_text(data = annotation_df[3:4, ], aes(x=x, y=y, label = text), hjust = "inward", color = RColorBrewer::brewer.pal(3, "Set1")[2])+
ylim(-1.5, 9.5)+
labs(title = "Ways to Segement EPSP Stim")
segmentation_demo_plt
save_plot(here("output", "plots", "segmentation_demo_plt.pdf"),
segmentation_demo_plt, base_height = 7, base_width = 10)
Proof of principle: Contrasting simulaiton with LIF to regression and VIR
M <- readABF_as_matrix2(path =
paste0(abf_dir_path, "190915_0019.abf"),
channels = "all",
relative.time = T)
temp <- as_tibble(M) %>%
janitor::clean_names(case = "upper_camel") %>%
mutate(Sweep = as.factor(Sweep)) %>%
group_by(Sweep) %>%
mutate(MinTime = min(Time, na.rm = T)) %>%
mutate(Time = Time - MinTime) %>%
filter(Sweep == 1) #%>%
# mutate(In7 = rollmean(In7, 17, fill = "extend"))
fm <- lm(In4 ~ In7, data = temp)
fm_newdata <- data.frame(In4 = NA,
In7 = temp$In7,
Time = temp$Time,
sweep = temp$Sweep)
fm_newdata$In4 <- predict(fm, newdata = fm_newdata)
temp <- temp[, c("Time", "In7", "In4")]
temp <- full_join(temp, rename(fm_newdata, lm = In4))
#TODO add in Rin / Tau programmatically
place_params_for_predict_voltage_trace(tau = 21.417860031127930, Rin = 8.535)
place_data_for_predict_voltage_trace(input.df = temp, Time.ch = "Time", Inj.ch = "In7")
run_predict_voltage_trace_py()
predicted_response <- retrieve_predicted_voltage()
t1 <- temp
t2 <- as_tibble(predicted_response)
t3 <- cbind(t1, t2[, "predicted"])
t3$Vinf <- (t3$In7*8.535)-49.5
t3$predicted <- t3$predicted-49.5
plt_left <-
t3 %>% gather(key = "Method", value = "mV", c("lm", "Vinf", "predicted")) %>%
ggplot(aes(Time, mV, color = Method))+
geom_path()+
geom_path(data = t3, aes(Time, In4), color = "black")+
theme_minimal()+
theme(legend.position = "")+
facet_grid(Method~.)
t3_subtracted <- t3
t3_subtracted$lm <- t3_subtracted$In4 - t3_subtracted$lm
t3_subtracted$Vinf <- t3_subtracted$In4 - t3_subtracted$Vinf
t3_subtracted$predicted <- t3_subtracted$In4 - t3_subtracted$predicted
t3_subtracted$In4 <- t3_subtracted$In4 - -49.5
plt_mid <-
t3_subtracted %>% gather(key = "Method", value = "mV", c("lm", "Vinf", "predicted")) %>%
ggplot(aes(Time, mV, color = Method))+
geom_path()+ #TODO It seems there is a nA that's causing early termination here?
geom_path(data = t3_subtracted, aes(Time, In4), color = "black")+
theme_minimal()+
theme(legend.position = "")+
facet_grid(Method~.)
plt_right <-
t3 %>% gather(key = "Method", value = "mV", c("lm", "Vinf", "predicted")) %>%
ggplot(aes(In4, mV, color = Method))+
geom_path()+
geom_abline(intercept = 0, slope = 1)+
theme_minimal()+
theme(legend.position = "")+
facet_grid(Method~.)
epsp_prediction_demo <- (plt_left | plt_mid | plt_right) + plot_layout(widths = c(3, 3, 2))
epsp_prediction_demo
save_plot(here("output", "plots", "epsp_prediction_demo.pdf"),
epsp_prediction_demo, base_height = 7, base_width = 10)
predict_epsp_responses <- function(trace = trace_list$EPSP,
tau = 19.8, #abs(mean(unlist(out$df[, "In9_Tau_est"]))),
rin = 20.0, #abs(mean(unlist(out$df[, "In9_R1"]))),
v_ch = "In9",
i_ch = "In12"
){
temp <- trace
place_params_for_predict_voltage_trace(tau = tau,
Rin = rin)
predicted_response_list <- map(seq(1, length(unique(temp$Sweep))), function(i){
place_data_for_predict_voltage_trace(input.df = temp[temp$Sweep == i, ],
Time.ch = "Time",
Inj.ch = i_ch)
run_predict_voltage_trace_py()
predicted_response <- retrieve_predicted_voltage()
pred_v <- predicted_response[predicted_response$Time >= 19.38, "predicted"]
meas_v <- temp[temp$Sweep == i & temp$Time >= 19.38, v_ch]
offset <- median(unlist(pred_v - meas_v), na.rm = T)
predicted_response$offset <- offset
return(predicted_response)
})
my_plts <- map(
seq_along(predicted_response_list),
function(i){
ggplot()+
geom_line(data = predicted_response_list[[i]], aes(Time, predicted-offset), color = "cornflowerblue")+
geom_line(data = temp[temp$Sweep == i, ], aes_string("Time", v_ch), color = "black")+
labs(x = "Seconds", y = "mV")+
coord_cartesian(y = c(-100, +100))
})
trace_plt <- cowplot::plot_grid(plotlist = my_plts)
return(list(diagnostic_plots = list(trace = trace_plt),
predictions = predicted_response_list))
}
out <- predict_epsp_responses()
out$diagnostic_plots$trace
out$predictions
# saveRDS(out$predictions,
# here("data", "predicted_voltage", paste(str_remove(input_file, ".abf"), "epsp", "predict_v.rds", sep = "-")))
# save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "epsp", "segment.pdf", sep = "-")),
# out$diagnostic_plots$trace, base_height = 7, base_width = 10)
danielkick/mRNA24hLC documentation built on Feb. 6, 2024, 2:18 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE) library(readxl) library(here) library(tidyverse) library(janitor) library(zoo) library(readABF) library(broom) library(minpack.lm) # getting Tau library(patchwork) library(ggsci) library(cowplot) library(patchwork) library(minpack.lm) # for use in fit_tau_1term_exp() theme_set(ggplot2::theme_minimal()) library(devtools) devtools::load_all()
abf_dir_path <- "D:/0000_Traces_mRNA24hLC/" abf_dir_path <- "./inst/extdata/ABFs/"
Handling different kinds of files (Proof of principle)
Load in data (ionic currents)
file_list <- list(FI = "190903_0025.abf", EPSP = "190903_0023.abf", GJCC = "190903_0021.abf", GJVC = "190903_0020.abf", HTK = "190903_0031.abf", HTKC = "190903_0034.abf", A = "190903_0030.abf", AC = "190903_0033.abf") trace_list <- map(seq_along(file_list), function(i){ temp <- readABF_as_matrix2(path = paste0(abf_dir_path, # "C:/Users/Daniel/Documents/Trace_Holding/", file_list[[i]]), channels = "all", relative.time = T) %>% as.data.frame() %>% janitor::clean_names(case = "upper_camel") %>% mutate(Sweep = as.factor(Sweep)) %>% group_by(Sweep) %>% mutate(MinTime = min(Time, na.rm = T)) %>% mutate(Time = Time - MinTime) %>% select(-MinTime) }) names(trace_list) <- names(file_list)
Show all
example_plots <- map(trace_list, function(i){ downsample_data(i, len = 5000) %>% ungroup() %>% gather(Ch, Value, c("In4","In7","In9","In12")) %>% mutate(Cell = case_when(Ch %in% c("In4","In7") ~ "Cell A", Ch %in% c("In9","In12") ~ "Cell B",)) %>% mutate(ChType = case_when(Ch %in% c("In4","In9") ~ "mV", Ch %in% c("In7","In12") ~ "nA",)) %>% mutate(Ch = factor(.$Ch, levels = c("In4","In7","In9","In12"))) %>% ggplot(aes(Time, Value, group = Sweep, color = Cell))+ geom_line()+ facet_wrap(ChType~Cell, scales = "free_y", ncol = 2)+ labs(y = "")+ theme(legend.position = "")+ scale_color_aaas() }) plot_grid(plotlist = example_plots)
GJVC
## GJVC: get Ig, Ileak==== process_gjvc_trace <- function(trace = trace_list$GJVC){ ### Setup i_leak <- trace %>% # filter(Time > 0.75) %>% mutate(In4 = as.numeric(In4), In7 = as.numeric(In7), In9 = as.numeric(In9), In12 = as.numeric(In12)) i_leak <- i_leak %>% #When are the step start times? # 1. In9: ~0.27575 # 2. In4: ~1.02532 mutate(Inj_in = ifelse(Time > 0.8, "In4", "In9"), Segment = case_when( Time >= (0.07575 - 0.02532) & Time <= (0.27575 - 0.02532) ~ "PreStep", # Inj In9 Time >= (0.27575+ 0.17468) & Time <= (0.27575+ 0.22468) ~ "Step", Time >= 0.8 & Time <= 1.0 ~ "PreStep", # Inj In4 Time >= 1.2 & Time <= 1.25 ~ "Step" )) plt_tagged_regions <- downsample_data(i_leak, len = 10000) %>% gather("key", "value", c("In4", "In7", "In9", "In12")) %>% mutate(Cell = case_when(key %in% c("In4", "In7") ~ "Cell1", key %in% c("In9", "In12") ~ "Cell2"), Aspect = case_when(key %in% c("In4", "In9") ~ "mV", key %in% c("In7", "In12") ~ "nA") ) %>% ggplot(aes(Time, value, color = Segment, group = Sweep))+ geom_path()+ # facet_wrap(key~., scales = "free") facet_grid(Aspect~Cell, scales = "free") i_leak_for_In4_to_In9 <- i_leak %>% filter(Inj_in == "In4") %>% group_by(Sweep, Segment, Inj_in) %>% mutate(In4Mean = mean(In4, na.rm = T), In7Mean = mean(In7, na.rm = T), In9Mean = mean(In9, na.rm = T), In12Mean =mean(In12, na.rm = T)) %>% ungroup() i_leak_for_In9_to_In4 <- i_leak %>% filter(Inj_in == "In9") %>% group_by(Sweep, Segment, Inj_in) %>% mutate(In4Mean = mean(In4, na.rm = T), In7Mean = mean(In7, na.rm = T), In9Mean = mean(In9, na.rm = T), In12Mean =mean(In12, na.rm = T)) %>% ungroup() ### Ig #### i_leak_for_In4_to_In9 <- i_leak_for_In4_to_In9 %>% select(Sweep, Segment, In4Mean, In7Mean, In9Mean, In12Mean) %>% filter(Segment != "NA") %>% distinct() %>% pivot_wider(names_from = "Segment", values_from = c("In4Mean", "In12Mean", "In9Mean", "In7Mean" )) %>% mutate(In4Mean_Step = In4Mean_Step - In4Mean_PreStep, In12Mean_Step = In12Mean_Step - In12Mean_PreStep, In9Mean_Step = In9Mean_Step - In9Mean_PreStep, In7Mean_Step = In7Mean_Step - In7Mean_PreStep ) %>% select(Sweep, In4Mean_Step, In12Mean_Step, In9Mean_Step, In7Mean_Step) %>% rename(In4Mean = In4Mean_Step, In12Mean = In12Mean_Step, In9Mean = In9Mean_Step, In7Mean = In7Mean_Step) %>% filter(abs(In4Mean) >= 4 | abs(In9Mean) >= 4 )# exclude no step condition i_leak_for_In9_to_In4 <- i_leak_for_In9_to_In4 %>% select(Sweep, Segment, In4Mean, In7Mean, In9Mean, In12Mean) %>% filter(Segment != "NA") %>% distinct() %>% pivot_wider(names_from = "Segment", values_from = c("In4Mean", "In12Mean", "In9Mean", "In7Mean" )) %>% mutate(In4Mean_Step = In4Mean_Step - In4Mean_PreStep, In12Mean_Step = In12Mean_Step - In12Mean_PreStep, In9Mean_Step = In9Mean_Step - In9Mean_PreStep, In7Mean_Step = In7Mean_Step - In7Mean_PreStep ) %>% select(Sweep, In4Mean_Step, In12Mean_Step, In9Mean_Step, In7Mean_Step) %>% rename(In4Mean = In4Mean_Step, In12Mean = In12Mean_Step, In9Mean = In9Mean_Step, In7Mean = In7Mean_Step) %>% filter(abs(In4Mean) >= 4 | abs(In9Mean) >= 4 )# exclude no step condition plt_ig_slopes <- ggplot()+ geom_smooth(data = i_leak_for_In4_to_In9, aes(In4Mean, In12Mean), method = lm, color = "red")+ geom_point( data = i_leak_for_In4_to_In9, aes(In4Mean, In12Mean), color = "red")+ geom_smooth(data = i_leak_for_In9_to_In4, aes(In9Mean, In7Mean), method = lm)+ geom_point( data = i_leak_for_In9_to_In4, aes(In9Mean, In7Mean))+ labs(x = "Difference in mV", y = "nA") In4_to_In9_Ig <- as.numeric(broom::tidy(lm(In12Mean ~ In4Mean, data = i_leak_for_In4_to_In9))[2, "estimate"]) In9_to_In4_Ig <- as.numeric(broom::tidy(lm(In7Mean ~ In9Mean, data = i_leak_for_In9_to_In4))[2, "estimate"]) ### Ileak #### In4_R <- as.numeric(broom::tidy(lm(In7Mean ~ In4Mean, data = i_leak_for_In4_to_In9))[2, "estimate"]) In9_R <- as.numeric(broom::tidy(lm(In12Mean ~ In9Mean, data = i_leak_for_In9_to_In4))[2, "estimate"]) In4_R <- In4_R^-1 In9_R <- In9_R^-1 return( list(diagnostic_plots = list( trace = plt_tagged_regions, fit = plt_ig_slopes), df = data.frame( In4_to_In9 = In4_to_In9_Ig, In9_to_In4 = In9_to_In4_Ig, In4_R1c = In4_R, In9_R1c = In9_R ) ) ) } out <- process_gjvc_trace(trace = trace_list$GJVC) out$diagnostic_plots$trace out$diagnostic_plots$fit out$df # # Save for review: # input_file <- "190903_0020.abf" # # write.csv(out$df, # here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjvc", "df.csv", sep = "-")) ) # # save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjvc", "trace.pdf", sep = "-")), # out$diagnostic_plots$trace, base_height = 7, base_width = 10) # # save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjvc", "fit.pdf", sep = "-")), # out$diagnostic_plots$fit, base_height = 7, base_width = 10)
GJCC
For reference, we use the following equations to determine the resistances of interest. $$r_{11}=\frac{v_1}{i_1}$$ $$r_{12}=\frac{v_2}{i_1}$$ $$r_1= \frac{r_{11}*r_{22} - r_{12}^2}{r_{22} - r_{12}}$$
$$r_c = \frac{r_{11}*r_{22} - r_{12}^2}{r_{12}}$$
## GJCC: get CC, R_C, C, Rin ==== process_gjcc_trace <- function(trace = trace_list$GJCC){ ## prep ==== temp <- trace temp_plt <- temp %>% # filter(Sweep == 4) %>% mutate(In4 = as.numeric(In4)) %>% mutate(In7 = as.numeric(In7)) # Show what the protocol looks like plt1 <- downsample_data(temp_plt, len = 1000) %>% gather(key, value, c("In4", "In9")) %>% ggplot(aes(Time, value, color = key, group = interaction(Sweep, key)))+ geom_path()+ scale_color_manual(values = RColorBrewer::brewer.pal(7, name = "PuOr")[c(1, 7)] ) plt2 <- downsample_data(temp_plt, len = 1000) %>% gather(key, value, c("In7", "In12")) %>% # mutate(key = factor(.$key, levels = c("In7", "In12"))) %>% ggplot(aes(Time, value, color = key, group = interaction(Sweep, key)))+ geom_path()+ scale_color_manual(values = RColorBrewer::brewer.pal(7, name = "PuOr")[c(7, 1)] )# Without flipping the factor level, get the right colors. plt_trace <- plt1/plt2 # annotate temp <- temp %>% mutate(Segment = case_when(Time >= 0 & Time <= 0.32820 ~ "A", #"PreStepIn12", Time >= 1.32836 & Time <= 1.82836 ~ "B", #"StepIn12", Time >= 2.57744 & Time <= 3.07744 ~ "C", #"PreStepIn7", Time >= 4.07760 & Time <= 4.57760 ~ "D", #"StepIn7" )) ## resistances/ coupling ==== # $$r_{11}=\frac{v_1}{i_1}$$ # $$r_{12}=\frac{v_2}{i_1}$$ # $$r_1= \frac{r_{11}*r_{22} - r_{12}^2}{r_{22} - r_{12}}$$ # $$r_c = \frac{r_{11}*r_{22} - r_{12}^2}{r_{12}}$$ temp_resist <- temp %>% group_by(Sweep, Segment) %>% filter(!is.na(Segment)) %>% summarise(MeanIn4 = mean(In4, na.rm = T), MeanIn7 = mean(In7, na.rm = T), MeanIn9 = mean(In9, na.rm = T), MeanIn12 = mean(In12, na.rm = T)) %>% ungroup() %>% pivot_wider(names_from = Segment, values_from = c("MeanIn4", "MeanIn7", "MeanIn9", "MeanIn12")) %>% mutate(In4_CC = MeanIn9_D/MeanIn4_D, #inj in In4 In4_R11 = MeanIn4_D/MeanIn7_D, In4_R12 = MeanIn9_D/MeanIn7_D, In9_CC = MeanIn4_B/MeanIn9_B, #inj in In9 In9_R11 = MeanIn9_B/MeanIn12_B, In9_R12 = MeanIn4_B/MeanIn12_B ) %>% mutate(mean_R12 = (In4_R12+In9_R12)/2) %>% mutate(In4_R1 = ((In4_R11*In9_R11) - (mean_R12^2)) / (In9_R11-mean_R12), In4_Rc = ((In4_R11*In9_R11) - (mean_R12^2)) / (mean_R12), In9_R1 = ((In9_R11*In4_R11) - (mean_R12^2)) / (In4_R11-mean_R12), In9_Rc = ((In9_R11*In4_R11) - (mean_R12^2)) / (mean_R12) ) %>% select(starts_with("In"), "Sweep") ## V_rest ==== temp_rmp <- temp %>% filter(Segment == "C") %>% group_by(Sweep) %>% summarise(In4_Vrest = min(In4, na.rm = T), In9_Vrest = min(In9, na.rm = T)) ## Tau ==== temp <- temp %>% mutate(Segment = case_when(Time >= 3.07744 & Time <= 4.47760 ~ "In4Step", #"PreStepIn7", Time >= 0.32820 & Time <= 1.72836 ~ "In9Step", #"StepIn7" )) plt1 <- downsample_data(temp, len = 1000) %>% gather(Key, Value, c("In4", "In9")) %>% ggplot(aes(Time, Value, color = Segment, group = interaction(Key, Sweep)))+ geom_path() plt2 <- downsample_data(temp, len = 1000) %>% ungroup() %>% gather(key, value, c("In7", "In12")) %>% mutate(key = factor(.$key, levels = c("In7", "In12"))) %>% ggplot(aes(Time, value, color = Segment, group = interaction(key, Sweep)))+ geom_path() plt_segments <- plt1/plt2 ### in4 temp_in4 <- temp %>% filter(Segment == "In4Step") %>% mutate(Time = Time - min(Time, na.rm = T)) test_fits <- map(unique(temp_in4$Sweep), function(i){ fit_tau_1term_exp(df = filter(temp_in4, Sweep == i), IV = "Time", DV = "In4") }) fits_in4 <- map(seq_along(transpose(test_fits)[[1]]), function(i){ df <- transpose(test_fits)[[1]][i][[1]] df <- df[df$term == "b", c("estimate", "std.error")] df$Sweep <- i return(df) }) %>% do.call(rbind, .) %>% rename(In4_Tau_est = estimate, In4_Tau_std.err = std.error) %>% mutate(Sweep = as.factor(Sweep)) test_fits_plts <- map(seq_along(test_fits), function(i){ ggplot()+ geom_line(data = as.data.frame(test_fits[[i]]$check_fit), aes(x = Time, y = Fit), color = "firebrick")+ geom_line(data = as.data.frame(test_fits[[i]]$check_fit), aes(x = Time, y = In4), alpha = 0.8)+ labs(subtitle = paste("tau =", as.character( round(test_fits[[i]]$fit[2, 2], digits = 3) ))) }) plt_tau_in4 <- cowplot::plot_grid(plotlist = test_fits_plts) ### in9 temp_in9 <- temp %>% filter(Segment == "In9Step") %>% mutate(Time = Time - min(Time, na.rm = T)) test_fits <- map(unique(temp_in9$Sweep), function(i){ fit_tau_1term_exp(df = filter(temp_in9, Sweep == i), IV = "Time", DV = "In9") }) fits_in9 <- map(seq_along(transpose(test_fits)[[1]]), function(i){ df <- transpose(test_fits)[[1]][i][[1]] df <- df[df$term == "b", c("estimate", "std.error")] df$Sweep <- i return(df) }) %>% do.call(rbind, .) %>% rename(In9_Tau_est = estimate, In9_Tau_std.err = std.error) %>% mutate(Sweep = as.factor(Sweep)) test_fits_plts <- map(seq_along(test_fits), function(i){ ggplot()+ geom_line(data = as.data.frame(test_fits[[i]]$check_fit), aes(x = Time, y = Fit), color = "firebrick")+ geom_line(data = as.data.frame(test_fits[[i]]$check_fit), aes(x = Time, y = In9), alpha = 0.8)+ labs(subtitle = paste("tau =", as.character( round(test_fits[[i]]$fit[2, 2], digits = 3) ))) }) plt_tau_in9 <- cowplot::plot_grid(plotlist = test_fits_plts) ## Condense ==== summary_df <- full_join(temp_resist, temp_rmp) %>% full_join(., fits_in4) %>% full_join(., fits_in9) return(list( diagnostic_plots = list( trace = plt_trace, segment = plt_segments, fit_In4 = plt_tau_in4, fit_In9 = plt_tau_in9 ), df = summary_df )) } out <- process_gjcc_trace(trace = trace_list$GJCC) out$diagnostic_plots$trace out$diagnostic_plots$segment out$diagnostic_plots$fit_In4 out$diagnostic_plots$fit_In9 out$df # # Save for review: # input_file <- "190903_0021.abf" # # write.csv(out$df, # here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjcc", "df.csv", sep = "-"))) # # save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjcc", "trace.pdf", sep = "-")), # out$diagnostic_plots$trace, base_height = 7, base_width = 10) # # save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjcc", "segment.pdf", sep = "-")), # out$diagnostic_plots$segment, base_height = 7, base_width = 10) # # save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjcc", "fit_in4.pdf", sep = "-")), # out$diagnostic_plots$fit_In4, base_height = 7, base_width = 10) # # save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "gjcc", "fit_in9.pdf", sep = "-")), # out$diagnostic_plots$fit_In9, base_height = 7, base_width = 10)
FI
temp <- trace_list$FI downsample_data(df = temp, len = 5000) %>% gather(key, value, c("In9", "In12")) %>% ungroup() %>% mutate(key = factor(.$key, levels = c("In9", "In12"))) %>% ggplot(aes(Time, value, group = Sweep))+ geom_line()+ facet_wrap(key~., scales = "free_y", ncol = 1) #r11 == ~15.3 #r1 == ~17.9 place_params_for_predict_voltage_trace(tau = 19.9, Rin = 15.3) # place_data_for_predict_voltage_trace(input.df = temp[temp$Sweep == 10, ], # Time.ch = "Time", # Inj.ch = "In12") # # run_predict_voltage_trace_py() # # predicted_response <- retrieve_predicted_voltage() predicted_response_list <- map(seq(1, length(unique(temp$Sweep))), function(i){ place_data_for_predict_voltage_trace(input.df = temp[temp$Sweep == i, ], Time.ch = "Time", Inj.ch = "In12") run_predict_voltage_trace_py() predicted_response <- retrieve_predicted_voltage() return(predicted_response) }) my_plts <- map( seq_along(predicted_response_list), function(i){ ggplot()+ geom_line(data = predicted_response_list[[i]], aes(Time, predicted-56.5))+ geom_line(data = temp[temp$Sweep == i, ], aes(Time, In9), color = "red")+ labs(x = "Seconds", y = "mV")+ coord_cartesian(y = c(-150, +120)) }) # cowplot::plot_grid(plotlist = my_plts) multi_plt <- my_plts[[11]] / my_plts[[11 - 2]] / my_plts[[1 + 2]] # # Save for review: # input_file <- "190903_0025.abf" # saveRDS(predicted_response_list, # here("data", "predicted_voltage", paste(str_remove(input_file, ".abf"), "fi", "predict_v.rds", sep = "-"))) # # save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "fi", "multi_plt.pdf", sep = "-")), # multi_plt, base_height = 7, base_width = 10) # test_load <- readRDS(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "fi", "predict_v.rds", sep = "-")))
predict_fi_responses <- function(trace = trace_list$FI, tau = 19.8, #abs(mean(unlist(out$df[, "In9_Tau_est"]))), rin = 20.0, #abs(mean(unlist(out$df[, "In9_R1"]))), v_ch = "In9", i_ch = "In12"){ # trace = trace_list$FI # tau = 19.8 # rin = 20.0 # v_ch = "In9" # i_ch = "In12" temp <- trace # downsample_data(df = temp, len = 10000) %>% # gather(key, value, c(v_ch, i_ch)) %>% # ungroup() %>% # mutate(key = factor(.$key, levels = c(v_ch, i_ch))) %>% # ggplot(aes(Time, value, group = Sweep))+ # geom_line()+ # facet_wrap(key~., scales = "free_y", ncol = 1) place_params_for_predict_voltage_trace(tau = tau, Rin = rin) predicted_response_list <- map(seq(1, length(unique(temp$Sweep))), function(i){ place_data_for_predict_voltage_trace(input.df = temp[temp$Sweep == i, ], Time.ch = "Time", Inj.ch = i_ch) run_predict_voltage_trace_py() predicted_response <- retrieve_predicted_voltage() pred_v <- predicted_response[predicted_response$Time >= 4, "predicted"] meas_v <- temp[temp$Sweep == i & temp$Time >= 4, v_ch] offset <- median(unlist(pred_v - meas_v), na.rm = T) predicted_response$offset <- offset return(predicted_response) }) my_plts <- map( seq_along(predicted_response_list), function(i){ ggplot()+ geom_line(data = predicted_response_list[[i]], aes(Time, predicted-offset), color = "red")+ geom_line(data = temp[temp$Sweep == i, ], aes(Time, In9), color = "black")+ labs(x = "Seconds", y = "mV")+ coord_cartesian(y = c(-150, +120)) }) return(list(diagnostic_plots = list(traces = my_plts), predictions = predicted_response_list)) } # my_plts[[8]] # # cowplot::plot_grid(plotlist = my_plts) # # multi_plt <- my_plts[[11]] / my_plts[[11 - 2]] / my_plts[[1 + 2]]
EPSP
(epsp) How could we slice up the trace?
## Here's the stimulus used in the protocol ==== # sweep duration should be 19.687 seconds epsp_stim <- readABF_as_matrix(here("inst", "extdata", "epsp_stim_files", "08 current injection.abf"), channels = "Axo1I2") epsp_stim <- as_tibble(epsp_stim) %>% mutate(Time = Time - min(Time, na.rm = T)) %>% rename(Stim = Axo1I2) shading_annotations <- data.frame( starts = c(0.40, 4.87, 10.74, 16.74), next_start = c(4.87, 10.74, 16.74, 19.685), equal_len = c(0.40, 4.87, 10.74, 16.74) + 2.945, on_end = c(2, 6.17, 12.54, 18.5), equal_on = c(0.40, 4.87, 10.74, 16.74) + 1.30 ) annotation_df <- data.frame(x = c(0.5, 0.5, 2.1, 1.8), y = seq(0.25, -1.25, length.out = 4), text = c("Period", "Min. Period", "Burst", "Min Burst")) segmentation_demo_plt <- ggplot()+ geom_vline(xintercept = c(shading_annotations$starts, shading_annotations$equal_len, shading_annotations$next_start, shading_annotations$equal_on), color = "gray")+ geom_path(data = epsp_stim, aes(Time, Stim+1.5))+ geom_rect(data = shading_annotations, aes(xmin = starts, xmax = next_start, ymin = 0, ymax = 0.5), alpha = 0.3, fill = RColorBrewer::brewer.pal(3, "Set1")[1], color = RColorBrewer::brewer.pal(3, "Set1")[1])+ geom_rect(data = shading_annotations, aes(xmin = starts, xmax = equal_len, ymin = -0.5, ymax = 0), alpha = 0.3, fill = RColorBrewer::brewer.pal(3, "Set1")[1], color = RColorBrewer::brewer.pal(3, "Set1")[1])+ geom_rect(data = shading_annotations, aes(xmin = starts, xmax = on_end, ymin = -1, ymax = -0.5), alpha = 0.3, fill = RColorBrewer::brewer.pal(3, "Set1")[2], color = RColorBrewer::brewer.pal(3, "Set1")[2])+ geom_rect(data = shading_annotations, aes(xmin = starts, xmax = equal_on, ymin = -1.5, ymax = -1), alpha = 0.3, fill = RColorBrewer::brewer.pal(3, "Set1")[2], color = RColorBrewer::brewer.pal(3, "Set1")[2])+ theme_classic()+ # geom_vline(xintercept = 12.54)+ geom_text(data = annotation_df[1:2, ], aes(x=x, y=y, label = text), hjust = "inward", color = RColorBrewer::brewer.pal(3, "Set1")[1])+ geom_text(data = annotation_df[3:4, ], aes(x=x, y=y, label = text), hjust = "inward", color = RColorBrewer::brewer.pal(3, "Set1")[2])+ ylim(-1.5, 9.5)+ labs(title = "Ways to Segement EPSP Stim") segmentation_demo_plt save_plot(here("output", "plots", "segmentation_demo_plt.pdf"), segmentation_demo_plt, base_height = 7, base_width = 10)
Proof of principle: Contrasting simulaiton with LIF to regression and VIR
M <- readABF_as_matrix2(path = paste0(abf_dir_path, "190915_0019.abf"), channels = "all", relative.time = T) temp <- as_tibble(M) %>% janitor::clean_names(case = "upper_camel") %>% mutate(Sweep = as.factor(Sweep)) %>% group_by(Sweep) %>% mutate(MinTime = min(Time, na.rm = T)) %>% mutate(Time = Time - MinTime) %>% filter(Sweep == 1) #%>% # mutate(In7 = rollmean(In7, 17, fill = "extend")) fm <- lm(In4 ~ In7, data = temp) fm_newdata <- data.frame(In4 = NA, In7 = temp$In7, Time = temp$Time, sweep = temp$Sweep) fm_newdata$In4 <- predict(fm, newdata = fm_newdata) temp <- temp[, c("Time", "In7", "In4")] temp <- full_join(temp, rename(fm_newdata, lm = In4)) #TODO add in Rin / Tau programmatically place_params_for_predict_voltage_trace(tau = 21.417860031127930, Rin = 8.535) place_data_for_predict_voltage_trace(input.df = temp, Time.ch = "Time", Inj.ch = "In7") run_predict_voltage_trace_py() predicted_response <- retrieve_predicted_voltage() t1 <- temp t2 <- as_tibble(predicted_response) t3 <- cbind(t1, t2[, "predicted"]) t3$Vinf <- (t3$In7*8.535)-49.5 t3$predicted <- t3$predicted-49.5 plt_left <- t3 %>% gather(key = "Method", value = "mV", c("lm", "Vinf", "predicted")) %>% ggplot(aes(Time, mV, color = Method))+ geom_path()+ geom_path(data = t3, aes(Time, In4), color = "black")+ theme_minimal()+ theme(legend.position = "")+ facet_grid(Method~.) t3_subtracted <- t3 t3_subtracted$lm <- t3_subtracted$In4 - t3_subtracted$lm t3_subtracted$Vinf <- t3_subtracted$In4 - t3_subtracted$Vinf t3_subtracted$predicted <- t3_subtracted$In4 - t3_subtracted$predicted t3_subtracted$In4 <- t3_subtracted$In4 - -49.5 plt_mid <- t3_subtracted %>% gather(key = "Method", value = "mV", c("lm", "Vinf", "predicted")) %>% ggplot(aes(Time, mV, color = Method))+ geom_path()+ #TODO It seems there is a nA that's causing early termination here? geom_path(data = t3_subtracted, aes(Time, In4), color = "black")+ theme_minimal()+ theme(legend.position = "")+ facet_grid(Method~.) plt_right <- t3 %>% gather(key = "Method", value = "mV", c("lm", "Vinf", "predicted")) %>% ggplot(aes(In4, mV, color = Method))+ geom_path()+ geom_abline(intercept = 0, slope = 1)+ theme_minimal()+ theme(legend.position = "")+ facet_grid(Method~.) epsp_prediction_demo <- (plt_left | plt_mid | plt_right) + plot_layout(widths = c(3, 3, 2)) epsp_prediction_demo save_plot(here("output", "plots", "epsp_prediction_demo.pdf"), epsp_prediction_demo, base_height = 7, base_width = 10)
predict_epsp_responses <- function(trace = trace_list$EPSP, tau = 19.8, #abs(mean(unlist(out$df[, "In9_Tau_est"]))), rin = 20.0, #abs(mean(unlist(out$df[, "In9_R1"]))), v_ch = "In9", i_ch = "In12" ){ temp <- trace place_params_for_predict_voltage_trace(tau = tau, Rin = rin) predicted_response_list <- map(seq(1, length(unique(temp$Sweep))), function(i){ place_data_for_predict_voltage_trace(input.df = temp[temp$Sweep == i, ], Time.ch = "Time", Inj.ch = i_ch) run_predict_voltage_trace_py() predicted_response <- retrieve_predicted_voltage() pred_v <- predicted_response[predicted_response$Time >= 19.38, "predicted"] meas_v <- temp[temp$Sweep == i & temp$Time >= 19.38, v_ch] offset <- median(unlist(pred_v - meas_v), na.rm = T) predicted_response$offset <- offset return(predicted_response) }) my_plts <- map( seq_along(predicted_response_list), function(i){ ggplot()+ geom_line(data = predicted_response_list[[i]], aes(Time, predicted-offset), color = "cornflowerblue")+ geom_line(data = temp[temp$Sweep == i, ], aes_string("Time", v_ch), color = "black")+ labs(x = "Seconds", y = "mV")+ coord_cartesian(y = c(-100, +100)) }) trace_plt <- cowplot::plot_grid(plotlist = my_plts) return(list(diagnostic_plots = list(trace = trace_plt), predictions = predicted_response_list)) } out <- predict_epsp_responses() out$diagnostic_plots$trace out$predictions # saveRDS(out$predictions, # here("data", "predicted_voltage", paste(str_remove(input_file, ".abf"), "epsp", "predict_v.rds", sep = "-"))) # save_plot(here("data", "data_audits", paste(str_remove(input_file, ".abf"), "epsp", "segment.pdf", sep = "-")), # out$diagnostic_plots$trace, base_height = 7, base_width = 10)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.