TFBSContent.Rmd

---
title: "TF Motif Enrichment Analysis — 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(dplyr)
library(ggplot2)
library(reshape2)
library(ggpubr)
library(ggsci)
library(gridExtra)
library(effectsize)

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)
```

Analysis of TF binding site content in sequences selected by active learning. Motif
scanning was performed with FIMO (MEME suite) using organism-specific databases and
thresholds:

- **Yeast**: JASPAR Fungi database, p-value threshold 0.0001
- **Human**: HOCOMOCO v13 database, FDR threshold 0.01

Custom nucleotide-frequency background files derived from each AL pool were supplied
to FIMO.  FIMO scanning is restricted to the AL pool sequences for filtering
consistency.

Enrichment is quantified with Cliff's delta (and the corresponding AUROC derived from
the Mann–Whitney U statistic) comparing TF motif counts between selected/control
sequences and the full AL pool.  A random-sampling baseline provides the null
expectation.

**Code to analyze and generate Fig5. E, F; Supplementary Fig. 13 E, F; Supplementary Fig. 14 A-D**

# Shared utility functions

```{r shared-functions}
#' Merge FIMO hits with a sequence table, keeping sequences with zero hits.
#'
#' @param data.tf Either a data.frame of pre-loaded FIMO hits (columns:
#'   sequence_name, motif_id) or a file path to raw FIMO output.
#' @param path.seq Path to a 3-column TSV (sequence_name, sequence, method).
#' @param tf.whitelist Character vector of TF IDs to keep when reading from
#'   file.  Ignored when data.tf is a data.frame (assumed already filtered).
#' @param motif_col Column name used for the TF identifier in the raw FIMO
#'   output ("motif_id" for HOCOMOCO, "motif_alt_id" for JASPAR).
#' @param split_motif_id Logical; if TRUE, split motif_id on "." and keep the
#'   first element (needed for HOCOMOCO IDs like "P53_HUMAN.H13MOA.0.A").
preprocess.tf <- function(data.tf, path.seq, tf.whitelist = NULL,
                          motif_col = "motif_id", split_motif_id = FALSE) {

  seq.df <- read.table(path.seq, header = FALSE, sep = "\t",
                        col.names = c("sequence_name", "sequence", "method"))
  seq.df$sequence <- NULL

  if (is.data.frame(data.tf)) {
    # data.tf already has columns: sequence_name, motif_id
    merged <- merge(seq.df, data.tf[, c("sequence_name", "motif_id")],
                    by = "sequence_name", all.x = TRUE)
    merged$motif_id <- ifelse(is.na(merged$motif_id), 0, merged$motif_id)
    return(merged)
  }

  # --- Read from file path ---
  tf.motifs <- read.table(data.tf, header = TRUE, sep = "\t")

  # Keep only sequence_name and the motif identifier column
  tf.motifs <- tf.motifs[, c("sequence_name", motif_col)]
  colnames(tf.motifs)[colnames(tf.motifs) == motif_col] <- "motif_id"

  if (split_motif_id) {
    tf.motifs$motif_id <- str_split_fixed(tf.motifs$motif_id, "\\.", 2)[, 1]
  }

  # Parse method and clean sequence_name (pipe-delimited: seqid|method)
  tf.motifs$method <- str_split_fixed(tf.motifs$sequence_name, "\\|", 2)[, 2]
  tf.motifs$sequence_name <- str_split_fixed(tf.motifs$sequence_name, "\\|", 2)[, 1]

  if (!is.null(tf.whitelist)) {
    tf.motifs <- subset(tf.motifs, motif_id %in% tf.whitelist)
  }

  merged <- merge(seq.df, tf.motifs[, c("sequence_name", "motif_id")],
                  by = "sequence_name", all.x = TRUE)
  merged$motif_id <- ifelse(is.na(merged$motif_id), 0, merged$motif_id)
  return(merged)
}


#' Count the total number of TF motif hits per sequence.
counts.tf <- function(tf.df) {
  counts.per.tf <- subset(tf.df, motif_id != 0) %>%
    group_by(sequence_name, method, motif_id) %>%
    summarise(counts = n(), .groups = "drop")

  zero.counts <- subset(tf.df, motif_id == 0)
  zero.counts$counts <- 0

  counts.per.tf <- rbind(counts.per.tf, zero.counts)

  tf.per.seq <- counts.per.tf %>%
    group_by(sequence_name, method) %>%
    summarise(counts = sum(counts), .groups = "drop")

  return(tf.per.seq)
}


