Chapter 1 Analysis for Figure 2 of “ASK1 inhibits browning of white adipose tissue in obesity”

# wrangling packages
library(here)
library(janitor)
library(readxl)
library(data.table)
library(stringr)

# analysis packages
library(emmeans)
library(car) # qqplot, spreadlevel
library(DHARMa)

# graphing packages
library(ggsci)
library(ggpubr)
library(ggforce)
library(cowplot)
library(lazyWeave) #pvalstring

here <- here::here()
data_path <- "data"

ggplotsci_path <- here::here("R", "ggplotsci.R")
source(ggplotsci_path)

Data source: ASK1 inhibits browning of white adipose tissue in obesity

This chunk assigns the path to the Excel data file for all panels of Figure 2. The data for each panel are in a single sheet in the Excel file named “Source Date_Figure 2”.

data_folder <- "ASK1 inhibits browning of white adipose tissue in obesity"
file_name <- "41467_2020_15483_MOESM4_ESM.xlsx"
file_path <- here::here(data_path, data_folder, file_name)
  
fig_2_sheet <- "Source Date_Figure 2"

1.1 useful functions

A function to import longitudinal data from Fig 2

# function to read in parts of 2b
import_fig_2_part <- function(range_2){
  fig_2_part <- read_excel(file_path,
                          sheet = fig_2_sheet,
                          range = range_2,
                          col_names = TRUE) %>%
    data.table()
  group <- colnames(fig_2_part)[1]
  setnames(fig_2_part, old = group, new = "treatment")
  fig_2_part[, treatment := as.character(treatment)] # this was read as logical
  fig_2_part[, treatment := group] # assign treatment group
  fig_2_part[, mouse_id := paste(group, .I)]
  return(fig_2_part)
}
# script to compute various area under the curves (AUC) using trapezoidal method
# le Floch's "incremental" auc substracts the baseline value from all points.
# This can create some elements with negative area if post-baseline values are less
# than baseline value.
# Some researchers "correct" this by setting any(y - ybar < 0 to zero. Don't do this.

auc <- function(x, y, method="auc", average = FALSE){
  # method = "auc", auc computed using trapezoidal calc
  # method = "iauc" is an incremental AUC of Le Floch
  # method = "pos_iauc" is a "positive" incremental AUC of Le Floch but not Wolever
  # method = "post_0_auc" is AUC of post-time0 values
  # if average then divide area by duration
  if(method=="iauc"){y <- y - y[1]}
  if(method=="pos_iauc"){y[y < 0] <- 0}
  if(method=="post_0_auc"){
    x <- x[-1]
    y <- y[-1]
  }
  n <- length(x)
  area <- 0
  for(i in 2:n){
    area <- area + (x[i] - x[i-1])*(y[i-1] + y[i])
  }
  value <- area/2
  if(average == TRUE){
    value <- value/(x[length(x)] - x[1])
  }
  return(value)
}
pal_nature_mod <- c(
  "#3DB7E9", # summer sky
  "#e69f00", # gamboge, squash, buttercup
  "#359B73", # ocean green
  "#2271B2", # honolulu blue
  "#f0e442", # holiday, 
  "#F748A5", # barbi pink
  "#d55e00" # bamboo
)

1.2 figure 2b – effect of ASK1 KO on growth (body weight)

1.2.1 figure 2b – import

range_list <- c("A21:N41", "A43:N56", "A58:N110", "A112:N170")
fig_2b_wide <- data.table(NULL)
for(range_i in range_list){
  part <- import_fig_2_part(range_i)
  fig_2b_wide <- rbind(fig_2b_wide,
                       part)
}

fig_2b <- melt(fig_2b_wide,
               id.vars = c("treatment", "mouse_id"),
               variable.name = "week",
               value.name = "body_weight")
fig_2b[, week := as.numeric(as.character(week))]
fig_2b[, c("ask1", "diet") := tstrsplit(treatment, " ", fixed=TRUE)]
fig_2b[, week_f := factor(week)]

1.2.2 figure 2b – exploratory plots

