Props_All.Rmd

---
title: "R Notebook"
output:
  html_document:
    df_print: paged
---


```{r setup_Props_All, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(MuMIn)
```


# Analysis of membrane properties from all animals

## Comparison of stellate cells with Wfs1+ cells

This analysis is currently elsewhere.


## Plot of stellate cell and pyramidal cell properties as a function of dorsoventral location

Gather variables into a single tidy data frame.
```{r}
plot_sc_wfs <- data.sc %>%
  select(vm:fi, dvlocmm) %>%
  mutate(classification = "SC")
plot_wfs <- data.wfs %>%
  select(vm:fi, dvlocmm, classification)
plot_sc_wfs <- bind_rows(plot_sc_wfs, plot_wfs) %>%
  gather("property", "value", vm:fi)

ggplot(plot_sc_wfs, aes(dvlocmm,value)) +
  geom_point(aes(colour = classification)) +
  facet_wrap(~property, scales = "free_y")

```

## Evaluate whether properties of pyramidal cells depend on dv location.

Nest pyramidal cells from wfs data
```{r}
data.wfs_r <- filter(data.wfs, classification == "Pyr") %>%
  select(-classification) %>%
  gather("property", "value", vm:fi) %>%
  group_by(property) %>%
  nest()
```


Fit linear model to each parameter as a function of location
```{r, warning=FALSE}
data.wfs_r <- data.wfs_r %>%
  mutate(mm_vsris = map(data, model_vsris)) %>%
  mutate(mm_vsris_null = map(data, model_vsris_null)) %>%
  mutate(linearmodel = map(data, linearmodel_to_fit))

data.wfs_r <- mixedmod_extract(data.wfs_r, mm_vsris)
data.wfs_r <- mixed_vs_linear_pchisqu(data.wfs_r, mm_vsris, linearmodel) 
data.wfs_r$mm_vsris_vslinear_pdiff_adj <- p.adjust(data.wfs_r$mm_vsris_vslinear_pdiff, method = "BH")

# Look at linear model output
select(data.wfs_r, property, linearmodel) %>%
  mutate(glance = map(linearmodel, broom::glance),
         adj.r2 = map(glance, ~.$adj.r.squared)) %>%
  select(property, adj.r2) %>%
  kableExtra::kable()

# Look at comparison between the linear model and trhe mixed model.
table_mixedvslinear(data.wfs_r, "mm_vsris_vslinear", "property")
```

Compare mixed model containing dvloc with null model
```{r, warning=FALSE}
data.wfs_r <- data.wfs_r %>%
  mutate(anova = map2(mm_vsris, mm_vsris_null, ~anova(.x,.y))) %>%
  mutate(tidy_anova = map(anova, broom::tidy)) %>%
  mutate(anova_p_val = map_dbl(tidy_anova, ~.$p.value[2]))

data.wfs_r <- data.wfs_r %>%
  mutate(anova_p_val_adj = p.adjust(data.wfs_r$anova_p_val, method = "BH"))
```



## Make summary table for the manuscript
```{r}
props_man_table <- c("property", "mm_vsris_gradient_slopes", "anova_p_val_adj", "mm_vsris_marginal.r2", "mm_vsris_conditional.r2", "mm_vsris_slope_min", "mm_vsris_slope_max", "mm_vsris_vslinear_pdiff_adj")
man_table <- unnest(data.wfs_r[props_man_table])

man_table$anova_p_val_adj <- format(man_table$anova_p_val_adj, digits = 3)
man_table$mm_vsris_vslinear_pdiff_adj <- format(man_table$mm_vsris_vslinear_pdiff_adj, digits = 3)
man_table$property <- c("Vm", "IR", "Sag", "Tm", "Res. frequency", "Res. magnitude", "Spike thresold", "Spike maximum", "Spike width", "Rheobase", "Spike AHP", "I-F slope")

(man_table <- man_table %>%
    knitr::kable(
      digits = 4,
      col.names = c("Property", "Slope", "p (slope)", "Marginal R2", "Conditional R2", "Slope (min)", "Slope (max)", "p (vs linear)")
    ) %>%
    kableExtra::kable_styling(bootstrap_options = c("striped", "hover")))

table_path <- paste0(getwd(), "/Tables/Wfs1_mm_props_vs_dv")
man_table %>% save_kable_jpg_html(table_path)

```



## Plot raw and transformed stellate cell properties as a function of location

