-
Notifications
You must be signed in to change notification settings - Fork 5
/
README.Rmd
351 lines (273 loc) · 15 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
fig.align = "center"
#out.width = "100%"
)
```
# memes
<!-- badges: start -->
[![Lifecycle: stable](https://img.shields.io/badge/lifecycle-stable-brightgreen.svg)](https://lifecycle.r-lib.org/articles/stages.html#stable)
[![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
[![Codecov test coverage](https://codecov.io/gh/snystrom/memes/branch/master/graph/badge.svg)](https://codecov.io/gh/snystrom/memes?branch=master)
![R-CMD-check-bioc](https://github.com/snystrom/memes/workflows/R-CMD-check-bioc/badge.svg)
![Bioconductor Build Status](https://bioconductor.org/shields/build/devel/bioc/memes.svg)
![Bioconductor Lifetime](https://bioconductor.org/shields/years-in-bioc/memes.svg)
<!-- badges: end -->
An R interface to the [MEME Suite](http://meme-suite.org/) family of tools,
which provides several utilities for performing motif analysis on DNA, RNA, and
protein sequences. memes works by detecting a local install of the MEME suite,
running the commands, then importing the results directly into R.
## Installation
### Bioconductor
memes is currently available on the Bioconductor `devel` branch:
```{r, eval=F}
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# The following initializes usage of Bioc devel
BiocManager::install(version='devel')
BiocManager::install("memes")
```
### Development Version (Github)
You can install the development version of memes from [GitHub](https://github.com/snystrom/memes) with:
```{r, eval=F}
if (!requireNamespace("remotes", quietly=TRUE))
install.packages("remotes")
remotes::install_github("snystrom/memes")
# To temporarily bypass the R version 4.1 requirement, you can pull from the following branch:
remotes::install_github("snystrom/memes", ref = "no-r-4")
```
### Docker Container
```{shell}
# Get development version from dockerhub
docker pull snystrom/memes_docker:devel
# the -v flag is used to mount an analysis directory,
# it can be excluded for demo purposes
docker run -e PASSWORD=<password> -p 8787:8787 -v <path>/<to>/<project>:/mnt/<project> snystrom/memes_docker:devel
```
## Detecting the MEME Suite
memes relies on a local install of the [MEME Suite](http://meme-suite.org/).
For installation instructions for the MEME suite, see the [MEME Suite
Installation Guide](http://meme-suite.org/doc/install.html?man_type=web).
memes needs to know the location of the `meme/bin/` directory on your local machine.
You can tell memes the location of your MEME suite install in 4 ways. memes
will always prefer the more specific definition if it is a valid path. Here they
are ranked from most- to least-specific:
1. Manually passing the install path to the `meme_path` argument of all memes functions
2. Setting the path using `options(meme_bin = "/path/to/meme/bin/")` inside your R script
3. Setting `MEME_BIN=/path/to/meme/bin/` in your `.Renviron` file
4. memes will try the default MEME install location `~/meme/bin/`
If memes fails to detect your install at the specified location, it will fall
back to the next option.
To verify memes can detect your MEME install, use `check_meme_install()` which
uses the search herirarchy above to find a valid MEME install. It will report
whether any tools are missing, and print the path to MEME that it sees. This can
be useful for troubleshooting issues with your install.
```{r check_install_works}
library(memes)
# Verify that memes detects your meme install
# (returns all green checks if so)
check_meme_install()
```
```{r check_install_fails}
# You can manually input a path to meme_path
# If no meme/bin is detected, will return a red X
check_meme_install(meme_path = 'bad/path')
```
## The Core Tools
| Function Name | Use | Sequence Input | Motif Input | Output |
|:-------------:|:----------------:|:--------------:|:-----------:|:-------------------------------------------------------|
| `runStreme()` | Motif Discovery (short motifs) | Yes | No | `universalmotif_df` |
| `runDreme()` | Motif Discovery (short motifs) | Yes | No | `universalmotif_df` |
| `runAme()` | Motif Enrichment | Yes | Yes | data.frame (optional: `sequences` column) |
| `runFimo()` | Motif Scanning | Yes | Yes | GRanges of motif positions |
| `runTomTom()` | Motif Comparison | No | Yes | `universalmotif_df` w/ `best_match_motif` and `tomtom` columns* |
| `runMeme()` | Motif Discovery (long motifs) | Yes | No | `universalmotif_df` |
\* **Note:** if `runTomTom()` is run using a `universalmotif_df`
the results will be joined with the `universalmotif_df` results as extra
columns. This allows easy comparison of *de-novo* discovered motifs with their
matches.
**Sequence Inputs** can be any of:
1. Path to a .fasta formatted file
2. `Biostrings::XStringSet` (can be generated from GRanges using `get_sequence()` helper function)
3. A named list of `Biostrings::XStringSet` objects (generated by `get_sequence()`)
**Motif Inputs** can be any of:
1. A path to a .meme formatted file of motifs to scan against
2. A `universalmotif` object, or list of `universalmotif` objects
3. A `runDreme()` results object (this allows the results of `runDreme()` to pass directly to `runTomTom()`)
4. A combination of all of the above passed as a `list()` (e.g. `list("path/to/database.meme", "dreme_results" = dreme_res)`)
**Output Types**:
`runDreme()`, `runStreme()`, `runMeme()` and `runTomTom()` return
`universalmotif_df` objects which are data.frames with special columns. The
`motif` column contains a `universalmotif` object, with 1 entry per row. The
remaining columns describe the properties of each returned motif. The following
column names are special in that their values are used when running
`update_motifs()` and `to_list()` to alter the properties of the motifs stored
in the `motif` column. Be careful about changing these values as these changes
will propagate to the `motif` column when calling `update_motifs()` or
`to_list()`.
- name
- altname
- family
- organism
- strand
- nsites
- bkgsites
- pval
- qval
- eval
memes is built around the [universalmotif package](https://www.bioconductor.org/packages/release/bioc/html/universalmotif.html)
which provides a framework for manipulating motifs in R. `universalmotif_df`
objects can interconvert between data.frame and `universalmotif` list format
using the `to_df()` and `to_list()` functions, respectively. This allows use of
`memes` results with all other Bioconductor motif packages, as `universalmotif`
objects can convert to any other motif type using `convert_motifs()`.
`runTomTom()` returns a special column: `tomtom` which is a `data.frame` of all
match data for each input motif. This can be expanded out using
`tidyr::unnest(tomtom_results, "tomtom")`, and renested with `nest_tomtom()`.
The `best_match_` prefixed columns returned by `runTomTom()` indicate values for
the motif which was the best match to the input motif.
## Quick Examples
### Motif Discovery with DREME
```{r}
suppressPackageStartupMessages(library(magrittr))
suppressPackageStartupMessages(library(GenomicRanges))
# Example transcription factor peaks as GRanges
data("example_peaks", package = "memes")
# Genome object
dm.genome <- BSgenome.Dmelanogaster.UCSC.dm6::BSgenome.Dmelanogaster.UCSC.dm6
```
The `get_sequence` function takes a `GRanges` or `GRangesList` as input and
returns the sequences as a `BioStrings::XStringSet`, or list of `XStringSet`
objects, respectively. `get_sequence` will name each fasta entry by the genomic
coordinates each sequence is from.
```{r}
# Generate sequences from 200bp about the center of my peaks of interest
sequences <- example_peaks %>%
resize(200, "center") %>%
get_sequence(dm.genome)
```
`runDreme()` accepts XStringSet or a path to a fasta file as input. You can use
other sequences or shuffled input sequences as the control dataset.
```{r}
# runDreme accepts all arguments that the commandline version of dreme accepts
# here I set e = 50 to detect motifs in the limited example peak list
# In a real analysis, e should typically be < 1
dreme_results <- runDreme(sequences, control = "shuffle", e = 50)
```
memes is built around the
[universalmotif](https://www.bioconductor.org/packages/release/bioc/html/universalmotif.html)
package. The results are returned in `universalmotif_df` format, which is an R data.frame that can seamlessly interconvert between data.frame and `universalmotif` format using `to_list()` to convert to `universalmotif` list format, and `to_df()` to convert back to data.frame format. Using `to_list()` allows using `memes` results with all `universalmotif` functions like so:
```{r}
library(universalmotif)
dreme_results %>%
to_list() %>%
view_motifs()
```
### Matching motifs using TOMTOM
Discovered motifs can be matched to known TF motifs using `runTomTom()`, which can accept as input a path to a .meme formatted file, a `universalmotif` list, or the results of `runDreme()`.
TomTom uses a database of known motifs which can be passed to the `database` parameter as a path to a .meme format file, or a `universalmotif` object.
Optionally, you can set the environment variable `MEME_DB` in `.Renviron` to a file on disk, or
the `meme_db` value in `options` to a valid .meme format file and memes will
use that file as the database. memes will always prefer user input to the
function call over a global variable setting.
```{r}
options(meme_db = system.file("extdata/flyFactorSurvey_cleaned.meme", package = "memes"))
m <- create_motif("CMATTACN", altname = "testMotif")
tomtom_results <- runTomTom(m)
```
```{r}
tomtom_results
```
### Using runDreme results as TOMTOM input
`runTomTom()` will add its results as columns to a `runDreme()` results data.frame.
```{r}
full_results <- dreme_results %>%
runTomTom()
```
### Motif Enrichment using AME
AME is used to test for enrichment of known motifs in target sequences. `runAme()`
will use the `MEME_DB` entry in `.Renviron` or `options(meme_db =
"path/to/database.meme")` as the motif database. Alternately, it will accept all
valid inputs similar to `runTomTom()`.
```{r}
# here I set the evalue_report_threshold = 30 to detect motifs in the limited example sequences
# In a real analysis, evalue_report_threshold should be carefully selected
ame_results <- runAme(sequences, control = "shuffle", evalue_report_threshold = 30)
ame_results
```
## Visualizing Results
`view_tomtom_hits` allows comparing the input motifs to the top hits from
TomTom. Manual inspection of these matches is important, as sometimes the top
match is not always the correct assignment. Altering `top_n` allows you to show
additional matches in descending order of their rank.
```{r}
full_results %>%
view_tomtom_hits(top_n = 1)
```
It can be useful to view the results from `runAme()` as a heatmap.
`plot_ame_heatmap()` can create complex visualizations for analysis of enrichment
between different region types (see vignettes for details). Here is a simple
example heatmap.
```{r, fig.height=3, fig.width=5}
ame_results %>%
plot_ame_heatmap()
```
# Scanning for motif occurances using FIMO
The FIMO tool is used to identify matches to known motifs. `runFimo` will return
these hits as a `GRanges` object containing the genomic coordinates of the motif
match.
```{r, fig.height=4, fig.width=3}
# Query MotifDb for a motif
e93_motif <- MotifDb::query(MotifDb::MotifDb, "Eip93F") %>%
universalmotif::convert_motifs()
# Scan for the E93 motif within given sequences
fimo_results <- runFimo(sequences, e93_motif, thresh = 1e-3)
# Visualize the sequences matching the E93 motif
plot_sequence_heatmap(fimo_results$matched_sequence)
```
## Importing Data from previous runs
memes also supports importing results generated using the MEME suite outside of
R (for example, running jobs on [meme-suite.org](meme-suite.org), or running on
the commandline). This enables use of preexisting MEME suite results with
downstream memes functions.
| MEME Tool | Function Name | File Type |
|:---------:|:-------------------:|:----------------:|
| Streme | `importStremeXML()` | streme.xml |
| Dreme | `importDremeXML()` | dreme.xml |
| TomTom | `importTomTomXML()` | tomtom.xml |
| AME | `importAme()` | ame.tsv* |
| FIMO | `importFimo()` | fimo.tsv |
| Meme | `importMeme()` | meme.txt |
\* `importAME()` can also use the "sequences.tsv" output when AME used `method = "fisher"`, this is optional.
# FAQs
### How do I use memes/MEME on Windows?
The MEME Suite does not currently support Windows, although it can be
installed under [Cygwin](https://www.cygwin.com/) or the [Windows Linux Subsytem](https://docs.microsoft.com/en-us/windows/wsl/install-win10) (WSL).
Please note that if MEME is installed on Cygwin or WSL, you must also run R
inside Cygwin or WSL to use memes.
An alternative solution is to use [Docker](https://www.docker.com/get-started)
to run a virtual environment with the MEME Suite installed. We provide a [memes docker container](https://github.com/snystrom/memes_docker)
that ships with the MEME Suite, R studio, and all `memes` dependencies
pre-installed.
# Citation
memes is a wrapper for a select few tools from the MEME Suite, which were
developed by another group. In addition to citing memes, please cite the MEME
Suite tools corresponding to the tools you use.
If you use `runDreme()` in your analysis, please cite:
Timothy L. Bailey, "DREME: Motif discovery in transcription factor ChIP-seq data", Bioinformatics, 27(12):1653-1659, 2011. [full text](https://academic.oup.com/bioinformatics/article/27/12/1653/257754)
If you use `runTomTom()` in your analysis, please cite:
Shobhit Gupta, JA Stamatoyannopolous, Timothy Bailey and William Stafford Noble, "Quantifying similarity between motifs", Genome Biology, 8(2):R24, 2007. [full text](http://genomebiology.com/2007/8/2/R24)
If you use `runAme()` in your analysis, please cite:
Robert McLeay and Timothy L. Bailey, "Motif Enrichment Analysis: A unified framework and method evaluation", BMC Bioinformatics, 11:165, 2010, doi:10.1186/1471-2105-11-165. [full text](http://www.biomedcentral.com/1471-2105/11/165)
If you use `runFimo()` in your analysis, please cite:
Charles E. Grant, Timothy L. Bailey, and William Stafford Noble, "FIMO: Scanning for occurrences of a given motif", Bioinformatics, 27(7):1017-1018, 2011. [full text](http://bioinformatics.oxfordjournals.org/content/early/2011/02/16/bioinformatics.btr064.full)
## Licensing Restrictions
The MEME Suite is free for non-profit use, but for-profit users should purchase a
license. See the [MEME Suite Copyright Page](http://meme-suite.org/doc/copyright.html) for details.