forked from BriVandrey/lec_to_mec
-
Notifications
You must be signed in to change notification settings - Fork 0
/
pharmacology.R
307 lines (245 loc) · 15.2 KB
/
pharmacology.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
# script used to analyse data reported in https://doi.org/10.1101/2022.08.25.505296
# data and analysis scripts can be found at www.github.com/BriVandrey/lec-to-mec-analysis
# to use, download files and set working directory (***) to relevant path for your machine on line 14
# outputs: statistical results are saved as .txt files, plots are saved as .eps files
# Contact: [email protected]
# import libraries
library(rstatix)
library(tidyverse)
library(gridExtra)
library(ggpubr)
#set working directory - add relevant path for your machine
wd = "~/Documents/lecmec manuscript/analysis"
#wd = " " # *** Set working directory here
setwd(wd)
#function to create and move into subdirectory
create_dir <- function(main_dir, sub_dir){
new_dir <- file.path(main_dir, sub_dir)
dir.create(new_dir, showWarnings=FALSE) #ignore if already exists
setwd(new_dir)
}
#-----------------------------------------------------------------------------------------------------------
#PLOTTING UTILITY-------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------
#function for producing lineplot with individual observations - for pharmacology data
pharma_lp <- function(df, col, filter, xlabel1, xlabel2, ylabel, filename)
{
# handle col variables
col2 <- enquo(col)
#filter for cell type
df2 <- filter(df, type==filter)
#Turn into a factor with the levels in the correct order
df2$drug <- as.character(df2$drug)
df2$drug <- factor(df2$drug, levels=unique(df2$drug))
#summarise data
sum <-group_by(df2, drug, cell) %>%
summarise(mean_col = mean(!!col2))
sum_mean_col <-group_by(sum, drug) %>%
get_summary_stats(type = "mean_se")
#lineplot
epsp_lp <- ggplot(sum_mean_col, aes(x=drug, y=mean))+
theme_classic()+
geom_line(data = sum, inherit.aes = FALSE, aes(x = drug, y = mean_col, group = cell), size =0.6, colour="grey74")+ #individual traces
geom_point(data = sum, inherit.aes = FALSE, aes(x = drug, y = mean_col, group = cell), size = 2, shape=111, colour="grey74") + #individual points
geom_point(size = 1) + #mean points
geom_line(aes(group=1), size = 0.6) + #mean trace
geom_errorbar(aes(ymax = mean+se, ymin = mean-se), width = 0.2, size = 0.6) + #error bars
scale_x_discrete(labels=c("Baseline", xlabel1, xlabel2)) +
xlab("") + ylab(ylabel)+
ylim(0,NA)
epsp_lp + theme(axis.text.x = element_text(size = 10, angle=-50, color='black'),
axis.text.y = element_text(size = 10, color='black'),
axis.title.y = element_text(size=12, margin= margin(t=0,r=0,b=0,l=0)))
ggsave(file=filename, width =35, height =45, units="mm") #save plot
}
# function to generate scatterplot with simple linear model
corr_plot <- function(df, x_col, y_col, x_label, y_label, filename){
# handle col variables
x_var <- enquo(x_col)
y_var <- enquo(y_col)
# plot data - option to colour code by cell type, uncommented
plt <- ggplot(df, aes(x=!!x_var, y=!!y_var)) +
geom_point(fill=NA, shape=1, size=2, na.rm=TRUE, show.legend=FALSE)+
#aes(colour=factor(type)) +
labs(title="", x=x_label, y=y_label) +
scale_color_brewer(palette="BuPu")+
xlim(-80, -50)+ #limits set for figures used in paper
ylim(0, 4.5)+
theme_classic()
plt + theme(axis.text.x = element_text(size = 12, angle=-40, color='black'),
axis.text.y = element_text(size = 12, color='black'),
axis.title.y = element_text(size=12, margin= margin(t=0,r=14,b=0,l=0)),
axis.title.x = element_text(size=12, margin= margin(t=14,r=0,b=0,l=0)),
legend.position = 'none') +
geom_smooth(method='glm')
ggsave(file=filename, width =40, height =45, units='mm') #save plot
}
# function to generate scatterplot with simple linear model - coloured by group
corr_plot_grouped <- function(df, x_col, y_col, x_label, y_label, filename){
# handle col variables
x_var <- enquo(x_col)
y_var <- enquo(y_col)
# plot data - option to colour code by cell type, uncommented
plt <- ggplot(df, aes(x=!!x_var, y=!!y_var)) +
geom_point(fill=NA, shape=1, size=2, na.rm=TRUE, show.legend=FALSE)+
aes(colour=factor(type)) +
labs(title="", x=x_label, y=y_label) +
#scale_color_brewer(palette="BuPu")+
xlim(-80, -50)+ #limits set for figures used in paper
ylim(0, 4.5)+
theme_classic()
plt + theme(axis.text.x = element_text(size = 12, angle=-40, color='black'),
axis.text.y = element_text(size = 12, color='black'),
axis.title.y = element_text(size=12, margin= margin(t=0,r=14,b=0,l=0)),
axis.title.x = element_text(size=12, margin= margin(t=14,r=0,b=0,l=0)),
legend.position = 'none') +
geom_smooth(method='glm')
ggsave(file=filename, width =40, height =45, units='mm') #save plot
}
#-----------------------------------------------------------------------------------------------------------
#ANALYSIS UTILITY-------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------
# function to perform Friedmans repeated measures test on average data
avg_fried <- function(df, cell_type, filename){
df <- filter(df, type == cell_type) #filter by cell-type
df<- as.matrix(df); df<- as.data.frame(df) # reformat data to correct type
f_epsp <- friedman_test(data=df, epsp~drug|cell)
eff_epsp<- friedman_effsize(data=df, epsp ~ drug|cell)
f_ipsp <- friedman_test(data=df, ipsp_abs~drug|cell)
eff_ipsp <- friedman_effsize(data=df, ipsp_abs ~ drug|cell)
df$epsp <- as.numeric(as.character(df$epsp)); df$ipsp_abs <- as.numeric(as.character(df$ipsp_abs)); # re-make numeric for pairwise tests
pw_epsp <- wilcox_test(data=df, epsp~drug, paired = TRUE, p.adjust.method = "bonferroni")
pw_ipsp <- wilcox_test(data=df, ipsp_abs~drug, paired = TRUE, p.adjust.method = "bonferroni")
df <- na.omit(df) # omit any empty rows for epsp halfwidth
df<- as.matrix(df); df<- as.data.frame(df) #reformat again
f_epsp_hw <- friedman_test(data=df, epsp_halfwidth~drug|cell)
eff_epsp_hw<- friedman_effsize(data=df, epsp_halfwidth ~ drug|cell)
df$epsp_halfwidth <- as.numeric(as.character(df$epsp_halfwidth))
pw_epsp_hw <- wilcox_test(data=df, epsp_halfwidth~drug, paired = TRUE, p.adjust.method = "bonferroni")
results <- list(f_epsp, eff_epsp, pw_epsp,
f_ipsp, eff_ipsp, pw_ipsp,
f_epsp_hw,eff_epsp_hw, pw_epsp_hw)
capture.output(results, file=filename)
}
# function to perform Friedmans repeated measures test on sweep data - epsp & ipsp amplitude
cell_by_cell_fried <- function (df, cell_name, filename){
df<- filter(df, cell==cell_name)
#friedmans and posthoc tests
f_epsp <- friedman_test(data=df, epsp~drug|pulse)
pw_epsp <- wilcox_test(data=df, epsp~drug, paired = TRUE, p.adjust.method = "bonferroni")
eff_epsp<- friedman_effsize(data=df, epsp ~ drug|pulse)
f_ipsp <- friedman_test(data=df, ipsp_abs~drug|pulse)
eff_ipsp <- friedman_effsize(data=df, ipsp_abs ~ drug|pulse)
pw_ipsp <- wilcox_test(data=df, ipsp_abs~drug, paired = TRUE, p.adjust.method = "bonferroni")
results <- list(f_epsp, eff_epsp, pw_epsp,
f_ipsp, eff_ipsp, pw_ipsp)
capture.output(results, file=filename)
}
# function to perform Friedmans repeated measures test on sweep data - halfwidth
cell_by_cell_fried_hw <- function (df, cell_name, filename){
df<- filter(df, cell==cell_name)
df<- na.omit(df)
df <- as.matrix(df); df <- as.data.frame(df)
#friedmans and posthoc tests
f_epsp_hw <- friedman_test(data=df, epsp_halfwidth~drug|pulse)
eff_epsp_hw<- friedman_effsize(data=df, epsp_halfwidth ~ drug|pulse)
df$epsp_halfwidth <- as.numeric(as.character(df$epsp_halfwidth)) # re-make numeric for pairwise tests
pw_epsp_hw <- wilcox_test(data=df, epsp_halfwidth~drug, paired = TRUE, p.adjust.method = "bonferroni")
results <- list(f_epsp_hw, eff_epsp_hw, pw_epsp_hw)
capture.output(results, file=filename)
}
# ---------------------------------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------------------------------
# import data from .csv files and subset for plotting
glut_df <- read.csv(paste(wd, "data/glut_pharma.csv", sep="/"), head=TRUE, sep=",") #NBQX and APV - AMPA receptor antags
gaba_df <- read.csv(paste(wd, "data/gaba_pharma.csv", sep="/"), head=TRUE, sep=",") #Gabazine and CGP - GABA receptor antags
gaba_df_pulses<- read.csv(paste(wd, "data/gaba_pharma_individual_traces.csv", sep="/"), head=TRUE, sep=",") # individual traces - GABA receptor antags
ttx_df <- read.csv(paste(wd, "data/ttx_pharma.csv", sep="/"), head=TRUE, sep=",") #ttx and ap4 - monosynaptic expts
#add absolute values for ipsp
glut_df$ipsp_abs <- abs(glut_df$ipsp)
gaba_df$ipsp_abs <- abs(gaba_df$ipsp)
gaba_df_pulses$ipsp_abs <- abs(gaba_df_pulses$ipsp)
ttx_df$ipsp_abs <- abs(ttx_df$ipsp)
#subset based on order of drug application
cgp <- subset(gaba_df, pharma=="gaba_cgpgz")#CGP -> Gabazine
gz <- subset(gaba_df, pharma=="gaba_gzcgp")#Gabazine -> CGP
#create directory for outputs
create_dir(wd, "pharmacology")
#-----------------------------------------------------------------------------------------------------------
#FIGURE 5 --------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------
# NBQX + APV - L2 SCs & L2 PCs - Panels A & B
pharma_lp(glut_df, epsp, "l2_sc", "NBQX", "APV", "EPSP (mV)", "l2_sc_nbqx_apv_epsp.eps")
pharma_lp(glut_df, epsp, "l2_pc", "NBQX", "APV", "EPSP (mV)", "l2_pc_nbqx_apv_epsp.eps")
pharma_lp(glut_df, ipsp_abs, "l2_sc", "NBQX", "APV", "IPSP (mV)", "l2_sc_nbqx_apv_ipsp.eps")
pharma_lp(glut_df, ipsp_abs, "l2_pc", "NBQX", "APV", "IPSP (mV)", "l2_pc_nbqx_apv_ipsp.eps")
# TTX & 4AP - L2 SCs & L2 PCs - Panels E & F
pharma_lp(ttx_df, epsp, "l2_sc", "TTX", "4-AP", "EPSP (mV)", "l2_sc_ttx_4ap_epsp.eps")
pharma_lp(ttx_df, epsp, "l2_pc", "TTX", "4-AP", "EPSP (mV)", "l2_pc_ttx_4ap_epsp.eps")
#-----------------------------------------------------------------------------------------------------------
#FIGURE 6 --------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------
# statistics using averages from multiple traces------------------------------------------------------------
cgp_nas_removed <- na.omit(cgp) # manually remove NAs - required for Friedmans
# friedmans tests
# stellate cells - gabazine, CGP
avg_fried(gz, 'l2_sc', 'friedmans_l2_sc_gzcgp.txt')
avg_fried(gz, 'l2_pc', 'friedmans_l2_pc_gzcgp.txt')
avg_fried(cgp_nas_removed, 'l2_sc', 'friedmans_l2_sc_cgpgz.txt')
# statistics using individual traces----------------------------------------------------------------------------
# cell by cell analysis ----------------------------------------------------------------------------------------
create_dir(getwd(), "pulses_stats")
#Friedmans repeated measures
# gabazine -> cgp, stellates, cell by cell
cell_by_cell_fried(gaba_df_pulses, "1237_S2C2", "1237_S2C2_fried.txt")
cell_by_cell_fried(gaba_df_pulses, "1241_S2C2", "1241_S2C2_fried.txt")
cell_by_cell_fried_hw(gaba_df_pulses, "1241_S2C2", "1241_S2C2_fried_halfwidth.txt")
cell_by_cell_fried(gaba_df_pulses, "1356_S5C1", "1356_S5C1_fried.txt")
cell_by_cell_fried_hw(gaba_df_pulses, "1356_S5C1", "1356_S5C1_fried_halfwidth.txt")
cell_by_cell_fried(gaba_df_pulses, "1355_S2C1", "1355_S2C1_fried.txt")
cell_by_cell_fried(gaba_df_pulses, "1354_S4C1", "1354_S4C1_fried.txt")
cell_by_cell_fried_hw(gaba_df_pulses, "1354_S4C1", "1354_S4C1_fried_halfwidth.txt")
# cgp -> gabazine, stellates, cell by cell
cell_by_cell_fried(gaba_df_pulses, "1309_S1C2", "1309_S1C2_fried.txt")
cell_by_cell_fried(gaba_df_pulses, "1309_S4C1", "1309_S4C1_fried.txt")
cell_by_cell_fried(gaba_df_pulses, "180321_S1C1", "180321_S1C1_fried.txt")
cell_by_cell_fried(gaba_df_pulses, "180321_S2C1", "180321_S2C1_fried.txt")
# gabazine -> cgp, pyramidal cells, cell by cell
cell_by_cell_fried(gaba_df_pulses, "1312_S2C1", "1312_S2C1_fried.txt")
cell_by_cell_fried(gaba_df_pulses, "1354_S3C1", "1354_S3C1_fried.txt")
cell_by_cell_fried_hw(gaba_df_pulses, "1354_S3C1", "1354_S3C1_fried_halfwidth.txt")
cell_by_cell_fried(gaba_df_pulses, "1355_S1C1", "1355_S1C1_fried.txt")
cell_by_cell_fried_hw(gaba_df_pulses, "1355_S1C1", "1355_S1C1_fried_halfwidth.txt")
cell_by_cell_fried(gaba_df_pulses, "1355_S4C1", "1355_S4C1_fried.txt")
cell_by_cell_fried_hw(gaba_df_pulses, "1355_S4C1", "1355_S4C1_fried_halfwidth.txt")
cell_by_cell_fried(gaba_df_pulses, "1356_S1C1", "1356_S1C1_fried.txt")
cell_by_cell_fried_hw(gaba_df_pulses, "1356_S1C1", "1356_S1C1_fried_halfwidth.txt")
# gabazine -> cgp, pyramidal cells, *** only one cell
cell_by_cell_fried(gaba_df_pulses, "1309_S2C1", "1309_S2C1_fried.txt")
cell_by_cell_fried_hw(gaba_df_pulses, "1309_S2C1", "1309_S2C1_fried_halfwidth.txt")
#--------------------------------------------------------------------------------------------------------------
setwd("..") # return to main pharmacology folder for plots
# lineplots
# Gabazine + CGP 55485 - L2 SCs & L2 PCs - Panels A & B
pharma_lp(gz, epsp, "l2_sc", "Gabazine", "CGP", "EPSP (mV)", "l2_sc_gz_cgp_epsp.eps")
pharma_lp(gz, epsp, "l2_pc", "Gabazine", "CGP", "EPSP (mV)", "l2_pc_gz_cgp_epsp.eps")
pharma_lp(gz, epsp_halfwidth, "l2_sc", "Gabazine", "CGP", "EPSP Hafwidth (ms)", "l2_sc_gz_cgp_epsp_hw.eps")
pharma_lp(gz, epsp_halfwidth, "l2_pc", "Gabazine", "CGP", "EPSP Hafwidth (ms)", "l2_pc_gz_cgp_epsp_hw.eps")
pharma_lp(gz, ipsp_abs, "l2_sc", "Gabazine", "CGP", "IPSP (mV)", "l2_sc_gz_cgp_ipsp.eps")
pharma_lp(gz, ipsp_abs, "l2_pc", "Gabazine", "CGP", "IPSP (mV)", "l2_pc_gz_cgp_ipsp.eps")
# CGP 55485 + Gabazine - L2 SCs & L2 PCs - Panels C & D
pharma_lp(cgp, epsp, "l2_sc", "CGP", "Gabazine", "EPSP (mV)", "l2_sc_cgp_gz_epsp.eps")
pharma_lp(cgp, epsp, "l2_pc", "CGP", "Gabazine", "EPSP (mV)", "l2_pc_cgp_gz_epsp.eps")
pharma_lp(cgp, epsp_halfwidth, "l2_sc", "CGP", "Gabazine", "EPSP Hafwidth (ms)", "l2_sc_cgp_gz_epsp_hw.eps")
pharma_lp(cgp, epsp_halfwidth, "l2_pc", "CGP", "Gabazine", "EPSP Hafwidth (ms)", "l2_pc_cgp_gz_epsp_hw.eps")
pharma_lp(cgp, ipsp_abs, "l2_sc", "CGP", "Gabazine", "IPSP (mV)", "l2_sc_cgp_gz_ipsp.eps")
pharma_lp(cgp, ipsp_abs, "l2_pc", "CGP", "Gabazine", "IPSP (mV)", "l2_pc_cgp_gz_ipsp.eps")
#-----------------------------------------------------------------------------------------------------------
#FIGURE 7 --------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------
# NBQX + APV - L1 & L2 INS - Panel E
pharma_lp(glut_df, epsp, "l1_in", "NBQX", "APV", "EPSP (mV)", "l1_in_nbqx_apv_epsp.eps")
pharma_lp(glut_df, epsp, "l2_in", "NBQX", "APV", "EPSP (mV)", "l2_in_nbqx_apv_epsp.eps")
# TTX & 4AP - L1 & L2 INS - Panel F
pharma_lp(ttx_df, epsp, "l1_in", "TTX", "4-AP", "EPSP (mV)", "l1_in_ttx_4ap_epsp.eps")
pharma_lp(ttx_df, epsp, "l2_in", "TTX", "4-AP", "EPSP (mV)", "l2_in_ttx_4ap_epsp.eps")