qplot(x = week,
      y = body_weight,
      data = fig_2b,
      color = treatment) +
  facet_grid(ask1~diet)

  1. no obvious outliers
  2. reduced growth rate at bigger size
qplot(x = week,
      y = body_weight,
      data = fig_2b,
      color = treatment) +
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

1. loess smooth. Growth in ASK1F/F + HFD clearly greater than other three treatment combinations.

1.3 Figure 2c – Effect of ASK1 KO on final body weight

1.3.1 Figure 2c – import

range_2c <- "A173:BD176"
y_cols <- c("ASK1F/F chow", "ASK1Δadipo chow", "ASK1F/F HFD", "ASK1Δadipo HFD")
fig_2c_import <- read_excel(file_path,
                     sheet = fig_2_sheet,
                     range = range_2c,
                     col_names = FALSE) %>%
  transpose(make.names=1) %>%
  data.table() %>%
  melt(measure.vars = y_cols,
       variable.name = "treatment",
       value.name = "body_weight_gain") %>%
  na.omit()
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * ...
fig_2c_import[, mouse_id := paste(treatment, .I, sep = "_"),
       by = treatment]

1.3.2 Figure 2c – check own computation of weight change v imported value

Note that three cases are missing from fig_2c import that are in fig_2b

# change colnames to char
fig_2c <- copy(fig_2b_wide)
weeks <- unique(fig_2b[, week])
setnames(fig_2c,
         old = colnames(fig_2c),
         new = c("treatment", paste0("week_", weeks), "mouse_id"))
fig_2c[, weight_gain := week_12 - week_0]
fig_2c <- fig_2c[, .SD, .SDcols = c("treatment",
                                "week_0",
                                "week_12",
                                "weight_gain")]
fig_2c[, mouse_id := paste(treatment, .I, sep = "_"),
       by = treatment]

fig_2c[, c("ask1", "diet") := tstrsplit(treatment, " ", fixed=TRUE)]

fig_2c_check <- merge(fig_2c,
                fig_2c_import,
                by = c("mouse_id"),
                all = TRUE)

#View(fig_2c_check)

1.3.3 Figure 2c – exploratory plots

qplot(x = treatment,
      y = weight_gain,
      data = fig_2c)

  • no obvious outliers
  • variance increases with mean, as expected from growth, suggests a multiplicative model. But start with simple lm.

1.3.4 Figure 2c – fit the model: m1 (lm)

fig_2c_m1 <- lm(weight_gain ~ week_0 + ask1*diet, data = fig_2c)

1.3.5 Figure 2c – check the model: m1

# check normality assumption
set.seed(1)
qqPlot(fig_2c_m1, id=FALSE)

spreadLevelPlot(fig_2c_m1, id=FALSE)

## 
## Suggested power transformation:  0.06419448
  • QQ indicates possible right skew but especially left side is squashed toward mean
  • spread-level indicates variance increases with mean

For p-value, this may not be too severe but for intervals, best to account for this. Try gamma with log link (which makes biological sense for growth)

1.3.6 Figure 2c – fit the model: m2 (gamma glm)

fig_2c_m2 <- glm(weight_gain ~ week_0 + ask1*diet,
          family = Gamma(link = "log"),
          data = fig_2c)

1.3.7 Figure 2c – check the model, m2

set.seed(1)
fig_2c_m2_sim  <-  simulateResiduals(fig_2c_m2,  n=250)
## Model family was recognized or set as continuous, but duplicate values were detected in the response. Consider if you are fitting an appropriate model.
plot(fig_2c_m2_sim) 

  • well behaved QQ and spread-level

1.3.8 Figure 2c – inference from the model

coef_table <- cbind(coef(summary(fig_2c_m2)),
                    exp(confint(fig_2c_m2))) %>%
  data.table(keep.rownames = TRUE)
