Mani_Timelapse_CF.Rmd

---
title: "Mani_Timelapse_CF"
author: "MANISHA"
date: "`r Sys.Date()`"
output: 
  html_document: 
    code_folding: hide
    fig_caption: yes
    fig_height: 6
    fig_width: 12
    fontsize: 10pt
    toc: true
    toc_float: true
    self_contained: false
---

# Load the config

In this file there is experiment that is analysed.
- Confocal imaging based.


This experiments was performed by Manisha in FV10i as first trial for live cell imaging.
apart.


```{r, eval=TRUE, message=FALSE, warning=FALSE}
library(magrittr)
library(dplyr)
library(ggplot2)
source("R/ggtheme_jack.R")
source("R/fun_multiplot.R")
ggplot2::theme_set(theme_jack())

layouts <- "TRUE"
cell_traces <- "TRUE"
fit_plots <- "TRUE"
Diagnosis <- "TRUE"
filtering <- "TRUE"
base_path <- c("../../../data")

exp_name <- c("20170528_THP1_WTBCG_PKNGBCG_LTR_TL_FV10_track_OUTPUT")

exp_read_out <- c("Cells")

dat <- readr::read_csv(sprintf("%s/%s/%s.csv",base_path, exp_name, exp_read_out))


dat <- dat[, !duplicated(colnames(dat))]



```

# Aim 

The experiment was DONE to determine the changes in the lysosomal number and intensities upon infection of WTBCG and PKNGBCG in THP1 cells with time.


There are three samples - uninfected, BCGinfected, PKNGinfected (Metadata_Inf_Status) 

- There are FIVE sites in each sample - 1,2,3 (Metadata_Sample)

- There are 79 timepoints for each site at  1 hour each (Metadata_Site)

- There are 3 Z slices for each site and each timepoint (Metadata_Time, Metadata_Zslice)

- The tracking id is for each cell (Trackobjects_Label_50)

We will sum up the values for z slices for each cell tracked.

Analysis Method -

1. We can fit the curves of each cell tracked - median fit or best fit
2. We can take the mean, median, stdev,mad etc for the lysosomal features in each condition and compare it.



## Join Layouts

```{r, message=FALSE, eval=layouts}

exp_layout <-
  readr::read_csv(
    sprintf("%s/%s/20170528_THP1_WTBCG_PKNGBCG_LTR_tl.csv", base_path, exp_name)
    ) %>%
  tbl_df()

dat$Metadata_Sample %<>% as.numeric()
dat$Metadata_Site %<>% as.numeric()

dat %<>%
  inner_join(exp_layout, .)

dat %<>%
  select(
    contains("Metadata"),
    contains("Children"),
    contains("AreaShape_Area"),
    contains("Label_50"),
    contains("Intensity_IntegratedIntensity"),
    contains("Object_Number")
    )

dat %<>% 
  select(
    -contains("Edge"), 
    -contains("Frame"), 
    -contains("Channel"),
    -contains("Series"),
    -contains("FileLocation")
    )


colnames(dat) <- gsub("_50", "", colnames(dat))

```

### Timepoints cleaning uup

```{r, message=FALSE, eval=filtering}

dat$Metadata_Time <- gsub("T001", "02", dat$Metadata_Time)
dat$Metadata_Time <- gsub("T002", "03", dat$Metadata_Time)
dat$Metadata_Time <- gsub("T003", "04", dat$Metadata_Time)
dat$Metadata_Time <- gsub("T004", "05", dat$Metadata_Time)
dat$Metadata_Time <- gsub("T005", "06", dat$Metadata_Time)
dat$Metadata_Time <- gsub("T006", "07", dat$Metadata_Time)
dat$Metadata_Time <- gsub("T007", "08", dat$Metadata_Time)
dat$Metadata_Time <- gsub("T008", "09", dat$Metadata_Time)
dat$Metadata_Time <- gsub("T009", "10", dat$Metadata_Time)
dat$Metadata_Time <- gsub("T010", "11", dat$Metadata_Time)
dat$Metadata_Time <- gsub("T011", "12", dat$Metadata_Time)

dat$Metadata_Time %<>% as.numeric()

dat %<>% filter(Metadata_Time < 13)

```



# Diagnostic plots without filtering

