---
title: "Tracking Analysis for Kathi"
author: "Julian Stopp"
date: '2022-02-04'
output:
  pdf_document: default
  html_document: default
updated: '2022-02-04'
---

```{r setup, include=FALSE}
# Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.

# Setup knitr with standard formattingg changes
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE, tidy = TRUE)

## Package Loading ## 
Packages = c("data.table", "tidyr", "plotly", "dplyr", "ggplot2");
invisible(lapply(Packages, library, character.only = TRUE))

## Empty workspace ##
rm(list=ls())

## Set working directory
knitr::opts_knit$set(root.dir = r"(.\Files)")
```

```{r Functions, include=FALSE}
# Function to randomly select specific groups and ALL their rows
sample_n_groups = function(grouped_df, size, replace = FALSE, weight=NULL) {
  grp_var <- grouped_df %>% 
    groups %>%
    unlist %>% 
    as.character
  random_grp <- grouped_df %>% 
    dplyr::summarise() %>% 
    dplyr::sample_n(size, replace, weight) %>% 
    dplyr::mutate(unique_id = 1:NROW(.))
  grouped_df %>% 
    dplyr::right_join(random_grp, by=grp_var) %>% 
    dplyr::group_by_(grp_var)
}

# Outlier Removal
Outlier_rem <- function(.dat, Col, probMin=.25, probMax=.75) {
  Q <- quantile(.dat[[Col]], probs=c(probMin, probMax), na.rm = FALSE)
  iqr <- IQR(.dat[[Col]])
  up <-  Q[2]+1.5*iqr  # Upper Range  
  low <- Q[1]-1.5*iqr  # Lower Range

return(subset(.dat, .dat[[Col]] > (Q[1] - 1.5*iqr) & .dat[[Col]] < (Q[2]+1.5*iqr)))
}

# Angle calculation taken from: https://www.r-bloggers.com/2007/05/calculate-turning-angles-and-step-lengths-from-location-data/
anglefun <- function(xx,yy,bearing=TRUE,as.deg=FALSE){
  ## calculates the compass bearing of the line between two points
  ## xx and yy are the differences in x and y coordinates between two points
  ## Options:
  ## bearing = FALSE returns +/- pi instead of 0:2*pi
  ## as.deg = TRUE returns degrees instead of radians
  c = 1
  if (as.deg){
    c = 180/pi
  }
  
  b<-sign(xx)
  b[b==0]<-1  #corrects for the fact that sign(0) == 0
  tempangle = b*(yy<0)*pi+atan(xx/yy)
  if(bearing){
    #return a compass bearing 0 to 2pi
    #if bearing==FALSE then a heading (+/- pi) is returned
    tempangle[tempangle<0]<-tempangle[tempangle<0]+2*pi
  }
  return(tempangle*c)
}
bearing.ta <- function(loc1,loc2,loc3,as.deg=FALSE){
  #browser()
  ## calculates the bearing and length of the two lines
  ##    formed by three points
  ## the turning angle from the first bearing to the
  ##    second bearing is also calculated
  ## locations are assumed to be in (X,Y) format.
  ## Options:
  ## as.deg = TRUE returns degrees instead of radians
  if (length(loc1) != 2 | length(loc2) != 2 | length(loc3) !=2){
    print("Locations must consist of either three vectors, length == 2,
or three two-column dataframes")
    return(NaN)
  }
  c = 1
  if (as.deg){
    c = 180/pi
  }
  
  locdiff1<-loc2-loc1
  locdiff2<-loc3-loc2
  bearing1<-anglefun(locdiff1[1],locdiff1[2],bearing=F)
  bearing2<-anglefun(locdiff2[1],locdiff2[2],bearing=F)
  if(is.data.frame(locdiff1)){
    dist1<-sqrt(rowSums(locdiff1^2))
    dist2<-sqrt(rowSums(locdiff2^2))
  }else{
    dist1<-sqrt(sum(locdiff1^2))
    dist2<-sqrt(sum(locdiff2^2))
  }
  
  ta=(bearing2-bearing1)
  
  ta[ta < -pi] = ta[ta < -pi] + 2*pi
  ta[ta > pi] = ta[ta > pi] - 2*pi
  #return(list(bearing1=unlist(bearing1*c),bearing2=unlist(bearing2*c), ta=unlist(ta*c),dist1=unlist(dist1),dist2=unlist(dist2)))
  return(ta=unlist(ta*c))
}
```

