Skip to contents

Calculate motif enrichment in a set of sequences

Source: R/motif_enrichment.R
motif_enrichment.Rd

motif_enrichment() calculates motif enrichment relative to a set of background sequences using Analysis of Motif Enrichment (AME) from memes.

Usage

motif_enrichment(
  peak_input,
  motif,
  genome_build,
  out_dir = tempdir(),
  verbose = FALSE,
  meme_path = NULL,
  ...
)

Arguments

peak_input

Either a path to the narrowPeak file or a GRanges peak object generated by read_peak_file().

motif

An object of class universalmotif.

genome_build

The genome build that the peak sequences should be derived from.

out_dir

Location to save the 0-order background file along with the AME output files.

verbose

A logical indicating whether to print verbose messages while running the function. (default = FALSE)

meme_path

path to "meme/bin/" (default: NULL). Will use default search behavior as described in check_meme_install() if unset.

...

Arguments passed on to memes::runAme

method

default: fisher (allowed values: fisher, ranksum, pearson, spearman, 3dmhg, 4dmhg)

sequences

logical(1) add results from sequences.tsv to sequences list column to returned data.frame. Valid only if method = "fisher". See AME outputs webpage for more information (Default: FALSE).

silent

whether to suppress stdout (default: TRUE), useful for debugging.

Value

A list containing a FIMO results data frame and a numeric referring to the proportion of peaks with a motif.

See also

runAme

Examples

data("CTCF_TIP_peaks", package = "MotifPeeker")
data("motif_MA1102.3", package = "MotifPeeker")

# \donttest{
    res <- motif_enrichment(
        peak_input = CTCF_TIP_peaks,
        motif = motif_MA1102.3,
        genome_build =
            BSgenome.Hsapiens.UCSC.hg38::BSgenome.Hsapiens.UCSC.hg38
    )
    print(res)
#> $tp
#> [1] 73.00 40.11
#> 
#> $fp
#> [1] 286.00  26.19
#> 
#> $positive_peaks
#> GRanges object with 73 ranges and 7 metadata columns:
#>                   seqnames            ranges strand |            name     score
#>                      <Rle>         <IRanges>  <Rle> |     <character> <numeric>
#>    TIPseq_peak_63    chr10 69801154-69801320      * |  TIPseq_peak_63       109
#>    TIPseq_peak_90    chr10 70559156-70559588      * |  TIPseq_peak_90       122
#>   TIPseq_peak_157    chr10 73529363-73529638      * | TIPseq_peak_157       106
#>   TIPseq_peak_167    chr10 73957451-73957745      * | TIPseq_peak_167        73
#>    TIPseq_peak_39    chr10 68988538-68988838      * |  TIPseq_peak_39       111
#>               ...      ...               ...    ... .             ...       ...
#>   TIPseq_peak_105    chr10 71270297-71270663      * | TIPseq_peak_105       138
#>    TIPseq_peak_51    chr10 69507836-69508269      * |  TIPseq_peak_51       206
#>   TIPseq_peak_169    chr10 73966072-73966376      * | TIPseq_peak_169        70
#>    TIPseq_peak_44    chr10 69232402-69232647      * |  TIPseq_peak_44        45
#>    TIPseq_peak_46    chr10 69392521-69392898      * |  TIPseq_peak_46       132
#>                   signalValue    pValue    qValue      peak    summit
#>                     <numeric> <numeric> <numeric> <integer> <integer>
#>    TIPseq_peak_63     7.19424   15.0082  10.96220        54  69801208
#>    TIPseq_peak_90     8.60956   16.5013  12.29430       312  70559468
#>   TIPseq_peak_157     8.79456   14.6750  10.66540       105  73529468
#>   TIPseq_peak_167     6.58212   11.0334   7.35935       125  73957576
#>    TIPseq_peak_39     9.10857   15.2180  11.14710       168  68988706
#>               ...         ...       ...       ...       ...       ...
#>   TIPseq_peak_105     9.39681  18.29610  13.87650       113  71270410
#>    TIPseq_peak_51    11.20730  25.79100  20.65960       119  69507955
#>   TIPseq_peak_169     6.33687  10.63470   7.00211       168  73966240
#>    TIPseq_peak_44     5.62082   7.95049   4.55375       110  69232512
#>    TIPseq_peak_46     9.25369  17.56690  13.23280       225  69392746
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths
#> 
# }