The motifmatchr package is designed for analyzing many sequences and many motifs to find which sequences contain which motifs. It uses the MOODS C++ library (developedby Pasi Rastas, Janne Korhonen, and Petri Martinmaki) internally for motif matching. The primary method of motifmatchr is matchMotifs, which takes in motif PWMs/PFMs and genomic ranges or sequences and returns either which ranges/sequences match which motifs or the positions of the matches.
Compared with alternative motif matching functions available in Bioconductor (e.g. matchPWM in Biostrings or searchSeq in TFBSTools), motifmatchr is designed specifically for the use case of determining whether many different sequences/ranges contain many different motifs. For example, when analyzing ChIP-seq or ATAC-seq data one might want to find what motifs in a collection of motifs like the JASPAR database are found in what peaks.
Example use case of motifmatchr with a set of peaks and a few motifs. For additional options for inputs & outputs, see remainder of vignette.
library(motifmatchr)
library(GenomicRanges)
library(SummarizedExperiment)
library(BSgenome)
# load some example motifs
data(example_motifs, package = "motifmatchr")
# Make a set of peaks
peaks <- GRanges(seqnames = c("chr1","chr2","chr2"),
ranges = IRanges(start = c(76585873,42772928,100183786),
width = 500))
# Get motif matches for example motifs in peaks
motif_ix <- matchMotifs(example_motifs, peaks, genome = "hg19")
motifMatches(motif_ix) # Extract matches matrix from result## 3 x 3 sparse Matrix of class "lgCMatrix"
## MA0599.1_KLF5 MA0107.1_RELA MA0137.3_STAT1
## [1,] | | .
## [2,] | . |
## [3,] | . |# Get motif positions within peaks for example motifs in peaks
motif_pos <- matchMotifs(example_motifs, peaks, genome = "hg19",
out = "positions") This method has two mandatory arguments:
Position weight matrices or position frequency matrices, stored in the PWMatrix, PFMatrix, PWMatrixList, or PFMatrixList objects from the TFBSTools package
Either a set of genomic ranges (GenomicRanges or RangedSummarizedExperiment object) or a set of sequences (either DNAStringSet, DNAString, or simple character vector)
If the second argument is a set of genomic ranges, a genome sequence is also required. By default BSgenome.Hsapiens.UCSC.hg19 is used — you will have to have installed BSgenome.Hsapiens.UCSC.hg19. If using another genome build, either the appropraiate BSgenome object for your species or a DNAStringSet or FaFile object for your species should be passed to the genome argument.
# using peaks
motif_ix_peaks <- matchMotifs(example_motifs, peaks, genome = "hg19")
# using SummarizedExperiment
example_SummarizedExperiment <-
SummarizedExperiment(assays = list(counts = matrix(1,
ncol = 4,
nrow = 3)),
rowRanges = peaks)
motif_ix_SummarizedExperiment <- matchMotifs(example_motifs,
example_SummarizedExperiment,
genome = "hg19")
all.equal(motifMatches(motif_ix_peaks),
motifMatches(motif_ix_SummarizedExperiment))## [1] TRUE# using BSgenomeViews
example_BSgenomeViews <- BSgenomeViews(BSgenome.Hsapiens.UCSC.hg19, peaks)
motif_ix_BSgenomeViews <- matchMotifs(example_motifs, example_BSgenomeViews) ## Warning in matrix(log2(bg(pwm)) - log2(bg_freqs), nrow = 4, ncol =
## length(pwm), : data length [12] is not a sub-multiple or multiple of the
## number of columns [10]
## Warning in matrix(log2(bg(pwm)) - log2(bg_freqs), nrow = 4, ncol =
## length(pwm), : data length [12] is not a sub-multiple or multiple of the
## number of columns [10]## Warning in matrix(log2(bg(pwm)) - log2(bg_freqs), nrow = 4, ncol =
## length(pwm), : data length [12] is not a sub-multiple or multiple of the
## number of columns [11]all.equal(motifMatches(motif_ix_peaks), motifMatches(motif_ix_BSgenomeViews))## [1] "Attributes: < Component \"i\": Numeric: lengths (6, 9) differ >"
## [2] "Attributes: < Component \"p\": Mean relative difference: 0.5 >"
## [3] "Attributes: < Component \"x\": Lengths (6, 9) differ (comparison on first 6 components) >"# using DNAStringSet
library(Biostrings)
library(BSgenome.Hsapiens.UCSC.hg19)
example_DNAStringSet <- getSeq(BSgenome.Hsapiens.UCSC.hg19, peaks)
motif_ix_DNAStringSet <- matchMotifs(example_motifs, example_DNAStringSet)
all.equal(motifMatches(motif_ix_peaks), motifMatches(motif_ix_DNAStringSet))## [1] TRUE# using character vector
example_character <- as.character(example_DNAStringSet)
motif_ix_character <- matchMotifs(example_motifs, example_character)
all.equal(motifMatches(motif_ix_peaks), motifMatches(motif_ix_character))## [1] TRUEIn determining motif matches, background nucleotide frequencies are used. By default the background frequencies are the nucleotide frequencies within the input sequence – to use alternate frequencies, supply the bg arument to match_pwms. If the input sequences are fairly short (as in our vignette example!), it probably makes sense to use other input frequencies. Here we show how to use even frequencies:
motif_ix <- matchMotifs(example_motifs, peaks, genome = "hg19", bg = "even") We can also choose to use the frequency from the genome. In this case:
motif_ix <- matchMotifs(example_motifs, peaks, genome = "hg19", bg = "genome") A genome must be specified if using bg = "genome"!
To specify frequencies per base pair, the order should be “A”,“C”,“G”, then “T”, or those nucleotides should be used as names in the vector.
motif_ix <- matchMotifs(example_motifs, peaks, genome = "hg19",
bg = c("A" = 0.2,"C" = 0.3, "G" = 0.3, "T" = 0.2)) PWMatrix objects have associated background frequencies that can be accessed using the bg function. If the supplied background frequencies to match_pwms do not match the frequencies in the input PWM, then the PWM is adjusted to refect the supplied background frequencies. The calculated score is based on this adjusted PWM and not the direct input PWM. To ensure the score is computed using the exact PWM input, simply make sure the background frequencies passed to matchMotifs match hose used for the input PWM and stored in bg slot of PWMatrix object.
motifmatchr expects input PWMs to use either natural logarithm or log 2. If the input is a PFM, the TFBSTools toPWM is used for making the PWM with the default psueodcounts of 0.8 and base 2 logarithm. For more control of the pseudocounts, simply use the toPWM function to convert your PFMs prior to calling matchMotifs.
library(TFBSTools)
example_pwms <- do.call(PWMatrixList,lapply(example_motifs, toPWM,
pseudocounts = 0.5))The default p-value cutoff is 0.00005. No adjustment for multiple comparisons is made in this p-value cutoff. This p-value threshold is used to determine a score threshold. # Outputs
The matchMotifs can return three possible outputs, depending on the out argument:
(Default, with out = matches) Boolean matrix indicating which ranges/sequences contain which motifs, stored as motifMatches in assays slot of SummarizedExperiment object. The motifMatches methods can be used to extract the boolean matrix. If either the subject argument is a GenomicRanges or RangedSummarizedExperiment object OR a ranges argument is provided, then a RangedSummarizedExeriment is returned rather than a SummarizedExperiment.
(out = scores) Same as (1) plus two additional assays – a matrix with the score of the high motif score within each range/sequence (score only reported if match present) and a matrix with the number of motif matches. The motifScores and motifCounts methods can be used to access these components.
(out = positions) A GenomicRangesList with the ranges of all matches within the input ranges/sequences. Note that if the subject argument is a character vector, DNAStringSet, or DNAString and a ranges argument corresponding to the ranges represented by the sequences is NOT provided, then a list of IRangesList objects will be returned instead with the relative positions withing the sequences.
motif_ix <- matchMotifs(example_motifs, peaks, genome = "hg19")
print(motif_ix)## class: RangedSummarizedExperiment
## dim: 3 3
## metadata(0):
## assays(1): motifMatches
## rownames: NULL
## rowData names(0):
## colnames(3): MA0599.1_KLF5 MA0107.1_RELA MA0137.3_STAT1
## colData names(1): namehead(motifMatches(motif_ix))## 3 x 3 sparse Matrix of class "lgCMatrix"
## MA0599.1_KLF5 MA0107.1_RELA MA0137.3_STAT1
## [1,] | | .
## [2,] | . |
## [3,] | . |motif_ix_scores <- matchMotifs(example_motifs, peaks,
genome = "hg19", out = "scores")
print(motif_ix_scores)## class: RangedSummarizedExperiment
## dim: 3 3
## metadata(0):
## assays(3): motifScores motifMatches motifCounts
## rownames: NULL
## rowData names(0):
## colnames(3): MA0599.1_KLF5 MA0107.1_RELA MA0137.3_STAT1
## colData names(1): namehead(motifMatches(motif_ix_scores))## 3 x 3 sparse Matrix of class "lgCMatrix"
## MA0599.1_KLF5 MA0107.1_RELA MA0137.3_STAT1
## [1,] | | .
## [2,] | . |
## [3,] | . |head(motifScores(motif_ix_scores))## 3 x 3 sparse Matrix of class "dgCMatrix"
## MA0599.1_KLF5 MA0107.1_RELA MA0137.3_STAT1
## [1,] 7.975145 13.55634 .
## [2,] 7.085762 . 16.29401
## [3,] 7.739609 . 14.81298head(motifCounts(motif_ix_scores))## 3 x 3 sparse Matrix of class "dgCMatrix"
## MA0599.1_KLF5 MA0107.1_RELA MA0137.3_STAT1
## [1,] 1 1 .
## [2,] 1 . 1
## [3,] 1 . 2motif_pos <- matchMotifs(example_motifs, peaks, genome = "hg19",
out = "positions")
print(motif_pos)## GRangesList object of length 3:
## $MA0599.1_KLF5
## GRanges object with 3 ranges and 1 metadata column:
## seqnames ranges strand | score
## <Rle> <IRanges> <Rle> | <numeric>
## [1] chr1 [ 76585928, 76585937] - | 7.9751454917072
## [2] chr2 [ 42773248, 42773257] - | 7.08576191714116
## [3] chr2 [100184023, 100184032] - | 7.73960919591018
##
## $MA0107.1_RELA
## GRanges object with 1 range and 1 metadata column:
## seqnames ranges strand | score
## [1] chr1 [76585933, 76585942] + | 13.5563372551788
##
## $MA0137.3_STAT1
## GRanges object with 3 ranges and 1 metadata column:
## seqnames ranges strand | score
## [1] chr2 [ 42773087, 42773097] - | 16.2940093372611
## [2] chr2 [100183852, 100183862] + | 14.8129754070887
## [3] chr2 [100183852, 100183862] - | 13.953566440506
##
## -------
## seqinfo: 2 sequences from an unspecified genome; no seqlengthsSys.Date()## [1] "2017-05-18"sessionInfo()## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.1 LTS
##
## 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
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] motifmatchr_0.99.7 BSgenome.Hsapiens.UCSC.hg19_1.4.0
## [3] BSgenome_1.42.0 rtracklayer_1.34.2
## [5] SummarizedExperiment_1.4.0 Biobase_2.34.0
## [7] GenomicRanges_1.26.4 GenomeInfoDb_1.10.3
## [9] TFBSTools_1.12.2 Biostrings_2.42.1
## [11] XVector_0.14.1 IRanges_2.8.2
## [13] S4Vectors_0.12.2 BiocGenerics_0.20.0
##
## loaded via a namespace (and not attached):
## [1] httr_1.2.1 pkgload_0.0.0.9000
## [3] VGAM_1.0-3 splines_3.3.1
## [5] R.utils_2.5.0 gtools_3.5.0
## [7] assertthat_0.2.0 Rsamtools_1.26.2
## [9] yaml_2.1.14 DirichletMultinomial_1.16.0
## [11] highlight_0.4.7 RSQLite_1.1-2
## [13] backports_1.0.5 lattice_0.20-34
## [15] digest_0.6.12.1 colorspace_1.3-2
## [17] htmltools_0.3.6 Matrix_1.2-7.1
## [19] R.oo_1.21.0 plyr_1.8.4
## [21] XML_3.98-1.7 devtools_1.12.0.9000
## [23] zlibbioc_1.20.0 purrr_0.2.2
## [25] xtable_1.8-2 GO.db_3.4.0
## [27] scales_0.4.1 whisker_0.3-2
## [29] BiocParallel_1.8.2 git2r_0.18.0
## [31] tibble_1.3.1 annotate_1.52.1
## [33] KEGGREST_1.14.1 ggplot2_2.2.1
## [35] withr_1.0.2 TFMPvalue_0.0.6
## [37] lazyeval_0.2.0 crayon_1.3.2
## [39] magrittr_1.5 memoise_1.1.0
## [41] poweRlaw_0.70.0 evaluate_0.10
## [43] R.methodsS3_1.7.1 RcppArmadillo_0.7.700.0.0
## [45] xml2_1.1.1 pkgbuild_0.0.0.9000
## [47] CNEr_1.10.2 BiocInstaller_1.24.0
## [49] tools_3.3.1 hms_0.3
## [51] stringr_1.2.0 munsell_0.4.3
## [53] callr_1.0.0.9000 AnnotationDbi_1.36.2
## [55] pkgdown_0.1.0.9000 caTools_1.17.1
## [57] rlang_0.1 grid_3.3.1
## [59] RCurl_1.95-4.8 bitops_1.0-6
## [61] rmarkdown_1.3 testthat_1.0.2
## [63] gtable_0.2.0 DBI_0.6-1
## [65] roxygen2_6.0.1 reshape2_1.4.2
## [67] R6_2.2.1 GenomicAlignments_1.10.1
## [69] knitr_1.15.1 seqLogo_1.40.0
## [71] commonmark_1.1 rprojroot_1.2
## [73] desc_1.1.0 readr_1.1.1
## [75] stringi_1.1.5 Rcpp_0.12.10
## [77] png_0.1-7