analyze_MLMSS.R
#import packages
library(tidyverse)
library(plyr)
library(doParallel)
library(reticulate)
library(tidyverse)
library(reshape2)
library(ggpol)
library(rstatix)
#functions for ggplot plot layout adapted from: https://rpubs.com/Koundy/71792
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")), ...)
}
#select data sets processed with "analysis script v2.R"
datasets <- c("/media/OIC-station2/Genetics/181005 BRCA2-Halo PCNA-iRFP IR tracking/",
"/media/OIC-station2/Genetics/181010 BRCA2-Halo PCNA-iRFP IR tracking/",
"/media/OIC-station2/Genetics/181012 BRCA2-Halo PCNA-iRFP IR tracking/")
segments <- list()
for(x in datasets){ #load files
load(file.path(x,"segments_all.Rdata"))
segments_all <- ldply(segments_all,.id = NA)
segments_all$.id <- NULL
segments[[basename(x)]] <- segments_all
}
segments_all <- ldply(segments)
rm(segments)
#create unique cellID
segments_all$cellID <- paste0(segments_all$cellID,segments_all$.id)
#source functions (python)for ML-MSS analysis
source('python/ML_py.R')
#segment tracks
segments_all <- ddply(segments_all, .variables = c("condition","cellID"), function(x){
ML_segment_tracks(x)
})
save(segments_all,file = "/home/maarten/Documents/DBD paper analysis/segments_all.Rdata")
ptm <- proc.time()
#initialize cluster
nodes <- detectCores()
cl <- makeCluster(nodes-12)
registerDoParallel(cl)
#make tracklet column
segments_all <- ddply(segments_all,.variables = c("condition","cellID"),.parallel = T,function(x){
ddply(x,.variables = c("track"), .parallel = F, function(x){
segment <- 1
tracklets <- vector(length=nrow(x))
tracklets[1] <- segment
for (i in 2:nrow(x)){
if(x$state[i]==x$state[i-1]){
tracklets[i] <-segment
} else {
segment <- segment+1
tracklets[i] <-segment
}
}
x$tracklet <- paste0(x$track,".",tracklets)
return(x)
})})
stopCluster(cl)
proc.time() - ptm
#calculate MSD and MSS
numPmsd <- 4
numPmss <- 4
minLen <- 6
p <- seq(from=0.5,to=6,length.out=12)
py$pixSize <- 0.100
py$t <- 0.032
source_python('python/getMSDandMSS_R.py')
MSD_MSS <- function(x){
if(nrow(x)>10){
out <- getMSDandMSS_R(x$X,x$Y)
return(tibble("D_ML"=out[[1]],"D_Smmss"=out[[2]]))
} else {
return(tibble("D_ML"=-1.0,"D_Smmss"=-1.0))
return(NA)
}
}
segs_nest <- segments_all %>%
select(cellID,condition,tracklet,X,Y) %>%
group_by(cellID,condition,tracklet) %>%
group_modify(~MSD_MSS(.x)) %>%
inner_join(y=segments_all,by=c("cellID","condition","tracklet"))
save(segs_nest,file = "/home/maarten/Documents/DBD paper analysis/segments_tbl.Rdata")
write_csv(segs_nest,path = "/home/maarten/Documents/DBD paper analysis/segments_tbl.txt")
load(file = "/home/maarten/Documents/DBD paper analysis/segments_tbl.Rdata")
segs_nest <- segs_nest %>%
filter(condition!="WT G10 cell cycle")
segs_nest$condition <- factor(segs_nest$condition,levels=c("WT G10 -IR","WT G10 +IR","dDBD E4 -IR","dDBD E4 +IR","dCTD A2 -IR","dCTD A2 +IR","dDBDdCTD F4 -IR","dDBDdCTD F4 +IR"))
#calculate statistics
x <- segs_nest %>%
filter(D_ML>0,state==0)%>%
filter(D_ML>0,grepl("-IR",condition))%>%
dplyr::distinct(condition,cellID,tracklet,.keep_all=T)%>%
group_by(condition,cellID)%>%
dplyr::summarise(mean=mean(D_ML)) %>%
ungroup()
x$condition <- droplevels(x$condition)
rstatix::pairwise_t_test(data = x,formula = mean ~ condition,paired = F)
#plotting the diffusion rate for fast fractions
p <- segs_nest %>%
filter(D_ML>0,state==0)%>%
filter(D_ML>0,grepl("-IR",condition))%>%
dplyr::distinct(condition,cellID,tracklet,.keep_all=T)%>%
group_by(condition,cellID)%>%
dplyr::summarise(mean=mean(D_ML)) %>%
group_by(condition)%>%
ggplot(aes(x=condition,y=mean, fill=condition))+geom_boxplot()+
scale_colour_Publication()+scale_fill_Publication()+theme_Publication(base_size=18)+ theme(legend.position = "none")+
xlab("")+ylab(expression(D[app]~mu~m^{2}/s))+ylim(0,3)
p
ggsave(p,filename = "plots/diffusionrate_fast.pdf",width = 8,height = 8,units = "in" )
#analyze -IR
segs_nest %>%
filter(D_ML>0,grepl("-IR",condition))%>%
dplyr::distinct(condition,cellID,tracklet,.keep_all=T)%>%
group_by(condition,cellID,state)%>%
dplyr::summarise(mean=mean(D_ML)) %>%
group_by(condition)%>%
ggplot(aes(x=condition,y=mean,fill=condition))+geom_boxplot(notch = F)+facet_wrap(.~state,scales = "free_y")+
scale_colour_Publication()+scale_fill_Publication()+theme_Publication(base_size=12)+ theme(legend.position = "none") + ylab(expression(D[app]~mu~m^{2}/s)) + xlab("")
fractions <- segs_nest %>%
filter(D_ML>0)%>%
dplyr::distinct(condition,cellID,tracklet,.keep_all=T)%>%
group_by(condition,cellID,state)%>%
dplyr::summarise(number=n())%>%
group_by(condition,cellID) %>%
dplyr::mutate(fraction=number/sum(number))
plt <- fractions %>%
filter(state==2) %>%
ggplot(aes(y=fraction,fill=condition))+geom_boxjitter(errorbar.draw = TRUE,jitter.height = 0, jitter.width = 0.015)+ylab("immobile fraction")+
scale_colour_Publication()+scale_fill_Publication()+theme_Publication(base_size=12)
plt
ggsave(plt,filename = "plots/immobile_fractions.pdf",width = 8,height = 8,units = "in" )
#plot diffusion rate histograms
p <- segs_nest %>%
filter(D_ML>0)%>%
dplyr::distinct(condition,cellID,tracklet,.keep_all=T)%>%
mutate(state_str=as.character(state))%>%
ggplot(aes(x=D_ML,fill=state_str,y=(..count..)/tapply(..count..,..PANEL..,sum)[..PANEL..]))+geom_histogram(position="identity",alpha=0.5)+scale_x_log10(limits=c(0.0001,10))+facet_wrap(.~condition,ncol=2)+
scale_colour_Publication()+scale_fill_Publication()+theme_Publication(base_size=12)+ theme(legend.position = "none")+ylab("relative frequency")+ xlab(expression(D[app]~mu~m^{2}/s)) +
geom_text(data=subset(y,state==0 ), aes(x=1., y=0.12, label=mean_fraction), colour="#c00000", inherit.aes=FALSE, parse=FALSE)+
geom_text(data=subset(y,state==0 ), aes(x=1., y=0.10, label=paste0("+/-",sd_fraction)), colour="#c00000", inherit.aes=FALSE, parse=FALSE)+
geom_text(data=subset(y,state==1 ), aes(x=0.08, y=0.06, label=mean_fraction), colour="#fdae61", inherit.aes=FALSE, parse=FALSE)+
geom_text(data=subset(y,state==1), aes(x=.08, y=0.04, label=paste0("+/-",sd_fraction)), colour="#fdae61", inherit.aes=FALSE, parse=FALSE)+
geom_text(data=subset(y,state==2 ), aes(x=0.003, y=0.09, label=mean_fraction), colour="#1f497d", inherit.aes=FALSE, parse=FALSE)+
geom_text(data=subset(y,state==2 ), aes(x=.003, y=0.07, label=paste0("+/-",sd_fraction)), colour="#1f497d", inherit.aes=FALSE, parse=FALSE)
p
ggsave(p,filename = "plots/diffusionhistograms.pdf",width = 8,height = 8,units = "in" )
