Introduction

A sequence motif is a short sequence conserved in multiple amino acid or nucleic acid sequences with biological significance such as DNA binding sites, where transcription factors (TFs) bind and catalytic sites in enzymes. Motifs are often represented as position count matrix (PCM)(Table 1), position frequency matrix (PFM)(Table 2) or position weight matrix (PWM) (reviewed in^1–3). Sequence logos, based on information theory or probability theory, have been used as a graphical representation of sequence motifs to depict the frequencies of bases or amino acids at each conserved position^4–8.

Table 1: Sample of PCM. Each row represents the counts of A, C, G and T, while each column is the counts for that nucleotide position.
A	1	1	23	22	1	2
C	19	16	0	0	2	4
G	2	4	0	1	19	15
T	1	2	0	0	1	2

Table 2: Sample of PFM. Each row represents the frequency of A, C, G and T, while each column is the frequency for that nucleotide position.
A	0.04	0.04	1	0.96	0.04	0.09
C	0.83	0.70	0	0.00	0.09	0.17
G	0.09	0.17	0	0.04	0.83	0.65
T	0.04	0.09	0	0.00	0.04	0.09

High-throughput experimental approaches have been developed to identify binding sequences for a large number of TFs in several organisms. The experimental methods include protein binding microarray (PBM)⁹, bacterial one-hybrid (B1H)¹⁰, systematic evolution of ligands by exponential enrichment (SELEX)^11,12, High-Throughput (HT) SELEX¹³, and selective microfluidics-based ligand enrichment followed by sequencing (SMiLE-seq)¹⁴. Several algorithms have been developed to discover motifs from the binding sequences, identified by the above experimental methods. For example, MEME Suite¹⁵ has been widely used to construct motifs for B1H and SELEX data, while Seed-and-Wobble (SW)¹⁶, Dream5¹⁷, and BEEML¹⁸ have been used to generate motifs for PBM data.

Visualization tools have been developed for representing a sequence motif as individual sequence logo, such as weblogo⁷and seqlogo¹⁹. To facilitate the comparison of multiple motifs, several motif alignment and clustering tools have been developed such as Tomtom²⁰, MatAlign²¹, and STAMP²². However, these tools are lacking flexibility in dispalying the similarities or differences among motifs.

The motifStack package²³ is designed for for the graphic representation of multiple motifs with their similarities or distances depicted. It provides the flexibility for users to customize the clustering and graphic parameters such as distance cutoff, the font family, colors, symbols, layout styles, and ways to highlight different regions or branches.

In this workshop, we will demonstrate the features and flexibility of motifStack for visualizing the alignment of multiple motifs in various styles. In addition, we will illustrate the utility of motifStack for providing insights into families of related motifs using several large collections of homeodomain (HD) DNA binding motifs from human, mouse and fly. Furthermore, we will show how motifStack facilitated the discovery that changing the computational method used to build motifs from the same binding data can lead to artificial segregation of motifs for identical TFs using the same motif alignment method.

Installation

BiocManager is used to install the released version of motifStack.

library(BiocManager)
install("motifStack")

To install the development version of motifStack, please try

install("jianhong/motifStack")

packageVersion("motifStack")

## [1] '1.29.5'

Please note that starting from version 1.28.0, motifStack requires cario (>=1.6) instead of ghostscript which is for a lower version of motifStack. The capabilities function could be used to check capablilities of cario for current R session. If motifStack version number is smaller than 1.28.0, the full path to the executable ghostscript must be available for R session.

In this tutorial, we will illustrate the functionalities of motifStack available in version >=1.29.4.

capabilities()["cairo"]

## cairo 
##  TRUE

Inputs of motifStack

The inputs of motifStack are one or more motifs represented as PCM or PFM. Motifs in PFM format are available in most of the motif databases such as TRANSFAC²⁴, JASPAR²⁵, CIS-BP²⁶ and HOCOMOCO²⁷. More databases could be found at http://meme-suite.org/db/motifs.

