---
title: "Motif comparisons and P-values"
date: 5 May 2026
author:
- Benjamin Jean-Marie Tremblay^[benjamin.tremblay@uwaterloo.ca]
bibliography: universalmotif.bib
abstract: >
  Two important but not often discussed topics with regards to motifs are motif comparisons and P-values. These are explored here, including implementation details and example use cases.
vignette: >
  %\VignetteIndexEntry{Motif comparisons and P-values}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output:
  bookdown::pdf_document2
---

```{r setup, echo=FALSE}
knitr::opts_chunk$set(collapse=TRUE, comment = "#>")
suppressPackageStartupMessages(library(universalmotif))
suppressPackageStartupMessages(library(Biostrings))
suppressMessages(suppressPackageStartupMessages(library(MotifDb)))
suppressMessages(suppressPackageStartupMessages(library(ggplot2)))
suppressMessages(suppressPackageStartupMessages(library(ggtree)))
suppressMessages(suppressPackageStartupMessages(library(dplyr)))
data(ArabidopsisPromoters)
data(ArabidopsisMotif)
motdb <- convert_motifs(MotifDb)
```

# Introduction

This vignette covers motif comparisons (including metrics, parameters and clustering) and P-values. For an introduction to sequence motifs, see the [introductory](IntroductionToSequenceMotifs.pdf) vignette. For a basic overview of available motif-related functions, see the [motif manipulation](MotifManipulation.pdf) vignette. For sequence-related utilities, see the [sequences](SequenceSearches.pdf) vignette.

# Motif comparisons

There are a couple of functions available in other Bioconductor packages which allow for motif comparison, such as `PWMSimilarity()` (`TFBSTools`) and `motifSimilarity()` (`PWMEnrich`). Unfortunately these functions are not designed for comparing large numbers of motifs. Furthermore they are restrictive in their option range. The `universalmotif` package aims to fix this by providing the `compare_motifs()` / `compare_motifs2()` function. Several other functions also make use of the core `compare_motifs()` / `compare_motifs2()` functionality, including `merge_motifs()` / `merge_motifs2()` and `view_motifs()` / `view_motifs2()`. (The `2` versions of these functions are newer, faster, simplified versions.) The following introduction of motif comparison is specifically concerning the original fully-featured `compare_motifs()`

## An overview of available comparison metrics

The `compare_motifs()` function has been written to allow comparisons using any of the following metrics:

* Euclidean distance (`EUCL`)
* Weighted Euclidean distance (`WEUCL`)
* Kullback-Leibler divergence (`KL`) [@kl; @roepcke]
* Hellinger distance (`HELL`) [@hellinger]
* Squared Euclidean distance (`SEUCL`)
* Manhattan distance (`MAN`)
* Pearson correlation coefficient (`PCC`)
* Weighted Pearson correlation coefficient (`WPCC`)
* Sandelin-Wasserman similarity (`SW`; or sum of squared distances) [@wasserman]
* Average log-likelihood ratio (`ALLR`) [@wang]
* Lower limit average log-likelihood ratio (`ALLR_LL`; minimum column score of -2) [@mahony]
* Bhattacharyya coefficient (`BHAT`) [@bhatt]

For clarity, here are the `R` implementations of these metrics:

```{r}
EUCL <- function(c1, c2) {
  sqrt( sum( (c1 - c2)^2 ) )
}

WEUCL <- function(c1, c2, bkg1, bkg2) {
  sqrt( sum( (bkg1 + bkg2) * (c1 - c2)^2 ) )
}

KL <- function(c1, c2) {
  ( sum(c1 * log(c1 / c2)) + sum(c2 * log(c2 / c1)) ) / 2
}

HELL <- function(c1, c2) {
  sqrt( sum( ( sqrt(c1) - sqrt(c2) )^2 ) ) / sqrt(2)
}

SEUCL <- function(c1, c2) {
  sum( (c1 - c2)^2 )
}

MAN <- function(c1, c2) {
  sum ( abs(c1 - c2) )
}

PCC <- function(c1, c2) {
  n <- length(c1)
  top <- n * sum(c1 * c2) - sum(c1) * sum(c2)
  bot <- sqrt( ( n * sum(c1^2) - sum(c1)^2 ) * ( n * sum(c2^2) - sum(c2)^2 ) )
  top / bot
}

WPCC <- function(c1, c2, bkg1, bkg2) {
  weights <- bkg1 + bkg2
  mean1 <- sum(weights * c1)
  mean2 <- sum(weights * c2)
  var1 <- sum(weights * (c1 - mean1)^2)
  var2 <- sum(weights * (c2 - mean2)^2)
  cov <- sum(weights * (c1 - mean1) * (c2 - mean2))
  cov / sqrt(var1 * var2)
}

SW <- function(c1, c2) {
  2 - sum( (c1 - c2)^2 )
}

ALLR <- function(c1, c2, bkg1, bkg2, nsites1, nsites2) {
  left <- sum( c2 * nsites2 * log(c1 / bkg1) )
  right <- sum( c1 * nsites1 * log(c2 / bkg2) )
  ( left + right ) / ( nsites1 + nsites2 )
}

BHAT <- function(c1, c2) {
  sum( sqrt(c1 * c2) )
}
```

