Introduction

A sequence logo, based on information theory, has been widely used as a graphical representation of sequence conservation (aka motif) in multiple amino acid or nucleic acid sequences. Sequence motif represents conserved characteristics such as DNA binding sites, where transcription factors bind, and catalytic sites in enzymes. Although many tools, such as seqlogo1, have been developed to create sequence motif and to represent it as individual sequence logo, software tools for depicting the relationship among multiple sequence motifs are still lacking. We developed a flexible and powerful open-source R/Bioconductor package, motifStack, for visualization of the alignment of multiple sequence motifs.

Import matrix

The importMatrix function is used to import motifs from files or convert XMatrix/XMatrixList object into motifStack compatable format.

convert motifs from XMatrixList

library(motifStack)
library(JASPAR2020)
motifs <- importMatrix(getMatrixSet(JASPAR2020, 
                                    list(species="Mus musculus")))
plot(motifs[[1]])
## Loading required namespace: Cairo
import data from PFMatrixList

import data from PFMatrixList

import motifs from files

The supported formats are “meme”, “transfac”, “jaspar”, “scpd”, “cisbp”, and “psam”.

RUNX1 <- importMatrix(system.file("extdata", "MA0002.1.jaspar",
                                  package = "motifStack",
                                  mustWork = TRUE))[[1]]
plot(RUNX1)
import data from file

import data from file

Examples of using motifStack

plot a DNA sequence logo with different fonts and colors

Users can select different fonts and colors to draw the sequence logo.

library(motifStack)
pcm <- read.table(file.path(find.package("motifStack"), 
                            "extdata", "bin_SOLEXA.pcm"))
pcm <- pcm[,3:ncol(pcm)]
rownames(pcm) <- c("A","C","G","T")
motif <- new("pcm", mat=as.matrix(pcm), name="bin_SOLEXA")
##pfm object
#motif <- pcm2pfm(pcm)
#motif <- new("pfm", mat=motif, name="bin_SOLEXA")
plot(motif)
Plot a DNA sequence logo with different fonts and colors

Plot a DNA sequence logo with different fonts and colors

#plot the logo with same height
plot(motif, ic.scale=FALSE, ylab="probability")
Plot a DNA sequence logo with different fonts and colors

Plot a DNA sequence logo with different fonts and colors

#try a different font and a different color group
motif@color <- colorset(colorScheme='basepairing')
plot(motif,font="serif")
Plot a DNA sequence logo with different fonts and colors

Plot a DNA sequence logo with different fonts and colors

plot sequence logo with markers

If you assign markers slot by a list of marker object, markers can be plotted in the figure. There are three type of markers, “rect”, “line” and “text”.

markerRect <- new("marker", type="rect", start=6, stop=7, gp=gpar(lty=2, fill=NA, col="orange"))
markerLine <- new("marker", type="line", start=2, stop=7, gp=gpar(lwd=2, col="red"))
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))
plot(motif)
Plot a DNA sequence logo with markers

Plot a DNA sequence logo with markers

change the x-axis labels

plot(motif, xaxis=paste0("pos", seq.int(7)+10))
Plot a DNA sequence logo with pre-defined xlabels

Plot a DNA sequence logo with pre-defined xlabels

plot sequence logo stack

To show multiple motifs on the same canvas as a sequence logo stack, the distance of motifs need to be calculated first. Previously, MotIV3::motifDistances ( R implementation of STAMP4) is used to calculate the distance. However, The MotIV package were dropped from Bioconductor 3_12. Currently, by default, R implementation of matalign is used. After alignment, users can use plotMotifLogoStack, plotMotifLogoStackWithTree or plotMotifStackWithRadialPhylog to draw sequence logos in different layouts. To make it easy to use, we integrated different functionalities into one workflow function named as motifStack.

library(motifStack)
#####Input#####
motifs<-importMatrix(dir(file.path(find.package("motifStack"),
                                   "extdata"),"pcm$", 
                         full.names = TRUE))

## plot stacks
motifStack(motifs, layout="stack", ncex=1.0)
Plot motifs with sequence logo stack style

Plot motifs with sequence logo stack style

rnaMotifs <- DNAmotifToRNAmotif(motifs)
names(rnaMotifs)
## [1] "bin_SOLEXA"   "fd64A_SOLEXA" "fkh_NAR"      "foxo_SOLEXA"  "FoxP_SOLEXA" 
## [6] "slp1_SOLEXA"  "slp2_SOLEXA"
motifStack(rnaMotifs, layout = "stack", 
           reorder=FALSE) ## we can also use reorder=FALSE to keep the order of input. 
Plot RNA motifs with sequence logo stack style

