Package 'universalmotif'

Add multi-letter information to a motif.

Description

If the original sequences are available for a particular motif, then they can be used to generate higher-order PPM matrices. See the "Motif import, export, and manipulation" vignette for more information.

Usage

add_multifreq(motif, sequences, add.k = 2:3, RC = FALSE,
  threshold = 0.001, threshold.type = "pvalue", motifs.perseq = 1,
  add.bkg = FALSE)

Arguments

motif	See convert_motifs() for acceptable formats. If the motif is not a universalmotif motif, then it will be converted.
sequences	XStringSet The alphabet must match that of the motif. If these sequences are all the same length as the motif, then they are all used to generate the multi-freq matrices. Otherwise scan_sequences() is first run to find the best sequence stretches within these.
add.k	numeric(1) The k-let lengths to add.
RC	logical(1) If TRUE, check reverse complement of the input sequences. Only available for DNA/RNA.
threshold	numeric(1) See details.
threshold.type	character(1) One of c('pvalue', 'qvalue', 'logodds', 'logodds.abs'). See details.
motifs.perseq	numeric(1) If scan_sequences() is run, then this indicates how many hits from each sequence is to be used.
add.bkg	logical(1) Indicate whether to add corresponding higher order background information to the motif. Can sometimes be detrimental when the input consists of few short sequences, which can increase the likelihood of adding zero or near-zero probabilities.

Details

See scan_sequences() for more info on scanning parameters.

At each position in the motif, then the probability of each k-let covering from the initial position to ncol - 1 is calculated. Only positions within the motif are considered: this means that the final k-let probability matrix will have ncol - 1 fewer columns. Calculating k-let probabilities for the missing columns would be trivial however, as you would only need the background frequencies. Since these would not be useful for scan_sequences() though, they are not calculated.

Currently add_multifreq() does not try to stay faithful to the default motif matrix when generating multifreq matrices. This means that if the sequences used for training are completely different from the actual motif, the multifreq matrices will be as well. However this is only really a problem if you supply add_multifreq() with a set of sequences of the same length as the motif. In this case add_multifreq() is forced to create the multifreq matrices from these sequences. Otherwise add_multifreq() will scan the input sequences for the motif and use the best matches to construct the multifreq matrices.

This 'multifreq' representation is only really useful within the universalmotif environment. Despite this, if you wish it can be preserved in text using write_motifs().

A note on motif size

The number of rows for each k-let matrix is n^k, with n being the number of letters in the alphabet being used. This means that the size of the k-let matrix can become quite large as k increases. For example, if one were to wish to represent a DNA motif of length 10 as a 10-let, this would require a matrix with 1,048,576 rows (though at this point if what you want is to search for exact sequence matches, the motif format itself is not very useful).

Value

A universalmotif object with filled multifreq slot. The bkg slot is also expanded with corresponding higher order probabilities if add.bkg = TRUE.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

scan_sequences(), convert_motifs(), write_motifs()

Examples

sequences <- create_sequences(seqlen = 10)
motif <- create_motif()
motif.trained <- add_multifreq(motif, sequences, add.k = 2:4)
## peek at the 2-let matrix:
motif.trained["multifreq"]$`2`

---

Arabidopsis motif in universalmotif format.

Description

Arabidopsis motif trained from ArabidopsisPromoters using MEME version 4. This motif was generated at the command line using the following command: ⁠meme promoters.fa -revcomp -nmotifs 3 -mod anr -dna⁠.

Usage

ArabidopsisMotif

Format

universalmotif

---

Arabidopsis promoters as a DNAStringSet.

Description

50 Arabidopsis promoters, each 1000 bases long. See the "Sequence manipulation and scanning" vignette for an example workflow describing extracting promoter sequences.

Usage

ArabidopsisPromoters

Format

DNAStringSet

---

Compare motifs.

Description

Compare motifs using one of the several available metrics. See the "Motif comparisons and P-values" vignette for detailed information. Note: for regular DNA/RNA motif comparison using PCC, see compare_motifs2() for a faster and simpler alternative.

Usage

compare_motifs(motifs, compare.to, db.scores, use.freq = 1,
  use.type = "PPM", method = "PCC", tryRC = TRUE, min.overlap = 6,
  min.mean.ic = 0.25, min.position.ic = 0, relative_entropy = FALSE,
  normalise.scores = FALSE, max.p = 0.01, max.e = 10, nthreads = 1,
  score.strat = "a.mean", output.report, output.report.max.print = 10)

Arguments

motifs	See convert_motifs() for acceptable motif formats.
compare.to	numeric If missing, compares all motifs to all other motifs. Otherwise compares all motifs to the specified motif(s).
db.scores	data.frame or DataFrame. See details.
use.freq	numeric(1). For comparing the multifreq slot.
use.type	character(1) One of 'PPM' and 'ICM'. The latter allows for taking into account the background frequencies if relative_entropy = TRUE. Note that 'ICM' is not allowed when method = c("ALLR", "ALLR_LL").
method	character(1) One of PCC, EUCL, SW, KL, ALLR, BHAT, HELL, SEUCL, MAN, ALLR_LL, WEUCL, WPCC. See details.
tryRC	logical(1) Try the reverse complement of the motifs as well, report the best score.
min.overlap	numeric(1) Minimum overlap required when aligning the motifs. Setting this to a number higher then the width of the motifs will not allow any overhangs. Can also be a number between 0 and 1, representing the minimum fraction that the motifs must overlap.
min.mean.ic	numeric(1) Minimum mean information content between the two motifs for an alignment to be scored. This helps prevent scoring alignments between low information content regions of two motifs. Note that this can result in some comparisons failing if no alignment passes the mean IC threshold. Use average_ic() to filter out low IC motifs to get around this if you want to avoid getting NAs in your output.
min.position.ic	numeric(1) Minimum information content required between individual alignment positions for it to be counted in the final alignment score. It is recommended to use this together with normalise.scores = TRUE, as this will help punish scores resulting from only a fraction of an alignment.
relative_entropy	logical(1) Change the ICM calculation affecting min.position.ic and min.mean.ic. See convert_type().
normalise.scores	logical(1) Favour alignments which leave fewer unaligned positions, as well as alignments between motifs of similar length. Similarity scores are multiplied by the ratio of aligned positions to the total number of positions in the larger motif, and the inverse for distance scores.
max.p	numeric(1) Maximum P-value allowed in reporting matches. Only used if compare.to is set.
max.e	numeric(1) Maximum E-value allowed in reporting matches. Only used if compare.to is set. The E-value is the P-value multiplied by the number of input motifs times two.
nthreads	numeric(1) Run compare_motifs() in parallel with nthreads threads. nthreads = 0 uses all available threads.
score.strat	character(1) How to handle column scores calculated from motif alignments. "sum": add up all scores. "a.mean": take the arithmetic mean. "g.mean": take the geometric mean. "median": take the median. "wa.mean", "wg.mean": weighted arithmetic/geometric mean. "fzt": Fisher Z-transform. Weights are the total information content shared between aligned columns.
output.report	character(1) Provide a filename for compare_motifs() to write an html ouput report to. The top matches are shown alongside figures of the match alignments. This requires the knitr and rmarkdown packages. (Note: still in development.)
output.report.max.print	numeric(1) Maximum number of top matches to print.

Details

Available metrics

The following metrics are available:

• Euclidean distance (EUCL) (Choi et al. 2004)

• Weighted Euclidean distance (WEUCL)

• Kullback-Leibler divergence (KL) (Kullback and Leibler 1951; Roepcke et al. 2005)

• Hellinger distance (HELL) (Hellinger 1909)

• 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 (Sandelin and Wasserman 2004)

• Average log-likelihood ratio (ALLR) (Wang and Stormo 2003)

• Lower limit ALLR (ALLR_LL) (Mahony et al. 2007)

• Bhattacharyya coefficient (BHAT) (Bhattacharyya 1943)

Comparisons are calculated between two motifs at a time. All possible alignments are scored, and the best score is reported. In an alignment scores are calculated individually between columns. How those scores are combined to generate the final alignment scores depends on score.strat.

See the "Motif comparisons and P-values" vignette for a description of the various metrics. Note that PCC, WPCC, SW, ALLR, ALLR_LL and BHAT are similarities; higher values mean more similar motifs. For the remaining metrics, values closer to zero represent more similar motifs.

Small pseudocounts are automatically added when one of the following methods is used: KL, ALLR, ALLR_LL, IS. This is avoid zeros in the calculations.

Calculating P-values

To note regarding p-values: P-values are pre-computed using the make_DBscores() function. If not given, then uses a set of internal precomputed P-values from the JASPAR2018 CORE motifs. These precalculated scores are dependent on the length of the motifs being compared. This takes into account that comparing small motifs with larger motifs leads to higher scores, since the probability of finding a higher scoring alignment is higher.

The default P-values have been precalculated for regular DNA motifs. They are of little use for motifs with a different number of alphabet letters (or even the multifreq slot).

Value

matrix if compare.to is missing; DataFrame otherwise. For the latter, function args are stored in the metadata slot.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Bhattacharyya A (1943). “On a measure of divergence between two statistical populations defined by their probability distributions.” Bulletin of the Calcutta Mathematical Society, 35, 99-109.

Choi I, Kwon J, Kim S (2004). “Local feature frequency profile: a method to measure structural similarity in proteins.” PNAS, 101, 3797-3802.

Hellinger E (1909). “Neue Begrundung der Theorie quadratischer Formen von unendlichvielen Veranderlichen.” Journal fur die reine und angewandte Mathematik, 136, 210-271.

Khan A, Fornes O, Stigliani A, Gheorghe M, Castro-Mondragon JA, van der Lee R, Bessy A, Cheneby J, Kulkarni SR, Tan G, Baranasic D, Arenillas DJ, Sandelin A, Vandepoele K, Lenhard B, Ballester B, Wasserman WW, Parcy F, Mathelier A (2018). “JASPAR 2018: update of the open-access database of transcription factor binding profiles and its web framework.” Nucleic Acids Research, 46, D260-D266.

Kullback S, Leibler RA (1951). “On information and sufficiency.” The Annals of Mathematical Statistics, 22, 79-86.

Itakura F, Saito S (1968). “Analysis synthesis telephony based on the maximum likelihood method.” In 6th International Congress on Acoustics, C-17.

Mahony S, Auron PE, Benos PV (2007). “DNA Familial Binding Profiles Made Easy: Comparison of Various Motif Alignment and Clustering Strategies.” PLoS Computational Biology, 3.

Pietrokovski S (1996). “Searching databases of conserved sequence regions by aligning protein multiple-alignments.” Nucleic Acids Research, 24, 3836-3845.

Roepcke S, Grossmann S, Rahmann S, Vingron M (2005). “T-Reg Comparator: an analysis tool for the comparison of position weight matrices.” Nucleic Acids Research, 33, W438-W441.

Sandelin A, Wasserman WW (2004). “Constrained binding site diversity within families of transcription factors enhances pattern discovery bioinformatics.” Journal of Molecular Biology, 338, 207-215.

Wang T, Stormo GD (2003). “Combining phylogenetic data with co-regulated genes to identify motifs.” Bioinformatics, 19, 2369-2380.

See Also

convert_motifs(), motif_tree(), view_motifs(), make_DBscores(), compare_motifs2(), view_motifs2()

Examples

motif1 <- create_motif(name = "1")
motif2 <- create_motif(name = "2")
motif1vs2 <- compare_motifs(c(motif1, motif2), method = "PCC")
## To get a dist object:
as.dist(1 - motif1vs2)

motif3 <- create_motif(name = "3")
motif4 <- create_motif(name = "4")
motifs <- c(motif1, motif2, motif3, motif4)
## Compare motif "2" to all the other motifs:
if (R.Version()$arch != "i386") {
compare_motifs(motifs, compare.to = 2, max.p = 1, max.e = Inf)
}

## If you are working with a large list of motifs and the mean.min.ic
## option is not set to zero, you may get a number of failed comparisons
## due to low IC. To filter the list of motifs to avoid these, use
## the average_ic() function to remove motifs with low average IC:
## Not run: 
library(MotifDb)
motifs <- convert_motifs(MotifDb)[1:100]
compare_motifs(motifs)
#> Warning in compare_motifs(motifs) :
#>   Some comparisons failed due to low IC
motifs <- motifs[average_ic(motifs) > 0.5]
compare_motifs(motifs)

## End(Not run)

---

Faster minimalist motif comparison.

Description

compare_motifs2() is a deliberately pared-down counterpart to compare_motifs(), with defaults and an algorithm that mirror those of the command-line tool yamtk. It uses the per-column Pearson correlation as its similarity metric, computes significance from a TOMTOM-style null PMF (either empirical, from the target database, or parametric, via a Dirichlet-Multinomial over a K=5 simplex grid), applies both Bonferroni and Benjamini-Hochberg adjustment, and returns either a long-format data.frame of the significant pairs or a square score / p-value / q-value matrix.

Usage

compare_motifs2(motifs, compare.to = NULL, qvalue = 0.1,
  min.overlap = 5L, RC = TRUE, null = c("empirical", "parametric"),
  bkg = NULL, matrix.out = c("score", "pvalue", "qvalue"), nthreads = 1)

Arguments

motifs	See convert_motifs() for accepted motif formats. DNA or RNA only.
compare.to	integer, character, or NULL. When NULL (default), all-pairs mode: every motif is compared against every other and a square matrix is returned (see matrix.out). When integer/character, the named/indexed motifs are treated as queries and the entire input is treated as the target database; a long-format data.frame of significant matches is returned.
qvalue	numeric(1). q-value cutoff for the long-format output. Default 0.1, matching ⁠yamtk cmp⁠. Ignored in matrix mode.
min.overlap	integer(1). Minimum overlap length (in columns) for a candidate alignment. Default 5.
RC	logical(1). If TRUE (default), also score query motifs in reverse-complement orientation and report the best of the two.
null	character(1). Null-PMF mode: "empirical" (default) or "parametric". See Details.
bkg	numeric(4) or NULL. Background base frequencies (A, C, G, T). If NULL, uniform c(.25, .25, .25, .25) is used. Critical for "parametric" (sets the Dirichlet prior); affects only pseudocount smoothing for "empirical".
matrix.out	character(1). In matrix mode, which metric to return: "score" (mean PCC; default), "pvalue", or "qvalue". Ignored in long-format mode.
nthreads	numeric(1). Number of threads. nthreads = 0 uses all available threads.

Details

Use compare_motifs() when you need any of the broader features: alternative similarity metrics (EUCL, ALLR, KL, …), score-aggregation strategies, IC-based filters, multifreq scoring, ICM/PPM dispatch, relative-entropy mode, or the JASPAR-derived db.scores lookup-table p-value path.

P-values are computed in two modes mirroring ⁠yamtk cmp⁠:

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 is O(N x T) where N is the number of target motifs and T their mean width. 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 * 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.

The reported score is the mean PCC over the overlap (always in $[-1, 1]$), not the raw sum.

Value

• Matrix mode (compare.to = NULL): a square numeric matrix indexed by motif name, holding the value selected by matrix.out.

• Long-format mode (compare.to provided): a data.frame with columns subject, subject.i, target, target.i, offset, strand, overlap, score, subject.consensus, target.consensus, pvalue, qvalue, filtered to qvalue <= qvalue and sorted by q-value ascending. Coordinates are 1-based; offset is the query column where target column 1 aligns.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Gupta S, Stamatoyannopoulos JA, Bailey TL, Noble WS (2007). "Quantifying similarity between motifs." Genome Biology, 8(2), R24.

