Unverified Commit 1c933ec1 authored by Maximilian von der Linde's avatar Maximilian von der Linde Committed by GitHub
Browse files

Add files via upload

parent 74fda16d
################## 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()
##### 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, ]
if (sequence$type[i] == "OptomotoL") or (sequence$type[i] == "OptomotoR") {
keeps = c("a_pos","fly")
period.data[[i]] <- temp[keeps] #list only position and torque data by period
}
else{
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")}
################## Some functions to import XML DTS data
require("XML")
####################### Functions for importing data ##########################
## Importing DTS data from an XML file
flyDataImport <- function(xml_name) {
## Import the data from the .xml file.
flyData <- xmlParse(xml_name)
flyDataXMLtop = xmlRoot(flyData)
#parse the metadata
URIs <- xmlToDataFrame(nodes=getNodeSet(flyData,"//metadata/URIs"))
experimenter <- xmlToDataFrame(nodes=getNodeSet(flyData,"//metadata/experimenter"))
fly <- xmlToDataFrame(nodes=getNodeSet(flyData,"//metadata/fly"))
experiment <- xmlToDataFrame(nodes=getNodeSet(flyData,"//metadata/experiment"))
#parse sequence data
NofPeriods = as.integer(xmlGetAttr(flyDataXMLtop[['sequence']], "periods"))
sequence <- xmlToDataFrame(nodes=getNodeSet(flyData,"//sequence/period"))
#parse time series meta-data
CSV_descriptor <- xmlToDataFrame(nodes=getNodeSet(flyData,"//timeseries/CSV_descriptor"))
variables <- xmlToDataFrame(nodes=getNodeSet(flyData,"//timeseries/variables/variable"))
#parse the time series raw data
rawdata <- read.table(text=xmlSApply(flyDataXMLtop[['timeseries']][['csv_data']], xmlValue), col.names=variables$type)
options(digits.secs = 3)
rawdata$date<-as.POSIXct(rawdata$time/1000, origin="1970-01-01", tz="UTC")
##list all data