Unverified Commit 5db7042e authored by Maximilian von der Linde's avatar Maximilian von der Linde Committed by GitHub
Browse files

Update plotdata.R

plotdata.R has been modified so that platform data can be processed and hopefully still works for Torquemeter.
Nevertheless there are many functions which have the wrong parameter or simply have been excluded for platform data (Performance indices!!)
parent bbdf89f2
################### Plotting functions for DTS data ######################
######## Plot torque and position data for each period ########
trq_pos_traces <- function(temp)
{
##modify a_Pos to introduce gaps between -180 and 180
pos <- temp$a_pos
for(p in 2:nrow(temp))
{
if (abs(as.numeric(pos[p]))>1950 && sign(as.numeric(pos[p]))!=sign(as.numeric(pos[p-1]))) {pos[p-1]="NA"}
}
par(mar=c(5, 4, 4, 4) + 0.1)
traces <- plot(x = temp$time, y = temp$a_pos, type = "n", axes=FALSE, xaxs = "i", yaxs = "i", ylim=c(-2047,2048), ylab = "",xlab="")
#shade quadrants
rect(temp$time[1],-1538,temp$time[nrow(temp)],-512, col = "grey95")
rect(temp$time[1],512,temp$time[nrow(temp)],1538, col = "grey95")
par(new=TRUE)
plot(x = temp$time, y = pos, type = "l", col="red3", xaxs = "i", yaxs = "i", ylim=c(-2047,2048), ylab = "position[arb.units]",xlab="time [ms]")
par(new=TRUE)
plot(x = temp$time, temp$torque, type = "l", col="blue", xaxs = "i", ylim=c(-500,500), main = paste("Fly Traces", flyname, "Period", i), axes=F, xlab=NA, ylab="")
lines(c(temp$time[1],temp$time[nrow(temp)]),c(0,0),type="l",lty=1,lwd=1, col="black")
axis(4)
mtext("torque [arb_units]", side = 4, line = 3)
return(traces)
}
dytraces <- function(rawdata)
{
traces <- dygraph(rawdata, main = paste("Time Traces", flyname)) %>%
dySeries("a_pos", label = "position", color = "darkred") %>%
dySeries("torque", axis = 'y2', color = "blue") %>%
dyAxis("x", drawGrid = FALSE) %>%
dyAxis("y", label = "Position [arb_units]", valueRange = c(-2048,2048)) %>%
dyOptions(gridLineColor = "lightgrey") %>%
dyAxis("y2", label = "Torque [arb_units]", independentTicks = TRUE, valueRange = c(-700,700), drawGrid = FALSE) %>%
dyOptions(includeZero = TRUE) %>%
dyRangeSelector()
return(traces)
}
############ Function to plot several plots into a single plot #########
multiplot <- function(..., plotlist = NULL, file, cols = 1, layout = NULL)
{
require(grid)
plots <- c(list(...), plotlist)
numPlots = length(plots)
if (is.null(layout)) {
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots == 1) {
print(plots[[1]])
} else {
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
for (i in 1:numPlots) {
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
############ Function to annotate plots after stats #########
#Wilcoxon annotations
wilcox.annotate <- function(boxes, wilcoxon)
{
annotate("text",
x=boxes,
y=1.1,
label=paste("p=",wilcoxon[boxes]))
}
#samplesizes annotations
samplesizes.annotate <- function(boxes, samplesizes)
{
annotate("text",
x=boxes,
y=-1.1,
label=paste("N=", samplesizes[boxes]))
}
########### Function to calculate False Positive Risk
#======================================================
calc.FPR = function(samplesizes,pval,prior,delta1){
sdiff=sqrt(1/samplesizes[1] + 1/samplesizes[2])
df=samplesizes[1]+samplesizes[2]-2
# Note FPR doesn't need calculation of power for p-equals case
#
#under H0, use central t distribution
tcrit=qt((1-pval/2),df,ncp=0)
x0=tcrit
y0=dt(x0,df,0)
#
# under H1 use non-central t distribution
ncp1=delta1/sdiff #non-centrality paramater
x1=x0 #tcrit
y1=dt(x1,df,ncp=ncp1)
# Calc false positive risk
p0=2*y0
p1=y1
FPR=((1-prior)*p0)/(((1-prior)*p0) + prior*p1)
FPR
output=c(FPR,x0,y0,x1,y1)
return(output)
################### Plotting functions for DTS data ######################
######## Plot torque and position data for each period ########
trq_pos_traces <- function(temp)
{
##modify a_Pos to introduce gaps between -180 and 180
pos <- temp$a_pos
for(p in 2:nrow(temp))
{
if (abs(as.numeric(pos[p]))>1950 && sign(as.numeric(pos[p]))!=sign(as.numeric(pos[p-1]))) {pos[p-1]="NA"}
}
par(mar=c(5, 4, 4, 4) + 0.1)
traces <- plot(x = temp$time, y = temp$a_pos, type = "n", axes=FALSE, xaxs = "i", yaxs = "i", ylim=c(-2047,2048), ylab = "",xlab="")
#shade quadrants
rect(temp$time[1],-1538,temp$time[nrow(temp)],-512, col = "grey95")
rect(temp$time[1],512,temp$time[nrow(temp)],1538, col = "grey95")
par(new=TRUE)
plot(x = temp$time, y = pos, type = "l", col="red3", xaxs = "i", yaxs = "i", ylim=c(-2047,2048), ylab = "position[arb.units]",xlab="time [ms]")
par(new=TRUE)
plot(x = temp$time, temp$torque, type = "l", col="blue", xaxs = "i", ylim=c(-500,500), main = paste("Fly Traces", flyname, "Period", i), axes=F, xlab=NA, ylab="")
lines(c(temp$time[1],temp$time[nrow(temp)]),c(0,0),type="l",lty=1,lwd=1, col="black")
axis(4)
mtext("torque [arb_units]", side = 4, line = 3)
return(traces)
}
dytraces <- function(rawdata)
{
traces <- dygraph(rawdata, main = paste("Time Traces", flyname)) %>%
dySeries("a_pos", label = "position", color = "darkred") %>%
dySeries("torque", axis = 'y2', color = "blue") %>%
dyAxis("x", drawGrid = FALSE) %>%
dyAxis("y", label = "Position [arb_units]", valueRange = c(-2048,2048)) %>%
dyOptions(gridLineColor = "lightgrey") %>%
dyAxis("y2", label = "Torque [arb_units]", independentTicks = TRUE, valueRange = c(-700,700), drawGrid = FALSE) %>%
dyOptions(includeZero = TRUE) %>%
dyRangeSelector()
return(traces)
}
dytraces_platform <- function(rawdata)
{
traces <- dygraph(rawdata, main = paste("Time Traces", flyname)) %>%
dySeries("a_pos", label = "position", color = "darkred") %>%
dySeries("fly", axis = 'y2', color = "blue") %>%
dyAxis("x", drawGrid = FALSE) %>%
dyAxis("y", label = "Position [arb_units]") %>%
dyOptions(gridLineColor = "lightgrey") %>%
dyAxis("y2", label = "Torque [arb_units]", independentTicks = TRUE, drawGrid = FALSE) %>%
dyOptions(includeZero = TRUE) %>%
dyRangeSelector()
return(traces)
}
############ Function to plot several plots into a single plot #########
multiplot <- function(..., plotlist = NULL, file, cols = 1, layout = NULL)
{
require(grid)
plots <- c(list(...), plotlist)
numPlots = length(plots)
if (is.null(layout)) {
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots == 1) {
print(plots[[1]])
} else {
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
for (i in 1:numPlots) {
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
############ Function to annotate plots after stats #########
#Wilcoxon annotations
wilcox.annotate <- function(boxes, wilcoxon)
{
annotate("text",
x=boxes,
y=1.1,
label=paste("p=",wilcoxon[boxes]))
}
#samplesizes annotations
samplesizes.annotate <- function(boxes, samplesizes)
{
annotate("text",
x=boxes,
y=-1.1,
label=paste("N=", samplesizes[boxes]))
}
########### Function to calculate False Positive Risk
#======================================================
calc.FPR = function(samplesizes,pval,prior,delta1){
sdiff=sqrt(1/samplesizes[1] + 1/samplesizes[2])
df=samplesizes[1]+samplesizes[2]-2
# Note FPR doesn't need calculation of power for p-equals case
#
#under H0, use central t distribution
tcrit=qt((1-pval/2),df,ncp=0)
x0=tcrit
y0=dt(x0,df,0)
#
# under H1 use non-central t distribution
ncp1=delta1/sdiff #non-centrality paramater
x1=x0 #tcrit
y1=dt(x1,df,ncp=ncp1)
# Calc false positive risk
p0=2*y0
p1=y1
FPR=((1-prior)*p0)/(((1-prior)*p0) + prior*p1)
FPR
output=c(FPR,x0,y0,x1,y1)
return(output)
}
\ No newline at end of file
################## An R-script to import DTS data and to plot it in several ways
library(ggplot2)
library(tidyr)
library(dygraphs)
library(grid)
library(reshape2)
library(dplyr)
## source the script with the functions needed for analysis
source("readXMLdatafile.R")
source("DTS_plotfunctions.R")
#create lists for collecting all single fly data by period
period.data <- list()
grouped.data <- list()
############# read file list and plot graphs for each file #############################
xml_list <- choose.files()
if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group have the identical experimental design
{
for (l in 1:length(xml_list))
{
xml_name=xml_list[l]
##### read the data with the corresponding function #######
singleflydata <- flyDataImport(xml_name)
##extract sequence meta-data
NofPeriods = singleflydata[[5]]
sequence <- singleflydata[[6]]
##extract fly meta-data
fly <- singleflydata[[3]]
flyname = fly$name[1]
##extract the rawdata
rawdata <- singleflydata[[9]]
#create plot lists
poshistos <- list()
trqhistos <- list()
##### analyze data from each period separately #######
for(i in 1:NofPeriods){
#save colors for later plotting
if(sequence$outcome[i]==0){sequence$color[i]="lightyellow"} else {sequence$color[i]="orange"}
if(sequence$outcome[i]==0){sequence$histocolor[i]="darkgreen"} else {sequence$histocolor[i]="orange"}
#only look at period data
temp <- rawdata[rawdata$period == i, ]
keeps = c("a_pos","torque")
period.data[[i]] <- temp[keeps] #list only position and torque data by period
if(sequence$type[i]!="optomotor")
{
## plot the torque and position time traces
filename = paste(flyname,"timetraces",i,".png", sep="_")
png(file = filename, width = 1920) # direct the following output to the image
trq_pos_traces(temp)
graphics.off()
}
#plot period histograms
#torque
trqhistos[[i]] <- ggplot(data=temp, aes_string(temp$torque)) +
geom_histogram(binwidth = 3, fill = sequence$histocolor[i]) +
labs(x="torque [arb units]", y="frequency") +
xlim(-600,600) +
ggtitle(paste(flyname, "Period", i))
#position
if(sequence$type[i]!="optomotor")
{
poshistos[[i]] <- ggplot(data=temp, aes_string(temp$a_pos)) +
geom_histogram(binwidth=10, fill = sequence$histocolor[i]) +
labs(x="position [arb units]", y="frequency") +
xlim(-2047,2048) +
ggtitle(paste(flyname, "Period", i))
}
##calculate PIs
##create data.frame for adding PIs if it doesn't eist, yet
if(!exists("PIprofile")){PIprofile <- data.frame(matrix(ncol = NofPeriods))}
#for position
if(sequence$type[i]!="optomotor")
{
t1 = sum(abs(temp$a_pos) >= 512 & abs(temp$a_pos) <= 1538)
t2 = nrow(temp)-t1
sequence$lambda[i] = (t1-t2)/(t1+t2)
if (sequence$contingency[i] == '1_3_Q'){sequence$lambda[i]=-sequence$lambda[i]}
if(l==1){PIprofile[1,i]=sequence$lambda[i]}
}
#for torque
#if(sequence$type[i]=="yt"||"sw_blue"||"sw_green")
#{}
} #for Number of Periods
########## analyze data for entire experiment ##############
## plot position and torque histograms ##
#torque
trqhistos[[NofPeriods+1]] <- ggplot(data=rawdata, aes_string(rawdata$torque)) +
geom_histogram(binwidth=3) +
labs(x="torque [arb units]", y="frequency") +
xlim(-600,600) +
ggtitle(paste(flyname, "total"))
#position
poshistos[[NofPeriods+1]] <- ggplot(data=rawdata, aes_string(rawdata$a_pos)) +
geom_histogram(binwidth=10) +
labs(x="position [arb units]", y="frequency") +
xlim(-2047,2048) +
ggtitle(paste(flyname, "total"))
##write histograms to file
#torque
filename = paste(flyname,"torque_histos.png", sep="_")
png(file = filename, width = 1000) # direct the following output to the image
multiplot(plotlist = trqhistos, cols=2)
dev.off()
#position
filename = paste(flyname,"position_histos.png", sep="_")
png(file = filename, width = 1000) # direct the following output to the image
multiplot(plotlist = poshistos, cols=2)
dev.off()
##plot PI bargraph
filename = paste(flyname,"performance.png", sep="_")
png(file = filename, width = 1000) # direct the following output to the image
barplot(sequence$lambda, main = paste("Performance Indices", flyname),
xlab="Time [periods]",
ylab="PI [rel. units]",
space=0,
col = sequence$color)
graphics.off()
##dyplot traces
traces <- dytraces(rawdata)
print(traces)
##move PIs to multi-experiment data.frame
if(l>1){PIprofile <- rbind2(PIprofile, as.vector(t(sequence$lambda)))}
##add period data to grouped data
grouped.data[[l]] <- period.data
} #for number of flies
########### plot graphs for all experiments #####################
##pool all data except "optomotor" by period
pooled.data<-list()
for(i in 1:NofPeriods)
{
period.data<-data.frame()
if(sequence$type[i]!="optomotor")
{
for (l in 1:length(xml_list))
{
period.data <- rbind(period.data, grouped.data[[l]][[i]])
}
}
pooled.data[[i]] <- period.data
}
## plot pooled position and torque histograms by period ##
for(i in 1:NofPeriods)
{
temp<-pooled.data[[i]]
#torque
trqhistos[[i]] <- ggplot(data=temp, aes_string(temp$torque)) +
geom_histogram(binwidth = 3, fill = sequence$histocolor[i]) +
labs(x="torque [arb units]", y="frequency") +
xlim(-600,600) +
ggtitle(paste("Period", i))
#position
if(sequence$type[i]!="optomotor")
{
poshistos[[i]] <- ggplot(data=temp, aes_string(temp$a_pos)) +
geom_histogram(binwidth=10, fill = sequence$histocolor[i]) +
labs(x="position [arb units]", y="frequency") +
xlim(-2047,2048) +
ggtitle(paste("Period", i))
}
}
## pool all torque and position data into single data.frame
all.data <- do.call(rbind, pooled.data)
## plot pooled histograms for all flies over all periods
#torque
trqhistos[[NofPeriods+1]] <- ggplot(data=all.data, aes_string(all.data$torque)) +
geom_histogram(binwidth=3) +
labs(x="torque [arb units]", y="frequency") +
xlim(-600,600) +
ggtitle("Pooled Torque Histogram")
#position
poshistos[[NofPeriods+1]] <- ggplot(data=all.data, aes_string(all.data$a_pos)) +
geom_histogram(binwidth=10) +
labs(x="position [arb units]", y="frequency") +
xlim(-2047,2048) +
ggtitle("Pooled Position Histogram")
###write pooled histograms to file###
#torque
filename = paste(flyname,"torque_histos.png", sep="_")
png(file = "pooled_torque_histos.png", width = 1000) # direct the following output to the image
multiplot(plotlist = trqhistos, cols=2)
dev.off()
#position
png(file = "pooled_position_histos.png", width = 1000) # direct the following output to the image
multiplot(plotlist = poshistos, cols=2)
dev.off()
###make colnames in PIprofile for plotting###
colnames(PIprofile) <- sprintf("PI%d", 1:NofPeriods)
# Plot box&dotplot with notches
ggplot(melt(PIprofile), aes(variable, value)) +
geom_boxplot(fill = unique(sequence$color), notch = TRUE, outlier.color=NA, width=0.8, size=0.6) +
geom_jitter(data = melt(PIprofile), aes(variable, value), position=position_jitter(0.3), cex=2, color="grey80") +
ggtitle(paste("PI Profile, N=",length(xml_list))) +
ylim(-1,1)+theme_light(base_size = 18) + theme(panel.grid.major.x = element_blank() ,panel.border = element_rect(size = 0.5, linetype = "solid", colour = "black", fill=NA)) +
theme(axis.text.y = element_text(size=18))+ ylab("PI [rel. units]") + theme(aspect.ratio=4/NofPeriods)
# plot violin plot
ggplot(melt(PIprofile), aes(variable, value)) +
geom_violin(fill = unique(sequence$color), adjust = .6) +
stat_summary(fun.y = mean, geom = "point", position = position_dodge(width = .9),
size = 6, shape = 4) +
ggtitle(paste("PI Profile, N=",length(xml_list))) +
ylim(-1,1)+theme_light(base_size = 18) + theme(panel.grid.major.x = element_blank() ,panel.border = element_rect(size = 0.5, linetype = "solid", colour = "black", fill=NA)) +
theme(axis.text.y = element_text(size=18))+ ylab("PI [rel. units]") + theme(aspect.ratio=4/NofPeriods)
# plot bar plot with sem
#compute summary statistics
error <- data.frame(period=integer(ncol(PIprofile)), mean=numeric(ncol(PIprofile)), sem=numeric(ncol(PIprofile)))
for(e in 1:ncol(PIprofile))
{
error$period[e]=e
error$mean[e]=mean(PIprofile[,e])
error$sem[e]=sd(PIprofile[,e]/sqrt(ncol(PIprofile)))
}
# plot graph
ggplot(error, aes(x=period, y=mean)) +
geom_bar(position=position_dodge(), stat="identity", fill = unique(sequence$color), colour="black") +
geom_errorbar(aes(ymin=mean-sem, ymax=mean+sem),
width=0,
size=1.5,
position=position_dodge(.9)) +
ggtitle(paste("PI Profile, N=",length(xml_list))) +
ylim(-1,1)+theme_light(base_size = 18) + theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank(), panel.border = element_rect(size = 0.5, linetype = "solid", colour = "black", fill=NA)) +
theme(axis.text.y = element_text(size=18))+ ylab("PI [rel. units]") + theme(aspect.ratio=4/NofPeriods)
} else {print("You have selected files with differing metadata")}
\ No newline at end of file
################## An R-script to import DTS data and to plot it in several ways
library(ggplot2)
library(tidyr)
library(dygraphs)
library(grid)
library(reshape2)
library(dplyr)
## source the script with the functions needed for analysis
source("readXMLdatafile.R")
source("DTS_plotfunctions.R")
#create lists for collecting all single fly data by period
period.data <- list()
grouped.data <- list()
############# read file list and plot graphs for each file #############################
xml_list <- choose.files()
if(MultiFlyDataVerification(xml_list)==TRUE) { # make sure all flies in a group have the identical experimental design
print("identical")
for (l in 1:length(xml_list))
{
xml_name=xml_list[l]
##### read the data with the corresponding function #######
singleflydata <- flyDataImport(xml_name)
##extract sequence meta-data
NofPeriods = singleflydata[[5]]
sequence <- singleflydata[[6]]
##extract fly meta-data
fly <- singleflydata[[3]]
flyname = fly$name[1]
##extract the rawdata
rawdata <- singleflydata[[9]]
#create plot lists
poshistos <- list()
trqhistos <- list()
flyhistos <- list()
##### analyze data from each period separately #######
for(i in 1:NofPeriods){
#save colors for later plotting
if(sequence$outcome[i]==0){sequence$color[i]="lightyellow"} else {sequence$color[i]="orange"}
if(sequence$outcome[i]==0){sequence$histocolor[i]="darkgreen"} else {sequence$histocolor[i]="orange"}
#only look at period data
temp <- rawdata[rawdata$period == i, ]