In addition, MotifDb package contains a collection of DNA-binding sequence motifs in PFM format. 9933 PFMs were collected from 11 sources: cisbp_1.02²⁶, FlyFactorSurvey²⁸, HOCOMOCO²⁷, HOMER²⁹, hPDI³⁰, JASPAR²⁵, Jolma2013³¹, ScerTF³², stamlab³³, SwissRegulon³⁴ and UniPROBE⁹.

library(MotifDb)
MotifDb[1]

## MotifDb object of length 1
## | Created from downloaded public sources: 2013-Aug-30
## | 1 position frequency matrices from 1 source:
## |         cisbp_1.02:    1
## | 1 organism/s
## |        Scerevisiae:    1
## Scerevisiae-cisbp_1.02-M0001_1.02

MotifDb[[1]]

##             1           2           3           4           5           6
## A 0.009615385 0.009615385 0.971153846 0.009615385 0.009615385 0.971153846
## C 0.009615385 0.009615385 0.009615385 0.009615385 0.971153846 0.009615385
## G 0.009615385 0.009615385 0.009615385 0.009615385 0.009615385 0.009615385
## T 0.971153846 0.971153846 0.009615385 0.971153846 0.009615385 0.009615385
##             7           8
## A 0.009615385 0.009615385
## C 0.971153846 0.009615385
## G 0.009615385 0.009615385
## T 0.009615385 0.971153846

The plotMotifLogo function can be used to plot sequence logo for motifs in PFM format (Figure 1). As default, the nucleotides in a sequence logo will be drawn proportional to their information content (IC)⁸.

library(motifStack)
plotMotifLogo(MotifDb[[1]])

Figure 1: Plot sequence logo for motif in PFM format.

A PFM can be formated as an pfm object which is the preferred format for motifStack (Figure 2).

matrix.bin <- query(MotifDb, "bin_SANGER")
motif <- new("pfm", mat=matrix.bin[[1]], name=names(matrix.bin)[1])
motifStack(motif)

Figure 2: Plot sequence logo for a pfm object.

Similary, a PCM can be converted to a pcm object easily (Figure 3).

path <- system.file("extdata", package = "motifStack")
pcm <- read.table(file.path(path, "bin_SOLEXA.pcm"))
pcm

##   V1 V2  V3   V4   V5   V6   V7  V8   V9
## 1  A  | 462    0 1068 1025 1068   0 1019
## 2  C  |  71   60    0   24    0 993    0
## 3  G  | 504    0    0    0    0  12    0
## 4  T  |  31 1008    0   19    0  63   49

pcm <- pcm[,3:ncol(pcm)]
rownames(pcm) <- c("A","C","G","T") # must have rownames
motif <- new("pcm", mat=as.matrix(pcm), name="bin_SOLEXA")
plot(motif)

Figure 3: Plot sequence logo for a pcm object.

In addition, importMatrix function can be used to import the motifs into a pcm or pfm object from files exported by TRANSFAC²⁴, JASPAR²⁵ or CIS-BP²⁶.

motifs <- importMatrix(dir(path, "*.pcm", full.names = TRUE))
motifStack(motifs)

Figure 4: Plot sequence logo for data imported by importMatrix

Plot single sequence logo

motifStack can be used to plot DNA/RNA sequence logo, affinity logo and amino acid sequence log. It provides the flexibility for users to customize the graphic parameters such as the font family, symbol colors and markers. Different from most of the sequence logo plotting tools, motifStack does not need to pre-define the path for plotting. It depends on grImport/2³⁵ package to import fonts into a path, which provides the power of motifStack more flexibility.

Plot a DNA/RNA/amino acid sequence logo

To plot an RNA sequence logo, simply change the rowname of the matrix from “T” to “U” (Figure 5).

rna <- pcm
rownames(rna)[4] <- "U"
RNA_motif <- new("pcm", mat=as.matrix(rna), name="RNA_motif")
plot(RNA_motif)

Figure 5: Plot an RNA sequence logo.

If the PFM of amino acid sequence motif is provided, motifStack can be used to draw amino acid sequence logo (Figure 6).

protein<-read.table(file.path(path,"cap.txt"))
protein<-t(protein[,1:20])
protein<-pcm2pfm(protein)
protein_motif<-new("pfm", mat=protein, name="CAP", 
                   color=colorset(alphabet="AA",colorScheme="chemistry"))
