---
title: "Myo10 Velocity Analysis"
date: "2024"
authors: Julia Shangguan, edited by Ron Rock
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(purrr)
```

# Functions

## Signal to molecules

This function will take a df that includes `signal` and `set` columns. 

```{r}
sig2mol <- function(df, signal_var) {
  # Hardcoded constants all in one place here
  totalsigS1 <- 1139752131
  totalsigS2 <- 903216652
  totalsigS3 <- 951213825

  totalmolS1 <- 317063.2 * 56
  totalmolS2 <- 652125.8 * 54 
  totalmolS3 <- 619201.1 * 58

  correction <- 0.9 # percent of Halotag that TMR-Ligand labels
  
  df |>
    mutate(mol = ({{ signal_var }} / correction) * case_when(
      set == 1 ~ totalmolS1/totalsigS1,
      set == 2 ~ totalmolS2/totalsigS2,
      set == 3 ~ totalmolS3/totalsigS3
    ))
}
```

# Create a function that will process each set individually to get velocities per trajectory

* TrackMate analysis files include “Track segment merging” and “Track segment splitting” events
* That means, sometimes two spots with the same cell identifier & track ID exist within the SAME FRAME. 
* Therefore, time_diff = 0 rows occurs whenever there is the same frame number in a row bc they refer to multiple spots of the same trajectory existing in the same frame. This results in an "Inf" value for speed
* As a result, all merge and split events are excluded

For speeds:
* All speed saved as microns per frame
* Each frame is 300 ms exposure
* To get nm/s, the "speed" column is multiplied by 1000, then divided by 0.3

```{r}
process_sets <- function(set) { # input is the dataframe of a dataset
  set |>
    arrange(Cell, TRACK_ID, FRAME) |>  # Make sure frames are in ascending order for each TrackID's cell trajectory
    group_by(Cell, TRACK_ID) |>
    mutate(mol_diff = mol - first(mol),
           displacement = sqrt((POSITION_X - lag(POSITION_X))^2 + 
                               (POSITION_Y - lag(POSITION_Y))^2),
           normalizedFRAME = dense_rank(FRAME),
           rel_time = normalizedFRAME*0.3,
           time_diff = FRAME - lag(FRAME),
           speed = (displacement / time_diff) * 1000/0.3,
           velocity_x = (POSITION_X - lag(POSITION_X)) / time_diff,
           velocity_y = (POSITION_Y - lag(POSITION_Y)) / time_diff,
           velocity_firstx = (POSITION_X - first(POSITION_X)) / time_diff,
           velocity_firsty = (POSITION_Y - first(POSITION_Y)) / time_diff) |>
    filter(!any(is.infinite(speed))) |> # Toss *whole tracks* that have merge/splits
    ungroup()
}
```

# Distribution of molecules per cell

Sanity check for the sig2mol function:

```{r}
live_mol <- sig2mol(read_csv('231109_totsig_allsets.csv', show_col_types = FALSE), sig)
```

```{r}
# Found this function online to have exponential formatting
# Modified so that zero displays as "0"

scientific_10 <- function(x) {
  formatted <- format(x, scientific = TRUE)
  formatted[x == 0] <- "0"
  parse(text = gsub("e\\+*", " %*% 10^", formatted))
}

ggplot(live_mol , aes(x = mol)) + 
  geom_histogram(color="black", fill="gray") +
  labs(y = 'Count',
       x = "Number of molecules (per cell)") +
  theme_bw() +
    theme(axis.line = element_line(color='black'),
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()) +
  scale_x_continuous(
    breaks = pretty(live_mol$mol, n = 6), 
    expand = expansion(mult = c(0, 0)), 
    label=scientific_10
    ) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  theme(text = element_text(size = 25), axis.text = element_text(size = 16)) +
         theme(plot.margin = margin(20, 20, 20, 20)) +
  coord_cartesian(xlim = c(0, 3450000)) 
```

```{r}

# Upper and lower bound estimates of the values

(min(live_mol$mol))
(median(live_mol$mol))
(max(live_mol$mol))

# Error estimates
lowbound_slope <- 20.38978 
upperbound_slope <- 28.25411
slope_true <- 24.322