Motif comparison involves comparing a single column from each motif individually, and adding up the scores from all column comparisons. Since this causes the score to be highly dependent on motif length, the scores can instead be averaged using the arithmetic mean, geometric mean, median, or Fisher Z-transform.

If you're curious as to how the comparison metrics perform, two columns can be compared individually using `compare_columns()`:

```{r}
c1 <- c(0.7, 0.1, 0.1, 0.1)
c2 <- c(0.5, 0.0, 0.2, 0.3)

compare_columns(c1, c2, "PCC")
compare_columns(c1, c2, "EUCL")
```

Note that some metrics do not work with zero values, and small pseudocounts are automatically added to motifs for the following:

* `KL`
* `ALLR`
* `ALLR_LL`

As seen in figure \ref{fig:fig1}, the distributions for random individual column comparisons tend to be very skewed. This is usually remedied when comparing the entire motif, though some metrics still perform poorly in this regard.

```{r,echo=FALSE,fig.wide=TRUE,fig.asp=1,fig.cap="\\label{fig:fig1}Distributions of scores from approximately 500 random motif and individual column comparisons"}
atm <- filter_motifs(motdb, organism = "Athaliana")
pool <- do.call(cbind, atm)@motif
pool <- pool + 0.01
metrics <- universalmotif:::COMPARE_METRICS

res <- vector("list", length(metrics))
names(res) <- metrics

res2 <- vector("list", length(metrics))
names(res2) <- metrics
res3 <- vector("list", length(metrics))
names(res3) <- metrics
res4 <- vector("list", length(metrics))
names(res4) <- metrics
res4 <- res4[!names(res4) %in% c("PCC", "ALLR", "ALLR_LL")]
res5 <- vector("list", length(metrics))
names(res5) <- metrics
res9 <- vector("list", length(metrics))
names(res9) <- metrics
res8 <- vector("list", length(metrics))
names(res8) <- metrics
res8 <- res8[!names(res8) %in% c("PCC", "ALLR", "ALLR_LL")]
res99 <- vector("list", length(metrics))
names(res99) <- metrics

y1 <- atm[sample(1:length(atm), 33)]
x1 <- pool[, sample(1:ncol(pool), 528)]
x2 <- pool[, sample(1:ncol(pool), 528)]

for (m in metrics) {
  res[[m]] <- numeric(528)
  for (i in 1:528) {
    res[[m]][i] <- compare_columns(x1[, i], x2[, i], m)
  }
  res2[[m]] <- as.numeric(as.dist(suppressWarnings(compare_motifs(y1, method=m, score.strat="sum", min.mean.ic=0))))
  res3[[m]] <- as.numeric(as.dist(suppressWarnings(compare_motifs(y1, method=m, score.strat="a.mean", min.mean.ic=0))))
  if (!m %in% c("PCC", "ALLR", "ALLR_LL")) {
    res4[[m]] <- as.numeric(as.dist(suppressWarnings(compare_motifs(y1, method=m, score.strat="g.mean", min.mean.ic=0))))
    res8[[m]] <- as.numeric(as.dist(suppressWarnings(compare_motifs(y1, method=m, score.strat="wg.mean", min.mean.ic=0))))
  }
  res5[[m]] <- as.numeric(as.dist(suppressWarnings(compare_motifs(y1, method=m, score.strat="median", min.mean.ic=0))))
  res9[[m]] <- as.numeric(as.dist(suppressWarnings(compare_motifs(y1, method=m, score.strat="wa.mean", min.mean.ic=0))))
  res99[[m]] <- as.numeric(as.dist(suppressWarnings(compare_motifs(y1, method=m, score.strat="fzt", min.mean.ic=0))))
}

l_2_df <- function(x) {

  d <- data.frame(key = character(), scores = numeric())

  for (i in seq_along(x)) {
    y <- x[[i]]
    y <- y[!is.na(y)]
    d <- rbind(d, data.frame(key = rep(names(x)[i], length(y)), scores = y))
  }

  d

}

res <- l_2_df(res)
res2 <- l_2_df(res2)
res3 <- l_2_df(res3)
res4 <- l_2_df(res4)
res5 <- l_2_df(res5)
res9 <- l_2_df(res9)
res8 <- l_2_df(res8)
res99 <- l_2_df(res99)

res$type <- "RawColumnScores"
res2$type <- "Sum"
res3$type <- "ArithMean"
res4$type <- "GeoMean"
res5$type <- "Median"
res9$type <- "WeightedArithMean"
res99$type <- "FisherZTrans"
res8$type <- "WeightedGeoMean"
res6 <- rbind(res, res3, res4, res5, res9, res8, res99)

dres <- res6 %>%
          group_by(key, type) %>%
          summarise(mean = mean(scores, na.rm = TRUE), .groups = "drop")

ggplot(res6, aes(x = scores, fill = type)) +
  geom_density(alpha = 0.3, adjust = 2) +
  geom_vline(aes(xintercept = 0), colour = "black", linetype = "dashed") +
  # geom_vline(data=dres, aes(xintercept=mean, colour = type)) +
  facet_wrap(key ~ ., ncol = 3, scales = "free") +
  theme_minimal() +
  theme(text = element_text(family = "Times"))
```