```{r, fig.width=10, fig.height=8, eval=Diagnosis}


p1 <- dat %>% ggplot(aes(Mean_Nuclei_AreaShape_Area)) + geom_density() + xlim(0,3500)+
  geom_vline(aes(xintercept=100)) +
  geom_vline(aes(xintercept=3250))
p2 <- dat %>% ggplot(aes(Intensity_IntegratedIntensity_origDAPI)) + geom_density() + xlim(0,1800)+ 
  geom_vline(aes(xintercept=200)) + 
  geom_vline(aes(xintercept=1200))
p3 <- dat %>% ggplot(aes(Children_BCG_Count))+geom_density()+ xlim(0,10)
p4 <- dat %>% ggplot(aes(Mean_BCG_AreaShape_Area))+geom_density()+ xlim(0,95)
p5 <- dat %>% ggplot(aes(Intensity_IntegratedIntensity_FilteredGCH))+geom_density()+ xlim(0,500)

p6 <- dat %>% ggplot(aes(Children_RCH_Count))+geom_density() + xlim(0,80)
p7 <- dat %>% ggplot(aes(Mean_RCH_AreaShape_Area)) + geom_density() + xlim(0,80)
p8 <- dat %>% ggplot(aes(Intensity_IntegratedIntensity_FilteredRCH)) + geom_density() + xlim(0,1000)
p9 <- dat %>% ggplot(aes(Intensity_IntegratedIntensity_origRCH))+geom_density()+ xlim(0,1500)

multiplot(p1,p2,p3,p4,p5,p6,p7,p8,p9, cols = 3)
```

# Filtering

```{r, message=FALSE, fig.width=5, fig.height=4, eval=filtering}

n1 <- nrow(dat)

dat %<>% 
  filter(Intensity_IntegratedIntensity_origDAPI > 200 & Intensity_IntegratedIntensity_origDAPI < 1200)
n2 <- nrow(dat)

dat %<>% 
  filter(Mean_Nuclei_AreaShape_Area < 3250 & Mean_Nuclei_AreaShape_Area > 100)

n3 <- nrow(dat)


n1
n2
n3

cells_filtered <- (n1-n3)*100/n1

sprintf("Percentage of cells filtered out: %s", cells_filtered)

```

# Diagnostic plots after filtering

```{r, fig.width=10, fig.height=8, eval=Diagnosis}
p1 <- dat %>% ggplot(aes(Mean_Nuclei_AreaShape_Area)) + geom_density() + xlim(0,3500)+
  geom_vline(aes(xintercept=100)) +
  geom_vline(aes(xintercept=3250))
p2 <- dat %>% ggplot(aes(Intensity_IntegratedIntensity_origDAPI)) + geom_density() + xlim(0,1800)+ 
  geom_vline(aes(xintercept=200)) + 
  geom_vline(aes(xintercept=1200))
p3 <- dat %>% ggplot(aes(Children_BCG_Count))+geom_density()+ xlim(0,10)
p4 <- dat %>% ggplot(aes(Mean_BCG_AreaShape_Area))+geom_density()+ xlim(0,95)
p5 <- dat %>% ggplot(aes(Intensity_IntegratedIntensity_FilteredGCH))+geom_density()+ xlim(0,500)

p6 <- dat %>% ggplot(aes(Children_RCH_Count))+geom_density() + xlim(0,80)
p7 <- dat %>% ggplot(aes(Mean_RCH_AreaShape_Area)) + geom_density() + xlim(0,80)
p8 <- dat %>% ggplot(aes(Intensity_IntegratedIntensity_FilteredRCH)) + geom_density() + xlim(0,1000)
p9 <- dat %>% ggplot(aes(Intensity_IntegratedIntensity_origRCH))+geom_density()+ xlim(0,1500)

multiplot(p1,p2,p3,p4,p5,p6,p7,p8,p9, cols = 3)
```

# labelling the infection status

```{r,  message=FALSE, eval=layouts}

unique(dat$Metadata_InfStatus)

dat$Metadata_InfStatus <- gsub("UNINF", "NI", dat$Metadata_InfStatus)

unique(dat$Metadata_InfStatus)

dat$Metadata_InfStatus[dat$Children_BCG_Count < "1" & dat$Metadata_InfStatus == "PKNG"] <- "Byst_PKNG"
dat$Metadata_InfStatus[dat$Children_BCG_Count < "1" & dat$Metadata_InfStatus == "WTBCG"] <- "Byst_WTBCG"

unique(dat$Metadata_InfStatus)

dat$Metadata_InfStatus <-
  factor(dat$Metadata_InfStatus, levels=c('NI','Byst_WTBCG','WTBCG', 'Byst_PKNG' , 'PKNG'))

```

# Grouping the zslices