coef_table[, Estimate:=exp(Estimate)]
knitr::kable(coef_table, digits = c(0,2,2,2,4,2,2))
rn Estimate Std. Error t value Pr(>|t|) 2.5 % 97.5 %
(Intercept) 10.87 0.35 6.83 0.0000 5.53 21.50
week_0 0.99 0.01 -0.84 0.3997 0.96 1.02
ask1ASK1Δadipo 1.03 0.11 0.24 0.8138 0.83 1.28
dietHFD 1.56 0.08 5.45 0.0000 1.33 1.82
ask1ASK1Δadipo:dietHFD 0.79 0.13 -1.93 0.0551 0.61 1.00
fig_2c_m2_emm <- emmeans(fig_2c_m2, specs = c("diet", "ask1"),
                         type = "response")
fig_2c_m2_pairs <- contrast(fig_2c_m2_emm,
                     method = "revpairwise",
                     simple = "each",
                     combine = TRUE,
                     adjust = "none") %>%
  summary(infer = TRUE)
fig_2c_m2_emm
##  diet ask1       response    SE  df asymp.LCL asymp.UCL
##  chow ASK1F/F        8.19 0.568 Inf      7.15      9.39
##  HFD  ASK1F/F       12.75 0.548 Inf     11.72     13.87
##  chow ASK1Δadipo     8.41 0.721 Inf      7.11      9.95
##  HFD  ASK1Δadipo    10.28 0.424 Inf      9.48     11.14
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale
fig_2c_m2_pairs
##  ask1       diet contrast             ratio     SE  df asymp.LCL asymp.UCL
##  ASK1F/F    .    HFD / chow           1.556 0.1263 Inf     1.328     1.825
##  ASK1Δadipo .    HFD / chow           1.222 0.1168 Inf     1.013     1.474
##  .          chow ASK1Δadipo / ASK1F/F 1.026 0.1127 Inf     0.828     1.273
##  .          HFD  ASK1Δadipo / ASK1F/F 0.806 0.0485 Inf     0.716     0.907
##  z.ratio p.value
##   5.453  <.0001 
##   2.097  0.0360 
##   0.236  0.8134 
##  -3.587  0.0003 
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## Tests are performed on the log scale
  • within ASK1 Cn, HFD mean is 1.6X chow mean
  • within ASK1 KO, HFD mean is 1.2X chow mean
  • within chow, ASK1 KO mean is 1.0X ASK1 Cn mean
  • within HFD, ASK1 KO mean is 0.86X ASK1 Cn mean
# same as interaction effect in coefficient table
contrast(fig_2c_m2_emm, interaction = "pairwise", by = NULL) %>%
  summary(infer = TRUE)
##  diet_pairwise ask1_pairwise        ratio     SE  df asymp.LCL asymp.UCL
##  chow / HFD    ASK1F/F / ASK1Δadipo 0.785 0.0982 Inf     0.614         1
##  z.ratio p.value
##  -1.934  0.0531 
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## Tests are performed on the log scale
  • the reduction in weight gain in the ASK1 KO mice compared to ASK1 CN is 0.785X. Notice that p > 0.05.

1.3.9 Figure 2c – plot the model

fig_2c_m2_emm_dt <- summary(fig_2c_m2_emm) %>%
  data.table
fig_2c_m2_pairs_dt <- data.table(fig_2c_m2_pairs)
fig_2c_m2_pairs_dt[ , p_pretty := pvalString(p.value)]

dodge_width <- 0.8 # separation between groups
# get x positions of brackets for p-values
# requires looking at table and mentally figuring out
# Chow is at x = 1 and HFD is at x = 2
fig_2c_m2_pairs_dt[, group1 := c(1-dodge_width/4,
                                 1+dodge_width/4,
                                 1-dodge_width/4,
                                 2-dodge_width/4)]
fig_2c_m2_pairs_dt[, group2 := c(2-dodge_width/4,
                                 2+dodge_width/4,
                                 1+dodge_width/4,
                                 2+dodge_width/4)]

