motifmatchr的作用就是分析众多的序列和众多的motifs, 从中找到哪个序列包含哪个motif. 它的核心函数就是matchMotifs，最大特点就是快，因为它用的是MOODS C++库用于motif匹配。
尽管Bioconductor上也有很多工具能够做motif匹配，比如说Biostrings::mathcPWM, TFBSTools::searchSeq，但是motifmatchr更适合于分析许多不同的序列包含许多不同的motif。例如，当分析ChIP-seq或者ATAC-seq数据时, 你可能会想知道在哪个peak里有哪种类型的motif.
R包安装和加载
library(motifmatchr)
library(GenomicRanges)
matchMotifs
motifmatchr的核心函数是matchMotifs，因此了解这个函数的输入数据是什么格式就行了。必须的两个输入是
- 位置权重矩阵(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矩阵, 下一步根据这些文件构建出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()读取构建成一个GRanges对象。因为是提供的GRanges对象，那么还需要额外设置一个参数genome, 利用Biostrings::readDNAStringSet()读取一个参考基因组序列就行了。
或者用bedtools先根据bed文件提取数据，然后利用Biostrings::readDNAStringSet()读取
示例数据
我们以包中提供的数据为例，进行演示
加载实例的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函数提取匹配矩阵
motifMatches(motif_ix)
输出结果是一个稀疏矩阵
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")
其他参数
除了必须的motif信息和你的序列信息输入外，还有一些其他的参数可以做一些设置。
- 背景核苷酸频率, 默认是
bg=subject, 也就是你的输入数据作为背景，也可设置成genome或even - P值: 用于判断匹配是否足够可信的参数，默认是0.00005，基本上不用修改
- 输出信息:
matchMotifs可以返回三种可能输出，matches, scores 和 positions
总的来说，这个R包是一个比较简单的工具，比较容易上手使用。
