---
title: "Data Analysis for JoVe"
output: 
  html_notebook: 
    fig_width: 15
    fig_height: 5
    fig_caption: yes
---

# Colour schemes used

```{r,include=FALSE}
library(ggplot2)
ivana_colors <- c(
  `cyan`        = "#00FFFF",
  `petrol`        = "#007e82",
  `EMBLgreen`      = "#6DA945",
  `yellow`        = "#E1E100",
  `orange`      = "#E18A33",
  `red`       = "#E10000",
  `dark red`     = "#B21217",
  `brown`     = "#B6734B",
  `light grey` = "#929297",
  `dark grey`  = "#666666",
  `black`      = "#000000")

ivana_colorsDark <- c(
  `cyan`        = "#00FFFF",
  `petrol`        = "#007E82",
  `EMBLgreen`      = "#6DA945",
  `yellow`        = "#E1E100",
  `orange`      = "#E18A33",
  `red`       = "#E10000",
  `dark red`     = "#B21217",
  `brown`     = "#B6734B",
  `light grey` = "#929297",
  `dark grey`  = "#666666",
  `white`      = "#FFFFFF")

#' Function to extract colors as hex codes
#' @scheme which colour scheme
#' @param ... Character names of 
#'
getColours <- function(scheme=ivana_colors,...) {
  cols <- c(...)

  if (is.null(cols))
    return (scheme)

  scheme[cols]
}



ivana_palettes <- list(
  `full`  = getColours(scheme=ivana_colors),

  `main`  = getColours(scheme=ivana_colors,'dark red', 'red',  'orange','EMBLgreen', 'cyan','petrol'),
  
  `main2`  = getColours(scheme=ivana_colors,'dark red', 'red',  'orange','yellow',"EMBLgreen", 'cyan','petrol'),
  
  `main2s`  = getColours(scheme=ivana_colors,'red',  'orange','yellow',"EMBLgreen", 'cyan'),
  
  `diverge`  = getColours(scheme=ivana_colors,'dark grey','yellow','dark red', 'cyan', 'red','petrol','orange','EMBLgreen','brown','dark grey','yellow','dark red', 'cyan', 'red','petrol','orange','EMBLgreen','brown'),
  
  `diverge2`  = getColours(scheme=ivana_colors,'yellow','dark red', 'cyan', 'red','petrol','orange','EMBLgreen','brown','dark gray',
                           'yellow','dark red', 'cyan', 'red','petrol','orange','EMBLgreen','brown','dark gray',
                           'yellow','dark red', 'cyan', 'red','petrol','orange','EMBLgreen','brown','dark gray',
                           'yellow','dark red', 'cyan', 'red','petrol','orange','EMBLgreen','brown','dark gray',
                           'yellow','dark red', 'cyan', 'red','petrol','orange','EMBLgreen','brown','dark gray',
                           'yellow','dark red', 'cyan', 'red','petrol','orange','EMBLgreen','brown','dark gray',
                           'yellow','dark red', 'cyan', 'red','petrol','orange','EMBLgreen','brown','dark gray',
                           'yellow','dark red', 'cyan', 'red','petrol','orange','EMBLgreen','brown','dark gray'),

  `hot`   = getColours(scheme=ivana_colors,"brown","dark red", 'red', "orange", "yellow"),

  `cold` = getColours(scheme=ivana_colors,"petrol", "EMBLgreen", "cyan", "light grey"),

  `grey`  = getColours(scheme=ivana_colors,'black', "dark grey","light grey")
)

#' Return function to interpolate a drsimonj color palette
#' @param palette e.g. list of colours
#' @param subpalette Character name of subpalette - a subset of colours from palette
#' @param reverse Boolean indicating whether the palette should be reversed
#' @param ... Additional arguments to pass to colorRampPalette()
#'
getPalettes <- function(pallete=ivana_palettes,subpalette = "main", reverse = FALSE, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("Package \"ggplot2\" needed for this function to work. Please install it.",
      call. = FALSE)}
  pal <- pallete[[subpalette]]

  if (reverse) pal <- rev(pal)
  colorRampPalette(pal, ...)

  
}

#' Color scale constructor for drsimonj colors
#' @param palette e.g. list of colours
#' @param subpalette Character name of subpalette - a subset of colours from palette
#' @param discrete Boolean indicating whether color aesthetic is discrete or not
#' @param reverse Boolean indicating whether the palette should be reversed
#' @param ... Additional arguments passed to discrete_scale() or
#'            scale_color_gradientn(), used respectively when discrete is TRUE or FALSE
#'
scale_color_KoehlerLab <- function(pallete=ivana_palettes,subpalette = "main", discrete = TRUE, reverse = FALSE, ...) {
  pal <- getPalettes(pallete=pallete,subpalette = subpalette, reverse = reverse)

  if (discrete) {
    discrete_scale("colour", paste0("KoehlerLab_", subpalette), palette = pal, ...)
  } else {
    scale_color_gradientn(colours = pal(256), ...)
  }
}

#' Fill scale constructor
#' @param palette e.g. list of colours
#' @param subpalette Character name of subpalette - a subset of colours from palette
#' @param discrete Boolean indicating whether color aesthetic is discrete or not
#' @param reverse Boolean indicating whether the palette should be reversed
#' @param ... Additional arguments passed to discrete_scale() or
#'            scale_fill_gradientn(), used respectively when discrete is TRUE or FALSE
#'
scale_fill_KoehlerLab <- function(pallete=ivana_palettes,subpalette = "main", discrete = TRUE, reverse = FALSE, ...) {
  pal <- getPalettes(pallete=pallete,subpalette = subpalette, reverse = reverse)

  if (discrete) {
    discrete_scale("fill", paste0("KoehlerLab_", subpalette), palette = pal, ...)
  } else {
    scale_fill_gradientn(colours = pal(256), ...)
  }
}




## Make my own theme

theme_white = function(base_size = 12, base_family = "",facetPanel="#6DA945") {
 
  theme_minimal(base_size = base_size, base_family = base_family) %+replace%
 
    theme(
      # Specify axis options
      axis.line = element_blank(),  
      axis.text.x = element_text(size = base_size*0.8, color = "black", lineheight = 0.9),  
      axis.text.y = element_text(size = base_size*0.8, color = "black", lineheight = 0.9),  
      axis.ticks = element_line(color = "black", size  =  0.2),  
      axis.title.x = element_text(size = base_size, color = "black", margin = margin(0, 10, 0, 0)),  
      axis.title.y = element_text(size = base_size, color = "black", angle = 90, margin = margin(0, 10, 0, 0)),  
      axis.ticks.length = unit(0.3, "lines"),   
      # Specify legend options
      legend.background = element_rect(color = NA, fill = "white"),  
      legend.key = element_rect(color = "black",  fill = "white"),  
      legend.key.size = unit(1.2, "lines"),  
      legend.key.height = NULL,  
      legend.key.width = NULL,      
      legend.text = element_text(size = base_size*0.8, color = "black"),  
      legend.title = element_text(size = base_size*0.8, face = "bold", hjust = 0, color = "black"),  
      legend.position = "right",  
      legend.text.align = NULL,  
      legend.title.align = NULL,  
      legend.direction = "vertical",  
      legend.box = NULL, 
      # Specify panel options
      panel.background = element_rect(fill = "white", color  =  NA),  
      panel.border = element_rect(fill = NA, color = "black"),  
      panel.grid.major = element_line(color = "light grey"),  
      panel.grid.minor = element_line(color = "light grey"),  
      panel.margin = unit(0.5, "lines"),   
      # Specify facetting options
      strip.background = element_rect(fill = facetPanel, color = "black"),  
      strip.text.x = element_text(size = base_size*0.8, color = "black"),  
      strip.text.y = element_text(size = base_size*0.8, color = "black",angle = -90),  
      # Specify plot options
      plot.background = element_rect(color = "white", fill = "white"),  
      plot.title = element_text(size = base_size*1.2, color = "black"),  
      plot.margin = unit(rep(1, 4), "lines")
      #Remove grid
      #panel.grid.major = element_blank(), panel.grid.minor = element_blank()
    )
  
    
}

```

