Memes: A motif analysis environment in R using tools from the MEME Suite (2024)

Identification of biopolymer motifs represents a key step in the analysis of biological sequences. The MEME Suite is a widely used toolkit for comprehensive analysis of biopolymer motifs; however, these tools are poorly integrated within popular analysis frameworks like the R/Bioconductor project, creating barriers to their use. Here we present memes, an R package that provides a seamless R interface to a selection of popular MEME Suite tools. memes provides a novel “data aware” interface to these tools, enabling rapid and complex discriminative motif analysis workflows. In addition to interfacing with popular MEME Suite tools, memes leverages existing R/Bioconductor data structures to store the multidimensional data returned by MEME Suite tools for rapid data access and manipulation. Finally, memes provides data visualization capabilities to facilitate communication of results. memes is available as a Bioconductor package at https://bioconductor.org/packages/memes, and the source code can be found at github.com/snystrom/memes.

Biologically active molecules such as DNA, RNA, and proteins are polymers composed of repeated subunits. For example, nucleotides are the subunits of DNA and RNA, and amino acids are the subunits of proteins. Functional properties of biopolymers are determined by short, recurring stretches of subunits known as motifs. Motifs can serve as binding sites between molecules, they can influence the structure of molecules, and they can contribute to enzymatic activities. For these reasons, motif analysis has become a key step in determining the function of biopolymers and in elucidating their roles in biological networks. The MEME Suite is a widely used compilation of tools used for identifying and analyzing motifs found in biological sequences. Here, we describe a new piece of software named “memes” that connects MEME Suite tools to R, the statistical analysis environment. By providing an interface between the MEME Suite and R, memes allows for improved motif analysis workflows and easy access to a wide range of existing data visualization tools, further expanding the utility of MEME Suite tools.

Table 1. A list of MEME Suite tools that are currently included in memes.

https://doi.org/10.1371/journal.pcbi.1008991.t001

Design & implementation

Data visualization

The MEME Suite provides a small set of data visualizations that have limited customizability. memes leverages the advantages of the R graphics environment to provide a wide range of data visualization options that are highly customizable. We describe two scenarios below.

The TomTom tool allows users to compare unknown motifs to a set of known motifs to identify the best match. Visual inspection of motif comparison data is a key step in assessing the quality and accuracy of a match [12]. The view_tomtom_hits() function allows users to compare query motifs with the list of potential matches as assigned by TomTom, similar to the commandline behavior. The force_best_match() function allows users to reassign the TomTom best match to a lower-ranked (but still statistically significant) match in order to highlight motifs with greater biological relevance (e.g. to skip over a transcription factor that is not expressed in the experimental sample, or when two motifs are matched equally well).

The AME tool searches for enrichment of motifs within a set of experimental sequences relative to a control set [15]. The meaningful result from this tool is the statistical parameter (for example, a p-value) associated with the significance of motif enrichment. However, AME does not provide a mechanism for visualizing these results.

The plot_ame_heatmap() function in memes returns a ggplot2 formatted heatmap of statistical significance of motif enrichment. If AME is used to examine motif content of multiple groups of sequences, the plot_ame_heatmap() function can also return a plot comparing motif significance within multiple groups. Several options exist to customize the heatmap values in order to capture different aspects of the output data. The ame_compare_heatmap_methods() function enables users to compare the distribution of values between samples in order to select a threshold that accurately captures the dynamic range of their results.

Containerized analysis maximizes availability and facilitates reproducibility

memes relies on a locally installed version of the MEME Suite which is accessed by the user’s R process. Although R is available on Windows, Mac, and Linux operating systems, the MEME Suite is incompatible with Windows, limiting its adoption by Windows users. Additional barriers also exist to installing a local copy of the MEME Suite on compatible systems, for example, the MEME Suite relies on several system-level dependencies whose installation can be difficult for novice users. Finally, some tools in the MEME Suite use python to generate shuffled control sequences for analysis, which presents a reproducibility issue as the random number generation algorithm changed between python2.7 and python3. The MEME Suite will build and install on both python2.7 and python3 systems, therefore without careful consideration of this behavior, the same code run on two systems may not produce identical results, even if using the same major version of the MEME Suite. In order to increase access to the MEME Suite on unsupported operating systems, and to facilitate reproducible motif analysis using memes, we have also developed a docker container with a preinstalled version of the MEME Suite along with an R/Bioconductor analysis environment including the most recent version of memes and its dependencies. As new container versions are released, they are version-tagged and stored to ensure reproducibility while allowing updates to the container.

De novo motif analysis

The DREME tool can be used to discover short, novel motifs present in a set of input sequences relative to a control set. The runDreme() command is the memes interface to the DREME tool. runDreme() syntax enables users to produce complex discriminative analysis designs using intuitive syntax. Examples of possible designs and their syntax are compared below.

runDreme(all_sequences, )

runDreme(sequences_by_response, )

runDreme(sequences_by_response$Increasing, sequences_by_response$Decreasing)

runDreme(sequences_by_response, )

runDreme(sequences_by_response, c(, ))

In order to discover motifs associated with dynamic changes to chromatin accessibility, we use open chromatin sites that do not change in accessibility (i.e. “static sites”) as the control set to discover motifs enriched in either increasing or decreasing sites.

