2_plot_data.Rmd

---
title: "Untitled"
output: html_document
editor_options: 
  markdown: 
    wrap: 72
---

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

```{r eval=FALSE, include=FALSE}
records <- list.files("./inst/extdata/")
records
```

```{r}
for(library_name in c(
  'tidyverse',
  'janitor',
  'broom',
  'cowplot',
  'ggh4x',
  'ggsignif',  
  'PNWColors',
  'agricolae'
  )){
  if(FALSE == require(library_name, character.only = TRUE)){
    print(paste0("installing '", library_name, "'"))
    install.packages(library_name)
  } 
  library(library_name, character.only = TRUE) 
}
```

## Load in metadata

The referenced files match naming of the xlsx file in teams. To have
plain text files rather than binary the sheets are stored as individual
csvs in `seeds in rootbot/`.

```{r}
if(!file.exists("./data/metadata.rds")){
#after drafting switch to above line
# if(TRUE){
  metadata <- read.csv("./inst/extdata/seeds in rootbot/Log.csv")
  # remove human readable but not computationally practical columns of text.
  metadata <- metadata[, c("X282.inbred",
                           "Rootbot.date", 
                           "scored", 
                           "source...200k.packets.from", 
                           "WW.cassette.", 
                           "WS.cassette.",
                           "Earliest_scored_plate")]
  
  metadata <- metadata %>% 
    rename("inbred" = "X282.inbred",
           "rootbot_date" = "Rootbot.date", 
           "scored" = "scored", 
           "source" = "source...200k.packets.from", 
           "ww_cassette" = "WW.cassette.", 
           "ws_cassette" = "WS.cassette.",
           "earliest_scored_plate" = "Earliest_scored_plate") 
  
  saveRDS(metadata, file = "./data/metadata.rds")
}else{
  metadata <- readRDS(file = "./data/metadata.rds")
}
```

## Load in records

```{r}
# if the aggregated dataframe (df.rds) exists load it otherwise build and save it
if(!file.exists("./data/df.rds")){
#TODO after drafting switch to above line
# if(TRUE){
  records <- list.files("./inst/extdata/plate_csvs")
  
  record_list <- list()
  for(record in records){
    # print(record) # for debugging
    entry <- read.csv(paste("./inst/extdata/plate_csvs", record, sep = "/"))
    entry["file"] <- record
    
    record_list[[length(record_list)+1]] <- entry
  }
  
  # sloppy and inefficent but functional.
  agg_df <- record_list[[1]]
  for(i in seq(2, length(record_list))){
    # agg_df <- full_join(agg_df, record_list[[i]])
    agg_df <- full_join(
      mutate(agg_df, 
             parent_name = as.character(parent_name),
             parent = as.character(parent)
             ),
      mutate(record_list[[i]], 
             parent_name = as.character(parent_name),
             parent = as.character(parent)
             )
    )
  }
  
  df <- janitor::clean_names(agg_df)
  # retain only primary root (should also rm int named lateral roots)
  df <- df[df$root_ontology == " Primary root", ]
  # discard .+behind
  df <- df[!(stringr::str_detect(df$root_name, '.+behind*')), ]
  # discard C49Adon't_use_   
  df <- df[!(stringr::str_detect(df$root_name, '.+don\'t')), ]
  
  df$root_name <- stringr::str_replace(df$root_name, 'seed', '')
  df$root_name <- stringr::str_replace(df$root_name, ' ', '')
  # Drop control seeds from experimental plates
  df <- df[!(stringr::str_detect(df$root_name, 'B73w.+')), ]
  df <- df[!(df$root_name %in% c('B73_WScml158q', 'B73M162W',  'B74wwA554')), ]

  df <- df[df$root_name != "", ]
  df$root_name <- stringr::str_replace(df$root_name, ' ', '')
  df$root_name <- stringr::str_replace(df$root_name, '_', '')
  
  
  # fix misnamed values  
  # Note: B73 plates appear to lack root numbers these will need to be imputed
  # or non-pairwise analysis will need to be used.
  df[df$root_name == 'C1231', 'root_name'] <- 'C123#1'
  df[df$root_name == 'B73', 'root_name'] <- 'B73#-9999'
  
  # split root_name into genotype/rep
  df <- separate(
    df,
    root_name,
    c("genotype", "replicate"),
    sep = "#",
    remove = TRUE,
    convert = FALSE,
    extra = "warn",
    fill = "warn"
  )
  
  # strip out leading space
  df$genotype <- str_remove_all(df$genotype, " ")
  
  saveRDS(df, file = "./data/df.rds")  
  
}else{
  df <- readRDS(file = "./data/df.rds")
}
```