pd <- position_dodge(width = dodge_width)
fig_2c_gg <- ggplot(data = fig_2c,
                    aes(x = diet,
                        y = weight_gain,
                        color = ask1)) +
  
  # points
  geom_sina(alpha = 0.5,
            position = pd) +
  
  # plot means and CI
  geom_errorbar(data = fig_2c_m2_emm_dt,
                aes(y = response,
                    ymin = asymp.LCL,
                    ymax = asymp.UCL,
                    color = ask1),
                width = 0,
                position = pd
  ) +
  
  geom_point(data = fig_2c_m2_emm_dt,
             aes(y = response,
                 color = ask1),
             size = 3,
             position = pd
  ) +
  
  # plot p-values (y positions are adjusted by eye)
  stat_pvalue_manual(fig_2c_m2_pairs_dt,
                     label = "p_pretty",
                     y.position=c(28.5, 31, 26, 26),
                     tip.length = 0.01) +
  
  # aesthetics
  ylab("Weight Gain") +
  scale_color_manual(values=pal_nature_mod,
                     name = NULL) +
  theme_pubr() +
  theme(legend.position="top") +
  theme(axis.title.x=element_blank()) +
  
  NULL

fig_2c_gg

1.3.10 Figure 2c – report

Results could be reported using either:

(This is inconsistent with plot, if using this, the plot should reverse what factor is on the x-axis and what factor is the grouping (color) variable) Mean weight gain in ASK1F/F mice on HFD was 1.56 (95% CI: 1.33, 1.82, \(p < 0.0001\)) times that of ASK1F/F mice on chow while mean weight gain in ASK1Δadipo mice on HFD was only 1.22 (95% CI: 1.01, 1.47, \(p = 0.036\)) times that of ASK1Δadipo mice on chow. This reduction in weight gain in ASK1Δadipo mice compared to ASK1F/F control mice was 0.79 times (95% CI; 0.61, 1.00, \(p = 0.0531\)).

(This is consistent with the plot in that its comparing difference in the grouping factor within each level of the factor on the x-axis) Mean weight gain in ASK1Δadipo mice on chow was trivially larger (1.03 times) than that in ASK1F/F mice on chow (95% CI: 0.83, 1.27, \(p = 0.81\)) while mean weight gain in ASK1Δadipo mice on HFD was smaller (0.81 times) than that in ASK1F/F control mice on HFD (95% CI: 0.72 , 0.91, \(p = 0.0003\)). This reduction in weight gain in ASK1Δadipo mice compared to ASK1Δadipo mice is 0.79 times (95% CI; 0.61, 1.00, \(p = 0.0531\)).

note to research team. The big difference in p-values between weight difference on chow and weight difference on HFD might lead one to believe there is a “difference in this difference”. Using a p-value = effect strategy, this is not supported.

1.4 Figure 2d – Effect of ASK1 KO on glucose tolerance (whole curve)

1.4.1 Figure 2d – Import

range_list <- c("A179:H189", "A191:H199", "A201:H214", "A216:H230")
fig_2d_wide <- data.table(NULL)
for(range_i in range_list){
  part <- import_fig_2_part(range_i)
  fig_2d_wide <- rbind(fig_2d_wide,
                       part)
}
fig_2d_wide[, c("ask1", "diet") := tstrsplit(treatment, " ", fixed=TRUE)]

# melt
fig_2d <- melt(fig_2d_wide,
               id.vars = c("treatment",
                           "ask1",
                           "diet",
                           "mouse_id"),
               variable.name = "time",
               value.name = "glucose")
fig_2d[, time := as.numeric(as.character(time))]

# for plot only (not analysis!)
shift <- 2
fig_2d[treatment == "ASK1F/F chow", time_x := time - shift*1.5]
fig_2d[treatment == "ASK1Δadipo chow", time_x := time - shift*.5]
fig_2d[treatment == "ASK1F/F HFD", time_x := time + shift*.5]
fig_2d[treatment == "ASK1Δadipo HFD", time_x := time + shift*1.5]

1.4.2 Figure 2d – exploratory plots

qplot(x = time_x, y = glucose, color = treatment, data = fig_2d) +
  geom_line(aes(group = mouse_id), alpha = 0.3)

