motifmatchr: 在R语言中分析peak中里是否有mot

2019-08-17 本文已影响9人 xuzhougeng

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包是一个比较简单的工具，比较容易上手使用。

我的个人博客 xuzhougeng.top

欢迎加入我的知识星球向我提问

知识星球