dreme_vs_static <- runDreme(sequences_by_response, control = "Static")

These data are readily visualized using the universalmotif::view_motifs() function. Visualization of the results reveals two distinct motifs associated with decreasing and increasing sites (Fig 1). In this analysis, the de-novo discovered motifs within each category appear visually similar to each other; however, the MEME Suite does not provide a mechanism to compare groups of motifs based on all pairwise similarity metrics. We utilized the universalmotif::compare_motifs() function to compute Pearson correlation coefficients for each set of motifs to quantitatively assess motif similarity (Fig 1).

Transcription factors often bind with sequence specificity to regulate activity of nearby genes. Therefore, comparison of de-novo motifs with known DNA binding motifs for transcription factors can be a key step in identifying transcriptional regulators. The TomTom tool is used to compare unknown motifs against a list of known motifs to identify matches, memes provides the runTomTom() function as an interface to this tool. By passing the results from runDreme() into runTomTom() and searching within a database of Drosophila transcription factor motifs, we can identify candidate transcription factors that may bind the motifs associated with increasing and decreasing chromatin accessibility.

dreme_vs_static_tomtom runTomTom(dreme_vs_static, , )

Using this approach, the results can be visualized using the view_tomtom_hits() function to visually inspect the matches assigned by TomTom, providing a simple way to assess the quality of the match. A representative plot of these data is shown in Fig 2, while the full results can be viewed in Table 2. This analysis revealed that the motifs associated with decreasing sites are highly similar to the E93 motif. Conversely, motifs found in increasing sites match the transcription factors Br, L(3)neo38, Gl, and Lola (Table 2). These data support the hypothesis that E93 binding to its motif may play a larger role in chromatin closing activity, while binding of the transcription factor br is associated with chromatin opening.

Motif enrichment testing using runAme

Discovery and matching of de-novo motifs is only one way to find candidate transcription factors within target sites. Indeed, in many instances, requiring that a motif is recovered de-novo is not ideal, as these approaches are less sensitive than targeted searches. Another approach, implemented by the AME tool, is to search for known motif instances in target sequences and test for their overrepresentation relative to a background set of sequences [15]. The runAme() function is the memes interface to the AME tool. It accepts a set of sequences as input and control sets, and it will perform enrichment testing using a provided motif database for occurrences of each provided motif.

A major limitation of this approach is that transcription factors containing similar families of DNA binding domain often possess highly similar motifs, making it difficult to identify the “true” binding factor associated with an overrepresented motif. Additionally, when searching for matches against a motif database, AME must account for multiple testing, therefore using a larger than necessary motif database can produce a large multiple testing penalty, limiting sensitivity of detection. One way to overcome these limitations is to limit the transcription factor motif database to include only motifs for transcription factors expressed in the sample of interest. Accounting for transcription factor expression during motif analysis has been demonstrated to increase the probability of identifying biologically relevant transcription factor candidates [14,16].

The universalmotif_df structure can be used to integrate expression data with a motif database to remove entries for transcription factors that are not expressed. To do so, we import a Drosophila transcription factor motif database generated by the Fly Factor Survey and convert to universalmotif_df format [17]. In this database, the altname column stores the gene symbol.

fly_factor <- universalmotif::read_meme() %>%

setNames(., purrr::map_chr(., ~{.x[]})) %>%

to_df()

## # A tibble: 3 x 16

##name altname family organism consensusalphabet strand icscore nsites

##<chr><chr><chr><chr><chr><chr><chr><dbl> <int>

## 1 ab_SANG… ab<NA><NA>BWNRCCAGGWMCN… DNA+-15.020

## 2 ab_SOLE… ab<NA><NA>NNNNHNRCCAGGW… DNA+-14.6446

## 3 Abd-A_F… abd-A<NA><NA>KNMATWAWDNA+-7.37 37

## # … with 7 more variables: bkgsites <int>, pval <dbl>, qval <dbl>, eval <dbl>,

## #type <chr>, bkg <named list>, motif <I<named list>>

Next, we import a pre-filtered list of genes expressed in a timecourse of Drosophila wing development.

wing_expressed_genes <- read.csv("data/wing_expressed_genes.csv")

Finally, we subset the motif database to only expressed genes using dplyr data.frame subsetting syntax (note that base R subsetting functions operate equally well on these data structures), then convert the data.frame back into universalmotif format using to_list(). This filtering step removes 43% of entries from the original database, greatly reducing the multiple-testing correction.

fly_factor_expressed <- fly_factor %>%

dplyr::filter(altname %in% wing_expressed_genes$symbol) %>%

to_list()

runAme() syntax is identical to runDreme() in that discriminative designs can be constructed by calling list entries by name. Data can be visualized using plot_ame_heatmap(), revealing that the E93 motif is strongly enriched at Decreasing sites, while l(3)neo38, lola, and br motifs are enriched in Increasing sites, supporting the de-novo discovery results (Fig 3). Additionally, AME detects several other transcription factor motifs that distinguish decreasing and increasing sites, providing additional clues to potential factors that also bind with E93 to affect chromatin accessibility.

ame_vs_static <- runAme(sequences_by_response, control = "Static", database = fly_factor_expressed)