## Comparison parameters

There are several key parameters to keep in mind when comparing motifs. Some of these are:

* `method`: one of the metrics listed previously
* `tryRC`: choose whether to try comparing the reverse complements of each motif as well
* `min.overlap`: limit the amount of allowed overhang between the two motifs
* `min.mean.ic`, `min.position.ic`: don't allow low IC alignments or positions to contribute to the final score
* `score.strat`: how to combine individual column scores in an alignment

See the following example for an idea as to how some of these settings impact scores:

```{r,echo=FALSE,fig.cap="\\label{fig:fig2}Example scores from comparing two motifs",fig.height=4,fig.width=5}
library(universalmotif)
library(MotifDb)

motifs <- convert_motifs(MotifDb)
motifs <- filter_motifs(motifs, altname = c("M0003_1.02", "M0004_1.02"))
# summarise_motifs(motifs)
view_motifs(motifs) +
  theme(text = element_text(family = "Times"))
try_all <- function(motifs, ...) {
  scores <- vector("list", 4)
  methods <- c("PCC", "WPCC", "EUCL", "SW", "KL", "ALLR",
               "BHAT", "HELL", "WEUCL", "SEUCL", "MAN", "ALLR_LL")
  for (i in seq_along(methods)) {
    scores[[1]][i] <- compare_motifs(motifs, method = methods[i])[1, 2]
    scores[[2]][i] <- compare_motifs(motifs, method = methods[i], normalise.scores = TRUE)[1, 2]
    scores[[3]][i] <- compare_motifs(motifs, method = methods[i], min.overlap = 99)[1, 2]
    scores[[4]][i] <- compare_motifs(motifs, method = methods[i], min.position.ic=0.25)[1, 2]
  }
  res <- data.frame(type = c("similarity", "similarity", "distance", 
                             "similarity", "distance", "similarity",
                             "similarity", "distance",
                             "distance", "distance", 
                             "distance", "similarity"),
                    method = methods,
                    default = scores[[1]],
                    normalised = scores[[2]],
                    # noRC = scores[[3]],
                    checkIC = scores[[4]])
  knitr::kable(res, format = "markdown", caption = "Comparing two motifs with various settings")
}
try_all(motifs)
```

Settings used in the previous table:

* normalised: `normalise.scores = TRUE`
* checkIC: `min.position.ic = 0.25`

## Comparison P-values

By default, `compare_motifs()` will compare all motifs provided and return a matrix. The `compare.to` argument will cause `compare_motifs()` to return P-values.

```{r}
library(universalmotif)
library(MotifDb)
motifs <- filter_motifs(MotifDb, organism = "Athaliana")

# Compare the first motif with everything and return P-values
head(compare_motifs(motifs, 1))
```

P-values are made possible by estimating distribution (usually the best fitting distribution for motif comparisons) parameters from randomised motif scores, then using the appropriate `stats::p*()` distribution function to return P-values. These estimated parameters are pre-computed with `make_DBscores()` and stored as `JASPAR2018_CORE_DBSCORES` and `JASPAR2018_CORE_DBSCORES_NORM`. Since changing any of the settings and motif sizes will affect the estimated distribution parameters, estimated parameters have been pre-computed for a variety of these. See `?make_DBscores` if you would like to generate your own set of pre-computed scores using your own parameters and motifs.

## A faster alternative: `compare_motifs2()`