plot(protein_motif)

Figure 6: Plot amino acid sequence logo.

motifStack can also be used to plot other types of logos (Figure 7).

m <- matrix(0, nrow = 10, ncol = 10, 
            dimnames = list(strsplit("motifStack", "")[[1]],
                            strsplit("motifStack", "")[[1]]))
for(i in seq.int(10)) m[i, i] <- 1
ms <- new("pfm", mat=m, name="motifStack")
plot(ms)

Figure 7: Plot a user-defined logo.

The sequence logo can be visulized with different fonts and colors (Figure 8).

motif$color <- colorset(colorScheme = 'blindnessSafe')
plot(motif, font="mono,Courier", fontface="plain")

Figure 8: Plot sequence logo with different fonts and colors.

Plot sequence logo with specified regions highlighted using the markers slot

The markers slot of a pfm or pcm object can be assigned by a list of marker objects. There are three type of markers, “rect”, “line” and “text”. If markers slot is assigned, markers can be plotted (Figure 9).

markerRect <- new("marker", type="rect", start=6, stop=7, 
                  gp=gpar(lty=2, fill=NA, col="orange", lwd=3))
markerRect2 <- new("marker", type="rect", start=2, 
                   gp=gpar(lty=3, fill=NA, col="green", lwd=2))
markerLine <- new("marker", type="line", start=2, stop=7,
                  gp=gpar(lwd=2, col="darkblue"))
markerText <- new("marker", type="text", start=c(1, 5), 
                  label=c("*", "core"), gp=gpar(cex=2, col="red"))
motif <- new("pcm", mat=as.matrix(pcm), name="bin_SOLEXA", 
             markers=c(markerRect, markerLine, markerText, markerRect2))
plot(motif)

Figure 9: Plot a sequence logo with selected regions highlighted.

Plot motifs with ggplot2

Sequence logo can also be plotted by geom_motif function which is compatible with ggplot2³⁶ (Figure 10).

len <- length(motifs)
motifs[[1]]$markers <- list(markerRect2, markerRect)
df <- data.frame(x=.5, y=(seq.int(len)-.5)/len, 
                 width=.75, height=1/(len+1))
df$motif <- motifs
library(ggplot2)
ggplot(df, aes(x=x, y=y, width=width, height=height, motif=motif)) +
  geom_motif(use.xy = TRUE) + theme_bw() + xlim(0, 1) + ylim(0, 1)

Figure 10: Plot sequence logo with ggplot2.

Plot multiple sequence logos

motifStack is designed to visulize the alignment of multiple motifs and facilitate the analysis of binding site diversity/conservation within families of TFs and evolution of binding sites among different species. There are many functions in motifStack such as plotMotifLogoStack, motifCloud, plotMotifStackWithRadialPhylog, motifPiles, and motifCircos for plotting motifs in different layouts, motifSignature for generating motif signatures by merging similar motifs. To make it easy to use, we created one workflow function called motifStack, which can be used to draw a single DNA, RNA or amino acid motif as a sequence logo, as well as to align motifs, generate motif signatures and plot aligned motifs as a stack, phylogeny tree, radial tree or cloud of sequence logos.

To align motifs, DNAmotifAlignment function implements Ungapped Smith–Waterman (local) alignment based on column comparison scores calculated by Kullback-Leibler distance³⁷. DNAmotifAlignment requires the inputs be sorted by the distance of motifs. To get the sorted motifs for alignment, distances can be calculated using motifDistance in motIV (R implementation of STAMP²²). Distances generated from MatAlign²¹ are also supported with a helper function ade4::newick2phylog³⁸.

Plot sequence logo stack

As shown in above, motifStack function can be used to plot motifs as a stack of sequence logos easily. By setting the layout to “tree” (Figure 11) or “phylog” (Figure 12), motifStack will plot the stack of sequence logos with a linear tree.

motifStack(motifs, layout = "tree")

Figure 11: Plot sequence logo stack with tree layout.

motifStack(motifs, layout="phylog", f.phylog=.15, f.logo=0.25)

Figure 12: Plot sequence logo stack with phylog layout.