Plot RNA motifs with sequence logo stack style

motif2 <- motif
motif2$mat <- motif$mat[, 5:12]
motif2$name <- "logo2"
psamMotifs <- list(motif, motif2)
motifStack(psamMotifs)
Plot affinity logos with sequence logo stack style

Plot affinity logos with sequence logo stack style

## plot stacks with hierarchical tree
motifStack(motifs, layout="tree")
Sequence logo stack with hierarchical cluster tree

Sequence logo stack with hierarchical cluster tree

## When the number of motifs is too much to be shown in a vertical stack, 
## motifStack can draw them in a radial style.
## random sample from MotifDb
library("MotifDb")
matrix.fly <- query(MotifDb, "Dmelanogaster")
motifs2 <- as.list(matrix.fly)
## use data from FlyFactorSurvey
motifs2 <- motifs2[grepl("Dmelanogaster\\-FlyFactorSurvey\\-",
                         names(motifs2))]
## format the names
names(motifs2) <- gsub("Dmelanogaster_FlyFactorSurvey_", "",
                       gsub("_FBgn\\d+$", "",
                            gsub("[^a-zA-Z0-9]","_",
                                 gsub("(_\\d+)+$", "", names(motifs2)))))
motifs2 <- motifs2[unique(names(motifs2))]
pfms <- sample(motifs2, 30)
## creat a list of object of pfm 
motifs2 <- mapply(pfms, names(pfms), FUN=function(.ele, .name){
  new("pfm",mat=.ele, name=.name)}, SIMPLIFY = FALSE)
## trim the motifs
motifs2 <- lapply(motifs2, trimMotif, t=0.4)
## setting colors
library(RColorBrewer)
color <- brewer.pal(10, "Set3")
## plot logo stack with radial style
motifStack(motifs2, layout="radialPhylog", 
           circle=0.3, cleaves = 0.2, 
           clabel.leaves = 0.5, 
           col.bg=rep(color, each=3), col.bg.alpha=0.3, 
           col.leaves=rep(color, each=3),
           col.inner.label.circle=rep(color, each=3), 
           inner.label.circle.width=0.05,
           col.outer.label.circle=rep(color, each=3), 
           outer.label.circle.width=0.02, 
           circle.motif=1.2,
           angle=350)
Plot motifs in a radial style when the number of motifs is too much to be shown in a vertical stack

Plot motifs in a radial style when the number of motifs is too much to be shown in a vertical stack

plot a sequence logo cloud

We can also plot a sequence logo cloud for DNA motifs.

## assign groups for motifs
groups <- rep(paste("group",1:5,sep=""), each=10)
names(groups) <- names(pfms)
## assign group colors
group.col <- brewer.pal(5, "Set3")
names(group.col)<-paste("group",1:5,sep="")
## create a list of pfm objects
pfms <- mapply(names(pfms), pfms, FUN=function(.ele, .pfm){
  new("pfm",mat=.pfm, name=.ele)}
               ,SIMPLIFY = FALSE)
## use matalign to calculate the distances of motifs
hc <- clusterMotifs(pfms)
## convert the hclust to phylog object
library(ade4)
phylog <- ade4::hclust2phylog(hc)
## reorder the pfms by the order of hclust
leaves <- names(phylog$leaves)
pfms <- pfms[leaves]
## extract the motif signatures
motifSig <- motifSignature(pfms, phylog, cutoffPval=0.0001, min.freq=1)
## draw the motifs with a tag-cloud style.
motifCloud(motifSig, scale=c(6, .5), 
           layout="rectangles", 
           group.col=group.col, 
           groups=groups, 
           draw.legend=TRUE)
Sequence logo cloud with rectangle packing layout

Sequence logo cloud with rectangle packing layout

motifCircos

We can also plot it with circos style. In circos style, we can plot two group of motifs and with multiple color rings.

## 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(), 30), 
                           sample(colors(), 30), 
                           sample(colors(), 30)),
            angle=350, motifScale="logarithmic")
Grouped sequence logo with circos style layout

Grouped sequence logo with circos style layout

motifPiles

We can also plot the motifs in pile style. In pile style, we can plot two group of motifs with multiple types of annotation, for example heatmap. The col.anno parameter should be set as a named list.