```{r}
metadata$inbred <- toupper(metadata$inbred)
df$genotype <- toupper(df$genotype)

mismatched_inbreds <- full_join(
  data.frame(
    inbred = unique(metadata$inbred),
    metadata = T),
  data.frame(
    inbred = unique(df$genotype),
    df = T))

mismatched_inbreds[is.na(mismatched_inbreds$metadata), 'metadata'] <- F
mismatched_inbreds[is.na(mismatched_inbreds$df), 'df'] <- F

mismatched_inbreds <- mismatched_inbreds %>% 
  mutate(match = case_when(metadata != df ~ F, 
                           metadata == df ~ T)) %>% 
  arrange(match, inbred)

write.csv(
  mismatched_inbreds, 
  "./missing_inbreds_not_in_data.csv")
```


## tidy up metadata

```{r}
metadata <- metadata %>% as_tibble()

metadata <- metadata %>% 
  # ignore non scored entries
  mutate(
    scored = case_when(
      scored == "yes" ~ TRUE,
      scored == "yes " ~ TRUE,
      scored == "no" ~ FALSE,
      scored == ""~ FALSE)
    ) %>% 
  filter(scored == TRUE)
  

metadata <- metadata %>% 
  # drop missing dates
  filter(!(rootbot_date %in% c("", " "))) %>% 
  # make into dates
  mutate(rootbot_date = lubridate::mdy(rootbot_date)) 
  
  
# Make tidy wrt cassette
metadata <- metadata %>% 
  pivot_longer(c("ww_cassette", "ws_cassette")) %>% 
  rename("condition" = "name",
         "cassette" = "value") %>% 
  mutate(condition = case_when(
    condition == "ww_cassette" ~ "ww",
    condition == "ws_cassette" ~ "ws"
  ))


metadata <- metadata %>% distinct()

metadata <- metadata %>% 
  mutate(uid = seq(1, nrow(metadata)))
```

## Tidy up data

```{r}
df <- as_tibble(df)

df["image_date"] <- stringr::str_extract(df$image, "\\d+") %>% lubridate::ymd()

df["image_datetime"] <- stringr::str_extract(df$image, "\\d+-\\d+") %>% lubridate::ymd_hms()

df["image_cassette"] <- stringr::str_extract(df$image, "\\d+$") %>% as.integer()
```


```{r}
df <- df %>% rename(
  inbred = genotype,
  cassette = image_cassette
)
df["in_df"] = T
metadata["in_meta"] = T

data <- full_join(metadata, df)

data[is.na(data$in_df), 'in_df'] = F
data[is.na(data$in_meta), 'in_meta'] = F

data %>% 
  group_by(in_df, in_meta) %>% 
  tally()
  
obs_w_missing_uid <- data %>% 
  filter(in_meta == F) %>% 
  select(-length, -replicate, -image_date, -image_datetime) %>% 
  distinct()
  
obs_w_missing_uid
write.csv(obs_w_missing_uid, './obs_w_missing_uid.csv')
```