Even though the layout of “phylog” allows more motifs to be plotted comparing to the layout of “tree”, both layouts only allow small sets of motifs to be readily compared. However, moderate throughput analysis of TFs has generated 100s to 1000s of motifs, which make the visualization of their relationships more challenging. To reduce the total number of motifs depicted, clusters of highly similar motifs can be merged with the motifSignature function using a distance threshold. The distance threshold is user-determined; a low threshold will merge only extremely similar motifs while a higher threshold can group motifs with lower similarity.

matrix.fly <- query(MotifDb, "FlyFactorSurvey")
motifs2 <- as.list(matrix.fly)
## format the name
names(motifs2) <- gsub("(_[\\.0-9]+)*_FBgn\\d+$", "", 
                       elementMetadata(matrix.fly)$providerName)
names(motifs2) <- gsub("[^a-zA-Z0-9]", "_", names(motifs2))

motifs2 <- motifs2[unique(names(motifs2))]
set.seed(1)
pfms <- sample(motifs2, 50)

## use MotIV to calculate the distances of motifs
jaspar.scores <- MotIV::readDBScores(file.path(find.package("MotIV"), 
                                               "extdata", 
                                               "jaspar2010_PCC_SWU.scores"))
d <- MotIV::motifDistances(lapply(pfms, pfm2pwm))
hc <- MotIV::motifHclust(d, method="average")
## convert the hclust to phylog object
phylog <- hclust2phylog(hc)
## reorder the pfms by the order of hclust
leaves <- names(phylog$leaves)
pfms <- pfms[leaves]
## create a list of pfm objects
pfms <- mapply(pfms, names(pfms), 
               FUN=function(.pfm, .name){
                 new("pfm",mat=.pfm, name=.name)})
## extract the motif signatures
motifSig <- motifSignature(pfms, phylog, groupDistance=0.01, min.freq=1)

Logos for these merged motifs can be plotted on the same dendrograms, using the relative size of the logos to reflect the number of TFs that contribute to the merged motif. Compared with the original dendrogram, the number of TFs with similar motifs and the relationships between the merged motifs can be more easily compared. The 50 motifs shown in the original image are now shown as only 28 motif clusters/signatures, with the members of each cluster clearly indicated. By setting the colors of dendrogram via parameter col.tree, the colors of lable of leaves via parameter col.leaves, and the colors of annotation bar via parameter col.anno, motifPiles function can mark each motif with multiple annotation colors (Figure 13).

## get the signatures from object of motifSignature
sig <- signatures(motifSig)
## get the group color for each signature
gpCol <- sigColor(motifSig)

library(RColorBrewer)
color <- brewer.pal(12, "Set3")
## plot the logo stack with pile style.
motifPiles(phylog=phylog, pfms=pfms, pfms2=sig, 
            col.tree=rep(color, each=5),
            col.leaves=rep(rev(color), each=5),
            col.pfms2=gpCol, 
            r.anno=c(0.02, 0.03, 0.04), 
            col.anno=list(sample(colors(), 50), 
                          sample(colors(), 50), 
                          sample(colors(), 50)),
            motifScale="logarithmic",
            plotIndex=TRUE,
            groupDistance=0.01)

Figure 13: Plot sequence logo stack by motifPiles.

Plot sequence logo cloud

To emphasize the weight of motif signatures, motif signatures can be plotted as a circular (Figure 14) or rectanglular word clouds with motif size representing the number of motifs that contributed to the signature. The larger the motif size, the larger the number of motifs within that motif signature.

## draw the motifs with a word-cloud style.
motifCloud(motifSig, scale=c(6, .5), layout="cloud")

Figure 14: Plot sequence logo cloud.

Plot sequence logo circos

The sequence logo cloud will loose the cluster information. Plotting motifs as sequence logo circos is a good solution to visualize motif signatures with dendrogram (Figure 15).

## plot the logo stack with radial style.
plotMotifStackWithRadialPhylog(phylog=phylog, pfms=sig, 
                              circle=0.4, cleaves = 0.3, 
                              clabel.leaves = 0.5, 
                              col.bg=rep(color, each=5), col.bg.alpha=0.3, 
                              col.leaves=rep(rev(color), each=5),
                              col.inner.label.circle=gpCol, 
                              inner.label.circle.width=0.03,
                              angle=350, circle.motif=1.2, 
                              motifScale="logarithmic")