* no obvious unplausible outliers but two mice w/ high values in “F/F HFD” * similar at time zero (initial effect is trivial)

use AUC conditional on time 0 glucose

qplot(x = time,
      y = glucose,
      data = fig_2d,
      color = treatment) +
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

1.4.3 Figure 2d – fit the model

1.4.4 Figure 2d – check the model

1.4.5 Figure 2d – inference

1.4.6 Figure 2d – plot the model

fig_2d_gg <- ggplot(data = fig_2d,
                    aes(x = time_x,
                        y = glucose,
                        color = ask1,
                        shape = diet)) +
  geom_point() +
  geom_line(aes(group = mouse_id), alpha = 0.3) +
  
  # aesthetics
  ylab("Blood glucose (mmol/L)") +
  xlab("Time (min)") +
  scale_x_continuous(breaks=c(0,15, 30, 45, 60, 90, 120)) +
  scale_color_manual(values=pal_nature_mod,
                     name = NULL) +
  scale_shape_manual(values=c(16, 1),
                     name = NULL) +
  
  theme_pubr() +
  theme(legend.position="top") +
  #theme(axis.title.x=element_blank()) +
  
  NULL

fig_2d_gg

1.5 Figure 2e – Effect of ASK1 KO on glucose tolerance (summary measure)

The researchers did create a table to import but this analysis uses the mean post-baseline glucose amount as the response instead of the area under the curve of over the full 120 minutes. This mean is computed as the post-baseline area under the curve divided by the duration of time of the post-baseline measures (105 minutes). This analysis will use fig_2d_wide since there is only one a single Y variable per mouse.

1.5.1 Figure 2e – message the data

# AUC of post-baseline values
# do this after melt as we don't need this in long format)
fig_2d_wide[, glucose_0 := get("0")]

times <- c(0, 15, 30, 45, 60, 90, 120)
time_cols <- as.character(times)
Y <- fig_2d_wide[, .SD, .SDcols = time_cols]
fig_2d_wide[, glucose_mean := apply(Y, 1, auc,
                           x=times,
                           method = "post_0_auc",
                           average = TRUE)]

1.5.2 Figure 2e – exploratory plots

qplot(x = treatment, y = glucose_mean, data = fig_2d_wide)

1.5.3 Figure 2e – fit the model

fig_2e_m1 <- lm(glucose_mean ~ glucose_0 + ask1*diet, data = fig_2d_wide)

1.5.4 Figure 2e – check the model

# check normality assumption
set.seed(1)
qqPlot(fig_2e_m1, id=FALSE)

spreadLevelPlot(fig_2e_m1, id=FALSE)

## 
## Suggested power transformation:  -0.4035073

1.5.5 Figure 2e – inference from the model

fig_2e_m1_coef <- coef(summary(fig_2e_m1))
fig_2e_m1_coef
##                          Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)             8.6835026  1.2556730  6.9154169 2.460045e-08
## glucose_0               1.2194720  0.2390919  5.1004329 8.592596e-06
## ask1ASK1Δadipo         -0.3488511  0.8759124 -0.3982716 6.925479e-01
## dietHFD                 4.2782121  0.7908612  5.4095613 3.185930e-06
## ask1ASK1Δadipo:dietHFD -2.7503448  1.1288320 -2.4364518 1.937783e-02
fig_2e_m1_emm <- emmeans(fig_2e_m1, specs = c("diet", "ask1"))
fig_2e_m1_pairs <- contrast(fig_2e_m1_emm,
                     method = "revpairwise",
                     simple = "each",
                     combine = TRUE,
                     adjust = "none") %>%
  summary(infer = TRUE)
