ISM.Rmd

---
title: "In-Silico Mutagenesis (ISM) — Yeast & Human"
author: "Emmanuel Cazottes"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    toc: true
    toc_float: true
    theme: united
    highlight: tango
    code_folding: show
    df_print: paged
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.width = 10,
  fig.height = 6,
  fig.align = "center"
)

library(stringr)
library(tidyverse)
library(reshape2)
library(ggpubr)
library(ggsci)
library(gridExtra)

publication_theme <- theme_pubr() +
  theme(
    axis.text = element_text(size = 7, color = "black"),
    axis.title = element_text(size = 10, face = "plain"),
    plot.title = element_text(size = 10, face = "plain", hjust = 0.5),
    legend.title = element_text(size = 10, face = "plain"),
    legend.text = element_text(size = 7),
    axis.line = element_line(size = 0.5, color = "black"),
    axis.ticks = element_line(size = 0.5, color = "black"),
    axis.ticks.length = unit(0.2, "cm"),
    legend.position = "right",
    legend.box.background = element_rect(color = "white"),
    legend.background = element_rect(fill = "white", color = NA),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_rect(fill = "white")
  )

theme_set(publication_theme)
```

Two complementary ISM analyses:

1. **Oracle ISM** — ISM scores for top-n-bottom control sequences only,
   comparing Rarely vs Often selected within each AL strategy (uncertainty /
   diversity).
2. **Full-pool ISM** — ISM scores for all AL pool sequences, with
   top-n-bottom annotations; compared to AL pool and a random baseline via
   Mann–Whitney U.


# Shared utility functions

```{r shared-functions}
#' Format a p-value for plotting (floor at 2.2e-16).
fmt_pval <- function(p) {
  ifelse(p <= 2.2e-16, "p<2.2e-16", paste0("p=", signif(p, 2)))
}


#' Rename raw control method labels using regex (works for both organisms).
rename_tnb_labels <- function(method_vec) {
  method_vec <- gsub("uncertainty_((1x60)|(3x20))_bottom\\d+", "Rarely uncertain", method_vec)
  method_vec <- gsub("uncertainty_((1x60)|(3x20))_top\\d+",    "Often uncertain",  method_vec)
  method_vec <- gsub("diversity_((1x60)|(3x20))_bottom\\d+",   "Rarely diverse",   method_vec)
  method_vec <- gsub("diversity_((1x60)|(3x20))_top\\d+",      "Often diverse",    method_vec)
  method_vec
}

LEVEL_ORDER <- c("AL pool", "Random selection",
                  "Rarely uncertain", "Often uncertain",
                  "Rarely diverse", "Often diverse")
INFO_ORDER  <- c("AL pool", "Random selection", "Rarely", "Often")

FILL_COLOURS <- c("Rarely" = "#F3756D", "Often" = "#1EBDC2",
                   "Random selection" = "white", "AL pool" = "white")
```


## Oracle ISM helpers

```{r oracle-helpers}
#' Load and annotate oracle ISM results.
#'
#' The oracle ISM file has columns: mean, max, group.  Group encodes the AL
#' strategy + pool size (e.g. "uncertainty_3x20", "diversity_1x60").  The
#' mapping from group order to human-readable labels is fixed by the upstream
#' pipeline.
load_oracle_ism <- function(path) {
  ism <- read.table(path, header = TRUE, sep = "\t")

  label_order <- c(rep(c("Rarely uncertain", "Rarely diverse"), 2),
                   rep(c("Often uncertain", "Often diverse"), 2))

  for (i in seq_along(unique(ism$group))) {
    ism[ism$group == unique(ism$group)[i], "al.method"] <- label_order[i]
  }

  ism$sel.seq <- factor(str_split_fixed(ism$al.method, " ", 2)[, 1],
                         levels = c("Rarely", "Often"))
  ism$round   <- str_split_fixed(ism$group, "_", 2)[, 2]
  ism$group   <- factor(str_split_fixed(ism$group, "_", 2)[, 1],
                         levels = c("uncertainty", "diversity"))
  ism$al.method <- factor(ism$al.method,
                           levels = c("Rarely uncertain", "Often uncertain",
                                      "Rarely diverse", "Often diverse"))
  ism
}


#' Mann–Whitney tests for oracle ISM (Rarely vs Often within each group).
#'
#' @param ism        Oracle ISM data.frame.
#' @param y.mean     Y position for mean p-value labels.
#' @param y.max      Y position for max p-value labels.
oracle_mw_test <- function(ism, y.mean, y.max) {
  map_dfr(c("3x20", "1x60"), function(r) {
    wt <- function(metric, al1, al2) {
      wilcox.test(subset(ism, round == r & al.method == al1)[[metric]],
                  subset(ism, round == r & al.method == al2)[[metric]])
    }

    data.frame(
      round   = r,
      group   = rep(c("diversity", "uncertainty"), 2),
      ism     = rep(c("mean", "mean", "max", "max"), 2),
      y.axis  = rep(c(y.mean, y.mean, y.max, y.max), 2),
      p.value = fmt_pval(c(
        wt("mean", "Rarely diverse",   "Often diverse")$p.value,
        wt("mean", "Rarely uncertain", "Often uncertain")$p.value,
        wt("max",  "Rarely diverse",   "Often diverse")$p.value,
        wt("max",  "Rarely uncertain", "Often uncertain")$p.value
      ))
    )
  })
}


#' Dodged violin/boxplot for oracle ISM.
plot_oracle_ism <- function(ism, mwt, round_filter, metric) {
  ylab.text <- ifelse(metric == "mean",
                       "Mean absolute ISM score", "Max ISM score")

  ggplot() +
    geom_violin(data = subset(ism, round == round_filter),
                aes(x = group, y = .data[[metric]], fill = sel.seq),
                width = 0.75, position = position_dodge(width = 0.5)) +
    geom_boxplot(data = subset(ism, round == round_filter),
                 aes(x = group, y = .data[[metric]], fill = sel.seq),
                 width = 0.1, outliers = FALSE,
                 position = position_dodge(width = 0.5)) +
    geom_text(data = subset(mwt, round == round_filter & ism == metric),
              aes(x = group, y = y.axis, label = p.value), size = 3) +
    xlab("") + ylab(ylab.text) + labs(fill = "") +
    theme(aspect.ratio = 1,
          plot.margin = margin(5, 5, 5, 5, "mm"),
          legend.position = "none")
}
```


## Full-pool ISM helpers

```{r pool-helpers}
#' Load full-pool ISM and merge with sequence metadata.
#'
#' @param path.ism    Path to the ISM results TSV (columns: seq, mean, max).
#' @param path.pool   Path to the AL pool sequence TSV.
#' @param path.tnb    Path to the top-n-bottom sequence TSV.
#' @param n.rnd       Number of random sequences to sample from the pool.
#' @param trim.seq    Optional: trim padding from sequences (yeast-specific).
#'                    A length-2 integer vector c(start, end) for str_sub,
#'                    or NULL to skip.
load_pool_ism <- function(path.ism, path.pool, path.tnb,
                           n.rnd, trim.seq = NULL) {

  ism.all <- read.table(path.ism, header = TRUE, sep = "\t")
  if (!is.null(trim.seq)) {
    ism.all$seq <- str_sub(ism.all$seq, trim.seq[1], trim.seq[2])
  }

  al.pool <- read.table(path.pool, header = FALSE, sep = "\t",
                          col.names = c("sequence_name", "seq", "method"))
  tnb     <- read.table(path.tnb, header = FALSE, sep = "\t",
                          col.names = c("sequence_name", "seq", "method"))

  ism.pool <- merge(al.pool, ism.all, by = "seq")
  ism.tnb  <- merge(tnb, ism.all, by = "seq")

  set.seed(42)
  ism.rnd <- ism.pool[sample(nrow(ism.pool), n.rnd), ]
  ism.rnd$method    <- "Random selection"
  ism.rnd$al.method <- "Random selection"
  ism.rnd$rounds    <- "Random selection"
  ism.rnd$info      <- "Random selection"

  # Annotate tnb
  ism.tnb$al.method <- rename_tnb_labels(ism.tnb$method)
  ism.tnb$rounds    <- str_split_fixed(ism.tnb$method, "_", 3)[, 2]
  ism.tnb$info      <- str_split_fixed(ism.tnb$al.method, " ", 2)[, 1]

  # Annotate pool
  ism.pool$al.method <- "AL pool"
  ism.pool$rounds    <- "AL pool"
  ism.pool$info      <- "AL pool"

  out <- rbind(ism.tnb, ism.pool, ism.rnd)
  out$al.method <- factor(out$al.method, levels = LEVEL_ORDER)
  out$info      <- factor(out$info, levels = INFO_ORDER)
  out
}


#' Mann–Whitney tests for full-pool ISM (each group vs AL pool).
pool_mw_test <- function(ism, y.mean, y.max) {
  ref.mean <- subset(ism, method == "ALpool")$mean
  ref.max  <- subset(ism, method == "ALpool")$max

  df.pval <- map_dfr(unique(ism$method), function(m) {
    sub <- subset(ism, method == m)
    mwt.mean <- wilcox.test(sub$mean, ref.mean)
    mwt.max  <- wilcox.test(sub$max, ref.max)

    data.frame(
      method         = m,
      al.method      = unique(sub$al.method),
      rounds         = unique(sub$rounds),
      info           = unique(sub$info),
      pvalue.mean    = fmt_pval(mwt.mean$p.value),
      statistic.mean = mwt.mean$statistic,
      pvalue.max     = fmt_pval(mwt.max$p.value),
      statistic.max  = mwt.max$statistic,
      y.mean         = y.mean,
      y.max          = y.max
    )
  })
  df.pval
}


#' Violin/boxplot for full-pool ISM.
plot_pool_ism <- function(ism, mwt, round_filter, metric, title) {
  ylab.text <- ifelse(metric == "mean",
                       "Mean absolute ISM score", "Max ISM score")
  pval.col <- ifelse(metric == "mean", "pvalue.mean", "pvalue.max")
  y.col    <- ifelse(metric == "mean", "y.mean", "y.max")

  ggplot(subset(ism, rounds %in% c(round_filter, "AL pool", "Random selection")),
         aes(x = al.method, y = .data[[metric]], fill = info)) +
    geom_violin(width = 1.25, outliers = FALSE) +
    geom_boxplot(width = 0.1, outliers = FALSE) +
    geom_text(data = mwt,
              aes(x = al.method, y = .data[[y.col]],
                  label = .data[[pval.col]]), size = 3) +
    scale_fill_manual(values = FILL_COLOURS) +
    xlab("") + ylab(ylab.text) + ggtitle(title) +
    theme(aspect.ratio = 0.75,
          axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
          plot.margin = margin(5, 5, 5, 5, "mm"),
          legend.position = "none")
}
```


# Oracle ISM {.tabset}

ISM scores computed on top-n-bottom control sequences only.  Mann–Whitney U
compares Rarely vs Often selected within each strategy.

## Data loading

```{r oracle-load}
oracle.yeast <- load_oracle_ism("yeast/sequence_analysis/yeast_ism_all.tsv")
oracle.human <- load_oracle_ism("human/sequence_analysis/human_ism_all.tsv")

cat("=== Yeast ===\n")
oracle.yeast %>% group_by(round, al.method) %>%
  summarise(avg_mean = mean(mean), median_mean = median(mean), .groups = "drop")

cat("\n=== Human ===\n")
summary(oracle.human$mean)
summary(oracle.human$max)
```

## Mann–Whitney tests

```{r oracle-mwt}
oracle.mwt.yeast <- oracle_mw_test(oracle.yeast, y.mean = 1.2, y.max = 9)
oracle.mwt.human <- oracle_mw_test(oracle.human, y.mean = 0.6, y.max = 6.5)
```

## Plots

```{r oracle-plots, fig.width=4, fig.height=3}
# Yeast
p.oracle.y.3x20.mean <- plot_oracle_ism(oracle.yeast, oracle.mwt.yeast, "3x20", "mean")
p.oracle.y.3x20.max  <- plot_oracle_ism(oracle.yeast, oracle.mwt.yeast, "3x20", "max")
p.oracle.y.1x60.mean <- plot_oracle_ism(oracle.yeast, oracle.mwt.yeast, "1x60", "mean")
p.oracle.y.1x60.max  <- plot_oracle_ism(oracle.yeast, oracle.mwt.yeast, "1x60", "max")

p.oracle.y.3x20.mean; p.oracle.y.3x20.max
p.oracle.y.1x60.mean; p.oracle.y.1x60.max

# Human
p.oracle.h.3x20.mean <- plot_oracle_ism(oracle.human, oracle.mwt.human, "3x20", "mean")
p.oracle.h.3x20.max  <- plot_oracle_ism(oracle.human, oracle.mwt.human, "3x20", "max")
p.oracle.h.1x60.mean <- plot_oracle_ism(oracle.human, oracle.mwt.human, "1x60", "mean")
p.oracle.h.1x60.max  <- plot_oracle_ism(oracle.human, oracle.mwt.human, "1x60", "max")

p.oracle.h.3x20.mean; p.oracle.h.3x20.max
p.oracle.h.1x60.mean; p.oracle.h.1x60.max
```

```{r oracle-save}
dir.create("fig", showWarnings = FALSE, recursive = TRUE)

save_oracle <- function(plots, prefix) {
  for (nm in names(plots)) {
    ggsave(paste0("fig/", prefix, ".ism.", nm, ".pdf"),
           plots[[nm]], height = 3, width = 4, dpi = 320)
  }
}

save_oracle(list(mean.3x20 = p.oracle.y.3x20.mean, max.3x20 = p.oracle.y.3x20.max,
                 mean.1x60 = p.oracle.y.1x60.mean, max.1x60 = p.oracle.y.1x60.max),
            "yeast")
save_oracle(list(mean.3x20 = p.oracle.h.3x20.mean, max.3x20 = p.oracle.h.3x20.max,
                 mean.1x60 = p.oracle.h.1x60.mean, max.1x60 = p.oracle.h.1x60.max),
            "human")
```


# Full-pool ISM {.tabset}

ISM scores for all AL pool sequences, with top-n-bottom annotations.
Mann–Whitney U compares each group to the AL pool.

## Data loading

```{r pool-load}
pool.ism.yeast <- load_pool_ism(
  path.ism  = "yeast/sequence_analysis/yeast_ism_all.tsv",
  path.pool = "yeast/sequences/yeast_pool.trim.tsv",
  path.tnb  = "yeast/sequences/yeast_top_n_bottom_1.trim.tsv",
  n.rnd     = 60000,
  trim.seq  = c(58, -14)
)

pool.ism.human <- load_pool_ism(
  path.ism  = "human/sequence_analysis/human_ism_all.tsv",
  path.pool = "human/sequences/pool.tsv",
  path.tnb  = "human/sequences/human_top_n_bottom_10.tsv",
  n.rnd     = 20000
)

cat("=== Yeast ===\n")
summary(pool.ism.yeast$mean); summary(pool.ism.yeast$max)

cat("\n=== Human ===\n")
summary(pool.ism.human$mean); summary(pool.ism.human$max)
```

## Mann–Whitney tests

```{r pool-mwt}
pool.mwt.yeast <- pool_mw_test(pool.ism.yeast, y.mean = 1.35, y.max = 8.7)
pool.mwt.human <- pool_mw_test(pool.ism.human, y.mean = 0.5,  y.max = 4.1)
```

## Plots

```{r pool-plots, fig.width=4, fig.height=4}
# Yeast
p.pool.y.3x20.mean <- plot_pool_ism(pool.ism.yeast, pool.mwt.yeast, "3x20", "mean", "Yeast — 3x20")
p.pool.y.3x20.max  <- plot_pool_ism(pool.ism.yeast, pool.mwt.yeast, "3x20", "max",  "Yeast — 3x20")
p.pool.y.1x60.mean <- plot_pool_ism(pool.ism.yeast, pool.mwt.yeast, "1x60", "mean", "Yeast — 1x60")
p.pool.y.1x60.max  <- plot_pool_ism(pool.ism.yeast, pool.mwt.yeast, "1x60", "max",  "Yeast — 1x60")

p.pool.y.3x20.mean; p.pool.y.3x20.max
p.pool.y.1x60.mean; p.pool.y.1x60.max

# Human
p.pool.h.3x20.mean <- plot_pool_ism(pool.ism.human, pool.mwt.human, "3x20", "mean", "Human — 3x20")
p.pool.h.3x20.max  <- plot_pool_ism(pool.ism.human, pool.mwt.human, "3x20", "max",  "Human — 3x20")
p.pool.h.1x60.mean <- plot_pool_ism(pool.ism.human, pool.mwt.human, "1x60", "mean", "Human — 1x60")
p.pool.h.1x60.max  <- plot_pool_ism(pool.ism.human, pool.mwt.human, "1x60", "max",  "Human — 1x60")

p.pool.h.3x20.mean; p.pool.h.3x20.max
p.pool.h.1x60.mean; p.pool.h.1x60.max
```

```{r pool-save}
save_pool <- function(plots, prefix) {
  for (nm in names(plots)) {
    ggsave(paste0("fig/", prefix, ".ism.all.", nm, ".pdf"),
           plots[[nm]], width = 4, height = 4, dpi = 320)
  }
}

save_pool(list(mean.3x20 = p.pool.y.3x20.mean, max.3x20 = p.pool.y.3x20.max,
               mean.1x60 = p.pool.y.1x60.mean, max.1x60 = p.pool.y.1x60.max),
          "yeast")
save_pool(list(mean.3x20 = p.pool.h.3x20.mean, max.3x20 = p.pool.h.3x20.max,
               mean.1x60 = p.pool.h.1x60.mean, max.1x60 = p.pool.h.1x60.max),
          "human")
```


```{r session-info}
sessionInfo()
```