analysis functions.R
msd_analyze_data <- function(directory,condition_list,framerate,n,fitzero,min_length,pixelsize,fitMSD,offset,max_tracks,track_file_name="tracks.simple.filtered.txt",extension=""){
segments_all <- list()
msd_fit_all <- list()
track_stats_all <- list()
for (i in 1:length(condition_list)){ #loop over the different conditions of the acquired data which is organized in sub-folders
segments <- list()
dir <- file.path(directory,condition_list[i])
filelist <- list.dirs(dir,full.names = T,recursive = F)
for(j in 1:length(filelist)){
tracks_simple <- read.csv(file.path(filelist[j],track_file_name),sep = "\t",header = F) #import tab separate file
names(tracks_simple) <- c("frame","X","Y","track","displacement","intensity","sigma","fit_error") #adjust column names
segments[[j]] <- data.frame(SEGMENT_STAT(tracks_simple),"cellID"=basename(dirname(file.path(filelist[j],track_file_name)))) #analyze the segments angle and displacements
}
track_msd <- TRACK_MSD(segments,n = n,framerate=framerate,truncate = FALSE,pxsize = pixelsize) #obtain MSD along multiple time intervals
save(track_msd,file=file.path(dir,"track_msd.Rdata"))
if(fitMSD){
tracks <- TRACK_MSD_fit(track_msd,n = n,fitzero = fitzero,framerate=framerate,offset=offset,pxsize = pixelsize) #fit linear curve to the MSD~t curve per track
save(tracks,file=file.path(dir,"msd_fit.Rdata"))
stats <- TRACK_STAT(x=segments)
save(stats,file=file.path(dir,"track_stats.Rdata"))
}
segments_all[[basename(dir)]] <- data.frame(ldply(segments),"condition"=basename(dir))
msd_fit_all[[basename(dir)]] <- data.frame(ldply(tracks),"condition"=basename(dir))
track_stats_all[[basename(dir)]] <- data.frame(ldply(stats),"condition"=basename(dir))
}
#save data to the folder
save(segments_all,file=file.path(directory,paste("segments_all_",extension,".Rdata")))
save(msd_fit_all,file=file.path(directory,paste("msd_fit_all_",extension,".Rdata")))
save(track_stats_all,file=file.path(directory,paste("track_stats_all_",extension,".Rdata")))
}
msd_histogram <- function(msd_fit_all,directory,name="",threshold=0.05,order=NULL,ymax=0.12,merge_var="cellID"){
library(plyr)
library(ggplot2)
library(reshape2)
library(ggpol)
theme_Publication <- function(base_size=14, base_family="sans") {
library(grid)
library(ggthemes)
(theme_foundation(base_size=base_size, base_family=base_family)
+ theme(plot.title = element_text(face = "bold",
size = rel(1.2), hjust = 0.5),
text = element_text(),
panel.background = element_rect(colour = NA),
plot.background = element_rect(colour = NA),
panel.border = element_rect(colour = NA),
axis.title = element_text(face = "bold",size = rel(1)),
axis.title.y = element_text(angle=90,vjust =2),
axis.title.x = element_text(vjust = -0.2),
axis.text = element_text(),
axis.line.x = element_line(colour="black"),
axis.line.y = element_line(colour="black"),
axis.ticks = element_line(),
panel.grid.major = element_line(colour="#ffffff"),
panel.grid.minor = element_blank(),
legend.key = element_rect(colour = NA),
legend.position = "bottom",
legend.direction = "horizontal",
legend.key.size= unit(0.4, "cm"),
legend.margin = unit(0.2, "cm"),
legend.title = element_text(face="italic"),
plot.margin=unit(c(10,5,5,5),"mm"),
strip.background=element_rect(colour="#ffffff",fill="#ffffff"),
strip.text = element_text(face="bold")
))
}
scale_fill_Publication <- function(...){
library(scales)
discrete_scale("fill","Publication",manual_pal(values = c("#c00000","#fdae61","#1f497d","#6599d9","#542788","#de77ae","#271d68","#6dc5aa")), ...)
}
scale_colour_Publication <- function(...){
library(scales)
discrete_scale("colour","Publication",manual_pal(values = c("#c00000","#fdae61","#1f497d","#6599d9","#542788","#de77ae","#217d68","#6dc5aa")), ...)
}
tracklength <- ldply(msd_fit_all,function(x){
nrow(x)
})
tracks <- ldply(msd_fit_all)
tracks$condition <-factor(tracks$condition, levels=names(msd_fit_all))
out <- ddply(tracks,.variables = c("condition",merge_var) ,function(x){
out <- hist(log10(x$D),breaks = seq(-7,2,0.1),plot = F)
return(out$counts/sum(out$counts))
})
n_tracks <- ddply(tracks,.variables = "condition",function(x){
nrow(x)
})
n_exp <- ddply(tracks,.variables = "condition",function(x){
length(unique(x$experiment))
})
n_cells <- ddply(tracks,.variables = "condition",function(x){
length(unique(x$cellID))
})
n_statistics <- data.frame(n_tracks$condition,"n_tracks" = n_tracks$V1,"n_cells" = n_cells$V1,"n_exp"=n_exp$V1)
write.table(n_statistics,file=file.path(directory,paste0(name,"_N_statistics.txt")),row.names = F,col.names = T)
means <- ddply(out,.variables = "condition",function(x){
colMeans(x[,-c(1,2)])
})
sds <- ddply(out,.variables = "condition",function(x){
apply(x[,-c(1,2)], 2, sd)
})
means <- melt(means)
sds <- melt(sds)
mids <- seq(-6.9,2,0.1)
mids <- 10^mids
histdata <- data.frame(".id"=means$condition,"x"=rep(mids,each=length(msd_fit_all)),'mean'=means$value,'se'=sds$value)
if(!is.null(order)){
histdata$.id <- factor(histdata$.id,levels(histdata$.id)[order])
}
#histdata$se <- 0
histdata <- na.omit(histdata)
data <- data.frame(histdata[,1:2])
names(data) <- c("x","y")
if(length(msd_fit_all)>2){
q1 <- ggplot(histdata, aes(x=x, y=mean,colour=.id,fill=.id)) +
geom_line(stat="identity")+
geom_ribbon(aes(ymin=mean-abs(se), ymax=mean+abs(se),colour=NULL), alpha=0.4)+
ylim(-0.005,ymax)+
# geom_smooth(stat="smooth",span = 0.2)+
theme(text = element_text(size=12),legend.position = "none")+
facet_wrap(~.id,nrow=2,dir="v") +scale_x_continuous(trans="log10")+
scale_colour_Publication()+scale_fill_Publication()+theme_Publication(base_size=12)+
theme(legend.position = "none")+ geom_vline(aes(xintercept=threshold),color = "red",linetype="dashed")+
xlab(expression(D[app]~mu~m^{2}/s))+ylab("")+theme(axis.line.x = element_line(color="black"),
axis.line.y = element_line(color="black"))
} else {
q1 <- ggplot(histdata, aes(x=x, y=mean,colour=.id,fill=.id)) +
geom_line(stat="identity")+
geom_ribbon(aes(ymin=mean-se, ymax=mean+se,colour=NULL), alpha=0.4)+
# geom_smooth(stat="smooth",span = 0.2)+
ylim(-0.005,ymax)+
theme(text = element_text(size=12),legend.position = "none")+
facet_wrap(~.id,nrow=1,dir="h") +scale_x_continuous(trans="log10")+
scale_colour_Publication()+scale_fill_Publication()+theme_Publication(base_size=12)+
theme(legend.position = "none")+ geom_vline(aes(xintercept=threshold),color = "red",linetype="dashed")+
xlab(expression(D[app]~mu~m^{2}/s))+ylab("")+theme(axis.line.x = element_line(color="black"),
axis.line.y = element_line(color="black"))
}
print(q1)
ggsave(filename = file.path(directory,paste0(name,"_D_histogram.pdf")),plot = q1,dpi=300,units="mm",height=100, width=150)
ggsave(filename = file.path(directory,paste0(name,"_D_histogram.png")),plot = q1+
theme(
plot.background = element_rect(fill = "transparent"),axis.line.x = element_line(colour="grey"),
axis.line.y = element_line(colour="grey"),
axis.ticks = element_line(colour="grey"),
axis.title.y = element_text(angle=90,vjust =2,colour="grey"),
axis.title.x = element_text(vjust = -0.2,colour="grey"),
axis.text = element_text(colour="grey")),bg = "transparent",dpi=400,limitsize = F,units="mm",height=100, width=150)
results <- llply(msd_fit_all,function(x) {
ddply(x,.variables = merge_var,function(x){
out <- table(x$D>threshold)/length(x$D)
# out <- c(out,x$cellID[1])
return(out)
})
}
)
results2 <- ldply(results)
results2$.id <- factor(results2$.id,names(msd_fit_all))
names(results2) <- c(".id","condition","immobile","mobile")
write.table(results2,file=file.path(directory,paste0(name,"_imm_fractions.txt")),row.names = F,col.names = T)
write.table(tracklength,file=file.path(directory,paste0(name,"_N_tracks.txt")),row.names = F,col.names = T)
if(!is.null(order)){
results2$.id <- factor(results2$.id,levels(histdata$.id)[order])
}
results2 <- na.omit(results2)
out <- ldply(results, function(x){
means <- mean(x[,2])
sems <- sd(x[,2])/sqrt(length(x[,2]))
return(data.frame("mean"=means,"se"=sems))
}
)
out$.id <- factor(out$.id,names(msd_fit_all))
#out$.id <- factor(out$.id,out$.id[c(1,2)])
#out$.id <- factor(out$.id,levels(histdata$.id)[c(1,2,3,4)])
out <- na.omit(out)
q2 <- ggplot(out, aes(x=.id, y=mean,fill=.id)) +
geom_bar(position=position_dodge(width = 1.1), stat="identity")+ geom_errorbar(aes(ymin=mean-se, ymax=mean+se),width=.2)+xlab("")+ylab("immobile fraction")+
scale_colour_Publication()+scale_fill_Publication()+theme_Publication(base_size=16)+theme(legend.position = "none")+
ylim(0,0.5)+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+ylab("immobile fraction")
print(q2)
ggsave(filename = file.path(directory,paste0(name,"_immobile_bar_graph.pdf")),plot = q2,units="mm",height=100, width=150)
ggsave(filename = file.path(directory,paste0(name,"_immobile_bar_graph.png")),plot = q2+
theme(
plot.background = element_rect(fill = "transparent"),axis.line.x = element_line(colour="grey"),
axis.line.y = element_line(colour="grey"),
axis.ticks = element_line(colour="grey"),
axis.title.y = element_text(angle=90,vjust =2,colour="grey"),
axis.title.x = element_text(vjust = -0.2,colour="grey"),
axis.text = element_text(colour="grey")),bg = "transparent",dpi=400,limitsize = F,units="mm",height=100, width=150)
if(!is.null(order)&&length(msd_fit_all)%%2==0){
datapoints <- list()
for (i in 1:length(msd_fit_all)){
if (i%%2==1){
one <- data.frame(results[[i]],'date'=ldply(strsplit(results[[i]]$cellID,split = "_"))[,1])
one_stat <- aggregate(one$'FALSE.'~one$date, data=one, FUN=function(x) c(mean=mean(x)))
two <- data.frame(results[[i+1]],'date'=ldply(strsplit(results[[i+1]]$cellID,split = "_"))[,1])
two_stat <- aggregate(two$'FALSE.'~two$date, data=two, FUN=function(x) c(mean=mean(x)))
merged <- (two_stat[,2]-one_stat[,2])/one_stat[,2]*100
datapoints[[ldply(strsplit(names(results),split = " "))[,1][i]]] <- merged
}
}
datapoints <- melt(datapoints)
datapoints <- aggregate(value~L1,data=datapoints, FUN=function(x) c(mean=mean(x)),simplify=T)
qx <- ggplot(datapoints,aes(x=L1,y=value,fill=L1))+geom_bar(stat = "identity")+
scale_colour_Publication()+scale_fill_Publication()+theme(legend.position = "none")+
xlab("")+ylab("% increase in immobilization")+ylim(c(0,100))
print(qx)
ggsave(filename = file.path(directory,paste0(name,"_immobile_increase.pdf")),plot = qx,)
}
#dotplot
q3 <- ggplot(results2,aes(x=factor(.id),y=immobile,fill=.id))+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_dotplot(binaxis = "y", stackdir = "center",dotsize = 0.7)+ xlab("")+
scale_colour_Publication()+scale_fill_Publication()+theme(legend.position = "none")+theme_Publication(base_size=12)+
ylim(0,0.5)+stat_summary(fun.data=mean_cl_normal,
geom="errorbar", color="black", width=0.3,size=1) +
stat_summary(fun.y=mean, geom="point", color="black")
print(q3)
ggsave(filename = file.path(directory,paste0(name,"_immobile_dot_plot.pdf")),plot = q3,units="mm",height=100, width=150)
ggsave(filename = file.path(directory,paste0(name,"_immobile_dot_plot.png")),plot = q3+
theme(
plot.background = element_rect(fill = "transparent"),axis.line.x = element_line(colour="grey"),
axis.line.y = element_line(colour="grey"),
axis.ticks = element_line(colour="grey"),
axis.title.y = element_text(angle=90,vjust =2,colour="grey"),
axis.title.x = element_text(vjust = -0.2,colour="grey"),
axis.text = element_text(colour="grey")),bg = "transparent",dpi=400,limitsize = F,units="mm",height=100, width=150)
#boxplot
q4 <- ggplot(results2,aes(x=factor(.id),y=immobile,fill=.id))+ theme(axis.text.x = element_text(angle = 90, hjust = 1))+
geom_boxjitter(errorbar.draw = TRUE,jitter.height = 0, jitter.width = 0.04)+ xlab("")+
scale_colour_Publication()+scale_fill_Publication()+theme(legend.position = "none")+theme_Publication(base_size=12)+
ylim(0,0.5)
print(q4)
ggsave(filename = file.path(directory,paste0(name,"_immobile_boxplot.pdf")),plot = q4,units="mm",height=100, width=150)
ggsave(filename = file.path(directory,paste0(name,"_immobile_boxplot.png")),plot = q4,bg = "transparent",dpi=400,limitsize = F,units="mm",height=100, width=75)
#qa <- grid.arrange(q1,q2,ncol=2)
}
merge_dataset <- function(dirs,conds,con_name){
if(length(dirs)!=length(conds)){
stop("dirs not equal to conds")
}
for(i in 1:length(dirs)){
load(file.path(dirs[i],"msd_fit_all.Rdata"))
date <- basename(dirs[i])
date <- strsplit(date," ")
date <- date[[1]][1]
msd_fit_all[[conds[i]]]$cellID <- paste0(date,"_",msd_fit_all[[conds[i]]]$cellID)
msd_fit_all[[conds[i]]]$experiment <- date
if (i==1) {
all <- msd_fit_all[[conds[i]]]
} else {
all <- rbind(all,msd_fit_all[[conds[i]]])
}
}
all$condition <- con_name
return(all)
}
sos_to_spoton_csv <- function(condition_list,directory,framerate,pixelsize){
condition_list <- list.dirs(directory,full.names = F,recursive = F)
#condition_list <- condition_list[c(4)]
#condition_list <- "SNAP-SiR"
segments_all <- list()
msd_fit_all <- list()
track_stats_all <- list()
#load data
for (i in 1:length(condition_list)){
segments <- list()
dir <- file.path(directory,condition_list[i])
filelist <- list.dirs(dir,full.names = T,recursive = F)
#filelist <- filelist[-grep("skip",x = filelist)]
total <- length(filelist)
# create progress bar
for(j in 1:total){
# Sys.sleep(0.1)
tracks_simple <- read.csv(file.path(filelist[j],"tracks.simple.filtered.txt"),sep = "\t",header = F)
spoton <- tracks_simple[c(1,1,4,2,3)]
spoton[,2] <- spoton[,2]/framerate/1000 #seconds
spoton[,4:5] <- spoton[,4:5]*pixelsize/1000 #um
spoton <- cbind(seq(0,nrow(spoton)-1),spoton)
names(spoton) <- c("","frame","t","trajectory","x","y")
write.csv(spoton,file = file.path(dirname(filelist[j]),paste(basename(dirname(filelist[j])),j,"tracks.spoton.txt",sep = "_")),row.names=FALSE,quote = FALSE)
}
}
}
merge_dataset_for_spoton <- function(dirs,conds,con_name,framerate,pixelsize,output){
if(length(dirs)!=length(conds)){
stop("dirs not equal to conds")
}
for(i in 1:length(dirs)){
segments <- list()
dir <- file.path(dirs[i],conds)
filelist <- list.dirs(dir,full.names = T,recursive = F)
date <- basename(dirs[i])
date <- strsplit(date," ")
date <- date[[1]][1]
#filelist <- filelist[-grep("skip",x = filelist)]
total <- length(filelist)
# create progress bar
for(j in 1:total){
# Sys.sleep(0.1)
tracks_simple <- read.csv(file.path(filelist[j],"tracks.simple.filtered.txt"),sep = "\t",header = F)
spoton <- tracks_simple[c(1,1,4,2,3)]
spoton[,2] <- spoton[,2]/framerate/1000 #seconds
spoton[,4:5] <- spoton[,4:5]*pixelsize/1000 #um
spoton <- cbind(seq(0,nrow(spoton)-1),spoton)
names(spoton) <- c("","frame","t","trajectory","x","y")
write.csv(spoton,file = file.path(output,paste(con_name,date,j,"tracks.spoton.txt",sep = "_")),row.names=FALSE,quote = FALSE)
}
}
}