Figure 15: Plot sequence logo circos.

motifCircos function provide the possibility to plot the motif signatures with their original sequence logos in the circos style (Figure 16).

## plot the logo stack with cirsoc style.
motifCircos(phylog=phylog, pfms=pfms, pfms2=sig, 
            col.tree.bg=rep(color, each=5), col.tree.bg.alpha=0.3, 
            col.leaves=rep(rev(color), each=5),
            col.inner.label.circle=gpCol, 
            inner.label.circle.width=0.03,
            col.outer.label.circle=gpCol, 
            outer.label.circle.width=0.03,
            r.rings=c(0.02, 0.03, 0.04), 
            col.rings=list(sample(colors(), 50), 
                           sample(colors(), 50), 
                           sample(colors(), 50)),
            angle=350, motifScale="logarithmic")

Figure 16: Plot sequence logo circos by motifCiros.

Plot fly data from different experimental methods

Following scripts will show how motifStack facilitated the discovery that changing the experiemntal method used to build motifs from the same binding data can lead to artificial segregation of motifs for identical TFs within a motif alignment. MatAlign²¹ was downloaded from http://stormo.wustl.edu/MatAlign/ and neighbor is a tool of phylip. MatAlign2tree.pl is a perl script for running MatAlign and neighbor to get the neighbor-joint phylogenetic tree for given motifs. Compare to STAMP, MatAlign alignment method is superior for discriminating between the closely related motifs²³.

Here we clustered motifs for a smaller set of fly HDs that were determined using four different experimental methods to measure DNA binding specificity, B1H, HT-SELEX, SMiLE-seq and PBM, to confirm that the motif generation method can impact clustering. As it shown in Figre 17 B1H motifs from different lab tend to cluster with each other (motifs 1,2, 37,38, 41-113, 132,133, 136-162, 166-171, 180-195, 206-236, 242,243) compared to PBM data (motifs ‘eve’ in blue). PMB data (motif 5-13, 18-21, 23-26) tend to have lower IC, which may be introduced by computation methods²³. HT-SELEX tend to have dimer (motif labels in red).

path <- system.file('libs', package = 'motifStackWorkshop')
pfmpath <- system.file('extdata', "pfms", package = 'motifStackWorkshop')
MatAlign2tree <- system.file('extdata', 'MatAlign2tree.pl', 
                             package = 'motifStackWorkshop')


pfmfiles <- dir(pfmpath, full.names = TRUE)
pfms <- lapply(pfmfiles, importMatrix, format="pfm", to="pfm")
pfms <- unlist(pfms)
## fix the name issue
names(pfms) <- gsub("\\.", "_", names(pfms))
pfms <- sapply(pfms, function(.ele){
  .ele@name <-gsub("\\.", "_", .ele@name)
  .ele
})
## get phylog
outpath="output"
matalign_path <- file.path(path, "matalign")
neighbor_path <- file.path(path, "neighbor")
system(paste("perl", MatAlign2tree, "--in", getwd(), "--pcmpath", pfmpath, 
             "--format pfm --ext pfm",
             "--out", outpath, "--matalign", matalign_path, 
             "--neighbor", neighbor_path, "--tree","UPGMA"))
newickstrUPGMA <- readLines(con=file.path(outpath, "NJ.matalign.distMX.nwk"))
unlink(outpath, recursive = TRUE)
phylog <- newick2phylog(newickstrUPGMA, FALSE)
leaves <- names(phylog$leaves)
motifs <- pfms[leaves]
groupDistance <- 2.5
trim <- 0.2
motifSig <- motifSignature(motifs, phylog, groupDistance=groupDistance, 
                           min.freq=1, trim=trim)
sig <- signatures(motifSig)
gpCol <- sigColor(motifSig)

##set methods color
colorSet <- c("B1H_Mathelier"="#40E0D0", "B1H_Zhu"="#008080", 
              "PBM_Busser"="brown", "PBM_Weirauch"="#F69156",
              "SELEX_Nitta"="blue", 
              "SMiLE_seq_Isakova"="#D900D9")