```{r Data Import - Positions & Area from imaris}
# Reads names of all position csv files in folder
data_to_read = list.files(pattern = "\\.csv$",
                          recursive = F, # reads subfolders
                          full.names = F)

# Reads all .tab files and creates a list
mydata <- data.table::rbindlist(sapply(data_to_read, data.table::fread, simplify = FALSE),
                    use.names = TRUE, idcol = "FileName")

# Convert to dataframe & reorganize & rename
cols = c("Global_ID", "TrackID", "Replicate") 
mydata <- mydata %>%
  dplyr::rename(
    TrackID = 'Track n°',
    FRAME = 'Slice n°',
    POSITION_X = X,
    POSITION_Y = Y
  ) %>%
  tidyr::separate(FileName, c("Genotype", "Replicate", "Rest"), sep = "_", extra ="merge") %>%
  tidyr::unite("Global_ID", TrackID, Genotype, Replicate, remove = F, sep = "_") %>%
  dplyr::select(-c("Pixel Value")) %>%
  dplyr::mutate_at(cols, factor) %>%
  dplyr::arrange(Global_ID)
```

```{r Data Manipulation}
# Calculate Speed & Gradient
## Correction term for time: we have 60 sec intervals, so we need to divide by 60 to get µm/sec
timeCor <- 3*60

pixelsize = 0.473

mydata <- mydata %>%
  dplyr::group_by(Global_ID) %>%
  dplyr::mutate(
    Time = FRAME * 3,
    
    FRAME_dif = FRAME-dplyr::lag(FRAME),
    
    xDis = POSITION_X - dplyr::lag(POSITION_X),
    yDis = POSITION_Y - dplyr::lag(POSITION_Y),
    Displacement = (sqrt((xDis^2)+(yDis^2))) * pixelsize,
    
    myNSpeed = Displacement/FRAME_dif/timeCor,
    IMSpeed = 
      if_else(row_number() == 1, dplyr::lead(Displacement) / ((dplyr::lead(FRAME) - FRAME) * timeCor),
              if_else(row_number() == n(), Displacement / ((FRAME - dplyr::lag(FRAME)) * timeCor), 
                      (dplyr::lead(Displacement) + Displacement) / ((dplyr::lead(FRAME) - dplyr::lag(FRAME)) * timeCor))),

# Correct Positions for Spider Plot
    POSITION_X_Norm = POSITION_X-dplyr::first(POSITION_X),
    POSITION_Y_Norm = POSITION_Y-dplyr::first(POSITION_Y),
     
# Calculate angles    
    Time_dif = Time-dplyr::lag(Time),
    
    POSITION_X_lag = lag(POSITION_X),
    POSITION_X_lead = lead(POSITION_X),
    POSITION_Y_lag = lag(POSITION_Y),
    POSITION_Y_lead = lead(POSITION_Y)
    ) %>% 
  dplyr::rowwise() %>%
  dplyr::mutate(
    directional_change = bearing.ta(c(POSITION_X_lag, POSITION_Y_lag), c(POSITION_X, POSITION_Y), c(POSITION_X_lead, POSITION_Y_lead), as.deg=T),
    directional_change_abs = abs(directional_change)
  ) %>%
  dplyr::ungroup()
```

```{r Summary}
mysum <- mydata %>%
  dplyr::group_by(Global_ID, Genotype) %>%
  dplyr::summarise(
    xTotDis = dplyr::last(POSITION_X) - dplyr::first(POSITION_X),
    yTotDis = dplyr::last(POSITION_Y) - dplyr::first(POSITION_Y),
    EucDis = sqrt((xTotDis^2)+(yTotDis^2)) * pixelsize,
    TotDis = sum(Displacement, na.rm = T),
    myNSpeed = median(myNSpeed, na.rm=T),
    myNSpeed_max = max(myNSpeed, na.rm=T),
    IMSpeed = median(IMSpeed, na.rm=T),
    IMSpeed_max = max(IMSpeed, na.rm=T),
    Velocity = median(Velocity, na.rm=T),
    Timepoints = n(),
    directional_change = median(directional_change, na.rm=T),
    directional_change_abs = median(directional_change_abs, na.rm=T) 
  ) %>%
  dplyr::mutate(
    Straightness = EucDis/TotDis
  ) %>%
  dplyr::ungroup()
```

