knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library("readxl")
library("tidyverse")
library("cowplot")
library("car")
library("SchulzLabTools")

if (Sys.info()[1] == "Darwin") {
  # on macos
  use.path <- "/Volumes/schulzlab/Data_Daniel/"
} else if (Sys.info()[1] == c("Windows")) {
  # on windows
  use.path <- "S:/Data_Daniel/"
} else {
  print(c("where should I look?"))
}

Initial data wrangling

Convenience functions:

Read in phase change data and manipulate it for merger.

#handle data with the old format (non data/metadata organized)
#i.e. 180207_phase_shift_reorganized
iclamp <- as.data.frame(read.csv(paste0(use.path, "180207_phase_shift_reorganized/data/iclamp2.csv")))
vclamp <- as.data.frame(read.csv(paste0(use.path, "180207_phase_shift_reorganized/data/vclamp2.csv")))
current <- as.data.frame(read.csv(paste0(use.path, "180207_phase_shift_reorganized/data/current2.csv")))
#1.  transform iclamp into something that can be merged with tecc
iclamp <- iclamp[,c("experiment", "file", "trace", "trace_start",
                    "IN4_mean", "IN6_mean", "IN7_mean", "IN9_mean", "IN11_mean", "IN12_mean", 
                    "IN4_base", "IN6_base", "IN7_base", "IN9_base", "IN11_base", "IN12_base", 
                    "filepath", "time", "phase_offset", "IN4_id", "IN9_id")]

#alter file col so that it has just the suffix of the filename (tevc$Recording)
iclamp$file <- gsub(pattern = ".+_", replacement = "", iclamp$file)
iclamp$file <- gsub(pattern = "\\.abf", replacement = "", iclamp$file)

#rename everything in `phase_offset` so that it can be converted to `condition`
iclamp[iclamp$phase_offset == 0, "phase_offset"] <- "Phase Offset 0 Degrees"
iclamp[iclamp$phase_offset == 45, "phase_offset"] <- "Phase Offset 45 Degrees"
iclamp[iclamp$phase_offset == 90, "phase_offset"] <- "Phase Offset 90 Degrees"
iclamp[iclamp$phase_offset == 180, "phase_offset"] <- "Phase Offset 180 Degrees"

#rename cols to match tevc
iclamp <- iclamp %>% 
  rename(Experiment = experiment) %>% rename(Recording = file) %>% 
  rename(Trace = trace) %>% rename(`Trace Start` = trace_start) %>% 
  rename(IN4_baseline = IN4_base) %>% rename(IN6_baseline = IN6_base) %>% 
  rename(IN7_baseline = IN7_base) %>% 
  rename(IN9_baseline = IN9_base) %>% rename(IN11_baseline = IN11_base) %>% 
  rename(IN12_baseline = IN12_base) %>% 
  rename(path = filepath) %>% rename(`Time Exposed` = time) %>% 
  rename(Condition = phase_offset) %>% 
  rename(IN4 = IN4_id) %>% rename(IN9 = IN9_id)

#add a Treatment col
iclamp$Treatment <- "Phase Shift"
#ready to rbind into dataset

# 2. transform vclamp into something that can be merged with tevc
#rearrange to get cols close to target str
vclamp <- vclamp[,c("experiment",
                    "file",
                    "trace",
                    "trace_start",
                    "IN4_mean",
                    "IN7_mean",
                    "IN9_mean",
                    "IN12_mean",
                    "IN4_base",
                    "IN7_base",
                    "IN9_base",
                    "IN12_base",
                    "filepath",
                    "time",
                    "phase_offset",
                    "IN4_id",
                    "IN9_id",

                    "injection",
                    "notes",
                    "X")]
#drop cols not in tevc (we'll infer injected cell later)
vclamp <- vclamp[, -which(names(vclamp) %in% c("injection", "notes", "X"))] 

#alter file col so that it has just the suffix of the filename (tevc$Recording)
vclamp$file <- gsub(pattern = ".+_", replacement = "", vclamp$file)
vclamp$file <- gsub(pattern = "\\.abf", replacement = "", vclamp$file)
vclamp$file <- as.character(vclamp$file)

#rename everything in `phase_offset` so that it can be converted to `condition`
vclamp[vclamp$phase_offset == 0, "phase_offset"] <- "Phase Offset 0 Degrees"
vclamp[vclamp$phase_offset == 45, "phase_offset"] <- "Phase Offset 45 Degrees"
vclamp[vclamp$phase_offset == 90, "phase_offset"] <- "Phase Offset 90 Degrees"
vclamp[vclamp$phase_offset == 180, "phase_offset"] <- "Phase Offset 180 Degrees"

#rename cols to match tevc
vclamp <- vclamp %>% 
  rename(Experiment = experiment) %>% rename(Recording = file) %>% 
  rename(Trace = trace) %>% rename(`Trace Start` = trace_start) %>% 
  rename(IN4_baseline = IN4_base) %>% rename(IN7_baseline = IN7_base) %>% 
  rename(IN9_baseline = IN9_base) %>% rename(IN12_baseline = IN12_base) %>% 
  rename(path = filepath) %>% rename(`Time Exposed` = time) %>% 
  rename(Condition = phase_offset) %>% 
  rename(IN4 = IN4_id) %>% rename(IN9 = IN9_id)

#add a Treatment col
vclamp$Treatment <- "Phase Shift"
#ready to rbind into dataset

#3. split current into two df that can be merged with `a` and `htk`
phase.a <- current[,c("experiment",
                      "a.file", 
                      "IN4.a.peak_mV", "IN4.a.peak_nA", 
                      "IN9.a.peak_mV", "IN9.a.peak_nA",
                      "IN4.a.end_mV", "IN4.a.end_nA", 
                      "IN9.a.end_mV", "IN9.a.end_nA", 
                      "time", 
                      "phase_offset", 
                      "IN4_id", 
                      "IN9_id")]

phase.htk <- current[,c("experiment",
                        "htk.file", 
                        "IN4.htk.peak_mV", "IN4.htk.peak_nA", 
                        "IN9.htk.peak_mV", "IN9.htk.peak_nA",
                        "IN4.htk.end_mV",  "IN4.htk.end_nA", 
                        "IN9.htk.end_mV", "IN9.htk.end_nA", 
                        "IN4.rin_used", 
                        "IN9.rin_used", 
                        "time", 
                        "phase_offset", 
                        "IN4_id", 
                        "IN9_id")]
#rename everything in `phase_offset` so that it can be converted to `condition`
phase.a[phase.a$phase_offset == 0, "phase_offset"] <- "Phase Offset 0 Degrees"
phase.a[phase.a$phase_offset == 45, "phase_offset"] <- "Phase Offset 45 Degrees"
phase.a[phase.a$phase_offset == 90, "phase_offset"] <- "Phase Offset 90 Degrees"
phase.a[phase.a$phase_offset == 180, "phase_offset"] <- "Phase Offset 180 Degrees"

phase.htk[phase.htk$phase_offset == 0, "phase_offset"] <- "Phase Offset 0 Degrees"
phase.htk[phase.htk$phase_offset == 45, "phase_offset"] <- "Phase Offset 45 Degrees"
phase.htk[phase.htk$phase_offset == 90, "phase_offset"] <- "Phase Offset 90 Degrees"
phase.htk[phase.htk$phase_offset == 180, "phase_offset"] <- "Phase Offset 180 Degrees"