```{r}
library(stringdist)
obs_w_missing_uid <- obs_w_missing_uid[, c("inbred", "image", "file", "in_meta", "in_df")] %>% 
  distinct()




inbred_google <- function(
    search_inbred = "RS229",
    n_matches = 10,
    possible_inbreds = unique(metadata$inbred)
){
 string_distances <- stringdist(search_inbred, possible_inbreds)

  distance_table <- 
  data.frame(
    inbred = possible_inbreds,
    dist   = string_distances)
  
  distance_table <- distance_table %>% 
    arrange(dist)
  
  best_matches <- distance_table[seq(1, n_matches), 
                 'inbred'] 
  return(best_matches)
}


n_matches = 10
obs_w_missing_uid %>% head()

obs_w_missing_uid$inbred


inbred_google(
  search_inbred = "RS229",
  n_matches = n_matches,
  possible_inbreds = unique(metadata$inbred)
)

obs_w_missing_uid[, ]

paste0("match_", as.character(1:n_matches))

# mask <- string_distances == min(string_distances)
# possible_inbreds[mask]




# ## Steps
# 
# If value in data but not in log file
# 
# for each missing value
# 
# 1.  Look for similar names in "seeds in rootbot.xlsx" Long
# 2.  Look at the photo and see if the name in that file is consistent
#     1.  Get genotypes for the same experiment to help locate correct packet
# 3.  Check downstairs to see if there is a typo in the packet
# 
# Example corrections:
# 
# tz18 == tzi8
# 
# gf12 == gt12
# 
# If value in log file but not in data
# 
# 1.  Look for names that are similar (something that looks like a typo could be another inbred line)

```


```{r}
# TODO ideally there should be no mismatched keys to be dropped
data <- data[(data$in_meta & data$in_df), ] %>% select(-in_df, -in_meta)

```


```{r}
# Group experiments by time
hours_threshold = 36


temp = data.frame(timestamp = sort(unique(c(data$rootbot_date, data$image_date))))

temp['next_timestamp'] <- c(as.Date('1999-01-01'), temp$timestamp[1:(nrow(temp)-1)])
temp['exp_group'] <- NA

ith_exp = 0
for(i in seq(1, nrow(temp))){
  if(difftime(temp[i, 'timestamp'], 
              temp[i, 'next_timestamp'], units = 'hours') > hours_threshold){
    ith_exp <- ith_exp +1
  }  
  temp[i, 'exp_group'] <- ith_exp
}

temp <- temp %>% select(-next_timestamp)

data <- full_join(data, rename(temp, image_date = timestamp))


temp <- data %>% 
  group_by(inbred, condition, cassette, uid, exp_group) %>% 
  summarise(
    image_datetime = image_datetime,
    min_datetime = min(image_datetime)) %>% 
  ungroup() %>% 
  distinct() %>% 
  mutate(delta_datetime = image_datetime - min_datetime) %>% 
  mutate(hours_elapsed = as.numeric(delta_datetime/(60*60)))

data <- full_join(data, temp)

```

```{r}
# apply time adjustments
# use Earliest_scored_plate column in log file
# todo anything without col starting at *** assume 6, 36
data[is.na(data$earliest_scored_plate), "earliest_scored_plate"] <- 0
data$hours_elapsed <- data$hours_elapsed + data$earliest_scored_plate

data
write.csv(data, './data.csv')
```


```{r}
data_summary_stats <- data %>% 
  select(length, inbred, condition, hours_elapsed) %>% 
  group_by(inbred, condition, hours_elapsed) %>% 
  summarise(min = min(   length, na.rm = TRUE),
            med = median(length, na.rm = TRUE),
            mean= mean(  length, na.rm = TRUE),
            max = max(   length, na.rm = TRUE),
            sd  = sd(    length, na.rm = TRUE),
            iqr = IQR(   length, na.rm = TRUE),
            )
write.csv(data_summary_stats, './data_summary_stats.csv')

# TODO replace with a more robust model (rank deficient) 
fm <- lm("length ~ hours_elapsed*condition*inbred", data)
res <- broom::tidy(fm)
res
# set up prediction df
temp <- data %>% select(inbred, condition) %>% distinct()
temp <- rbind(mutate(temp, hours_elapsed = 6),
              mutate(temp, hours_elapsed = 36))
temp['length_yHat'] <- predict(fm, temp)

write.csv(temp, './data_summary_preds.csv')
```


## merge
```{r}
# quick confirmation that the merge worked as expected
ggplot(data, aes(x = rootbot_date))+geom_density()+geom_point(aes(y=0))+theme_minimal()

ggplot(data, aes(x = hours_elapsed))+geom_density()+geom_point(aes(y=0))+theme_minimal()
```