```{r,  message=FALSE, eval=layouts}

meta_cols <- dat %>% select(contains("Metadata")) %>% names()
feat_cols <- dat %>% select(-contains("Metadata"), -TrackObjects_Label) %>% names()

dat_mean <-
  dat %>%
  group_by(
    Metadata_Sample,
    Metadata_Site,
    Metadata_Time,
    Metadata_InfStatus,
    TrackObjects_Label) %>%
  dplyr::summarize_each_(funs(mean(., na.rm=TRUE)),
                         vars = feat_cols)

dat_sum <-
  dat %>%
  dplyr::group_by(
    Metadata_Sample,
    Metadata_Site,
    Metadata_Time,
    Metadata_InfStatus,
    TrackObjects_Label) %>%
  dplyr::summarize_each_(dplyr::funs(sum(., na.rm=TRUE)),
                         vars = feat_cols)

```

# Filtering Single track label cells

```{r,  message=FALSE, eval=layouts}


fat <- 
  dat %>% 
  group_by( 
  Metadata_Sample,
  Metadata_Site,
  Metadata_Time,
  Metadata_InfStatus,
  TrackObjects_Label
  ) %>% count() %>%
  filter(n == 1) %>%
  select(-n)

dat1 <- inner_join(fat,dat)

```




##Plotting data

# Single cell traces

Imaging done per 1 hour for 12 hours

```{r, fig.width=12, fig.height=16, warning=FALSE, eval=cell_traces}

ggplot(dat1,
       aes(
         x = Metadata_Time,
         y = Intensity_IntegratedIntensity_origRCH,
         color = Metadata_Zslice,
         group = TrackObjects_Label)) +
  geom_line(alpha=0.5) + theme(legend.position="none") +
  facet_grid(Metadata_Site + Metadata_Zslice ~ Metadata_InfStatus)

```



# Ploting changes in lysotracker content in cell over time

```{r, message=FALSE, warning=FALSE, fig.width=8, fig.height=8, eval=TRUE}
feat_cols <- c("Intensity_IntegratedIntensity_origRCH")

dat_med <-
  dat1 %>%
  dplyr::group_by(
    Metadata_Time,
    Metadata_InfStatus) %>%
  dplyr::summarize_each_(dplyr::funs(mean(., na.rm=TRUE)),
                         vars = feat_cols)

p1 <-
  ggplot(dat1,
         aes(
           x = Metadata_Time,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Fitted Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Time (in hours)")
p2 <-
  ggplot(dat_med,
         aes(
           x = Metadata_Time,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Median Collapsed Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Time (in hours)")

multiplot(p1, p2, cols=1)

######

feat_cols <- c("Intensity_IntegratedIntensity_FilteredRCH")



dat_med <-
  dat1 %>%
  dplyr::group_by(
    Metadata_Time,
    Metadata_InfStatus) %>%
  dplyr::summarize_each_(dplyr::funs(mean(., na.rm=TRUE)),
                         vars = feat_cols)

p1 <-
  ggplot(dat1,
         aes(
           x = Metadata_Time,
           y = (Intensity_IntegratedIntensity_FilteredRCH),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Fitted Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Time (in hours)")
p2 <-
  ggplot(dat_med,
         aes(
           x = Metadata_Time,
           y = (Intensity_IntegratedIntensity_FilteredRCH),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Median Collapsed Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Time (in hours)")

multiplot(p1, p2, cols=1)

######

feat_cols <- c("Children_RCH_Count")

dat_med <-
  dat1 %>%
  dplyr::group_by(
    Metadata_Time,
    Metadata_InfStatus) %>%
  dplyr::summarize_each_(dplyr::funs(mean(., na.rm=TRUE)),
                         vars = feat_cols)

p1 <-
  ggplot(dat1,
         aes(
           x = Metadata_Time,
           y = (Children_RCH_Count),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Fitted Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Count of lysotracker",
    x="Time (in hours)")
p2 <-
  ggplot(dat_med,
         aes(
           x = Metadata_Time,
           y = (Children_RCH_Count),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Median Collapsed Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Count of lysotracker",
    x="Time (in hours)")

multiplot(p1, p2, cols=1)




```



# Ploting changes in lysotracker content in cell over time without unexposed

