Updated Figure S1 to include signal traces

Analysis/01_Analysis_Dataset1.R

@@ -599,9 +599,8 @@ plot_bayes_mima <- function(data, groups) {
           panel.grid.minor = element_blank(),
           axis.line = element_line(colour = "black"),
           axis.title = element_blank(),
-          strip.text.x = element_text(size = 16),
           strip.background = element_rect(colour="NA", fill="NA"),
-          strip.text.y = element_blank(),
+          strip.text = element_blank(),
           panel.border = element_blank(),
           text = element_text(size=16),
           axis.text = element_text(colour="black"))
@@ -621,9 +620,8 @@ plot_bayes_mima <- function(data, groups) {
           panel.grid.minor = element_blank(),
           axis.line = element_line(colour = "black"),
           axis.title=element_blank(),
-          strip.text.x = element_text(size = 16),
           strip.background = element_rect(colour="NA", fill="NA"),
-          strip.text.y = element_blank(),
+          strip.text = element_blank(),
           panel.border = element_blank(),
           text = element_text(size=16),
           axis.text=element_text(colour="black"))
@@ -645,9 +643,8 @@ plot_bayes_mima <- function(data, groups) {
           axis.title=element_blank(),
           plot.tag = element_text(size = 14, hjust = 0.5, vjust = 0.5),
           plot.tag.position = c(0.02, 1 - 0.03),
-          strip.text.x = element_text(size = 16),
           strip.background = element_rect(colour="NA", fill="NA"),
-          strip.text.y = element_blank(),
+          strip.text = element_blank(),
           panel.border = element_blank(),
           text = element_text(size=16),
           axis.text=element_text(colour="black"))