#' Draw nb.sample random subsets of size nb.rnd from a counts data.frame.
random.sampling <- function(df, nb.rnd, nb.sample) {
  lapply(seq_len(nb.sample), function(i) {
    idx <- sample(nrow(df), nb.rnd)
    out <- df[idx, ]
    out$method <- paste0("RandomSampling_", i)
    out
  })
}


#' Compute Cliff's delta for each element of a list vs a reference.
compute.cliff.delta <- function(counts.list, counts.ref) {
  cd <- lapply(counts.list, function(x) {
    cliffs_delta(x$counts, counts.ref$counts)
  })
  cd <- do.call("rbind", cd)
  cd$method <- vapply(counts.list, function(x) unique(x$method), character(1))
  cd
}


#' Summarise random-baseline Cliff's delta into a single row.
summarise.rnd.cliff <- function(cliff.d.rnd.raw) {
  data.frame(
    r_rank_biserial = mean(cliff.d.rnd.raw$r_rank_biserial),
    CI              = 0.95,
    CI_low          = min(cliff.d.rnd.raw$r_rank_biserial),
    CI_high         = max(cliff.d.rnd.raw$r_rank_biserial),
    al.method       = "Random selection"
  )
}


#' Parse model architecture and AL method from the combined method label.
annotate.cliff.delta <- function(cd) {
  cd$model.arch <- toupper(str_extract(cd$method, "(attn|rnn|cnn)"))
  cd$al.method  <- str_remove(cd$method, "_(attn|rnn|cnn)")
  cd$model.arch <- factor(cd$model.arch, levels = c("RNN", "CNN", "ATTN"))
  cd
}


#' Split a Cliff's delta table into 3x20 and 1x60 subsets.
split.by.pool <- function(cd) {
  cd.60k  <- cd[grep("60k", cd$method), ]
  cd.60k$al.method <- gsub("_60k", "", cd.60k$al.method)
  cd.3x20 <- cd[grep("60k", cd$method, invert = TRUE), ]
  cd.3x20$al.method <- gsub("_3x20", "", cd.3x20$al.method)
  list(x3x20 = cd.3x20, x1x60 = cd.60k)
}


#' Rename control-sequence labels for plotting.
rename.controls <- function(cd, suffix_pattern) {
  cd$al.method <- gsub(paste0("diversity_", suffix_pattern, "_bottom\\d+"),
                        "Rarely diverse", cd$al.method)
  cd$al.method <- gsub(paste0("diversity_", suffix_pattern, "_top\\d+"),
                        "Often diverse", cd$al.method)
  cd$al.method <- gsub(paste0("uncertainty_", suffix_pattern, "_bottom\\d+"),
                        "Rarely uncertain", cd$al.method)
  cd$al.method <- gsub(paste0("uncertainty_", suffix_pattern, "_top\\d+"),
                        "Often uncertain", cd$al.method)
  cd$al.method <- factor(cd$al.method,
                          levels = c("Rarely uncertain", "Often uncertain",
                                     "Rarely diverse", "Often diverse"))
  cd
}


#' Cliff's delta dot-plot with CI and a random-baseline reference line.
plot.cliff.delta <- function(cd, rnd.ref, title,
                              ylimits = c(-0.4, 0.4), ystep = 0.1,
                              facet_by_arch = TRUE) {
  p <- ggplot(cd, aes(al.method, y = r_rank_biserial)) +
    geom_point(size = 2) +
    geom_linerange(aes(ymin = CI_low, ymax = CI_high),
                   width = 0.5, colour = "black") +
    geom_hline(yintercept = rnd.ref, linetype = 2) +
    scale_y_continuous(limits = ylimits,
                       breaks = round(seq(ylimits[1], ylimits[2], by = ystep), 2),
                       expand = expansion(mult = 0.05)) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
          axis.text = element_text(size = 7),
          plot.margin = margin(5, 5, 5, 5, "mm")) +
    labs(x = "", y = "Cliff's delta", title = title)

  if (facet_by_arch) p <- p + facet_wrap(~model.arch, scales = "free_x")
  if (!facet_by_arch) p <- p + theme(aspect.ratio = 1 / 0.75)
  p
}


