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

Adaptation of HTM_DTS_project.R

In this patch HTM_DTS_project.R has been adapted to work for platform data and torque data (already tested!).
The two R Markdown files are a fundamental part of the R file and have been updated too. 

Features still missing are:
- In project.Rmd: Two graphs per group with histograms of all periods for left turning arena from all flys and of all periods for right turning arena from all flys (only for platform data of course)

- general finetuning of the visual presentation in the output files if platform data is used
parent bbdf89f2
################## An R-script to read YAML DTS project files, visualize and statistically evaluate data ################## An R-script to read YAML DTS project files, visualize and statistically evaluate data
library(ggplot2) library(ggplot2)
library(tidyr) library(tidyr)
library(dygraphs) library(dygraphs)
library(grid) library(grid)
library(reshape2) library(reshape2)
library(dplyr) library(dplyr)
library(gridExtra) library(gridExtra)
library(yaml) library(yaml)
library(ggsignif) library(ggsignif)
library(effsize) library(effsize)
library(pwr) library(pwr)
library(BayesFactor) library(BayesFactor)
library(genefilter) library(genefilter)
library(seewave) library(seewave)
library(lubridate) library(lubridate)
library(rmarkdown) library(rmarkdown)
library(markdown) library(markdown)
library(knitr) library(knitr)
## source the script with the functions needed for analysis ## source the script with the functions needed for analysis
source("readXMLdatafile.R") source("readXMLdatafile.R")
source("DTS_plotfunctions.R") source("DTS_plotfunctions.R")
start.wd = getwd() start.wd = getwd()
## read single YAML project file ## read single YAML project file
project.file <- file.choose() project.file <- file.choose()
project.path = dirname(project.file) project.path = dirname(project.file)
project.data<-yaml.load_file(project.file) project.data<-yaml.load_file(project.file)
#make sure the evaluations are written in a subfolder of the data folder #make sure the evaluations are written in a subfolder of the data folder
evaluation.path = paste(project.path,"evaluations", sep = "/") evaluation.path = paste(project.path,"evaluations", sep = "/")
dir.create(evaluation.path, showWarnings = FALSE) dir.create(evaluation.path, showWarnings = FALSE)
setwd(evaluation.path) setwd(evaluation.path)
#how many experimental groups will need to be evaluated #how many experimental groups will need to be evaluated
NofGroups=lengths(project.data["resources"]) NofGroups=lengths(project.data["resources"])
grouped.trqhistos <- list() #Torque histograms for group in a list of length NofPeriods grouped.trqhistos <- list() #Torque histograms for group in a list of length NofPeriods
grouped.poshistos <- list() #Position histograms for group in a list of length NofPeriods grouped.poshistos <- list() #Position histograms for group in a list of length NofPeriods
grouped.PIprofiles <- list() #PIProfile data frames in a list of length NofGroups grouped.flyhistos <- list()
grouped.periods <- list() #Period designs in a list of length NofGroups grouped.flyhistosL <- list() #Platform position histograms for group in a list of length NofPeriods
grouped.spectra <- list() #Power spectra in a list of length NofGroups grouped.flyhistosR <- list() #Platform position histograms for group in a list of length NofPeriods
grouped.PIprofiles <- list() #PIProfile data frames in a list of length NofGroups
for(x in 1:NofGroups) grouped.periods <- list() #Period designs in a list of length NofGroups
{ grouped.spectra <- list() #Power spectra in a list of length NofGroups
xml_list = paste(project.path, project.data[["resources"]][[x]][["data"]], sep = "/")
for(x in 1:NofGroups)
#create/empty lists for collecting all single fly data by period {
period.data <- list() #data grouped by period xml_list = paste(project.path, project.data[["resources"]][[x]][["data"]], sep = "/")
grouped.data <- list() #total data grouped
speclist <- list() #create/empty lists for collecting all single fly data by period
period.data <- list() #data grouped by period
#start evaluating grouped.data <- list() #total data grouped
if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group have the identical experimental design speclist <- list()
{
for (l in 1:length(xml_list)) #start evaluating
{ if(MultiFlyDataVerification(xml_list)==TRUE) # make sure all flies in a group have the identical experimental design
xml_name=xml_list[l] {
for (l in 1:length(xml_list))
##### read the data with the corresponding function ####### {
singleflydata <- flyDataImport(xml_name) xml_name=xml_list[l]
##extract sequence meta-data ##### read the data with the corresponding function #######
NofPeriods = singleflydata[[5]] singleflydata <- flyDataImport(xml_name)
sequence <- singleflydata[[6]]
samplerate = as.numeric(as.character(singleflydata[[4]]$sample_rate)) ##extract sequence meta-data
NofPeriods = singleflydata[[5]]
##extract fly meta-data sequence <- singleflydata[[6]]
fly <- singleflydata[[3]] samplerate = as.numeric(as.character(singleflydata[[4]]$sample_rate))
flyname = fly$name[1]
##extract fly meta-data
##extract experiment meta-data fly <- singleflydata[[3]]
experimenter <- singleflydata[[2]] flyname = fly$name[1]
experiment <- singleflydata[[4]]
##extract experiment meta-data
##extract the rawdata experimenter <- singleflydata[[2]]
rawdata <- singleflydata[[9]] experiment <- singleflydata[[4]]
#create/empty plot lists ##extract the rawdata
poshistos <- list() rawdata <- singleflydata[[9]]
trqhistos <- list()
#create/empty plot lists
#### call RMarkdown for single fly evaluations ################################################ poshistos <- list()
rmarkdown::render('b:/GitHub/DTSevaluations/single_fly.Rmd', trqhistos <- list()
output_file = paste(flyname,"descr_anal.html", sep="_"), flyhistos <- list() #still needed?
output_dir = evaluation.path) flyhistosR <- list()
#### end RMarkdown for single fly evaluations ################################################ flyhistosL <- list()
#vectors for pooled platform periods
##move PIs to multi-experiment data.frame optomotoL_flydata <- c()
if(l>1){PIprofile <- rbind2(PIprofile, as.vector(t(sequence$lambda)))} optomotoR_flydata <- c()
##add period data to grouped data #### call RMarkdown for single fly evaluations ################################################
grouped.data[[l]] <- period.data single_fly_path <- paste(start.wd, "single_fly.Rmd", sep = "/")
rmarkdown::render(single_fly_path,
output_file = paste(flyname,"descr_anal.html", sep="_"),
} #for number of flies in xml_list output_dir = evaluation.path)
#### end RMarkdown for single fly evaluations ################################################
########### plot graphs for all experiments #####################
##move PIs to multi-experiment data.frame
##pool all data except "optomotor" by period if(l>1 & isTRUE(sequence$type[1]=="fs"||sequence$type[1]=="color")){
PIprofile <- rbind2(PIprofile, as.vector(t(sequence$lambda)))
pooled.data<-list() }
for(i in 1:NofPeriods) ##add period data to grouped data
{ grouped.data[[l]] <- period.data
period.data<-data.frame()
if(sequence$type[i]!="optomotor")
{ } #for number of flies in xml_list
for (l in 1:length(xml_list))
{ ########### plot graphs for all experiments #####################
period.data <- rbind(period.data, grouped.data[[l]][[i]])
} ##pool all data except "optomotor" by period
}
pooled.data[[i]] <- period.data pooled.data<-list()
} #for number of periods
for(i in 1:NofPeriods)
## plot pooled position and torque histograms by period ## {
period.data<-data.frame()
for(i in 1:NofPeriods) if(sequence$type[i]!="optomotor")
{ {
temp<-pooled.data[[i]] for (l in 1:length(xml_list))
{
#torque period.data <- rbind(period.data, grouped.data[[l]][[i]])
trqhistos[[i]] <- ggplot(data=temp, aes_string(temp$torque)) + }
geom_histogram(binwidth = 3, fill = sequence$histocolor[i]) + }
labs(x="torque [arb units]", y="frequency") + pooled.data[[i]] <- period.data
xlim(-600,600) + } #for number of periods
ggtitle(paste("Period", i))
## plot pooled position and torque or platform histograms by period ##
#position
if(sequence$type[i]=="fs"||sequence$type[i]=="color"||sequence$type[i]=="optomotor") for(i in 1:NofPeriods)
{ {
poshistos[[i]] <- ggplot(data=temp, aes_string(temp$a_pos)) + temp<-pooled.data[[i]]
geom_histogram(binwidth=10, fill = sequence$histocolor[i]) + if(sequence$type[i]=="fs"||sequence$type[i]=="color"||sequence$type[i]=="optomotor") {
labs(x="position [arb units]", y="frequency") + #torque
xlim(-2047,2048) + trqhistos[[i]] <- ggplot(data=temp, aes_string(temp$torque)) +
ggtitle(paste("Period", i)) geom_histogram(binwidth = 3, fill = sequence$histocolor[i]) +
} labs(x="torque [arb units]", y="frequency") +
} xlim(-600,600) +
ggtitle(paste("Period", i))
## pool all torque and position data into single data.frame
#position
all.data <- do.call(rbind, pooled.data) 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") +
## plot pooled histograms for all flies over all periods xlim(-2047,2048) +
ggtitle(paste("Period", i))
#torque } else if (sequence$type[i]=="OptomotoL"||sequence$type[i]=="OptomotoR") {
trqhistos[[NofPeriods+1]] <- ggplot(data=all.data, aes_string(all.data$torque)) + #platform
geom_histogram(binwidth=3) + flyhistos[[i]] <- ggplot(data=temp, aes_string(temp$fly)) +
labs(x="torque [arb units]", y="frequency") + geom_histogram(binwidth = 0.0175, fill = sequence$histocolor[i]) +
xlim(-600,600) + labs(x="fly [arb units]", y="frequency") +
ggtitle("Pooled Torque Histogram") xlim(-5,2) +
ggtitle(paste("Period", i))
#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") + ## pool all torque/platform and position data into single data.frame
xlim(-2047,2048) +
ggtitle("Pooled Position Histogram") all.data <- do.call(rbind, pooled.data)
} else {print("You have selected files with differing metadata")}
## plot pooled histograms for all flies over all periods
###Collect the data from each group in their respective lists and empty the variables if(sequence$type[1]=="fs"||sequence$type[1]=="color"||sequence$type[1]=="optomotor") {
#torque
#Period sequence design meta-data trqhistos[[NofPeriods+1]] <- ggplot(data=all.data, aes_string(all.data$torque)) +
grouped.periods[[x]] = periods geom_histogram(binwidth=3) +
labs(x="torque [arb units]", y="frequency") +
#Torque Histograms xlim(-600,600) +
grouped.trqhistos[[x]] = trqhistos #add torque histograms to list of grouped histograms ggtitle("Pooled Torque Histogram")
trqhistos <- list() #empty torque histograms
#position
#Position Histograms poshistos[[NofPeriods+1]] <- ggplot(data=all.data, aes_string(all.data$a_pos)) +
grouped.poshistos[[x]] = poshistos #add torque histograms to list of grouped position histograms geom_histogram(binwidth=10) +
poshistos <- list() #empty list of position histograms labs(x="position [arb units]", y="frequency") +
xlim(-2047,2048) +
#PI data ggtitle("Pooled Position Histogram")
colnames(PIprofile) <- sprintf("PI%d", 1:NofPeriods) #make colnames in PIprofile } else if (sequence$type[1]=="OptomotoL"||sequence$type[1]=="OptomotoR") {
grouped.PIprofiles[[x]] = PIprofile #add PIprofile to list of grouped PIs #platform
PIprofile <- PIprofile[0,] #empty PIprofile flyhistos[[NofPeriods+1]] <- ggplot(data=all.data, aes_string(all.data$fly)) +
geom_histogram(binwidth=0.0175) +
#Power spectra labs(x="fly [arb units]", y="frequency") +
spectemp <- do.call(cbind, speclist) #combine all power spectra xlim(-5,2) +
colnames(spectemp)[1]<-"freq" #label the first x-axis as frequency ggtitle("Pooled Platform Position Histogram")
spectemp$freq <- spectemp$freq*1000 #convert kHz to Hz }
spectemp <- spectemp[, -grep("x", colnames(spectemp))] #drop all x-axes exept the one now labelled "freq"
spectemp[length(spectemp)+1] <- rowMeans(spectemp[, grep("y", colnames(spectemp))]) #calculate the mean power spectrum in the group } else {print("You have selected files with differing metadata")}
spectemp[length(spectemp)+1] <- rowSds(spectemp[, grep("y", colnames(spectemp))])/sqrt(length(project.data[["resources"]][[x]][["data"]])) #calculate the standard deviation in the group
spectemp[, grep("y", colnames(spectemp))] <- NULL #drop all raw data for summary data ###Collect the data from each group in their respective lists and empty the variables
spectemp$group <- as.factor(rep(paste(project.data[["resources"]][[x]][["name"]], ", N=", length(project.data[["resources"]][[x]][["data"]]), sep = ""), nrow(spectemp))) #add grouping variable for later plotting
colnames(spectemp)[2] <- "mean" #Period sequence design meta-data
colnames(spectemp)[3] <- "sd" grouped.periods[[x]] = periods
grouped.spectra[[x]] = spectemp #save group mean/sd
#Torque Histograms
} #for nofGroups grouped.trqhistos[[x]] = trqhistos #add torque histograms to list of grouped histograms
trqhistos <- list() #empty torque histograms
####################################################
######## plots and statistical evaluations ######## #Platform Histograms
################################################### grouped.flyhistos[[x]] = flyhistos #add platform histograms to list of grouped p_position histograms
flyhistos <- list() #empty list of p_position histograms
###### Plots ###### grouped.flyhistosL[[x]] = flyhistosL #add platform histograms to list of grouped p_position histograms for left turning arena
flyhistosL <- list()
###generate important variables for later plotting and annotation grouped.flyhistosR[[x]] = flyhistosR #add platform histograms to list of grouped p_position histograms for right turning arena
colorrange = project.data[["statistics"]][["color-range"]] flyhistosR <- list()
boxcolors = c(colorrange[1:NofGroups])
boxes<-c(1:NofGroups) #Position Histograms
grouped.poshistos[[x]] = poshistos #add torque histograms to list of grouped position histograms
###### if there are more than two groups, try to pool the data into 'experimental' and 'control' poshistos <- list() #empty list of position histograms
if(NofGroups>2){}
#PI data
###### continue for all projects with two groups if(sequence$type[1]=="fs"||sequence$type[1]=="color"||sequence$type[1]=="optomotor") {
colnames(PIprofile) <- sprintf("PI%d", 1:NofPeriods) #make colnames in PIprofile
#### call RMarkdown for project evaluations ################################################ grouped.PIprofiles[[x]] = PIprofile #add PIprofile to list of grouped PIs
rmarkdown::render('b:/GitHub/DTSevaluations/project.Rmd', PIprofile <- PIprofile[0,] #empty PIprofile
output_file = paste(project.data$name,"html", sep = "."), }
output_dir = evaluation.path)
#### end RMarkdown for project evaluations ################################################# #Power spectra
spectemp <- do.call(cbind, speclist) #combine all power spectra
colnames(spectemp)[1]<-"freq" #label the first x-axis as frequency
setwd(start.wd) spectemp$freq <- spectemp$freq*1000 #convert kHz to Hz
\ No newline at end of file spectemp <- spectemp[-grep("x", colnames(spectemp))] #drop all x-axes exept the one now labelled "freq"
spectemp[length(spectemp)+1] <- rowMeans(spectemp[, grep("y", colnames(spectemp))]) #calculate the mean power spectrum in the group
spectemp[length(spectemp)+1] <- rowSds(spectemp[, grep("y", colnames(spectemp))])/sqrt(length(project.data[["resources"]][[x]][["data"]])) #calculate the standard deviation in the group
spectemp[, grep("y", colnames(spectemp))] <- NULL #drop all raw data for summary data
spectemp$group <- as.factor(rep(paste(project.data[["resources"]][[x]][["name"]], ", N=", length(project.data[["resources"]][[x]][["data"]]), sep = ""), nrow(spectemp))) #add grouping variable for later plotting
colnames(spectemp)[2] <- "mean"
colnames(spectemp)[3] <- "sd"
grouped.spectra[[x]] = spectemp #save group mean/sd
} #for nofGroups
####################################################
######## plots and statistical evaluations ########
###################################################
###### Plots ######
###generate important variables for later plotting and annotation
#check if statistcs are included, if not generate color range manually
if (!is.null(project.data[["statistics"]])) {
colorrange = project.data[["statistics"]][["color-range"]]
} else {colorrange = c("khaki", "olivedrab3", "cornflowerblue", "goldenrod1", "indianred1", "plum3")}
boxcolors = c(colorrange[1:NofGroups])
boxes<-c(1:NofGroups)
###### if there are more than two groups, try to pool the data into 'experimental' and 'control'
if(NofGroups>2){}
###### continue for all projects with two groups
#### call RMarkdown for project evaluations ################################################
project_path <- paste(start.wd, "project.Rmd", sep = "/")
rmarkdown::render(project_path,
output_file = paste(project.data$name,"html", sep = "."),
output_dir = evaluation.path)
#### end RMarkdown for project evaluations #################################################
setwd(start.wd)
This diff is collapsed.
<style type="text/css"> <style type="text/css">
body .main-container { body .main-container {
max-width: 1800px !important; max-width: 1800px !important;
} }
</style> </style>
```{r setup, include=FALSE} ```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(echo = TRUE)
``` ```
```{r, packages, eval=FALSE, echo=FALSE} ```{r, packages, eval=FALSE, echo=FALSE}
# load packages # load packages
library(ggplot2) library(ggplot2)
library(tidyr) library(tidyr)
library(dygraphs) library(dygraphs)
library(grid) library(grid)
library(reshape2) library(reshape2)
library(dplyr) library(dplyr)
library(gridExtra) library(gridExtra)
library(yaml) library(yaml)
library(ggsignif) library(ggsignif)
library(effsize) library(effsize)
library(pwr) library(pwr)
library(BayesFactor) library(BayesFactor)
library(genefilter) library(genefilter) # if not available for newest version: setRepositories(addURLs = c(CRANxtras = "http://www.stats.ox.ac.uk/pub/RWin"))
library(seewave) library(seewave)
## source the script with the functions needed for analysis ## source the script with the functions needed for analysis
source("readXMLdatafile.R") source("readXMLdatafile.R")
source("DTS_plotfunctions.R") source("DTS_plotfunctions.R")
``` ```
<center><h1>Single Fly Descriptive Data Evaluation Sheet</h1></center> <center><h1>Single Fly Descriptive Data Evaluation Sheet</h1></center>
<center><h2>for time series experiment `r paste(flyname)`</h2></center> <center><h2>for time series experiment `r paste(flyname)`</h2></center>
### 1. Metadata: ### 1. Metadata:
**Date and time of the experiment:** `r paste(experiment$dateTime)` **Date and time of the experiment:** `r paste(experiment$dateTime)`
**Experimenter:** `r paste(experimenter$firstname, experimenter$lastname, sep = " ")`; **ORCID:** [`r paste(experimenter$orcid)`](http://orcid.org/`r paste(experimenter$orcid)`) **Experimenter:** `r paste(experimenter$firstname, experimenter$lastname, sep = " ")`; **ORCID:** [`r paste(experimenter$orcid)`](http://orcid.org/`r paste(experimenter$orcid)`)
**Experiment description:** `r paste(experiment$description)` **Experiment description:** `r paste(experiment$description)`
**Experiment duration:** `r paste(experiment$duration)`s **Experiment duration:** `r paste(experiment$duration)`s
**Sampling rate:** `r paste(experiment$sample_rate)`Hz; **Arena type:** `r paste(experiment$arena_type)`; **Torquemeter type:** `r paste(experiment$meter_type)` **Sampling rate:** `r paste(experiment$sample_rate)`Hz; **Arena type:** `r paste(experiment$arena_type)`; **Torquemeter type:** `r paste(experiment$meter_type)`
**Fly description:** `r paste(fly$description)`; **FlyBase ID:** [`r paste(fly$flybase)`](http://flybase.org/reports/`r paste(fly$flybase)`) **Fly description:** `r paste(fly$description)`; **FlyBase ID:** [`r paste(fly$flybase)`](http://flybase.org/reports/`r paste(fly$flybase)`)
### 2. Experimental Design: ### 2. Experimental Design:
```{r expdesign, eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE, fig.align = "center", fig.height = 3, fig.width = 9} ```{r expdesign, eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE, fig.align = "center", fig.height = 3, fig.width = 9}
#plot table of period sequence #plot table of period sequence
periods=t(sequence) periods=t(sequence)
colnames(periods)<-sprintf("Period %s",1:NofPeriods)
grid.table(periods) # check if number of periods in metadata is true
``` if (ncol(periods) != NofPeriods) {
NofPeriods = ncol(periods)
### 3. Interactive Time Traces: print("sequence$type metadata faulty, function will still be working correctly")
}
```{r dytraces, eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE, fig.align = "center", fig.height = 4, fig.width = 12} colnames(periods)<-sprintf("Period %s",1:NofPeriods)
##dyplot traces grid.table(periods)
dytraces(rawdata) ```
``` ### 3. Interactive Time Traces:
### 4. Time Traces by Period: ```{r dytraces, eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE, fig.align = "center", fig.height = 4, fig.width = 12}
```{r PeriodTraces, eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE, fig.align = "center", fig.height = 6, fig.width = 10} ##dyplot traces
if (sequence$type[1] == "OptomotoR" | sequence$type[1] == "OptomotoL") {
for(i in 1:NofPeriods){ dytraces_platform(rawdata)
} else {
#save colors for later plotting dytraces(rawdata)
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 ### 4. Time Traces by Period:
temp <- rawdata[rawdata$period == i, ]
keeps = c("a_pos","torque") ```{r PeriodTraces, eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE, fig.align = "center", fig.height = 6, fig.width = 10}
period.data[[i]] <- temp[keeps] #list only position and torque data by period
for(i in 1:NofPeriods){
if(sequence$type[i]=="fs"||sequence$type[i]=="color"||sequence$type[i]=="optomotor")
{ #save colors for later plotting
## plot the torque and position time traces if(sequence$outcome[i]==0){sequence$color[i]="lightyellow"} else {sequence$color[i]="orange"}
trq_pos_traces(temp) if(sequence$outcome[i]==0){sequence$histocolor[i]="darkgreen"} else {sequence$histocolor[i]="orange"}
}
#only look at period data
#Calculate period histograms temp <- rawdata[rawdata$period == i, ]
#torque
trqhistos[[i]] <- ggplot(data=temp, aes_string(temp$torque)) + if (sequence$type[i] == "OptomotoL" || sequence$type[i] == "OptomotoR") {
geom_histogram(binwidth = 3, fill = sequence$histocolor[i]) + keeps = c("a_pos","fly")
labs(x="torque [arb units]", y="frequency") + } else {
xlim(-600,600) + keeps = c("a_pos","torque")
ggtitle(paste(flyname, "Period", i)) }
period.data[[i]] <- temp[keeps] #list only position and torque/platform data by period
#position