For DNA or RNA motif comparison that does not need the alternative
similarity metrics (EUCL, KL, ALLR, and so on), the score-aggregation
strategies, the IC-based filters, multifreq scoring, or the
JASPAR-derived `db.scores` lookup-table P-value path, the lighter
`compare_motifs2()` is generally much faster than `compare_motifs()`,
and it parallelises rather better across queries. It uses the
per-column Pearson correlation as its only similarity metric, computes
p-values from a TOMTOM-style empirical or parametric null PMF (with
both Bonferroni and Benjamini-Hochberg adjustment), and returns either
a square matrix or a long-format `data.frame`. Its default surface is
aligned with the command-line tool
[yamtk](https://github.com/bjmt/yamtk).

The function has eight formals (`motifs`, `compare.to`, `qvalue`,
`min.overlap`, `RC`, `null`, `bkg`, `matrix.out`, `nthreads`); see
`?compare_motifs2` for the details. There are two null-PMF modes:

- **`null = "empirical"`** (default): per-query column, score against
  every target-database column, histogram into 51 PCC bins, convolve
  to overlap length $L$. Cost per query scales with database size;
  best for databases with more than a few hundred total columns.
- **`null = "parametric"`**: enumerate the 56 compositions of A/C/G/T
  on a $K=5$ simplex grid, weight each by a Dirichlet-Multinomial PMF
  under $\alpha = 4 \times \mathtt{bkg}$, score against the query
  column. Per-query cost is constant in database size; best for small
  databases or when histogram sparsity would distort the empirical
  estimate.

```{r, eval = FALSE}
## Score matrix (default: matrix.out = "score" returns mean PCC,
## always in [-1, 1])
m <- compare_motifs2(motifs)

## Significance matrix (p-values or BH q-values across the
## query-vs-database axis)
mp <- compare_motifs2(motifs, matrix.out = "pvalue")
mq <- compare_motifs2(motifs, matrix.out = "qvalue")

## Long-format: motif 1 against the rest, q-value <= 0.05
hits <- compare_motifs2(motifs, compare.to = 1, qvalue = 0.05)
```

When `compare_motifs()` is called with arguments that map cleanly onto
`compare_motifs2()`'s feature set (PCC + sum strategy, the default IC
filters, no multifreq, no HTML report, no `db.scores`, and a DNA/RNA
alphabet), it emits a one-line hint pointing you at the faster
function. (Silence it with
`options(universalmotif.suggest.compare_motifs2 = FALSE)`.)

### Speed comparison vs `compare_motifs()` and `yamtk cmp` {.unnumbered}

The table below records the wall-clock time for the three functions
on a single shared fixture: the first 50 (or 200) motifs of
HOCOMOCOv11, compared against themselves, with RC = TRUE,
min.overlap = 5, a uniform background, and the empirical null PMF. The
`compare_motifs2()` numbers are shown for both `matrix.out = "score"`
(which skips the p-value path entirely, since the best alignment is
chosen by score rather than by p-value) and `matrix.out = "pvalue"`
(the full significance computation, which is equivalent work to
`yamtk cmp`). Each cell shows `nthreads = 1 / nthreads = 4`. The
benchmark script lives at `benchmarks/compare_motifs2_vignette.log` in
the package repository, and the values are baked in here rather than
recomputed at vignette build time. These benchmarks were run on an
MBP M5 under R 4.4.1.

| Workload | yamtk¹ | `compare_motifs()` | `cm2()` score² | `cm2()` pvalue |
|---|---:|---:|---:|---:|
| 50 × 50 (2 500 pairs) | 0.028 / 0.012 | 0.598 / 0.194 | 0.014 / 0.004 | 0.285 / 0.082 |
| 200 × 200 (40 000 pairs) | 0.114 / 0.038 | 9.537 / 2.590 | 0.224 / 0.064 | 1.673 / 0.470 |

`cm2()` is shorthand for `compare_motifs2()`; the `score` / `pvalue`
labels indicate the value of the `matrix.out` argument.

¹ `yamtk cmp` is a dedicated C CLI tool; threading is by query motif.
The `compare_motifs2()` algorithm is a port of yamtk's.

² Score-only matrix mode skips the per-pair p-value lookup, which
otherwise dominates runtime on dense fixtures. This is safe because
the best alignment is selected by score, not p-value; p-values are
computed only as a post-hoc significance check on the already-chosen
alignment.

To summarise:

- `yamtk cmp` is the fastest of the three at every scale; the
  `compare_motifs2()` port stays within ~2× at 4 threads for score-only
  output and within ~12× for full p-value output.
- `compare_motifs2()` is 5-40× faster than `compare_motifs()` at
  4 threads on the same set, with the larger speedups in score-only
  mode. `compare_motifs()`'s p-value path uses pre-computed JASPAR
  score-distribution tables (`JASPAR2018_CORE_DBSCORES`) and a different
  null model, so the numbers are not directly comparable as
  "significance" but are comparable as "wall-clock to a similarity
  matrix".
- Threading scales 3-3.5× at 4 threads in all three universalmotif
  modes.
