figure-03-50k-sim.Rmd

# Figure 3

> **_NOTE:_** This plot omits the non-human primate CFU and flowcytometry data
> scatter plot layer shown in the publication because we do not have permission
> to redistribute that data.

```{r imports}
#| cache: false
#| message: false
library(readr)     # read_csv
library(dplyr)
library(snakecase) # to_snake_case
library(tidyr)     # pivot_longer
library(forcats)   # fct
library(stringr)   # str_remove
library(ggplot2)
library(magick)    # image_read
library(grid)      # rasterGrob
library(ggpp)      # geom_grob
library(patchwork)

source("R/config.R")
```

```{r prepare-inputs-temporal}
#| cache: true
file_stat <-
  file.path(dir_data,
            "2024-07-19-A-gr-50k",
            "output-lhs-stat-cols-35.csv.gz")

df_runs <-
  ## Read in the data with minimal changes.
  read_csv(file_stat,
           col_types = c(exp = "c", time = "i"),
           show_col_types = FALSE) %>%
  ## Convert ticks to days.
  mutate(time = time / ticks_per_day)
df_runs

## Labels to remap the outputs for plotting.
outputs <-
  c(mac_scaled_ring = "Macrophages",
    tc_scaled_ring = "T cells",
    tot_mtb_in_gran_scaled_ring = "CFU")
df_std <-
  df_runs %>%
  ## Remove spaces and dots from column names.
  rename_with(to_snake_case) %>%
  ## Classify by Mtb level.
  group_by(exp) %>%
  mutate(classif =
         case_when(
           tot_mtb_in_gran_scaled_ring[time == max(time)] == 0
           ~ "Sterile",
           (tot_mtb_in_gran_scaled_ring[time == 154] -
            tot_mtb_in_gran_scaled_ring[time == 105]) <
             1.1 * tot_mtb_in_gran_scaled_ring[time == 105] |
           tot_mtb_in_gran_scaled_ring[time == 154] < 1e4
           ~ "Controlling",
           .default = "Uncontrolled"
         ) %>%
         fct(levels = c("Uncontrolled", "Sterile", "Controlling")),
         .before = time) %>%
  ungroup() %>%
  ## Sum T cell counts.
  mutate(exp = fct(exp),
         tc_scaled_ring = c(tgam_scaled_ring +
                              tcyt_scaled_ring +
                              treg_scaled_ring)) %>%
  ## Subset to relevant data.
  dplyr::select(exp:time, ends_with("ring")) %>%
  ## Convert to long format for ploting.
  pivot_longer(-c(exp:time), names_to = "output") %>%
  dplyr::filter(output %in% names(outputs)) %>%
  ## Use labels.
  mutate(output = outputs[output])
df_std
df_std |>
  distinct(exp, classif) |>
  count(classif) |>
  mutate(percent = 100 * n / sum(n),
         classif = fct_relevel(classif, "Sterile", after = Inf)) |>
  arrange(classif)
```

```{r prepare-inputs-proprietary}
#| cache: true
#| dependson: -1
## Files containing proprietary calibration data.  These a may or may not be
## present when generating the plot.
dir_data_proprietary <-
  file.path("~/immunology/GR-ABM-ODE/simulation/scripts",
            "calibration/abc/src/gransimcal/data")
file_cfu <-
  file.path(dir_data_proprietary, "CombinedCFU_updateUponMoreData.csv")
file_flow <-
  file.path(dir_data_proprietary, "ModelCalibrationCountLung_2018.csv")
## Read in reference CFU and flowcytometry datasets.
if (all(file.exists(file_cfu, file_flow))) {
  df_cfu <-
    read_csv(file_cfu, col_select = c(3, 1, 2), show_col_types = FALSE) %>%
    rename_with(to_snake_case) %>%
    ## Remove single outlier diameter.
    dplyr::filter(granuloma_diameter_mm >= 0.5) %>%
    select(-granuloma_diameter_mm) %>%
    pivot_longer(-time, names_to = "output") %>%
    ## Ignore sterilization.
    dplyr::filter(value > 0)
  df_flow <-
    read_csv(file_flow, show_col_types = FALSE) %>%
    ## Map columns to GranSim agents.
    dplyr::rename(time = Timepoint,
                  tc_scaled_ring = "CD3") %>%
    mutate(mac_scaled_ring = `CD11b+` - `CD11b+/Calp+`) %>%
    ## Remove untreated control.
    dplyr::filter(time > 0) %>%
    dplyr::select(time, any_of(names(outputs))) %>%
    pivot_longer(-time, names_to = "output")
}
```

```{r plot-temporal}
#| cache: true
#| dependson:
#| - -2
#| - -1
#| fig.width: 6.5
#| fig.height: 3
colors <- palette.colors(4, palette = "Okabe-Ito")[-1]
plot_time <-
  df_std %>%
  mutate(time = time / 7) %>%
  ggplot(aes(x = time,
             y = value)) +
  facet_grid(~ output,
             scales = "free_y") +
  labs(x = "Week (post-infection)",
       y = "Cell counts",
       color = "Classification") +
  theme_bw() +
  theme(legend.position = "top",
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(0, -10, -10, -10)) +
  ## The subsequent geom_line(linewidth = 0.1, ...) setting will make the
  ## legend harder to see, so reset the legend to the default linewidth.
  guides(color = guide_legend(override.aes = list(linewidth = 1))) +
  scale_y_log10(labels = scales::label_log()) +
  ## Reorder the legend to match how the classifications appear in the plot.
  scale_color_manual(breaks = c("Uncontrolled", "Controlling", "Sterile"),
                     values = colors) +
  geom_line(linewidth = 0.1,
            aes(group = exp,
                color = classif))
## Overlay the proprietary data if it is accessible.
if (all(c("df_cfu", "df_flow") %in% ls())) {
  plot_time <-
    plot_time +
    geom_point(data =
                 bind_rows(df_cfu, df_flow) %>%
                 mutate(output = outputs[output],
                        time = time / 7),
               size = 0.5)
}
plot_time
```

```{r plot-spatial}
#| cache: true
#| fig.width: 6.5
#| fig.height: 2.3
files_spatial_meta <-
  file.path(dir_data,
            "2025-03-02-A-matches-spatial-x25-each-fig-2024-07-19-A-gr-50k",
            "runlist-metadata.csv")
files_spatial <-
  file.path(dir_data,
            "2025-03-02-A-matches-spatial-x25-each-fig-2024-07-19-A-gr-50k",
            "exp*",
            "exp*",
            "exp*umpng") |>
  Sys.glob()
files_thumbnail <-
  file.path(dir_data,
            "2025-03-02-A-matches-spatial-x25-each-fig-2024-07-19-A-gr-50k",
            "images",
            "*.png") |>
  Sys.glob()

## Plot grviz-lung rendered PNG images with legend thumbnails.
plot_space <-
  read_csv(files_spatial_meta, col_types = c(classif = "f", ntile = "i")) |>
  right_join(
    tibble(path = files_spatial,
           ## grviz-lung fixes the calibration bar size as 1/3rd the image
           ## width.
           cal =
             files_spatial |>
             basename() |>
             str_extract("(\\d+[.]\\d+)um", group = 1) |>
             as.numeric(),
           exp =
             path |>
             basename() |>
             str_extract("^exp(\\d+)", group = 1) |>
             str_c(".1") |>
             as.numeric()),
    by = "exp") |>
  mutate(classif = fct_relabel(classif, function(x) {
    sprintf("%s\nCFU = %s", x, format(round(CFU), trim = TRUE, big.mark = ","))
  })) |>
  ggplot() +
  lims(x = c(0, 512),
       y = c(0, 512)) +
  theme_bw() +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank()) +
  coord_fixed(expand = FALSE) +
  facet_wrap(~ classif) +
  geom_grob(aes(x = 256,
                y = 256,
                label = img),
            data =
              . %>%
              mutate(img =
                       path |>
                       lapply(image_read) |>
                       lapply(rasterGrob)),
            vp.width = 1,
            vp.height = 1) +
  geom_line(data =
              . %>%
              filter(classif == classif[1]) %>%
              group_by(classif) %>%
              reframe(x =
                        ## grviz-lung fixes the calibration bar to 0.3 of the
                        ## image width.  Re-normalize to 500um.
                        -c(0, 0.3 * cal * (500 / cal)) +
                        quantile(c(0, 512), 0.95),
                      y = quantile(c(0, 512), 0.05)),
            aes(x = x, y = y, group = classif),
            color = "white")
plot_space

names <-
  c(## Macropahges.
    "Resting" = "mac-resting.png",
    "Activated" = "mac-active.png",
    "Infected" = "mac-infected.png",
    "Chronically~infected" =  "mac-chronically-infected.png",
    ## T cells.
    "Cytotoxic" = "tcyt.png",
    "Regulatory" = "treg.png",
    "IFN~gamma~producing" = "tgam.png",
    ## Mtb.
    "Mtb" = "mtb-non-square.png",
    ## Grid.
    "Caseation" = "caseation.png")
stopifnot(all(names %in% basename(files_thumbnail)))
## Plot legend tumbnails.
df_thumb <-
  tibble(cell = names(names),
         basename = setNames(names, NULL)) |>
  left_join(tibble(path = files_thumbnail,
                   basename = basename(path)),
            by = "basename") |>
  mutate(group = case_when(
    str_detect(basename, "^mac-") ~ "Macrophages",
    str_detect(basename, "^t.{3}[.]") ~ "T cells",
    str_detect(basename, "^mtb") ~ "CFU",
    .default = "Grid"),
    .before = 1L) |>
  mutate(image =
           path |>
           lapply(image_read) |>
           lapply(image_modulate, brightness = 200) |>
           lapply(rasterGrob)) |>
  select(-basename,
         -path)
layout <- "
A#FJ
BDGK
CEHL
##IM
####
"
plot_titles <-
  lapply(c(1, 5, 8, 9), function(i) {
    wrap_elements(panel = textGrob(df_thumb[i, "group"],
                                   x = unit(0, "npc"),
                                   just = "left",
                                   gp = gpar(fontsize = 8.5)),
                  clip = FALSE) +
      theme(plot.tag = element_blank())
  })
plot_thumbnails <-
  lapply(1:9, function(i) {
    ggplot() +
      geom_grob(aes(x = 0L,
                    y = 0L,
                    label = image),
                data = df_thumb[i, "image"],
                vp.width = 1,
                vp.height = 1) +
      labs(y = parse(text = pull(df_thumb[i, "cell"]))) +
      coord_fixed() +
      theme(plot.tag = element_blank())
  })
plot_space_legend <-
  ## Macrophages.
  plot_titles[[1]] +
  plot_thumbnails[[1]] +
  plot_thumbnails[[2]] +
  plot_thumbnails[[3]] +
  plot_thumbnails[[4]] +
  ## T cells.
  plot_titles[[2]] +
  plot_thumbnails[[5]] +
  plot_thumbnails[[6]] +
  plot_thumbnails[[7]] +
  ## CFU.
  plot_titles[[3]] +
  plot_thumbnails[[8]] +
  ## Grid.
  plot_titles[[4]] +
  plot_thumbnails[[9]] +
  plot_layout(design = layout) &
  ## Place y-axis on right per https://stackoverflow.com/a/54828826
  scale_y_continuous(position = "right",
                     sec.axis = dup_axis()) &
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y.left = element_blank(),
        axis.title.y.right = element_text(angle = 0,
                                          hjust = 0,
                                          vjust = 0.5,
                                          size = 8.5),
        plot.margin = unit(rep(0, 4), "cm"))
plot_space_legend
```

```{r figure-03-50k-sim}
#| cache: true
#| dependson:
#| - -2
#| - -1
#| fig.width: 6.5
#| fig.height: 6
(plot_time & theme(plot.tag.position = c(0.03, 0.88))) /
  plot_spacer() /
  (plot_space & theme(plot.tag.position = c(0.03, 0.95))) /
  plot_space_legend +
  plot_layout(heights = c(1, 0.1, 1, 0.4)) +
  plot_annotation(tag_levels = "A") &
  theme(plot.margin = margin(0.1))
```

```{r session-info}
#| cache: false
sessionInfo()
```