#rename cols to match a/htk
phase.a <- phase.a %>% 
  rename(Experiment = experiment) %>% rename(Recording = a.file) %>% 
  rename(IN4_mV_peak = IN4.a.peak_mV) %>% rename(IN7_nA_peak = IN4.a.peak_nA) %>% 
  rename(IN9_mV_peak = IN9.a.peak_mV) %>% rename(IN12_nA_peak = IN9.a.peak_nA) %>% 
  rename(IN4_mV_end = IN4.a.end_mV) %>% rename(IN7_nA_end = IN4.a.end_nA) %>% 
  rename(IN9_mV_end = IN9.a.end_mV) %>% rename(IN12_nA_end = IN9.a.end_nA) %>% 
  rename(`Time Exposed` = time) %>% 
  rename(Condition = phase_offset) %>% 
  rename(IN4 = IN4_id) %>% rename(IN9 = IN9_id)

phase.htk <- phase.htk %>% 
  rename(Experiment = experiment) %>% rename(Recording = htk.file) %>% 
  rename(IN4_mV_peak = IN4.htk.peak_mV) %>% rename(IN7_nA_peak = IN4.htk.peak_nA) %>% 
  rename(IN9_mV_peak = IN9.htk.peak_mV) %>% rename(IN12_nA_peak = IN9.htk.peak_nA) %>% 
  rename(IN4_mV_end = IN4.htk.end_mV) %>% rename(IN7_nA_end = IN4.htk.end_nA) %>% 
  rename(IN9_mV_end = IN9.htk.end_mV) %>% rename(IN12_nA_end = IN9.htk.end_nA) %>% 
  rename(IN4_rin = IN4.rin_used) %>% rename(IN9_rin = IN9.rin_used) %>% 
  rename(`Time Exposed` = time) %>% 
  rename(Condition = phase_offset) %>% 
  rename(IN4 = IN4_id) %>% rename(IN9 = IN9_id)

#ensure that Recording is char
phase.a$Recording <- as.character(phase.a$Recording)
phase.htk$Recording <- as.character(phase.htk$Recording)

#add a Treatment col
phase.a$Treatment <- "Phase Shift"
phase.htk$Treatment <- "Phase Shift"
#ready to rbind into datasets

Prepare inverse, aberrant, and controls for merger.

inverted.wave <- condense_data_to_lists(dir.name = "180117_inverted_wave", treatment.name = "Inverted Wave")
ctrl.tea <- condense_data_to_lists(dir.name = "180110_Blocker_washon_hold_to_ctrl", treatment.name = "Control Wave + TEA")
aberrant.wave <- condense_data_to_lists(dir.name = "180205_TEA_wave_in_control", treatment.name = "Aberrant Wave")
ctrl.silent <- condense_data_to_lists(dir.name = "180209_silent_time_control", treatment.name = "Silent Control")
small.phase.angle <- condense_data_to_lists(dir.name = "180308_Small_Phase_Angle", treatment.name = "Small Phase Angle")

Merge by data type

#3. combine into a single df based on data type (e.g. tevc_a)
tecc    <- rbind(iclamp, inverted.wave[[1]], ctrl.tea[[1]], aberrant.wave[[1]], ctrl.silent[[1]], small.phase.angle[[1]])
tevc    <- rbind(vclamp, inverted.wave[[2]], ctrl.tea[[2]], aberrant.wave[[2]], ctrl.silent[[2]], small.phase.angle[[2]])
a       <- rbind(phase.a, inverted.wave[[3]], ctrl.tea[[3]], aberrant.wave[[3]], ctrl.silent[[3]], small.phase.angle[[3]])
htk     <- rbind(phase.htk, inverted.wave[[4]], ctrl.tea[[4]], aberrant.wave[[4]], ctrl.silent[[4]], small.phase.angle[[4]])
htk.ls  <- rbind(inverted.wave[[5]], ctrl.tea[[5]], aberrant.wave[[5]], ctrl.silent[[5]], small.phase.angle[[5]])

Control current subset of data ----------------------------------------------

if (TRUE == FALSE){
  possible.treatments <- c("Phase Shift","Inverted Wave","Control Wave + TEA","Aberrant Wave","Silent Control","Small Phase Angle")

  current.treatment <- possible.treatments[1]

  tecc <- tecc[tecc$Treatment == current.treatment,]
  tevc <- tevc[tevc$Treatment == current.treatment,]
  a <- a[a$Treatment == current.treatment,]
  htk <- htk[htk$Treatment == current.treatment,]
  htk.ls <- htk.ls[htk.ls$Treatment == current.treatment,]
}

Work with current subset of data --------------------------------------------

Cleaning: Drop non-standard timepoints and experiments that don't go to completion.

tecc <- tecc[(tecc$IN4 == "LC4" | tecc$IN4 == "LC5") &
       (tecc$IN4 == "LC4" | tecc$IN4 == "LC5"),]
# Drop timepoints that were only collected in early experiments
tecc <- tecc[tecc$`Time Exposed` != 10 &
               tecc$`Time Exposed` != 30 &
               tecc$`Time Exposed` != 50, ]

tevc <- tevc[tevc$`Time Exposed` != 10 &
               tevc$`Time Exposed` != 30 &
               tevc$`Time Exposed` != 50, ]

a <- a[a$`Time Exposed` != 10 &
               a$`Time Exposed` != 30 &
               a$`Time Exposed` != 50, ]

htk <- htk[htk$`Time Exposed` != 10 &
               htk$`Time Exposed` != 30 &
               htk$`Time Exposed` != 50, ]

tevc <- tecc[!(tecc$Treatment == "Phase Shift" & tecc$`Time Exposed` > 90) &
       !(tecc$Treatment == "Small Phase Angle" & tecc$`Time Exposed` > 160) &
       !(tecc$Treatment == "Inverted Wave" & tecc$`Time Exposed` > 90) &
       !(tecc$Treatment == "Control Wave + TEA" & tecc$`Time Exposed` > 90) &
       !(tecc$Treatment == "Aberrant Wave" & tecc$`Time Exposed` > 90) &
       !(tecc$Treatment == "PhaseShift" & tecc$`Time Exposed` > 90) &
       !(tecc$Treatment == "Silent Control" & tecc$`Time Exposed` > 90),]

tevc <- tevc[!(tevc$Treatment == "Phase Shift" & tevc$`Time Exposed` > 90) &
       !(tevc$Treatment == "Small Phase Angle" & tevc$`Time Exposed` > 160) &
       !(tevc$Treatment == "Inverted Wave" & tevc$`Time Exposed` > 90) &
       !(tevc$Treatment == "Control Wave + TEA" & tevc$`Time Exposed` > 90) &
       !(tevc$Treatment == "Aberrant Wave" & tevc$`Time Exposed` > 90) &
       !(tevc$Treatment == "PhaseShift" & tevc$`Time Exposed` > 90) &
       !(tevc$Treatment == "Silent Control" & tevc$`Time Exposed` > 90),]