- For new code that does not need any of `compare_motifs()`'s
  advanced features, try `compare_motifs2()`.

# Merging motifs

The `universalmotif` package offers two ways to collapse similar
motifs into consensus motifs: `merge_motifs()` / `merge_similar()`,
and `merge_motifs2()` / `merge_similar2()` (the latter pair are
newer, built on `compare_motifs2()`, and use significance-based
clustering).

## A faster alternative: `merge_motifs2()` / `merge_similar2()`

For simple DNA/RNA motif merging tasks:

- `merge_motifs2(motifs)` aligns every motif against the highest-IC
  anchor via `compare_motifs2()` and column-averages them in a single
  shared frame.

- `merge_similar2(motifs, qvalue)` clusters motifs by statistical
  significance: two motifs are linked if their pairwise q-value
  meets the user's `qvalue` cutoff. Clustering is by connected
  components on the resulting graph (deterministic, no linkage-method
  to choose). Each component is collapsed via `merge_motifs2()`.

```{r}
library(universalmotif)
m1 <- create_motif("TTGACATA", name = "a")
m2 <- create_motif("CTTGACAT", name = "b")
m3 <- create_motif("TGACATAT", name = "c")
m4 <- create_motif("GGGCCCCC", name = "unrelated")

## Single merged consensus across 3 related motifs:
merge_motifs2(list(m1, m2, m3))

## Cluster + merge: returns a shorter list, with similar motifs collapsed.
merge_similar2(list(m1, m2, m3, m4), qvalue = 0.05)

## Or just the cluster assignments (no merging):
merge_similar2(list(m1, m2, m3, m4), qvalue = 0.05, return.clusters = TRUE)
```

When a `merge_motifs()` or `merge_similar()` call uses only features
that are also available to `merge_motifs2()` / `merge_similar2()`, a
one-line hint is emitted pointing you at the leaner function. (Silence
it with `options(universalmotif.suggest.merge_motifs2 = FALSE)` and
`options(universalmotif.suggest.merge_similar2 = FALSE)`.)

# Motif trees with ggtree

## Using `motif_tree()`

Additionally, this package introduces the `motif_tree()` function for generating basic tree-like diagrams for comparing motifs. This allows for a visual result from `compare_motifs()`. All options from `compare_motifs()` are available in `motif_tree()`. This function uses the `ggtree` package and outputs a `ggplot` object (from the `ggplot2` package), so altering the look of the trees can be done easily after `motif_tree()` has already been run. There is also a new `motif_tree2()`, which is based on `compare_motifs2()`.

```{r}
library(universalmotif)
library(MotifDb)

motifs <- filter_motifs(MotifDb, family = c("AP2", "B3", "bHLH", "bZIP",
                                            "AT hook"))
motifs <- motifs[sample(seq_along(motifs), 100)]
tree <- motif_tree(motifs, layout = "daylight", linecol = "family")

## Make some changes to the tree in regular ggplot2 fashion:
# tree <- tree + ...

tree
```

## Using `compare_motifs()` and `ggtree()`

While `motif_tree()` works as a quick and convenient tree-building function, it can be inconvenient when more control is required over tree construction. For this purpose, the following code goes through how exactly `motif_tree()` generates trees.

```{r}
library(universalmotif)
library(MotifDb)
library(ggtree)
library(ggplot2)

motifs <- convert_motifs(MotifDb)
motifs <- filter_motifs(motifs, organism = "Athaliana")
motifs <- motifs[sample(seq_along(motifs), 25)]

## Step 1: compare motifs

comparisons <- compare_motifs(motifs, method = "PCC", min.mean.ic = 0,
                              score.strat = "a.mean")

## Step 2: create a "dist" object

# The current metric, PCC, is a similarity metric
comparisons <- 1 - comparisons

comparisons <- as.dist(comparisons)

# We also want to extract names from the dist object to match annotations
labels <- attr(comparisons, "Labels")

## Step 3: get the comparisons ready for tree-building

# The R package "ape" provides the necessary "as.phylo" function
comparisons <- ape::as.phylo(hclust(comparisons))

## Step 4: incorporate annotation data to colour tree lines

family <- sapply(motifs, function(x) x["family"])
family.unique <- unique(family)

# We need to create a list with an entry for each family; within each entry
# are the names of the motifs belonging to that family
family.annotations <- list()
for (i in seq_along(family.unique)) {
  family.annotations <- c(family.annotations,
                          list(labels[family %in% family.unique[i]]))
}
names(family.annotations) <- family.unique

# Now add the annotation data:
comparisons <- ggtree::groupOTU(comparisons, family.annotations)

## Step 5: draw the tree

tree <- ggtree(comparisons, aes(colour = group), layout = "rectangular") +
          theme(legend.position = "bottom", legend.title = element_blank())

## Step 6: add additional annotations

# If we wish, we can add additional annotations such as tip labelling and size

# Tip labels:
tree <- tree + geom_tiplab()

# Tip size:
tipsize <- data.frame(label = labels,
                      icscore = sapply(motifs, function(x) x["icscore"]))

tree <- tree %<+% tipsize + geom_tippoint(aes(size = icscore))

```

