DREME原理和安装使用方法

  • A+
所属分类:Bioinformatics

写在前面

文献中常用的有DREME和HOMER,这次先搞定DREME,下次再写HOMER。

使用MEME套件中的DREME,用于鉴定meRIP-Seq数据中peak的motifmotif是序列中反复出现的由几个碱基组成的小片段,可能有特定的生物学意义。不过看meRIP-seq的文章里,这块的分析在结果中一般体现为一句话,交代找到的显著motif与已有报道一致,或者说多个发育阶段鉴定的保守motif是什么。

DREME的原理(http://meme-suite.org/doc/dreme-tutorial.html?man_type=web):

  • 输入文件:positive sequences和negative sequences(背景序列,如果没有提供,会根据positive sequences自动生成);
  • 寻找positive sequences中的3-8bp(默认长度)片段, 计算每个片段在positive和negative中出现的次数,通过Fisher's exact test 确定显著性,按照P-value排序成regular expression motifs.
  • 选取前100个motif,将其中1个碱基替换成ambiguity code,重新计算p-value;重复该过程,使得每个碱基都可以有ambiguity code,得到generalized RE motifs.(相当于consensus sequences)
  • 针对上一步得到的generalized RE motifs,计算每个motif在positive和negative中出现的次数,通过Fisher's exact test 确定显著性,按照设定的E-value值输出motif及相应的频率矩阵。

实战

01 安装

下载软件:http://meme-suite.org/doc/download.html

#install
mkdir meme
tar zxf meme_5.0.1.tar.gz 
cd meme_5.0.1
./configure --prefix=/root/cc/biosoft/bin/meme --with-url=http://meme-suite.org --enable-build-libxml2 --enable-build-libxslt
make
make test
make install

#加入PATH(~/.bash_profile)
PATH=/root/cc/biosoft/bin/meme/bin:$PATH
source ~/.bash_profile
#test
dreme 

#根据提示需要安装imagemagick
yum install ImageMagick
convert -version

02 实战
输入文件:fasta格式的序列文件

dreme -p input.fa -oc outDir -dna -png -eps -e 0.05 -t 18000 -m 20 
#-dna, -rna : 默认是dna
#-png, -eps : 输出图片格式
# -e, -t, -m : 是设置软件结束运行的参数,不设置的话,默认在e-value>e时停止,否则会一直运行下去。
# -e: e-value,0.05(default)  一直search,直到下一个motif 的e-value >e
# -t: 运行时间 一直search,直到达到这个时间
# -m: 找到的motif数目 一直search,直到找到的motif数目达到这个值
#还可以设置Core Motif Width:--mink 最小值;--maxk 最大值

设置为-dna时,可以识别U(等同于T),但输出LOGO时使用T.

设置为-rna时,仍然能识别T(等同于U),但输出LOGO时使用U.

同样的序列和参数,设置为-dna-rna,找到的motif是不一样的。Why? 没找到说明

是不是灰常简单?嘻嘻(#.#)

03 LOGO编辑

DREME生成的xml文件中,有每个motif的PWM matrix,里面有motif的每个位置上每种碱基的probability。可以用来直接画LOGO,方便的工具有如下几个。

seqLogo可以调整横纵坐标、碱基高度按比例/绝对大小展示等。

source("https://bioconductor.org/biocLite.R")
biocLite("seqLogo")
require(seqLogo)
mFile <- read.table("inp")
mFile
#  V1       V2 V3 V4 V5       V6
#1  0 0.000000  1  0  1 0.000000
#2  0 0.000000  0  1  0 0.000000
#3  1 0.907213  0  0  0 0.636825
#4  0 0.092787  0  0  0 0.363175

p<-makePWM(m)
seqLogo(p)  
seqLogo(p,ic.scale=FALSE)

seqLogo(p)

seqLogo(p,ic.scale=FALSE)

不足:目前seqLogo只能识别DNA alphabet,无法识别RNA,Protein

motifStack用途更广泛,这里仅用于把motif中的T替换成U。

需要安装ghostscript,然后将安装路径设置为环境变量R_GSCMD.

比如我安装的是64位的,安装在C:\Program Files\gs\gs9.23,按如下设置:

Sys.setenv(R_GSCMD=file.path("C:", "Program Files", "gs", 
                             "gs9.23", "bin", "gswin64c.exe"))

将矩阵中的T换成U。

source("https://bioconductor.org/biocLite.R")
biocLite("motifStack")
require(motifStack)
Sys.setenv(R_GSCMD=file.path("C:", "Program Files", "gs", 
                             "gs9.23", "bin", "gswin64c.exe"))

pcm=read.table("inp")
#0 0.000000  1  0  1 0.000000
#0 0.000000  0  1  0 0.000000
#1 0.907213  0  0  0 0.636825
#0 0.092787  0  0  0 0.363175

rownames(pcm) <- c("A","C","G","U")
motif <- new("pcm", mat=as.matrix(pcm), name="KD")
plot(motif,ncex=2,xaxis=F)  #ncex: fond size of name; xaxis 不显示横坐标

T--> U

04 Note

meRIP-seq的数据,分析得到的peak有时会跨区域。比如两个相邻的exon都是peak,但是得到的是一个包含这两个exon的区域,那么,提取序列时,应该提取的是真实的peak区域,即相应的exon区域。如果是exomePeak的结果,生成的BED12文件里包含block的信息,可以利用bedtools getfasta -split直接提取相应block的序列,再“串联”起来;如果是普通的BED文件(chr start end),那么就需要结合exon的坐标(比如GTF,GFF等注释文件),自己生成相应的BED12文件了。

avatar

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: