In mikeod38/dauergut: Scripts for rict-1 dauer paper

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
pathname = here::here()
library(magrittr)
library(knitr)
library(dplyr)
library(ggplot2)
library(ggbeeswarm)
library(beeswarm)
library(ggthemes)
library(viridis)
library(lme4)
library(schoRsch)
library(DHARMa)
library(multcomp)
library(tidyr)
library(rstanarm)
library(ggrepel)
library(rstan)
library(rstanarm)
library(MASS)
devtools::install_github("mikeod38/dauergut")
devtools::install_github("rasmusab/bayesian_first_aid")
library(dauergut)
library(plotly)


dauer_ANOVA <- . %>% lm(data = ., formula = pct ~ genotype)

daf-28, but not daf-7, is strongly reduced by temperature{.tabset}

3A

strains = "N2"
#strains<-c("N2","ft7")
days<-c("11_23_16", "12_6_16")
#days <- "11_23_16" # 12_6 was likely L2s based on the timing

d7GFP<-read.csv(file.path(pathname, "extdata","2A_3A_daf7GFP.csv")) %>%
  filter(mean!=4095 & genotype %in% strains & date %in% days & pheromone == 0 & food == "OP50") %>%
  mutate(genotype = factor(genotype, levels=strains), ID = as.character(ID)) %>%
  separate(ID, c("ID.A", "ID.B"), sep = ":", extra = "drop") %>% 
  mutate(genoID = as.factor(paste(date, genotype, ID.A, sep = ":")), cell.norm = mean)

df <- d7GFP %>% group_by(genotype, date, neuron, temp, genoID) %>% summarise(cell.norm = mean(cell.norm))

#plot
p<-df %>% ggplot(aes(x=temp, y=cell.norm)) +
  list(add.scatter(),add.median(width = .5),add.quartiles()) + 
  labs(
  title = "daf-7 expression does not decrease at high temperatures",
  x = "temperature (ºC)",
  y = "mean intensity") +
  stat_smooth(method="lm",linetype="dashed", size = 0.2, colour = "grey", alpha = 0.05) +
  add.n() +
  theme_classic() + 
  theme(axis.title.x = element_text(size=16),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 15),
        strip.text.x = element_text(size=20))
p

library(lsmeans)
log.tran<-make.tran(type="genlog", param = c(0.5,10))
lm.ASI.lin <- with(log.tran, lm(linkfun(cell.norm) ~ temp, data = df))

Figure 3A daf-7 expression does not strongly decrease at high temperatures. Quantification of daf-7p::gfp fluorescence in the ASI,neurons in L1 larvae grown at the indicated temperatures. Each data point is the average GFP fluorescence of a single animal (2 neurons per animal); numbers in parentheses below indicate the total number of animals examined. Horizontal lines are the median; error bars are quartiles. Dashed lines indicate slope of linear regression as a function of temperature.

library(sjPlot)
sjt.lm(lm.ASI.lin, depvar.labels = "log(GFP expression)", show.fstat = TRUE)

3B

tempdays<-c("1", "4", "6", "8", "10")
foods <- "OP50"
strains <- "N2"
d28<-read.csv(file.path(pathname,'extdata','3BEF_daf-28_GFP_rict-1.csv')) %>%
separate(ID, c("ID.A", "ID.B"), sep = ":", extra = "drop") %>%
  subset(pheromone == 0) %>% 
  mutate(day = as.factor(day), genotype = factor(genotype, levels = strains),
         genoID = interaction(date,neuron,temp,genotype,ID.A))

df <- d28 %>% subset(day %in% tempdays & food %in% foods & mean !=4095 & genotype %in% strains) %>% group_by(date,temp,genotype,neuron,genoID) %>% summarise(cell.norm = mean(mean)) %>%
  data.frame() %>% mutate(group.id = interaction(genotype, neuron)) # take mean of each worm

## plot
(p<-df %>% ggplot(aes(x=temp, y=cell.norm)) +
   list(add.scatter(),add.median(width = .5),add.quartiles()) + 
    theme_classic() +
  facet_grid(.~neuron) +
  scale_color_viridis(discrete = TRUE) +
  labs(title = "daf-28 expression is controlled by temperature",
    x = "temperature",
    y = "mean fluorescence",
    colour = "Temperature") + add.n() +
   stat_smooth(method="lm",linetype="dashed", size = 0.2, colour = "grey", alpha = 0.05) +
   theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size= 20),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 15),
        strip.text.x = element_text(size=20),
        panel.spacing.x=unit(2, "lines")))


lm.ASI.lin <- with(log.tran, lm(linkfun(cell.norm) ~ temp, data = subset(df, neuron == "ASI")))
#lm.ASI.poly.2 <- update(lm.ASI.lin, formula = . ~ . - temp + poly(temp,2,raw=TRUE))
#anova(lm.ASI.lin,lm.ASI.poly.2) #quadratic fits better but not enough temp points to model
#anova(lm.ASI.lin)

lm.ASJ.lin <- with(log.tran, lm(linkfun(cell.norm) ~ temp, data = subset(df, neuron == "ASJ")))
#anova(lm.ASJ.lin)

Figure 3B
daf-28 expression strongly decreases at high temperatures. Quantification of daf-28p::gfp fluorescence in the ASI (left) and ASJ (right) neurons in L1 larvae grown at the indicated temperatures. Each data point is the average GFP fluorescence of a single animal (2 neurons per animal); numbers in parentheses below indicate the total number of animals examined. Horizontal lines are the median; error bars are quartiles. Dashed lines indicate slope of linear regression as a function of temperature.

library(sjPlot)
sjt.lm(lm.ASI.lin, lm.ASJ.lin, depvar.labels = c("ASI", "ASJ"), show.fstat = TRUE)

3C

FISH <- read.csv(file.path(pathname, 'extdata', '3C_daf7_28_FISH.csv')) %>%
  separate(Label, c("group","number"), sep = "-") %>%
  mutate(number = case_when(
    is.na(number) ~ "01",
    TRUE ~ number
  )) %>% 
  separate(group, c("measure", "genotype", "probe", "temp", "tube"), sep = "_") %>%
  mutate(animal = interaction(genotype,probe,temp,tube,number))

background <- FISH %>% filter(type == "background") %>% 
  select(animal,Mean) %>%
  rename(Mean.background = Mean)

df <- filter(FISH, type == "cell") %>% 
  left_join(.,background) %>% mutate(cell.norm = Mean - Mean.background)

daf7 <- filter(df, probe == "daf7Cy5")

#simple lm
linmod<-lm(data=daf7, formula=cell.norm~temp)
#check for outliers
daf7 <- dauergut::flag_outliers(linmod, daf7, threshold = 4, noplot = TRUE) %>% filter(outlier.status == FALSE)


contrast <- daf7 %$% t.test(cell.norm ~ temp)
contrast <- data.frame(p.value = contrast$p.value) %>% dauergut::prange()
plot.contrasts <- c("",contrast$prange)

#linmod<-lm(data=daf7[daf7$outlier.status == FALSE,], formula=cell.norm~temp)
#Fstat not sig#
fit <- BayesianFirstAid::bayes.t.test(daf7[daf7$temp == "25", "cell.norm"], 
                               daf7[daf7$temp == "27", "cell.norm"])
fit_sum <- fit$stats %>% 
  data.frame() %>% 
  select(mean,q2.5.,q97.5.,q25.,q75.) %>%
  tibble::rownames_to_column('parameter') %>%
  dplyr::filter(., grepl("mu",parameter)) %>%
  dplyr::filter(., !grepl("diff",parameter)) %>%
  rename(lower.CL = q2.5., 
         upper.CL = q97.5.,
         lower.25 = q25.,
         upper.75 = q75.) %>%
  mutate(temp = as.factor(c(25,27)),
         x.pos = c(1.3,2.3))

p <- ggplot(daf7, aes(x=temp)) + #x-layer
    theme_my +
    geom_errorbar(data=fit_sum, aes(x=x.pos,y=mean, ymin=lower.CL, ymax=upper.CL),
                  width=0,colour ="grey", lwd=0.15) + #90% cred int for 95% one-sided H0
    geom_errorbar(data=fit_sum, aes(x=x.pos,y=mean, ymin = lower.25, ymax = upper.75),
                  width=0,colour = "darkgrey", lwd = 0.15+0.7) + #75% cred interval
    geom_segment(data = fit_sum, aes(x = x.pos-0.005*nrow(fit_sum),
                                   y = mean, xend = x.pos+0.005*nrow(fit_sum), 
                                   yend = mean), colour = "darkgrey") + 
    scale_x_discrete(labels=function(x) sub(" ","\n",x,fixed=TRUE)) +
    stat_summary(aes(x=temp, y=800), geom="text", label=plot.contrasts, show.legend = TRUE, size=4) + # pvalues
    theme(axis.text.x = element_text(size = 12),
          axis.text.y = element_text(size = 12),
          axis.line = element_line(size=0.2),
          axis.title = element_text(size=16)) +
  geom_quasirandom(aes(y=cell.norm),colour = "#990000", cex=1,
                           width = 0.02,size=0.3,
                           method = 'smiley') +
          stat_summary(aes(y=cell.norm),fun.y = median, 
                       fun.ymin = function(z) {quantile(z,0.25)}, 
                       fun.ymax = function(z) {quantile(z,0.75)},
                       geom = "errorbar", width = 0.15, lwd = 0.15) +
          stat_summary(aes(y=cell.norm),fun.y = median, 
                       fun.ymin = median, 
                       fun.ymax = median,
                       geom = "crossbar", width = 0.25, lwd = 0.35) +
          labs(y = "normalized expression",x = "temperature")

p

Figure 3C Quantification of daf-7 (3C) or daf-28 (3D) mRNA levels in ASI assessed via fluorescent in situ hybridization in animals grown to L1 molt stage at the indicated temperatures. Expression was normalized by subtracting mean background pixel values for each image. Each dot is the mean fluorescence intensity in a single ASI neuron (see Materials and Methods). Horizontal thick bar indicates median. Error bars are quartiles. Light gray thin and thick vertical bars at right indicate Bayesian 95% and 75% credible intervals, respectively. *** P<0.001 relative to animals grown at 25ºC (Welch’s t-test). Representative image of daf-28 mRNA FISH in ASI neurons of animals grown at the indicated temperatures are shown on right in D. Scale bar: 5 μm.

3D

daf28 <- filter(df, probe == "daf28Cy5")

#simple lm
linmod<-lm(data=daf28, formula=cell.norm~temp)
#check for outliers
daf28 <- dauergut::flag_outliers(linmod, daf28, threshold = 4, noplot = TRUE) %>%
  filter(outlier.status == FALSE)

contrast <- daf28 %$% t.test(cell.norm ~ temp)
contrast <- data.frame(p.value = contrast$p.value) %>% dauergut::prange()
plot.contrasts <- c("",contrast$prange)

#linmod<-lm(data=daf28[daf28$outlier.status == FALSE,], formula=cell.norm~temp)
#Fstat not sig#
fit <- BayesianFirstAid::bayes.t.test(daf28[daf28$temp == "25", "cell.norm"], 
                               daf28[daf28$temp == "27", "cell.norm"])
fit_sum <- fit$stats %>% 
  data.frame() %>% 
  select(mean,q2.5.,q97.5.,q25.,q75.) %>%
  tibble::rownames_to_column('parameter') %>%
  dplyr::filter(., grepl("mu",parameter)) %>%
  dplyr::filter(., !grepl("diff",parameter)) %>%
  rename(lower.CL = q2.5., 
         upper.CL = q97.5.,
         lower.25 = q25.,
         upper.75 = q75.) %>%
  mutate(temp = as.factor(c(25,27)),
         x.pos = c(1.3,2.3))