#' AUROC dot-plot with a random-baseline reference line.
plot.auroc <- function(auroc.df, rnd.auroc, title, y_col,
                        colour_col = NULL, facet_by_arch = FALSE) {
  mapping <- if (!is.null(colour_col)) {
    aes(al.method, y = .data[[y_col]], colour = .data[[colour_col]])
  } else {
    aes(al.method, y = .data[[y_col]])
  }

  p <- ggplot(auroc.df, mapping) +
    geom_point(size = 2) +
    geom_hline(yintercept = rnd.auroc, linetype = 2) +
    scale_y_continuous(limits = c(0, 1),
                       breaks = round(seq(0, 1, by = 0.1), 2),
                       expand = expansion(mult = 0.05)) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
          axis.text = element_text(size = 7),
          plot.margin = margin(5, 5, 5, 5, "mm")) +
    labs(x = "", y = "TFBS content enrichment (AUROC)", title = title)

  if (facet_by_arch) p <- p + facet_wrap(~model.arch, scales = "free_x")
  if (!facet_by_arch) p <- p + theme(aspect.ratio = 1 / 0.75, legend.position = "none")
  p
}
```


# Yeast — data loading & TF filtering {.tabset}

Yeast TF identifiers are obtained from the sacCer3 gene annotation; only TFs
with gene symbols present in the annotation are retained.

```{r yeast-tf-filtering}
yeast.gtf <- read.table("yeast/reference/sacCer3.ncbiRefSeq.gtf.gz",
                          header = FALSE, sep = "\t")
yeast.genes <- unique(gsub(";", "",
  str_split_fixed(yeast.gtf$V9, " ", Inf)[, ncol(str_split_fixed(yeast.gtf$V9, " ", Inf))]
))
# Simpler extraction matching the original logic:
yeast.genes <- unique(gsub(";", "",
  vapply(str_split(yeast.gtf$V9, " "), function(x) x[length(x)], character(1))
))
```

FIMO was run per-chromosome; read all result files and concatenate.

```{r yeast-fimo-import}
fimo.files <- list.files("yeast/sequence_analysis/FIMO_ALpool_jaspar/",
                          pattern = "fimo.tsv", recursive = TRUE,
                          full.names = TRUE)

fimo.al.yeast <- lapply(fimo.files, function(x) {
  df <- read.table(x, header = TRUE, sep = "\t")
  df <- subset(df, motif_alt_id %in% yeast.genes)
  df <- df[, c("sequence_name", "motif_alt_id")]
  colnames(df)[2] <- "motif_id"
  df$method <- str_split_fixed(df$sequence_name, "\\|", 2)[, 2]
  df$sequence_name <- str_split_fixed(df$sequence_name, "\\|", 2)[, 1]
  df
})
fimo.al.yeast <- do.call("rbind", fimo.al.yeast)
```

```{r yeast-preprocess}
tf.motifs.al.yeast       <- preprocess.tf(fimo.al.yeast,
                                           "yeast/sequences/yeast_pool.trim.tsv")
tf.motifs.tnb.yeast      <- preprocess.tf(fimo.al.yeast,
                                           "yeast/sequences/yeast_top_n_bottom_1.trim.tsv")
tf.motifs.tnb.yeast$method <- gsub("1x60", "60k", tf.motifs.tnb.yeast$method)

tf.motifs.selected.yeast <- preprocess.tf(fimo.al.yeast,
                                           "yeast/sequences/master_output.tsv")
tf.motifs.selected.yeast$method <- str_remove(
  tf.motifs.selected.yeast$method, "(_[1-5])"
)

# Remove duplicated architecture pairs (keep RNN copy only)
duplicated.methods <- c("rnn-cnn_attn", "diff_arch_attn",
                        "rnn-cnn_60k_attn", "diff_arch_60k_attn",
                        "rnn-cnn_cnn", "diff_arch_cnn",
                        "rnn-cnn_60k_cnn", "diff_arch_60k_cnn")
tf.motifs.selected.yeast <- subset(tf.motifs.selected.yeast,
                                    !(method %in% duplicated.methods))
```


# Human — data loading & TF filtering {.tabset}

Human TF identifiers are obtained by cross-referencing GTEx Whole Blood
expression (TPM >= 1) with the TF masterlist from Lambert et al. to map gene
symbols to UniProt IDs used in HOCOMOCO.  Filtering for K562-relevant TFs.

```{r human-tf-filtering}
gtex <- read.table(
  "human/reference/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct.gz",
  header = TRUE, sep = "\t", skip = 2
)

