Add files via upload

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()
+
+
+

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)
+

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)
+