p <- ggplot(daf28, aes(x = temp)) + # x-layer
  theme_my +
  geom_errorbar(
    data = fit_sum, aes(x = x.pos, y = mean, ymin = lower.CL, ymax = upper.CL),
    width = 0, colour = "grey", lwd = 0.15
  ) + # 90% cred int for 95% one-sided H0
  geom_errorbar(
    data = fit_sum, aes(x = x.pos, y = mean, ymin = lower.25, ymax = upper.75),
    width = 0, colour = "darkgrey", lwd = 0.15 + 0.7
  ) + # 75% cred interval
  geom_segment(data = fit_sum, aes(
    x = x.pos - 0.005 * nrow(fit_sum),
    y = mean, xend = x.pos + 0.005 * nrow(fit_sum),
    yend = mean
  ), colour = "darkgrey") +
  scale_x_discrete(labels = function(x) sub(" ", "\n", x, fixed = TRUE)) +
  coord_cartesian(ylim = c(0,320)) +
  stat_summary(aes(x = temp, y = 300), geom = "text", label = plot.contrasts, show.legend = TRUE, size = 4) + # pvalues
  theme(
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    axis.line = element_line(size = 0.2),
    axis.title = element_text(size = 16)
  ) +
  geom_quasirandom(
    aes(y = cell.norm), colour = "#990000", cex = 1,
    width = 0.02, size = 0.3,
    method = "smiley"
  ) +
  stat_summary(
    aes(y = cell.norm), fun.y = median,
    fun.ymin = function(z) {
      quantile(z, 0.25)
    },
    fun.ymax = function(z) {
      quantile(z, 0.75)
    },
    geom = "errorbar", width = 0.15, lwd = 0.15
  ) +
  stat_summary(
    aes(y = cell.norm), fun.y = median,
    fun.ymin = median,
    fun.ymax = median,
    geom = "crossbar", width = 0.25, lwd = 0.35
  ) +
  labs(y = "normalized expression", x = "temperature")

p

Figure 3D Quantification of daf-7 (3C) or daf-28 (3D) mRNA levels in ASI assessed via fluorescent in situ hybridization in animals grown to L1 molt stage at the indicated temperatures. Expression was normalized by subtracting mean background pixel values for each image. Each dot is the mean fluorescence intensity in a single ASI neuron (see Materials and Methods). Horizontal thick bar indicates median. Error bars are quartiles. Light gray thin and thick vertical bars at right indicate Bayesian 95% and 75% credible intervals, respectively. *** P<0.001 relative to animals grown at 25ºC (Welch’s t-test). Representative image of daf-28 mRNA FISH in ASI neurons of animals grown at the indicated temperatures are shown on right in D. Scale bar: 5 μm.

3E-F

days = as.factor(8:10) # days containing N2 and ft7 27º data
temps = c("25","27")
strains = c("N2", "ft7")
foods = "OP50"
d28<-read.csv(file.path(pathname,'extdata', '3BEF_daf-28_GFP_rict-1.csv')) %>%
separate(ID, c("ID.A", "ID.B"), sep = ":", extra = "drop") %>%
  subset(mean!=4095 & food == foods &
           temp %in% temps &
           genotype %in% strains &
           pheromone == 0 &
           day %in% days) %>%
  mutate(day = as.factor(day), genotype = factor(genotype, levels = strains),
         genoID = interaction(date,genotype,ID.A), temp = factor(temp, levels = temps))

df <- d28 %>% group_by(date,temp,genotype,neuron,genoID) %>% summarise(cell.norm = mean(mean)) %>%
  data.frame() %>% mutate(group.id = interaction(genotype, temp)) %>% arrange(temp, genotype, neuron) # take mean 

# anova(with(log.tran,lm(linkfun(cell.norm) ~ temp, data = df[df$neuron == "ASI",])),with(log.tran,lm(linkfun(cell.norm) ~ temp+genotype, data = df[df$neuron == "ASI",])))
# # for ASI, F = 4.0071, 1Df p = 0.0481 for genotype
# anova(with(log.tran,lm(linkfun(cell.norm) ~ temp, data = df[df$neuron == "ASJ",])),with(log.tran,lm(linkfun(cell.norm) ~ temp+genotype, data = df[df$neuron == "ASJ",])))
# # for ASJ, F = 38.082, 1Df, p < 0.001 for genotype
# # no significant interaction of temp with genotype for either neuron

lm.ASI <- with(log.tran, lm(linkfun(cell.norm) ~ temp + genotype, data = df[df$neuron == "ASI",]))
lm.ASJ <- with(log.tran, lm(linkfun(cell.norm) ~ temp + genotype, data = df[df$neuron == "ASJ",]))