tf.masterlist <- read.table("human/reference/tf_masterlist.tsv",
                             header = TRUE, sep = "\t", fill = TRUE)
tf.masterlist <- tf.masterlist[grep("HUMAN", tf.masterlist$curated.uniprot_id), ]
tf.masterlist$curated.uniprot_id <- str_split_fixed(
  tf.masterlist$curated.uniprot_id, "_", 2
)[, 1]
tf.masterlist <- tf.masterlist[, c("curated.uniprot_id", "auto.gene_symbol")]

gtex.merged <- merge(gtex[, c("Description", "Whole.Blood")],
                      tf.masterlist,
                      by.x = "Description", by.y = "auto.gene_symbol")
human.selected.TF <- subset(gtex.merged, Whole.Blood >= 1)$curated.uniprot_id
```

```{r human-preprocess}
tf.motifs.al.human <- preprocess.tf(
  "human/sequence_analysis/hocomoco_ALpool_FDR0.01/fimo.tsv",
  "human/sequences/pool.tsv",
  tf.whitelist    = human.selected.TF,
  motif_col       = "motif_id",
  split_motif_id  = TRUE
)

tf.motifs.selected.human <- preprocess.tf(tf.motifs.al.human,
                                           "human/sequences/master_output.tsv")
tf.motifs.selected.human$method <- str_remove(
  tf.motifs.selected.human$method, "(_[1-5])"
)
tf.motifs.selected.human$method <- gsub("(_1x60)", "",
                                         tf.motifs.selected.human$method)

tf.motifs.tnb.human <- preprocess.tf(
  tf.motifs.al.human,
  "human/sequences/human_top_n_bottom_10.tsv"
)
tf.motifs.tnb.human$method <- gsub("1x60", "60k", tf.motifs.tnb.human$method)

tf.motifs.selected.human <- subset(tf.motifs.selected.human,
                                    !(method %in% duplicated.methods))
```


# TF motif counts per sequence

```{r motif-counts}
# Yeast
counts.selected.yeast <- counts.tf(tf.motifs.selected.yeast)
counts.tnb.yeast      <- counts.tf(tf.motifs.tnb.yeast)
counts.al.yeast       <- counts.tf(tf.motifs.al.yeast)

# Human
counts.selected.human <- counts.tf(tf.motifs.selected.human)
counts.tnb.human      <- counts.tf(tf.motifs.tnb.human)
counts.al.human       <- counts.tf(tf.motifs.al.human)
```


# Cliff's delta {.tabset}

Cliff's delta measures the degree to which TF motif counts in a given set of
sequences tend to be higher (or lower) than in the full AL pool.

## Random baseline

```{r random-baseline}
set.seed(42)

# Yeast — 60 000 sequences x 15 replicates
counts.rnd.yeast <- random.sampling(counts.al.yeast, 60000, 15)
cliff.d.rnd.raw.yeast <- compute.cliff.delta(counts.rnd.yeast, counts.al.yeast)
cliff.d.rnd.yeast <- summarise.rnd.cliff(cliff.d.rnd.raw.yeast)

cat("Yeast random baseline Cliff's delta:\n")
summary(cliff.d.rnd.raw.yeast$r_rank_biserial)
cat("SD:", sd(cliff.d.rnd.raw.yeast$r_rank_biserial), "\n")

# Human — 20 000 sequences x 15 replicates
counts.rnd.human <- random.sampling(counts.al.human, 20000, 15)
cliff.d.rnd.raw.human <- compute.cliff.delta(counts.rnd.human, counts.al.human)
cliff.d.rnd.human <- summarise.rnd.cliff(cliff.d.rnd.raw.human)

cat("\nHuman random baseline Cliff's delta:\n")
summary(cliff.d.rnd.raw.human$r_rank_biserial)
cat("SD:", sd(cliff.d.rnd.raw.human$r_rank_biserial), "\n")
```

## Selected sequences

```{r cliff-delta-selected}
# --- Yeast ---
counts.selected.yeast.list <- counts.selected.yeast %>%
  group_by(method) %>%
  group_split()