# Required packages

```{r, echo=FALSE}
library(dplyr)
library(tidyr)
library("ggplot2") #
#library("grid")
library("gridExtra")#
library("scales")#
library("cowplot")#
library("minpack.lm")
library("boot")#
library("ellipse")#
library("cluster")
library("R.matlab")
library("mvoutlier")
library("data.table")
library("pls")
library(stringr)
library(mclust)

library(purrr)
```

# Graphs for JoVe 2022

## Figure 5: The Fourier ring correlation (FRC) resolution of single molecule localization microscopy images depends on the distance of the imaged z-plane to the plane of the coverslip. 

### Import and prepare data

Input data was prepared in SMAP: a modular super-resolution microscopy analysis platform for SMLM data (Ries 2020, Nat Methods 17, 870–872 (2020). https://doi.org/10.1038/s41592-020-0938-1)

```{r}
setwd("Z:/Lab_common/Manuscripts/2022_JoVE/SupplementaryFiled")
data <- read.table('220531_FRCcurves.csv',sep = ',',header = TRUE, fileEncoding = 'UTF-8-BOM') #211130_HIM-6-V5locs
# because sml files have additional chars in filename

descriptors <- read.table('SMLMDataDescriptors.csv',sep = ',',header = TRUE)

# Combine descriptors into the big FRC data table from SMAP/Matlab

df <- inner_join(data,descriptors,by=c("folder","image")) 
```


### Main


```{r}

p<-ggplot(data=df,mapping=aes(x = spatialfreq,y = curve, color =piezo.distance.from,group=piezo.distance.from))+
  stat_smooth(aes(group=piezo.distance.from),method = "loess",span=0.1,se=FALSE)+
  geom_hline(yintercept = 1/7)+facet_wrap('dye')+ylab('FRC')+xlab('Spatial frequency (nm^-1)')+scale_color_KoehlerLab(pallete = ivana_palettes, subpalette = 'main2s',discrete=F,reverse = T)+theme_bw()
p#ylim(c(0.1,0.20))+xlim(c(0.01,0.035))
```

```{r}
pdf(file=paste(Sys.Date(),"_FRC-curves-plotSmooth.pdf",sep = ""),useDingbats = TRUE,width = 12, height =6 )
p
dev.off()
```

### Inset
#### Get summary table

```{r}

summaryD <- df %>% group_by(folder,image,stage,disection,sample,piezo.distance.from,dye)%>%summarise(err=mean(errors),frc=mean(values))

summaryD
```


#### Write out summary table

```{r}
write.table(summaryD,paste(Sys.Date(),"_FRC-summary.csv",sep = ""),
              append = FALSE,sep = ',',row.names = FALSE)
```

#### Plot Figure 5 insets

```{r, fig.height=3, fig.width=6}
plotCORR<-ggplot(summaryD, mapping=aes(x=frc, y=piezo.distance.from, color=piezo.distance.from))+geom_point(size=2)+
  theme_bw()+facet_wrap(~dye)+scale_color_KoehlerLab(pallete = ivana_palettes, subpalette = 'main2s',discrete = F,reverse = T)

plotCORR
```


```{r}
pdf(file=paste(Sys.Date(),"_FRC-distance-vs-frcWithErrNEW.pdf",sep = ""),width = 6, height =3)
plotCORR
dev.off()
```



## Figure 6:  Fluorescence Lifetime Imaging (FLIM) enhanced Tau-Stimulated emission depletion microscopy can resolve two localization bands of both HTP-3 and C-terminus of SYP-5 (SYP-5::HA). 

### B

#### Import data and align the line profiles

```{r}
fit2Gaussian<-function(xA,yA,protein){
  if(protein=='HTP-3'){
    lb=c(-Inf,-Inf,-200,-Inf)
    ub=c(Inf,Inf,200,Inf)
  }else{
    lb=c(-Inf,-Inf,-80,-Inf)
    ub=c(Inf,Inf,80,Inf)
  }
  
  result<-nlsLM(y ~ (C/(w*sqrt(2*pi)) * exp(-0.5*( (x-(a-b))/w)**2) +
                       C/(w*sqrt(2*pi)) * exp(-0.5*( (x-(a+b))/w)**2) ),
                data= data.frame(x=xA,y=yA),start=list(C=50, a=0, b=20, w=50),algorithm="LM")
  out<-coef(result)
  return(out)
}


data <- read.table('FrontalViews.csv',sep = ',',header = TRUE) 
data<-data%>%ungroup()%>%group_by(image,roi)%>%mutate(ctr=(Distance-weighted.mean(Distance,Gray_Value/max(Gray_Value)))*1000,normGV=Gray_Value/max(Gray_Value))%>%group_by(image,roi,protein)%>%mutate(ctrM=ctr-fit2Gaussian(ctr,normGV,unique(protein))[2])

```

#### Plotting

```{r}
stedLP<-ggplot(data, mapping=aes(x=ctrM,y=normGV, color=protein))+
  geom_line(data=data[data$protein=='HTP-3',],aes(color=paste(protein,image,roi,sep='-')),alpha=.3,size=1)+
  geom_line(data=data[data$protein=='SYP-5::HA',],aes(color=paste(protein,image,roi,sep='-')),alpha=.3,size=1)+
  theme_bw()+theme(legend.position = 'none')+scale_color_KoehlerLab(pallete = ivana_palettes, subpalette = 'main',reverse = F)+
  stat_smooth(aes(group=protein,color=protein),size=1.5,se =F,span=.2)+scale_x_continuous(breaks = round(seq(-200, 200, by = 50),1))+xlim(c(-200,200))+
  xlab("Distance (nm)")+ylab("Normalized intensity")

stedLP
```

```{r}
pdf(file=paste(Sys.Date(),"_STED-LineProfilesUseDingbatsFin.pdf",sep = ""),width = 6, height =6,useDingbats = T)
stedLP+xlim(c(-200,200))
dev.off()
```

#### Fit distances and summarise across all line profiles


```{r}


data%>%group_by(image,roi,protein)%>%summarise(b=fit2Gaussian(ctr,normGV,unique(protein))[3],
                                               a=fit2Gaussian(ctr,normGV,unique(protein))[2],
                                               w=fit2Gaussian(ctr,normGV,unique(protein))[4],
                                               C=fit2Gaussian(ctr,normGV,unique(protein))[1])%>%ungroup()%>%group_by(protein)%>%
  summarise(n=n(),across(c('a','b','w',"C"), list(min = min, max = max,
                                            sd=sd,var=var,
                                            mean=mean,median=median)))

```



### C

#### Import data and align the line profiles

```{r}
fit2Gaussian<-function(xA,yA,protein){
  if(protein=='HTP-3'){
    lb=c(-Inf,-Inf,-200,-Inf)
    ub=c(Inf,Inf,200,Inf)
  }else{
    lb=c(-Inf,-Inf,-80,-Inf)
    ub=c(Inf,Inf,80,Inf)
  }
  
  result<-nlsLM(y ~ (C/(w*sqrt(2*pi)) * exp(-0.5*( (x-(a-b))/w)**2) +
                       C/(w*sqrt(2*pi)) * exp(-0.5*( (x-(a+b))/w)**2) ),
                data= data.frame(x=xA,y=yA),start=list(C=50, a=0, b=20, w=50),algorithm="LM")
  out<-coef(result)
  return(out)
}


data <- read.table('SideViews.csv',sep = ',',header = TRUE) 
data<-data%>%ungroup()%>%group_by(image,roi)%>%mutate(ctr=(Distance-weighted.mean(Distance,Gray_Value/max(Gray_Value)))*1000,normGV=Gray_Value/max(Gray_Value))%>%group_by(image,roi,protein)%>%mutate(ctrM=ctr-fit2Gaussian(ctr,normGV,unique(protein))[2])

```

#### Plotting

```{r}
stedLP<-ggplot(data, mapping=aes(x=ctrM,y=normGV, color=protein))+
  geom_line(data=data[data$protein=='HTP-3',],aes(color=paste(protein,image,roi,sep='-')),alpha=.3,size=1)+
  geom_line(data=data[data$protein=='SYP-5::HA',],aes(color=paste(protein,image,roi,sep='-')),alpha=.3,size=1)+
  theme_bw()+theme(legend.position = 'none')+scale_color_KoehlerLab(pallete = ivana_palettes, subpalette = 'main',reverse = F)+xlim(c(-200,200))+
  xlab("Distance (nm)")+ylab("Normalized intensity")

stedLP
```

```{r}
pdf(file=paste(Sys.Date(),"_STED-SideViews-LineProfilesUseDingbatsFin.pdf",sep = ""),width = 6, height =6,useDingbats = T)
stedLP+xlim(c(-200,200))
dev.off()
```

#### Fit distances for SYP-5::HA and HTP-3

```{r}

data%>%group_by(image,roi,protein)%>%summarise(b=fit2Gaussian(ctr,normGV,unique(protein))[3],
                                               a=fit2Gaussian(ctr,normGV,unique(protein))[2],
                                               w=fit2Gaussian(ctr,normGV,unique(protein))[4],
                                               C=fit2Gaussian(ctr,normGV,unique(protein))[1])

```