fig_2e_m1_emm
##  diet ask1       emmean    SE df lower.CL upper.CL
##  chow ASK1F/F      14.7 0.587 40     13.5     15.9
##  HFD  ASK1F/F      18.9 0.520 40     17.9     20.0
##  chow ASK1Δadipo   14.3 0.659 40     13.0     15.7
##  HFD  ASK1Δadipo   15.9 0.493 40     14.9     16.8
## 
## Confidence level used: 0.95
fig_2e_m1_pairs
##  ask1       diet contrast             estimate    SE df lower.CL upper.CL
##  ASK1F/F    .    HFD - chow              4.278 0.791 40    2.680     5.88
##  ASK1Δadipo .    HFD - chow              1.528 0.824 40   -0.138     3.19
##  .          chow ASK1Δadipo - ASK1F/F   -0.349 0.876 40   -2.119     1.42
##  .          HFD  ASK1Δadipo - ASK1F/F   -3.099 0.715 40   -4.544    -1.65
##  t.ratio p.value
##   5.410  <.0001 
##   1.853  0.0712 
##  -0.398  0.6925 
##  -4.335  0.0001 
## 
## Confidence level used: 0.95

1.5.6 Figure 2e – plot the model

fig_2e_gg <- gg_mean_error(data = fig_2d_wide,
                          fit = fig_2e_m1,
                          fit_emm = fig_2e_m1_emm,
                          fit_pairs = fig_2e_m1_pairs,
                          x_col = "diet",
                          y_col = "glucose_mean",
                          g_col = "ask1",
                          wrap_col = NULL,
                          x_label = "none",
                          y_label = "Post-baseline glucose (mmol per l)",
                          g_label = "",
                          dots = "sina",
                          dodge_width = 0.8,
                          adjust = 0.5,
                          p_show = c(3, 4, 2,1),
                          p_pos = c(1,1,2,3))
fig_2e_gg

1.6 Figure 2f – Effect of ASK1 on glucose infusion rate

1.6.1 Figure 2f – import

range_2f <- "A239:I240"
treatment_levels <- c("ASK1F/F", "ASK1Δadipo")
fig_2f <- read_excel(file_path,
                     sheet = fig_2_sheet,
                     range = range_2f,
                     col_names = FALSE) %>%
  transpose(make.names=1) %>%
  data.table() %>%
  melt(measure.vars = treatment_levels,
       variable.name = "treatment",
       value.name = "glucose_infusion_rate") %>%
  na.omit()

fig_2f[, treatment := factor(treatment, treatment_levels)]

1.6.2 Figure 2f – exploratory plots

1.6.3 Figure 2f – fit the model

fig_2f_m1 <- lm(glucose_infusion_rate ~ treatment, data = fig_2f)

1.6.4 Figure 2f – check the model

1.6.5 Figure 2f – inference

fig_2f_m1_coef <- summary(fig_2f_m1) %>%
  coef()
fig_2f_m1_emm <- emmeans(fig_2f_m1, specs = "treatment")
fig_2f_m1_pairs <- contrast(fig_2f_m1_emm,
                            method = "revpairwise") %>%
  summary(infer = TRUE)
fig_2f_m1_pairs
##  contrast             estimate  SE df lower.CL upper.CL t.ratio p.value
##  ASK1Δadipo - ASK1F/F     18.9 6.3 12     5.18     32.6 3.000   0.0111 
## 
## Confidence level used: 0.95

1.6.6 Figure 2f – plot the model

fig_2f_m1_emm_dt <- summary(fig_2f_m1_emm) %>%
  data.table
fig_2f_m1_pairs_dt <- data.table(fig_2f_m1_pairs)
fig_2f_m1_pairs_dt[ , p_pretty := pvalString(p.value)]
fig_2f_m1_pairs_dt[, group1 := 1]
fig_2f_m1_pairs_dt[, group2 := 2]