cd.selected.yeast <- compute.cliff.delta(counts.selected.yeast.list,
                                          counts.al.yeast)
cd.selected.yeast <- annotate.cliff.delta(cd.selected.yeast)
cd.sel.yeast      <- split.by.pool(cd.selected.yeast)

# Align 1x60 AL method labels with 3x20 (same underlying methods)
cd.sel.yeast$x1x60$al.method <- cd.sel.yeast$x3x20$al.method

cat("Yeast selected — Cliff's delta summary:\n")
summary(cd.selected.yeast$r_rank_biserial)

# --- Human ---
counts.selected.human.list <- counts.selected.human %>%
  ungroup() %>%
  group_by(method) %>%
  group_split()

cd.selected.human <- compute.cliff.delta(counts.selected.human.list,
                                          counts.al.human)
cd.selected.human <- annotate.cliff.delta(cd.selected.human)
cd.sel.human      <- split.by.pool(cd.selected.human)

cat("\nHuman selected — Cliff's delta summary:\n")
summary(cd.selected.human$r_rank_biserial)
```

## Control sequences (top-n-bottom)

```{r cliff-delta-controls}
process.controls <- function(counts.al, counts.tnb, cliff.d.rnd) {
  counts.ctrl <- rbind(counts.al, counts.tnb) %>%
    group_by(method) %>%
    group_split()

  cd.ctrl <- compute.cliff.delta(counts.ctrl, counts.al)
  cd.ctrl$al.method <- vapply(counts.ctrl, function(x) unique(x$method), character(1))
  cd.ctrl$al.method <- gsub("1x60", "60k", cd.ctrl$al.method)

  # Keep only the columns shared with the random-baseline summary before binding
  shared.cols <- c("r_rank_biserial", "CI", "CI_low", "CI_high", "al.method")
  cd.ctrl <- rbind(cd.ctrl[, shared.cols], cliff.d.rnd[, shared.cols])

  cd.3x20 <- cd.ctrl[grep("(60k)|(ALpool)|(Random selection)",
                            cd.ctrl$al.method, invert = TRUE), ]
  cd.3x20 <- rename.controls(cd.3x20, "3x20")

  cd.60k <- cd.ctrl[grep("60k", cd.ctrl$al.method), ]
  cd.60k <- rename.controls(cd.60k, "60k")

  list(x3x20 = cd.3x20, x1x60 = cd.60k)
}

cd.ctrl.yeast <- process.controls(counts.al.yeast, counts.tnb.yeast,
                                   cliff.d.rnd.yeast)
cd.ctrl.human <- process.controls(counts.al.human, counts.tnb.human,
                                   cliff.d.rnd.human)
```


# Plots — Cliff's delta {.tabset}

## Selected sequences

```{r plots-cliff-selected, fig.width=6, fig.height=3}
# Yeast
p.cd.sel.yeast.3x20 <- plot.cliff.delta(
  cd.sel.yeast$x3x20,
  rnd.ref  = c(cliff.d.rnd.yeast$CI_low, cliff.d.rnd.yeast$CI_high),
  title    = "Yeast — 3x20",
  ylimits  = c(-0.6, 0.6), ystep = 0.2
)
p.cd.sel.yeast.1x60 <- plot.cliff.delta(
  cd.sel.yeast$x1x60,
  rnd.ref  = c(cliff.d.rnd.yeast$CI_low, cliff.d.rnd.yeast$CI_high),
  title    = "Yeast — 1x60",
  ylimits  = c(-0.6, 0.6), ystep = 0.2
)
p.cd.sel.yeast.3x20
p.cd.sel.yeast.1x60
```

```{r plots-cliff-selected-human, fig.width=6, fig.height=3}
# Human
p.cd.sel.human.3x20 <- plot.cliff.delta(
  cd.sel.human$x3x20,
  rnd.ref  = cliff.d.rnd.human$r_rank_biserial,
  title    = "Human — 3x20",
  ylimits  = c(-0.3, 0.3), ystep = 0.1
)
p.cd.sel.human.1x60 <- plot.cliff.delta(
  cd.sel.human$x1x60,
  rnd.ref  = cliff.d.rnd.human$r_rank_biserial,
  title    = "Human — 1x60",
  ylimits  = c(-0.3, 0.3), ystep = 0.1
)
p.cd.sel.human.3x20
p.cd.sel.human.1x60
```

## Control sequences

```{r plots-cliff-controls, fig.width=5, fig.height=3.5}
# Yeast
p.cd.ctrl.yeast.3x20 <- plot.cliff.delta(
  cd.ctrl.yeast$x3x20,
  rnd.ref  = c(cliff.d.rnd.yeast$CI_low, cliff.d.rnd.yeast$CI_high),
  title    = "Yeast — 3x20",
  ylimits  = c(-0.2, 0.2), ystep = 0.1, facet_by_arch = FALSE
)
p.cd.ctrl.yeast.1x60 <- plot.cliff.delta(
  cd.ctrl.yeast$x1x60,
  rnd.ref  = c(cliff.d.rnd.yeast$CI_low, cliff.d.rnd.yeast$CI_high),
  title    = "Yeast — 1x60",
  ylimits  = c(-0.2, 0.2), ystep = 0.1, facet_by_arch = FALSE
)
p.cd.ctrl.yeast.3x20
p.cd.ctrl.yeast.1x60