```{r}
data %>% 
  filter(inbred == "C123") %>% 
  mutate(delta_datetime = delta_datetime/(60*60)) %>% # rescale from seconds to hours
  ggplot(aes(delta_datetime, length, color = condition))+
  geom_point()+
  geom_smooth(method = "lm", aes(fill = condition), alpha = 0.3)+
  # facet_grid(.~condition)
  labs(x = "Hours Elapsed", y = "Root Length (mm)", title = "Genotype: C123")
  
```

```{r}
# where ws is longer than ww
# slope or abs?


M <- data %>% 
  filter(inbred == "C123") %>% 
  mutate(delta_datetime = delta_datetime/(60*60)) 


fm <- lm("length ~ delta_datetime*condition", M)
res <- broom::tidy(fm)

# is there a significant difference?
res[res$term == "delta_datetime:conditionww", "p.value"] < 0.05

# ws has a more positive slope than ww
res[res$term == "delta_datetime:conditionww", "estimate"] < 0
```

```{r}
# find all inbreds without both conditions:
rm_geno <- data %>% 
  group_by(inbred) %>% 
  select(condition) %>% 
  distinct() %>% 
  tally() %>% 
  filter(n <2)
  
data <- data[!(data$inbred %in% rm_geno$inbred), ]


res_list <- purrr::map(unique(data$inbred), function(e){
  M <- data %>% 
  filter(inbred == e) %>% 
  mutate(delta_datetime = delta_datetime/(60*60)) 

  fm <- lm("length ~ delta_datetime*condition", M)
  res <- broom::tidy(fm)
  
})

res_list_sig <- purrr::map(res_list, function(res){
  # is there a significant difference?
  unlist(res[res$term == "delta_datetime:conditionww", "p.value"] )
})


res_list_est <- purrr::map(res_list, function(res){
  # ws has a more positive slope than ww
  unlist(res[res$term == "delta_datetime:conditionww", "estimate"])
})

res_overview <- data.frame(
  list_idx = seq(1, length(res_list)),
  p_value = unlist(res_list_sig),
  interact = unlist(res_list_est)
)
```

```{r}
res_overview %>%
  filter(!is.na(p_value)) %>% 
  # filter(p_value < 0.05) %>% 
  group_by(p_value < 0.05, interact < 0) %>% 
  tally()
```

# visualize those with ws favored

```{r}
list_idxs <- res_overview %>% filter(#p_value < 0.05, 
                                     interact < 0) %>% select(list_idx) %>% unlist()

i = 1
res_list[[list_idxs[i]]]


unique(data$inbred)[list_idxs[i]]


plts <- purrr::map(list_idxs,
  # c(15, #17, 18, 
  #   22, 32, 54, 
  #   #83, 
  #   99, 106),
                   function(ii){
  geno <- unique(data$inbred)[ii]
  
  data %>% 
  filter(inbred == geno) %>% 
  mutate(delta_datetime = delta_datetime/(60*60)) %>% # rescale from seconds to hours
  mutate(condition = factor(condition, levels = c("ww", "ws"))) %>% 
  ggplot(aes(delta_datetime, length, color = condition))+
  geom_point()+
  geom_smooth(method = "lm", aes(fill = condition), alpha = 0.3)+
  # facet_grid(.~condition)
  labs(x = "Hours Elapsed", y = "Root Length (cm)", title = paste0("Genotype: ", geno))+
  theme_bw()+
  theme(legend.position = "bottom")+
  ggsci::scale_color_lancet()+
  ggsci::scale_fill_lancet()
})


cowplot::plot_grid(plotlist = plts)


```

```{r eval=FALSE, include=FALSE}
# we could also consider an interaction model (but really we should used a mixed model)

fm <- lm("length ~ delta_datetime*condition*inbred", data)
res <- broom::tidy(fm)

# library(lme4)

# lmer(length ~ genotype*group , random = ~1|rootid)

```

![Example of Genotype Plot](reports/figures/33-16.png)