fig_2f_gg <- ggplot(data = fig_2f,
                    aes(x = treatment,
                        y = glucose_infusion_rate,
                        color = treatment)) +
  
  # points
  geom_sina(alpha = 0.5,
            position = pd) +
  
  # plot means and CI
  geom_errorbar(data = fig_2f_m1_emm_dt,
                aes(y = emmean,
                    ymin = lower.CL,
                    ymax = upper.CL,
                    color = treatment),
                width = 0,
                position = pd,
                color = "black"
  ) +
  
  geom_point(data = fig_2f_m1_emm_dt,
             aes(y = emmean,
                 color = treatment),
             size = 3,
             position = pd,
             color = "black"
  ) +
  
 # plot p-values (y positions are adjusted by eye)
  stat_pvalue_manual(fig_2f_m1_pairs_dt,
                     label = "p_pretty",
                     y.position=c(95),
                     tip.length = 0.01) +
  
  # aesthetics
  ylab("Glucose infusion rate") +
  scale_color_manual(values=pal_nature_mod,
                     name = NULL) +
  theme_pubr() +
  theme(legend.position="none") +
  theme(axis.title.x=element_blank()) +
  
  NULL

fig_2f_gg

1.7 Figure 2g

1.7.1 Figure 2g – import

range_list <- c("A244:G247", "A250:H253")
# import ASK1F/F
fig_2g_1 <- read_excel(file_path,
                     sheet = fig_2_sheet,
                     range = "A244:G247",
                     col_names = FALSE) %>%
  transpose(make.names=1) %>%
  data.table()
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * ...
fig_2g_1[, treatment := "ASK1F/F"]

# import ASK1Δadipo
fig_2g_2 <- read_excel(file_path,
                     sheet = fig_2_sheet,
                     range = "A250:H253",
                     col_names = FALSE) %>%
  transpose(make.names=1) %>%
  data.table()
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * ...
fig_2g_2[, treatment := "ASK1Δadipo"]

# combine
fig_2g <- rbind(fig_2g_1, fig_2g_2)

1.7.2 Figure 2g – exploratory plots

1.7.3 Figure 2g – fit the model

# a more sophisticated would be a mixed model to dampen noise
fig_2g_m1_ingWAT <- lm(ingWAT ~ treatment, data = fig_2g)
fig_2g_m1_epiWAT <- lm(epiWAT ~ treatment, data = fig_2g)
fig_2g_m1_Muscle <- lm(Muscle ~ treatment, data = fig_2g)
fig_2g_m1_BAT <- lm(BAT ~ treatment, data = fig_2g)

1.7.4 Figure 2g – check the model

1.7.5 Figure 2g – inference

fig_2g_infer <- function(m1){
  m1_emm <- emmeans(m1, specs = "treatment")
  m1_pairs <- contrast(m1_emm,
                              method = "revpairwise") %>%
    summary(infer = TRUE)
  return(list(emm = m1_emm,
              pairs = m1_pairs))
}

fig_2g_m1_emm_dt <- data.table(NULL)
fig_2g_m1_pairs_dt <- data.table(NULL)
m1_list <- list(fig_2g_m1_ingWAT, 
             fig_2g_m1_epiWAT,
             fig_2g_m1_Muscle,
             fig_2g_m1_BAT)
y_cols <- c("ingWAT", "epiWAT", "Muscle", "BAT")
for(i in 1:length(y_cols)){
  m1_infer <- fig_2g_infer(m1_list[[i]])

  m1_emm_dt <- summary(m1_infer$emm) %>%
    data.table
  fig_2g_m1_emm_dt <- rbind(fig_2g_m1_emm_dt,
                            data.table(tissue = y_cols[i], m1_emm_dt))
  m1_pairs_dt <- m1_infer$pairs %>%
    data.table
  fig_2g_m1_pairs_dt <- rbind(fig_2g_m1_pairs_dt,
                            data.table(tissue = y_cols[i], m1_pairs_dt))
}
fig_2g_m1_pairs_dt
##    tissue             contrast  estimate        SE df     lower.CL
## 1: ingWAT ASK1Δadipo - ASK1F/F  3.595000  1.468289 10   0.32344725
## 2: epiWAT ASK1Δadipo - ASK1F/F  1.390238  0.669957 11  -0.08432738
## 3: Muscle ASK1Δadipo - ASK1F/F  2.694048  5.675468 11  -9.79757382
## 4:    BAT ASK1Δadipo - ASK1F/F 33.855000 28.715230  7 -34.04572935
##      upper.CL   t.ratio    p.value
## 1:   6.866553 2.4484273 0.03435010
## 2:   2.864804 2.0751153 0.06222096
## 3:  15.185669 0.4746829 0.64429728
## 4: 101.755729 1.1789911 0.27691810