## Plotting motifs alongside trees

Unfortunately, the `universalmotif` package does not provide any function to easily plot motifs as part of trees (as is possible via the `motifStack` package). However, it can be done (somewhat roughly) by plotting a tree and a set of motifs side by side. In the following example, the `cowplot` package is used to glue the two plots together, though other packages which perform this function are available.

```{r}
library(universalmotif)
library(MotifDb)
library(cowplot)

## Get our starting set of motifs:
motifs <- convert_motifs(MotifDb[1:10])

## Get the tree: make sure it's a horizontal type layout
tree <- motif_tree(motifs, layout = "rectangular", linecol = "none")

## Now, make sure we order our list of motifs to match the order of tips:
mot.names <- sapply(motifs, function(x) x["name"])
names(motifs) <- mot.names
new.order <- tree$data$label[tree$data$isTip]
new.order <- rev(new.order[order(tree$data$y[tree$data$isTip])])
motifs <- motifs[new.order]

## Plot the two together (finessing of margins and positions may be required):
plot_grid(nrow = 1, rel_widths = c(1, -0.15, 1),
  tree + xlab(""), NULL,
  view_motifs(motifs, names.pos = "right") +
    ylab(element_blank()) +
    theme(
      axis.line.y = element_blank(),
      axis.ticks.y = element_blank(),
      axis.text.y = element_blank(),
      axis.text = element_text(colour = "white")
    )
)
```

# Motif P-values

Motif P-values are not usually discussed outside of the bioinformatics literature, but are actually quite a challenging topic. To illustrate this, consider the following example motif:

```{r}
library(universalmotif)

m <- matrix(c(0.10,0.27,0.23,0.19,0.29,0.28,0.51,0.12,0.34,0.26,
              0.36,0.29,0.51,0.38,0.23,0.16,0.17,0.21,0.23,0.36,
              0.45,0.05,0.02,0.13,0.27,0.38,0.26,0.38,0.12,0.31,
              0.09,0.40,0.24,0.30,0.21,0.19,0.05,0.30,0.31,0.08),
            byrow = TRUE, nrow = 4)
motif <- create_motif(m, alphabet = "DNA", type = "PWM")
motif
```

Let us then use this motif with `scan_sequences()`:

```{r}
data(ArabidopsisPromoters)

res <- scan_sequences(motif, ArabidopsisPromoters, verbose = 0,
  calc.pvals = FALSE, threshold = 0.8, threshold.type = "logodds")
head(res)
```

Now let us imagine that we wish to rank these matches by P-value. First, we must calculate the match probabilities:

```{r}
## One of the matches was CTCTAGAGAC, with a score of 5.869 (max possible = 6.531)

bkg <- get_bkg(ArabidopsisPromoters, 1)
bkg <- structure(bkg$probability, names = bkg$klet)
bkg
```
Now, use these to calculate the probability of getting CTCTAGAGAC.
```{r}
hit.prob <- bkg["A"]^3 * bkg["C"]^3 * bkg["G"]^2 * bkg["T"]^2
hit.prob <- unname(hit.prob)
hit.prob
```

Calculating the probability of a single match was easy, but then comes the challenging part: calculating the probability of all possible matches with a score higher than 5.869, and then summing these. This final sum then represents the probability of finding a match which scores at least 5.869. One way is to list all possible sequence combinations, then filtering based on score; however, this "brute force" approach is unreasonable for all but the smallest of motifs.

## The dynamic programming algorithm for calculating P-values and scores

Instead of trying to find and calculate the probabilities of all matches with a score equal to or higher than the query score, one can use a dynamic programming algorithm to generate a much smaller distribution of probabilities for the possible range of scores using set intervals. This method is implemented by the FIMO tool [@grant11]. The theory behind it is also explained in @gupta07, though the purpose of the algorithm is for motif comparison instead of motif P-values (however it is the same algorithm). The basic concept will also be briefly explained here.

