motifmatchr的作用就是分析眾多的序列和眾多的motifs, 從中找到哪個序列包含哪個motif. 它的核心函數(shù)就是matchMotifs，最大特點就是快，因為它用的是MOODS C++庫用于motif匹配。

盡管Bioconductor上也有很多工具能夠做motif匹配，比如說Biostrings::mathcPWM, TFBSTools::searchSeq，但是motifmatchr更適合于分析許多不同的序列包含許多不同的motif。例如，當(dāng)分析ChIP-seq或者ATAC-seq數(shù)據(jù)時, 你可能會想知道在哪個peak里有哪種類型的motif.

R包安裝和加載

library(motifmatchr)
library(GenomicRanges)

matchMotifs

motifmatchr的核心函數(shù)是matchMotifs，因此了解這個函數(shù)的輸入數(shù)據(jù)是什么格式就行了。必須的兩個輸入是

• 位置權(quán)重矩陣(position weight matrices, PWM)或位置頻率矩陣(position frequency matrices, PFM), 保存在PWMatrix, PFMatrix, PWMatrixList或PFMatrixList
• 一組基因組范圍(GenomicRanges或RangedSUmmarizedExperiment對象)或一組序列(DNAStringSet, DNAString 或 簡單的字符串向量)

PWM或PFM矩陣

Motif的PWM或PFM矩陣可以從JASPAR, CIS-BP下載。

例如，我從CIS-BP下載擬南芥Arabidopsis_thaliana_2019_08_17_2_32_am.zip，解壓縮后里有一個pwm_all_motifs文件夾，里面的文本文件存放的就是我所需要的PWM矩陣, 下一步根據(jù)這些文件構(gòu)建出matchMotifs的輸入對象

motif_dir <- "Arabidopsis_thaliana_2019_08_17_2_32_am/pwms_all_motifs"
PWList <- list()
for (file in list.files(motif_dir, pattern = ".txt")){
  df <- read.table(file.path(motif_dir, file), 
                   header = T,
                   row.names = 1)
  mt <- as.matrix(df)
  if (nrow(mt) ==0) next
  motif_id <- substr(file, 1,6)
  PWList[[motif_id]] <- PWMatrix(ID = motif_id, profileMatrix = t(mt))
}

PWMatrixLists <- do.call(PWMatrixList,PWList)

對于JASPAR，我們有更加方便的提取方法

library(JASPAR2018)
species <- "Arabidopsis thaliana"
collection <- "CORE"
opts <- list()
opts["species"] <- species
opts["collection"] <- collection
out  <- TFBSTools::getMatrixSet(JASPAR2018, opts)

if (!isTRUE(all.equal(TFBSTools::name(out), names(out)))) 
  names(out) <- paste(names(out), TFBSTools::name(out), 
                      sep = "_")
motif <- out

JASPAR的motif比較可靠，因此motif相對比較少。

給定范圍或序列

這部分的輸入就比較容易獲取，你可以提供MACS2的輸出BED，利用rtracklayer::import.bed()讀取構(gòu)建成一個GRanges對象。因為是提供的GRanges對象，那么還需要額外設(shè)置一個參數(shù)genome, 利用Biostrings::readDNAStringSet()讀取一個參考基因組序列就行了。

或者用bedtools先根據(jù)bed文件提取數(shù)據(jù)，然后利用Biostrings::readDNAStringSet()讀取

示例數(shù)據(jù)

我們以包中提供的數(shù)據(jù)為例，進行演示

加載實例的motif和GRanges對象

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

獲取motif在peak中的位置

motif_ix <- matchMotifs(example_motifs, peaks, genome = "hg19")

用motifMatches函數(shù)提取匹配矩陣

motifMatches(motif_ix)

輸出結(jié)果是一個稀疏矩陣

3 x 3 sparse Matrix of class "lgCMatrix"
     MA0599.1_KLF5 MA0107.1_RELA MA0137.3_STAT1
[1,]             |             |              .
[2,]             |             .              |
[3,]             |             .              |

其中的.就表示存在motif匹配。

我們還可以提取motif在peak中的位置

# Get motif positions within peaks for example motifs in peaks 
motif_ix <- matchMotifs(example_motifs, peaks, genome = "hg19",
                         out = "positions") 

其他參數(shù)

除了必須的motif信息和你的序列信息輸入外，還有一些其他的參數(shù)可以做一些設(shè)置。

• 背景核苷酸頻率, 默認(rèn)是bg=subject, 也就是你的輸入數(shù)據(jù)作為背景，也可設(shè)置成genome或even
• P值: 用于判斷匹配是否足夠可信的參數(shù)，默認(rèn)是0.00005，基本上不用修改
• 輸出信息: matchMotifs可以返回三種可能輸出，matches, scores 和 positions

總的來說，這個R包是一個比較簡單的工具，比較容易上手使用。