@@ -661,9 +658,9 @@ plot_bayes_mima <- function(data, groups) {
 ### Figure S1 #####
 Fly_mean <- plot_bayes_mima(data = data2, groups = c(creeks, 'Bee', 'PID'))
 
-FigS1A <- egg::ggarrange(Fly_mean$fen, Fly_mean$bee, Fly_mean$pid, clip = FALSE,
-                      labels = c('', '', expression(x ~ 10^3)), label.args = list(vjust=2),
-                      heights = 5, nrow = 1, ncol=3, widths = c(5, 1, 1), debug = FALSE)
+FigS1B <- egg::ggarrange(Fly_mean$fen, Fly_mean$bee, Fly_mean$pid, clip = FALSE,
+                      labels = c('', '', expression(x ~ 10^3)), label.args = list(vjust=0),
+                      nrow = 1, ncol=3, widths = c(5, 1, 1), heights = c(6), debug = FALSE)
 
 y.grob <- textGrob(paste("Response strength (mV)"), 
                    gp=gpar(fontsize=18), rot=90)
@@ -672,7 +669,7 @@ x.grob <- textGrob("Pulse duration (ms)",
                    gp=gpar(fontsize=18))
 
 #add to plot
-FigS1A <- grid.arrange(arrangeGrob(FigS1A, left = y.grob, bottom = x.grob))
+FigS1B <- grid.arrange(arrangeGrob(FigS1B, left = y.grob, bottom = x.grob))
 
 
 
@@ -947,6 +944,215 @@ Fig1C <- egg::ggarrange(traces300$traceW, traces300$traceB, traces300$traceP,
                         heights = 5, nrow = 1, ncol=3, widths = c(5, 1, 1), debug = FALSE)
 
 
+########################
+
+data <- data2
+data <- data %>%
+  filter(Group %in% c(creeks, "Bee", "PID")) %>%
+  filter(Odor %in% c("2-heptanone"), !Frequency %in% 10, !PulseDuration %in% 300)
+
+data <- droplevels(data)
+
+data$Odor <- factor(data$Odor, levels = levels(data$Odor),
+                    labels = c("2-heptanone"))
+
+# calculate the mean over all odor-pulseduration configurations
+d_mean <- data %>%
+  group_by(Group, Winged, Odor, PulseDuration) %>%
+  summarise_at(.funs = median, .vars=paste0("t", 1:5000))
+
+
+# calculate sample size (n) over all odor-pulseduration configurations
+d_n <- data %>%
+  group_by(Group, Winged, Odor, PulseDuration) %>%
+  summarise(N=n()) %>%
+  group_by(Group, Winged) %>%
+  summarise(N=unique(N))
+
+d_max <- d_mean %>%
+  group_by(Group, Odor, PulseDuration) %>%
+  summarise(Max=max(t400))
+
+
+# create d_all and add the sample size
+d_all <- d_mean
+
+d_all$PulseDuration <- factor(d_all$PulseDuration, levels = c('15', '30', '150'))
+d_all$PulseDur <-  factor(d_all$PulseDuration)
+
+
+data_long <- gather(d_all, time, Voltage, t1:t5000, factor_key=TRUE)
+data_long$time <- (as.numeric(gsub("t", "", data_long$time)) - 200)/1000
+
+
+
+
+
+print(d_n)
+
+geom.text.size <- 5.2
+data_long$Winged <- factor(data_long$Winged,
+                           levels = levels(data_long$Winged),
+                           labels = c("Full-winged", "Wing-reduced", "none"))
+
+
+data_segm1 <- data.frame(y=c(0, 0), yend=c(0.1, 0), x=c(Inf, 2), xend=c(Inf, 3),
+                         PulseDuration = factor(150, levels = c(15, 30, 150)),
+                         Group=factor(c("Lug", 'Mt Burns'), levels = creeks),
+                         Odor=c("2-heptanone"))
+label1 <- data.frame(x=c(Inf, 3), y=c(0.05, 0.02),
+                     text=c("0.1 mV", '1 s'),
+                     PulseDuration= factor(150, levels = c(15, 30, 150)),
+                     Group=factor(c("Lug", 'Mt Burns'), levels = creeks),
+                     Odor=c("2-heptanone"))
+
+fw_text <- data.frame(x = 0, y = 0.02,
+                      text=c("full-winged"),
+                      PulseDuration= factor(15, levels = c(15, 30, 150)),
+                      Group=factor("Mt Burns", levels = creeks),
+                      Odor=c("2-heptanone"))
+wr_text <- data.frame(x = 0, y = 0.05,
+                      text=c("wing-reduced"),
+                      PulseDuration = factor(15, levels = c(15, 30, 150)),
+                      Group=factor("Mt Burns", levels = creeks),
+                      Odor=c("2-heptanone"))
+
+pulseDur_text <- data.frame(x = 0.25, y = 0.1,
+                      text=c("15 ms", "30 ms", "150 ms"),
+                      PulseDuration = factor(c("15", "30", "150"), levels = c("15", "30", "150")),
+                      Group=factor("Black Jacks", levels = creeks),
+                      Odor=c("2-heptanone"))
+
+data_long$Group <- factor(data_long$Group, levels = c(creeks, 'Bee', 'PID'))
+
+
+col <- c(fw_l, wr_l)
+
+Fen <- ggplot(data_long %>% filter(Group %in% creeks), aes(x=time, y=Voltage)) + 
+  facet_grid(PulseDuration ~ Group, switch = "y") +
+  geom_rect(mapping=aes(xmin=0, xmax=as.numeric(as.character(PulseDuration))/1000,
+                        ymin=-Inf, ymax=Inf),
+            fill="grey90", alpha=1) +
+  scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
+  geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
+  geom_line(aes(colour = Winged)) +
+  theme_bw() +
+  labs(y = ' ') +
+  theme(legend.position = "none",
+        panel.grid.major = element_blank(),
+        panel.grid.minor = element_blank(),
+        axis.line = element_blank(),
+        axis.ticks = element_blank(),
+        axis.title.x = element_blank(),
+        panel.border = element_blank(),
+        plot.margin = unit(c(5.5, 5.5, 5.5, 5.5), "points"),
+        strip.background = element_rect(colour="NA", fill="NA"),
+        text = element_text(size=18),
+        aspect.ratio = ar * 3,
+        strip.text.y.left = element_blank(),
+        strip.text.x = element_text(size = 20),
+        axis.text=element_blank()) +
+  geom_segment(data=data_segm1,
+               aes(x=x, y=y, yend=yend, xend=xend)) +
+  geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
+  geom_text(data = wr_text, aes(x=x, y=y, label=text), size=6, hjust=0, color=wr_l) +
+  geom_text(data = fw_text, aes(x=x, y=y, label=text), size=6, hjust=0, color=fw_l) +
+  geom_text(data = pulseDur_text, aes(x=x, y=y, label=text), size=5, hjust=0, color='black') +
+  # geom_text(data = Diffs, aes(x=0.4, y=Max + 0.025, label=Star), size=geom.text.size, hjust=0.5) +
+  coord_cartesian(expand = TRUE,
+                  clip = 'off')
+Fen
+
+col <- c('blue')
+data_segm1 <- data.frame(y=0, yend=1, x=Inf, xend=Inf,
+                         PulseDuration = factor(150, levels = c(15, 30, 150)),
+                         Group=factor("Bee"),
+                         Odor=c("2-heptanone"))
+label1 <- data.frame(x=Inf, y=0.5,
+                     text=c("1 mV"),
+                     PulseDuration= factor(150, levels = c(15, 30, 150)),
+                     Group=factor("Bee"),
+                     Odor=c("2-heptanone"))
+
+
+Bee <- ggplot(data_long %>% filter(Group %in% 'Bee'), aes(x=time, y=Voltage)) + 
+  facet_grid(PulseDuration ~ Group, switch = "y") +
+  geom_rect(mapping=aes(xmin=0, xmax=as.numeric(as.character(PulseDuration))/1000,
+                        ymin=-Inf, ymax=Inf),
+            fill="grey90", alpha=1) +
+  scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
+  geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
+  geom_line(aes(colour = Winged)) +
+  theme_bw() +
+  theme(legend.position = "none",
+        panel.grid.major = element_blank(),
+        panel.grid.minor = element_blank(),
+        axis.line = element_blank(),
+        axis.title = element_blank(),
+        plot.margin = unit(c(5.5, 5.5, 5.5, 5.5), "points"),
+        strip.background = element_rect(colour="NA", fill="NA"),
+        panel.border = element_blank(),
+        text = element_text(size=18),
+        axis.ticks = element_blank(),
+        aspect.ratio = ar * 3,
+        strip.text.y = element_blank(),
+        strip.text.x = element_text(size = 20),
+        axis.text=element_blank()) +
+  geom_segment(data=data_segm1,
+               aes(x=x, y=y, yend=yend, xend=xend)) +
+  geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
+  # geom_text(data = Diffs, aes(x=0.4, y=Max + 0.025, label=Star), size=geom.text.size, hjust=0.5) +
+  coord_cartesian(expand = TRUE,
+                  clip = 'off')
+Bee
+
+data_segm1 <- data.frame(y=0, yend=3, x=Inf, xend=Inf,
+                         PulseDuration = factor(150, levels = c(15, 30, 150)),
+                         Group=factor("PID"),
+                         Odor=c("2-heptanone"))
+label1 <- data.frame(x=Inf, y=1.5,
+                     text=c("3 V"),
+                     PulseDuration= factor(150, levels = c(15, 30, 150)),
+                     Group=factor("PID"),
+                     Odor=c("2-heptanone"))
+
+PID <- ggplot(data_long %>% filter(Group %in% 'PID'), aes(x=time, y=Voltage)) + 
+  facet_grid(PulseDuration ~ Group, switch = "y") +
+  geom_rect(mapping=aes(xmin=0, xmax=as.numeric(as.character(PulseDuration))/1000,
+                        ymin=-Inf, ymax=Inf),
+            fill="grey90", alpha=1) +
+  scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
+  geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
+  geom_line(aes(colour = Winged)) +
+  theme_bw() +
+  theme(legend.position = "none",
+        panel.grid.major = element_blank(),
+        panel.grid.minor = element_blank(),
+        axis.line = element_blank(),
+        axis.title = element_blank(),
+        plot.margin = unit(c(5.5, 5.5, 5.5, 5.5), "points"),
+        strip.background = element_rect(colour="NA", fill="NA"),
+        panel.border = element_blank(),
+        text = element_text(size=18),
+        axis.ticks = element_blank(),
+        aspect.ratio = ar * 3,
+        strip.text.y = element_blank(),
+        strip.text.x = element_text(size = 20),
+        axis.text=element_blank()) +
+  geom_segment(data=data_segm1,
+               aes(x=x, y=y, yend=yend, xend=xend)) +
+  geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
+  # geom_text(data = Diffs, aes(x=0.4, y=Max + 0.025, label=Star), size=geom.text.size, hjust=0.5) +
+  coord_cartesian(expand = TRUE,
+                  clip = 'off')
+PID
+
+
+FigS1A <- egg::ggarrange(Fen, Bee, PID,
+                         heights = 5, nrow = 1, ncol=3, widths = c(5, 1, 1), debug = FALSE)
+
+########################
+
 
 #### Function to plot parameters for Response Dynamic:
 
@@ -1126,43 +1332,31 @@ Times <- data.frame("Odor"=times_tst_on$Var1, "on"=times_tst_on$on, "off"=times_
 
 Times <- Times[!(Times$Odor=="Blank"), ]
 
-plot_freq <- function(data, group, times, type="mean") {
-  
+plot_freq <- function(data, group, times) {
   
   data <- data %>%
-    filter(PulseDuration==50,
-           Odor!="Blank",
-           Odor!="Stonefly")
+    filter(PulseDuration == 50,
+           Odor == "2-heptanone")
   
   data <- droplevels(data)
   
   d_mean <- data %>%
     filter(Group %in% group) %>%
-    group_by(Winged, Odor, PulseDuration) %>%
-    summarise_at(.funs = mean, .vars=paste0("t", 1:5000))
+    group_by(Group, Winged, Odor, PulseDuration) %>%
+    summarise_at(.funs = median, .vars=paste0("t", 1:5000))
   
   # calculate the standard error 
   d_se <- data %>%
     filter(Group %in% group) %>%
-    group_by(Winged, Odor, PulseDuration) %>%
+    group_by(Group, Winged, Odor, PulseDuration) %>%
     summarise_at(.funs = se, .vars=paste0("t", 1:5000))
   
-  d_lower_quan <- data %>%
-    filter(Group %in% group) %>%
-    group_by(Winged, Odor, PulseDuration) %>%
-    summarise_at(.funs = quantile, probs=0.25, .vars=paste0("t", 1:5000))
-  
-  d_upper_quan <- data %>%
-    filter(Group %in% group) %>%
-    group_by(Winged, Odor, PulseDuration) %>%
-    summarise_at(.funs = quantile, probs=0.75, .vars=paste0("t", 1:5000))
-  
   # calculate sample size (n) 
   d_n <- data %>%
     filter(Group %in% group) %>%
-    group_by(Winged, Odor, PulseDuration) %>%
+    group_by(Group, Winged, Odor, PulseDuration) %>%
     summarise(N=n()) %>%
-    group_by(Winged) %>%
+    group_by(Group, Winged) %>%
     summarise(N=unique(N))
   d_n <- unite(d_n, col = "Winged", sep = " (n = ")
   d_n$t <- ")"
@@ -1173,75 +1367,154 @@ plot_freq <- function(data, group, times, type="mean") {
   
   
   data_se <- gather(d_se, time, se, t1:t5000, factor_key=TRUE)
-  data_upper <- gather(d_upper_quan, time, upp, t1:t5000, factor_key=TRUE)
-  data_lower <- gather(d_lower_quan, time, low, t1:t5000, factor_key=TRUE)
-  
-  
-  d_all$Odor <- factor(d_all$Odor, levels = c("2-heptanone", "1-octanol", 
-                                              "Propionic acid", "2-butanone"),
-                       labels = c("2-heptanone", "1-octanol", 
-                                  "Propionic acid", "2-butanone"))
-  times <- times[!times$Odor=="Stonefly", ]
-  times <- times[!times$Odor=="Blank", ]
+
   
+  d_all$Odor <- factor(d_all$Odor)
+  times <- times[times$Odor == "2-heptanone", ]
   
   data_long <- gather(d_all, time, Voltage, t1:t5000, factor_key=TRUE)
   data_long$time <- (as.numeric(gsub("t", "", data_long$time)) - 200)/1000
   data_long$SE <- data_se$se
-  data_long$upper <- data_upper$upp
-  data_long$lower <- data_lower$low
-  
-  
-  data_long$Winged <- factor(data_long$Winged,
-                             labels = t(d_n))
-  
+
   data_long$PulseDuration <- factor(data_long$PulseDuration, levels = c("50"), labels = c("10 Hz"))
   
   
-  if ( "Bee" %in% group | "Fenestrata-Bee" %in% group) {
-    col <- c("#000000")
-    
-  } else {
-    col <- c("#FF00FF", "#00DC00")
-  }
-  print(d_n)
+  col <- c("grey", "#00DC00")
+  
+  data_segm1 <- data.frame(y=c(0, -0.02), yend=c(0.05, -0.02), x=c(Inf, 1), xend=c(Inf, 2),
+                           PulseDuration = factor(150, levels = c(15, 30, 150)),
+                           Group=factor(c("Six Mile", 'Mt Burns')),
+                           Odor=c("2-heptanone"))
+  label1 <- data.frame(x=c(Inf, 2), y=c(0.035, -0.035),
+                       text=c("0.05 mV", '1 s'),
+                       PulseDuration= factor(150, levels = c(15, 30, 150)),
+                       Group=factor(c("Six Mile", 'Mt Burns')),
+                       Odor=c("2-heptanone"))
   
+  ## For Lug, we stimultaed with a different stimulus protocol were the stimulation
+  ## with several pulses in a specific frequency was not 10 Hz and the total length
+  ## o f stimulation was less than 3 seconds. We exclude these data here.
+  dat <- data_long %>% filter(Group == creeks, Group != 'Lug')
+  dat$Group <- factor(dat$Group, levels = c(creeks))
   
-  if (type=="mean") {
-    ggplot(data_long, aes(x=time, y=Voltage)) + 
+  W <- ggplot(dat, aes(x=time, y=Voltage)) + 
       geom_rect(data=times, inherit.aes=FALSE, fill = "grey95",
-                aes(xmin=on, xmax=off, ymin=-Inf,ymax=Inf)) +
-      geom_ribbon(aes(ymin=Voltage + SE, ymax=Voltage - SE, fill=Winged), color=NA, alpha=0.3) +
+                aes(xmin=on, xmax=off, ymin=-Inf, ymax=Inf)) +
+      # geom_ribbon(aes(ymin=Voltage + SE, ymax=Voltage - SE, fill=Winged), color=NA, alpha=0.3) +
       geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
       geom_line(aes(colour = Winged)) +
       scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
-      facet_grid(Odor ~ Winged, scales = "free_y") +
+      facet_grid(. ~ Group, scales = "free_y", drop = FALSE) +
       scale_y_continuous(breaks = pretty_breaks(n = 3)) +
       theme_bw() +
-      labs(x="Time (s)", y="Voltage (mV)") +
-      guides(fill = FALSE) +
+    labs(y = " ") +
+      guides(fill = 'none') +
+    geom_segment(data=data_segm1,
+                 aes(x=x, y=y, yend=yend, xend=xend)) +
+    geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
+      coord_cartesian(xlim =c(0, 3.5)) +
       theme(panel.grid.major = element_blank(),
             panel.grid.minor = element_blank(),
             legend.position = "none",
             strip.background = element_rect(colour="NA", fill="NA"),
-            panel.border = element_blank(),
+            strip.text.y = element_blank(),
             legend.title = element_blank(),
-            axis.line = element_line("black"),
-            strip.text.y = element_text(angle = 0, hjust = 0),
-            text = element_text(size=15))  } else {
-              
-            }
+            axis.line = element_blank(),
+            axis.ticks = element_blank(),
+            axis.title.x = element_blank(),
+            axis.title.y = element_text(size = 35),
+            panel.border = element_blank(),
+            axis.text=element_blank(),
+            strip.text = element_blank(),
+            text = element_text(size=18)) 
+  W
+  
+  col <- c("blue")
+  data_segm1 <- data.frame(y=0, yend=1, x=Inf, xend=Inf,
+                           PulseDuration = factor(150, levels = c(15, 30, 150)),
+                           Group=factor("Bee"),
+                           Odor=c("2-heptanone"))
+  label1 <- data.frame(x=Inf, y=0.5,
+                       text=c("1 mV"),
+                       PulseDuration= factor(150, levels = c(15, 30, 150)),
+                       Group=factor("Bee"),
+                       Odor=c("2-heptanone"))
+  
+  B <-  ggplot(data_long[data_long$Group %in% 'Bee', ], aes(x=time, y=Voltage)) + 
+    geom_rect(data=times, inherit.aes=FALSE, fill = "grey95",
+              aes(xmin=on, xmax=off, ymin=-Inf,ymax=Inf)) +
+    # geom_ribbon(aes(ymin=Voltage + SE, ymax=Voltage - SE, fill=Winged), color=NA, alpha=0.3) +
+    geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
+    geom_line(aes(colour = Winged)) +
+    scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
+    facet_grid(Odor ~ Group, scales = "free_y") +
+    scale_y_continuous(breaks = pretty_breaks(n = 3)) +
+    theme_bw() +
+    coord_cartesian(xlim =c(0, 3.5)) +
+    # guides(fill = 'none') +
+    geom_segment(data=data_segm1,
+                 aes(x=x, y=y, yend=yend, xend=xend)) +
+    geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
+    theme(panel.grid.major = element_blank(),
+          panel.grid.minor = element_blank(),
+          legend.position = "none",
+          strip.background = element_rect(colour="NA", fill="NA"),
+          panel.border = element_blank(),
+          legend.title = element_blank(),
+          axis.line = element_blank(),
+          strip.text = element_blank(),
+          axis.text=element_blank(),
+          axis.ticks = element_blank(),
+          axis.title = element_blank(),
+          text = element_blank()) 
+  B
+  
+  
+  data_segm1 <- data.frame(y=0, yend=1, x=3, xend=3,
+                           PulseDuration = factor(150, levels = c(15, 30, 150)),
+                           Group=factor("PID"),
+                           Odor=c("2-heptanone"))
+  label1 <- data.frame(x=3, y=0.5,
+                       text=c("1 V"),
+                       PulseDuration= factor(150, levels = c(15, 30, 150)),
+                       Group=factor("PID"),
+                       Odor=c("2-heptanone"))
+  
+  P <- ggplot(data_long[data_long$Group %in% 'PID', ], aes(x=time, y=Voltage)) + 
+    geom_rect(data = times, inherit.aes = FALSE, fill = "grey95",
+              aes(xmin = on, xmax=off, ymin = -Inf, ymax = Inf)) +
+    # geom_ribbon(aes(ymin=Voltage + SE, ymax=Voltage - SE, fill=Winged), color=NA, alpha=0.3) +
+    geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
+    geom_line(aes(colour = Winged)) +
+    scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
+    facet_grid(Odor ~ Group, scales = "free_y") +
+    scale_y_continuous(breaks = pretty_breaks(n = 3)) +
+    theme_bw() +
+    coord_cartesian(xlim =c(0, 3.5)) +
+    geom_segment(data=data_segm1,
+                 aes(x=x, y=y, yend=yend, xend=xend)) +
+    geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
+    # guides(fill = 'none') +
+    theme(panel.grid.major = element_blank(),
+          panel.grid.minor = element_blank(),
+          legend.position = "none",
+          strip.background = element_rect(colour="NA", fill="NA"),
+          panel.border = element_blank(),
+          legend.title = element_blank(),
+          axis.ticks = element_blank(),
+          axis.text=element_blank(),
+          axis.line = element_blank(),
+          axis.title = element_blank(),
+          strip.text = element_blank(),
+          text = element_blank()) 
+  P
   
+  list(W = W, B = B, P = P)
   
 }
 
-freq_lug <- plot_freq(data = data2, group = "Lug", times = Times, type="mean")
-freq_sixmile <- plot_freq(data2, "Six Mile", Times, type="mean")
-freq_blackjack <- plot_freq(data2, "Black Jacks", Times, type="mean")
-freq_whisky <- plot_freq(data2, "Whiskey", Times, type="mean")
-freq_burns <- plot_freq(data2, "Mt Burns", Times, type="mean")
-freq_bee <- plot_freq(data2, "Bee", Times, type="mean")
-freq_pid <- plot_freq(data2, "PID", Times, type="mean")
+
+FigS1C_new <- plot_freq(data2, c(creeks, 'Bee', 'PID'), Times)
 
 
 ## Power spectrum density analysis ####
@@ -1340,7 +1613,13 @@ plot_freq <- function(data) {
   
   odors <- data.frame("Odor"=levels(data$Odor))
   
-  W <- ggplot(data %>% filter(Group %in% creeks),
+  ## For Lug, we stimultaed with a different stimulus protocol were the stimulation
+  ## with several pulses in a specific frequency was not 10 Hz and the total length
+  ## o f stimulation was less than 3 seconds. We exclude these data here.
+  dat <- data %>% filter(Group %in% creeks & Group != 'Lug')
+  dat$Group <- factor(dat$Group, levels = c(creeks))
+  
+  W <- ggplot(dat,
               aes(x=x, y=log10(y))) + 
     geom_rect(mapping=aes(xmin=9.5, xmax=10.5,
                           ymin=-Inf, ymax=Inf),
@@ -1348,7 +1627,7 @@ plot_freq <- function(data) {
     geom_ribbon(aes(ymin=log10(y + se), ymax=log10(y - se), fill=Winged),
                 color=NA, alpha=0.3) +
     geom_line(aes(colour = Winged)) +
-    facet_grid(. ~ Group) +
+    facet_grid(. ~ Group, drop = FALSE) +
     scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
     theme_bw() +
     scale_y_continuous(breaks = pretty_breaks(n=3)) +
@@ -1422,17 +1701,21 @@ plot_freq <- function(data) {
 p5 <- plot_freq(data = spec_sel); p5
 #create common x and y labels
 
-FigS1B <- egg::ggarrange(p5$traceW, p5$traceB, p5$traceP,
+FigS1D <- egg::ggarrange(p5$traceW, p5$traceB, p5$traceP,
                          bottom = textGrob("Frequency (Hz)", gp = gpar(fontsize=16), vjust = 0.5),
                          left = textGrob(expression(paste("Power (", mV^2/Hz, ") [log10]")),
                                          gp = gpar(fontsize=18), vjust = 0.5, rot = 90),
                          ncol = 3, nrow=1,
                          label.args = list(gp=gpar(fontface="bold", fontsize=18), hjust=0, vjust=1),
                          widths = c(5,1,1),
-                         debug=FALSE)
-FigS1B
-
+                         debug = FALSE)
+FigS1D
 
+FigS1C <- egg::ggarrange(FigS1C_new$W, FigS1C_new$B, FigS1C_new$P,
+                         ncol = 3, nrow=1, heights = 1,
+                         widths = c(5, 1, 1),
+                         debug = FALSE)
+FigS1C
 
 Fig1AB <- cowplot::plot_grid(Fig1A, Fig1B, 
                              align = "hv", axis = "trbl",
@@ -1446,7 +1729,6 @@ cowplot::plot_grid(Fig1AB, Fig1C, nrow = 2, rel_heights = c(1, 1.3),
 
 
 ggsave("Figure1.pdf", path = fig_path, width = 12, height = 10)
-ggsave("Figure1.jpg", path = fig_path, width = 12, height = 10)
 
 
 ### Fig. 2
@@ -1455,12 +1737,10 @@ cowplot::plot_grid(Fig2A, Fig2B, nrow = 2,
 
 
 ggsave("Figure2.pdf", path = fig_path, width = 9, height = 12)
-ggsave("Figure2.jpg", path = fig_path, width = 9, height = 12)
 
 
 ### Fig. S1
-cowplot::plot_grid(FigS1A, NULL, FigS1B, nrow = 3, rel_heights = c(1.3,0.1, 1),
-                   labels = c('A', NULL, 'B'), label_size = 18)
+cowplot::plot_grid(FigS1A, FigS1B, NULL, FigS1C, FigS1D, nrow = 5, rel_heights = c(1.2, 1, 0.05, 0.3, 0.8),
+                   labels = c('A', 'B', '', 'C', 'D'), label_size = 18)
 
-ggsave("FigureS1.pdf", path = fig_path, width = 14, height = 8)
-ggsave("FigureS1.jpg", path = fig_path, width = 14, height = 8)
+ggsave("FigureS1.pdf", path = fig_path, width = 14, height = 14)