Tremblay BJM (2026). yamtk: Yet Another Motif ToolKit. https://github.com/bjmt/yamtk.

See Also

compare_motifs(), scan_sequences2(), view_motifs2()

Examples

library(universalmotif)
set.seed(1)
motifs <- lapply(1:6, function(i)
  create_motif(paste(sample(c("A","C","G","T"), 8, replace = TRUE),
                     collapse = ""), name = paste0("M", i)))
## Matrix of mean PCC scores
m <- compare_motifs2(motifs)
round(m, 2)
## Long-format hits at q <= 0.5 with motif 1 as the query
compare_motifs2(motifs, compare.to = 1, qvalue = 0.5)

---

Convert motif class.

Description

Allows for easy transfer of motif information between different classes as defined by other Bioconductor packages. This function is also used by nearly all other functions in this package, so any motifs of a compatible class can be used without needing to be converted beforehand.

Usage

convert_motifs(motifs, class = "universalmotif-universalmotif")

## S4 method for signature 'AsIs'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'list'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'universalmotif'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'MotifList'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'TFFMFirst'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'PFMatrix'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'PWMatrix'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'ICMatrix'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'XMatrixList'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'pwm'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'pcm'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'pfm'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'PWM'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'Motif'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

## S4 method for signature 'matrix'
convert_motifs(motifs,
  class = "universalmotif-universalmotif")

Arguments

motifs	Single motif object or list. See details.
class	character(1) Desired motif class. Input as 'package-class'. If left empty, defaults to 'universalmotif-universalmotif'. (See details.)

Details

Input

The following packge-class combinations can be used as input:

• MotifDb-MotifList

• TFBSTools-PFMatrix

• TFBSTools-PWMatrix

• TFBSTools-ICMatrix

• TFBSTools-PFMatrixList

• TFBSTools-PWMatrixList

• TFBSTools-ICMatrixList

• TFBSTools-TFFMFirst

• seqLogo-pwm

• motifStack-pcm

• motifStack-pfm

• PWMEnrich-PWM

• motifRG-Motif

• universalmotif-universalmotif

• matrix

Output

The following package-class combinations can be output:

• MotifDb-MotifList

• TFBSTools-PFMatrix

• TFBSTools-PWMatrix

• TFBSTools-ICMatrix

• TFBSTools-TFFMFirst

• seqLogo-pwm

• motifStack-pcm

• motifStack-pfm

• PWMEnrich-PWM

• Biostrings-PWM (type = 'log2prob')

• rGADEM-motif

• universalmotif-universalmotif (the default, no need to specify this)

Note: MotifDb-MotifList output was a later addition to convert_motifs(). As a result, to stay consistent with previous behaviour most functions will always convert MotifDb-MotifList objects to a list of universalmotif motifs, even if other formats would be simply returned as is (e.g. for other formats, filter_motifs() will return the input format; for MotifDb-MotifList, a list of universalmotif objects will be returned).

Value

Single motif object or list.

Methods (by class)

• convert_motifs(AsIs): Generate an error to remind users to run to_list() instead of using the column from to_df() directly.

• convert_motifs(list): Convert a list of motifs.

• convert_motifs(universalmotif): Convert a universalmotif object.

• convert_motifs(MotifList): Convert MotifList motifs. (MotifDb)

• convert_motifs(TFFMFirst): Convert TFFMFirst motifs. (TFBSTools)

• convert_motifs(PFMatrix): Convert PFMatrix motifs. (TFBSTools)

• convert_motifs(PWMatrix): Convert PWMatrix motifs. (TFBSTools)

• convert_motifs(ICMatrix): Convert ICMatrix motifs. (TFBSTools)

• convert_motifs(XMatrixList): Convert XMatrixList motifs. (TFBSTools)

• convert_motifs(pwm): Convert pwm motifs. (seqLogo)

• convert_motifs(pcm): Convert pcm motifs. (motifStack)

• convert_motifs(pfm): Convert pfm motifs. (motifStack)

• convert_motifs(PWM): Convert PWM motifs. (PWMEnrich)

• convert_motifs(Motif): Convert Motif motifs. (motifRG)

• convert_motifs(matrix): Create motif from matrices.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Bembom O (2018). seqLogo: Sequence logos for DNA sequence alignments. R package version 1.46.0.

Droit A, Gottardo R, Robertson G, Li L (2014). rGADEM: de novo motif discovery. R package version 2.28.0.

Mercier E, Gottardo R (2014). MotIV: Motif Identification and Validation. R package version 1.36.0.

Ou J, Wolfe SA, Brodsky MH, Zhu LJ (2018). “motifStack for the analysis of transcription factor binding site evolution.” Nature Methods, 15, 8-9. doi: 10.1038/nmeth.4555.

Shannon P, Richards M (2018). MotifDb: An Annotated Collection of Protein-DNA Binding Sequence Motifs. R package version 1.22.0.

Stojnic R, Diez D (2015). PWMEnrich: PWM enrichment analysis. R package version 4.16.0.

Tan G, Lenhard B (2016). “TFBSTools: an R/Bioconductor package for transcription factor binding site analysis.” Bioinformatics, 32, 1555-1556. doi: 10.1093/bioinformatics/btw024.

Yao Z (2012). motifRG: A package for discriminative motif discovery, designed for high throughput sequencing dataset. R package version 1.24.0.

Examples

# Convert from universalmotif:
jaspar <- read_jaspar(system.file("extdata", "jaspar.txt",
                                  package = "universalmotif"))
if (requireNamespace("motifStack", quietly = TRUE)) {
  jaspar.motifstack.pfm <- convert_motifs(jaspar, "motifStack-pfm")
}

# Convert from another class to universalmotif:
if (requireNamespace("TFBSTools", quietly = TRUE)) {
library(TFBSTools)
data(MA0003.2)
motif <- convert_motifs(MA0003.2)

# Convert from another class to another class
if (requireNamespace("PWMEnrich", quietly = TRUE)) {
  motif <- convert_motifs(MA0003.2, "PWMEnrich-PWM")
}

# The 'convert_motifs' function is embedded in the rest of the universalmotif
# functions: non-universalmotif class motifs can be used
MA0003.2.trimmed <- trim_motifs(MA0003.2)
# Note: if the motif object going in has information that the
# 'universalmotif' class can't hold, it will be lost
}

---

Convert universalmotif type.

Description

Switch between position count matrix (PCM), position probability matrix (PPM), position weight matrix (PWM), information count matrix (ICM), and contribution weight matrix (CWM) types. See the "Introduction to sequence motifs" vignette for details. Please also note that type conversion occurs implicitly throughout the universalmotif package, so there is generally no need to perform this manual conversion. Also please be aware that the message concerning pseudocount-adjusting motifs can be disabled via options(pseudocount.warning=FALSE).

Usage

convert_type(motifs, type, pseudocount, nsize_correction = FALSE,
  relative_entropy = FALSE)

Arguments

motifs	See convert_motifs() for acceptable formats.
type	character(1) One of c('PCM', 'PPM', 'PWM', 'ICM'). 'CWM' is only valid as the source type of a conversion; convert_type(other, "CWM") errors because contribution scores cannot be recovered from a probabilistic matrix.
pseudocount	numeric(1) Correction to be applied to prevent -Inf from appearing in PWM matrices. If missing, the pseudocount stored in the universalmotif 'pseudocount' slot will be used.
nsize_correction	logical(1) If true, the ICM at each position will be corrected to account for small sample sizes. Only used if relative_entropy = FALSE.
relative_entropy	logical(1) If true, the ICM will be calculated as relative entropy. See details.

Details

PCM

Position count matrix (PCM), also known as position frequency matrix (PFM). For n sequences from which the motif was built, each position is represented by the numbers of each letter at that position. In theory all positions should have sums equal to n, but not all databases are this consistent. If converting from another type to PCM, column sums will be equal to the 'nsites' slot. If empty, 100 is used.

PPM

Position probability matrix (PPM), also known as position frequency matrix (PFM). At each position, the probability of individual letters is calculated by dividing the count for that letter by the total sum of counts at that position (letter_count / position_total). As a result, each position will sum to 1. Letters with counts of 0 will thus have a probability of 0, which can be undesirable when searching for motifs in a set of sequences. To avoid this a pseudocount can be added ((letter_count + pseudocount) / (position_total + pseudocount)).

PWM

Position weight matrix (PWM; Stormo et al. (1982)), also known as position-specific weight matrix (PSWM), position-specific scoring matrix (PSSM), or log-odds matrix. At each position, each letter is represented by it's log-likelihood (log2(letter_probability / background_probility)), which is normalized using the background letter frequencies. A PWM matrix is constructed from a PPM. If any position has 0-probability letters to which pseudocounts were not added, then the final log-likelihood of these letters will be -Inf.

ICM

Information content matrix (ICM; Schneider and Stephens 1990). An ICM is a PPM where each letter probability is multiplied by the total information content at that position. The information content of each position is determined as: totalIC - Hi, where the total information totalIC

totalIC <- log2(alphabet_length), and the Shannon entropy (Shannon 1948) for a specific position (Hi)