1.7.6 Figure 2g – plot the model

# melt fig_2g

fig_2g_long <- melt(fig_2g,
                    id.vars = "treatment",
                    variable.name = "tissue",
                    value.name = "glucose_uptake")

# change name of ASK1Δadipo label
fig_2g_long[treatment == "ASK1Δadipo", treatment := "ASK1-/-adipo"]
fig_2g_m1_emm_dt[treatment == "ASK1Δadipo", treatment := "ASK1-/-adipo"]

fig_2g_m1_pairs_dt[ , p_pretty := pvalString(p.value)]
fig_2g_m1_pairs_dt[, group1 := 1]
fig_2g_m1_pairs_dt[, group2 := 2]

fig_2g_plot <- function(tissue_i,
                        y_lab = FALSE, # title y-axis?
                        g_lab = FALSE  # add group label?
                        ){
  y_max <- max(fig_2g_long[tissue == tissue_i, glucose_uptake], na.rm=TRUE)
  y_min <- min(fig_2g_long[tissue == tissue_i, glucose_uptake], na.rm=TRUE)
  y_pos <- y_max + (y_max-y_min)*.05
  
  gg <- ggplot(data = fig_2g_long[tissue == tissue_i],
                      aes(x = treatment,
                          y = glucose_uptake,
                          color = treatment)) +
    
    # points
    geom_sina(alpha = 0.5,
              position = pd) +
    
    # plot means and CI
    geom_errorbar(data = fig_2g_m1_emm_dt[tissue == tissue_i],
                  aes(y = emmean,
                      ymin = lower.CL,
                      ymax = upper.CL,
                      color = treatment),
                  width = 0,
                  position = pd
    ) +
    
    geom_point(data = fig_2g_m1_emm_dt[tissue == tissue_i],
               aes(y = emmean,
                   color = treatment),
               size = 3,
               position = pd
    ) +
    
    # plot p-values (y positions are adjusted by eye)
    stat_pvalue_manual(fig_2g_m1_pairs_dt[tissue == tissue_i],
                       label = "p_pretty",
                       y.position=c(y_pos),
                       tip.length = 0.01) +
    
    # aesthetics
    ylab("Glucose Uptake") +
    scale_color_manual(values=pal_nature_mod,
                       name = NULL) +
    ggtitle(tissue_i)+
    
    theme_pubr() +
    theme(legend.position="top") +
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank()) +
    
    NULL
  
  if(y_lab == FALSE){
    gg <- gg + theme(axis.title.y = element_blank())
  }
  if(g_lab == FALSE){
    gg <- gg + theme(legend.position="none")
  }
  
  return(gg)
}

y_cols <- c("ingWAT", "epiWAT", "Muscle", "BAT")
legend <- get_legend(fig_2g_plot("ingWAT", y_lab = TRUE, g_lab = TRUE))
gg1 <- fig_2g_plot("ingWAT", y_lab = TRUE, )
gg2 <- fig_2g_plot("epiWAT")
gg3 <- fig_2g_plot("Muscle")
gg4 <- fig_2g_plot("BAT")
top_gg <- plot_grid(gg1, gg2, gg3, gg4,
                       nrow=1,
                       rel_widths = c(1.15, 1, 1.05, 1.1)) # by eye
plot_grid(top_gg, legend,
          nrow=2,
          rel_heights = c(1, 0.1))

1.8 Figure 2h

1.9 Figure 2i

1.10 Figure 2j

1.11 Exercises

  1. Using the code chunks for Figure 2f, complete the analysis for Figure 2i. This should simply require copying the chunks from 2f and swapping out the correct names (and the range of the Excel file!) of variables. At a minimum you will need to run
  • the setup chunk with the libraries
  • the chunk defining the path of the data folder and the sheet with Figure 2
  • the chunk defining the pal_nature_mod palette in the “useful functions” section