Nothing
#' Create P-value databases.
#'
#' Generate data used by [compare_motifs()] for P-value calculations. By default,
#' [compare_motifs()] uses an internal database based on the JASPAR2018 core motifs
#' \insertCite{jaspar}{universalmotif}. Parameters for distributions are
#' are estimated for every combination of motif `widths`.
#'
#' @param db.motifs `list` Database motifs.
#' @param shuffle.db `logical(1)` Deprecated. Does nothing.
#' generate random motifs with [create_motif()].
#' @param shuffle.k `numeric(1)` See [shuffle_motifs()].
#' @param shuffle.method `character(1)` See [shuffle_motifs()].
#' @param rand.tries `numeric(1)` Approximate number of comparisons
#' to perform for every combination of `widths`.
#' @param widths `numeric` Motif widths to use in P-value database calculation.
#' @param progress `logical(1)` Show progress.
#' @param nthreads `numeric(1)` Run [compare_motifs()] in parallel with `nthreads`
#' threads. `nthreads = 0` uses all available threads.
#'
#' @return A `DataFrame` with score distributions for the
#' input database. If more than one [make_DBscores()] run occurs (i.e. args
#' `method`, `normalise.scores` or `score.strat` are longer than 1), then
#' the function args are included in the `metadata` slot.
#'
#' @details
#' See [compare_motifs()] for more info on comparison parameters.
#'
#' To replicate the internal \pkg{universalmotif} DB scores, run
#' [make_DBscores()] with the default settings. Note that this will be
#' a slow process.
#'
#' Arguments `widths`, `method`, `normalise.scores` and `score.strat` are
#' vectorized; all combinations will be attempted.
#'
#' @examples
#' \dontrun{
#' library(MotifDb)
#' motifs <- convert_motifs(MotifDb[1:100])
#' scores <- make_DBscores(motifs, method = "PCC")
#' compare_motifs(motifs, 1:100, db.scores = scores)
#' }
#'
#' @references
#' \insertRef{jaspar}{universalmotif}
#'
#' @seealso [compare_motifs()]
#' @author Benjamin Jean-Marie Tremblay, \email{b2tremblay@@uwaterloo.ca}
#' @inheritParams compare_motifs
#' @export
make_DBscores <- function(db.motifs,
method = c("PCC", "EUCL", "SW", "KL", "WEUCL",
"ALLR", "BHAT", "HELL", "WPCC", "SEUCL",
"MAN", "ALLR_LL"),
shuffle.db = TRUE,
shuffle.k = 3, shuffle.method = "linear",
rand.tries = 1000, widths = 5:30,
min.position.ic = 0,
normalise.scores = c(FALSE, TRUE), min.overlap = 6,
min.mean.ic = 0.25, progress = TRUE,
nthreads = 1, tryRC = TRUE,
score.strat = c("sum", "a.mean", "g.mean", "median",
"wa.mean", "wg.mean", "fzt")) {
args <- as.list(environment())
if (length(method) > 1 || length(normalise.scores) > 1 || length(score.strat) > 1) {
# Need to vectorize through all three args and combine into data.frames, one per
# method. Finally, all data.frames are kept in a list with an additional entry
# for function args.
out <- vector("list", length(method) + 1)
names(out) <- c(method, "args")
mc <- 1
total <- length(method) * length(normalise.scores) * length(score.strat)
if (progress) t1 <- Sys.time()
for (m in method) {
out[[m]] <- DataFrame(subject = integer(), target = integer(),
paramA = numeric(), paramB = numeric(),
method = character(), normalised = logical(),
strat = character(), distribution = character())
for (norm in normalise.scores) {
for (strat in score.strat) {
if (progress) message("[", mc, "/", total, "]", " method=\"", m,
"\" normalise.scores=", norm, " score.strat=\"",
strat, "\" ", appendLF = FALSE)
mc <- mc + 1
if (strat %in% c("g.mean", "wg.mean") && m %in% c("ALLR", "ALLR_LL", "PCC")) {
if (progress)
message("\n > Skipping: (w)g.mean not allowed with ALLR/ALLR_LL/PCC\n")
next
}
if (progress) start <- Sys.time()
tmp <- make_DBscores(db.motifs, m, shuffle.db, shuffle.k, shuffle.method,
rand.tries, widths, min.position.ic,
norm, min.overlap, min.mean.ic,
progress, nthreads, tryRC, strat)
out[[m]] <- rbind(out[[m]], tmp)
if (progress) stop <- Sys.time()
if (progress) message(" > ", format(difftime(stop, start)), "\n")
}
}
}
if (progress) t2 <- Sys.time()
if (progress) message(" *** Total runtime: ", format(difftime(t2, t1)), " ***")
out <- do.call(rbind, out)
out <- out[order(as.vector(out$method), as.vector(out$normalised),
as.vector(out$strat), as.vector(out$distribution),
as.vector(out$subject), as.vector(out$target)), ]
out[, 1] <- Rle(out[, 1])
out[, 5] <- Rle(out[, 5])
out[, 6] <- Rle(out[, 6])
out[, 7] <- Rle(out[, 7])
out[, 8] <- Rle(out[, 8])
out@metadata <- args[-1]
return(out)
}
# param check --------------------------------------------
method <- match.arg(method, COMPARE_METRICS)
char_check <- check_fun_params(list(method = args$method,
shuffle.method = args$shuffle.method,
score.strat = args$score.strat),
c(0, 1, 1), logical(), TYPE_CHAR)
num_check <- check_fun_params(list(shuffle.k = args$shuffle.k,
rand.tries = args$rand.tries,
min.overlap = args$min.overlap,
min.mean.ic = args$min.mean.ic,
nthreads = args$nthreads,
min.position.ic = args$min.position.ic),
numeric(), logical(), TYPE_NUM)
logi_check <- check_fun_params(list(shuffle.db = args$shuffle.db,
progress = args$progress,
normalise.scores = args$normalise.scores,
tryRC = args$tryRC),
numeric(), logical(), TYPE_LOGI)
all_checks <- c(char_check, num_check, logi_check)
if (length(all_checks) > 0) stop(all_checks_collapse(all_checks))
#---------------------------------------------------------
# having min.mean.ic > 0 can really mess with scores sometimes
if (score.strat %in% c("g.mean", "wg.mean") && method %in%
c("ALLR", "ALLR_LL", "PCC"))
stop(wmsg("'g.mean'/'wg.mean' is not allowed for methods which can generate negative ",
"values: ALLR, ALLR_LL, PCC"))
if (!score.strat %in% c("sum", "a.mean", "g.mean", "median", "wa.mean",
"wg.mean", "fzt"))
stop("'score.strat' must be one of 'sum', 'a.mean', 'g.mean', 'median', ",
"'wa.mean', 'wg.mean', 'fzt'")
rand.mots <- make_DBscores_motifs(db.motifs, widths, rand.tries, shuffle.k)
comps <- get_comparisons(widths)
total <- nrow(comps)
totry <- DataFrame(subject = comps[, 1], target = comps[, 2],
paramA = rep(NA_real_, total),
paramB = rep(NA_real_, total),
method = rep(method, total),
normalised = rep(normalise.scores, total),
strat = rep(score.strat, total),
distribution = rep("logistic", total))
if (progress) {
print_pb(0)
}
counter <- 1
for (i in seq_len(total)) {
subject <- rand.mots[[as.character(totry$subject[i])]]
target <- rand.mots[[as.character(totry$target[i])]]
tmpall <- c(subject, target)
mtmp <- lapply(tmpall, function(x) x@motif)
comps <- get_comp_indices(seq_along(subject), length(tmpall))
comps <- comps[!comps[, 2] %in% seq_along(subject), ]
res <- compare_motifs_cpp(mtmp, comps[, 1] - 1, comps[, 2] - 1,
method, min.overlap, tryRC,
get_bkgs(tmpall), 1, FALSE, min.mean.ic,
normalise.scores, nthreads,
min.position.ic, get_nsites(tmpall),
score.strat)
if (length(unique(res)) == 1)
stop(wmsg("failed to estimate distribution due to uniform random scores ",
"at comparison: ", as.vector(totry$subject[i]), " - ",
as.vector(totry$target[i]), " ; ",
"perhaps too few reference motifs"))
res[res == min_max_doubles()$min | res == min_max_doubles()$max] <- NA_real_
res <- res[!is.na(res)]
if (length(res) <= 1)
stop(wmsg("failed to obtain scores due to low motif IC at comparison: ",
as.vector(totry$subject[i]), " - ", as.vector(totry$target[i]),
" ; perhaps too few reference motifs or too many low IC motifs"))
d <- find_best_dist(res)
totry$paramA[i] <- d$paramA
totry$paramB[i] <- d$paramB
totry$distribution[i] <- d$dist
if (progress) update_pb(counter, total)
counter <- counter + 1
}
totry <- totry[order(totry$method, totry$normalised, totry$strat,
totry$distribution, totry$subject, totry$target), ]
totry[, 1] <- Rle(totry[, 1])
totry[, 5] <- Rle(totry[, 5])
totry[, 6] <- Rle(totry[, 6])
totry[, 7] <- Rle(totry[, 7])
totry[, 8] <- Rle(totry[, 8])
rownames(totry) <- NULL
totry
}
find_best_dist <- function(x) {
d <- c("normal", "logistic", "weibull")
pvals <- numeric(length(d))
# NORMAL
n <- tryCatch(suppressWarnings(fitdistr(x, "normal")), error = function(e) FALSE)
if (isFALSE(n)) {
pvals[1] <- -1
} else {
pvals[1] <- suppressWarnings(ks.test(x, "pnorm", n$estimate["mean"],
n$estimate["sd"])$p)
}
# LOGISTIC
l <- tryCatch(suppressWarnings(fitdistr(x, "logistic")), error = function(e) FALSE)
if (isFALSE(l)) {
pvals[2] <- -1
} else {
pvals[2] <- suppressWarnings(ks.test(x, "plogis", l$estimate["location"],
l$estimate["scale"])$p)
}
# WEIBULL
w <- tryCatch(suppressWarnings(fitdistr(x, "weibull")), error = function(e) FALSE)
if (isFALSE(w)) {
pvals[3] <- -1
} else {
pvals[3] <- suppressWarnings(ks.test(x, "pweibull", w$estimate["shape"],
w$estimate["scale"])$p)
}
dist <- d[which.max(pvals)[1]]
switch(dist,
"normal" = list(dist = "normal", paramA = n$estimate["mean"],
paramB = n$estimate["sd"]),
"logistic" = list(dist = "logistic", paramA = l$estimate["location"],
paramB = l$estimate["scale"]),
"weibull" = list(dist = "weibull", paramA = w$estimate["shape"],
paramB = w$estimate["scale"]))
}
get_comparisons <- function(widths) {
out <- do.call(rbind, combn(widths, 2, simplify = FALSE))
rbind(out, matrix(c(widths, widths), ncol = 2))
}
make_DBscores_motifs <- function(motifs, widths, rand.tries, shuffle.k) {
mot.mats <- lapply(motifs, function(x) x@motif)
mot.bkg <- get_bkgs(motifs)
mot.nsites <- get_nsites(motifs)
alph <- unique(vapply(motifs, function(x) x@alphabet, character(1)))
if (length(alph) > 1) stop("all motifs must share the same alphabet")
comb.mats <- do.call(cbind, mot.mats)
if (ncol(comb.mats) < max(widths)) {
while (ncol(comb.mats) < max(widths)) {
comb.mats <- cbind(comb.mats, comb.mats)
}
}
if (shuffle.k == 1) comb.mats <- comb.mats[, sample(seq_len(ncol(comb.mats)))]
else {
old.order <- as.character(seq_len(ncol(comb.mats)))
new.order <- shuffle_linear(old.order, shuffle.k, mode = 2)
comb.mats <- comb.mats[, as.numeric(new.order)]
}
nmots <- round(sqrt(rand.tries))
out <- vector("list", length(widths))
names(out) <- as.character(widths)
for (i in seq_along(widths)) {
j <- sample.int(length(mot.bkg), 1)
out[[i]] <- lapply(seq_len(nmots),
function(x) sample_motif_(comb.mats, widths[i],
mot.bkg[[j]], mot.nsites[[j]],
alph))
}
out
}
sample_motif_ <- function(pool, size, bkg, nsites, alph) {
n <- ncol(pool)
n <- sample.int(n - size + 1, 1)
mat <- pool[, seq(n, n + size - 1), drop = FALSE]
create_motif(mat, alphabet = alph, bkg = bkg, nsites = nsites)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
make_DBscores_v1 <- function(db.motifs, method, shuffle.db = TRUE,
shuffle.k = 3, shuffle.method = "linear",
shuffle.leftovers = "asis", rand.tries = 1000,
normalise.scores = TRUE, min.overlap = 6,
min.mean.ic = 0, progress = TRUE, BP = FALSE) {
# param check --------------------------------------------
args <- as.list(environment())
char_check <- check_fun_params(list(method = args$method,
shuffle.method = args$shuffle.method,
shuffle.leftovers = args$shuffle.leftovers),
numeric(), logical(), TYPE_CHAR)
num_check <- check_fun_params(list(shuffle.k = args$shuffle.k,
rand.tries = args$rand.tries,
min.overlap = args$min.overlap,
min.mean.ic = args$min.mean.ic),
numeric(), logical(), TYPE_NUM)
logi_check <- check_fun_params(list(shuffle.db = args$shuffle.db,
progress = args$progress, BP = args$BP,
normalise.scores = args$normalise.scores),
numeric(), logical(), TYPE_LOGI)
all_checks <- c(char_check, num_check, logi_check)
if (length(all_checks) > 0) stop(all_checks_collapse(all_checks))
#---------------------------------------------------------
db.motifs <- convert_motifs(db.motifs)
db.ncols <- vapply(db.motifs, function(x) ncol(x@motif), numeric(1))
if (shuffle.db) {
rand.mots <- shuffle_motifs(db.motifs, k = shuffle.k,
method = shuffle.method)
if (length(rand.mots) != rand.tries) {
if (length(rand.mots) < rand.tries) {
while (length(rand.mots) < rand.tries) {
more.rand.mots <- shuffle_motifs(db.motifs, k = shuffle.k,
method = shuffle.method)
rand.mots <- c(rand.mots, more.rand.mots)
}
}
if (length(rand.mots) > rand.tries) {
rand.mots <- rand.mots[sample(seq_along(rand.mots), rand.tries)]
}
}
} else {
rand.mots <- lapply(seq_len(rand.tries),
function(x) create_motif(sample.int(26, 1) + 4))
}
rand.ncols <- vapply(rand.mots, function(x) ncol(x@motif), numeric(1))
totry <- expand.grid(list(subject = sort(unique(rand.ncols)),
target = sort(unique(db.ncols))))
totry$mean <- rep(NA, nrow(totry))
totry$sd <- rep(NA, nrow(totry))
res <- vector("list", nrow(totry))
if (progress) print_pb(0)
for (i in seq_len(nrow(totry))) {
tmp1 <- db.motifs[totry[i, 2] == db.ncols]
tmp2 <- rand.mots[totry[i, 1] == rand.ncols]
res[[i]] <- compare_motifs(c(tmp2, tmp1), seq_along(tmp2), method = method,
min.overlap = min.overlap, min.mean.ic = min.mean.ic,
max.e = Inf, max.p = Inf,
normalise.scores = normalise.scores)$score
totry$mean[i] <- mean(res[[i]])
totry$sd[i] <- sd(res[[i]])
if (progress) update_pb(i, nrow(totry))
}
totry$method <- rep(method, nrow(totry))
totry$normalised <- rep(normalise.scores, nrow(totry))
totry
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
make_DBscores_v3 <- function(db.motifs, rand.tries = 10000, widths = 1:100,
progress = TRUE, BP = FALSE) {
mots <- lapply(db.motifs, function(x) x@motif)
cols <- do.call(cbind, mots)
if (progress) message("Calculating scores")
scores <- lapply_(COMPARE_METRICS, function(x) calc_scores(cols, x),
BP = BP, PB = progress)
names(scores) <- COMPARE_METRICS
if (progress) message("Calculating distribution parameters")
params <- lapply_(scores, function(x) get_params(x, rand.tries, widths),
BP = BP, PB = progress)
}
get_params <- function(scores, n, widths) {
means <- numeric(length(widths))
sds <- numeric(length(widths))
for (i in seq_along(widths)) {
ans <- vapply(seq_len(n),
function(x) sum(sample(scores, widths[i])),
numeric(1))
means[i] <- mean(ans)
sds[i] <- sd(ans)
}
names(means) <- widths
names(sds) <- widths
list(mean = means, sd = sds)
}
calc_scores <- function(cols, method) {
columns <- switch(method,
KL = columns + 0.01,
ALLR = columns + 0.01,
ALLR_LL = columns + 0.01,
columns
)
out <- numeric(ncol(columns)^2 / 2 - ncol(columns) / 2)
counter <- 0
for (i in seq_len(ncol(columns))) {
for (j in seq_len(ncol(columns))[-seq_len(i)]) {
counter <- counter + 1
out[counter] <- compare_columns(columns[, i], columns[, j], method)
}
}
out
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.