(min(live_mol$mol)* lowbound_slope/slope_true)
(median(live_mol$mol)* lowbound_slope/slope_true)
(max(live_mol$mol)* lowbound_slope/slope_true)

(min(live_mol$mol)* upperbound_slope/slope_true)
(median(live_mol$mol)* upperbound_slope/slope_true)
(max(live_mol$mol)* upperbound_slope/slope_true)
```

# Load spot movement datasets


```{r}
# Function to read in a set of experiments, and to attach the set number to the
# combined dataframe
# See R4DS 2nd ed,  26.3 (esp. 26.3.4)
read_sets <- function(filenames) {
  dfs <- filenames |>
    set_names(dirname) |>
    map(\(file) read_csv(file, show_col_types = FALSE)) |>
    bind_rows(.id = "set") |> #diff function than ron's
    mutate(set = parse_number(set),
           Cell = str_extract(`...1`, "Set(\\d+)Cell(\\d+)"))
}

# Now read each type
NspeedresultAll <- read_sets(
  c("Set1/240219Set1_nascentCells.csv", 
    "Set2/240219Set2_nascentCells.csv",
    "Set3/240219Set3_nascentCells.csv")
  ) |>
  sig2mol(TOTAL_INTENSITY_CH1) |>
  process_sets()

GspeedresultAll <- read_sets(
  c("Set1/240219Set1_growth_another_Cells.csv", 
    "Set2/240219Set2_growth_another_Cells.csv",
    "Set3/240219Set3_growth_another_Cells.csv")
  ) |>
  sig2mol(TOTAL_INTENSITY_CH1) |>
  process_sets()

RspeedresultAll <- read_sets(
  c("Set1/240219Set1_retrograde_Cells.csv", 
    "Set2/240219Set2_retrograde_Cells.csv",
    "Set3/240219Set3_retrograde_Cells.csv")
  ) |>
  sig2mol(TOTAL_INTENSITY_CH1) |>
  process_sets()
```

# Filtering spot movements that have strange behavior

Done here before any analysis / plotting, so everything uses the same sets

## Nascent filopodium:
	• Set 2, cell 1, 298
		○ 757210 --> looks like the spot glitches? Like a diff spot came into view at the same site
	• Set 3, cell 1, 1058
		○ 2596829 --> during time that filopodium curls and suddenly jerks
    ○ 2587619--> during time that filopodium curls and suddenly jerks

```{r}
NspeedresultAll_ed <- NspeedresultAll |> 
  filter(Cell != "nascent_Set2Cell1" | TRACK_ID != 298 | ID != 757210)
```

## Growth from existing (2nd phase) filopodium:

	• Set 2, cell 11, 2464 
		○ 2000213 --> not real signal; background noise being tracked
		○ 2011539 --> not real signal; background noise being tracked
	• Set 1, cell 1, 1602
		○ 5080436 --> looks like the spot glitches in trajectory; wrong thing temporally tracked
    ○ 5100633 --> looks like the spot glitches in trajectory; wrong thing temporally tracked

```{r}
GspeedresultAll_ed <- GspeedresultAll |> 
  filter(Cell != "growth_another_Set2Cell11" | TRACK_ID != 2464  | ID != 2000213) |>
  filter(Cell != "growth_another_Set2Cell11" | TRACK_ID != 2464  | ID != 2011539) |>
  filter(Cell != "growth_another_Set1Cell1" | TRACK_ID != 1602  | ID != 5080436 ) |>
  filter(Cell != "growth_another_Set1Cell1" | TRACK_ID != 1602  | ID != 5100633 ) 
```

## Returning trajectories:
	• Set 2, cell 1, 183
		○ 639643 --> seems fine
		○ 576594 --> seems fine
	• Set 2, cell 2 part 2, 49
		○ 159322 --> not real signal; background noise being tracked
		○ 34288 --> spot jumps, but seems like real signal
	• Set 2, cell 2 part 2, 22
		○ 38962 --> seems fine
    ○ 891589--> seems fine

```{r}
RspeedresultAll_ed <- RspeedresultAll |> 
  filter(Cell != "retrograde_Set2Cell2" | TRACK_ID != 49  | ID != 159322) 