## plot the logo stack with heatmap.
df <- data.frame(A=runif(n = 30), B=runif(n = 30), C=runif(n = 30), D=runif(n = 30))
map2col <- function(x, pal){
  rg <- range(x)
  pal[findInterval(x, seq(rg[1], rg[2], length.out = length(pal)+1), 
                   all.inside = TRUE)]
}
dl <- lapply(df, map2col, pal=heat.colors(10))
## alignment of the pfms, this step will make the motif logos occupy 
## more space. Users can skip this alignment to see the difference.
pfmsAligned <- DNAmotifAlignment(pfms)
## plot motifs
motifPiles(phylog=phylog, pfms=pfmsAligned, 
            col.tree=rep(color, each=5),
            col.leaves=rep(rev(color), each=5),
            col.pfms2=gpCol, 
            r.anno=rep(0.02, length(dl)), 
            col.anno=dl,
            motifScale="logarithmic",
            plotIndex=TRUE,
            groupDistance=10)
Grouped sequence logo with a heatmap

Grouped sequence logo with a heatmap

plot motifs with d3.js

Interactive plot can be generated using browseMotifs function which leverages the d3.js library. All motifs on the plot are draggable and the plot can be easily exported as a Scalable Vector Graphics (SVG) file.

browseMotifs(pfms = pfms, phylog = phylog, layout="tree", yaxis = FALSE, baseWidth=6, baseHeight = 15)

Plot the motifs in radialPhylog layout.

browseMotifs(pfms = pfms, phylog = phylog, layout="radialPhylog", yaxis = FALSE, xaxis = FALSE, baseWidth=6, baseHeight = 15)

docker container for motifStack

Docker container allows software to be packaged into containers which can be run in any platform using a virtual machine called boot2docker. To ease the installation of motifStack and its depencies, we have created a docker image containing all the components needed to run motifStack. Users can download the motifStack docker image using the following code snippet.

cd ~ ## in windows, please try cd c:\\ Users\\ username
docker pull jianhong/motifstack:latest
mkdir tmp4motifstack ## this will be the share folder for your host and container.
docker run -ti --rm -v ${PWD}/tmp4motifstack:/volume/data jianhong/motifstack:latest bash
  In motifstack:latest docker
    1  cd /volume/data
    2  git clone https://github.com/jianhong/motifStack.documentation.git
    3  cd motifStack.documentation/
    4  cp /usr/bin/matalign app/matalign-v4a
    5  cp /usr/bin/phylip/neighbor app/neighbor.app/Contents/MacOS/neighbor
    6  R cmd -e "rmarkdown::render('suppFigure2.Rmd')"
    7  R cmd -e "rmarkdown::render('suppFigure6.Rmd')"

You will see the test.pdf file in the folder of tmp4motifstack.

plot motifs with ggplot2

motifs could be plotted by geom_motif function.

pcm <- read.table(file.path(find.package("motifStack"), 
                            "extdata", "bin_SOLEXA.pcm"))
pcm <- pcm[,3:ncol(pcm)]
rownames(pcm) <- c("A","C","G","T")
markerRect <- new("marker", type="rect", start=6, stop=7, gp=gpar(lty=2, fill=NA, col="orange"))
markerLine <- new("marker", type="line", start=3, stop=5, gp=gpar(lwd=2, col="red"))
markerText <- new("marker", type="text", start=1, label="*", gp=gpar(cex=2, col="red"))
motif <- new("pcm", mat=as.matrix(pcm), name="bin_SOLEXA", 
             markers=c(markerRect, markerLine, markerText))
pfm <- pcm2pfm(motif)
df <- data.frame(xmin=c(.25, .25), ymin=c(.25, .75), xmax=c(.75, .75), ymax=c(.5, 1), 
                 fontfamily=c("serif", "mono"), fontface=c(2, 1))
df$motif <- list(pfm, pfm)

library(ggplot2)

ggplot(df, aes(xmin=xmin, ymin=ymin, xmax=xmax, ymax=ymax, motif=motif, 
               fontfamily=fontfamily, fontface=fontface)) + 
    geom_motif() + theme_bw() + ylim(0, 1) + xlim(0, 1)

df <- data.frame(x=.5, y=c(.25, .75), width=.5, height=.25, 
                 fontfamily=c("serif", "mono"), fontface=c(2, 1))
df$motif <- list(pfm, pfm)

ggplot(df, aes(x=x, y=y, width=width, height=height, motif=motif, 
               fontfamily=fontfamily, fontface=fontface)) + 
    geom_motif(use.xy=TRUE) + theme_bw() + ylim(0, 1) + xlim(0, 1)

Session Info