a <- a[!(a$Treatment == "Phase Shift" & a$`Time Exposed` > 90) &
       !(a$Treatment == "Small Phase Angle" & a$`Time Exposed` > 160) &
       !(a$Treatment == "Inverted Wave" & a$`Time Exposed` > 90) &
       !(a$Treatment == "Control Wave + TEA" & a$`Time Exposed` > 90) &
       !(a$Treatment == "Aberrant Wave" & a$`Time Exposed` > 90) &
       !(a$Treatment == "PhaseShift" & a$`Time Exposed` > 90) &
       !(a$Treatment == "Silent Control" & a$`Time Exposed` > 90),]

htk <- htk[!(htk$Treatment == "Phase Shift" & htk$`Time Exposed` > 90) &
       !(htk$Treatment == "Small Phase Angle" & htk$`Time Exposed` > 160) &
       !(htk$Treatment == "Inverted Wave" & htk$`Time Exposed` > 90) &
       !(htk$Treatment == "Control Wave + TEA" & htk$`Time Exposed` > 90) &
       !(htk$Treatment == "Aberrant Wave" & htk$`Time Exposed` > 90) &
       !(htk$Treatment == "PhaseShift" & htk$`Time Exposed` > 90) &
       !(htk$Treatment == "Silent Control" & htk$`Time Exposed` > 90),]

Prepare Two electrode current clamp df

Prep work: condense to median values

#Silent Control, Aberrant Wave, and Control Wave + TEA:
#IN7_mean and IN12_mean show POSITIVE values instead of negative. This is a problem in the xlsx which since I copied directly from clampfit means the error is coming from there.
#Inverted Wave
#IN7_mean and IN12_mean include values around POSITIVE and NEGATIVE 2. If something chaned in clampfit, it may have been over the course of working with these files.
#Phase Shift
#Only negative values. No problems here.

#to solve for the above, we'll take all injections with an absolute value above threshfold and set them as negative. Since this is current clamp anything positive should be dropped (from depricated protocols that stepped to positive voltages) or noise anyway

threshold.nA <- -0.9
#fix positive entries
tecc[abs(tecc$IN7_mean) > abs(threshold.nA), "IN7_mean"] <- 
  (abs(tecc[abs(tecc$IN7_mean) > abs(threshold.nA), "IN7_mean"]) * -1)

tecc[abs(tecc$IN12_mean) > abs(threshold.nA), "IN12_mean"] <- 
  (abs(tecc[abs(tecc$IN12_mean) > abs(threshold.nA), "IN12_mean"]) * -1)
#Looks like the same problem exists with IN4_mean and IN9_mean: postive and negative values. Setting everything to negative to fix this. 
tecc[, "IN4_mean"] <- (abs(tecc[, "IN4_mean"]) * -1)
tecc[, "IN9_mean"] <- (abs(tecc[, "IN9_mean"]) * -1)

#drop rows where neither electrode was injecting current
tecc <- tecc[tecc$IN7_mean < threshold.nA | tecc$IN12_mean < threshold.nA, ]
tecc$IN4 <- as.character(tecc$IN4)
tecc$IN9 <- as.character(tecc$IN9)

tecc[,"Inj_Cell"] <- "LC0"
tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "Inj_Cell"] <- as.character(tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "IN4"])
tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "Inj_Cell"] <- as.character(tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "IN9"])
tecc <- check_start_end(input.df = tecc, start.time = 0, end.time = 60)invisible(readline(prompt="Press [enter] to continue"))  #https://stackoverflow.com/questions/15272916/how-to-wait-for-a-keypress-in-r
#Differs from take 3 in that this iterates over the time column instead of the recording col.
tick <- Sys.time()

aa <- find_median_deflections(input.df = tecc[tecc$Treatment == "Phase Shift",])
bb <- find_median_deflections(input.df = tecc[tecc$Treatment == "Inverted Wave",])
cc <- find_median_deflections(input.df = tecc[tecc$Treatment == "Control Wave + TEA",])
dd <- find_median_deflections(input.df = tecc[tecc$Treatment == "Aberrant Wave",])
ee <- find_median_deflections(input.df = tecc[tecc$Treatment == "Silent Control",])
ff <- find_median_deflections(input.df = tecc[tecc$Treatment == "Small Phase Angle",])

output.df <- full_join(aa, bb)
output.df <- full_join(output.df, cc)
output.df <- full_join(output.df, dd)
output.df <- full_join(output.df, ee)
output.df <- full_join(output.df, ff)

print(Sys.time()-tick)
#41 sec
tecc <- output.df

Calculate Resistances per Bennett 1966

Apparant Cell resistances: $$r_{11}=\frac{v_1}{i_1}$$

#when inj IN7...
tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "r11"] <- tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "IN4_mean"] / tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "IN7_mean"]

#when inj IN12...
tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "r11"] <- tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "IN9_mean"] / tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "IN12_mean"]

Transfer resistances: $$r_{12}=\frac{v_2}{i_1}$$

#when inj IN7...
tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "r12"] <- tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "IN9_mean"] / tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "IN7_mean"]

#when inj IN12...
tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "r12"] <- tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "IN4_mean"] / tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "IN12_mean"]

$$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}}$$

tecc <- solve_for_resistances(input.df = tecc)
tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "Coupling Coefficient"] <- 
  tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "IN9_mean"] /
  tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "IN4_mean"]

tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "Coupling Coefficient"] <- 
  tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "IN4_mean"] / 
  tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "IN9_mean"]


tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "Input Resistance"] <- 
  tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "IN4_mean"] / 
  tecc[abs(tecc$IN7_mean) > abs(tecc$IN12_mean), "IN7_mean"]

tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "Input Resistance"] <- 
  tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "IN9_mean"] / 
  tecc[abs(tecc$IN7_mean) < abs(tecc$IN12_mean), "IN12_mean"]
tecc$r11.0 <- 0
tecc$r12.0 <- 0
tecc$r1.0 <- 0
tecc$rc.0 <- 0
tecc$`Coupling Coefficient.0` <- 0
tecc$`Input Resistance.0` <- 0

tecc$interact <- interaction(tecc$Experiment, tecc$Inj_Cell)
tick <- Sys.time()
for(INTERACT in unique(tecc$interact)){
  #INTERACT =  "170623b.LC5"
  #INTERACT = "180110.LC4"
  #print(INTERACT)
  tecc[tecc$interact ==  INTERACT, c("r11.0", "r12.0", "r1.0", "rc.0", "Coupling Coefficient.0", "Input Resistance.0")] <- lapply(tecc[tecc$interact ==  INTERACT & tecc$`Time Exposed` == 0, c("r11", "r12", "r1", "rc", "Coupling Coefficient", "Input Resistance")], mean)
}
print(Sys.time() - tick)
tecc$Condition <- factor(tecc$Condition, c("Phase Offset 0 Degrees", "11.2_degrees", "22.5_degrees", "Phase Offset 45 Degrees" , "Phase Offset 90 Degrees" ,"Phase Offset 180 Degrees","Inverted Wave","TEA_ControlWave", "TEA_Silent","Aberrant Depol" ,"hold_at_-53mV"))

Prepare Two electrode voltage clamp df

Voltage Clamp Measurements of Coupling Conductance

tevc <- tevc[(tevc$IN4 == "LC4" | tevc$IN4 == "LC5") &
       (tevc$IN4 == "LC4" | tevc$IN4 == "LC5"),]
