diff --git a/four_bumps_analysis.R b/four_bumps_analysis.R
new file mode 100644
index 0000000..7a01c79
--- /dev/null
+++ b/four_bumps_analysis.R
@@ -0,0 +1,137 @@
+#install package to read matrices
+install.packages("R.matlab")
+
+#load R.matlab package and load the data
+library(R.matlab)
+bumps_data <- readMat("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/four_bumps_bsl.mat")
+peak1<- t(bumps_data$col.peaks[,1])
+peak2<- t(bumps_data$col.peaks[,2])
+peak3<- t(bumps_data$col.peaks[,3])
+peak4<- t(bumps_data$col.peaks[,4])
+
+#file_names contains both the filenames and the layer
+file_names_data <- read.csv("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/filenames.csv", header = TRUE,sep = ',')
+
+
+col_bumps_data <- t(cbind(peak1,peak2,peak3,peak4))
+
+#Create bump index column
+bmp_name <- c("peak1","peak2","peak3","peak4")
+bumpnb <- matrix(rep(bmp_name, each =74), ncol=4)
+col_bumpnb <- t(cbind(t(bumpnb[,1]),t(bumpnb[,2]),t(bumpnb[,3]),t(bumpnb[,4])))
+
+#create channel index
+channel_idx = t(t(rep(1:74, times =4)))
+
+#data frame
+
+org_data = data.frame(channel_idx,col_bumpnb, col_bumps_data, file_names_data)
+
+#anova
+install.packages("tidyverse")
+library(tidyverse)
+anovaPeaks = aov( col_bumps_data ~ col_bumpnb + Error(channel_idx), org_data)
+summary(anovaPeaks)
+
+#lmer
+library(lmerTest)
+library(car)
+library(lme4)
+linearPeaks = lmer(col_bumps_data ~ col_bumpnb +(1|channel_idx), org_data)
+summary(linearPeaks)
+  
+
+
+#posthoc
+library(lsmeans)
+lsmip(linearPeaks, col_bumpnb) #got an error in this one 
+
+   #other way for posthocs
+library(multcompView)
+CLD(lsmeans(linearPeaks, c("col_bumpnb")))
+
+   #more detailed way for posthoc
+CLD(emmeans(linearPeaks, ~col_bumpnb ))
+library(multcomp)
+posthoc = emmeans(linearPeaks, "col_bumpnb")
+summary(as.glht(update(pairs(posthoc), by = NULL)), test = adjusted("free")) #got an error
+
+    #more simple method #couldn't get it to work, need the package
+posthoc(linearPeaks, alpha = 0.05,test = "pairwise", base = 1, alternative = "two.sided",
+    statistic = "Wald",
+    padj.method = "free")
+
+## lmer accounting for the layers
+linearPeaks = lmer(col_bumps_data ~ col_bumpnb*layer +(1|channel_idx), org_data)
+summary(linearPeaks)
+CLD(lsmeans(linearPeaks, c("col_bumpnb", "layer")))
+CLD(emmeans(linearPeaks, ~col_bumpnb|layer ))
+library(multcomp)
+posthoc = emmeans(linearPeaks, "col_bumpnb", by = "layer")
+summary(as.glht(update(pairs(posthoc), by = NULL)))#, test = adjusted("free")) #got an error
+
+
+
+# Plot the data as a function of layers
+org_data$layer <- factor(org_data$layer, 
+                                 levels = c("P", "M", "K"),
+                                 labels = c("P", "M", "K"))
+org_data$col_bumpnb <- factor(org_data$col_bumpnb, 
+                                 levels = c("peak1", "peak2","peak3","peak4"),
+                                 labels = c("peak1", "peak2","peak3","peak4"))
+library("ggpubr")
+library("colorspace")
+png('all_layers_plot.png')
+ggboxplot(org_data, x = "col_bumpnb", y = "col_bumps_data", color = "layer")
+          palette = c("#00AFBB", "#E7B800","#E495A5","#ABB065")
+dev.off()
+
+##other plots
+   # organizing each leayer's data
+p_data = org_data$col_bumps_data[org_data$layer=="P"]
+peak_nb = org_data$col_bumpnb[org_data$layer=="P"]
+p_idx = org_data$channel_idx[org_data$layer=="P"]
+p_lay = data.frame(peak_nb, p_data, p_idx)
+
+m_data = org_data$col_bumps_data[org_data$layer=="M"]
+mpeak_nb = org_data$col_bumpnb[org_data$layer=="M"]
+m_idx = org_data$channel_idx[org_data$layer=="M"]
+m_lay = data.frame(mpeak_nb, m_data, m_idx)
+
+k_data = org_data$col_bumps_data[org_data$layer=="K"]
+kpeak_nb = org_data$col_bumpnb[org_data$layer=="K"]
+k_idx = org_data$channel_idx[org_data$layer=="K"]
+k_lay = data.frame(kpeak_nb, k_data, k_idx)
+
+   #comparisons
+my_comparisons = list(c("peak1", "peak2"), c("peak1", "peak3"), c("peak1","peak4"), c("peak2","peak3"), c("peak2", "peak4"), c("peak3","peak4"))
+
+
+png('p_layer_plot.png')
+dodge = position_dodge(0.6)
+p_lay %>%
+ggplot(aes(x = peak_nb, y = p_data, na.rm = TRUE),add = "dotplot")+
+ geom_boxplot(width = .3, position = dodge, alpha = .5) +
+  #scale_fill_viridis(discrete = TRUE) +
+scale_y_continuous(limits = c(0, 150), breaks = seq(0,150, by=10)) +
+stat_compare_means(comparisons = my_comparisons, data = p_lay, label = "p.signif", hide.ns = FALSE, method = "t.test")+
+ theme_bw() +
+  theme(
+    plot.title   = element_text(color = "black", size = 15, face = "bold", hjust=0.5),
+    axis.title.x = element_text(color = "black", size = 13),
+    axis.title.y = element_text(color = "black", size = 13),
+    panel.border = element_blank(),
+    axis.text    = element_text(size  = 9), 
+    #legend.position = "none",
+    #strip.background = element_rect(colour="gray94", fill="gray94")
+  )+
+
+          labs(
+          title = "P layer peaks",
+          x = "Peak Number",
+          y = "Spike rate (spike/s)" )
+          #palette = c("#00AFBB", "#E7B800","#E495A5","#ABB065")
+dev.off()
+
+ 
+
diff --git a/individual_channels_analysis.R b/individual_channels_analysis.R
new file mode 100644
index 0000000..b2865f0
--- /dev/null
+++ b/individual_channels_analysis.R
@@ -0,0 +1,116 @@
+library(R.matlab)
+library(lmerTest)
+library(car)
+library(lme4)
+
+
+filenames <-list.files("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels/", pattern=NULL, all.files=FALSE,
+    full.names=FALSE)
+#peaknames <- t(matrix(rep(labels, each =1), ncol=74))
+
+## treat each data file one by one
+
+bumps_data <- readMat("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels/160517_I_p01_uclust1_cinterocdrft_stab_fft_sigatatDE50_NDE0.mat")
+
+peak1<- t(bumps_data$all.pks[1,])
+peak2<- t(bumps_data$all.pks[2,])
+peak3<- t(bumps_data$all.pks[3,])
+peak4<- t(bumps_data$all.pks[4,])
+col_bumps_data <- t(cbind(peak1,peak2,peak3,peak4))
+
+#lbpksdat <- data.frame(col_bumps_data)
+#colnames(lbpksdat) <- labels
+
+labels <- c("P1","P2","P3","P4")
+bumpnb <- matrix(rep(labels, each =length(bumps_data$all.pks[1,])), ncol=4)
+p_label <- t(cbind(t(bumpnb[,1]),t(bumpnb[,2]),t(bumpnb[,3]),t(bumpnb[,4])))
+
+#create table
+org_data = data.frame(p_label, col_bumps_data)
+
+#fit linear model
+linearPeaks = lm(col_bumps_data ~ p_label, org_data)
+summary(linearPeaks)
+
+##
+## create a for loop to treat all the data at once
+filenames <-list.files("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels/", pattern=NULL, all.files=FALSE,
+data_filename_list <- Sys.glob("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels/*.mat")
+
+require(pbkrtest)
+pvals <- array(0, c(length(data_filename_list),4))
+tvals <- 0
+df.KR <- 0
+p.KR <- 0
+for(i in 1:length(data_filename_list)){
+
+filename <- filenames[i]
+bumps_data <- readMat(data_filename_list[i])
+labels <- c("P1","P2","P3","P4")
+peak1<- t(bumps_data$all.pks[1,])
+peak2<- t(bumps_data$all.pks[2,])
+peak3<- t(bumps_data$all.pks[3,])
+peak4<- t(bumps_data$all.pks[4,])
+col_bumps_data <- t(cbind(peak1,peak2,peak3,peak4))
+
+#create labels
+labels <- c("P1","P2","P3","P4")
+bumpnb <- matrix(rep(labels, each =length(bumps_data$all.pks[1,])), ncol=4)
+p_label <- t(cbind(t(bumpnb[,1]),t(bumpnb[,2]),t(bumpnb[,3]),t(bumpnb[,4])))
+
+#create trial index
+trial_idx = t(t(rep(1:length(bumps_data$all.pks[1,]), times =4)))
+
+#create table
+org_data = data.frame(trial_idx, p_label, col_bumps_data)
+
+##fit linear model
+linearPeaks = lmer(col_bumps_data ~ p_label +(1|trial_idx), org_data)
+linsummary <- summary(linearPeaks)
+
+##work in progress to extract values from the table
+#sink("linsummary.txt")
+#linsummary
+#sink()
+#linsummary.string = readChar("linsummary.txt",20000)
+#splitstring = strsplit(linsummary.string,"\r\n")
+#interstring = strsplit(splitstring[[1]][19]," ")
+#Lines contain TABs, so the p-value is the eighth element along:
+#pvalue = as.numeric(interstring[[1]][8])
+#store t-stats and pvalues in a vector
+
+coefs <- data.frame(coef(summary(linearPeaks)))
+df.KR <- get_ddf_Lb(linearPeaks, fixef(linearPeaks))
+p.KR <- 2 * (1 - pt(abs(coefs$t.value), df.KR))
+pvals[i,] <- p.KR
+
+#for(j in 1:4){ 
+#if you wanna consider the distribution normal (not recommended if small sample)
+#Vcov <- vcov(linearPeaks, useScale = FALSE) 
+#betas <- fixef(linearPeaks) 
+#se <- sqrt(diag(Vcov)) 
+#zval <- betas / se 
+#pvals[i,] <- 2 * pnorm(abs(zval), lower.tail = FALSE)
+#pvals[i,] <- 2 * (1 - pnorm(abs(coefs$t.value)))
+#}
+
+##save output in a mat file
+#filename2 <- paste("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels/lmer_results",filename, ".mat", sep = "")
+#writeMat(filename2, table = linsummary)
+
+print(p.KR) 
+#print(filename)
+print(linsummary)
+}
+#print(tvals)
+#print(df.KR)
+#print(p.KR)
+print(pvals)
+#save pvalues as a .mat file
+filename3 <- paste("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels/lmer_results", ".mat", sep = "")
+writeMat(filename3, pvalues = pvals)
+
+#save pvalues as a .csv file
+filename4 <- paste("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels/lmer_results", ".csv", sep = "")
+write.csv(pvals, file = filename4)
+
diff --git a/individual_channels_analysis_KRpvalcorr.R b/individual_channels_analysis_KRpvalcorr.R
new file mode 100644
index 0000000..eec6058
--- /dev/null
+++ b/individual_channels_analysis_KRpvalcorr.R
@@ -0,0 +1,77 @@
+library(R.matlab)
+library(lmerTest)
+library(car)
+library(lme4)
+
+require(pbkrtest)
+
+
+##
+## create a for loop to treat all the data at once
+filenames <-list.files("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels_peakadj2/", pattern=NULL, all.files=FALSE,  full.names=FALSE)
+data_filename_list <- Sys.glob("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels_peakadj2/*.mat")
+pvals <- array(0, c(length(data_filename_list),4))
+tvals <- 0
+df.KR <- 0
+p.KR <- 0
+for(i in 1:length(data_filename_list)){
+
+filename <- filenames[i]
+bumps_data <- readMat(data_filename_list[i])
+labels <- c("P1","P2","P3","P4")
+peak1<- t(bumps_data$all.pks[1,])
+peak2<- t(bumps_data$all.pks[2,])
+peak3<- t(bumps_data$all.pks[3,])
+peak4<- t(bumps_data$all.pks[4,])
+col_bumps_data <- t(cbind(peak1,peak2,peak3,peak4))
+
+#create labels
+labels <- c("P1","P2","P3","P4")
+bumpnb <- matrix(rep(labels, each =length(bumps_data$all.pks[1,])), ncol=4)
+p_label <- t(cbind(t(bumpnb[,1]),t(bumpnb[,2]),t(bumpnb[,3]),t(bumpnb[,4])))
+
+#create trial index
+trial_idx = t(t(rep(1:length(bumps_data$all.pks[1,]), times =4)))
+
+#create table
+org_data = data.frame(trial_idx, p_label, col_bumps_data)
+
+##fit linear model
+linearPeaks = lmer(col_bumps_data ~ p_label +(1|trial_idx), org_data)
+linsummary <- summary(linearPeaks)
+
+coefs <- data.frame(coef(summary(linearPeaks)))
+df.KR <- get_ddf_Lb(linearPeaks, fixef(linearPeaks))
+p.KR <- 2 * (1 - pt(abs(coefs$t.value), df.KR))
+pvals[i,] <- p.KR
+
+#for(j in 1:4){ 
+#if you wanna consider the distribution normal (not recommended if small sample)
+#Vcov <- vcov(linearPeaks, useScale = FALSE) 
+#betas <- fixef(linearPeaks) 
+#se <- sqrt(diag(Vcov)) 
+#zval <- betas / se 
+#pvals[i,] <- 2 * pnorm(abs(zval), lower.tail = FALSE)
+#pvals[i,] <- 2 * (1 - pnorm(abs(coefs$t.value)))
+#}
+
+##save output in a mat file
+#filename2 <- paste("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels/lmer_results",filename, ".mat", sep = "")
+#writeMat(filename2, table = linsummary)
+
+print(p.KR) 
+#print(filename)
+#print(linsummary)
+}
+#print(tvals)
+#print(df.KR)
+#print(p.KR)
+print(pvals)
+#save pvalues as a .mat file
+filename3 <- paste("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels_peakadj2/lmer_results/lmer_results", ".mat", sep = "")
+writeMat(filename3, pvalues = pvals)
+
+#save pvalues as a .csv file
+filename4 <- paste("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels_peakadj2/lmer_results/lmer_results", ".csv", sep = "")
+write.csv(pvals, file = filename4)
+
diff --git a/individual_channels_analysis_troughs_KRpvalcorr.R b/individual_channels_analysis_troughs_KRpvalcorr.R
new file mode 100644
index 0000000..1d0df5c
--- /dev/null
+++ b/individual_channels_analysis_troughs_KRpvalcorr.R
@@ -0,0 +1,75 @@
+library(R.matlab)
+library(lmerTest)
+library(car)
+library(lme4)
+
+require(pbkrtest)
+
+
+##
+## create a for loop to treat all the data at once
+filenames <-list.files("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels_troughadj/", pattern=NULL, all.files=FALSE,  full.names=FALSE)
+data_filename_list <- Sys.glob("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels_troughadj/*.mat")
+pvals <- array(0, c(length(data_filename_list),3))
+tvals <- 0
+df.KR <- 0
+p.KR <- 0
+for(i in 1:length(data_filename_list)){
+
+filename <- filenames[i]
+bumps_data <- readMat(data_filename_list[i])
+labels <- c("T1","T2","T3")
+trough1<- t(bumps_data$all.trghs[1,])
+trough2<- t(bumps_data$all.trghs[2,])
+trough3<- t(bumps_data$all.trghs[3,])
+col_bumps_data <- t(cbind(trough1,trough2,trough3))
+
+#create labels
+labels <- c("T1","T2","T3")
+bumpnb <- matrix(rep(labels, each =length(bumps_data$all.trghs[1,])), ncol=3)
+t_label <- t(cbind(t(bumpnb[,1]),t(bumpnb[,2]),t(bumpnb[,3])))
+
+#create trial index
+trial_idx = t(t(rep(1:length(bumps_data$all.trghs[1,]), times =3)))
+
+#create table
+org_data = data.frame(trial_idx, t_label, col_bumps_data)
+
+##fit linear model
+linearPeaks = lmer(col_bumps_data ~ t_label +(1|trial_idx), org_data)
+linsummary <- summary(linearPeaks)
+
+coefs <- data.frame(coef(summary(linearPeaks)))
+df.KR <- get_ddf_Lb(linearPeaks, fixef(linearPeaks))
+p.KR <- 2 * (1 - pt(abs(coefs$t.value), df.KR))
+pvals[i,] <- p.KR
+
+#for(j in 1:4){ 
+#if you wanna consider the distribution normal (not recommended if small sample)
+#Vcov <- vcov(linearPeaks, useScale = FALSE) 
+#betas <- fixef(linearPeaks) 
+#se <- sqrt(diag(Vcov)) 
+#zval <- betas / se 
+#pvals[i,] <- 2 * pnorm(abs(zval), lower.tail = FALSE)
+#pvals[i,] <- 2 * (1 - pnorm(abs(coefs$t.value)))
+#}
+
+##save output in a mat file
+#filename2 <- paste("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels/lmer_results",filename, ".mat", sep = "")
+#writeMat(filename2, table = linsummary)
+
+print(p.KR) 
+#print(filename)
+#print(linsummary)
+}
+#print(tvals)
+#print(df.KR)
+#print(p.KR)
+print(pvals)
+#save pvalues as a .mat file
+filename3 <- paste("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels_troughadj/lmer_results/lmer_results", ".mat", sep = "")
+writeMat(filename3, pvalues = pvals)
+
+#save pvalues as a .csv file
+filename4 <- paste("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels_troughadj/lmer_results/lmer_results", ".csv", sep = "")
+write.csv(pvals, file = filename4)
diff --git a/plots_peaks_power.R b/plots_peaks_power.R
new file mode 100644
index 0000000..f1138d0
--- /dev/null
+++ b/plots_peaks_power.R
@@ -0,0 +1,321 @@
+library(R.matlab)
+library("ggpubr")
+library("colorspace")
+bumps_data <- readMat("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/individual_channels_peakadj2/mean_peak1_peak4.mat")
+
+parts_data <- readMat("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/power_spectrum/part1_part2_norm_power.mat")
+
+#file_names contains both the filenames and the layer
+file_names_data <- read.csv("C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/filenames.csv", header = TRUE,sep = ',', na.strings="")
+
+peak1<- t(t(bumps_data$peak1.peak4[,1]))
+peak4<- t(t(bumps_data$peak1.peak4[,2]))
+
+
+
+#create channel index
+channel_idx = t(t(rep(1:74, times =1)))
+
+#data frame
+
+org_data = data.frame(channel_idx,peak1, peak4, file_names_data)
+
+
+library(tidyr)
+#no_na_data <- drop_na(org_data)
+
+row.has.na <- apply(org_data, 1, function(x){any(is.na(x))})
+org_data.filt <- org_data[!row.has.na,]
+
+##add normalized diff between peak 1 and peak 4 for colour gradient in plot multiplied by .045 for them to be between -1 and 1
+#org_data.filt$diff <- ((org_data.filt$peak1 - org_data.filt$peak4)-min(org_data.filt$peak1 - org_data.filt$peak4))/(max(org_data.filt$peak1 - org_data.filt$peak4)-min(org_data.filt$peak1 - org_data.filt$peak4))
+org_data.filt$diff <- ((org_data.filt$peak1 - org_data.filt$peak4))/(min(org_data.filt$peak4 - org_data.filt$peak1))
+#org_data.filt$diff <- (org_data.filt$peak4 - org_data.filt$peak1)*.045
+#max_val <- max(abs(org_data.filt$diff))
+
+
+## plot each cell class independently
+
+
+library(reshape2)
+library(ggplot2)
+
+## plot the P cell class
+peak1 = org_data.filt$peak1[org_data.filt$layer=="P"]
+peak4 = org_data.filt$peak4[org_data.filt$layer=="P"]
+p_idx =org_data.filt$channel_idx[org_data.filt$layer=="P"]
+p_diff = org_data.filt$diff[org_data.filt$layer == "P"]
+p_peakdata = data.frame(p_idx,peak1,peak4, p_diff)
+
+# convert from wide to long
+plot_dat <- melt(p_peakdata, id.var=c('p_idx','p_diff'))
+
+# plot
+svg('C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/plots/p_class_plot_peak1_peak4_normnegpos.svg')
+ggplot(plot_dat) + 
+  # simple lines
+  geom_line(aes(x=paste0('dots ', variable), y=value, group=p_idx, color=p_diff)) +
+  scale_colour_gradient(low="#FF3366")+
+  #scale_colour_gradientn(colours=heat.colors(50, alpha=.7, rev = TRUE))+
+  # box plots and jitter points, with modified x value
+  geom_violin(aes(x=variable, y=value,fill=variable), width=0.6, trim=F)+
+  scale_fill_manual(values=c("#FF3366", "#FF6699"))+
+  geom_boxplot(aes(x=variable, y=value),width=0.1) +
+  geom_jitter(alpha =.5, aes(x=paste0('dots ', variable), y=value, color=p_diff), position = position_jitter(width = .02)) +
+
+  # specify x value order
+ scale_x_discrete(limits=c('peak1', 'dots peak1', 'dots peak4',
+                           'peak4'))+
+scale_y_continuous(limits = c(0, 150), breaks = seq(0,150, by=10)) +
+ theme_bw() +
+  theme(
+    plot.title   = element_text(color = "black", size = 15, face = "bold", hjust=0.5),
+    axis.title.x = element_text(color = "black", size = 13),
+    axis.title.y = element_text(color = "black", size = 13),
+    panel.border = element_blank(),
+    axis.text    = element_text(size  = 9),
+    axis.line = element_line(colour = "black", size = 1) 
+     )+
+
+          labs(
+          title = "P cell class comparison: peak1 vs peak4 spiking activity",
+          x = "Peak Number",
+          y = "Spike rate (spikes/s)" )
+dev.off()
+
+## plot the M cells class
+peak1 = org_data.filt$peak1[org_data.filt$layer=="M"]
+peak4 = org_data.filt$peak4[org_data.filt$layer=="M"]
+p_idx =org_data.filt$channel_idx[org_data.filt$layer=="M"]
+p_diff = org_data.filt$diff[org_data.filt$layer == "M"]
+m_peakdata = data.frame(p_idx,peak1,peak4,p_diff)
+
+# convert from wide to long
+plot_dat <- melt(m_peakdata, id.var=c('p_idx','p_diff'))
+
+# plot
+svg('C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/plots/m_class_plot_peak1_peak4_normnegpos.svg')
+ggplot(plot_dat) + 
+  # simple lines
+  geom_line(aes(x=paste0('dots ', variable), y=value, group=p_idx, color=p_diff)) +
+  scale_colour_gradient(low="#333333")+
+  # box plots and jitter points, with modified x value
+  geom_violin(aes(x=variable, y=value,fill=variable), width=0.6, trim=F)+
+  scale_fill_manual(values=c("#555555", "#999999"))+
+  geom_boxplot(aes(x=variable, y=value),width=0.1) +
+  geom_jitter(alpha =.5, aes(x=paste0('dots ', variable), y=value, color=p_diff), position = position_jitter(width = .02)) +
+
+  # specify x value order
+ scale_x_discrete(limits=c('peak1', 'dots peak1', 'dots peak4',
+                           'peak4'))+
+scale_y_continuous(limits = c(-20, 150), breaks = seq(-20,150, by=10)) +
+ theme_bw() +
+  theme(
+    plot.title   = element_text(color = "black", size = 15, face = "bold", hjust=0.5),
+    axis.title.x = element_text(color = "black", size = 13),
+    axis.title.y = element_text(color = "black", size = 13),
+    panel.border = element_blank(),
+    axis.text    = element_text(size  = 9),
+    axis.line = element_line(colour = "black", size = 1) 
+     )+
+
+          labs(
+          title = "M cell class comparison: peak1 vs peak4 spiking activity",
+          x = "Peak Number",
+          y = "Spike rate (spikes/s)" )
+dev.off()
+
+
+## plot the K cells class
+peak1 = org_data.filt$peak1[org_data.filt$layer=="K"]
+peak4 = org_data.filt$peak4[org_data.filt$layer=="K"]
+p_idx =org_data.filt$channel_idx[org_data.filt$layer=="K"]
+p_diff = org_data.filt$diff[org_data.filt$layer == "K"]
+k_peakdata = data.frame(p_idx,peak1,peak4,p_diff)
+
+# convert from wide to long
+plot_dat <- melt(k_peakdata, id.var=c('p_idx','p_diff'))
+
+# plot
+svg('C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/plots/k_class_plot_peak1_peak4_5chan_normposneg.svg')
+ggplot(plot_dat) + 
+  # simple lines
+  geom_line(aes(x=paste0('dots ', variable), y=value, group=p_idx, color=p_diff)) +
+  scale_colour_gradient(low="#9CCC00")+
+  # box plots and jitter points, with modified x value
+  geom_violin(aes(x=variable, y=value,fill=variable), width=0.6, trim=F)+
+  scale_fill_manual(values=c("#9CCC00", "#CCFF00"))+
+  geom_boxplot(aes(x=variable, y=value),width=0.1) +
+  geom_jitter(alpha =.5, aes(x=paste0('dots ', variable), y=value, color=p_diff), position = position_jitter(width = .02)) +
+
+  # specify x value order
+ scale_x_discrete(limits=c('peak1', 'dots peak1', 'dots peak4',
+                           'peak4'))+
+scale_y_continuous(limits = c(0, 120), breaks = seq(0,120, by=10)) +
+ theme_bw() +
+  theme(
+    plot.title   = element_text(color = "black", size = 15, face = "bold", hjust=0.5),
+    axis.title.x = element_text(color = "black", size = 13),
+    axis.title.y = element_text(color = "black", size = 13),
+    panel.border = element_blank(),
+    axis.text    = element_text(size  = 9),
+    axis.line = element_line(colour = "black", size = 1) 
+     )+
+
+          labs(
+          title = "K cell class comparison: peak1 vs peak4 spiking activity",
+          x = "Peak Number",
+          y = "Spike rate (spikes/s)" )
+dev.off()
+
+
+#### power plots
+
+window1<- t(t(parts_data$parts[,1]))
+window2<- t(t(parts_data$parts[,2]))
+#create channel index
+channel_idx = t(t(rep(1:74, times =1)))
+
+#data frame
+
+org_data = data.frame(channel_idx,window1, window2, file_names_data)
+
+
+library(tidyr)
+#no_na_data <- drop_na(org_data)
+
+row.has.na <- apply(org_data, 1, function(x){any(is.na(x))})
+org_data.filt <- org_data[!row.has.na,]
+org_data.filt$diff <- ((org_data.filt$window1 - org_data.filt$window2))/(min(org_data.filt$window2 - org_data.filt$window1))
+
+## plot the P cell class
+window1 = org_data.filt$window1[org_data.filt$layer=="P"]
+window2 = org_data.filt$window2[org_data.filt$layer=="P"]
+p_idx =org_data.filt$channel_idx[org_data.filt$layer=="P"]
+p_diff = org_data.filt$diff[org_data.filt$layer == "P"]
+
+p_powerdata = data.frame(p_idx,window1,window2,p_diff)
+
+# convert from wide to long
+plot_dat <- melt(p_powerdata, id.var=c('p_idx','p_diff'))
+
+# plot
+svg('C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/power_spectrum/plots/p_class_plot_part1_part2_normposneg.svg')
+ggplot(plot_dat) + 
+  # simple lines
+  geom_line(aes(x=paste0('dots ', variable), y=value, group=p_idx, color=p_diff)) +
+  scale_colour_gradient(low="#FF3366")+
+  # box plots and jitter points, with modified x value
+  geom_violin(aes(x=variable, y=value,fill=variable), width=0.6, trim=F)+
+  scale_fill_manual(values=c("#FF3366", "#FF6699"))+
+  geom_boxplot(aes(x=variable, y=value),width=0.1) +
+  geom_jitter(alpha =.5, aes(x=paste0('dots ', variable), y=value, color=p_diff), position = position_jitter(width = .02)) +
+
+  # specify x value order
+ scale_x_discrete(limits=c('window1', 'dots window1', 'dots window2',
+                           'window2'))+
+scale_y_continuous(limits = c(0, 1), breaks = seq(0,1, by=.1)) +
+ theme_bw() +
+  theme(
+    plot.title   = element_text(color = "black", size = 15, face = "bold", hjust=0.5),
+    axis.title.x = element_text(color = "black", size = 13),
+    axis.title.y = element_text(color = "black", size = 13),
+    panel.border = element_blank(),
+    axis.text    = element_text(size  = 9),
+    axis.line = element_line(colour = "black", size = 1) 
+     )+
+
+          labs(
+          title = "P cell class comparison: window1 vs window2 mean power",
+          x = "Window Number",
+          y = "Power at 4Hz (normalized)" )
+dev.off()
+
+
+## plot the M cell class
+window1 = org_data.filt$window1[org_data.filt$layer=="M"]
+window2 = org_data.filt$window2[org_data.filt$layer=="M"]
+p_idx =org_data.filt$channel_idx[org_data.filt$layer=="M"]
+p_diff = org_data.filt$diff[org_data.filt$layer == "M"]
+m_powerdata = data.frame(p_idx,window1,window2,p_diff)
+
+# convert from wide to long
+plot_dat <- melt(m_powerdata, id.var=c('p_idx','p_diff'))
+
+# plot
+svg('C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/power_spectrum/plots/m_class_plot_part1_part2_normposneg.svg')
+ggplot(plot_dat) + 
+  # simple lines
+  geom_line(aes(x=paste0('dots ', variable), y=value, group=p_idx, color=p_diff)) +
+  scale_colour_gradient(low="#333333")+
+  # box plots and jitter points, with modified x value
+  geom_violin(aes(x=variable, y=value,fill=variable), width=0.6, trim=F)+
+  scale_fill_manual(values=c("#555555", "#999999"))+
+  geom_boxplot(aes(x=variable, y=value),width=0.1) +
+  geom_jitter(alpha =.5, aes(x=paste0('dots ', variable), y=value, color=p_diff), position = position_jitter(width = .02)) +
+
+  # specify x value order
+ scale_x_discrete(limits=c('window1', 'dots window1', 'dots window2',
+                           'window2'))+
+scale_y_continuous(limits = c(0, 1), breaks = seq(0,1, by=.1)) +
+ theme_bw() +
+  theme(
+    plot.title   = element_text(color = "black", size = 15, face = "bold", hjust=0.5),
+    axis.title.x = element_text(color = "black", size = 13),
+    axis.title.y = element_text(color = "black", size = 13),
+    panel.border = element_blank(),
+    axis.text    = element_text(size  = 9),
+    axis.line = element_line(colour = "black", size = 1) 
+     )+
+
+          labs(
+          title = "M cell class comparison: window1 vs window2 mean power",
+          x = "Window Number",
+          y = "Power at 4Hz (normalized)" )
+dev.off()
+
+
+## plot the K cell class
+window1 = org_data.filt$window1[org_data.filt$layer=="K"]
+window2 = org_data.filt$window2[org_data.filt$layer=="K"]
+p_idx =org_data.filt$channel_idx[org_data.filt$layer=="K"]
+p_diff = org_data.filt$diff[org_data.filt$layer == "K"]
+k_powerdata = data.frame(p_idx,window1,window2,p_diff)
+
+# convert from wide to long
+plot_dat <- melt(k_powerdata, id.var=c('p_idx','p_diff'))
+
+# plot
+svg('C:/Users/maier/Documents/LGN_data/single_units/inverted_power_channels/good_single_units_data_4bumps_more/power_spectrum/plots/k_class_plot_part1_part2_5chan_normposdiff.svg')
+ggplot(plot_dat) + 
+  # simple lines
+  geom_line(aes(x=paste0('dots ', variable), y=value, group=p_idx, color=p_diff)) +
+  scale_colour_gradient(low="#9CCC00")+
+  # box plots and jitter points, with modified x value
+  geom_violin(aes(x=variable, y=value,fill=variable), width=0.6, trim=F)+
+  scale_fill_manual(values=c("#9CCC00", "#CCFF00"))+
+  geom_boxplot(aes(x=variable, y=value),width=0.1) +
+  geom_jitter(alpha =.5, aes(x=paste0('dots ', variable), y=value, color=p_diff), position = position_jitter(width = .02)) +
+
+  # specify x value order
+ scale_x_discrete(limits=c('window1', 'dots window1', 'dots window2',
+                           'window2'))+
+scale_y_continuous(limits = c(0, 1), breaks = seq(0,1, by=.1)) +
+ theme_bw() +
+  theme(
+    plot.title   = element_text(color = "black", size = 15, face = "bold", hjust=0.5),
+    axis.title.x = element_text(color = "black", size = 13),
+    axis.title.y = element_text(color = "black", size = 13),
+    panel.border = element_blank(),
+    axis.text    = element_text(size  = 9),
+    axis.line = element_line(colour = "black", size = 1) 
+     )+
+
+          labs(
+          title = "K cell class comparison: window1 vs window2 mean power",
+          x = "Window Number",
+          y = "Power at 4Hz (normalized)" )
+dev.off()
+
+
+