contrasts.ASI<-summary(lsmeans(lm.ASI, pairwise ~ genotype | temp)$contrasts) %>% dauergut::prange()
contrasts.ASJ<-summary(lsmeans(lm.ASJ, pairwise ~ genotype | temp)$contrasts) %>% dauergut::prange()

stan.ASJ<-with(log.tran, stan_lmer(linkfun(cell.norm) ~ group.id + (1|genoID), data = df[df$neuron == "ASJ",], adapt_delta = 0.99, iter = 4000))

stan.ASI<-with(log.tran, stan_lmer(linkfun(cell.norm) ~ group.id + (1|genoID), data = df[df$neuron == "ASI",], adapt_delta = 0.99, iter = 4000))

mixed.ASI <- stan.ASI %>% getStan_CIs(type = "log", group = "temp", base=10, compare_gt = TRUE)
mixed.ASJ <- stan.ASJ %>% getStan_CIs(type = "log", group = "temp", base=10, compare_gt = TRUE)

labels <- c("WT", "rict-1(ft7)")


#### plot ASI ###
mixed <- mixed.ASI
p1 <- df %>% dplyr::filter(neuron == "ASI") %>% ggplot(aes(x = genotype, y = cell.norm)) + 
  geom_quasirandom(aes(y=cell.norm, pch = temp),colour = "#339900", cex=1,
                           width = 0.075,size=0.3,
                           method = 'smiley') + 
  list(add.median(width=0.25),add.quartiles()) + facet_grid(.~temp) +
  stat_summary(aes(x=as.numeric(as.factor(genotype)) + 0.3, y=0),
                   fun.data = fun_length, geom = "text", size = 3) +
  add.Bayes.CI() +
  coord_cartesian(ylim=c(0,4095)) +
  labs(title = "ASI",
    x = "genotype",
    y = "mean fluorescence (AU)",
    colour = "Temperature") +
   theme_classic() + 
  theme(axis.text.x = element_text(size = 12, face = "italic"),
          axis.text.y = element_text(size = 12),
          axis.line = element_line(size=0.2),
          axis.title = element_text(size=16)) +
  scale_x_discrete(labels = labels) +
  guides(pch=FALSE)

#### plot ASJ ###
mixed<-mixed.ASJ
p2 <- df %>% dplyr::filter(neuron == "ASJ") %>% ggplot(aes(x = genotype, y = cell.norm)) + 
  geom_quasirandom(aes(y=cell.norm, pch = temp),colour = "#339900", cex=1,
                           width = 0.075,size=0.3,
                           method = 'smiley') + 
  list(add.median(width=0.25),add.quartiles()) + facet_grid(.~temp) +
  stat_summary(aes(x=as.numeric(as.factor(genotype)) + 0.3, y=0),
                   fun.data = fun_length, geom = "text", size = 3) +
  add.Bayes.CI() +
  coord_cartesian(ylim=c(0,4095))+
  labs(title = "ASJ",
    x = "genotype",
    y = "mean fluorescence (AU)",
    colour = "Temperature") +
  theme_classic() + 
  theme(axis.text.x = element_text(size = 12, face = "italic"),
          axis.text.y = element_text(size = 12),
          axis.line = element_line(size=0.2),
          axis.title = element_text(size=16)) +
  scale_x_discrete(labels = labels) +
  guides(pch=FALSE)

wzxhzdk:11

wzxhzdk:12

**Figure 3E-F** _daf-28_ expression is reduced in _rict-1_ mutants at multiple temperatures. Quantification of _daf-28p::gfp_ fluorescence in the ASI (left) and ASJ (right) neurons in L1 larvae grown at the indicated temperatures. Each data point is the average GFP fluorescence of a single animal (2 neurons per animal); numbers in parentheses below indicate the total number of animals examined. Horizontal lines are the median; error bars are quartiles.

library(sjPlot)
sjt.lm(lm.ASI, lm.ASJ, depvar.labels = c("ASI", "ASJ"), show.fstat = TRUE)

mikeod38/dauergut documentation built on May 30, 2019, 7:16 p.m.