tevc <- check_start_end(input.df = tevc, start.time = 0, end.time = 60)
invisible(readline(prompt="Press [enter] to continue"))
tevc <- find_and_fix_switched_cols(input.df = tevc)
#Ensure channel cell ids are characters, not factors
tevc$IN4 <- as.character(tevc$IN4)
tevc$IN9 <- as.character(tevc$IN9)
tic <- Sys.time()
outputs <- process_tevc(input.df = tevc)

cowplot::ggsave("val.tevc.all.png", plot = cowplot::plot_grid(plotlist = outputs[[2]]), 
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 12, 
                height = 8.83,
                dpi = 300, limitsize = TRUE)

output.df <- outputs[[1]]

print(Sys.time() - tic)
tevc[tevc$Inj_Cell == "LC4", "Direction"] <- "LC4 -> LC5"
tevc[tevc$Inj_Cell == "LC5", "Direction"] <- "LC5 -> LC4"

tevc[, "interact"] <- interaction(tevc$Experiment, tevc$Direction)
#Experiment 170817b has a slope greater than 9. Dropping this experiment's data
tevc <- tevc[tevc$Experiment != "170817b",]
tevc <- tevc[!is.na(tevc),]

Prepare A type df

a <- a[(a$IN4 == "LC4" | a$IN4 == "LC5") &
       (a$IN4 == "LC4" | a$IN4 == "LC5"),]
a <- check_start_end(input.df = a, start.time = 0, end.time = 60)
invisible(readline(prompt="Press [enter] to continue"))
#read in data
#drop entries that don't go the full lenght of the experiment

a.peak <- sep_peak_end(input.df = a, current.type = "a")[[1]] %>% 
  simplify_df() %>% 
  stack_channels()

a.end <- sep_peak_end(input.df = a, current.type = "a")[[2]] %>%
  simplify_df() %>% 
  stack_channels()

a.peak$interact <- interaction(a.peak$Experiment, a.peak$Cell)
a.end$interact <- interaction(a.end$Experiment, a.end$Cell)

Prepare HTK df

htk <- htk[(htk$IN4 == "LC4" | htk$IN4 == "LC5") &
       (htk$IN4 == "LC4" | htk$IN4 == "LC5"),]
htk <- check_start_end(input.df = htk, start.time = 0, end.time = 60)
invisible(readline(prompt="Press [enter] to continue"))
htk.peak <- sep_peak_end(input.df = htk, current.type = "htk")[[1]] %>% simplify_df() %>% stack_channels()

htk.end <- sep_peak_end(input.df = htk, current.type = "htk")[[2]] %>%
  simplify_df() %>% 
  stack_channels()

htk.peak$interact <- interaction(htk.peak$Experiment, htk.peak$Cell)
htk.end$interact <- interaction(htk.end$Experiment, htk.end$Cell)

Prepare HTK LS df

htk.ls <- htk.ls[(htk.ls$IN4 == "LC4" | htk.ls$IN4 == "LC5") &
       (htk.ls$IN4 == "LC4" | htk.ls$IN4 == "LC5"),]
htk.ls <- check_start_end(input.df = htk.ls, start.time = 0, end.time = 60)
invisible(readline(prompt="Press [enter] to continue"))
htk.ls.peak <- sep_peak_end(input.df = htk.ls, current.type = "htk.ls")[[1]] %>% simplify_df() %>% stack_channels()

htk.ls.end <- sep_peak_end(input.df = htk.ls, current.type = "htk.ls")[[2]] %>%
  simplify_df() %>% 
  stack_channels()

htk.ls.peak$interact <- interaction(htk.ls.peak$Experiment, htk.ls.peak$Cell)
htk.ls.end$interact <- interaction(htk.ls.end$Experiment, htk.ls.end$Cell)
#saveRDS(tecc, file = "S:/Data_Daniel/activity_dept_change_project/out/processed.data/tecc01.rds")
#saveRDS(tevc, file = "S:/Data_Daniel/activity_dept_change_project/out/processed.data/tevc01.rds")
#saveRDS(a, file = "S:/Data_Daniel/activity_dept_change_project/out/processed.data/a01.rds")
#saveRDS(htk, file = "S:/Data_Daniel/activity_dept_change_project/out/processed.data/htk01.rds")
#saveRDS(htk.ls, file = "S:/Data_Daniel/activity_dept_change_project/out/processed.data/htkls01.rds")

Restore point ---------------------------------------------------------------

#tecc <- readRDS(file = "S:/Data_Daniel/activity_dept_change_project/out/processed.data/tecc01.rds")
#tevc <- readRDS(file = "S:/Data_Daniel/activity_dept_change_project/out/processed.data/tevc01.rds")
#a <- readRDS(file = "S:/Data_Daniel/activity_dept_change_project/out/processed.data/a01.rds")
#htk <- readRDS(file = "S:/Data_Daniel/activity_dept_change_project/out/processed.data/htk01.rds")
#htk.ls <- readRDS(file = "S:/Data_Daniel/activity_dept_change_project/out/processed.data/htkls01.rds")

Figures ---------------------------------------------------------------------

Method Validation ==========================================================

How are Input Resistance and Membrane Resistance Related?

ggplot(tecc[tecc$`Time Exposed` == 0, ], aes(x = `Input Resistance`, y = r1))+
  geom_abline(intercept = 0, slope = 1, linetype = "dotted")+
  geom_point(aes(color = Treatment))+
  geom_smooth(method = lm, se = FALSE, linetype = "dashed")+
  labs(x = "Input Resistance",  y = "Membrane Resistance")+
  theme(legend.position = "bottom")+
  xlim(0,20)+
  ylim(0,20)+
  geom_text(x = 10, y = 1, label = lm_eqn(m = lm(r1 ~ `Input Resistance`, data = tecc[tecc$`Time Exposed` == 0, ])
                                          ), parse = TRUE)+
  labs(title = "Membrane Resistance vs Input Resistance: All Data")

cowplot::ggsave("val.r1vr11.all.png", plot = ggplot2::last_plot(), 
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.3, 
                height = 6.2,
                dpi = 300, limitsize = TRUE)

ggplot(tecc[tecc$`Time Exposed` == 0 & tecc$`Input Resistance` < 13, ], aes(x = `Input Resistance`, y = r1))+
  geom_abline(intercept = 0, slope = 1, linetype = "dotted")+
  geom_point(aes(color = Treatment))+
  geom_smooth(method = lm, se = FALSE, linetype = "dashed")+
  labs(x = "Input Resistance",  y = "Membrane Resistance")+
  theme(legend.position = "bottom")+
  xlim(0,20)+
  ylim(0,20)+
  geom_text(x = 10, y = 1, label = lm_eqn(m = lm(r1 ~ `Input Resistance`, data = tecc[tecc$`Time Exposed` == 0 & tecc$`Input Resistance` < 13, ])
                                          ), parse = TRUE)+
  labs(title = "Membrane Resistance vs Input Resistance \nOutlier Removed")

cowplot::ggsave("val.r1vr11.most.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.3, 
                height = 6.2,
                dpi = 300, limitsize = TRUE)

How much overlap is there between the input resistances of each group?

ggplot(tecc[tecc$`Time Exposed` == 0, ], aes(x = Condition, y = `Input Resistance`, fill = Condition))+
  geom_boxplot()+
  geom_point()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "")+
  labs(title = "Input Resistance By Group")