methods <- factor(sub("^.*?__", "", leaves))
lev <- levels(methods)
levels(methods) <- colorSet[lev]
methods <- as.character(methods)

leaveNames <- gsub("__.*$", "", leaves)
leaveCol <- ifelse(leaveNames %in% c("unc_4", "hbn", "bsh", "Optix", 
                                     "al", "PHDP", "CG32532", "CG11294"),
                   "#D55E00", 
                   ifelse(leaveNames == "eve", "#0072B2", "#000000"))

## calculate average of top 8 position information content for each motif
icgp <- sapply(sapply(motifs, getIC), function(.ele)
  mean(sort(.ele, decreasing=TRUE)[1:min(8, length(.ele))]))
icgp.ranges <- range(icgp)
icgp <- cut(icgp, 10, labels=colorRampPalette(c("#009E73", "#000000", "#D55E00"))(10))
icgp.image <- as.raster(matrix(rev(levels(icgp)), ncol=1))
icgp <- as.character(icgp)

motifCircos(phylog=phylog, pfms=DNAmotifAlignment(motifs), 
            col.tree.bg=methods, col.tree.bg.alpha=.3, 
            r.rings=c(.05, .1), col.rings=list(icgp, rep("white", length(icgp))), 
            col.inner.label.circle=gpCol, inner.label.circle.width=0.03,
            labels.leaves=leaveNames, col.leaves=leaveCol, 
            cleaves=.1, r.tree=1.3, r.leaves=.2, clabel.leaves=.7, 
            motifScale="logarithmic", angle=358, groupDistance=groupDistance,
            plotIndex=TRUE, IndexCex=.7, groupDistanceLineCol='#D55E00')

legend(-2.4, 2.4, legend=names(colorSet), 
       fill=highlightCol(colorSet, .3), 
       border="white", lty=NULL, bty = "n")

text(2.05, 2.3, labels="Information Content")
rasterImage(icgp.image, 2, 1.8, 2.04, 2.2, interpolate = FALSE)
text(c(2.05, 2.05), c(1.8, 2.2), labels=format(icgp.ranges, digits=3), adj=0)

Figure 17: Fly data from different methods.

References

1. Das, M. K. & Dai, H.-K. A survey of dna motif finding algorithms. in BMC bioinformatics 8, S21 (BioMed Central, 2007).

2. Sandve, G. K., Abul, O., Walseng, V. & Drabløs, F. Improved benchmarks for computational motif discovery. BMC bioinformatics 8, 193 (2007).

3. Tompa, M. et al. Assessing computational tools for the discovery of transcription factor binding sites. Nature biotechnology 23, 137 (2005).

4. Schneider, T. D. & Stephens, R. M. Sequence logos: A new way to display consensus sequences. Nucleic acids research 18, 6097–6100 (1990).

5. Schneider, T. D. in Methods in enzymology 274, 445–455 (Elsevier, 1996).

6. Schneider, T. D. Strong minor groove base conservation in sequence logos implies dna distortion or base flipping during replication and transcription initiation. Nucleic acids research 29, 4881–4891 (2001).

7. Crooks, G. E., Hon, G., Chandonia, J.-M. & Brenner, S. E. WebLogo: A sequence logo generator. Genome research 14, 1188–1190 (2004).

8. Workman, C. T. et al. EnoLOGOS: A versatile web tool for energy normalized sequence logos. Nucleic acids research 33, W389–W392 (2005).

9. Newburger, D. E. & Bulyk, M. L. UniPROBE: An online database of protein binding microarray data on protein–dna interactions. Nucleic acids research 37, D77–D82 (2008).

10. Meng, X., Brodsky, M. H. & Wolfe, S. A. A bacterial one-hybrid system for determining the dna-binding specificity of transcription factors. Nature biotechnology 23, 988 (2005).

11. Ellington, A. D. & Szostak, J. W. In vitro selection of rna molecules that bind specific ligands. nature 346, 818 (1990).

12. Tuerk, C. & Gold, L. Systematic evolution of ligands by exponential enrichment: RNA ligands to bacteriophage t4 dna polymerase. Science 249, 505–510 (1990).