```{r, message=FALSE, warning=FALSE, fig.width=8, fig.height=8, eval=TRUE}
feat_cols <- c("Intensity_IntegratedIntensity_origRCH")

dat_med <-
  dat1 %>%
  dplyr::group_by(
    Metadata_Time,
    Metadata_InfStatus,
    Metadata_Zslice) %>%
  dplyr::summarize_each_(dplyr::funs(mean(., na.rm=TRUE)),
                         vars = feat_cols)

p1 <-
  ggplot(dat1 %>% filter(Metadata_InfStatus !="NI"),
         aes(
           x = Metadata_Time,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) +
  geom_smooth() +
  labs(
    title="Fitted Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Time (in hours)")
p2 <-
  ggplot(dat_med %>% filter(Metadata_InfStatus !="NI"),
         aes(
           x = Metadata_Time,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Median Collapsed Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Time (in hours)")

multiplot(p1, p2, cols=1)


######

feat_cols <- c("Intensity_IntegratedIntensity_FilteredRCH")



dat_med <-
  dat1 %>%
  dplyr::group_by(
    Metadata_Time,
    Metadata_InfStatus) %>%
  dplyr::summarize_each_(dplyr::funs(mean(., na.rm=TRUE)),
                         vars = feat_cols)

p1 <-
  ggplot(dat1 %>% filter(Metadata_InfStatus !="NI"),
         aes(
           x = Metadata_Time,
           y = (Intensity_IntegratedIntensity_FilteredRCH),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Fitted Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Time (in hours)")
p2 <-
  ggplot(dat_med %>% filter(Metadata_InfStatus !="NI"),
         aes(
           x = Metadata_Time,
           y = (Intensity_IntegratedIntensity_FilteredRCH),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Median Collapsed Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Time (in hours)")

multiplot(p1, p2, cols=1)

######

feat_cols <- c("Children_RCH_Count")

dat_med <-
  dat1 %>%
  dplyr::group_by(
    Metadata_Time,
    Metadata_InfStatus) %>%
  dplyr::summarize_each_(dplyr::funs(mean(., na.rm=TRUE)),
                         vars = feat_cols)

p1 <-
  ggplot(dat1 %>% filter(Metadata_InfStatus !="NI"),
         aes(
           x = Metadata_Time,
           y = (Children_RCH_Count),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Fitted Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Count of lysotracker",
    x="Time (in hours)")
p2 <-
  ggplot(dat_med %>% filter(Metadata_InfStatus !="NI"),
         aes(
           x = Metadata_Time,
           y = (Children_RCH_Count),
           color = Metadata_InfStatus)) +
  geom_smooth()+
  labs(
    title="Median Collapsed Kinetics of lysotracker after 0 hpi",
    y="Whole Cell Int-Count of lysotracker",
    x="Time (in hours)")

multiplot(p1, p2, cols=1)


```


# Ploting changes in lysotracker content in cell at time intervals of hrs as boxplot
```{r, message=FALSE, warning=FALSE, fig.width=14, fig.height=9, eval=TRUE}

feat_cols <- c("Intensity_IntegratedIntensity_origRCH")



time_point <- c("2" , "4", "6", "8", "10", "12")

n1 <- ggplot(dat1 %>% filter(Metadata_InfStatus !="NI" && Metadata_Time == time_point[1]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 2 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")


n2 <- ggplot(dat1 %>% filter(Metadata_InfStatus !="NI" && Metadata_Time == time_point[2]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 4 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")

n3 <- ggplot(dat1 %>% filter(Metadata_InfStatus !="NI" && Metadata_Time == time_point[3]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 6 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")

n4 <- ggplot(dat1 %>% filter(Metadata_InfStatus !="NI" && Metadata_Time == time_point[4]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 8 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")

n5 <- ggplot(dat1 %>% filter(Metadata_InfStatus !="NI" && Metadata_Time == time_point[5]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 10 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")

n6 <- ggplot(dat1 %>% filter(Metadata_InfStatus !="NI" && Metadata_Time == time_point[6]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 12 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")



multiplot(n1,n2,n3,n4,n5,n6, cols=3)

```


# Ploting changes in lysotracker content in cell at time intervals of hrs as boxplot
```{r, message=FALSE, warning=FALSE, fig.width=14, fig.height=9, eval=TRUE}

feat_cols <- c("Intensity_IntegratedIntensity_origRCH")



time_point <- c("2" , "4", "6", "8", "10", "12")

n1 <- ggplot(dat1 %>% filter(Metadata_Time == time_point[1]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 2 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")


n2 <- ggplot(dat1 %>% filter(Metadata_Time == time_point[2]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 4 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")

n3 <- ggplot(dat1 %>% filter(Metadata_Time == time_point[3]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 6 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")

n4 <- ggplot(dat1 %>% filter(Metadata_Time == time_point[4]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 8 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")

n5 <- ggplot(dat1 %>% filter(Metadata_Time == time_point[5]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 10 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")

n6 <- ggplot(dat1 %>% filter(Metadata_Time == time_point[6]),
         aes(
           x = Metadata_InfStatus,
           y = (Intensity_IntegratedIntensity_origRCH),
           color = Metadata_InfStatus)) + ylim(150,1250) +
  geom_boxplot()+
  facet_grid( ~ Metadata_Time) +
  labs(
    title="lysotracker intensity at 12 hpi",
    y="Whole Cell Int-Intensity of lysotracker",
    x="Infection status at Time (in hours)")


multiplot(n1,n2,n3,n4,n5,n6, cols=3)

```