Background nucleotide frequencies

In 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.

Log base and pseudocounts

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))

P value

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:

1. (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.

2. (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.

3. (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): name

head(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): name

head(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.81298

head(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             .              2

motif_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.97515
##   [2]     chr2   42773248-42773257      - |   7.08576
##   [3]     chr2 100184023-100184032      - |   7.73961
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths
## 
## $MA0107.1_RELA
## GRanges object with 1 range and 1 metadata column:
##       seqnames            ranges strand |     score
##          <Rle>         <IRanges>  <Rle> | <numeric>
##   [1]     chr1 76585933-76585942      + |   13.5563
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths
## 
## $MA0137.3_STAT1
## GRanges object with 3 ranges and 1 metadata column:
##       seqnames              ranges strand |     score
##          <Rle>           <IRanges>  <Rle> | <numeric>
##   [1]     chr2   42773087-42773097      - |   16.2940
##   [2]     chr2 100183852-100183862      + |   14.8130
##   [3]     chr2 100183852-100183862      - |   13.9536
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths