expression.Rmd

---
title: "Regulatory Activity of Selected Sequences — 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)
```

Comparison of regulatory activity (expression / log2 fold-change) between
control sequences (top-n-bottom for diversity and uncertainty) and the AL
pool.  A random subsample of 100 000 AL pool sequences provides a visual
baseline.

Expression data sources:

- **Yeast**: GPRA data from de Boer et al.
- **Human**: MPRA data from Gosai et al. (K562 only)


# Shared utility functions

```{r shared-functions}
#' Read a 3-column sequence TSV and return IDs + method.
load.sequences <- function(path) {
  seq.df <- read.table(path, header = FALSE, sep = "\t",
                        col.names = c("IDs", "sequence", "method"))
  seq.df$sequence <- NULL
  seq.df
}


#' Rename raw control method labels to human-readable names.
#'
#' Uses regex so the same function works for both yeast (suffix _bottom1 /
#' _top101) and human (suffix _bottom10 / _top1010).
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
}


#' Mann–Whitney U test of each group vs the AL pool reference.
#'
#' @param data     Data.frame with columns `al.method` (factor) and `value`.
#' @param ref.label  Label of the reference group (default "AL pool").
#' @return Data.frame with method, formatted p-value, and test statistic.
mann_whitney_vs_pool <- function(data, ref.label = "AL pool") {
  ref.vals <- subset(data, al.method == ref.label)$value

  map_dfr(levels(data$al.method), function(lvl) {
    grp <- subset(data, al.method == lvl)$value
    mwt <- wilcox.test(grp, ref.vals)
    data.frame(
      al.method = lvl,
      pvalue    = ifelse(mwt$p.value <= 2.2e-16,
                          "p<2.2e-16",
                          paste0("p=", signif(mwt$p.value, 2))),
      statistic = mwt$statistic
    )
  })
}


#' Violin + boxplot of expression by method, with Mann–Whitney p-values.
#'
#' @param data       Data.frame with `al.method` (factor) and `value`.
#' @param pval.df    Data.frame from mann_whitney_vs_pool.
#' @param ylab.text  Y-axis label.
#' @param pval.y     Y position for the p-value text.
#' @param title      Plot title.
#' @param violin.adjust  Bandwidth adjustment for the violin (default 1).
plot_expression_violin <- function(data, pval.df, ylab.text, pval.y,
                                    title = NULL, violin.adjust = 1) {
  pval.df$y.axis <- pval.y

  ggplot(data, aes(x = al.method, y = value)) +
    geom_violin(width = 1.25, outliers = FALSE, adjust = violin.adjust) +
    geom_boxplot(width = 0.1, outliers = FALSE) +
    geom_text(data = pval.df,
              aes(x = al.method, y = y.axis, label = pvalue), size = 3) +
    xlab("") + ylab(ylab.text) +
    { if (!is.null(title)) 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"),
          axis.text = element_text(size = 7),
          axis.title = element_text(size = 10))
}


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


# Yeast {.tabset}

## Data loading

```{r yeast-load}
# Expression data: headerless 2-column file (sequence, expression)
yeast.expr <- read.table("yeast/reference/GPRA_Rafi2024.txt", header = FALSE,
                          sep = "\t", col.names = c("Seq", "expression"))
yeast.expr$IDs <- seq(0, nrow(yeast.expr) - 1)
yeast.expr$Seq <- NULL

# Sequences
yeast.alpool <- load.sequences("yeast/sequences/yeast_pool.trim.tsv")

yeast.alpool$al.method <- "AL pool"

yeast.tnb <- load.sequences("yeast/sequences/yeast_top_n_bottom_1.trim.tsv")
yeast.tnb$al.method <- rename_tnb_labels(yeast.tnb$method)

# Split by pool size
yeast.tnb$rounds <- ifelse(grepl("1x60", yeast.tnb$method), "1x60k", "3x20k")
yeast.tnb.3x20   <- subset(yeast.tnb, rounds == "3x20k")
yeast.tnb.1x60   <- subset(yeast.tnb, rounds == "1x60k")
```

## Merge with expression & assemble

```{r yeast-merge}
yeast.tnb.3x20.expr <- merge(yeast.tnb.3x20, yeast.expr, by = "IDs")
yeast.tnb.1x60.expr <- merge(yeast.tnb.1x60, yeast.expr, by = "IDs")
yeast.alpool.expr   <- merge(yeast.alpool,    yeast.expr, by = "IDs")

set.seed(42)
yeast.rnd.expr <- yeast.alpool.expr[sample(nrow(yeast.alpool.expr), 100000), ]
yeast.rnd.expr$al.method <- "Random selection"

# Rename value column for the shared plotting function
yeast.data.3x20 <- rbind(yeast.tnb.3x20.expr[,-4], yeast.rnd.expr, yeast.alpool.expr)
yeast.data.3x20$value <- yeast.data.3x20$expression
yeast.data.3x20$al.method <- factor(yeast.data.3x20$al.method, levels = LEVEL_ORDER)

yeast.data.1x60 <- rbind(yeast.tnb.1x60.expr[,-4], yeast.rnd.expr, yeast.alpool.expr)
yeast.data.1x60$value <- yeast.data.1x60$expression
yeast.data.1x60$al.method <- factor(yeast.data.1x60$al.method, levels = LEVEL_ORDER)
```

## Plots

```{r yeast-plots, fig.width=4, fig.height=4}
yeast.pval.3x20 <- mann_whitney_vs_pool(yeast.data.3x20)
yeast.pval.1x60 <- mann_whitney_vs_pool(yeast.data.1x60)

p.yeast.3x20 <- plot_expression_violin(
  yeast.data.3x20, yeast.pval.3x20,
  ylab.text = "Expression", pval.y = 18, violin.adjust = 3
)
p.yeast.1x60 <- plot_expression_violin(
  yeast.data.1x60, yeast.pval.1x60,
  ylab.text = "Expression", pval.y = 18, violin.adjust = 3
)

p.yeast.3x20
p.yeast.1x60
```

```{r yeast-save}
dir.create("fig", showWarnings = FALSE, recursive = TRUE)
ggsave("fig/yeast.boxplot.expression.3x20.pdf", p.yeast.3x20,
       width = 4, height = 4, dpi = 320)
ggsave("fig/yeast.boxplot.expression.1x60.pdf", p.yeast.1x60,
       width = 4, height = 4, dpi = 320)
```


# Human {.tabset}

## Data loading

```{r human-load}
# MPRA data from Gosai et al. — keep IDs and K562 columns only
human.mpra <- read.table("human/reference/41586_2024_8070_MOESM4_ESM.txt",
                          header = TRUE, sep = "\t")
human.mpra <- human.mpra[, c("IDs", "K562_log2FC")]  # IDs + K562-related columns

# Sequences
human.alpool <- load.sequences("human/sequences/pool.tsv")
human.alpool$al.method <- "AL pool"

human.tnb <- load.sequences("human/sequences/human_top_n_bottom_10.tsv")
human.tnb$al.method <- rename_tnb_labels(human.tnb$method)

# Split by pool size
human.tnb.3x20 <- human.tnb[grep("1x60", human.tnb$method, invert = TRUE), ]
human.tnb.1x60 <- human.tnb[grep("1x60", human.tnb$method), ]
```

## Merge with expression & assemble

```{r human-merge}
human.tnb.3x20.mpra <- merge(human.tnb.3x20, human.mpra, by = "IDs")
human.tnb.1x60.mpra <- merge(human.tnb.1x60, human.mpra, by = "IDs")
human.alpool.mpra   <- merge(human.alpool,    human.mpra, by = "IDs")

set.seed(42)
human.rnd.mpra <- human.alpool.mpra[sample(nrow(human.alpool.mpra), 100000), ]
human.rnd.mpra$al.method <- "Random selection"

# Rename value column for the shared plotting function
human.data.3x20 <- rbind(human.tnb.3x20.mpra, human.rnd.mpra, human.alpool.mpra)
human.data.3x20$value <- human.data.3x20$K562_log2FC
human.data.3x20$al.method <- factor(human.data.3x20$al.method, levels = LEVEL_ORDER)

human.data.1x60 <- rbind(human.tnb.1x60.mpra, human.rnd.mpra, human.alpool.mpra)
human.data.1x60$value <- human.data.1x60$K562_log2FC
human.data.1x60$al.method <- factor(human.data.1x60$al.method, levels = LEVEL_ORDER)
```

## Plots

```{r human-plots, fig.width=4, fig.height=4}
human.pval.3x20 <- mann_whitney_vs_pool(human.data.3x20)
human.pval.1x60 <- mann_whitney_vs_pool(human.data.1x60)

p.human.3x20 <- plot_expression_violin(
  human.data.3x20, human.pval.3x20,
  ylab.text = "log2FC(K562)", pval.y = 10
)
p.human.1x60 <- plot_expression_violin(
  human.data.1x60, human.pval.1x60,
  ylab.text = "log2FC(K562)", pval.y = 10
)

p.human.3x20
p.human.1x60
```

```{r human-save}
ggsave("fig/human.boxplot.expression.log2FCK562.3x20.pdf", p.human.3x20,
       width = 4, height = 4, dpi = 320)
ggsave("fig/human.boxplot.expression.log2FCK562.1x60.pdf", p.human.1x60,
       width = 4, height = 4, dpi = 320)
```


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