cowplot::ggsave("val.r11.box.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.5, 
                height = 6,
                dpi = 300, limitsize = TRUE)

ggplot(tecc[tecc$`Time Exposed` == 0, ], aes(x = `Input Resistance`, fill = Condition))+
  geom_histogram(bins = 17)+
  labs(title = "Input Resistance is Unimodal")

cowplot::ggsave("val.r11.hist.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 4.85,
                dpi = 300, limitsize = TRUE)

Membrane resistance across conditions

ggplot(tecc[tecc$`Time Exposed` == 0 & tecc$r1 > 0, ], aes(x = Condition, y = r1, fill = Condition))+
  geom_boxplot()+
  geom_point()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "")+
  labs(title = "Membrane Resistance By Group")

cowplot::ggsave("val.r1.box.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.5, 
                height = 6,
                dpi = 300, limitsize = TRUE)

Transfer resistance across conditions

ggplot(tecc[tecc$`Time Exposed` == 0, ], aes(x = Condition, y = r12, fill = Condition))+
  geom_boxplot()+
  geom_point()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "")+
  labs(title = "Transfer Resistance By Group")

cowplot::ggsave("val.r12.box.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.5, 
                height = 6,
                dpi = 300, limitsize = TRUE)

Coupling coefficient across conditions

ggplot(tecc[tecc$`Time Exposed` == 0, ], aes(x = Condition, y = `Coupling Coefficient`, fill = Condition))+
  geom_boxplot()+
  geom_point()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "")+
  labs(title = "Coupling Coefficient By Group")

cowplot::ggsave("val.cc.box.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.5, 
                height = 6,
                dpi = 300, limitsize = TRUE)

Coupling resistance across conditions

ggplot(tecc[tecc$`Time Exposed` == 0 & tecc$rc>0, ], aes(x = Condition, y = rc, fill = Condition))+
  geom_boxplot()+
  geom_point()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "")+
  labs(title = "Coupling Resistance By Group")

cowplot::ggsave("val.rc.box.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.5, 
                height = 6,
                dpi = 300, limitsize = TRUE)

Is input resistance of one LC predictive of the other's?

rect <- tecc[tecc$`Time Exposed` == 0,c("Experiment", "Inj_Cell", "Condition", "r11","r12","r1","rc","Coupling Coefficient")]
inj.lc5 <- rect[rect$Inj_Cell == "LC5",]
inj.lc4 <- rect[rect$Inj_Cell == "LC4",]
names(inj.lc5) <- c("Experiment","Inj_Cell","Condition",
                    "LC5.r11","LC5->LC4.r12","LC5.r1",
                    "LC5->LC4.rc","LC5->LC4.Coupling Coefficient")
names(inj.lc4) <- c("Experiment","Inj_Cell","Condition",
                    "LC4.r11","LC4->LC5.r12","LC4.r1",
                    "LC4->LC5.rc","LC4->LC5.Coupling Coefficient")
rect <- cbind(inj.lc4, inj.lc5)

# [1] "Experiment"                    "Inj_Cell"                     
# [3] "Condition"                     "LC4.r11"                      
# [5] "LC4->LC5.r12"                  "LC4.r1"                       
# [7] "LC4->LC5.rc"                   "LC4->LC5.Coupling Coefficient"
# [9] "LC5.r11"                       "LC5->LC4.r12"                 
#[11] "LC5.r1"                        "LC5->LC4.rc"                  
#[13] "LC5->LC4.Coupling Coefficient"
ggplot(rect, aes(x = LC4.r11, y = LC5.r11))+
  geom_point(aes(color = Condition))+
  geom_abline(slope = 1, intercept = 0, linetype = 2)+
  geom_smooth(method = lm)+
  ylim(0,8)+
  xlim(0,8)+
  geom_text(x = 4, y = 0, label = lm_eqn(m = lm(LC5.r11 ~ LC4.r11, data = rect)
                                          ), parse = TRUE)+
  labs(title = "LC5 Input Resistance by LC4 Input Resistance")+
  theme(legend.position = "")