sessionInfo()
## R version 4.4.0 beta (2024-04-15 r86425)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    grid      stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-3   JASPAR2020_0.99.10   TFBSTools_1.42.0    
##  [4] ggplot2_3.5.1        ade4_1.7-22          MotifDb_1.46.0      
##  [7] Biostrings_2.72.0    XVector_0.44.0       GenomicRanges_1.56.0
## [10] GenomeInfoDb_1.40.0  IRanges_2.38.0       S4Vectors_0.42.0    
## [13] BiocGenerics_0.50.0  motifStack_1.48.0    knitr_1.46          
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.2.2                   bitops_1.0-7               
##   [3] rlang_1.1.3                 magrittr_2.0.3             
##   [5] matrixStats_1.3.0           compiler_4.4.0             
##   [7] RSQLite_2.3.6               png_0.1-8                  
##   [9] vctrs_0.6.5                 reshape2_1.4.4             
##  [11] stringr_1.5.1               pwalign_1.0.0              
##  [13] pkgconfig_2.0.3             crayon_1.5.2               
##  [15] fastmap_1.1.1               labeling_0.4.3             
##  [17] splitstackshape_1.4.8       caTools_1.18.2             
##  [19] utf8_1.2.4                  Rsamtools_2.20.0           
##  [21] rmarkdown_2.26              tzdb_0.4.0                 
##  [23] pracma_2.4.4                UCSC.utils_1.0.0           
##  [25] DirichletMultinomial_1.46.0 bit_4.0.5                  
##  [27] xfun_0.43                   zlibbioc_1.50.0            
##  [29] cachem_1.0.8                CNEr_1.40.0                
##  [31] jsonlite_1.8.8              blob_1.2.4                 
##  [33] highr_0.10                  DelayedArray_0.30.0        
##  [35] BiocParallel_1.38.0         jpeg_0.1-10                
##  [37] parallel_4.4.0              R6_2.5.1                   
##  [39] bslib_0.7.0                 stringi_1.8.3              
##  [41] rtracklayer_1.64.0          jquerylib_0.1.4            
##  [43] Rcpp_1.0.12                 SummarizedExperiment_1.34.0
##  [45] base64enc_0.1-3             R.utils_2.12.3             
##  [47] readr_2.1.5                 Matrix_1.7-0               
##  [49] tidyselect_1.2.1            abind_1.4-5                
##  [51] yaml_2.3.8                  codetools_0.2-20           
##  [53] curl_5.2.1                  lattice_0.22-6             
##  [55] tibble_3.2.1                plyr_1.8.9                 
##  [57] withr_3.0.0                 KEGGREST_1.44.0            
##  [59] Biobase_2.64.0              evaluate_0.23              
##  [61] pillar_1.9.0                BiocManager_1.30.22        
##  [63] MatrixGenerics_1.16.0       generics_0.1.3             
##  [65] grImport2_0.3-1             RCurl_1.98-1.14            
##  [67] hms_1.1.3                   munsell_0.5.1              
##  [69] scales_1.3.0                BiocStyle_2.32.0           
##  [71] gtools_3.9.5                xtable_1.8-4               
##  [73] glue_1.7.0                  seqLogo_1.70.0             
##  [75] tools_4.4.0                 TFMPvalue_0.0.9            
##  [77] BiocIO_1.14.0               data.table_1.15.4          
##  [79] BSgenome_1.72.0             annotate_1.82.0            
##  [81] GenomicAlignments_1.40.0    XML_3.99-0.16.1            
##  [83] Cairo_1.6-2                 poweRlaw_0.80.0            
##  [85] AnnotationDbi_1.66.0        colorspace_2.1-0           
##  [87] GenomeInfoDbData_1.2.12     restfulr_0.0.15            
##  [89] cli_3.6.2                   fansi_1.0.6                
##  [91] S4Arrays_1.4.0              dplyr_1.1.4                
##  [93] gtable_0.3.5                R.methodsS3_1.8.2          
##  [95] sass_0.4.9                  digest_0.6.35              
##  [97] SparseArray_1.4.0           rjson_0.2.21               
##  [99] htmlwidgets_1.6.4           R.oo_1.26.0                
## [101] memoise_2.0.1               htmltools_0.5.8.1          
## [103] lifecycle_1.0.4             httr_1.4.7                 
## [105] GO.db_3.19.1                bit64_4.0.5                
## [107] MASS_7.3-60.2

Reference

1. Bembom, O. SeqLogo: Sequence logos for dna sequence alignments. R package version 1.5.4 (2006).

2. Foat, B. C., Morozov, A. V. & Bussemaker, H. J. Statistical mechanical modeling of genome-wide transcription factor occupancy data by matrixreduce. Bioinformatics 22, e141–e149 (2006).

3. Mercier, E. & Gottardo, R. MotIV: Motif identification and validation. R package version 1.10.0 (2010).

4. S, M. & PV, B. STAMP: A web tool for exploring dna-binding motif similarities. Nucleic Acids Res. 35(Web Server issue), W253–W258 (2007).