# Human
p.cd.ctrl.human.3x20 <- plot.cliff.delta(
  cd.ctrl.human$x3x20,
  rnd.ref  = cliff.d.rnd.human$r_rank_biserial,
  title    = "Human — 3x20",
  ylimits  = c(-0.4, 0.4), ystep = 0.1, facet_by_arch = FALSE
)
p.cd.ctrl.human.1x60 <- plot.cliff.delta(
  cd.ctrl.human$x1x60,
  rnd.ref  = cliff.d.rnd.human$r_rank_biserial,
  title    = "Human — 1x60",
  ylimits  = c(-0.4, 0.4), ystep = 0.1, facet_by_arch = FALSE
)
p.cd.ctrl.human.3x20
p.cd.ctrl.human.1x60
```


# Save Cliff's delta plots

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

# Selected
ggsave("fig/yeast.tfmotifs.selected.cliff.d.3x20.pdf",
       plot = p.cd.sel.yeast.3x20, height = 3, width = 3, dpi = 320)
ggsave("fig/yeast.tfmotifs.selected.cliff.d.60k.pdf",
       plot = p.cd.sel.yeast.1x60, height = 3, width = 3, dpi = 320)
ggsave("fig/human.tfmotifs.selected.cliff.d.3x20.pdf",
       plot = p.cd.sel.human.3x20, height = 3, width = 3, dpi = 320)
ggsave("fig/human.tfmotifs.selected.cliff.d.60k.pdf",
       plot = p.cd.sel.human.1x60, height = 3, width = 3, dpi = 320)

# Controls
ggsave("fig/yeast.tfmotifs.controls.cliff.d.3x20.pdf",
       plot = p.cd.ctrl.yeast.3x20, height = 3.5, width = 2.5, dpi = 320)
ggsave("fig/yeast.tfmotifs.controls.cliff.d.1x60.pdf",
       plot = p.cd.ctrl.yeast.1x60, height = 3.5, width = 2.5, dpi = 320)
ggsave("fig/human.tfmotifs.controls.cliff.d.3x20.pdf",
       plot = p.cd.ctrl.human.3x20, height = 3.5, width = 2.5, dpi = 320)
ggsave("fig/human.tfmotifs.controls.cliff.d.1x60.pdf",
       plot = p.cd.ctrl.human.1x60, height = 3.5, width = 2.5, dpi = 320)
```


# TFBS enrichment — AUROC {.tabset}