cowplot::ggsave("val.r11.lc5vlc4.dir.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

Is there directionality in membrane resistance?

ggplot(rect, aes(x = LC4.r1, y = LC5.r1))+
  geom_point(aes(color = Condition))+
  geom_abline(slope = 1, intercept = 0, linetype = 2)+
  geom_smooth(method = lm)+
  ylim(0,15)+
  xlim(0,15)+
  geom_text(x = 7.5, y = 0, label = lm_eqn(m = lm(LC5.r1 ~ LC4.r1, data = rect)
                                          ), parse = TRUE)+
  labs(title = "LC5 Membrane Resistance by LC4 Membrane Resistance")+
  theme(legend.position = "")

cowplot::ggsave("val.r1.lc5vlc4.dir.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

Is there directionality in transfer resistance?

ggplot(rect, aes(x = `LC4->LC5.r12`, y = `LC5->LC4.r12`))+
  geom_point(aes(color = Condition))+
  geom_abline(slope = 1, intercept = 0, linetype = 2)+
  geom_smooth(method = lm)+
  ylim(0,7)+
  xlim(0,7)+
  geom_text(x = 3.5, y = 0, label = lm_eqn(m = lm(`LC5->LC4.r12`~`LC4->LC5.r12`, data = rect)
                                          ), parse = TRUE)+
  labs(title = "LC4->LC5 Transfer Resistance by \nLC5->LC4 Transfer Resistance")+
  theme(legend.position = "")

cowplot::ggsave("val.r12.lc5vlc4.dir.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

Is there directionality in coupling resistance?

ggplot(rect, aes(x = `LC4->LC5.rc`, y = `LC5->LC4.rc`))+
  geom_point(aes(color = Condition))+
  geom_abline(slope = 1, intercept = 0, linetype = 2)+
  geom_smooth(method = lm)+
  ylim(0,11)+
  xlim(0,11)+
  geom_text(x = 5.5, y = 0, label = lm_eqn(m = lm(`LC5->LC4.rc`~`LC4->LC5.rc`, data = rect)
                                          ), parse = TRUE)+
  labs(title = "LC4->LC5 Coupling Resistance by \nLC5->LC4 Coupling Resistance")+
  theme(legend.position = "")

cowplot::ggsave("val.rc.lc5vlc4.dir.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

Is there directionality in coupling coefficient?

ggplot(rect, aes(x = `LC4->LC5.Coupling Coefficient`, y = `LC5->LC4.Coupling Coefficient`))+
  geom_point(aes(color = Condition))+
  geom_abline(slope = 1, intercept = 0, linetype = 2)+
  geom_smooth(method = lm)+
  ylim(0,1)+
  xlim(0,1)+
  geom_text(x = 0.5, y = 0, label = lm_eqn(m = lm(`LC5->LC4.Coupling Coefficient`~`LC4->LC5.Coupling Coefficient`, data = rect)
                                          ), parse = TRUE)+
  labs(title = "LC4->LC5 Coupling Coefficient \nby LC5->LC4 Coupling Coefficient")+
  theme(legend.position = "")

cowplot::ggsave("val.cc.lc5vlc4.dir.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

How close a online and offline leak subtracted measures of HTK?

merge_htk_methods <- function(input.df = htk.peak,
                              input.df.ls = htk.ls.peak){
  #input.df = htk.peak
  #input.df.ls = htk.ls.peak

  input.df <- input.df[, !(names(input.df) %in% c("Recording"))]
  input.df.ls <- input.df.ls[, !(names(input.df.ls) %in% c("Recording"))]

  #input.df$online <- FALSE
  #input.df.ls$online <- TRUE

  input.df$Experiment <- as.character(input.df$Experiment)
  input.df$IN4 <- as.character(input.df$IN4)
  input.df$IN9 <- as.character(input.df$IN9)

  input.df$target.mV <- ((round((input.df$mV)/10))*10)
  input.df.ls$target.mV <- ((round((input.df.ls$mV)/10))*10)

  input.df <- input.df %>% rename(manual.nA = "nA")
  input.df <- input.df %>% rename(manual.mV = "mV")

  input.df.ls <- input.df.ls %>% rename(online.mV = "mV")
  input.df.ls <- input.df.ls %>% rename(online.nA = "nA")

  #dplyr::full_join(input.df, input.df.ls) %>% head()

  #full_join(input.df, input.df.ls, by = c("Experiment", "Time Exposed", "Cell", "target.mV")) %>% head()

  return(dplyr::full_join(input.df, input.df.ls))
}

htk.comparison <- merge_htk_methods(input.df = htk.peak,
                                      input.df.ls = htk.ls.peak)

All the data: peak

ggplot(htk.comparison, aes(x = manual.mV, y = online.mV))+
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
  geom_point(aes(color = interact), alpha = 0.4)+
  geom_smooth(method = lm)+
  geom_text(x = -30, y = -80, label = lm_eqn(m = lm(`online.mV`~`manual.mV`, data = htk.comparison)), parse = TRUE)+
  theme(legend.position = "")+
  ylim(-80,20)+
  xlim(-80,20)+
  labs(title = "No Systematic Difference in Voltage Control")

cowplot::ggsave("val.man.pn.htk.mv.peak.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

ggplot(htk.comparison, aes(x = manual.nA, y = online.nA))+
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
  geom_point(aes(color = interact), alpha = 0.4)+
  geom_smooth(method = lm)+
  geom_text(x = 225, y = 0, label = lm_eqn(m = lm(`online.nA`~`manual.nA`, data = htk.comparison)), parse = TRUE)+
  ylim(0,450)+
  xlim(0,450)+
  theme(legend.position = "")+
  labs(title = "Measures of HTK Transient Differ Based on \nLeak Subtraction Method")

cowplot::ggsave("val.man.pn.htk.na.peak.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

Broken down by target voltage: peak

voltages <- c(-80, -70, -60, -50, -40, -30, -20, -10,  0, 10, 20)
compare.htk.methods <- list()
for (i in 1:length(voltages)){
  compare.htk.methods[[i]] <- 
    ggplot(htk.comparison[htk.comparison$target.mV == voltages[i], ], aes(x = manual.nA, y = online.nA))+
    geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
    geom_point(aes(color = interact))+
    geom_smooth(method = lm, color = "red")+
    labs(title = paste("At", as.character(voltages[i]), "mV"))+
    theme(legend.position = "")

}
cowplot::plot_grid(plotlist = compare.htk.methods)

cowplot::ggsave("val.man.pn.htk.multi.peak.mv.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.72, 
                height = 6.23,
                dpi = 300, limitsize = TRUE)

Looking at the differences: peak

htk.comparison[,"nA.diff"] <- (htk.comparison[,"online.nA"] - htk.comparison[,"manual.nA"])

lower.sd <- (mean(na.omit(htk.comparison$nA.diff)) - sd(na.omit(htk.comparison$nA.diff)))
upper.sd <- (mean(na.omit(htk.comparison$nA.diff)) + sd(na.omit(htk.comparison$nA.diff)))
q1 <- as.numeric(summary(htk.comparison$nA.diff)[2])
q2 <- as.numeric(summary(htk.comparison$nA.diff)[3])
q3 <- as.numeric(summary(htk.comparison$nA.diff)[5])

htk.comparison$target.mV <- as.factor(htk.comparison$target.mV)
ggplot(htk.comparison, aes(x = nA.diff, fill = target.mV))+
  geom_histogram(bins = 100)+
  geom_segment(aes(x = lower.sd,
                   y = -15,
                   xend = upper.sd, 
                   yend = -15), 
                   color = "red",
                   size = 1
               )+
  geom_segment(aes(x = mean(na.omit(htk.comparison$nA.diff)),
                   y = -10,
                   xend = mean(na.omit(htk.comparison$nA.diff)), 
                   yend = -20), 
                   color = "red",
                   size = 1
               )+
    geom_segment(aes(x = q1,
                   y = -5,
                   xend = q3, 
                   yend = -5), 
                   color = "blue",
                   size = 1
               )+
  geom_segment(aes(x = q2,
                   y = 0,
                   xend = q2, 
                   yend = -10), 
                   color = "blue",
                   size = 1)+
  labs(title = "Online less Offline nA: \nQuartiles Denoted in Blue \nStandard Deviation in Red")+
  theme(legend.position = "right")

cowplot::ggsave("val.man.pn.htk.hist.peak.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

All the data: end

htk.comparison <- merge_htk_methods(input.df = htk.end,
                                      input.df.ls = htk.ls.end)
ggplot(htk.comparison, aes(x = manual.mV, y = online.mV))+
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
  geom_point(aes(color = interact), alpha = 0.4)+
  geom_smooth(method = lm)+
  theme(legend.position = "")+
  ylim(-80,20)+
  xlim(-80,20)+
  labs(title = "No Systematic Difference in Voltage Control")+
  geom_text(x = -30, y = -75, label = lm_eqn(m = lm(`online.mV`~`manual.mV`, data = htk.comparison)), parse = TRUE)

cowplot::ggsave("val.man.pn.htk.mv.end.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

ggplot(htk.comparison, aes(x = manual.nA, y = online.nA))+
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
  geom_point(aes(color = interact), alpha = 0.4)+
  geom_smooth(method = lm)+
  ylim(0,450)+
  xlim(0,450)+
  theme(legend.position = "")+
  labs(title = "Measures of HTK Transient Differ Based on \nLeak Subtraction Method")+
  geom_text(x = 225, y = 0, label = lm_eqn(m = lm(`online.nA`~`manual.nA`, data = htk.comparison)), parse = TRUE)

cowplot::ggsave("val.man.pn.htk.na.end.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

Broken down by target voltage: end

voltages <- c(-80, -70, -60, -50, -40, -30, -20, -10,  0, 10, 20)
compare.htk.methods <- list()
for (i in 1:length(voltages)){
  compare.htk.methods[[i]] <- 
    ggplot(htk.comparison[htk.comparison$target.mV == voltages[i], ], aes(x = manual.nA, y = online.nA))+
    geom_abline(slope = 1, intercept = 0, linetype = "dashed")+
    geom_point(aes(color = interact))+
    geom_smooth(method = lm, color = "red")+
    labs(title = paste("At", as.character(voltages[i]), "mV"))+
    theme(legend.position = "")
}
cowplot::plot_grid(plotlist = compare.htk.methods)

cowplot::ggsave("val.man.pn.htk.multi.end.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.72, 
                height = 6.23,
                dpi = 300, limitsize = TRUE)

Looking at the differences: end

htk.comparison[,"nA.diff"] <- (htk.comparison[,"online.nA"] - htk.comparison[,"manual.nA"])

lower.sd <- (mean(na.omit(htk.comparison$nA.diff)) - sd(na.omit(htk.comparison$nA.diff)))
upper.sd <- (mean(na.omit(htk.comparison$nA.diff)) + sd(na.omit(htk.comparison$nA.diff)))
q1 <- as.numeric(summary(htk.comparison$nA.diff)[2])
q2 <- as.numeric(summary(htk.comparison$nA.diff)[3])
q3 <- as.numeric(summary(htk.comparison$nA.diff)[5])

htk.comparison$target.mV <- as.factor(htk.comparison$target.mV)
ggplot(htk.comparison, aes(x = nA.diff, fill = target.mV))+
  geom_histogram(bins = 100)+
  geom_segment(aes(x = lower.sd,
                   y = -15,
                   xend = upper.sd, 
                   yend = -15), 
                   color = "red",
                   size = 1
               )+
  geom_segment(aes(x = mean(na.omit(htk.comparison$nA.diff)),
                   y = -10,
                   xend = mean(na.omit(htk.comparison$nA.diff)), 
                   yend = -20), 
                   color = "red",
                   size = 1
               )+
    geom_segment(aes(x = q1,
                   y = -5,
                   xend = q3, 
                   yend = -5), 
                   color = "blue",
                   size = 1
               )+
  geom_segment(aes(x = q2,
                   y = 0,
                   xend = q2, 
                   yend = -10), 
                   color = "blue",
                   size = 1)+
  labs(title = "Online less Offline nA: \nQuartiles Denoted in Blue \nStandard Deviation in Red")+
  theme(legend.position = "right")

cowplot::ggsave("val.man.pn.htk.hist.end.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 6.44, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

Visualize Treatments =======================================================

1. Phase Angle

All groups

TREATMENT <- "Phase Shift"

return.plts <- list()
# Un-normalized ---------------------------------------------------------------
return.plts[[1]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = (rc),
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "rc",
    title = paste(TREATMENT, "Coupling Resistance")
  ) +
  theme_Publication()+
  facet_grid(~Condition)

return.plts[[2]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "CC",
    title = paste(TREATMENT, "Coupling Coefficient")
  ) +
  theme_Publication()+
  facet_grid(~Condition)

return.plts[[3]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Input Resistance`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Input Resistance",
    title = paste(TREATMENT, "Input Resistance")
  ) +
  theme_Publication()+
  facet_grid(~Condition)

# Normalized ------------------------------------------------------------------
return.plts[[4]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = (rc / rc.0),
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized rc",
    title = paste(TREATMENT, "Normalized Coupling Resistance")
  ) +
  theme_Publication()+
  facet_grid(~Condition)

return.plts[[5]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient` / `Coupling Coefficient.0`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized CC",
    title = paste(TREATMENT, "Normalized Coupling Coefficient")
  ) +
  theme_Publication()+
  facet_grid(~Condition)

return.plts[[6]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Input Resistance` / `Input Resistance.0`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized Input Resistance",
    title = paste(TREATMENT, "Normalized Input Resistance")
  ) +
  theme_Publication()+
  facet_grid(~Condition)

# Group interaction style plots -----------------------------------------------

return.plts[[7]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient` / `Coupling Coefficient.0`,
    color = Condition, group = Condition
  )) +
  geom_point(position = position_jitter(width = 0.1), size = 2) +
  stat_summary(fun.y = mean, geom = "point", pch = 12, size = 5) +
  stat_summary(fun.y = mean, geom = "line", size = 1) +
  labs(color = "Phase \nAngle") +
  labs(title = "Interaction of Treatment:Time \non Coupling Coefficient (normalized)", x = "Time in Minutes", y = "Coupling Coefficient") +
  theme_Publication() +
  ylim(0, 2)

return.plts[[8]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT & tecc$`Time Exposed` < 90, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient`,
    color = Condition, group = interaction(Condition, `Time Exposed`)
  )) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2) +
  labs(title = "Coupling Coefficients over Time", x = "Time in Minutes", y = "Coupling Coefficient") +
  facet_grid(~Condition) +
  theme_Publication() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.position = "") +
  ylim(0, 1)

# -----------------------------------------------------------------------------
save_fig_list(input.list = return.plts,
                          title.prefix = "phaseshift",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))
ggplot(tevc[tevc$Treatment == "Phase Shift" & tevc$`Time Exposed` < 90,], aes(x = `Time Exposed`, y = slopes, color = Condition, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")+
  theme(legend.position = "")+
  facet_grid(~Condition)

cowplot::ggsave("phaseshift.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)

2. Small Phase Angles

#show_tecc_plots(TREATMENT = "Small Phase Angle")
TREATMENT <- "Small Phase Angle"

return.plts <- list()
# Un-normalized ---------------------------------------------------------------
return.plts[[1]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = (rc),
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "rc",
    title = paste(TREATMENT, "Coupling Resistance")
  ) +
  theme_Publication()

return.plts[[2]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "CC",
    title = paste(TREATMENT, "Coupling Coefficient")
  ) +
  theme_Publication()

return.plts[[3]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = `Input Resistance`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Input Resistance",
    title = paste(TREATMENT, "Input Resistance")
  ) +
  theme_Publication()

# Normalized ------------------------------------------------------------------
return.plts[[4]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = (rc / rc.0),
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized rc",
    title = paste(TREATMENT, "Normalized Coupling Resistance")
  ) +
  theme_Publication()

return.plts[[5]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient` / `Coupling Coefficient.0`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized CC",
    title = paste(TREATMENT, "Normalized Coupling Coefficient")
  ) +
  theme_Publication()

return.plts[[6]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = `Input Resistance` / `Input Resistance.0`,
    color = Condition,
    group = interact
  )) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1, linetype = 2) +
  labs(
    y = "Normalized Input Resistance",
    title = paste(TREATMENT, "Normalized Input Resistance")
  ) +
  theme_Publication()

# Group interaction style plots -----------------------------------------------

return.plts[[7]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = `Time Exposed`,
    y = `Coupling Coefficient` / `Coupling Coefficient.0`,
    color = interact, group = Condition
  )) +
  geom_point(position = position_jitter(width = 0.1), size = 2) +
  stat_summary(fun.y = mean, geom = "point", pch = 12, size = 5) +
  stat_summary(fun.y = mean, geom = "line", size = 1) +
  labs(color = "Phase \nAngle") +
  labs(title = "Interaction of Treatment:Time \non Coupling Coefficient (normalized)", x = "Time in Minutes", y = "Coupling Coefficient") +
  theme_Publication() +
  ylim(0, 2)

return.plts[[8]] <-
  ggplot(tecc[tecc$Treatment == TREATMENT, ], aes(
    x = Condition,
    y = `Coupling Coefficient`,
    color = interact, group = Condition
  )) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2) +
  labs(title = "Coupling Coefficients over Time", x = "Time in Minutes", y = "Coupling Coefficient") +
  facet_grid(~`Time Exposed`) +
  theme_Publication() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.position = "") +
  ylim(0, 1)
# -----------------------------------------------------------------------------
save_fig_list(input.list = return.plts,
                          title.prefix = "smallphase",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))
ggplot(tevc[tevc$Treatment == "Small Phase Angle",], aes(x = `Time Exposed`, y = slopes, color = interact, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")

cowplot::ggsave("smallphase.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)
return.plts <- standard_ionic_plts(current.treatment = "Small Phase Angle",
                                title.prefix = "A Peak",
                                input.df = a.peak[a.peak$Treatment == "Small Phase Angle",])

save_fig_list(input.list = return.plts,
                          title.prefix = "smallphase.a.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts(current.treatment = "Small Phase Angle",
                                title.prefix = "A End",
                                input.df = a.end[a.end$Treatment == "Small Phase Angle",])

save_fig_list(input.list = return.plts,
                          title.prefix = "smallphase.end.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts(current.treatment = "Small Phase Angle",
                                title.prefix = "HTK Peak",
                                input.df = htk.peak[htk.peak$Treatment == "Small Phase Angle",])

save_fig_list(input.list = return.plts,
                          title.prefix = "smallphase.htk.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts(current.treatment = "Small Phase Angle",
                                title.prefix = "HTK End",
                                input.df = htk.end[htk.end$Treatment == "Small Phase Angle",])

save_fig_list(input.list = return.plts,
                          title.prefix = "smallphase.htk.end",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

3. Inverted Wave

inv.plts <- show_tecc_plots(TREATMENT = "Inverted Wave")
save_fig_list(input.list = inv.plts,
                          title.prefix = "inv",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))
ggplot(tevc[tevc$Treatment == "Inverted Wave",], aes(x = `Time Exposed`, y = slopes, color = interact, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")

cowplot::ggsave("inv.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)
return.plts <- standard_ionic_plts_1c(current.treatment = "Inverted Wave",
                                title.prefix = "A Peak",
                                input.df = a.peak[a.peak$Treatment == "Inverted Wave",])

save_fig_list(input.list = return.plts,
                          title.prefix = "inv.a.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Inverted Wave",
                                title.prefix = "A End",
                                input.df = a.end[a.end$Treatment == "Inverted Wave",])

save_fig_list(input.list = return.plts,
                          title.prefix = "inv.a.end",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Inverted Wave",
                                title.prefix = "HTK Peak",
                                input.df = htk.peak[htk.peak$Treatment == "Inverted Wave",])

save_fig_list(input.list = return.plts,
                          title.prefix = "inv.htk.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Inverted Wave",
                                title.prefix = "HTK End",
                                input.df = htk.end[htk.end$Treatment == "Inverted Wave",])

save_fig_list(input.list = return.plts,
                          title.prefix = "inv.htk.end",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

4. Aberrant Wave

ab.plts <- show_tecc_plots(TREATMENT = "Aberrant Wave")
save_fig_list(input.list = ab.plts,
                          title.prefix = "ab",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))
ggplot(tevc[tevc$Treatment == "Aberrant Wave",], aes(x = `Time Exposed`, y = slopes, color = interact, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")

cowplot::ggsave("ab.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)
return.plts <- standard_ionic_plts_1c(current.treatment = "Aberrant Wave",
                                title.prefix = "A Peak",
                                input.df = a.peak[a.peak$Treatment == "Aberrant Wave",])

save_fig_list(input.list = return.plts,
                          title.prefix = "ab.a.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Aberrant Wave",
                                title.prefix = "A End",
                                input.df = a.end[a.end$Treatment == "Aberrant Wave",])

save_fig_list(input.list = return.plts,
                          title.prefix = "ab.a.end",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Aberrant Wave",
                                title.prefix = "HTK Peak",
                                input.df = htk.peak[htk.peak$Treatment == "Aberrant Wave",])

save_fig_list(input.list = return.plts,
                          title.prefix = "ab.htk.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Aberrant Wave",
                                title.prefix = "HTK End",
                                input.df = htk.end[htk.end$Treatment == "Aberrant Wave",])

save_fig_list(input.list = return.plts,
                          title.prefix = "ab.htk.end",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

5. Silent

sh.plts <- show_tecc_plots(TREATMENT = "Silent Control")
save_fig_list(input.list = sh.plts,
                          title.prefix = "sh",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))
ggplot(tevc[tevc$Treatment == "Silent Control",], aes(x = `Time Exposed`, y = slopes, color = interact, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")

cowplot::ggsave("sh.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)
return.plts <- standard_ionic_plts_1c(current.treatment = "Silent Control",
                                title.prefix = "A Peak",
                                input.df = a.peak[a.peak$Treatment == "Silent Control",])

save_fig_list(input.list = return.plts,
                          title.prefix = "sh.a.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Silent Control",
                                title.prefix = "A End",
                                input.df = a.end[a.end$Treatment == "Silent Control",])

save_fig_list(input.list = return.plts,
                          title.prefix = "sh.a.end",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Silent Control",
                                title.prefix = "HTK Peak",
                                input.df = htk.peak[htk.peak$Treatment == "Silent Control",])

save_fig_list(input.list = return.plts,
                          title.prefix = "sh.htk.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Silent Control",
                                title.prefix = "HTK End",
                                input.df = htk.end[htk.end$Treatment == "Silent Control",])

save_fig_list(input.list = return.plts,
                          title.prefix = "sh.htk.end",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

6. Control in TEA

ctea.plts <- show_tecc_plots(TREATMENT = "Control Wave + TEA")
save_fig_list(input.list = ctea.plts,
                          title.prefix = "ctea",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))
ggplot(tevc[tevc$Treatment == "Control Wave + TEA",], aes(x = `Time Exposed`, y = slopes, color = interact, group = interaction(Experiment,Inj_Cell)))+
  geom_point()+
  geom_line()+
  labs(y = "Ratio of post / pre-synaptic current", title = "Change Over Time in the Small Phase Angle Group")

cowplot::ggsave("ctea.tevc.png", plot = ggplot2::last_plot(),
                device = NULL, path = paste0(
                  use.path, "180315_act_dept_report/figures/"), 
                scale = 1, 
                width = 9.86, 
                height = 5.86,
                dpi = 300, limitsize = TRUE)
return.plts <- standard_ionic_plts_1c(current.treatment = "Control Wave + TEA",
                                title.prefix = "A Peak",
                                input.df = a.peak[a.peak$Treatment == "Control Wave + TEA",])

save_fig_list(input.list = return.plts,
                          title.prefix = "ctea.a.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Control Wave + TEA",
                                title.prefix = "A End",
                                input.df = a.end[a.end$Treatment == "Control Wave + TEA",])

save_fig_list(input.list = return.plts,
                          title.prefix = "ctea.a.end",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Control Wave + TEA",
                                title.prefix = "HTK Peak",
                                input.df = htk.peak[htk.peak$Treatment == "Control Wave + TEA",])

save_fig_list(input.list = return.plts,
                          title.prefix = "ctea.htk.peak",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))

return.plts <- standard_ionic_plts_1c(current.treatment = "Control Wave + TEA",
                                title.prefix = "HTK End",
                                input.df = htk.end[htk.end$Treatment == "Control Wave + TEA",])

save_fig_list(input.list = return.plts,
                          title.prefix = "ctea.htk.end",
                          file.type = ".png",
                          save.path = paste0(use.path, "180315_act_dept_report/figures/"))


danielkick/esynvmod documentation built on May 17, 2019, 7:02 p.m.