For each individual position-letter score in the PWM, the chance of getting that score from the respective background probability of that letter is added to the intervals in which getting that specific score could allow the final score to land. Once this probability distribution is generated, it can be converted to a cumulative distribution and re-used for any input P-value/score to output the equivalent score/P-value. For P-value inputs, it finds the specific score interval where the accompanying P-value in the cumulative distribution is smaller than or equal to it, then reports the score of the previous interval. For score inputs, the scores are rounded to the nearest interval in the cumulative distribution and the accompanying P-value retrieved. The major advantages of this method include only looking for the probabilities of the range of scores with a set interval, cutting down on needing to find the probabilities of all actual possible scores (and thus increasing performance by several orders of magnitude for larger/higher-order motifs), and being able to re-use the distribution for any number of query P-value/scores. Although this method involves rounding off scores to allow a small set interval, in practice within the `universalmotif` package it offers the same maximum possible level of accuracy as the exhaustive method (described in the next section) as motif PWMs are always internally rounded to a thousandth of a decimal place for speed. This leaves as the only downside the inability to allow non-finite values to exist in the PWM (e.g. from zero-probabilities) since then a known range with set intervals could not possibly be created.

Going back to our example, we can see this in action using the `motif_pvalue()` function:

```{r}
res <- res[1:6, ]
pvals <- motif_pvalue(motif, res$score, bkg.probs = bkg)
res2 <- data.frame(motif=res$motif,match=res$match,pval=pvals)[order(pvals), ]
knitr::kable(res2, digits = 22, row.names = FALSE, format = "markdown")
```

To illustrate that we can also do the inverse of this calculation:

```{r}
res$score
motif_pvalue(motif, pvalue = pvals, bkg.probs = bkg)
```

You may occasionally see slight errors at the last couple of digits. These are generally unavoidable due to the internal rounding mechanisms of the `universalmotif` package.

Let us consider more examples, such as the following larger motif:

```{r}
data(ArabidopsisMotif)
ArabidopsisMotif
```

Using the `motif_range()` utility, we can get an idea of the possible range of scores:

```{r}
motif_range(ArabidopsisMotif)
```

We can use these ranges to confirm our cumulative distribution of P-values:

```{r}
(pvals2 <- motif_pvalue(ArabidopsisMotif, score = motif_range(ArabidopsisMotif)))
```

And again, going back to scores from these P-values:

```{r}
motif_pvalue(ArabidopsisMotif, pvalue = pvals2)
```

As a note: if you ever provide scores which are outside the possible ranges, then you will get the following behaviour:

```{r}
motif_pvalue(ArabidopsisMotif, score = c(-200, 100))
```

We can also use this function for the higher-order `multifreq` motif representation.

```{r}
data(examplemotif2)
examplemotif2["multifreq"]["2"]
motif_range(examplemotif2, use.freq = 2)
motif_pvalue(examplemotif2, score = 15, use.freq = 2)
motif_pvalue(examplemotif2, pvalue = 0.00001, use.freq = 2)
```

Feel free to use this function with any alphabets, such as amino acid motifs or even made up ones!

```{r}
(m <- create_motif(alphabet = "QWERTY"))
motif_pvalue(m, pvalue = c(1, 0.1, 0.001, 0.0001, 0.00001))
```

## The branch-and-bound algorithm for calculating P-values from scores