```

## Number of cells

```{r}
length(unique(NspeedresultAll_ed$Cell))

length(unique(GspeedresultAll_ed$Cell))

length(unique(RspeedresultAll_ed$Cell))
```

## Track vs. cell listing

```{r}
NspeedresultAll_ed  |> select(Cell, TRACK_ID) |> unique()

GspeedresultAll_ed  |> select(Cell, TRACK_ID) |> unique()

RspeedresultAll_ed  |> select(Cell, TRACK_ID) |> unique()
```


# Starting molecules histograms (4A-C)

Start with a plotting function:
```{r}
plot_init_mols_histograms <- function(df) {
  init_mols <- df |>
    arrange(Cell, TRACK_ID, FRAME) |>
    group_by(Cell, TRACK_ID) |>
    slice(1:2) |> # calculate the average of the first two values
    summarise(mol = mean(mol))

  cat("\nminimum      = ", min(init_mols$mol), "\n")
  cat("median       = ", median(init_mols$mol), "\n")
  cat("maximum      = ", max(init_mols$mol), "\n\n")
  cat("minimum (LB) = ", (min(init_mols$mol))* lowbound_slope/slope_true, "\n")
  cat("median  (LB) = ", (median(init_mols$mol))* lowbound_slope/slope_true, "\n")
  cat("maximum (LB) = ", (max(init_mols$mol))* lowbound_slope/slope_true, "\n\n")
  cat("minimum (UB) = ", (min(init_mols$mol))* upperbound_slope/slope_true, "\n")
  cat("median  (UB) = ", (median(init_mols$mol))* upperbound_slope/slope_true, "\n")
  cat("maximum (UB) = ", (max(init_mols$mol))* upperbound_slope/slope_true, "\n")
  
  ggplot(init_mols, aes(x = mol)) + 
    geom_histogram(color="black", fill="gray") +
    labs(y = 'Count',
         x = "Myo10 molecules") +
    theme_bw() +
    theme(# axis.line = element_line(color='black'),
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()) +
    scale_x_continuous(# breaks = scales::pretty_breaks(n = 6), 
                       expand = expansion(mult = c(0.0, 0.0))) + 
    scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
    coord_cartesian(xlim = c(0, 1500)) #+
#    theme(text = element_text(size = 20), axis.text = element_text(size = 20)) +
#    theme(plot.margin = margin(20, 40, 20, 20))
}
```

## Nascent filopodium

```{r}
p_imol_n <- plot_init_mols_histograms(NspeedresultAll_ed)
p_imol_n
```

## Growth from existing (2nd phase) filopodium

```{r}
p_imol_g <- plot_init_mols_histograms(GspeedresultAll_ed)
p_imol_g
```

## Returning filopodium

```{r}
p_imol_r <- plot_init_mols_histograms(RspeedresultAll_ed)
p_imol_r
```


# Molecules vs. time traces 

Start with a plotting function:
```{r}
plot_mol_vs_time <- function(df) {
  ggplot(df, aes(x = rel_time, y = mol_diff)) + 
    geom_line(aes(color = as.factor(interaction(Cell, TRACK_ID))), size = 0.3, alpha = 0.3) +
    geom_smooth(method = "gam", size = 1, data = filter(df, rel_time < 100)) +
    labs(y = expression(atop("Myo10 molecules", paste("(", Delta, " from frame 1)"))) ,
         x = "Time (s)") +
    theme_bw() +
    theme(
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = "none" ) +
    scale_y_continuous(expand = expansion(mult = c(0, 0.05)))
}

```

## Nascent filopodium

```{r}
p_mvt_n <- plot_mol_vs_time(NspeedresultAll_ed)
p_mvt_n
```



## Growth from existing (2nd phase) filopodium

```{r}
p_mvt_g <- plot_mol_vs_time(GspeedresultAll_ed)
p_mvt_g
```

## Returning filopodium

```{r}
p_mvt_r <- plot_mol_vs_time(RspeedresultAll_ed)
p_mvt_r
```

## t-test on last point of trajectories

```{r}
NspeedresultAll_ed |>
  group_by(TRACK_ID, Cell) |>
  summarize(last_pt = last(mol_diff)) |>
  pull(last_pt) |>
  t.test()