13. Hoinka, J. et al. Large scale analysis of the mutational landscape in ht-selex improves aptamer discovery. Nucleic acids research 43, 5699–5707 (2015).

14. Isakova, A. et al. SMiLE-seq identifies binding motifs of single and dimeric transcription factors. Nature methods 14, 316 (2017).

15. Bailey, T. L. et al. MEME suite: Tools for motif discovery and searching. Nucleic acids research 37, W202–W208 (2009).

16. Badis, G. et al. Diversity and complexity in dna recognition by transcription factors. Science 324, 1720–1723 (2009).

17. Weirauch, M. T. et al. Evaluation of methods for modeling transcription factor sequence specificity. Nature biotechnology 31, 126 (2013).

18. Zhao, Y. & Stormo, G. D. Quantitative analysis demonstrates most transcription factors require only simple models of specificity. Nature biotechnology 29, 480 (2011).

19. Bembom, O. SeqLogo: Sequence logos for dna sequence alignments. R package version 1,

20. Gupta, S., Stamatoyannopoulos, J. A., Bailey, T. L. & Noble, W. S. Quantifying similarity between motifs. Genome biology 8, R24 (2007).

21. Zhao, G. et al. Conserved motifs and prediction of regulatory modules in caenorhabditis elegans. G3: Genes, Genomes, Genetics 2, 469–481 (2012).

22. Mahony, S. & Benos, P. V. STAMP: A web tool for exploring dna-binding motif similarities. Nucleic acids research 35, W253–W258 (2007).

23. Ou, J., Wolfe, S. A., Brodsky, M. H. & Zhu, L. J. MotifStack for the analysis of transcription factor binding site evolution. Nature methods 15, 8 (2018).

24. Wingender, E. et al. TRANSFAC: An integrated system for gene expression regulation. Nucleic acids research 28, 316–319 (2000).

25. Sandelin, A., Alkema, W., Engström, P., Wasserman, W. W. & Lenhard, B. JASPAR: An open-access database for eukaryotic transcription factor binding profiles. Nucleic acids research 32, D91–D94 (2004).

26. Weirauch, M. T. et al. Determination and inference of eukaryotic transcription factor sequence specificity. Cell 158, 1431–1443 (2014).

27. Kulakovskiy, I. V. et al. HOCOMOCO: Towards a complete collection of transcription factor binding models for human and mouse via large-scale chip-seq analysis. Nucleic acids research 46, D252–D259 (2017).

28. Zhu, L. J. et al. FlyFactorSurvey: A database of drosophila transcription factor binding specificities determined using the bacterial one-hybrid system. Nucleic acids research 39, D111–D117 (2010).

29. Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and b cell identities. Molecular cell 38, 576–589 (2010).

30. Xie, Z., Hu, S., Blackshaw, S., Zhu, H. & Qian, J. HPDI: A database of experimental human protein–dna interactions. Bioinformatics 26, 287–289 (2009).

31. Jolma, A. et al. DNA-binding specificities of human transcription factors. Cell 152, 327–339 (2013).

32. Spivak, A. T. & Stormo, G. D. ScerTF: A comprehensive database of benchmarked position weight matrices for saccharomyces species. Nucleic acids research 40, D162–D168 (2011).

33. Neph, S. et al. Circuitry and dynamics of human transcription factor regulatory networks. Cell 150, 1274–1286 (2012).

34. Pachkov, M., Erb, I., Molina, N. & Van Nimwegen, E. SwissRegulon: A database of genome-wide annotations of regulatory sites. Nucleic acids research 35, D127–D131 (2006).

35. Murrell, P. & others. Importing vector graphics: The grImport package for r. Journal of Statistical Software 30, 1–37 (2009).

36. Wickham, H. Ggplot2: Elegant graphics for data analysis. (Springer, 2016).

37. Roepcke, S., Grossmann, S., Rahmann, S. & Vingron, M. T-reg comparator: An analysis tool for the comparison of position weight matrices. Nucleic acids research 33, W438–W441 (2005).

38. Dray, S., Dufour, A.-B. & others. The ade4 package: Implementing the duality diagram for ecologists. Journal of statistical software 22, 1–20 (2007).