The alternative to the dynamic programming algorithm is to exhaustively find all actual possible hits with a score equal to or greater than the input score. Generally there is no advantage to solving this exhaustively, with the exception that it allows non-finite values to be present (i.e., zero-probability letters which were not pseudocount-adjusted during the calculation of the PWM). A few algorithms have been proposed to make solving this problem exhaustively more efficient, but the method adopted by the `universalmotif` package is that of @pvalues. The authors propose using a branch-and-bound^[https://en.wikipedia.org/wiki/Branch_and_bound] algorithm (with a few tricks) alongside a certain approximation. Briefly: motifs are first reorganized so that the highest scoring positions and letters are considered first in the branch-and-bound algorithm. Then, motifs past a certain width (in the original paper, 10) are split in sub-motifs. All possible combinations are found in these sub-motifs using the branch-and-bound algorithm, and P-values calculated for the sub-motifs. Finally, the P-values are combined.

The `motif_pvalue()` function modifies this process slightly by allowing the size of the sub-motifs to be specified via the `k` parameter; and additionally, whereas the original implementation can only calculate P-values for motifs with a maximum of 17 positions (and motifs can only be split in at most two), the `universalmotif` implementation allows for any length of motif to be used (and motifs can be split any number of times). Changing `k` allows one to decide between speed and accuracy; smaller `k` leads to faster but worse approximations, and larger `k` leads to slower but better approximations. If `k` is equal to the width of the motif, then the calculation is _exact_. It is important to note however that this is still a computationally intensive task for larger motifs unless it is broken up into several sub-motifs, though at this point significant accuracy is lost due to the high level of approximation.

Now, let us return to our original example, and this time for the branch-and-bound algorithm set `method = "exhaustive"`:

```{r}
res <- res[1:6, ]
pvals <- motif_pvalue(motif, res$score, bkg.probs = bkg, method = "e")
res2 <- data.frame(motif=res$motif,match=res$match,pval=pvals)[order(pvals), ]
knitr::kable(res2, digits = 22, row.names = FALSE, format = "markdown")
```

The default `k` in `motif_pvalue()` is 8. I have found this to be a good trade-off between speed and P-value correctness.

To demonstrate the effect that `k` has on the output P-value, consider the following (and also note that for this motif `k = 10` represents an exact calculation):

```{r}
scores <- c(-6, -3, 0, 3, 6)
k <- c(2, 4, 6, 8, 10)
out <- data.frame(k = c(2, 4, 6, 8, 10),
                  score.minus6 = rep(0, 5),
                  score.minus3 = rep(0, 5),
                  score.0 = rep(0, 5),
                  score.3 = rep(0, 5),
                  score.6 = rep(0, 5))

for (i in seq_along(scores)) {
  for (j in seq_along(k)) {
    out[j, i + 1] <- motif_pvalue(motif, scores[i], k = k[j], bkg.probs = bkg,
      method = "e")
  }
}

knitr::kable(out, format = "markdown", digits = 10)
```

For this particular motif, while the approximation worsens slightly as `k` decreases, it is still quite accurate when the number of motif subsets is limited to two. Usually, you should only have to worry about `k` for longer motifs (such as those sometimes generated by `MEME`), where the number of sub-motifs increases.

## The random subsetting algorithm for calculating scores from P-values

Similarly to calculating P-values, exact scores can be calculated from small motifs, and approximate scores from larger motifs using subsetting. When an exact calculation is performed, all possible scores are extracted and a quantile function extracts the appropriate score. For approximate calculations, the overall set of scores is approximated several times by randomly adding up all possible scores from each `k` subset before a quantile function is used.

Starting from a set of P-values and setting `method = "exhaustive"`:

```{r}
bkg <- c(A=0.25, C=0.25, G=0.25, T=0.25)
pvals <- c(0.1, 0.01, 0.001, 0.0001, 0.00001)
scores <- motif_pvalue(motif, pvalue = pvals, bkg.probs = bkg, k = 10,
  method = "e")

scores.approx6 <- motif_pvalue(motif, pvalue = pvals, bkg.probs = bkg, k = 6,
  method = "e")
scores.approx8 <- motif_pvalue(motif, pvalue = pvals, bkg.probs = bkg, k = 8,
  method = "e")

pvals.exact <- motif_pvalue(motif, score = scores, bkg.probs = bkg, k = 10,
  method = "e")

pvals.approx6 <- motif_pvalue(motif, score = scores, bkg.probs = bkg, k = 6,
  method = "e")
pvals.approx8 <- motif_pvalue(motif, score = scores, bkg.probs = bkg, k = 8,
  method = "e")

res <- data.frame(pvalue = pvals, score = scores,
                  pvalue.exact = pvals.exact,
                  pvalue.k6 = pvals.approx6,
                  pvalue.k8 = pvals.approx8,
                  score.k6 = scores.approx6,
                  score.k8 = scores.approx8)
knitr::kable(res, format = "markdown", digits = 22)
```

Starting from a set of scores:

```{r}
bkg <- c(A=0.25, C=0.25, G=0.25, T=0.25)
scores <- -2:6
pvals <- motif_pvalue(motif, score = scores, bkg.probs = bkg, k = 10,
  method = "e")

scores.exact <- motif_pvalue(motif, pvalue = pvals, bkg.probs = bkg, k = 10,
  method = "e")

scores.approx6 <- motif_pvalue(motif, pvalue = pvals, bkg.probs = bkg, k = 6,
  method = "e")
scores.approx8 <- motif_pvalue(motif, pvalue = pvals, bkg.probs = bkg, k = 8,
  method = "e")

pvals.approx6 <- motif_pvalue(motif, score = scores, bkg.probs = bkg, k = 6,
  method = "e")
pvals.approx8 <- motif_pvalue(motif, score = scores, bkg.probs = bkg, k = 8,
  method = "e")

res <- data.frame(score = scores, pvalue = pvals,
                  pvalue.k6 = pvals.approx6,
                  pvalue.k8 = pvals.approx8,
                  score.exact = scores.exact,
                  score.k6 = scores.approx6,
                  score.k8 = scores.approx8)
knitr::kable(res, format = "markdown", digits = 22)
```
As you may have noticed, results from exact calculations are not _quite_ exact. This is due to the `universalmotif` package rounding off values internally for speed.

# Session info {.unnumbered}

```{r sessionInfo, echo=FALSE}
sessionInfo()
```

# References {.unnumbered}