GspeedresultAll_ed |>
  group_by(TRACK_ID, Cell) |>
  summarize(last_pt = last(mol_diff)) |>
  pull(last_pt) |>
  t.test()

RspeedresultAll_ed |>
  group_by(TRACK_ID, Cell) |>
  summarize(last_pt = last(mol_diff)) |>
  pull(last_pt) |>
  t.test()
```

# Spot speed vs. molecules plots 

Start with a plotting function:
```{r}
plot_speed_vs_mols <- function(df) { 
  cat("\nminimum = ", min(na.omit(df$speed)), " nm/s\n")
  cat("median  = ", median(na.omit(df$speed)), " nm/s\n")
  cat("maximum = ", max(na.omit(df$speed)), " nm/s\n\n")
  
  ggplot(df, aes(x = mol, y = speed)) + 
    geom_point(aes(color = interaction(Cell, TRACK_ID)), alpha = 0.1) +
    geom_smooth() +
    xlim(0, 4000) + ylim(0, 1000) +
    labs(y = 'Speed (nm/s)',
         x = "Myo10 molecules") +
    theme_bw() +
    theme(
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()) +
    guides(color = FALSE, size = FALSE) 
}
```

## Nascent filopodium

```{r}
summary(lm(speed ~ mol, data = NspeedresultAll_ed))
p_svm_n <- plot_speed_vs_mols(NspeedresultAll_ed)
p_svm_n
```

## Growth from existing (2nd phase) filopodium

```{r}
summary(lm(speed ~ mol, data = GspeedresultAll_ed))
p_svm_g <- plot_speed_vs_mols(GspeedresultAll_ed)
p_svm_g
```


## Returning filopodium

```{r}
summary(lm(speed ~ mol, data = RspeedresultAll_ed))
p_svm_r <- plot_speed_vs_mols(RspeedresultAll_ed)
p_svm_r
```

## Spearman's Rho test for speed vs. molecules

```{r}
cor.test(NspeedresultAll_ed$mol, NspeedresultAll$speed, exact=FALSE,
         method = "s", alternative = "less")

cor.test(GspeedresultAll_ed$mol, GspeedresultAll$speed, 
         method = "s", alternative = "less")

cor.test(RspeedresultAll_ed$mol, RspeedresultAll$speed, 
         method = "s", alternative = "less")
```

# How many trajectories analyzed in velocity plots?

```{r}
# Tracks
NspeedresultAll_ed %>% filter(is.na(speed) == FALSE) %>% select(Cell, TRACK_ID) %>% unique()

GspeedresultAll_ed %>% filter(is.na(speed) == FALSE) %>% select(Cell, TRACK_ID) %>% unique()

RspeedresultAll_ed %>% filter(is.na(speed) == FALSE) %>% select(Cell, TRACK_ID) %>% unique()

# Cells
NspeedresultAll_ed %>% filter(is.na(speed) == FALSE) %>% select(Cell) %>% unique()

GspeedresultAll_ed %>% filter(is.na(speed) == FALSE) %>% select(Cell) %>% unique()

RspeedresultAll_ed %>% filter(is.na(speed) == FALSE) %>% select(Cell) %>% unique()
```

# Final figure

```{r fig.width=10, fig.asp=0.8}
library(patchwork)

  # first column
  (p_imol_n + ggtitle("Initiation")) + 
  (p_imol_g + ggtitle("2nd Phase")) + 
  (p_imol_r + ggtitle("Retraction")) +
  # second column
  p_mvt_n + p_mvt_g + p_mvt_r +
  # third column
  p_svm_n + p_svm_g + p_svm_r +
  # some annotations that didn't look as good as the above
#  grid::textGrob('Initiation', rot = 90) +
#  grid::textGrob('2nd phase elongation', rot = 90) +
#  grid::textGrob('Retraction', rot = 90) +
  plot_layout(nrow = 3, byrow = FALSE) +
  plot_annotation(tag_levels = 'A')

#ggsave("Fig 4.pdf")
#ggsave("Fig 4.png")
```