```{r eval=FALSE}

for(current_inbred in unique(data$inbred)){

  temp <- data[data$inbred == current_inbred, ]
  
  # tweak group ordering for plotting
  temp$condition <- factor(temp$condition, c('ww', 'ws'))
  temp
  
  
  
  if( length(unique(temp$delta_datetime)) == 1 ){
  
    
    fm_ph <- agricolae::HSD.test(fm, trt = c('condition'))
    fm_ph <- fm_ph$groups
    fm_ph$comparison <- rownames(fm_ph)
    fm_ph$comparison <- stringr::str_replace(fm_ph$comparison, ":", ".")
    
  } else if( length(unique(temp$condition)) == 1 ){
  
    
    fm_ph <- agricolae::HSD.test(fm, trt = c('delta_datetime'))
    fm_ph <- fm_ph$groups
    fm_ph$comparison <- rownames(fm_ph)
    fm_ph$comparison <- stringr::str_replace(fm_ph$comparison, ":", ".")
    
  } else{
    fm <- lm(length ~ delta_datetime*condition, temp)
    # annotations
    fm_tab <- broom::tidy(fm)
    fm_tab    
    
    fm_ph <- agricolae::HSD.test(fm, trt = c('delta_datetime', 'condition'))
    fm_ph <- fm_ph$groups
    fm_ph$comparison <- rownames(fm_ph)
    fm_ph$comparison <- stringr::str_replace(fm_ph$comparison, ":", ".")
  }
  
  # set up base of the figure. Pass in aes to each geom so that we can use 
  # different columns for annotation
  fig <- temp %>% 
    ggplot()+
    ggdist::stat_halfeye(
      aes(interaction(delta_datetime, condition), length, fill = condition),
      adjust = .5, 
      width = .6, 
      .width = 0, 
      justification = -.2, 
      point_colour = NA
    ) + 
    geom_boxplot(
      aes(interaction(delta_datetime, condition), length, fill = condition),
      width = .15, 
      outlier.shape = NA
    ) +
    ## add justified jitter from the {gghalves} package
    gghalves::geom_half_point(
          aes(interaction(delta_datetime, condition), length, fill = condition),
          side = "l", 
      range_scale = .4, 
      alpha = .3
    ) +
    scale_x_discrete(guide = "axis_nested")+
    xlab("")+
    ylab("Length (mm)")+
    labs(title = paste('inbred: ', current_inbred))+
    ggthemes::theme_clean()+
    ggthemes::scale_fill_few()+
    theme(legend.position = 'bottom')
  
  
  
  # check if there's more than one condition and add in the relevant comparisons 
  if( !((length(unique(temp$condition)) == 1) & (length(unique(temp$delta_datetime)) == 1)) ){
  
    y_pos_start <- max(temp$length)
    y_pos_step <- 1 #y_pos_start / 10
    # select each hsd condition. Then for each letter in in the conditioning code add in a 
    # label with a letter specific y axis
    for(i in seq_along(unique(fm_ph$comparison)) ){
      fm_ph_subset = fm_ph[fm_ph$comparison == unique(fm_ph$comparison)[i] , ]
      
      # this is a messy way to get the comparisons at different locations.
      hsd_letters <- strsplit(fm_ph_subset$groups, split = "")[[1]]
      for(hsd_letter in hsd_letters){
        y_increment = seq_along(letters)[hsd_letter == letters]
  
        
        dat = data.frame(
            x=fm_ph_subset$comparison,
            y=c(y_pos_start+(y_pos_step*(y_increment))),
            annotation=toupper(fm_ph_subset$groups))
  
        fig <- fig+geom_label(
          data=dat,
          aes(x=x, y=y, label=annotation),
          fontface="bold",
          fill = pnw_palette("Sailboat",4,type="continuous")[y_increment],
          color = pnw_palette("Sailboat",4,type="continuous")[y_increment]
          )
  
        fig <- fig+geom_text(
          data=dat,
          aes(x=x, y=y, label=annotation),
          fontface="bold"
          )
      }
    }    
  }
  ggsave(paste0('./reports/figures/',current_inbred,'.tiff'), fig)
  ggsave(paste0('./reports/figures/',current_inbred,'.png'), fig)
}
```