⁠Hi <- -sum(sapply(alphabet_frequencies, function(x) x * log(2))⁠.

As a result, the total sum or height of each position is representative of it's sequence conservation, measured in the unit 'bits', which is a unit of energy (Schneider 1991; see https://fr-s-schneider.ncifcrf.gov/logorecommendations.html for more information). However not all programs will calculate information content the same. Some will 'correct' the total information content at each position using a correction factor as described by Schneider et al. (1986). This correction can applied by setting nsize_correction = TRUE, however it will only be applied if the 'nsites' slot is not empty. This is done using TFBSTools:::schneider_correction (Tan and Lenhard 2016). As such, converting from an ICM to which some form of correction has been applied will result in a PCM/PPM/PWM with slight inaccuracies.

Another method of calculating information content is calculating the relative entropy, also known as Kullback-Leibler divergence (Kullback and Leibler 1951). This accounts for background frequencies, which can be useful for genomes with a heavy imbalance in letter frequencies. For each position, the individual letter frequencies are calculated as letter_freq * log2(letter_freq / bkg_freq). When calculating information content using Shannon entropy, the maximum content for each position will always be log2(alphabet_length). This does not hold for information content calculated as relative entropy. Please note that conversion from ICM assumes the information content was not calculated as relative entropy.

CWM

Contribution weight matrix (CWM), as produced by attribution methods such as TF-MoDISco (Shrikumar et al. 2018) and DeepLIFT (Shrikumar, Greenside, and Kundaje 2017). Each cell holds a signed real-valued contribution score: how much the given letter at the given position drives a model's prediction. Unlike the four probabilistic types above, CWMs are not normalised (columns do not sum to 1 or to nsites), values may be negative (the model penalises that letter at that position), and there is no canonical maximum height. CWM is therefore only ever a source type for convert_type(): convert_type(cwm, "PPM") column-normalises absolute values via the TF-MoDISco-lite convention ⁠ppm[i,j] = |cwm[i,j]| / sum_i |cwm[i,j]|⁠ (columns whose absolute column sum is zero collapse to a uniform distribution). convert_type(cwm, "PWM" | "ICM" | "PCM") routes through CWM -> PPM first and then applies the standard ⁠PPM -> *⁠ step, so CWM -> PWM produces a real log-odds matrix (log2(ppm / bkg)), not a re-labelled CWM. The reverse direction (PPM, PWM, ICM, or PCM -> CWM) errors with a clear message: contribution scores cannot be inferred from a probability distribution, so a CWM must originate from an attribution method (read via read_meme() with CWM = TRUE, or built with create_motif(matrix, type = "CWM")).

Value

See convert_motifs() for possible output motif objects.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Kullback S, Leibler RA (1951). “On information and sufficiency.” The Annals of Mathematical Statistics, 22, 79-86.

Nishida K, Frith MC, Nakai K (2009). “Pseudocounts for transcription factor binding sites.” Nucleic Acids Research, 37, 939-944.

Schneider TD, Stormo GD, Gold L, Ehrenfeucht A (1986). “Information content of binding sites on nucleotide sequences.” Journal of Molecular Biology, 188, 415-431.

Schneider TD, Stephens RM (1990). “Sequence Logos: A New Way to Display Consensus Sequences.” Nucleic Acids Research, 18, 6097-6100.

Schneider TD (1991). “Theory of Molecular Machines. II. Energy Dissipation from Molecular Machines.” Journal of Theoretical Biology, 148, 125-137.

Shannon CE (1948). “A Mathematical Theory of Communication.” Bell System Technical Journal, 27, 379-423.

Shrikumar A, Greenside P, Kundaje A (2017). “Learning important features through propagating activation differences.” Proceedings of the 34th International Conference on Machine Learning, 70, 3145-3153.

Shrikumar A, Tian K, Avsec Ž, Shcherbina A, Banerjee A, Sharmin M, Nair S, Kundaje A (2018). “TF-MoDISco v0.4.4.2-alpha: Technical Note.” arXiv:1811.00416.

Stormo GD, Schneider TD, Gold L, Ehrenfeucht A (1982). “Use of the Perceptron algorithm to distinguish translational initiation sites in E. coli.” Nucleic Acids Research, 10, 2997-3011.

Tan G, Lenhard B (2016). “TFBSTools: an R/Bioconductor package for transcription factor binding site analysis.” Bioinformatics, 32, 1555-1556. doi: 10.1093/bioinformatics/btw024.

See Also

convert_motifs()

Examples

jaspar.pcm <- read_jaspar(system.file("extdata", "jaspar.txt",
                                      package = "universalmotif"))

## The motifs pseudocounts are 1: these will be used in the PCM->PPM
## calculation
jaspar.pwm <- convert_type(jaspar.pcm, type = "PPM")

## Setting pseudocount to 0 will prevent any correction from being
## applied to PPM/PWM matrices, overriding the motifs own pseudocounts
jaspar.pwm <- convert_type(jaspar.pcm, type = "PWM", pseudocount = 0)

## CWM -> PPM via |cwm| / sum (the TF-MoDISco convention):
cwm.mat <- matrix(c( 0.10, -0.20,  0.80, -0.10,
                    -0.30,  0.90, -0.20,  0.10,
                     0.70, -0.10, -0.10, -0.40,
                    -0.20, -0.10, -0.10,  0.90),
                  nrow = 4, byrow = FALSE,
                  dimnames = list(c("A","C","G","T"), NULL))
cwm <- create_motif(cwm.mat, type = "CWM", name = "modisco_example")
convert_type(cwm, "PPM")

---

Create a motif.

Description

Create a motif from a set of sequences, a matrix, or generate a random motif. See the "Motif import, export and manipulation" vignette for details.

Usage

create_motif(input, alphabet, type = "PPM", name = "motif",
  pseudocount = 0, bkg, nsites, altname, family, organism, bkgsites, strand,
  pval, qval, eval, extrainfo, add.multifreq)

## S4 method for signature 'missing'
create_motif(input, alphabet, type = "PPM",
  name = "motif", pseudocount = 0, bkg, nsites, altname, family, organism,
  bkgsites, strand, pval, qval, eval, extrainfo, add.multifreq)

## S4 method for signature 'numeric'
create_motif(input, alphabet, type = "PPM",
  name = "motif", pseudocount = 0, bkg, nsites, altname, family, organism,
  bkgsites, strand, pval, qval, eval, extrainfo, add.multifreq)

## S4 method for signature 'character'
create_motif(input, alphabet, type = "PPM",
  name = "motif", pseudocount = 0, bkg, nsites, altname, family, organism,
  bkgsites, strand, pval, qval, eval, extrainfo, add.multifreq)

## S4 method for signature 'matrix'
create_motif(input, alphabet, type = "PPM",
  name = "motif", pseudocount = 0, bkg, nsites, altname, family, organism,
  bkgsites, strand, pval, qval, eval, extrainfo, add.multifreq)

## S4 method for signature 'DNAStringSet'
create_motif(input, alphabet, type = "PPM",
  name = "motif", pseudocount = 0, bkg, nsites, altname, family, organism,
  bkgsites, strand, pval, qval, eval, extrainfo, add.multifreq)

## S4 method for signature 'RNAStringSet'
create_motif(input, alphabet, type = "PPM",
  name = "motif", pseudocount = 0, bkg, nsites, altname, family, organism,
  bkgsites, strand, pval, qval, eval, extrainfo, add.multifreq)

## S4 method for signature 'AAStringSet'
create_motif(input, alphabet, type = "PPM",
  name = "motif", pseudocount = 0, bkg, nsites, altname, family, organism,
  bkgsites, strand, pval, qval, eval, extrainfo, add.multifreq)

## S4 method for signature 'BStringSet'
create_motif(input, alphabet, type = "PPM",
  name = "motif", pseudocount = 0, bkg, nsites, altname, family, organism,
  bkgsites, strand, pval, qval, eval, extrainfo, add.multifreq)

Arguments

input	character, numeric, matrix, XStringSet, or missing.
alphabet	character(1) One of c('DNA', 'RNA', 'AA'), or a combined string representing the letters. If no alphabet is provided then it will try and guess the alphabet from the input.
type	character(1) One of c('PCM', 'PPM', 'PWM', 'ICM', 'CWM'). 'CWM' is intended for matrix input only. See convert_type() for the CWM details and the available conversions.
name	character(1) Motif name.
pseudocount	numeric(1) Correction to be applied to prevent -Inf from appearing in PWM matrices. Defaults to 0.
bkg	numeric A vector of probabilities, each between 0 and 1. If higher order backgrounds are provided, then the elements of the vector must be named. If unnamed, then the order of probabilities must be in the same order as the alphabetically sorted sequence alphabet.
nsites	numeric(1) Number of sites the motif was constructed from. If blank, then create_motif() will guess the appropriate number if possible. To prevent this, provide nsites = numeric().
altname	character(1) Alternate motif name.
family	character(1) Transcription factor family.
organism	character(1) Species of origin.
bkgsites	numeric(1) Total number of sites used to find the motif.
strand	character(1) Whether the motif is specific to a certain strand. Acceptable strands are '+', '-', and '+-' (to represent both strands). Note that '-+' and '*' can also be provided to represent both strands, but the final strand in the universalmotif object will be set to '+-'.
pval	numeric(1) P-value associated with motif.
qval	numeric(1) Adjusted P-value associated with motif.
eval	numeric(1) E-value associated with motif.
extrainfo	character Any other extra information, represented as a named character vector.
add.multifreq	numeric If the motif is created from a set of sequences, then the add_multifreq() function can be run at the same time (with RC = FALSE).

Details

The aim of this function is provide an easy interface to creating universalmotif motifs, as an alternative to the default class constructor (i.e. new('universalmotif', name=...)). See examples for potential use cases.

Note: when generating random motifs, the nsites slot is also given a random value.

Note: type = "CWM" is only supported when input is a matrix. CWMs originate from contribution-attribution methods (TF-MoDISco, DeepLIFT) and cannot be derived from sequences, consensus strings, or random generation. The input matrix is stored as-is (signed real values, no column-sum constraint); no auto-normalisation is applied. See convert_type() for the full CWM semantics.

See the examples section for more info on motif creation.

Value

universalmotif object.

Methods (by class)

• create_motif(missing): Create a random motif of length 10.

• create_motif(numeric): Create a random motif with a specified length.

• create_motif(character): Create motif from a consensus string.

• create_motif(matrix): Create motif from a matrix.

• create_motif(DNAStringSet): Create motif from a DNAStringSet.

• create_motif(RNAStringSet): Create motif from a RNAStringSet.

• create_motif(AAStringSet): Create motif from a AAStringSet.

• create_motif(BStringSet): Create motif from a BStringSet.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

convert_type(), add_multifreq(), create_sequences(), shuffle_motifs().

Examples

##### create motifs from a single string

# Motif is by default generated as a PPM: change final type as desired
DNA.motif <- create_motif("TATAWAW")
DNA.motif <- create_motif("TATAWAW", type = "PCM")

# Nsites will be set to the number of input sequences unless specified or
# a single string is used as input
DNA.motif <- create_motif("TTTTTTT", nsites = 10)

# Ambiguity letters can be used:
DNA.motif <- create_motif("TATAWAW")
DNA.motif <- create_motif("NNVVWWAAWWDDN")

# Be careful about setting nsites when using ambiguity letters!
DNA.motif <- create_motif("NNVVWWAAWWDDN", nsites = 1)

RNA.motif <- create_motif("UUUCCG")

# 'create_motif' will try to detect the alphabet type; this can be
# unreliable for AA and custom alphabets as DNA and RNA alphabets are
# detected first
AA.motif <- create_motif("AVLK", alphabet = "AA")

custom.motif <- create_motif("QWER", alphabet = "QWER")
# Specify custom alphabet
custom.motif <- create_motif("QWER", alphabet = "QWERASDF")

###### Create motifs from multiple strings of equal length

DNA.motif <- create_motif(c("TTTT", "AAAA", "AACC", "TTGG"), type = "PPM")
DNA.motif <- create_motif(c("TTTT", "AAAA", "AACC", "TTGG"), nsites = 20)
RNA.motif <- create_motif(c("UUUU", "AAAA", "AACC", "UUGG"), type = "PWM")
AA.motif <- create_motif(c("ARNDCQ", "EGHILK", "ARNDCQ"), alphabet = "AA")
custom.motif <- create_motif(c("POIU", "LKJH", "POIU", "CVBN"),
                             alphabet = "POIULKJHCVBN")

# Ambiguity letters are only allowed for single consensus strings: the
# following fails
## Not run: 
create_motif(c("WWTT", "CCGG"))
create_motif(c("XXXX", "XXXX"), alphabet = "AA")

## End(Not run)

##### Create motifs from XStringSet objects

library(Biostrings)

DNA.set <- DNAStringSet(c("TTTT", "AAAA", "AACC", "TTGG"))
DNA.motif <- create_motif(DNA.set)
RNA.set <- RNAStringSet(c("UUUU", "AACC", "UUCC"))
RNA.motif <- create_motif(RNA.set)
AA.set <- AAStringSet(c("VVVLLL", "AAAIII"))
AA.motif <- create_motif(AA.set)

# Custom motifs can be created from BStringSet objects
B.set <- BStringSet(c("QWER", "ASDF", "ZXCV", "TYUI"))
custom.motif <- create_motif(B.set)

##### Create motifs with filled 'multifreq' slot

DNA.motif.k2 <- create_motif(DNA.set, add.multifreq = 2)

##### Create motifs from matrices

mat <- matrix(c(1, 1, 1, 1,
                2, 0, 2, 0,
                0, 2, 0, 2,
                0, 0, 0, 0),
                nrow = 4, byrow = TRUE)
DNA.motif <- create_motif(mat, alphabet = "DNA")
RNA.motif <- create_motif(mat, alphabet = "RNA", nsites = 20)
custom.motif <- create_motif(mat, alphabet = "QWER")

# Specify custom alphabet
custom.motif <- create_motif(mat, alphabet = "QWER")

# Alphabet can be detected from rownames
rownames(mat) <- DNA_BASES
DNA.motif <- create_motif(mat)
rownames(mat) <- c("Q", "W", "E", "R")
custom.motif <- create_motif(mat)

# Matrices can also be used as input
mat.ppm <- matrix(c(0.1, 0.1, 0.1, 0.1,
                    0.5, 0.5, 0.5, 0.5,
                    0.1, 0.1, 0.1, 0.1,
                    0.3, 0.3, 0.3, 0.3),
                    nrow = 4, byrow = TRUE)

DNA.motif <- create_motif(mat.ppm, alphabet = "DNA", type = "PPM")

# Contribution weight matrices (signed, no column-sum constraint).
# Typically these come from TF-MoDISco / DeepLIFT pipelines and are
# imported via `read_meme(..., CWM = TRUE)`, but can be built from
# a matrix directly:
cwm.mat <- matrix(c( 0.10, -0.20,  0.80, -0.10,
                    -0.30,  0.90, -0.20,  0.10,
                     0.70, -0.10, -0.10, -0.40,
                    -0.20, -0.10, -0.10,  0.90),
                  nrow = 4, byrow = FALSE,
                  dimnames = list(c("A","C","G","T"), NULL))
CWM.motif <- create_motif(cwm.mat, type = "CWM", name = "modisco_example")

##### Create random motifs

# These are generated as PPMs with 10 positions

DNA.motif <- create_motif()
RNA.motif <- create_motif(alphabet = "RNA")
AA.motif <- create_motif(alphabet = "AA")
custom.motif <- create_motif(alphabet = "QWER")

# The number of positions can be specified

DNA.motif <- create_motif(5)

# If the background frequencies are not provided, they are generated
# using `rpois`; positions are created using `rdirichlet(1, bkg)`.
# (calling `create_motif()` creates motifs with an average
# positional IC of 1)

DNA.motif <- create_motif(bkg = c(0.3, 0.2, 0.2, 0.3))
DNA.motif <- create_motif(10, bkg = c(0.1, 0.4, 0.4, 0.1))

---

Create random sequences.

Description

Generate random sequences from any set of characters, represented as XStringSet objects.

Usage

create_sequences(alphabet = "DNA", seqnum = 100, seqlen = 100, freqs,
  nthreads = 1, rng.seed = sample.int(10000, 1))

Arguments

alphabet	character(1) One of c('DNA', 'RNA', 'AA'), or a string of characters to be used as the alphabet.
seqnum	numeric(1) Number of sequences to generate.
seqlen	numeric(1) Length of random sequences.
freqs	numeric A named vector of probabilities. The length of the vector must be the power of the number of letters in the sequence alphabet. Probabilities can only be provided for a single size k.
nthreads	numeric(1) Run create_sequences() in parallel with nthreads threads. nthreads = 0 uses all available threads. Note that no speed up will occur for jobs with seqnum = 1.
rng.seed	numeric(1) Set random number generator seed. Since sequence creation can occur simultaneously in multiple threads using C++, it cannot communicate with the regular R random number generator state and thus requires an independent seed. Each individual sequence creation instance is given the following seed: rng.seed * index. The default is to pick a random number as chosen by sample(), which effectively is making create_sequences() dependent on the R RNG state.

Value

XStringSet The returned sequences are unnamed.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

create_motif(), shuffle_sequences()

Examples

## Create DNA sequences with slightly increased AT content:
sequences <- create_sequences(freqs = c(A=0.3, C=0.2, G=0.2, T=0.3))
## Create custom sequences:
sequences.QWER <- create_sequences("QWER")
## You can include non-alphabet characters are well, even spaces:
sequences.custom <- create_sequences("!@#$ ")

---

Deduplicate overlapping motif hits.

Description

Given a set of motif hits, typically the output of scan_sequences() or scan_sequences2(), or a GRanges / data.frame from any other source, drop overlapping hits within each ⁠(sequence, motif, strand)⁠ group using a greedy algorithm: cluster overlapping hits, then within each cluster keep the hit with the best priority and, after that, any other hit that does not overlap an already-kept hit. This is the same algorithm implemented by yamtk.

Usage

dedup_hits(hits, by = "pvalue", reverse = FALSE, ignore.strand = FALSE,
  ignore.motif = FALSE, seq = "sequence", motif = "motif",
  strand = "strand")

Arguments

hits	data.frame or GRanges. For a data.frame, must have columns start, end, and the columns named by by, seq, motif, strand. For a GRanges, the seqnames/start/end/ strand slots are used; the motif and priority columns are read from mcols().
by	character(1). Name of the priority column. The semantic is auto-detected: a column named pvalue (or anything starting with pval) is treated as lower-is-better; anything else (e.g. score) is treated as higher-is-better. Use reverse = TRUE to flip.
reverse	logical(1). Flip the priority direction.
ignore.strand	logical(1). If TRUE, + and - hits compete with each other at overlapping coordinates.
ignore.motif	logical(1). If TRUE, hits from different motifs compete with each other at overlapping coordinates.
seq	character(1). Name of the sequence-id column in a data.frame input. Default "sequence". Ignored for GRanges.
motif	character(1). Name of the motif-id column in a data.frame input. Default "motif". Ignored for GRanges.
strand	character(1). Name of the strand column in a data.frame input. Default "strand". Ignored for GRanges.

Details

Note that this differs from scan_sequences()'s built-in no.overlaps argument, which collapses every connected-overlap cluster to a single survivor (hierarchical clustering with cutree). The greedy approach used here keeps legitimately distinct hits even when an intermediate noisy hit bridges them, which is usually what users want.

Coordinate convention: rows must have start <= end. Overlaps are inclusive: two hits overlap if a$start <= b$end & b$start <= a$end.

Value

The input hits, subsetted to the kept rows in the original order.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Tremblay BJM (2026). yamtk: Yet Another Motif ToolKit. https://github.com/bjmt/yamtk.

See Also

scan_sequences2(), scan_sequences()

Examples

library(universalmotif)
motifs <- list(create_motif("TATAAA", name = "M1"),
               create_motif("CACGTG", name = "M2"))
seqs <- create_sequences(seqnum = 5, seqlen = 400, rng.seed = 1)
hits <- scan_sequences2(motifs, seqs, pvalue = 5e-2,
                        return.granges = FALSE)
deduped <- dedup_hits(hits)
nrow(hits)
nrow(deduped)

---

Enrich for input motifs in a set of sequences.

Description

Given a set of target and background sequences, test if the input motifs are significantly enriched in the targets sequences relative to the background sequences. See the "Sequence manipulation and scanning" vignette. Note: for regular DNA/RNA motif enrichment, see enrich_motifs2() for a faster and simpler alternative.

Usage

enrich_motifs(motifs, sequences, bkg.sequences, max.p = 0.001,
  max.q = 0.001, max.e = 0.001, qval.method = "fdr", threshold = 1e-04,
  threshold.type = "pvalue", verbose = 0, RC = TRUE, use.freq = 1,
  shuffle.k = 2, shuffle.method = "euler", return.scan.results = FALSE,
  nthreads = 1, rng.seed = sample.int(10000, 1), motif_pvalue.k = 8,
  use.gaps = TRUE, allow.nonfinite = FALSE, warn.NA = TRUE,
  no.overlaps = TRUE, no.overlaps.by.strand = FALSE,
  no.overlaps.strat = "score", respect.strand = FALSE,
  motif_pvalue.method = c("dynamic", "exhaustive"),
  scan_sequences.qvals.method = c("BH", "fdr", "bonferroni"),
  mode = c("total.hits", "seq.hits"), pseudocount = 1)

Arguments

motifs	See convert_motifs() for acceptable motif formats.
sequences	XStringSet Sequences to scan. Alphabet should match motif.
bkg.sequences	XStringSet Optional. If missing, shuffle_sequences() is used to create background sequences from the input sequences.
max.p	numeric(1) P-value threshold.
max.q	numeric(1) Adjusted P-value threshold. This is only useful if multiple motifs are being enriched for.
max.e	numeric(1). The E-value is calculated by multiplying the P-value with the number of input motifs times two (McLeay and Bailey 2010).
qval.method	character(1) See stats::p.adjust().
threshold	numeric(1) See details.
threshold.type	character(1) One of c('pvalue', 'qvalue', 'logodds', 'logodds.abs'). See details.
verbose	numeric(1) 0 for no output, 4 for max verbosity.
RC	logical(1) If TRUE, check reverse complement of the input sequences. Only available for DNA/RNA.
use.freq	numeric(1) The default, 1, uses the motif matrix (from the motif['motif'] slot) to search for sequences. If a higher number is used, then the matching k-let matrix from the motif['multifreq'] slot is used. See add_multifreq().
shuffle.k	numeric(1) The k-let size to use when shuffling input sequences. Only used if no background sequences are input. See shuffle_sequences().
shuffle.method	character(1) One of c('euler', 'markov', 'linear'). See shuffle_sequences().
return.scan.results	logical(1) Return output from scan_sequences(). For large jobs, leaving this as FALSE can save a small amount time by preventing construction of the complete results data.frame from scan_sequences().
nthreads	numeric(1) Run scan_sequences() in parallel with nthreads threads. nthreads = 0 uses all available threads. Note that no speed up will occur for jobs with only a single motif and sequence.
rng.seed	numeric(1) Set random number generator seed. Since shuffling can occur simultaneously in multiple threads using C++, it cannot communicate with the regular R random number generator state and thus requires an independent seed. Each individual sequence in an XStringSet object will be given the following seed: rng.seed * index. See shuffle_sequences().
motif_pvalue.k	numeric(1) Control motif_pvalue() approximation. See motif_pvalue().
use.gaps	logical(1) Set this to FALSE to ignore motif gaps, if present.
allow.nonfinite	logical(1) If FALSE, then apply a pseudocount if non-finite values are found in the PWM. Note that if the motif has a pseudocount greater than zero and the motif is not currently of type PWM, then this parameter has no effect as the pseudocount will be applied automatically when the motif is converted to a PWM internally. This value is set to FALSE by default in order to stay consistent with pre-version 1.8.0 behaviour. A message will be printed if a pseudocount is applied. To disable this, set options(pseudocount.warning=FALSE).
warn.NA	logical(1) Whether to warn about the presence of non-standard letters in the input sequence, such as those in masked sequences.
no.overlaps	logical(1) Remove overlapping hits from the same motifs. Overlapping hits from different motifs are preserved. Please note that the current implementation of this feature can add significantly to the run time for large inputs.
no.overlaps.by.strand	logical(1) Whether to discard overlapping hits from the opposite strand (TRUE), or to only discard overlapping hits on the same strand (FALSE).
no.overlaps.strat	character(1) One of c("score", "order"). The former option keeps the highest scoring overlapping hit (and the first of these within ties), and the latter simply keeps the first overlapping hit.
respect.strand	logical(1) If motifs are DNA/RNA, then setting this option to TRUE will make scan_sequences() only scan the strands of the input sequences as indicated in the motif strand slot.
motif_pvalue.method	character(1) One of c("dynamic", "exhaustive"). Algorithm used for calculating P-values. The "exhaustive" method involves finding all possible motif matches at or above the specified score using a branch-and-bound algorithm, which can be computationally intensive (Hartman et al., 2013). Additionally, the computation must be repeated for each hit. The "dynamic" method calculates the distribution of possible motif scores using a much faster dynamic programming algorithm, and can be recycled for multiple scores (Grant et al., 2011). The only disadvantage is the inability to use allow.nonfinite = TRUE. See motif_pvalue() for details.
scan_sequences.qvals.method	character(1) One of c("fdr", "BH", "bonferroni"). The method for calculating adjusted P-values for individual motif hits. These are described in depth in the Sequence Searches vignette.
mode	character(1) One of c("total.hits", "seq.hits"). The former enriches for the total count of motif hits across all sequences, whereas the latter only counts motif hits once per sequence (useful for cases where there are many small sequences).
pseudocount	integer(1) Add a pseudocount to the motif hit counts when performing the Fisher test.

Details

To find enriched motifs, scan_sequences() is run on both target and background sequences. stats::fisher.test() is run to test for enrichment.

See scan_sequences() for more info on scanning parameters.

Value

DataFrame Enrichment results in a DataFrame. Function args and (optionally) scan results are stored in the metadata slot.

Author(s)

Benjamin Jean-Marie Tremblay [email protected]

References

McLeay R, Bailey TL (2010). “Motif Enrichment Analysis: A unified framework and method evaluation.” BMC Bioinformatics, 11.

See Also

scan_sequences(), shuffle_sequences(), add_multifreq(), motif_pvalue(), enrich_motifs2(), scan_sequences2(), match_bkg()

Examples

data(ArabidopsisPromoters)
data(ArabidopsisMotif)
if (R.Version()$arch != "i386") {
enrich_motifs(ArabidopsisMotif, ArabidopsisPromoters, threshold = 0.01)
}

---

Fast minimalist motif enrichment.

Description

enrich_motifs2() is a deliberately pared-down counterpart to enrich_motifs(), with a default surface that mirrors the command-line tool yamtk. It exposes a single p-value cutoff for the hits, a single q-value cutoff for the results, and two test modes ("seqs" and "sites", both Fisher's exact), and it leans on scan_sequences2() under the hood to scan the target and background sequences. The p-value adjustment is hard-coded to Benjamini-Hochberg.

Usage

enrich_motifs2(motifs, sequences, bkg.sequences = NULL, pvalue = 1e-04,
  qvalue = 0.1, test = c("seqs", "sites"), RC = TRUE, shuffle.k = 2L,
  rng.seed = sample.int(1e+06, 1), pseudocount = 1L, nthreads = 1)

Arguments

motifs	See convert_motifs() for accepted motif formats. DNA or RNA only; amino-acid and custom alphabet motifs are rejected.
sequences	XStringSet. Target sequences.
bkg.sequences	XStringSet or NULL. Background sequences. If NULL (default), target sequences are shuffled k-let-conserving via shuffle_sequences() with k = shuffle.k, matching both enrich_motifs() and ⁠yamtk enr⁠ default behaviour.
pvalue	numeric(1). P-value cutoff passed to scan_sequences2() for reporting individual hits. Default 1e-4.
qvalue	numeric(1). Result-level q-value cutoff. Motifs whose BH-adjusted enrichment p-value exceeds this are filtered out. Default 0.1.
test	character(1). One of "seqs" (default) or "sites". "seqs" runs Fisher's exact on per-sequence hit presence (does this sequence contain >=1 hit?), equivalent to enrich_motifs()'s mode = "seq.hits". "sites" runs Fisher's exact on per-position rates (how many hits per scannable position?), equivalent to mode = "total.hits".
RC	logical(1). If TRUE (default), scan both strands.
shuffle.k	integer(1). K-let size for the background shuffle (only used when bkg.sequences = NULL). Default 2L.
rng.seed	integer(1). RNG seed for background shuffling. Default sample.int(1e6, 1); set explicitly for reproducible runs.
pseudocount	integer(1). Pseudocount added uniformly to each cell of the Fisher 2x2 table; useful when one side has zero counts. Default 1L (matches yamtk's ⁠-p 1⁠ default).
nthreads	numeric(1). Number of threads passed to scan_sequences2() and shuffle_sequences(). nthreads = 0 uses all available threads.

Details

Use enrich_motifs() when you need any of: q-value adjustment methods other than BH; multifreq / gapped motifs; multiple-testing E-values; the respect.strand, allow.nonfinite, or non-pvalue threshold.type machinery; or amino-acid motifs.

Internally, scan_sequences2() is run once on the target set and once on the background set; the resulting hit tables are aggregated per motif into a 2x2 contingency table whose entries depend on test:

"seqs"

a = target sequences with >=1 hit; b = background sequences with >=1 hit; c = target sequences with no hits; d = background sequences with no hits.

"sites"

a = total hits in target; b = total hits in background; c = scannable target positions minus a; d = scannable background positions minus b. Scannable positions are computed per motif as sum(pmax(width(seqs) - motif_width + 1, 0)), doubled if RC = TRUE.

Each cell has pseudocount added before stats::fisher.test() is run with alternative = "greater". Q-values are computed via p.adjust(method = "BH") across all input motifs (not just the ones that pass the filter), then rows with qvalue > qvalue are dropped.

Value

A data.frame with one row per significant motif, sorted by pvalue ascending. Columns:

• motif: motif name

• motif.i: 1-based index into the input motifs

• consensus: IUPAC consensus of the motif

• target.seq.n: number of target sequences scanned

• target.seq.hits: number of target sequences containing >=1 hit

• target.site.hits: total hit count across all target sequences

• bkg.seq.n: number of background sequences scanned

• bkg.seq.hits: number of background sequences containing >=1 hit

• bkg.site.hits: total hit count across all background sequences

• enrichment: fold-enrichment (rate ratio; floor of 0.5/n in the denominator to avoid division by zero), per yamtk's formula

• log2.enrichment: log2(enrichment)

• pvalue: Fisher's exact one-sided ("greater") p-value

• qvalue: Benjamini-Hochberg adjusted p-value across motifs

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Tremblay BJM (2026). yamtk: Yet Another Motif ToolKit. https://github.com/bjmt/yamtk.

McLeay R, Bailey TL (2010). "Motif Enrichment Analysis: A unified framework and method evaluation." BMC Bioinformatics, 11.

See Also

enrich_motifs(), scan_sequences2(), shuffle_sequences(), match_bkg(), motif_pvalue()

Examples

library(universalmotif)
data(ArabidopsisPromoters)
data(ArabidopsisMotif)
if (R.Version()$arch != "i386") {
  enrich_motifs2(ArabidopsisMotif, ArabidopsisPromoters,
                 pvalue = 1e-3, qvalue = 0.5, rng.seed = 1)
}

---

Example motif in universalmotif format.

Description

A simple DNA motif. To recreate this motif: create_motif("TATAWAW", nsites = numeric())

Usage

examplemotif

Format

universalmotif

---

Another example motif in universalmotif format.

Description

A simple DNA motif with a non-empty multifreq slot. To recreate to this motif: add_multifreq(examplemotif, DNAStringSet(rep(c("CAAAACC", "CTTTTCC"), 3)))

Usage

examplemotif2

Format

universalmotif

---

Filter a list of motifs.

Description

Filter motifs based on the contents of available universalmotif slots. If the input motifs are not of universalmotif, then they will be converted for the duration of the filter_motifs() operation.

Usage

filter_motifs(motifs, name, altname, family, organism, width, alphabet, type,
  icscore, nsites, strand, pval, qval, eval, extrainfo)

Arguments

motifs	list See convert_motifs() for acceptable formats.
name	character Keep motifs by names.
altname	character Keep motifs by altnames.
family	character Keep motifs by family.
organism	character Keep motifs by organism.
width	numeric(1) Keep motifs with minimum width.
alphabet	character Keep motifs by alphabet.
type	character Keep motifs by type.
icscore	numeric(1) Keep motifs with minimum total IC.
nsites	numeric(1) Keep motifs with minimum number of target sites.
strand	character Keeps motifs by strand.
pval	numeric(1) Keep motifs by max P-value.
qval	numeric(1) Keep motifs by max Q-value.
eval	numeric(1) Keep motifs by max E-val.
extrainfo	character Named character vector of items that must be present in motif extrainfo slots.

Value

list Motifs. An attempt will be made to preserve the original class, see convert_motifs() for limitations.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

Examples

## By minimum IC:
jaspar <- read_jaspar(system.file("extdata", "jaspar.txt",
                                  package = "universalmotif"))
jaspar.ic3 <- filter_motifs(jaspar, icscore = 3)

## Starting from version 1.10.0 of the universalmotif package, one
## could instead make use of the universalmotif_df structure:
jaspar.ic3 <- jaspar |> to_df() |> subset(icscore > 3) |> to_list()

## By organism:
## Not run: 
library(MotifDb)
motifs <- convert_motifs(MotifDb)
motifs <- filter_motifs(motifs, organism = c("Athaliana", "Mmusculus"),
                        extrainfo = c("dataSource" = "cisbp_1.02"))
## Or:
motifs <- convert_motifs(MotifDb) |> to_df() |>
  subset(organism %in% c("Athaliana", "Mmusculus") &
    dataSource == "cisbp_1.02") |> to_list()

## End(Not run)

---

Polygon coordinates for plotting letters.

Description

DataFrame of polygon coordinates used by view_motifs() for plotting letters. It was generated using the createPolygons function from the gglogo package for the font Roboto Medium.

Usage

fontDFroboto

Format

DataFrame

---

Calculate sequence background.

Description

For a set of input sequences, calculate the overall sequence background for any k-let size. For very large sequences DNA and RNA sequences (in the billions of bases), please be aware of the much faster and more efficient Biostrings::oligonucleotideFrequency(). get_bkg() can still be used in these cases, though it may take several seconds or minutes to calculate the results (depending on requested k-let sizes).

Usage

get_bkg(sequences, k = 1:3, as.prob = NULL, pseudocount = 0,
  alphabet = NULL, to.meme = NULL, RC = FALSE, list.out = NULL,
  nthreads = 1, merge.res = TRUE, window = FALSE, window.size = 0.1,
  window.overlap = 0)

Arguments

sequences	XStringSet Input sequences. Note that if multiple sequences are present, the results will be combined into one (unless merge.res = FALSE).
k	integer Size of k-let. Background can be calculated for any k-let size.
as.prob	Deprecated.
pseudocount	integer(1) Add a count to each possible k-let. Prevents any k-let from having 0 or 1 probabilities.
alphabet	character(1) Provide a custom alphabet to calculate a background for. If NULL, then standard letters will be assumed for DNA, RNA and AA sequences, and all unique letters found will be used for BStringSet type sequences. Note that letters which are not a part of the standard DNA/RNA/AA alphabets or in the provided alphabet will not be counted in the totals during probability calculations.
to.meme	If not NULL, then get_bkg() will return the sequence background in MEME Markov Background Model format. Input for this argument will be used for cat(..., file = to.meme) within get_bkg(). See http://meme-suite.org/doc/bfile-format.html for a description of the format.
RC	logical(1) Calculate the background of the reverse complement of the input sequences as well. Only valid for DNA/RNA.
list.out	Deprecated.
nthreads	numeric(1) Run get_bkg() in parallel with nthreads threads. nthreads = 0 uses all available threads. Note that no speed up will occur for jobs with only a single sequence.
merge.res	logical(1) Whether to merge results from all sequences or return background data for individual sequences.
window	logical(1) Determine background in windows.
window.size	numeric Window size. If a number between 0 and 1 is provided, the value is calculated as the number multiplied by the sequence length.
window.overlap	numeric Overlap between windows. If a number between 0 and 1 is provided, the value is calculated as the number multiplied by the sequence length.

Value

If to.meme = NULL, a DataFrame with columns klet, count, and probability. If merge.res = FALSE, there will be an additional sequence column. If window = TRUE, there will be an additional start and stop columns.

If to.meme is not NULL, then NULL is returned, invisibly.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Bailey TL, Elkan C (1994). “Fitting a mixture model by expectation maximization to discover motifs in biopolymers.” Proceedings of the Second International Conference on Intelligent Systems for Molecular Biology, 2, 28-36.

See Also

create_sequences(), scan_sequences(), shuffle_sequences()

Examples

## Compare to Biostrings version
library(Biostrings)
seqs.DNA <- create_sequences()
bkg.DNA <- get_bkg(seqs.DNA, k = 3)
bkg.DNA2 <- oligonucleotideFrequency(seqs.DNA, 3, 1, as.prob = FALSE)
bkg.DNA2 <- colSums(bkg.DNA2)
all(bkg.DNA$count == bkg.DNA2)

## Create a MEME background file
get_bkg(seqs.DNA, k = 1:3, to.meme = stdout(), pseudocount = 1)

## Non-DNA/RNA/AA alphabets
seqs.QWERTY <- create_sequences("QWERTY")
bkg.QWERTY <- get_bkg(seqs.QWERTY, k = 1:2)

---

Implant sampled motif instances into sequences at known positions.

Description

implant_motifs() overwrites bases in an existing XStringSet with samples drawn column-by-column from one or more motif PPMs at chosen positions, returning the modified sequences together (optionally) with a ground-truth data.frame of where each implant was placed. This produces benchmark-grade positive sequences for testing scan_sequences2(), motif_finder(), enrich_motifs2(), motif_peaks(), or any other discovery / scanning pipeline against a known answer.

Usage

implant_motifs(motifs, sequences, n.per.seq = 1L, rate = NULL,
  positions = NULL, strand = c("both", "+"), centre.bias = 1L,
  min.spacing = 0L, max.retries = 100L, return.indices = FALSE)

Arguments

motifs	See convert_motifs() for accepted motif formats. DNA or RNA only, amino-acid and custom alphabet motifs are rejected. A single motif or a list of motifs.
sequences	XStringSet. Host sequences to plant into. Must share the motif alphabet.
n.per.seq	integer(1). Mode A. Number of instances to plant in each sequence. Default 1L.
rate	numeric(1) or NULL. Mode B. Per-bp Poisson rate. If non-NULL, overrides n.per.seq. Default NULL.
positions	list or NULL. Mode C. List of length length(sequences); element i is an integer vector of 1-based start positions. Overrides n.per.seq and rate. Default NULL.
strand	"both" or "+". Whether to also draw from the reverse complement. Default "both".
centre.bias	integer(1). Irwin-Hall N. 1L (default) is uniform; larger values cluster implants near the centre of each sequence. Ignored in positions mode.
min.spacing	integer(1). Minimum gap (in bp) between the end of one implant and the start of the next within the same sequence. 0L (default) allows abutting but not overlapping implants. Ignored in positions mode.
max.retries	integer(1). Per-implant cap on placement attempts before giving up (with a warning at the end if any targets fell short). Default 100L.
return.indices	logical(1). If FALSE (default), return only the modified XStringSet. If TRUE, return a list with the modified sequences plus a data.frame of per-implant ground-truth metadata.

Details

Three insertion modes (mutually exclusive; later non-NULL argument wins):

1. Fixed count (n.per.seq): plant exactly n.per.seq instances in every sequence. Matches HOMER-style ground-truth simulations.

2. Poisson rate (rate): per-sequence count drawn from rpois(1, rate * width(seq)). Density scales naturally with sequence length.

3. Explicit positions (positions): a list of length length(sequences); element i is a 1-based integer vector of start positions for that sequence. Bypasses centre.bias, min.spacing, and max.retries.

For each implant the function (a) picks a motif uniformly at random from motifs, (b) chooses a strand (uniform +/- if strand = "both"), (c) picks a start position (uniform or centre- biased via Irwin-Hall, retried up to max.retries times if it collides with an already-placed implant within min.spacing), (d) samples one instance from the motif's PPM column-by-column, reverse-complements if needed, and (e) overwrites that span in the target sequence.

Value

If return.indices = FALSE: an XStringSet of the same length, width, and class as sequences, with bases overwritten at the implant positions. Names are preserved. If return.indices = TRUE: a list with two elements – sequences (the modified XStringSet, as above) and indices (a data.frame with columns sequence.i (1-based), motif.i (1-based index into motifs), start (1-based), width, strand ("+" / "-"), and planted (the actual implanted string, already reverse-complemented when strand == "-")).

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

sample_sites(), create_sequences(), scan_sequences2(), motif_finder(), motif_peaks(), enrich_motifs2()

Examples

## Not run: 
library(universalmotif)
data(ArabidopsisMotif)
seqs  <- create_sequences("DNA", seqnum = 100, seqlen = 500)
set.seed(1)
res   <- implant_motifs(ArabidopsisMotif, seqs, n.per.seq = 1,
                        return.indices = TRUE)
out   <- res$sequences
truth <- res$indices
## Recover-rate check
hits <- scan_sequences(ArabidopsisMotif, out, threshold = 0.85,
                       threshold.type = "logodds")
## Density check
peaks <- motif_peaks(out, motif = ArabidopsisMotif)
plot_motif_peaks(peaks)

## End(Not run)

---

JASPAR2018 CORE database scores

Description

For use with compare_motifs(). The precomputed scores allow for fast P-value estimation. These scores were generated using make_DBscores() with the JASPAR2018 CORE motif set. The scores are organized in a DataFrame. In this DataFrame is the location and scale of scores resulting from a statistical distribution using the the comparisons of JASPAR2018 CORE motifs with randomized motifs of the specified subject and target motif length. Created using make_DBscores() from universalmotif v1.4.0. The parameters used can be seen via S4Vectors::metadata(JASPAR2018_CORE_DBSCORES).

Usage

JASPAR2018_CORE_DBSCORES

Format

DataFrame with function args in the metadata slot.

---

Create P-value databases.

Description

Generate data used by compare_motifs() for P-value calculations. By default, compare_motifs() uses an internal database based on the JASPAR2018 core motifs (Khan et al. 2018). Parameters for distributions are are estimated for every combination of motif widths.

Usage

make_DBscores(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"))

Arguments

db.motifs	list Database motifs.
method	character(1) One of PCC, EUCL, SW, KL, ALLR, BHAT, HELL, SEUCL, MAN, ALLR_LL, WEUCL, WPCC. See details.
shuffle.db	logical(1) Deprecated. Does nothing. generate random motifs with create_motif().
shuffle.k	numeric(1) See shuffle_motifs().
shuffle.method	character(1) See shuffle_motifs().
rand.tries	numeric(1) Approximate number of comparisons to perform for every combination of widths.
widths	numeric Motif widths to use in P-value database calculation.
min.position.ic	numeric(1) Minimum information content required between individual alignment positions for it to be counted in the final alignment score. It is recommended to use this together with normalise.scores = TRUE, as this will help punish scores resulting from only a fraction of an alignment.
normalise.scores	logical(1) Favour alignments which leave fewer unaligned positions, as well as alignments between motifs of similar length. Similarity scores are multiplied by the ratio of aligned positions to the total number of positions in the larger motif, and the inverse for distance scores.
min.overlap	numeric(1) Minimum overlap required when aligning the motifs. Setting this to a number higher then the width of the motifs will not allow any overhangs. Can also be a number between 0 and 1, representing the minimum fraction that the motifs must overlap.
min.mean.ic	numeric(1) Minimum mean information content between the two motifs for an alignment to be scored. This helps prevent scoring alignments between low information content regions of two motifs. Note that this can result in some comparisons failing if no alignment passes the mean IC threshold. Use average_ic() to filter out low IC motifs to get around this if you want to avoid getting NAs in your output.
progress	logical(1) Show progress.
nthreads	numeric(1) Run compare_motifs() in parallel with nthreads threads. nthreads = 0 uses all available threads.
tryRC	logical(1) Try the reverse complement of the motifs as well, report the best score.
score.strat	character(1) How to handle column scores calculated from motif alignments. "sum": add up all scores. "a.mean": take the arithmetic mean. "g.mean": take the geometric mean. "median": take the median. "wa.mean", "wg.mean": weighted arithmetic/geometric mean. "fzt": Fisher Z-transform. Weights are the total information content shared between aligned columns.

Details

See compare_motifs() for more info on comparison parameters.

To replicate the internal 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.

Value

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.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Khan A, Fornes O, Stigliani A, Gheorghe M, Castro-Mondragon JA, van der Lee R, Bessy A, Cheneby J, Kulkarni SR, Tan G, Baranasic D, Arenillas DJ, Sandelin A, Vandepoele K, Lenhard B, Ballester B, Wasserman WW, Parcy F, Mathelier A (2018). “JASPAR 2018: update of the open-access database of transcription factor binding profiles and its web framework.” Nucleic Acids Research, 46, D260-D266.

See Also

compare_motifs()

Examples

## Not run: 
library(MotifDb)
motifs <- convert_motifs(MotifDb[1:100])
scores <- make_DBscores(motifs, method = "PCC")
compare_motifs(motifs, 1:100, db.scores = scores)

## End(Not run)

---

Sample composition-matched background sequences from a universe.

Description

For each sequence in sequences, match_bkg() draws one or more background sequences from universe that share its GC fraction and length. This is the HOMER-style binned matching used by motif enrichment tools to control for sequence-composition biases (GC, length) that would otherwise inflate enrichment estimates relative to shuffle-only nulls. Matching can be extended to additional external per-sequence covariates (e.g. ChIP signal strength, conservation, mappability) via covariates, which simply add more binned axes to the grid.

Usage

match_bkg(sequences, universe = NULL, genome = NULL, exclude = NULL,
  n.candidates = NULL, chromosomes = NULL, n.per.target = 1L,
  n.bins.gc = 20L, n.bins.length = 10L, covariates = NULL,
  universe.covariates = NULL, n.bins.covariates = 10L,
  bin.type.covariates = "quantile", unique = TRUE,
  return.indices = FALSE)

Arguments

sequences	XStringSet. Target (positive) sequences. DNA or RNA only – GC matching is undefined for other alphabets.
universe	XStringSet or NULL. Pool of candidate background sequences. Must share sequences's alphabet and be at least n.per.target * length(sequences) long when unique = TRUE. Mutually exclusive with genome: supply exactly one of the two.
genome	A BSgenome object, or NULL. When supplied, the universe is built internally by sampling n.candidates random genomic windows whose widths match the target width distribution, optionally excluding any window that overlaps exclude. This is the typical shortcut for ChIP-seq style backgrounds, where the universe would otherwise be a manual sample() + subsetByOverlaps() + getSeq() boilerplate. Requires the BSgenome, GenomicRanges, and IRanges packages.
exclude	GRanges or NULL. Only honoured when genome is non-NULL. Random windows that overlap any range in exclude are dropped before sampling. Typical use: pass the target peak set so the background doesn't sample from peak regions.
n.candidates	integer(1) or NULL. Only honoured when genome is non-NULL. Number of random windows to sample before any exclude filtering. Default NULL picks max(10000, length(sequences) * 10).
chromosomes	character or NULL. Only honoured when genome is non-NULL. Names of seqlevels to sample from. Default NULL uses GenomeInfoDb::standardChromosomes(genome) when available (typically the autosomes plus the sex chromosomes), falling back to all seqlevels of the genome.
n.per.target	integer(1). Number of background sequences to draw per target. Default 1L.
n.bins.gc	integer(1). Number of equal-width GC bins. Default 20L.
n.bins.length	integer(1). Number of quantile-based length bins. Default 10L.
covariates	NULL, or a data.frame / matrix / numeric vector of external per-sequence covariates to additionally match on, with one row per sequence in sequences and one numeric column per covariate. A bare numeric vector is treated as a single covariate named "covariate1". Each covariate becomes an extra binned matching axis. Column names must not be "i", "gc" or "length" (these collide with the return.indices columns). Covariate values must be finite (no NA/NaN/Inf). Mutually exclusive with genome: internally-sampled genomic windows have no covariate values, so covariate matching requires an explicit universe paired with universe.covariates. Default NULL.
universe.covariates	NULL, or the covariate table for the universe sequences. Required (and only used) when covariates is supplied. Must have one row per universe sequence and the same column names, in the same order, as covariates. Default NULL.
n.bins.covariates	integer(1) applied to every covariate, or an integer vector named by covariate column (naming every covariate). Number of bins per covariate axis; each must be ⁠>= 2⁠. Default 10L.
bin.type.covariates	character. How to bin each covariate axis: "quantile" (default, breaks from the universe-covariate quantiles, like length) or "equalwidth" (equal-width breaks over the combined target/universe range, like GC). A single value applies to every covariate, or supply a character vector named by covariate column. Quantile binning is rank-based, hence scale- and shift-invariant; no covariate normalisation is needed.
unique	logical(1). If TRUE (default), each universe sequence is drawn at most once. Falls back to with-replacement sampling (with a warning) if the universe is too small to honour uniqueness.
return.indices	logical(1). If TRUE, return a data.frame of per-pair match indices and composition values instead of the matched XStringSet. Default FALSE.

Details

Algorithm (HOMER-style):

1. Compute per-sequence GC fraction and length for sequences and universe (plus any supplied covariates).

2. Bin the universe into a grid: GC bins equal-width over [0, 1], length bins quantile-based on the universe widths, and one extra axis per covariate (quantile- or equal-width binned, see bin.type.covariates).

3. For each target, locate its joint bin and sample n.per.target universe sequences from that bin (without replacement when unique = TRUE). If the bin is empty or undersized, expand outward in Manhattan-distance rings.

Use the result as a bkg.sequences argument to enrich_motifs2() or motif_finder() for a composition-controlled null. shuffle_sequences() is a complementary alternative that randomises each input in place (preserves k-let composition); match_bkg() instead samples real sequences from a larger pool.

Value

If return.indices = FALSE (default): an XStringSet (same subclass as sequences) of length length(sequences) * n.per.target. The first n.per.target entries correspond to sequences[[1]], the next to sequences[[2]], etc. If return.indices = TRUE: a data.frame with columns target.i, target.gc, target.length, universe.i, universe.gc, universe.length, followed (when covariates is supplied) by a ⁠target.<name>⁠ / ⁠universe.<name>⁠ pair for each covariate.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Heinz S, Benner C, Spann N, et al. (2010). "Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities." Molecular Cell, 38(4):576-589. doi:10.1016/j.molcel.2010.05.004.

See Also

shuffle_sequences(), enrich_motifs2(), motif_finder(), plot_match_bkg()

Examples

## Not run: 
library(universalmotif)
data(ArabidopsisPromoters)
target   <- ArabidopsisPromoters[1:10]
universe <- ArabidopsisPromoters[11:50]
bkg <- match_bkg(target, universe, n.per.target = 2)
## Use as a null background for enrichment / discovery
# enrich_motifs2(motifs, target, bkg.sequences = bkg)
# motif_finder(target, bkg.sequences = bkg)
plot_match_bkg(target, bkg)

## Additionally match on an external covariate (e.g. a ChIP signal
## score aligned to each sequence). Covariate matching needs an
## explicit universe + universe.covariates.
sig.t <- runif(length(target))
sig.u <- runif(length(universe))
idx <- match_bkg(target, universe,
                 covariates          = data.frame(signal = sig.t),
                 universe.covariates = data.frame(signal = sig.u),
                 return.indices      = TRUE)
plot_match_bkg(indices = idx)   # GC, length and signal panels

## End(Not run)

---

Merge motifs.

Description

Aligns the motifs using compare_motifs(), then averages the motif PPMs. Currently the multifreq slot, if filled in any of the motifs, will be dropped. Only 0-order background probabilities will be kept. Motifs are merged one at a time, starting with the first entry in the list. For merging regular DNA/RNA motifs, see the simpler and faster merge_motifs2().

Usage

merge_motifs(motifs, method = "ALLR", use.type = "PPM", min.overlap = 6,
  min.mean.ic = 0.25, tryRC = TRUE, relative_entropy = FALSE,
  normalise.scores = FALSE, min.position.ic = 0, score.strat = "sum",
  new.name = NULL)

Arguments

motifs	See convert_motifs() for acceptable motif formats.
method	character(1) One of PCC, EUCL, SW, KL, ALLR, BHAT, HELL, SEUCL, MAN, ALLR_LL, WEUCL, WPCC. See details.
use.type	character(1) One of 'PPM' and 'ICM'. The latter allows for taking into account the background frequencies if relative_entropy = TRUE. Note that 'ICM' is not allowed when method = c("ALLR", "ALLR_LL").
min.overlap	numeric(1) Minimum overlap required when aligning the motifs. Setting this to a number higher then the width of the motifs will not allow any overhangs. Can also be a number between 0 and 1, representing the minimum fraction that the motifs must overlap.
min.mean.ic	numeric(1) Minimum mean information content between the two motifs for an alignment to be scored. This helps prevent scoring alignments between low information content regions of two motifs. Note that this can result in some comparisons failing if no alignment passes the mean IC threshold. Use average_ic() to filter out low IC motifs to get around this if you want to avoid getting NAs in your output.
tryRC	logical(1) Try the reverse complement of the motifs as well, report the best score.
relative_entropy	logical(1) Change the ICM calculation affecting min.position.ic and min.mean.ic. See convert_type().
normalise.scores	logical(1) Favour alignments which leave fewer unaligned positions, as well as alignments between motifs of similar length. Similarity scores are multiplied by the ratio of aligned positions to the total number of positions in the larger motif, and the inverse for distance scores.
min.position.ic	numeric(1) Minimum information content required between individual alignment positions for it to be counted in the final alignment score. It is recommended to use this together with normalise.scores = TRUE, as this will help punish scores resulting from only a fraction of an alignment.
score.strat	character(1) How to handle column scores calculated from motif alignments. "sum": add up all scores. "a.mean": take the arithmetic mean. "g.mean": take the geometric mean. "median": take the median. "wa.mean", "wg.mean": weighted arithmetic/geometric mean. "fzt": Fisher Z-transform. Weights are the total information content shared between aligned columns.
new.name	character(1), NULL Instead of collapsing existing names (if NULL), assign a new one manually for the merged motif.

Details

See compare_motifs() for more info on comparison parameters.

If using a comparison metric where 0s are not allowed (KL, ALLR, ALLR_LL, IS), then pseudocounts will be added internally. These pseudocounts are only used for comparison and alignment, and are not used in the final merging step.

Note: score.strat = "a.mean" is NOT recommended, as merge_motifs() will not discriminate between two alignments with equal mean scores, even if one alignment is longer than the other.

Value

A single motif object. See convert_motifs() for available formats.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

compare_motifs(), merge_motifs2(), merge_similar(), trim_motifs()

Examples

## Not run: 
library(MotifDb)
merged.motif <- merge_motifs(MotifDb[1:5])

## End(Not run)

m1 <- create_motif("TTAAACCCC", name = "1")
m2 <- create_motif("AACC", name = "2")
m3 <- create_motif("AACCCCGG", name = "3")
view_motifs(merge_motifs(c(m1, m2, m3)))

---

Merge a list of motifs into a single consensus motif.

Description

merge_motifs2() is a faster minimalist counterpart to merge_motifs() that aligns motifs using compare_motifs2()'s alignment finder and averages their position-probability columns in a single shared coordinate frame.

Usage

merge_motifs2(motifs, min.overlap = 5L, RC = TRUE, new.name = NULL,
  weighted = FALSE, nthreads = 1)

Arguments

motifs	list of motifs. See convert_motifs() for accepted formats. DNA or RNA only.
min.overlap	integer(1). Minimum overlap (in columns) for an alignment to be considered. Default 5L. Forwarded to compare_motifs2().
RC	logical(1). If TRUE (default), also test reverse-complement alignments. Motifs whose best alignment hits the - strand of the anchor are reverse-complemented before averaging, so the merged motif is always returned in the anchor's orientation.
new.name	character(1) or NULL. Name for the merged motif. If NULL (default), set to paste0("merged_", paste(names, collapse = "+")).
weighted	logical(1). If TRUE, weight each input motif's column contribution by its ⁠@nsites⁠ slot. Default FALSE (simple unweighted mean). Motifs whose ⁠@nsites⁠ is missing or zero fall back to weight 1 in either mode.
nthreads	numeric(1). Number of threads passed to compare_motifs2(). nthreads = 0 uses all available threads.

Details

For each non-anchor motif, compare_motifs2() gives the best alignment offset, strand, and overlap against the anchor. Each aligned motif (reverse-complemented if its best alignment is - strand) is then placed in a unified column coordinate frame whose origin coincides with the anchor's first column. The output's column at frame position p is the (optionally ⁠@nsites⁠-weighted) mean of the PPM columns from all input motifs that cover p; positions covered by only some motifs are still emitted (no automatic IC-trimming). Run trim_motifs() on the result if you want to shrink low-information edges.

The merged motif's ⁠@bkg⁠ is the unweighted mean of the input motifs' ⁠@bkg⁠ vectors (or uniform if missing). ⁠@nsites⁠ is the sum of input ⁠@nsites⁠.

Value

A single universalmotif S4 object of type "PPM".

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

merge_motifs(), merge_similar2(), compare_motifs2(), trim_motifs()

Examples

library(universalmotif)
m1 <- create_motif("TTGACATA", name = "a")
m2 <- create_motif("CTTGACAT", name = "b")
m3 <- create_motif("TGACATAT", name = "c")
merged <- merge_motifs2(list(m1, m2, m3))
merged

---

Identify and merge similar motifs within a collection of motifs (or simply cluster motifs).

Description

Given a list of motifs, merge_similar() will identify similar motifs with compare_motifs(), and merge similar ones with merge_motifs(). For merging regular DNA/RNA motifs, see the simpler and faster merge_similar2().

Usage

merge_similar(motifs, threshold = 0.95, threshold.type = "score.abs",
  method = "PCC", use.type = "PPM", min.overlap = 6, min.mean.ic = 0,
  tryRC = TRUE, relative_entropy = FALSE, normalise.scores = FALSE,
  min.position.ic = 0, score.strat.compare = "a.mean",
  score.strat.merge = "sum", nthreads = 1, return.clusters = FALSE)

Arguments

motifs	See convert_motifs() for acceptable motif formats.
threshold	numeric(1) The minimum (for similarity metrics) or maximum (for distance metrics) threshold score for merging.
threshold.type	character(1) Type of score used for thresholding. Currently unused.
method	character(1) One of PCC, EUCL, SW, KL, BHAT, HELL, SEUCL, MAN, WEUCL, WPCC. See compare_motifs(). (The ALLR and ALLR_LL methods cannot be used for distance matrix construction.)
use.type	character(1) One of 'PPM' and 'ICM'. The latter allows for taking into account the background frequencies if relative_entropy = TRUE. Note that 'ICM' is not allowed when method = c("ALLR", "ALLR_LL").
min.overlap	numeric(1) Minimum overlap required when aligning the motifs. Setting this to a number higher then the width of the motifs will not allow any overhangs. Can also be a number between 0 and 1, representing the minimum fraction that the motifs must overlap.
min.mean.ic	numeric(1) Minimum mean information content between the two motifs for an alignment to be scored. This helps prevent scoring alignments between low information content regions of two motifs. Note that this can result in some comparisons failing if no alignment passes the mean IC threshold. Use average_ic() to filter out low IC motifs to get around this if you want to avoid getting NAs in your output.
tryRC	logical(1) Try the reverse complement of the motifs as well, report the best score.
relative_entropy	logical(1) Change the ICM calculation affecting min.position.ic and min.mean.ic. See convert_type().
normalise.scores	logical(1) Favour alignments which leave fewer unaligned positions, as well as alignments between motifs of similar length. Similarity scores are multiplied by the ratio of aligned positions to the total number of positions in the larger motif, and the inverse for distance scores.
min.position.ic	numeric(1) Minimum information content required between individual alignment positions for it to be counted in the final alignment score. It is recommended to use this together with normalise.scores = TRUE, as this will help punish scores resulting from only a fraction of an alignment.
score.strat.compare	character(1) The score.strat parameter used by compare_motifs(). For clustering purposes, the "sum" option cannot be used.
score.strat.merge	character(1) The score.strat parameter used by merge_motifs(). As discussed in merge_motifs(), the "sum" option is recommended over "a.mean" to maximize the overlap between motifs.
nthreads	numeric(1) Run compare_motifs() in parallel with nthreads threads. nthreads = 0 uses all available threads.
return.clusters	logical(1) Return the clusters instead of merging.

Details

See compare_motifs() for more info on comparison parameters, and merge_motifs() for more info on motif merging.

Value

See convert_motifs() for available output formats.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

compare_motifs(), merge_motifs(), merge_similar2()

Examples

## Not run: 
library(MotifDb)
motifs <- filter_motifs(MotifDb, family = "bHLH")[1:50]
length(motifs)
motifs <- merge_similar(motifs)
length(motifs)

## End(Not run)

---

Cluster similar motifs by significance and merge each cluster.

Description

merge_similar2() is a faster minimalist counterpart to merge_similar() that builds a pairwise similarity-significance graph using compare_motifs2()'s empirical or parametric null q-values, finds connected components on that graph, and collapses each component (cluster) into a single merged motif via merge_motifs2().

Usage

merge_similar2(motifs, qvalue = 0.01, null = c("empirical", "parametric"),
  min.overlap = 5L, RC = TRUE, bkg = NULL, weighted = FALSE,
  return.clusters = FALSE, nthreads = 1)

Arguments

motifs	list of motifs. DNA or RNA only.
qvalue	numeric(1). BH-adjusted q-value cutoff for linking two motifs in the significance graph. Default 0.01.
null	character(1). "empirical" (default) or "parametric", forwarded to compare_motifs2().
min.overlap	integer(1). Forwarded to compare_motifs2() / merge_motifs2(). Default 5L.
RC	logical(1). If TRUE (default), test reverse-complement alignments.
bkg	numeric(4) or NULL. Background frequencies, forwarded to compare_motifs2(). NULL uses uniform.
weighted	logical(1). Forwarded to merge_motifs2(). Default FALSE.
return.clusters	logical(1). If TRUE, return a universalmotif_df of cluster assignments instead of merged motifs. Default FALSE.
nthreads	numeric(1). Threads forwarded to compare_motifs2() and merge_motifs2(). nthreads = 0 uses all available threads.

Details

Unlike merge_similar() (which builds a distance matrix from raw similarity scores and applies hierarchical clustering at an absolute score cutoff) this function clusters on statistical significance: two motifs are linked if their pairwise q-value meets the user's qvalue cutoff. This matches the STAMP / TOMTOM-clustering semantics ("group all motifs that are significantly similar"), is deterministic, has no linkage-method choice to expose, and lets the user reason about the cutoff in q-value units rather than abstract distance.

The significance graph is built from a symmetrised q-value matrix: Qsym[i,j] = pmin(Q[i,j], Q[j,i]). compare_motifs2()'s per-query BH adjustment makes its native q-value matrix non-symmetric; the min- symmetrisation defines a pair as linked if either direction meets the cutoff (more permissive of the two, which is appropriate when the goal is to group similar motifs).

Clustering is by connected components on the resulting graph, implemented with a small union-find. Each component (size >= 2) becomes a single merged motif via merge_motifs2(); singletons (size == 1) pass through unchanged.

Value

If return.clusters = FALSE (default): a list of universalmotif S4 objects with similar motifs collapsed (the list is shorter than the input whenever any merging occurred). Singletons pass through unchanged.

If return.clusters = TRUE: a universalmotif_df (as produced by to_df()) describing the input motifs – one row per input motif, with the motif column holding the motif objects and name their names – augmented with two extra columns: motif.i (1-based input index) and cluster (1-based cluster id; motifs in the same cluster share an id).

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

merge_similar(), merge_motifs2(), compare_motifs2()

Examples

## Not run: 
library(universalmotif)
## A clear cluster of 3 related motifs + 1 unrelated motif
m1 <- create_motif("TTGACATA", name = "a")
m2 <- create_motif("CTTGACAT", name = "b")
m3 <- create_motif("TGACATAT", name = "c")
m4 <- create_motif("GGGCCCCC", name = "unrelated")
merged <- merge_similar2(list(m1, m2, m3, m4), qvalue = 0.05)
length(merged)  # 2: one merged + the unrelated singleton

## End(Not run)

---

Find significantly co-occurring motif pairs in a set of sequences.

Description

motif_coocc() tests every motif pair for over-co-occurrence across a set of sequences using a one-sided Fisher's exact test on the 2x2 contingency table of per-sequence presence/absence. The result is the long-format pair table, BH-corrected across all tested pairs. Optionally, when max.distance is non-NULL, the function also reports two descriptive spatial columns, both.clustered (number of co-occurring sequences with a within-max.distance (A, B) hit pair) and median.distance (the median nearest-pair spacing), so users can flag heterodimer-like arrangements. The Fisher p-value itself is always computed on the unfiltered 2x2 (the spatial filter is informational, not a test).

Usage

motif_coocc(motifs, sequences = NULL, hits = NULL, n.sequences = NULL,
  pvalue = 1e-04, RC = TRUE, max.distance = NULL, min.coocc = 1L,
  pseudocount = 0L, self.pairs = FALSE, nthreads = 1L)

Arguments

motifs	See convert_motifs() for accepted motif formats. A single motif or a list of motifs. Required in both paths to provide motif names and to dimension the pair table.
sequences	XStringSet or NULL. If supplied, the function scans internally via scan_sequences2(); DNA/RNA only.
hits	data.frame, GRanges, or NULL. Precomputed hit table with motif.i and sequence.i columns. If supplied, sequences is ignored and n.sequences must be provided. Any alphabet accepted on this path.
n.sequences	integer(1) or NULL. Total number of host sequences (needed to fill the "neither" cell of the 2x2 when using the precomputed-hits path).
pvalue	numeric(1). Per-hit p-value threshold for the internal scan. Default 1e-4. Ignored on the hit-table path.
RC	logical(1). Scan the reverse complement too? Default TRUE. Ignored when a hits table is provided.
max.distance	integer(1) or NULL. If non-NULL, two extra descriptive columns are added to the output: both.clustered (subset of co-occurring sequences with a within- max.distance (A, B) hit pair) and median.distance (median of those nearest-pair distances). The Fisher p-value is unchanged; it always tests "do A and B occur in the same sequences more often than chance?" on the unfiltered 2x2. Requires a start column on the hit table.
min.coocc	integer(1). Pairs with fewer than min.coocc co-occurring sequences are tagged with NA p- and q-values (not tested). Default 1L (only test pairs that co-occur at least once).
pseudocount	numeric(1). Added to every cell of the 2x2 contingency before the Fisher test (matches the convention from enrich_motifs2()). Default 0L.
self.pairs	logical(1). If TRUE, include ⁠(i, i)⁠ rows (testing whether motif i co-occurs with itself). Usually only meaningful for motifs that bind as homodimers. Default FALSE.
nthreads	integer(1). Threads for the internal scan (when sequences is supplied). The pairwise Fisher loop is single- threaded R code regardless.

Details

Two input paths are supported:

1. Internal scan (default): pass motifs + sequences and motif_coocc() will call scan_sequences2() with the given pvalue and RC arguments to build the hit table. DNA/RNA only (the restriction comes from scan_sequences2()).

2. Precomputed hits: pass motifs + hits + n.sequences. hits is a data.frame or GRanges with motif.i and sequence.i columns (the native return shape of scan_sequences() / scan_sequences2()), optionally start when max.distance is set. This path accepts any motif alphabet; once a hit table exists, co-occurrence is pure set arithmetic on integer indices and no sequence bytes are read.

Value

data.frame with columns:

• motif_a, motif_b: motif names (1-based indexing into the supplied motifs list; deduplicated names if duplicates exist).

• a_only, b_only, both, neither: 2x2 contingency cells.

• odds_ratio: conditional MLE odds ratio from fisher.test.

• pvalue: one-sided Fisher's exact p-value (alternative = "greater"). NA when both < min.coocc.

• qvalue: BH-corrected q-value across all tested pairs.

• both.clustered (spatial mode only): subset of both where a within-max.distance (A, B) hit pair exists.

• median.distance (spatial mode only): median of those nearest-pair distances; NA if no clustered pair. Rows are sorted by q-value ascending; NA-pvalue rows go last.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Heinz S, Benner C, Spann N, et al. (2010). "Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities." Molecular Cell, 38(4):576-589. doi:10.1016/j.molcel.2010.05.004. (HOMER's pairwise-Fisher co-occurrence convention.)

Whitington T, Frith MC, Johnson J, Bailey TL (2011). "Inferring transcription factor complexes from ChIP-seq data." Nucleic Acids Research, 39(15):e98. doi:10.1093/nar/gkr341. (SpaMo, the spatial-distance motif-pair analysis.)

See Also

scan_sequences2(), enrich_motifs2(), motif_peaks()

Examples

## Not run: 
library(universalmotif)
data(ArabidopsisMotif)
seqs <- create_sequences("DNA", seqnum = 200, seqlen = 500)
## Toy: implant two copies of the same motif jointly in all seqs.
set.seed(1)
imp <- implant_motifs(list(ArabidopsisMotif, ArabidopsisMotif),
                      seqs, n.per.seq = 1, return.indices = TRUE)
co <- motif_coocc(list(ArabidopsisMotif, ArabidopsisMotif),
                  imp$sequences, pvalue = 1e-3)
co

## End(Not run)

---

Discover de novo motifs in a set of sequences.

Description

motif_finder() is a minimalist de novo motif discovery function, whose defaults mirror the command-line tool yamtk. The pipeline works through a user-controlled range of motif widths; at each width it enumerates the over-represented k-mer seeds (using a Fisher's exact test on per-sequence presence against a background set), aligns the Hamming-1 neighbours of the best seed to form a PPM, refines that PPM over two re-scan passes, accepts the motif if its Fisher's exact p-value passes stop.pvalue, and then masks the covered positions and repeats, up to nmotifs times per width. Once every width is done, the motifs are deduplicated across widths by coverage overlap, BH-adjusted, IC-trimmed at the flanks, and returned in p-value order.

Usage

motif_finder(sequences, bkg.sequences = NULL, min.width = 6L,
  max.width = 15L, nmotifs = 10L, hit.pvalue = 1e-04,
  stop.pvalue = 0.001, qvalue = 0.001, dedup.overlap = 0.5, RC = TRUE,
  shuffle.k = 2L, rng.seed = sample.int(1e+06, 1), pseudocount = 1L,
  nthreads = 1)

Arguments

sequences	XStringSet. DNA or RNA target sequences.
bkg.sequences	XStringSet or NULL. Background sequences. If NULL (default), target sequences are shuffled k-let-conserving via shuffle_sequences() with k = shuffle.k, matching ⁠yamtk me⁠ default behaviour.
min.width	integer(1). Minimum motif width. Default 6L.
max.width	integer(1). Maximum motif width. Default 15L. Hard ceiling 30.
nmotifs	integer(1). Maximum number of motifs to discover per width before moving on. Default 10L.
hit.pvalue	numeric(1). P-value threshold for scoring a single hit during scanning. Default 1e-4.
stop.pvalue	numeric(1). Per-motif Fisher's exact p-value below which a candidate motif is accepted. Default 1e-3.
qvalue	numeric(1). Final BH-adjusted q-value cutoff. Motifs with qvalue above this are dropped. Default 1e-3.
dedup.overlap	numeric(1). Cross-width dedup overlap fraction. When the intersection of two motifs' positional coverage divided by the smaller motif's coverage exceeds this, the higher-p-value motif is dropped. Default 0.5.
RC	logical(1). If TRUE (default), scan both strands.
shuffle.k	integer(1). K-let size for the background shuffle (only used when bkg.sequences = NULL). Default 2L.
rng.seed	integer(1). RNG seed for background shuffling. Default sample.int(1e6, 1); set explicitly for reproducible runs.
pseudocount	integer(1). Pseudocount used when converting per-position counts to log-odds scores. Default 1L.
nthreads	numeric(1). Number of threads. Parallelism is at the width level. nthreads = 0 uses all available threads.

Details

Algorithm and defaults are a faithful port of yamtk, whose own design is in turn based on the STREME algorithm (Bailey 2021): seed enumeration via word counting with per-sequence Fisher's exact ranking, iterative PPM refinement on positive sequences, per-motif Fisher's exact significance against a shuffled or user-supplied background, and position masking between discoveries.

Motifs with qvalue > qvalue are dropped from the result. To see all discovered motifs regardless of significance, set qvalue = 1.

Value

A universalmotif_df (see to_df()) with one row per discovered motif, sorted by pvalue ascending. The standard to_df() columns (name, altname, family, organism, consensus, alphabet, type, nsites, pval, eval, motif, ...) are populated from each discovered motif. The yamtk-specific stats are carried as additional columns:

• rank: 1-based discovery rank

• width: motif width (post IC-trim)

• seqs_pos: positive sequences with >=1 hit

• seqs_neg: background sequences with >=1 hit

• sites_pos: total target hits

• sites_neg: total background hits

• n_pos: number of positive sequences

• n_neg: number of background sequences

• pvalue: Fisher's exact one-sided p-value

• qvalue: Benjamini-Hochberg q-value

Use to_list() on the returned object to recover a plain list of universalmotif S4 objects.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Bailey TL (2021). "STREME: accurate and versatile sequence motif discovery." Bioinformatics, 37(18), 2834-2840. doi:10.1093/bioinformatics/btab203.

Tremblay BJM (2026). yamtk: Yet Another Motif ToolKit. https://github.com/bjmt/yamtk.

See Also

scan_sequences2(), enrich_motifs2(), shuffle_sequences(), create_motif(), to_df(), to_list()

Examples

## Not run: 
library(universalmotif)
library(Biostrings)
set.seed(1)
planted <- DNAString("TTGACATA")
seqs <- create_sequences(seqnum = 100, seqlen = 200, rng.seed = 1)
## Plant the motif at random positions in 80% of sequences:
for (i in seq_len(80)) {
  pos <- sample.int(200 - 8, 1)
  seqs[[i]][pos:(pos + 7)] <- planted
}
motifs <- motif_finder(seqs, qvalue = 1, rng.seed = 1)
to_list(motifs)[[1]]

## End(Not run)

---

CentriMo-style positional enrichment of motif hits.

Description

motif_peaks() tests whether motif hits cluster non-uniformly along input sequences. It implements an analytical CentriMo-style test (Bailey & Machanick 2012): under the null that hit positions are uniformly distributed over ⁠[1, seq.length]⁠, the count of hits falling inside a candidate window follows a binomial distribution. For each motif and each candidate window, a one-sided binomial p-value is computed; the best-scoring window per motif is reported, Bonferroni-corrected over the number of windows tested, then BH- adjusted across motifs.

Usage

motif_peaks(hits, seq.length = NULL, mode = c("central", "local"),
  qvalue = 0.1, min.window = 10L, max.window = NULL, window.step = 10L,
  position.step = 5L, nthreads = 1)

Arguments

hits	Motif hit table from scan_sequences() or scan_sequences2(). Accepted as either a data.frame or a GRanges. Required columns / metadata: motif, score, start, and either end (scan_sequences2) or stop (scan_sequences). For data.frame input, a sequence / sequence.i column identifies which sequence each hit belongs to. For GRanges input, the sequence is taken from seqnames().
seq.length	integer(1). Common length of the input sequences (in bases). Required for data.frame input. For GRanges input, taken from seqlengths(hits) if not provided.
mode	character(1). "central" (default) tests only centre-of-sequence windows; "local" also varies the window centre across the sequence.
qvalue	numeric(1). BH-adjusted q-value cutoff for the final reported motifs. Default 0.1. Set to 1 to return every motif.
min.window, max.window	integer(1). Bounds of the window widths to test. max.window = NULL (default) uses floor(seq.length / 2). min.window defaults to 10L.
window.step	integer(1). Step between consecutive window widths. Default 10L.
position.step	integer(1). Step between consecutive window centres in mode = "local". Ignored in mode = "central". Default 5L.
nthreads	numeric(1). Number of threads. nthreads = 0 uses all available cores. Default 1.

Details

Two modes are supported. mode = "central" (default) tests only windows centred on seq.length / 2, which is the appropriate hypothesis for centred input sequences (e.g. ChIP-seq peak summits). mode = "local" additionally varies the window centre, scanning the whole sequence for any positionally-enriched region; this is more permissive and uses a larger Bonferroni correction.

All sequences are assumed to be of equal length; the function does not enforce centring on a biological reference, but its results are only meaningful if such centring is in place upstream.

Value

A data.frame with one row per motif passing the q-value cutoff, sorted by pvalue ascending. Columns: motif, motif.i, mode, seq.length, nhits, best.window, best.center, hits.in, hits.out, expected.in, enrichment, log2.enrichment, pvalue, qvalue, centers (list-column of integer vectors of hit centres, one entry per sequence after best-hit-per-sequence dedup; consumed by plot_motif_peaks()).

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Bailey TL, Machanick P (2012). "Inferring direct DNA binding from ChIP-seq." Nucleic Acids Research, 40(17):e128. doi:10.1093/nar/gks433.

See Also

scan_sequences2(), plot_motif_peaks()

Examples

## Not run: 
library(universalmotif)
library(Biostrings)
set.seed(1)
seqs <- create_sequences(seqnum = 200, seqlen = 500, rng.seed = 1)
planted <- DNAString("TTGACATA")
for (i in seq_len(150)) {
  pos <- 246L + sample.int(8, 1) - 1L
  subseq(seqs[[i]], start = pos, width = 8) <- planted
}
m <- create_motif("TTGACATA", name = "test")
hits <- scan_sequences2(m, seqs, pvalue = 1e-3, return.granges = TRUE)
motif_peaks(hits)

## End(Not run)

---

Motif P-value and scoring utility

Description

For calculating P-values and logodds scores from P-values for any number of motifs.

Usage

motif_pvalue(motifs, score, pvalue, bkg.probs, use.freq = 1, k = 8,
  nthreads = 1, rand.tries = 10, rng.seed = sample.int(10000, 1),
  allow.nonfinite = FALSE, method = c("dynamic", "exhaustive"))

Arguments

motifs	See convert_motifs() for acceptable motif formats.
score	numeric, list Get a P-value for a motif from a logodds score. See details for an explanation of how to vectorize the calculation for method = "dynamic".
pvalue	numeric, list Get a logodds score for a motif from a P-value. See details for an explanation of how to vectorize the calculation for method = "dynamic".
bkg.probs	numeric, list A vector background probabilities. If supplying individual background probabilities for each motif, a list of such vectors. If missing, retrieves the background from the motif bkg slot. Note that this option is only used when method = "dynamic", or when method = "exhaustive" and providing a P-value and returning a score; for the inverse, the motifs are first converted to PWMs via convert_type(), which uses the motif bkg slot for background adjustment.
use.freq	numeric(1) By default uses the regular motif matrix; otherwise uses the corresponding multifreq matrix. Max is 3 when method = "exhaustive".
k	numeric(1) For speed, scores/P-values can be approximated after subsetting the motif every k columns when method = "exhaustive". If k is a value equal or higher to the size of input motif(s), then the calculations are exact. The default, 8, is recommended to those looking for a good tradeoff between speed and accuracy for jobs requiring repeated calculations. Note that this is ignored when method = "dynamic", as subsetting is not required.
nthreads	numeric(1) Run motif_pvalue() in parallel with nthreads threads. nthreads = 0 uses all available threads. Currently only applied when method = "exhaustive".
rand.tries	numeric(1) When ncol(motif) < k and method = "exhaustive", an approximation is used. This involves randomly approximating the overall motif score distribution. To increase accuracy, the distribution is approximated rand.tries times and the final scores averaged. Note that this is ignored when method = "dynamic", as subsetting is not required.
rng.seed	numeric(1) In order to allow motif_pvalue() to perform C++ level parallelisation, it must work independently from R. This means it cannot communicate with R to get/set the R RNG state. To get around this, the RNG seed used by the C++ function can be set with rng.seed. To make sure each thread gets a different seed however, the seed is multiplied with the iteration count. For example: when working with two motifs, the second motif gets the following seed: rng.seed * 2. The default is to pick a random number as chosen by sample(), which effectively makes motif_pvalue() dependent on the R RNG state. Note that this is ignored when method = "dynamic", as the random subsetting is only used for method = "exhaustive".
allow.nonfinite	logical(1) If FALSE, then apply a pseudocount if non-finite values are found in the PWM. Note that if the motif has a pseudocount greater than zero and the motif is not currently of type PWM, then this parameter has no effect as the pseudocount will be applied automatically when the motif is converted to a PWM internally. Note that this option is incompatible with method = "dynamic". A message will be printed if a pseudocount is applied. To disable this, set options(pseudocount.warning=FALSE).
method	character(1) One of c("dynamic", "exhaustive"). Algorithm used for calculating P-values. The "exhaustive" method involves finding all possible motif matches at or above the specified score using a branch-and-bound algorithm, which can be computationally intensive (Hartman et al., 2013). Additionally, the computation must be repeated for each hit. The "dynamic" method calculates the distribution of possible motif scores using a much faster dynamic programming algorithm, and can be recycled for multiple scores (Grant et al., 2011). The only disadvantage is the inability to use allow.nonfinite = TRUE.

Details

Regarding vectorization

A note regarding vectorizing the calculation when method = "dynamic" (no vectorization is possible with method = "exhaustive"): to avoid performing the P-value/score calculation repeatedly for individual motifs, provide the score/pvalue arguments as a list, with each entry corresponding to the scores/P-values to be calculated for the respective motifs provided to motifs. If you simply provide a list of repeating motifs and a single numeric vector of corresponding input scores/P-values, then motif_pvalue() will not vectorize. See the Examples section.

The dynamic method

One of the algorithms available to motif_pvalue() to calculate scores or P-values is the dynamic programming algorithm used by FIMO (Grant et al., 2011). In this method, a small range of possible scores from the possible miminum and maximum is created and the cumulative probability of each score in this distribution is incrementally calculated using the logodds scores and the background probabilities. This distribution of scores and associated P-values can be used to calculate P-values or scores for any input, any number of times. This method scales well with large motifs, and multifreq representations. The only downside is that it is incompatible with allow.nonfinite = TRUE, as this would not allow for the creation of the initial range of scores. Although described for a different purpose, the basic premise of the dynamic programming algorithm is also described in Gupta et al. (2007).

The exhaustive method

Calculating P-values exhaustively for motifs can be very computationally intensive. This is due to how P-values must be calculated: for a given score, all possible sequences which score equal or higher must be found, and the probability for each of these sequences (based on background probabilities) summed. For a DNA motif of length 10, the number of possible unique sequences is 4^10 = 1,048,576. Finding all possible sequences higher than a given score can be done very efficiently and quickly with a branch-and-bound algorithm, but as the motif length increases even this calculation becomes impractical. To get around this, the P-value calculation can be approximated.

In order to calculate P-values for longer motifs, this function uses the approximation proposed by Hartmann et al. (2013), where the motif is subset, P-values calculated for the subsets, and finally combined for a total P-value. The smaller the size of the subsets, the faster the calculation; but also, the bigger the approximation. This can be controlled by setting k. In fact, for smaller motifs (< 13 positions) calculating exact P-values can be done individually in reasonable time by setting k = 12.

To calculate a score from a P-value, all possible scores are calculated and the (1 - pvalue) * 100 nth percentile score returned. When k < ncol(motif), the complete set of scores is instead approximated by randomly adding up all possible scores from each subset. Note that this approximation can actually be potentially quite expensive at times and even slower than the exact version; for jobs requiring lots of repeat calculations, a bit of benchmarking beforehand can be useful to find the optimal settings.

Please note that bugs are more likely to occur when using the exhaustive method, as the algorithm contains several times more code compared to the dynamic method. Unless you have a strong need to use allow.nonfinite = TRUE then avoid using this method.

Value

numeric, list A vector or list of vectors of scores/P-values.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Grant CE, Bailey TL, Noble WS (2011). "FIMO: scanning for occurrences of a given motif." Bioinformatics, 27, 1017-1018.

Gupta S, Stamatoyannopoulos JA, Bailey TL, Noble WS (2007). "Quantifying similarity between motifs." Genome Biology, 8, R24.

Hartmann H, Guthohrlein EW, Siebert M, Soding SLJ (2013). “P-value-based regulatory motif discovery using positional weight matrices.” Genome Research, 23, 181-194.

See Also

get_matches(), get_scores(), motif_range(), motif_score(), prob_match(), prob_match_bkg(), score_match()

Examples

if (R.Version()$arch != "i386") {

## P-value/score calculations are performed using the PWM version of the
## motif
data(examplemotif)

## Get a minimum score based on a P-value
motif_pvalue(examplemotif, pvalue = 0.001)

## Get the probability of a particular sequence hit
motif_pvalue(examplemotif, score = 0)

## The calculations can be performed for multiple motifs
motif_pvalue(c(examplemotif, examplemotif), pvalue = c(0.001, 0.0001))

## Compare score thresholds and P-value:
scores <- motif_score(examplemotif, c(0.6, 0.7, 0.8, 0.9))
motif_pvalue(examplemotif, scores)

## Calculate the probability of getting a certain match or better:
TATATAT <- score_match(examplemotif, "TATATAT")
TATATAG <- score_match(examplemotif, "TATATAG")
motif_pvalue(examplemotif, TATATAT)
motif_pvalue(examplemotif, TATATAG)

## Get all possible matches by P-value:
get_matches(examplemotif, motif_pvalue(examplemotif, pvalue = 0.0001))

## Vectorize the calculation for multiple motifs and scores/P-values:
m <- create_motif()
motif_pvalue(c(examplemotif, m), list(1:5, 2:3))
## The non-vectorized equivalent:
motif_pvalue(
  c(rep(list(examplemotif), 5), rep(list(m), 2)), c(1:5, 2:3)
)
}

---

Get the reverse complement of a DNA or RNA motif.

Description

For any motif, change the motif slot to it's reverse complement. If the multifreq slot is filled, then it is also applied. No other slots are affected.

Usage

motif_rc(motifs, ignore.alphabet = FALSE)

Arguments

motifs	See convert_motifs() for acceptable formats
ignore.alphabet	logical(1) If TRUE, then motif_rc() throws an error when it detects a non-DNA/RNA motif. If FALSE, it will proceed regardless.

Value

See convert_motifs() for available output formats.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

Examples

jaspar <- read_jaspar(system.file("extdata", "jaspar.txt",
                                  package = "universalmotif"))
jaspar.rc <- motif_rc(jaspar)

---

Generate ggplot2 motif trees with ggtree.

Description

For more powerful motif tree functions, see the motifStack package. The motif_tree() function compares motifs with compare_motifs() to create a distance matrix, which is used to generate a phylogeny. This can be plotted with ggtree::ggtree(). The purpose of this function is simply to combine the compare_motifs() and ggtree::ggtree() steps into one. For more control over tree creation, it is recommend to do these steps separately. See the "Motif comparisons and P-values" vignette for such a workthrough. This function requires the ape and ggtree packages to be installed separately. For a leaner alternative using compare_motifs2(), see motif_tree2().

Usage

motif_tree(motifs, layout = "circular", linecol = "family",
  labels = "none", tipsize = "none", legend = TRUE,
  branch.length = "none", db.scores, method = "EUCL", use.type = "PPM",
  min.overlap = 6, min.position.ic = 0, tryRC = TRUE, min.mean.ic = 0,
  relative_entropy = FALSE, progress = FALSE, nthreads = 1,
  score.strat = "a.mean", ...)

Arguments

motifs	list, dist See convert_motifs() for available formats. Alternatively, the resulting comparison matrix from compare_motifs() (run as.dist(results) beforehand; if the comparison was performed with a similarity metric, make sure to convert to distances first).
layout	character(1) One of c('rectangular', 'slanted', 'fan', 'circular', 'radial', 'equal_angle', 'daylight'). See ggtree::ggtree().
linecol	character(1) universalmotif slot to use to colour lines (e.g. 'family'). Not available for dist input (see examples for how to add it manually). See ggtree::ggtree().
labels	character(1) universalmotif slot to use to label tips (e.g. 'name'). For dist input, only 'name' is available. See ggtree::ggtree().
tipsize	character(1) universalmotif slot to use to control tip size (e.g. 'icscore'). Not available for dist input (see examples for how to add it manually). See ggtree::ggtree().
legend	logical(1) Show legend for line colour and tip size. See ggtree::ggtree().
branch.length	character(1) If 'none', draw a cladogram. See ggtree::ggtree().
db.scores	data.frame See compare_motifs().
method	character(1) One of PCC, EUCL, SW, KL, ALLR, BHAT, HELL, SEUCL, MAN, ALLR_LL, WEUCL, WPCC. See details.
use.type	character(1)c('PPM', 'ICM')⁠. The latter allows for taking into account the background frequencies (only if ⁠relative_entropy = TRUE'). See compare_motifs().
min.overlap	numeric(1) Minimum overlap required when aligning the motifs. Setting this to a number higher then the width of the motifs will not allow any overhangs. Can also be a number between 0 and 1, representing the minimum fraction that the motifs must overlap.
min.position.ic	numeric(1) Minimum information content required between individual alignment positions for it to be counted in the final alignment score. It is recommended to use this together with normalise.scores = TRUE, as this will help punish scores resulting from only a fraction of an alignment.
tryRC	logical(1) Try the reverse complement of the motifs as well, report the best score.
min.mean.ic	numeric(1) Minimum mean information content between the two motifs for an alignment to be scored. This helps prevent scoring alignments between low information content regions of two motifs. Note that this can result in some comparisons failing if no alignment passes the mean IC threshold. Use average_ic() to filter out low IC motifs to get around this if you want to avoid getting NAs in your output.
relative_entropy	logical(1) Change the ICM calculation affecting min.position.ic and min.mean.ic. See convert_type().
progress	logical(1) Show message regarding current step.
nthreads	numeric(1) Run compare_motifs() in parallel with nthreads threads. nthreads = 0 uses all available threads.
score.strat	character(1) How to handle column scores calculated from motif alignments. "sum": add up all scores. "a.mean": take the arithmetic mean. "g.mean": take the geometric mean. "median": take the median. "wa.mean", "wg.mean": weighted arithmetic/geometric mean. "fzt": Fisher Z-transform. Weights are the total information content shared between aligned columns.
...	ggtree params. See ggtree::ggtree().

Details

See compare_motifs() for more info on comparison parameters.

Value

ggplot object.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

References

Wickham H (2009). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-0-387-98140-6, <URL: http://ggplot2.org>.

Yu G, Smith D, Zhu H, Guan Y, Lam TT (2017). “ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data.” Methods in Ecology and Evolution, 8, 28-36. doi: 10.1111/2041-210X.12628.

See Also

motifStack::motifStack(), compare_motifs(), ggtree::ggtree(), ggplot2::ggplot(), motif_tree2()

Examples

jaspar <- read_jaspar(system.file("extdata", "jaspar.txt",
                                  package = "universalmotif"))
if (requireNamespace("ggtree", quietly = TRUE)) {
jaspar.tree <- motif_tree(jaspar, linecol = "none", labels = "name",
                          layout = "rectangular")
}

## Not run: 
## When inputting a dist object, the linecol and tipsize options are
## not available. To add these manually:

library(MotifDb)
library(ggtree)
library(ggplot2)

motifs <- filter_motifs(MotifDb, organism = "Athaliana")[1:50]
comparison <- compare_motifs(motifs, method = "PCC", score.strat = "a.mean")
comparison <- as.dist(1 - comparison)
mot.names <- attr(comparison, "Labels")
tree <- motif_tree(comparison)

annotations <- data.frame(label = mot.names,
                          icscore = sapply(motifs, function(x) x["icscore"]),
                          family = sapply(motifs, function(x) x["family"]))

tree <- tree %<+% annotations +
          geom_tippoint(aes(size = icscore)) +
          aes(colour = family) +
          theme(legend.position = "right",
                legend.title = element_blank())

## End(Not run)

---

Generate ggtree-based motif trees using compare_motifs2().

Description

motif_tree2() is the leaner counterpart of motif_tree(): it builds the distance matrix via compare_motifs2() (mean Pearson correlation, with built-in p-value / q-value machinery) instead of compare_motifs(). The distance for hclust() is derived from the mean Pearson correlation matrix as (1 - score) / 2, mapping similarity in ⁠[-1, 1]⁠ to a symmetric distance in ⁠[0, 1]⁠. The tree-construction and rendering steps (hclust() followed by ape::as.phylo() followed by ggtree::ggtree()) match motif_tree() exactly, so the visual conventions, the linecol / labels / tipsize slot-mapping arguments, and the dist-input short-circuit all behave the same.

Usage

motif_tree2(motifs, layout = "circular", linecol = "family",
  labels = "none", tipsize = "none", legend = TRUE,
  branch.length = "none", min.overlap = 6, tryRC = TRUE, nthreads = 1,
  progress = FALSE, ...)

Arguments

motifs	list or dist. See convert_motifs() for accepted motif formats. Alternatively, a pre-built dist object (e.g. the result of as.dist((1 - compare_motifs2(...)) / 2)) skips the comparison step entirely and goes straight to tree construction.
layout	character(1). One of c('rectangular', 'slanted', 'fan', 'circular', 'radial', 'equal_angle', 'daylight'). Passed through to ggtree::ggtree().
linecol	character(1). Motif slot used to colour the branches (e.g. "family"). Not used when motifs is a dist.
labels	character(1). Motif slot used to label the tips (e.g. "name"). For dist input, only "name" is meaningful.
tipsize	character(1). Motif slot used to size the tips (e.g. "icscore"). Not used when motifs is a dist.
legend	logical(1). Show legend for line colour and tip size.
branch.length	character(1). If "none", draw a cladogram. Passed through to ggtree::ggtree().
min.overlap	integer(1). Minimum overlap in columns for the pairwise alignment to count. Forwarded to compare_motifs2(). Default 6.
tryRC	logical(1). Also test reverse-complement alignments. Forwarded to compare_motifs2(). Default TRUE.
nthreads	numeric(1). Threads passed to compare_motifs2(). nthreads = 0 uses all available threads.
progress	logical(1). Print progress messages.
...	Additional arguments passed through to ggtree::ggtree().

Details

The PCC-to-distance conversion d = (1 - score) / 2 is symmetric (the comparison matrix from compare_motifs2() is symmetric in matrix mode), so as.dist() and hclust() accept it directly. Anti-correlated motif pairs land at the maximum distance of 1; identical motifs land at 0.

For more control over tree construction, run compare_motifs2() directly with matrix.out = "score", convert to a dist object yourself, and pass that dist to motif_tree2(). The function detects the dist input and skips the comparison step.

Value

A ggplot (ggtree) object.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

motif_tree(), compare_motifs2(), merge_similar2(), ggtree::ggtree()

Examples

## Not run: 
library(universalmotif)
jaspar <- read_jaspar(system.file("extdata", "jaspar.txt",
                                  package = "universalmotif"))
if (requireNamespace("ggtree", quietly = TRUE)) {
  motif_tree2(jaspar, linecol = "none", labels = "name",
              layout = "rectangular")
}

## Equivalent two-step form, useful when you want to inspect the
## distance matrix before drawing the tree:
if (requireNamespace("ggtree", quietly = TRUE)) {
  score.mat <- compare_motifs2(jaspar, matrix.out = "score")
  d <- as.dist((1 - score.mat) / 2)
  motif_tree2(d, labels = "name")
}

## End(Not run)

---

Visualise composition match between target and background sequences.

Description

Companion to match_bkg(). Overlays density curves for the target and matched-background sets in a faceted ggplot, useful for visually verifying that the binned matching landed where it was supposed to. Called with two XStringSets it plots GC fraction and sequence length; called with the indices data.frame from match_bkg(..., return.indices = TRUE) it plots every matched axis, including any covariates.

Usage

plot_match_bkg(sequences = NULL, bkg = NULL, by = NULL, bins = 30L,
  indices = NULL)

Arguments

sequences	XStringSet or NULL. The target sequences passed to match_bkg(). Ignored when indices is supplied.
bkg	XStringSet or NULL. The matched background sequences returned by match_bkg(). Ignored when indices is supplied.
by	character or NULL. Axes to plot. For the sequences/bkg form, a subset of c("gc", "length") (default both). For the indices form, a subset of the available axis names ("gc", "length", and each covariate name); NULL (default) plots them all.
bins	integer(1). Ignored when geom = "density"; kept for API parity with other v2 plotting helpers. Default 30L.
indices	data.frame or NULL. The output of match_bkg(..., return.indices = TRUE). When supplied, every ⁠target.<axis>⁠ / ⁠universe.<axis>⁠ column pair (apart from the index column) becomes a density panel, so covariate matches can be inspected alongside GC and length. Default NULL.

Value

A ggplot object.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

match_bkg()

---

Plot the position distribution of motif hits.

Description

Companion to motif_peaks(). Returns a ggplot faceted by motif, showing a histogram of best-hit centre positions and a shaded rectangle marking the most-enriched window.

Usage

plot_motif_peaks(peaks, motifs = NULL, ncol = 2L, bins = 50L,
  fill.window = "#ffe4a3")

Arguments

peaks	data.frame returned by motif_peaks(). Must carry the centers list-column.
motifs	character or NULL. Subset to these motifs (by name). NULL (default) plots every row of peaks.
ncol	integer(1). Number of facet columns. Default 2.
bins	integer(1). Histogram bin count. Default 50.
fill.window	character(1). Colour of the best-window shade. Default "#ffe4a3" (a pale yellow).

Value

A ggplot object.

Author(s)

Benjamin Jean-Marie Tremblay, [email protected]

See Also

motif_peaks()

---