Oxidative Stress

Analyzing the effect of Condition (Hypoxia vs Normoxia) on a Oxidative Stress gene panel

source(here::here("src", "setup.R"), echo = FALSE)

## Since downlit doesn't seem to link packages without an explicit & visible library() call anymore ...

library(here)
library(pipebind)

library(dplyr)
library(tidyr)
library(stringr)
library(purrr)
library(ggplot2)

library(glmmTMB)
library(parameters)
library(insight)
library(performance)
library(emmeans)

library(DT)
library(plotly)

## Data

supplementary_data <- load_supplementary_data()

OS_data <- load_pcr_data("OS")

1 Clean Data

Here’s the data from our Oxidative Stress gene panel, which has been cleaned by:

• Removing missing DCq values
• Removing genes missing either N or IH Condition
• Removing genes having less than 3 data points in either the N or IH Condition

2 Model fitting & diagnostics

2.1 Model fitting

Let’s define the model we will fit to each Gene’s DCq data.

Here, we will fit a Linear Mixed Effect Model, with a random effect per Experiment, to account for potential clustering:

\[ \begin{aligned} DCq_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1}(Condition), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for Experiment j = 1,} \dots \text{,J} \end{aligned} \]

OS_model <- function(data) {
    form <- "DCq ~ Condition"

    # If the data contains an Experiment variable with more than two levels, add a random effect by Experiment
    if ("Experiment" %in% colnames(data) && n_distinct(data$Experiment) >= 2) {
        form <- str_c(form, " + (1 | Experiment)")
    }

    glmmTMB(
        as.formula(form),
        family = gaussian("identity"),
        data = data,
        contrasts = list(Condition = "contr.sum")
    )
}

Now, let’s fit said model to each Gene’s data, for a given Stage and Layer:

compute_fold_change <- function(mod) {
    return(
        get_data(mod) |>
            select(Condition, DCq) |>
            pivot_wider(
                names_from = Condition,
                values_from = DCq,
                values_fn = mean
            ) |>
            summarize(Fold = 2**(-1 * (IH - N))) |>
            pull(Fold)
    )
}

(OS_data$models <- OS_data$clean |>
    group_split(Stage, Layer, Gene) |>
    map_dfr(
        \(d) {
            summarize(
                d,
                Mod = pick(Experiment, Condition, DCq) |> OS_model() |> list(),
                .by = c(Stage, Layer, Gene)
            )
        },
        .progress = "Fitting models:"
    ) |>
    filter(!has_na_coefs(Mod)) |> # Removing models that did not fit properly
    mutate(Fold = map_dbl(Mod, compute_fold_change)) |> # Adding the Fold change
    select(Stage, Layer, Gene, Fold, Mod))

2.2 Model diagnostics

Here are the diagnostic plots for a random sample of the fitted models:

3 Model analysis

For each model we fit, we can then extract the CI and p.value for the relevant contrasts, and use those to establish if a Gene was up or down-regulated:

get_emmeans_data <- function(mod) {
    return(
        emmeans(mod, specs = "Condition", type = "response") |>
            contrast(method = "pairwise", adjust = "none", infer = TRUE) |>
            as.data.frame() |>
            pivot_wider(names_from = contrast, values_from = estimate) |>
            select(last_col(), LCB = lower.CL, UCB = upper.CL, p.value) |>
            mutate(across(where(is.character), \(x) na_if(x, "NaN")))
    )
}

get_regulation_type <- function(fold, p_value) {
    case_when(
        p_value <= .05 & fold < 1 ~ regulation_type$DOWNREG,
        p_value <= .05 & fold > 1 ~ regulation_type$UPREG,
        is.na(p_value) | is.na(fold) ~ NA_character_,
        .default = regulation_type$NOT_REG
    )
}

(OS_data$predictions <- OS_data$models |>
    group_split(Stage, Layer, Gene) |>
    map_dfr(
        \(d) mutate(d, get_emmeans_data(Mod[[1]])),
        .progress = "Extracting model predictions:"
    ) |>
    filter(!is.na(p.value)) |>
    mutate(Expression = get_regulation_type(Fold, p.value)) |>
    select(
        Stage,
        Layer,
        Gene,
        Fold,
        Expression,
        matches("-|/"),
        LCB,
        UCB,
        p.value
    ))

3.1 Gene regulation overview

We can get a general overview of which Gene are up or down-regulated through a Sunburst plot, stacked by Stage and Layer.

make_suburst_plot <- function(dat, layers, tooltips = NULL, colors = NULL, plot_options = list()) {
    if (!is.null(colors)) {
        colors <- set_names(colors, \(x) x |> str_replace_all(" ", "\n") |> str_replace_all("and\n", "and "))
    }

    .make_sunburst_data <- function(dat, layers_part) {
        if (all(layers == layers_part)) {
            sun_dat <- summarize(dat, across(any_of(c(layers_part, tooltips)), first), .by = all_of(layers_part))
        } else {
            sun_dat <- summarize(dat, across(any_of(layers_part), first), .by = all_of(layers_part))
        }

        return(
            bind_cols(
                transmute(sun_dat, ids = pick(any_of(layers_part)) |> pmap_chr(\(...) str_c(..., sep = ">>"))),
                transmute(
                    sun_dat,
                    parents = pick(any_of(layers_part)) |>
                        select(-last_col()) |>
                        pmap_chr(\(...) str_c(..., sep = ">>")) |>
                        bind(x, if (is_empty(x)) "" else x)
                ),
                transmute(sun_dat, labels = pick(any_of(layers_part)) |> select(last_col()) |> pull(1))
            ) |>
                bind(
                    x,
                    if (all(layers == layers_part) && !is.null(tooltips)) {
                        bind_cols(
                            x,
                            transmute(
                                sun_dat,
                                hovertext = pmap_chr(
                                    pick(all_of(tooltips)),
                                    \(...) {
                                        str_c(
                                            names(c(...)),
                                            map_if(list(...), is.numeric, \(x) round(x, 3)) |>
                                                map_if(is.factor, as.character),
                                            sep = ": ",
                                            collapse = "\n"
                                        )
                                    }
                                )
                            )
                        )
                    } else {
                        x
                    }
                ) |>
                mutate(labels = labels |> str_replace_all(" ", "\n") |> str_replace_all("and\n", "and "))
        )
    }

    params <- map_dfr(seq_along(layers), \(x) .make_sunburst_data(dat, layers[1:x])) |>
        mutate(bg_color = map_chr(labels, \(x) colors[[x]] %||% NA))

    rlang::exec(
        plotly::plot_ly,
        !!!select(params, -bg_color),
        type = "sunburst",
        branchvalues = "total",
        hoverinfo = "text",
        sort = FALSE,
        marker = list(colors = ~ params$bg_color),
        !!!plot_options
    )
}

OS_data$predictions |>
    filter(Expression != regulation_type$NOT_REG) |>
    left_join(supplementary_data$gene_data$OS, join_by(Gene)) |>
    make_suburst_plot(
        layers = c("Stage", "Layer", "Expression", "Gene"),
        tooltips = c("Gene", "Pathway", "Fold", "p.value"),
        colors = sunburst_pcr_colors,
        plot_options = list(insidetextorientation = 'radial')
    )

3.2 Gene regulation timeline

To get a better idea of how each Gene’s regulation changes through time, we can plot a timeline of their expression, split by Layer and Pathway, with a color coding of the Effect of this regulation.

make_fold_timeline_plot <- function(
    dat,
    facet_rows = "Pathway",
    trans = "identity",
    color_by = NULL,
    colors = colors_effect,
    size_boost = 1
) {
    origin <- do.call(trans, list(1))

    dat <- (dat |>
        mutate(Fold_trans = do.call(trans, list(Fold))) |>
        mutate(
            Fold_Amp = ifelse(
                max(Fold_trans, na.rm = TRUE) - min(Fold_trans, na.rm = TRUE) != 0,
                max(Fold_trans, na.rm = TRUE) - min(Fold_trans, na.rm = TRUE),
                mean(Fold_trans, na.rm = TRUE)
            ) *
                0.1,
            .by = all_of(c(facet_rows, "Stage"))
        ))

    timeline <- (ggplot(dat) +
        {
            if (is.null(color_by)) aes(x = Gene, color = Fold >= 1) else aes(x = Gene, color = .data[[color_by]])
        } +
        geom_linerange(aes(ymax = Fold_trans), ymin = origin, linewidth = 2 + (size_boost * 0.5)) +
        geom_hline(yintercept = origin, linewidth = 0.3, linetype = "dotted") +
        geom_text(
            aes(
                label = str_c(round(Fold, 2), stars.pval(p.value) |> str_replace(fixed("."), ""), sep = " "),
                y = ifelse(Fold_trans > origin, Fold_trans + Fold_Amp, Fold_trans - Fold_Amp),
                hjust = ifelse(Fold > 1, 0, 1)
            ),
            vjust = 0.5,
            angle = 0,
            size = 2 + (size_boost * 0.25),
            check_overlap = TRUE
        ) +
        scale_color_manual(" ", values = colors) +
        scale_y_continuous(breaks = c(0, 1, 2, 3), expand = expansion(mult = 1.01 * (1 + (size_boost / 100)))) +
        scale_x_discrete(expand = expansion(add = 1 * size_boost), limits = \(x) rev(x)) +
        labs(
            x = "",
            y = ifelse(trans != "identity", str_glue("Fold Change *({trans} scale)*"), "Fold Change")
        ) +
        coord_flip() +
        facet_grid(
            vars(.data[[facet_rows]]),
            vars(Stage),
            scales = "free_y",
            space = "free_y",
            labeller = label_wrap_gen(width = 12, multi_line = TRUE)
        ) +
        {
            if (!is.null(color_by)) guides(color = guide_legend(title = color_by))
        } +
        theme(
            legend.position = ifelse(is.null(color_by), "none", "bottom"),
            axis.text.x = element_blank(),
            axis.title.x = element_markdown(size = 9),
            axis.text.y = element_text(size = 7),
            strip.text = element_text(size = 5 * size_boost),
            plot.title = element_markdown(size = 9, face = "plain", vjust = 1, hjust = 0.5)
        ))

    return(timeline)
}

render_OS_timeline <- function(dat, group, size_boost = 1.2) {
    if (group[[1]][1] == "Whole") {
        dat |>
            group_by(Figure) |>
            group_map(\(d, g) {
                render_OS_timeline(
                    d,
                    str_glue("Whole Cerebellum - {first(g[[1]])}"),
                    1.5
                )
            })
    } else {
        group_name <- group[[1]][1]
        if (group_name == "PC") {
            group_name <- "Purkinje Cells"
        }

        if (group_name == "Cellular Response") {
            cell_resp_levels <- c(
                "Cell Death and Protection",
                "Apoptotic Pathways",
                "Inflammation Pathways",
                "Autophagy and Mitophagy"
            )
            dat <- mutate(dat, Pathway = factor(Pathway, levels = cell_resp_levels))
        }

        plot <- make_fold_timeline_plot(
            dat,
            facet_rows = "Pathway",
            trans = "log",
            color_by = "Effect",
            size_boost = size_boost
        )

        width <- 2 + n_distinct(dat$Stage)
        height <- dat |>
            group_by(Pathway) |>
            group_map(\(d, g) n_distinct(d$Gene) * 0.1 + 1) |>
            flatten_dbl() |>
            sum()

        template_md <- c(
            '### `r group_name`',
            '```{r}',
            '#| echo: false',
            '#| fig-width: !expr width',
            '#| fig-height: !expr height',
            'plot',
            '```'
        )

        knitr::knit_child(text = template_md, envir = rlang::env(), quiet = TRUE)
    }
}

OS_timelines <- (OS_data$predictions |>
    left_join(supplementary_data$gene_data$OS, join_by(Gene)) |>
    filter(p.value <= .05) |>
    mutate(
        Effect = case_when(
            str_detect(Expression, "Downregulated") &
                Effect == "Beneficial" ~ "Deleterious",
            str_detect(Expression, "Downregulated") &
                Effect == "Deleterious" ~ "Beneficial",
            .default = Effect
        )
    ) |>
    select(
        Stage,
        Layer,
        Gene,
        Fold,
        p.value,
        Expression,
        Effect,
        Pathway,
        Figure
    ) |>
    group_by(Layer) |>
    group_map(render_OS_timeline) |>
    unlist())

Purkinje Cells

Whole Cerebellum - Cellular Response

Whole Cerebellum - Oxidative Stress