AUROC is derived from the Mann–Whitney U statistic via the relationship
AUROC = (Cliff's delta + 1) / 2.

```{r auroc-compute}
build.auroc.controls <- function(cd.ctrl.3x20, cd.ctrl.60k, cd.rnd) {
  auroc.ctrl <- data.frame(
    al.method  = cd.ctrl.3x20$al.method,
    auroc.3x20 = (cd.ctrl.3x20$r_rank_biserial + 1) / 2,
    auroc.1x60 = (cd.ctrl.60k$r_rank_biserial  + 1) / 2,
    freq       = factor(rep(c("Rarely", "Often"), 2),
                        levels = c("Rarely", "Often"))
  )
  auroc.rnd <- data.frame(
    al.method = cd.rnd$al.method,
    auroc     = (cd.rnd$r_rank_biserial + 1) / 2
  )
  list(controls = auroc.ctrl, rnd = auroc.rnd)
}

build.auroc.selected <- function(cd.sel.3x20, cd.sel.60k) {
  data.frame(
    auroc.3x20 = (cd.sel.3x20$r_rank_biserial + 1 )/ 2,
    auroc.1x60 = (cd.sel.60k$r_rank_biserial   + 1) / 2,
    method     = cd.sel.3x20$method,
    model.arch = cd.sel.3x20$model.arch,
    al.method  = cd.sel.3x20$al.method
  )
}

auroc.yeast.ctrl <- build.auroc.controls(cd.ctrl.yeast$x3x20,
                                          cd.ctrl.yeast$x1x60,
                                          cliff.d.rnd.yeast)
auroc.human.ctrl <- build.auroc.controls(cd.ctrl.human$x3x20,
                                          cd.ctrl.human$x1x60,
                                          cliff.d.rnd.human)

auroc.yeast.sel <- build.auroc.selected(cd.sel.yeast$x3x20,
                                         cd.sel.yeast$x1x60)
auroc.human.sel <- build.auroc.selected(cd.sel.human$x3x20,
                                         cd.sel.human$x1x60)
```

## Control sequences

```{r auroc-controls-plots, fig.width=5, fig.height=3.5}
# Yeast
p.auroc.ctrl.yeast.3x20 <- plot.auroc(
  auroc.yeast.ctrl$controls, auroc.yeast.ctrl$rnd$auroc,
  title = "Yeast — 3x20", y_col = "auroc.3x20", colour_col = "freq"
)
p.auroc.ctrl.yeast.1x60 <- plot.auroc(
  auroc.yeast.ctrl$controls, auroc.yeast.ctrl$rnd$auroc,
  title = "Yeast — 1x60", y_col = "auroc.1x60", colour_col = "freq"
)
p.auroc.ctrl.yeast.3x20
p.auroc.ctrl.yeast.1x60

# Human
p.auroc.ctrl.human.3x20 <- plot.auroc(
  auroc.human.ctrl$controls, auroc.human.ctrl$rnd$auroc,
  title = "Human — 3x20", y_col = "auroc.3x20", colour_col = "freq"
)
p.auroc.ctrl.human.1x60 <- plot.auroc(
  auroc.human.ctrl$controls, auroc.human.ctrl$rnd$auroc,
  title = "Human — 1x60", y_col = "auroc.1x60", colour_col = "freq"
)
p.auroc.ctrl.human.3x20
p.auroc.ctrl.human.1x60
```

## Selected sequences

```{r auroc-selected-plots, fig.width=6, fig.height=3}
# Yeast
p.auroc.sel.yeast.3x20 <- plot.auroc(
  auroc.yeast.sel, auroc.yeast.ctrl$rnd$auroc,
  title = "Yeast — 3x20", y_col = "auroc.3x20", facet_by_arch = TRUE
)
p.auroc.sel.yeast.1x60 <- plot.auroc(
  auroc.yeast.sel, auroc.yeast.ctrl$rnd$auroc,
  title = "Yeast — 1x60", y_col = "auroc.1x60", facet_by_arch = TRUE
)
p.auroc.sel.yeast.3x20
p.auroc.sel.yeast.1x60

# Human
p.auroc.sel.human.3x20 <- plot.auroc(
  auroc.human.sel, auroc.human.ctrl$rnd$auroc,
  title = "Human — 3x20", y_col = "auroc.3x20", facet_by_arch = TRUE
)
p.auroc.sel.human.1x60 <- plot.auroc(
  auroc.human.sel, auroc.human.ctrl$rnd$auroc,
  title = "Human — 1x60", y_col = "auroc.1x60", facet_by_arch = TRUE
)
p.auroc.sel.human.3x20
p.auroc.sel.human.1x60
```


# Save AUROC plots

```{r save-auroc-plots}
dir.create("fig")
# Controls
ggsave("fig/yeast.tfmotifs.controls.auroc.3x20.pdf",
       plot = p.auroc.ctrl.yeast.3x20, height = 3.5, width = 2.5, dpi = 320)
ggsave("fig/yeast.tfmotifs.controls.auroc.1x60.pdf",
       plot = p.auroc.ctrl.yeast.1x60, height = 3.5, width = 2.5, dpi = 320)
ggsave("fig/human.tfmotifs.controls.auroc.3x20.pdf",
       plot = p.auroc.ctrl.human.3x20, height = 3.5, width = 2.5, dpi = 320)
ggsave("fig/human.tfmotifs.controls.auroc.1x60.pdf",
       plot = p.auroc.ctrl.human.1x60, height = 3.5, width = 2.5, dpi = 320)

# Selected
ggsave("fig/yeast.tfmotifs.selected.auroc.3x20.pdf",
       plot = p.auroc.sel.yeast.3x20, height = 3, width = 3, dpi = 320)
ggsave("fig/yeast.tfmotifs.selected.auroc.60k.pdf",
       plot = p.auroc.sel.yeast.1x60, height = 3, width = 3, dpi = 320)
ggsave("fig/human.tfmotifs.selected.auroc.3x20.pdf",
       plot = p.auroc.sel.human.3x20, height = 3, width = 3, dpi = 320)
ggsave("fig/human.tfmotifs.selected.auroc.60k.pdf",
       plot = p.auroc.sel.human.1x60, height = 3, width = 3, dpi = 320)
```


# Export sequences with TF counts

```{r export-counts}
dir.create("files4Justin/selected", showWarnings = FALSE, recursive = TRUE)

export.seq.counts <- function(path.pool, path.tnb, path.selected,
                               counts.al, counts.tnb,
                               counts.selected, organism) {

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

  # AL pool
  al.out <- merge(pool.seq[, c("sequence_name", "sequence")],
                   counts.al, by = "sequence_name")
  al.out <- unique(al.out)

  # Top-n-bottom controls
  tnb.out <- merge(tnb.seq[, c("sequence_name", "sequence")],
                    counts.tnb, by = "sequence_name")
  tnb.out <- unique(tnb.out)

  write.table(al.out,
              paste0("files4Justin/", organism, ".alPool.TFcounts.tsv"),
              sep = "\t", col.names = TRUE, row.names = FALSE, quote = FALSE)
  write.table(tnb.out,
              paste0("files4Justin/", organism, ".TopnBottom.TFcounts.tsv"),
              sep = "\t", col.names = TRUE, row.names = FALSE, quote = FALSE)

  # Selected sequences (one file per method)
  if (!is.null(path.selected)) {
    sel.seq <- read.table(path.selected, header = FALSE, sep = "\t",
                           col.names = c("sequence_name", "sequence", "method"))
    sel.list <- counts.selected %>% ungroup() %>% group_split(method)
    for (df in sel.list) {
      out <- merge(sel.seq[, c("sequence_name", "sequence")], df,
                   by = "sequence_name")
      out <- unique(out)
      write.table(out,
                  paste0("files4Justin/selected/", organism, ".",
                         unique(df$method), ".tsv"),
                  sep = "\t", col.names = TRUE, row.names = FALSE,
                  quote = FALSE)
    }
  }
}

# Yeast
export.seq.counts(
  path.pool     = "sequences/yeast_pool.trim.tsv",
  path.tnb      = "sequences/yeast_top_n_bottom_1.trim.tsv",
  path.selected = "sequences/master_output.tsv",
  counts.al     = counts.al.yeast,
  counts.tnb    = counts.tnb.yeast,
  counts.selected = counts.selected.yeast,
  organism      = "yeast"
)

# Human
export.seq.counts(
  path.pool     = "sequences/pool.tsv",
  path.tnb      = "sequences/top_n_bottom/human_top_n_bottom_10.tsv",
  path.selected = "sequences/master_output.tsv",
  counts.al     = counts.al.human,
  counts.tnb    = counts.tnb.human,
  counts.selected = counts.selected.human,
  organism      = "human"
)
```


# Summary statistics

```{r summary-stats}
cat("=== Yeast — control sequence TF count summaries ===\n")
counts.ctrl.yeast <- rbind(counts.al.yeast, counts.tnb.yeast)
counts.ctrl.yeast %>%
  group_by(method) %>%
  group_split() %>%
  lapply(function(x) { cat(unique(x$method), "\n"); print(summary(x$counts)) })

cat("\n=== Human — control sequence TF count summaries ===\n")
counts.ctrl.human <- rbind(counts.al.human, counts.tnb.human)
counts.ctrl.human %>%
  group_by(method) %>%
  group_split() %>%
  lapply(function(x) { cat(unique(x$method), "\n"); print(summary(x$counts)) })
```

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