```{r Plotting}
mysum %>%
  Outlier_rem("IMSpeed", 0.05, 0.95) %>%
  plot_ly(x=~Genotype, y=~IMSpeed, color = ~Genotype,
          type="box")
mysum %>%
  Outlier_rem("myNSpeed", 0.05, 0.95) %>%
  plot_ly(x=~Genotype, y=~myNSpeed, color = ~Genotype,
          type="box")

mysum %>%
  Outlier_rem("Straightness", 0.05, 0.95) %>%
  plot_ly(x=~Genotype, y=~Straightness, color = ~Genotype,
          type="box")

mydata %>%
  tidyr::drop_na() %>%
  Outlier_rem("directional_change_abs", 0.05, 0.95) %>%
  plot_ly(x=~Genotype, y=~directional_change_abs, color = ~Genotype,
          type="box")
          
```{r Trajectory}
ggplotly(
  ggplot(data = mydata, aes(x = POSITION_X, y = POSITION_Y, color = Global_ID, group = 1,
                                text = paste(Global_ID))) +
    geom_path() +
    geom_point(data =  mydata %>%
                 dplyr::group_by(Global_ID) %>%
                 dplyr::slice(which.max(Time)) %>%
                 droplevels(),
               shape = 1, stroke = 0.2, size = 1.6, alpha = 1, colour = "black") +
    geom_point(data =  mydata %>%
                 dplyr::group_by(Global_ID) %>%
                 dplyr::slice(which.min(Time)) %>%
                 droplevels(),
               shape = 1, stroke = 0.2, size = 1, alpha = 1, colour = "black") +
      scale_y_reverse() +
    facet_wrap(~Genotype),
  tooltip = c("text")
) %>%
  layout(
    title = "<b>Cell Trajectories</b>",
    showlegend = F)
```

```{r Vector}
myVdata <- mydata %>%
  dplyr::select(Global_ID, Replicate, Genotype,
                POSITION_X, POSITION_Y, Time) %>%
  dplyr::group_by(Global_ID) %>%
  dplyr::arrange(Global_ID, Time) %>%
  dplyr::mutate(
    xDis = POSITION_X - dplyr::lag(POSITION_X),
    yDis = POSITION_Y - dplyr::lag(POSITION_Y),
    Displacement = sqrt((xDis^2)+(yDis^2)),
    
    # Definition of  velocity vector - According to Michael Riedl
    xVec = xDis/(dplyr::lead(Time) - Time),
    yVec = yDis/(dplyr::lead(Time) - Time),
    
    #sumi = (xVec * dplyr::lag(xVec, 1)) + (yVec * dplyr::lag(yVec, 1))
    )

myVsum <- myVdata %>%
  dplyr::group_by(Global_ID, Replicate, Genotype) %>%
  dplyr::summarise(
    X_Pos = dplyr::first(POSITION_X),
    Y_Pos = dplyr::first(POSITION_Y),
    xavgDis = median(xVec, na.rm=T),
    yavgDis = median(yVec, na.rm=T),
    Loc = "Start"
  )

myVsum = myVsum %>%
  group_by(Global_ID, Replicate, Genotype) %>%
  dplyr::summarise(
    X_Pos = X_Pos+xavgDis,
    Y_Pos = Y_Pos+yavgDis,
    Loc = "End"
  ) %>%
  bind_rows(myVsum,.) %>%
  arrange(Global_ID)

ggplot(myVsum, aes(x=X_Pos, y=Y_Pos, color=Global_ID)) +
  geom_point(size = 3, alpha=0.1) +
  geom_path(size=1, arrow = arrow(type = "open", angle = 30, length = unit(0.1, "inches"))) +
  facet_wrap(~Genotype, scales = "free") +
  theme(legend.position = "none")
```

```