Plot experimentally measured properties as a function of location.
```{r}
plot_all_raw_data <- data.sc %>%
  select(vm:fi, dvlocmm) %>%
  gather("property", "value", vm:fi)

ggplot(plot_all_raw_data, aes(dvlocmm,value)) +
  geom_point() +
  facet_wrap(~property, scales = "free_y")
```


Plot transformed properties as a function of location.
```{r}
plot_all_tprops_data <- data.sc_TP %>%
  select(vm:fi, dvlocmm) %>%
  gather("property", "value", vm:fi)

ggplot(plot_all_tprops_data, aes(dvlocmm,value)) +
  geom_point() +
  facet_wrap(~property, scales = "free_y")
```

Plot transformed properties as a function of transformed location.
```{r}
plot_all_tall_data <- data.sc_TP %>%
  select(vm:fi, dvlocmm1) %>%
  gather("property", "value", vm:fi)

ggplot(plot_all_tall_data, aes(dvlocmm1,value)) +
  geom_point() +
  facet_wrap(~property, scales = "free_y")
```


## Plot Wfs1 pyramidal cell data as a function of location
```{r}
plot_all_raw_data_wfs <- data.wfs %>%
  filter(classification == "Pyr") %>%
  select(vm:fi, dvlocmm, id) %>%
  gather("property", "value", vm:fi)

ggplot(plot_all_raw_data_wfs, aes(dvlocmm,value)) +
  geom_point(aes(colour = id)) +
  facet_wrap(~property, scales = "free_y", labeller = labeller(property = c(ahp = "AHP min. (mV)", fi = "F-I (Hz / pA)", ir = "IR (MΩ)", resf = "Res F (Hz)", resmag = "Res. mag.", rheo = "Rheobase (pA)", sag = "Sag", spkhlf = "Spike h-w (ms)", spkmax = "Spike max. (mV)", spkthr = "Spike thres. (mV)", tau = "Tm (ms)", vm = "Vrest (mV)"))) +
  theme_classic() +
  theme(legend.position = "bottom")

```


## Plot of SC properties as a function of mediolateral position

The problem with this plot is that apparent differences could be due to differential sampling of dorsal and ventral locations.
```{r}
props_labels <- c(
  ahp = "AHP (mV)",
  fi = "F-I (Hz / pA)",
  ir = "IR (MΩ)",
  resf = "Res F (Hz)",
  resmag = "Res mag",
  rheo = "Rheobase (pA)",
  sag = "Sag",
  spkhlf = "Spike HW (ms)",
  spkmax = "Spike max (mV)",
  spkthr = "Spike threshold (pA)",
  tau = "Tm (ms)",
  vm = "Vm (mV)"
)

data.sc_ml_plot <- filter(data.sc, mlpos != "NA") %>%
  gather("property", "value", vm:fi)

data.sc_ml_plot <- transform(data.sc_ml_plot, property = factor(property, levels = c("vm", "ir", "sag", "tau", "resf", "resmag", "spkthr", "spkmax", "spkhlf", "rheo", "ahp", "fi")))

(ml_boxplots <- ggplot(data.sc_ml_plot, aes(mlpos, value)) +
  geom_boxplot() +
    scale_y_continuous(name="") +
    scale_x_discrete(name="") +
  facet_wrap(vars(property), scales = "free", labeller = labeller(property = props_labels)) +
    theme(strip.text.x = element_text(size = 8),
          axis.text = element_text(size = 9)))
```

Plot properties as a function of dorsoventral location but colour code for ML position.
```{r}
(ml_dv_plot <- ggplot(data.sc_ml_plot, aes(dvloc, value)) +
  geom_point(aes(colour = mlpos)) +
  scale_x_continuous(name="DV location (µm)", limits=c(0, 2500), breaks = c(0, 1000, 2000)) +
  scale_y_continuous(name="") +
  labs(colour = "") +
  facet_wrap(vars(property), scales = "free", labeller = labeller(property = props_labels)) +
  theme(text = element_text(size=9),
        strip.text.x = element_text(size = 8),
        axis.text = element_text(size = 9)))
```

Save the figure
```{r}
ggsave("Figures/mlpos.png", height = 150, width = 180, units = "mm")
```


```{r}
remove(list = c("plot_all_raw_data", "plot_all_raw_data_wfs", "plot_wfs", "plot_sc_wfs"))
```