Motif预测

Motif预测
在许多生物学研究过程中，我们需要对真核和原核生物的转录调控，蛋白质结构活性位点，以及DNA、RNA的酶切位点进行识别。

在这些实际过程中，我们通常会碰到下面两个问题：
（1）我们通常需要从未知的数据尽可能多的发掘有用的信息，也就是说，在探索实验过程中，由于科研工作者一开始并不知道哪些信息是我们真正“有用”的，在没有其他太多的信息可以利用的时候，一个很自然的想法就是我们是否能知道哪些因素是这些数据共有部分，也就是我们通常所谓的数据的公共特征提取，或者叫motif的预测或模式识别，这里的“模式”（motif）可以简单的理解为特定数据的共同特征。

（2）另一方面，假设当我们已经知道了某个特征，我们需要把具有这些特征的数据都收集起来，那么，我们需要在大量数据中挑选符合我们要求地数据，这个过程就是所谓的“模式匹配”。

“模式识别”和“模式匹配”是生物信息学辅助实验生物学的一个重要手段。

关于这些模式的具体形式可以是相当广泛的。

从广义上讲，即使是多序列比对寻找保守区、蛋白质结构预测等问题，也可以纳入模式识别的范畴。

但下面我们要介绍是以围绕转录因子结合位点（TFBS）为代表的一类模式识别的程序，是大家传统意义上的模式识别问题。

这一类模式（motif）的最大的特点就是大家共有特征较短，一般一个TFBS位点的长度在5-20bp左右，而且信号比较灵活多变，所以这些调控元件在通常情况下不适合直接使用多序列联配的方式来寻找，而是需要一些专门的算法来解决这个问题。

下面我们介绍几个较著名的模式（motif）预测软件。

7.3.1 MEME/MAST系统
MEME和MAST是由T.L.Bailey、Charles Elkan和Bill Noble合作开发的一套搜索motif的程序组合套件。

是目前生物信息学领域，对motif预测方面最著名的程序之一。

这个组件的两个程序分别执行的是motif的预测和搜索两个不同的功能。

其中MEME的全称是Multi EM for motif elicitation，是一个基于EM算法的一个motif 预测程序。

而MAST的全称是Motif Alignment & search Tool，是一个motif搜索组件。

下载：
这个软件包可以从MEME的官方网站免费下载，我们这里使用的最新版本是version 3.5.4，其源代码可以从如下链接下载：
/downloads/meme_3.5.4.tar.gz
安装：
这里介绍在Linux/unix下的安装步骤：
1、把源代码文件上载本地的Linux/Unix指定目录:
/mnt/disk1/motif_workplace/
2、解压：
$ gunzip -c meme_3.5.4.tar.gz | tar xvf -
$ cd meme_3.5.4
3、配置：通常情况下，可以使用默认参数。

如果需要指定安装路径，可以使用参数--prefix=your-install-path
$ ./configure或$ ./configure --prefix=your-install-path
如果配置成功，会出现以下画面：
========================
Configuration parameters
========================
Install path: your-install-path
Install UID: usrid
Version: 3.5.4
C compiler: gcc
C compiler flags: -Wall -DUNIX -D__USE_FIXED_PROTOTYPES__ -O3 Linker: /usr/bin/ld
Special Libs: -lm
MPICC
MPIINC
MPILIBDIR
MPIFLAGS -DPARALLEL -O
MPIRUN
URL: http://usr_home/meme
PERL: /usr/bin/perl
MEME_LOGS: ${prefix}/LOGS
Run the following commands to compile, test and install meme: make
make test
make install
4、编译和安装：
在命令行键入make开始编译，
$ make
当编译完成后，键入命令：make test，测试，
$ make test
当返回提示：
Meme test for
---------------------------------
Dataset Model result
---------------------------------
crp0 oops OK
crp0 zoops OK
crp0 tcm OK
......
PASS: runcheck
==================
All 1 tests passed
==================
......
说明测试成功，可以键入命令：make install进行安装。

$ make install
其中测试画面中出现的三个模型oops、zoops和tcm是MEMmotif预测中最重要的选择参数，分别表示对任意给定的motif来讲，每条序列当且仅当出现一次、每条序列至多有一个以及每条序列出现两个或两个以上三种情况。

对于一个给定的序列，我们可以根据我们需要研究的问题，选择这三个模型中的一个。

关于meme参数的选择将在后面参数一部分详细讨论。

使用：
当我们在命令行键入meme命令，可以看到meme的参数列表和帮助文档。

$ meme
USAGE:
meme <dataset> [optional arguments]
<dataset> file containing sequences in FASTA format
[-h] print this message
[-dna] sequences use DNA alphabet
[-protein] sequences use protein alphabet
[-mod oops|zoops|anr] distribution of motifs
[-nmotifs <nmotifs>] maximum number of motifs to find
[-evt <ev>] stop if motif E-value greater than <evt>
[-nsites <sites>] number of sites for each motif
[-minsites <minsites>] minimum number of sites for each motif
[-maxsites <maxsites>] maximum number of sites for each motif
[-wnsites <wnsites>] weight on expected number of sites
[-w <w>] motif width
[-minw <minw>] minumum motif width
[-maxw <maxw>] maximum motif width
[-nomatrim] do not adjust motif width using multiple alignment [-wg <wg>] gap opening cost for multiple alignments
[-ws <ws>] gap extension cost for multiple alignments
[-noendgaps] do not count end gaps in multiple alignments
[-bfile <bfile>] name of background Markov model file
[-revcomp] allow sites on + or - DNA strands
[-pal] force palindromes (requires -dna)
[-maxiter <maxiter>] maximum EM iterations to run
[-distance <distance>] EM convergence criterion
[-prior dirichlet|dmix|mega|megap|addone]
type of prior to use
[-b <b>] strength of the prior
[-plib <plib>] name of Dirichlet prior file
[-spfuzz <spfuzz>] fuzziness of sequence to theta mapping
[-spmap uni|pam] starting point seq to theta mapping type
[-cons <cons>] consensus sequence to start EM from
[-text] output in text format (default is HTML)
[-maxsize <maxsize>] maximum dataset size in characters
[-nostatus] do not print progress reports to terminal
[-p <np>] use parallel version with <np> processors
[-time <t>] quit before <t> CPU seconds consumed
[-sf <sf>] print <sf> as name of sequence file
... ...
当我们在命令行键入mast命令，可以看到mast的参数列表和帮助文档。

$ mast
USAGE:
mast <mfile> [optional arguments ...]
<mfile> file containing motifs to use; may be a MEME output
file or a file with the format given below
[<database>] or
[-d <database>] database to search with motifs or
[-stdin] read database from standard input;
Default: reads database specified inside <mfile>
[-c <count>] only use the first <count> motifs
[-a <alphabet>] <mfile> is assumed to contain motifs in the
format output by make_logodds
and <alphabet> is their alphabet; -d <database>
or -stdin must be specified when this option is used
[-stdout] print output to standard output instead of file
[-text] output in text (ASCII) format;
(default: hypertext (HTML) format)
[-sep] score reverse complement DNA strand as a separate sequence [-norc] do not score reverse complement DNA strand
[-dna] translate DNA sequences to protein
[-comp] adjust p-values and E-values for sequence composition
[-rank <rank>] print results starting with <rank> best (default: 1)
[-smax <smax>] print results for no more than <smax> sequences (default: all) [-ev <ev>] print results for sequences with E-value < <ev> (default: 10) [-mt <mt>] show motif matches with p-value < mt (default: 0.0001)
[-w] show weak matches (mt<p-value<mt*10) in angle brackets [-bfile <bfile>] read background frequencies from <bfile>
[-seqp] use SEQUENCE p-values for motif thresholds
(default: use POSITION p-values)
[-mf <mf>] print <mf> as motif file name
[-df <df>] print <df> as database name
[-minseqs <minseqs>] lower bound on number of sequences in db
[-mev <mev>]+ use only motifs with E-values less than <mev>
[-m <m>]+ use only motif(s) number <m> (overrides -mev)
[-diag <diag>] nominal order and spacing of motifs
[-best] include only the best motif in diagrams
[-remcorr] remove highly correlated motifs from query
[-brief] brief output--do not print documentation
[-b] print only sections I and II
[-nostatus] do not print progress report
[-hit_list] print machine-readable list of all hits only; implies-text
... ...
下面我们开始逐一解释这些基本参数的含义以及在实际应用中需要注意哪些问题。

MEME参数的说明：
1）基本参数篇
我们可以看到，meme程序的参数很多，为了更好的说明问题，meme把参数大致地分几类，我们按照meme文档里的分类模式，简要说明：
1. <数据文件>：meme输入文件格式是fasta格式，主要形式如下：
>seq1
GDIFYPGYCPDVKPVNDFDLSAFAGAWHEIAK
>seq2
GDMFCPGYCPDVKPVGDFDLSAFAGAWHELAK
>seq3
QKVAGTWYSLAMAASDISLLDAQSAPLRVYVEELKPTPEGDLEILLQKW
可以是DNA序列的fasta文件，也可以蛋白质序列的fasta格式文件。

2.可选的几个参数类：①字母表类；②分布类；③搜索类；④系统类
①字母表参数类：
-dna/-protein，说明输入序列类型是DNA还是蛋白质。

②分布模型类：
-mod<string>有3个可选的值，就是前面提到的3个可选模型：
oops：(One Occurrence Per Sequence) 意思是每种motif在每条序列中只都出现一次。

zoops：(Zero or One Occurrence Per Sequence)，意思是每种motif在中至多出现一次。

anr：(Any Number of Repetitions)，意思是每种motif可以出现任意数目。

这个其实是前面提到的tcm模型和zoops功能的一个组合。

③搜索参数类
MEME预测的motif的过程其实是一个寻找给定目标函数最优的过程，通过目标函数计算每个可能的motif的对数似然比，近似用E值表示，所以候选种子的E值越小，表示这个种子是真正的motif的可能性越大，所以结果文件会按E 值从小到大的顺序输出。

所以下面的一些有关输入输出的参数是可以根据我们实际需要，进行调整。

-nmotifs <n>：选择输出motif个数n，表示meme的运行结果文件输出n个不同motif，这个参数的默认值是1。

-evt <p>：E-value的阈值，当这个参数为p，表示我们只输出E值比p小的那些motif。

所谓E值是一个描述motif模型可信度的一个标志，一般而言，E值越小，结果可信度越高。

这个参数的默认值是无穷大，也就是在默认状态下，这个参数不起作用。

-nsite<n>：表示每个motif出现的期望。

-minsites <n>：表示每个motif出现的期望的下限。

-maxsites <n>：表示每个motif出现的期望的上限。

这3个参数不是相互独立的，只能取-nsite和其他两个参数不能同时选择。

当我们取-nsite时，只有那个给定数目的种子作为程序候选的motif进行运算。

而-minsites和-maxsites参数则给出这个期望数目的一个下限和上限约束的范围。

只有落在这个范围内的种子才给予考虑。

注意，在OOPS参数下，这些参数的设置为无效设置。

而在其他两个模型参数的条件下，默认值为
-minsites：sqrt（序列数）
-maxsites：当ZOOPS时，默认为序列数n
当anr时，默认为min（序列数×5，50）。

-wnsites <n>: 表示nsite初始的权重。

它是控制种子满足nsites或minsites 和maxsites限制的一个贡献值，这个权值是[0,1)的一个值，这个值越大，表示motif倾向与满足限制的可能越大。

默认值是用0.8。

-w <n>：表示motif种子的宽度。

-minw <n>：表示motif种子宽度的下限
-maxw <n>：表示motif种子宽度的上限
也就是说，当-w给定，程序只尝试宽度为-w的种子，否则，尝试motif种子宽度在-minw和-maxw之间的值。

默认的-minw是8 -maxw是50。

-text：输出文本格式（默认的输出格式是HTML文件）
-maxsize <maxsize>：最大的数据集大小
-sf <sf>：打印<sf>作为序列文件名
[-nostatus]：不要在终端输出打印报告
[-time <t>]：运行在<t>时间前自动退出程序
MEME高级参数技巧篇
-nomatrim：这个参数的选用令MEME跳过用多序列联配截短的步骤。

-wg <a>：空位罚分，默认值是11
-ws <a>：空格罚分，默认值是1
-noendgaps：对最后的空位不罚分，默认值是对最后的空位罚分。

-bfile <bfile>：背景概率分布文件。

-revcomp：翻转序列，即motif可以在本身序列上，也可以在补链上
-pal：考虑回文结构
-maxiter <n>：表示EM迭代次数，默认是50次
-distance <a>：表示收敛规则，表示两个连续的频率矩阵的欧几里得距离小于这个值是，跳出迭代。

默认值是0.001。

-prior <string>：先验概率模型参数
dirichlet：简单Dirichlet先验，作为-dna和-alph的默认先验。

dmix：混合Dirichlets先验，-protein的默认分布。

mega：方差非常小的混合Dirichlets分布，方差除以数据集变换。

megap：在最后一次迭代用dmix参数，其他使用mega参数。

addone：给每一个观察值+1。

（laplace 法则）。

-b <a>：先验模型的强度。

<a> = 0：表示使用固有的先验强度。

默认值通常：
0.01：对dirichlet先验，<a>=0.01。

0 if：对于dmix而言，默认<a>=0。

-plib <string>：存放Dirichlet先验的文件，文件格式如prior30.plib所示。

-spfuzz <a>：作模糊操作具体如spmap；
-spmap <string>：使用映射函数，
uni：增加一个先验，并把一个子串转为theta的估计，这里的theta是EM算法中的一个估计参数。

默认的参数是-spfuzz <a>：0.5
pam：把PAM矩阵的子串转化为theta的估计。

默认的参数-spfuzz <a>：120 (PAM 120)
-cons <string>：废除采样的起始点，使用给定的字符串作种子。

这个参数在一直特定motif的基本形式时会很有用。

④系统参数：
最后几个参数是和大型机系统有关，MEME程序可以在服务器上支持并行运算，因此最后几个参数只有在一定的服务器环境下可以使用。

[-p <np>]：用<np>个CPU并行运算
以上是MEME各参数的说明，对于参数的选择，需要大家根据具体问题灵活调整，而对于这里提到的高级参数部分，如果大家对MEME的概率模型不熟悉的话，建议使用默认参数即可。

例1：假设我们有5条水稻的TSS上游调控序列，序列文件rice9311_sample.fa，我们需要寻找他们的motif，并把结果保存在rice_meme.html里面。

在命令行中输入命令如下：
$ meme rice9311_sample.fa -dna -nmotifs 4 -mod zoops -minw 5 -maxw 15 >rice_meme_htmlFormat.html
<网页形式见rice_meme.html>
或
$ meme rice9311_sample.fa -dna -nmotifs 4 -mod zoops -minw 5 -maxw 15 -text>rice_meme_txtFormat.out
<文本格式>
下表列出rice_meme_txtFormat.out的一部分文件如下，我们将分段解析结果文件包含的内容，从整体上看，MEME的结果文件主要由文件头、基本数据信息、输入命令信息、motif结果条目以及motif综合信息五大模块组成。

下面我们以例1运行的结果为例，解析各区段的作用。

（1）头文件：
我们可以看到，MEME结果的文件头可以分两大部分，一部分包含MEME 版本信息当前版本是（v3.5.4）以及官方主页。

********************************************************************* MEME - Motif discovery tool
********************************************************************* MEME version 3.5.4 (Release date: 3.5.4)
For further information on how to interpret these results or to get
a copy of the MEME software please access .
This file may be used as input to the MAST algorithm for searching
sequence databases for matches to groups of motifs. MAST is available
for interactive use and downloading at .
********************************************************************* 头文件的另一部分是关于文献引用的说明，在论文中如果使用MEME的结果，可以应用下面给出的T.Bailey1994的文献。

希望使用MEME的读者能在论文中正确的引用MEME的文献。

*********************************************************************
REFERENCE
********************************************************************* If you use this program in your research, please cite:
Timothy L. Bailey and Charles Elkan,
"Fitting a mixture model by expectation maximization to discover
motifs in biopolymers", Proceedings of the Second International
Conference on Intelligent Systems for Molecular Biology, pp. 28-36,
AAAI Press, Menlo Park, California, 1994.
********************************************************************* （2）数据集基本信息。

这段数据保存了输入数据集的一些基本的信息。

包括输入序列文件名和字符表以及每条输入序列的序列名称, 序列的权重，和序列长度。

其中这里的序列权重是表示这条序列在预测motif中发挥作用大小的一个度量。

********************************************************************* TRAINING SET
********************************************************************* DATAFILE= rice9311_sample.fa
ALPHABET= ACGT
Sequence name Weight Length Sequence name Weight Length
------------- ------ ------ ------------- ------ ------
BGK02457 1.0000 600 BGK03149 1.0000 600
BGK04166 1.0000 600 BGK01655 1.0000 600
BGK03537 1.0000 600 BGK02838 1.0000 600
BGK04759 1.0000 600
********************************************************************* 我们可以看到在这里我们的输入序列文件名：rice9311_sample.fa。

因为输入是单纯的DNA序列，序列中并没有出现DNA的简并形式，所以这里的字符表是{A,C,G,T}，我们还可以看到我们输入的序列名称，BGKXXXX，一般情况下，序列的权重是1，现在我们的输入序列的长度、都是600bp。

（3）命令行参数。

这一部分主要打印运行MEME时的人工参数以及程序默认的参数。

********************************************************************* COMMAND LINE SUMMARY
********************************************************************* This information can also be useful in the event you wish to report a problem with the MEME software.
command: meme rice9311_sample.fa -dna -nmotifs 4 -mod zoops -minw 5 -maxw 15 -dir
/home/goodgood
model: mod= zoops nmotifs= 4 evt= inf object
function= E-value of product of p-values
width: minw= 5 maxw= 15 minic= 0.00
width: wg= 11 ws= 1 endgaps= yes
nsites: minsites= 2 maxsites= 7 wnsites= 0.8
theta: prob= 1 spmap= uni spfuzz= 0.5
em: prior= dirichlet b= 0.01 maxiter= 50 distance=
1e-05
data: n= 4200 N= 7
strands: +
sample: seed= 0 seqfrac= 1
Letter frequencies in dataset:
A 0.267 C 0.277 G 0.211 T 0.245
Background letter frequencies (from dataset with add-one prior applied):
A 0.267 C 0.277 G 0.211 T 0.245
********************************************************************* 我们可以看到上表列出的参数中,在我们命令中，我们使用的是比较常用的参数，除了文件名和序列类型（-dna/-protein）以外，我们还使用了模型参数-mod zoops，在这里值得注意的是，这个模型参数的指定是MEME结果能否与生物学意义紧密联系的重要参数之一。

一般通常没有特殊假设的情况下使用ZOOPS模型比较合适。

总之，恰当模型的选择将是我们解决问题捷径。

另一个十分重要的参数是[-revcomp]，这个参数是用来指定程序是否让序列的补链上也参与motif预测。

因为这个参数的选择是和我们的数据集以及我们期望解决的生物学问题紧密关联，所以这个参数是不可以默认的。

在我们这个例子的问题中，由于我们的数据是水稻基因5'端上游启动子序列，所以我们并不需要补链信息，所以我们没有使用-revcomp参数。

关于其他的程序参数如何选择，很多时候并没太严格的规定，希望大家能在使用中仔细揣摩，因为参数的选择问题有些时候更象是一门艺术。

结果：
这一区段主要是MEME的预测得到的motif的形式和描述和评估。

在这个部分，结果文件可以根据找到的motif的E值由小到大的顺序，显示我们在参数中指定输出的最显著的若干的motif。

在这个区段中，MEME返回了多种形式的motif，并且给出了每个motif的统计显著性的E值和P值。

对于这个区段里列出每个motif，都会包含motif的基本信息和描述、位置p 值的排序、块的图解和motif的块形式、motif的位置特异的打分矩阵和位置特异概率矩阵形式以及正则表达式形式。

[1] motif 1 的基本信息。

********************************************************************* MOTIF 1 width = 13 sites = 7 llr = 88 E-value = 7.9e-001
********************************************************************* 这里给出这个motif的最基本信息，包括motif的宽度、位点的匹配数目、对数似然比(llr)和这个motif的E值。

这里llr的基本公式可以用llr = log (Pr(sites | motif) / Pr(sites | back))计算，其中Pr(sites|motif)条件概率，表示在motif条件下匹配的概率；Pr(sites | back)表示随机背景条件下出现位点匹配的概率。

这里的llr 值和E值都是反映这个motif是否是真正motif的指标，其中llr越大，E值越小，则说明这个motif越可靠。

[2]motif 1 的描述
----------------------------------------------------------------------------
Motif 1 Description
----------------------------------------------------------------------------
Simplified A a143a:a:4:69:
pos.-specific C :13:::::::::1
probability G :637:a:966419
matrix T :1:::::1:4:::
bits 2.2 *
2.0 *
1.8 * ***
1.6 * **** *
Information 1.3 * ***** **
content 1.1 * **********
(18.0 bits) 0.9 * **********
0.7 * **********
0.4 *************
0.2 *************
0.0 -------------
Multilevel AGAGAGAGGGAAG
consensus CA ATG
sequence G
----------------------------------------------------------------------------
从上面的描述片断，我们可以看到三种不同形式的motif的表示形式反映motif不同位置的保守性，一种是用简单的频数概率矩阵，一种是用motif各位置信息熵值，熵值越大，说明这个位点的保守性越高。

我们看到motif从左往右的第5位、第6位和第7位的信息熵较大分别是是1.8、2.2和1.8，从频度矩阵上也可以看到这几个位置最保守。

这里motif表示的第三种形式是多水平的一致序列形式，可以直观的看到那些位置的保守碱基是哪些。

[3]位置P值的排序
这一部分列出了这个motif在序列上的匹配位置，并根据它的匹配的P值大小，按P值从小到大的顺序排列序列，也就是根据匹配的可靠程度，排列序列，反映序列和motif的关系。

这里的P值可以理解为假阳性率，P值越小，结果越可靠。

----------------------------------------------------------------------------
Motif 1 sites sorted by position p-value
----------------------------------------------------------------------------
Sequence name Start P-value Site
------------- ----- --------- -------------
BGK03149 547 1.21e-08 CTAGCGAGCG AGAGAGAGGGAAG GGTTGCGACT
BGK04759 510 3.58e-08 GGATACAGGT AGAGAGAGGTGAG AAGGCAGTGG
BGK04166 68 9.35e-07 CAATAGAGTT ATAGAGAGATAAG AGAAGAGGTA
BGK02457 558 1.40e-06 TTTCGTGAGC AGGAAGAGGGGGG GGGGGGGGGG
BGK03537 102 1.62e-06 ACACCAAGTG ACCGAGAGATAAG GCTCATACAG
BGK01655 278 2.51e-06 GACTAGACGG AGCAAGATGGAAG TACACAGTCA
BGK02838 315 4.48e-06 TTAAATAAAA AAGGAGAGAGGAC TGACAATTTG
----------------------------------------------------------------------------
[4] motif块的图解
显示motif在各序列上的位置，并在图上显示，直观地反映motif序列的分布。

----------------------------------------------------------------------------
Motif 1 block diagrams
----------------------------------------------------------------------------
SEQUENCE NAME POSITION P-VALUE MOTIF DIAGRAM
------------- ---------------- -------------
BGK03149 1.2e-08 546_[1]_41
BGK04759 3.6e-08 509_[1]_78
BGK04166 9.4e-07 67_[1]_520
BGK02457 1.4e-06 557_[1]_30
BGK03537 1.6e-06 101_[1]_486
BGK01655 2.5e-06 277_[1]_310
BGK02838 4.5e-06 314_[1]_273
----------------------------------------------------------------------------
例如：BGK03149序列上的546_[1]_41表示motif [1]是夹在序列的从5’到3'的第546位开始到从3’到5'的第41 位之间的长度恰好是13nt的位置。

[5] motif的块表示
像多序列联配那样，motif位点的匹配也可以对齐位置，成为一个块，这是motif的联配形式，也可以反映位置的保守性，在某些情况下，我们也可以使用这种motif块形式来。

比如有个简单的motif logo画图程序weblogo就可以根据块形式motif输入，绘制漂亮的motif logo。

----------------------------------------------------------------------------
Motif 1 in BLOCKS format
----------------------------------------------------------------------------
BL MOTIF 1 width=13 seqs=7
BGK03149 ( 547) AGAGAGAGGGAAG 1
BGK04759 ( 510) AGAGAGAGGTGAG 1
BGK04166 ( 68) ATAGAGAGATAAG 1
BGK02457 ( 558) AGGAAGAGGGGGG 1
BGK03537 ( 102) ACCGAGAGATAAG 1
BGK01655 ( 278) AGCAAGATGGAAG 1
BGK02838 ( 315) AAGGAGAGAGGAC 1
//
----------------------------------------------------------------------------
[6]位置特异打分矩阵
----------------------------------------------------------------------------
Motif 1 position-specific scoring matrix
----------------------------------------------------------------------------
log-odds matrix: alength= 4 w= 13 n= 4116 bayes= 9.19722 E= 7.9e-001 190 -945 -945 -945
-90 -95 143 -78
68 4 44 -945
10 -945 176 -945
190 -945 -945 -945
-945 -945 224 -945
190 -945 -945 -945
-945 -945 202 -78
68 -945 143 -945
-945 -945 143 81
110 -945 102 -945
168 -945 -56 -945
-945 -95 202 -945
----------------------------------------------------------------------------
MEME 中的位置特异打分矩阵的分值的计算公式如下：
另x是{A,C,G,T}中任意一个碱基，ns是位点的总匹配数，这里ns=7，这里是n_{ix}是在位点的i列上x出现的次数，另位点第i位为碱基x的分记作score(i,x)有，若n_{ix}不为0，则score（i,x）=100*log2(n_{ix}/n/p(x))。

否则，score(i,x)=100*log2(0.01/ns)。

所以，当i=2为例，score（1,A）=100*log(（1/7）/0.267)约等于-90。

[7]位置特异概率矩阵
----------------------------------------------------------------------------
Motif 1 position-specific probability matrix
----------------------------------------------------------------------------
letter -probability matrix: alength= 4 w= 13 nsites= 7 E= 7.9e-001 1.000000 0.000000 0.000000 0.000000
0.142857 0.142857 0.571429 0.142857
0.428571 0.285714 0.285714 0.000000
0.285714 0.000000 0.714286 0.000000
1.000000 0.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
1.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.857143 0.142857
0.428571 0.000000 0.571429 0.000000
0.000000 0.000000 0.571429 0.428571
0.571429 0.000000 0.428571 0.000000
0.857143 0.000000 0.142857 0.000000
0.000000 0.142857 0.857143 0.000000
----------------------------------------------------------------------------
位置特异的概率矩阵的计算十分简单，motif的第i位为碱基x的概率：P(i,x)=n_{ix}/ns
[8]motif的正则表达式形式
为了方面人们使用脚本大规模处理数据，MEME还给出简单的正则表达式形式，正则表达式是计算机领域字符串匹配一个常用的形式，它有一套明确的表示方式来处里各种模式。

----------------------------------------------------------------------------
Motif 1 regular expression
----------------------------------------------------------------------------
AG[ACG][GA]AGAG[GA][GT][AG]AG
----------------------------------------------------------------------------
这里[...]表示方括号内部的几个字符都是允许的，以AG[ACG]A这个正则表达式为例，它的第三位字符可以是A或C或G，也就是说AGAA或AGCA或AGGA 这三个字串都是符合上AG[ACG]A的语法。

Time 19.47 secs.
... ...
（这里省略其他的motif的类似的描述片段）
... ...
motif结果综合信息：
最后给出一个结果文件包含的各个motif在序列上的位置的图，这一区段可以直观地反映各motif之间的位置关系可以作为最基本组合分析的第一个环节。

********************************************************************* SUMMARY OF MOTIFS
********************************************************************* ----------------------------------------------------------------------------
Combined block diagrams: non-overlapping sites with p-value < 0.0001
----------------------------------------------------------------------------
SEQUENCE NAME COMBINED P-VALUE MOTIF DIAGRAM
------------- ---------------- -------------
BGK02457 2.00e-10 79_[2(3.26e-07)]_51_[3(2.97e-09)]_397_[1(1.40e-06)]_30
BGK03149 1.46e-13
31_[1(7.50e-05)]_55_[2(6.96e-10)]_[1(2.05e-05)]_15_[1(8.43e-05)]_246_[4(7.63e-05 )]_80_[3(1.83e-07)]_43_[1(1.21e-08)]_41
BGK04166 5.37e-09
22_[4(7.63e-05)]_38_[1(9.35e-07)]_202_[1(5.22e-06)]_63_[3(8.08e-07)]_100_[2(2.1 0e-07)]_112
BGK01655 1.53e-05 117_[4(7.63e-05)]_153_[1(2.51e-06)]_16_[3(3.89e-07)]_279 BGK03537 4.20e-08 9_[2(4.70e-07)]_77_[1(1.62e-06)]_229_[4(7.63e-05)]_189_[3(2.11e-06)]_46
BGK02838 2.95e-08 158_[2(5.16e-05)]_25_[2(2.52e-08)]_101_[1(4.48e-06)]_189_[3(5.10e-06)]_69
BGK04759 6.89e-07 155_[4(7.63e-05)]_100_[3(1.22e-06)]_232_[1(3.58e-08)]_78
----------------------------------------------------------------------------
********************************************************************* 从上面的区段可以看到这一部分有几个结果数据、除序列名以外，还包括给定序列的组合P值，以及一个motif的分布示意图。

以BGK02457为例，组合P 值是 2.00e-10，表示在随机背景模型下生成符合这种形式的组合的概率是2.00e-10，关于motif示意图。

79_[2(3.26e-07)]_51_[3(2.97e-09)]_397_[1(1.40e-06)]_30。

这个示意图的意思是BGK02457序列从5’到3'方向经过79nt后，有motif 2的匹配,图上用表示，（3.26e-07）表示的是这个motif的位置p值。

见前面motif描述中位置p值部分。

********************************************************************* Stopped because nmotifs = 4 reached.
********************************************************************* CPU: shangao
********************************************************************* 在这里，我们简单地解释了MEME文本格式的结果文件的主要项目，除文本形式以外，MEME还提供更漂亮的HTML格式的结果，其结构和文本格式几乎一样，这里不再重复解释。

有兴趣的读者可以自行比较。

对于这里使用的完整的结果可以见附件：
rice_meme_txtFormat.out
rice_meme_htmlFormat.html
如果需要对结果文件各项更详细的解释，可以参看MEME的网上文档：/meme/meme-output-example.html#explanation
（1）根据例1的数据集，比较MEME双链motif预测和单链motif预测的结果差异，各举一个实际例子（提示：用MEME参数-revcomp）
（2）利用MEME挑选合适的模型寻找人和小鼠的可能的转录因子结合位点（提示：promoter序列可以从EPD下载，core promoter区大约是位于TSS上游150，到TSS下游50）。

EPD地址：ftp://ftp.epd.unil.ch/pub/databases/epd/ （3）调整各种概率模型（ZOOPS，OOPS，…），以及不同的字串长度，找到P值最小的motif，并尝试评估不同概率模型对结果的影响。

Motif的搜索：MAST
接着，我们可以利用MAST的搜索工具对预测的motif进行搜索，最简单的请况是使用命令。

我们对MAST常用其他参数说明，具体含义见下表：
mast命令的基本格式：mast <mfile> [optional arguments ...]
参数如下：
<mfile>：保存有motif的文件，可以是meme的结果文件，也可以下面提到的某些格式文件。

[<database>]或[-d <database>]：搜索motif的数据集。

[-stdin]：从标准输入中读取数据库。

默认是从meme结果文件作预测的那个数据集。

[-a <alphabet>]：字符集。

[-stdout]：使用标准输出输出结果。

[-text]：输出文本格式（默认是html格式输出）。

[-sep]：让DNA互补链作为单独的序列进行打分。

[-norc]：不对DNA互补链进行打分。

[-dna]：把DNA序列翻译成蛋白。

[-comp]：根据序列组成调整p值和E值。

[-rank <rank>]：打印最好的<rank>个结果（默认值：1）
[-smax <smax>]：打印至多<smax>条序列（默认为全部序列）
[-ev <ev>]：打印E值小于<ev>的序列（默认的<ev>是10）
[-mt <mt>]：显示p值小于<mt>的motif （默认值：0.0001）
[-w]：在尖括号显示弱motif（mt<p-value<mt*10）in angle brackets
[-bfile <bfile>]：从<bfile>文件中读取背景频率
[-seqp]：使用序列的p值作为motif的阈值。

（default: use POSITION p-values）[-mf <mf>]：打印<mf>作为motif文件名
[-df <df>]：打印<df>作为数据库名
[-minseqs <minseqs>]：数据库中的序列数的下界
[-mev <mev>]+：只使用E值小于<mev>的motif
[-m <m>]+：使用motif数目为<m>（忽略-mev）
[-diag <diag>]：motif的次序和间隔
[-best]：只在图表中包含最好的motif
[-remcorr]：去除高度相关的motif
[-brief]：简洁输出--不打印文档
[-b]：只打印第一节和第二节
[-nostatus]：不打印程序报告
[-hit_list]：打印所有匹配的机器可读列表；包括-text
实例2：
用实例1中给定的水稻数据集，搜索MEME预测的motif。

基本命令如下：
$ mast rice_meme_txtFormat.out rice9311_sample.fa -brief -text
<文本格式>
我们把前面MEME的结果文件输入MAST。

和MEME的结果文件类似，MAST结果可以也可以分文件头、数据、motif、高分序列、motif图和序列注释五部分。

（1）头文件。

和MEME一样，MAST结果的头文件，交代了MAST版本必要的版本信息、官方主页和文献引用。

********************************************************************* MAST - Motif Alignment and Search Tool
********************************************************************* MAST version 3.5.4 (Release date: 3.5.4)
For further information on how to interpret these results or to get
a copy of the MAST software please access .
********************************************************************* REFERENCE
********************************************************************* If you use this program in your research, please cite:
Timothy L. Bailey and Michael Gribskov,
"Combining evidence using p-values: application to sequence homology searches", Bioinformatics, 14(48-54), 1998.
********************************************************************* （2）数据库和motif。

这个部分主要介绍用于搜索motif的数据集的一些基本信息，计算给出motif的预处理信息。

********************************************************************* DATABASE AND MOTIFS
********************************************************************* DATABASE rice9311_sample.fa (nucleotide)
Last updated on Tue Jan 30 22:24:44 2007
Database contains 7 sequences, 4200 residues
Scores for positive and reverse complement strands are combined.
MOTIFS rice_meme_txtBriefFormat.out (nucleotide)
MOTIF WIDTH BEST POSSIBLE MATCH
----- ----- -------------------
1 13 AGAGAGAGGGAAG
2 15 TGGTTTTTAGGTGGA
3 15 CTCGCTTCCCCCGCA
4 7 AGAAAAA
PAIRWISE MOTIF CORRELATIONS:
MOTIF 1 2 3
----- ----- ----- -----
2 0.18
3 0.06 0.13
4 0.57 0.33 0.11
No overly similar pairs (correlation > 0.60) found.
Random model letter frequencies (from non-redundant database):
A 0.281 C 0.222 G 0.229 T 0.267
********************************************************************* 这里可以看到我们的数据集是rice9311_sample.fa，有7条序列，共4200nt 的残基。

我们看到用MEME的输出结果rice_meme_txtBriefFormat.out作为MAST 的输入文件，并列出了有4个motif用于搜索。

关于motif的预处理步骤，其实是计算几个motif之间的相关性。

最后还给出随机背景的概率分布。

(3) 高分序列。

这部分结果返回一个按E值由小到大序列排列，E值越小，说明匹配的显著性越高。

********************************************************************* SECTION I: HIGH-SCORING SEQUENCES
********************************************************************* SEQUENCE NAME DESCRIPTION E-VALUE LENGTH
------------- ----------- -------- ------ ------
BGK03149 3.8e-11 600
BGK02457 1.1e-08 600
BGK04166 9.2e-07 600
BGK02838 9.4e-07 600
BGK03537 6.5e-06 600
BGK04759 4.4e-05 600
BGK01655 0.00049 600
********************************************************************* （4）motif示意图。

和MEME类似，MAST的motif示意图标注了所有motif 在每条序列上匹配的位置。

********************************************************************* SECTION II: MOTIF DIAGRAMS
********************************************************************* SEQUENCE NAME E-VALUE MOTIF DIAGRAM
------------- -------- -------------
BGK03149 3.8e-11 99_[+2]_[+1]_361_[+3]_43_[+1]_41
BGK02457 1.1e-08
79_[+2]_51_[+3]_193_[+3]_133_[-1]_22_[-1]_8_[+1]_7_[-3]_8
BGK04166 9.2e-07 67_[+1]_202_[+1]_63_[+3]_100_[+2]_112
BGK02838 9.4e-07 198_[+2]_101_[+1]_178_[-1]_82
BGK03537 6.5e-06 9_[+2]_77_[+1]_425_[+3]_46
BGK04759 4.4e-05 262_[+3]_32_[-1]_187_[+1]_78
BGK01655 0.00049 65_[-1]_199_[+1]_16_[+3]_279
********************************************************************* 以BGK03149为例，99_[+2]_[+1]_361_[+3]_43_[+1]_41，这里的数字99表示以5’到3'的99nt距离，[+2][+1][+3]表示motif 1、2和3,+号表示正链匹配、-号表示反链匹配。

(5)序列注释。

MAST不光是给出各motif的位置信息，还能对每个序列每个motif匹配的位置给出详细的注释信息。

从下面的结果我们可以看到，MAST不光给出了序列P值和E值，还给出每个motif匹配单独的位置P值。

********************************************************************* SECTION III: ANNOTATED SEQUENCES
********************************************************************* BGK03149
LENGTH = 600 COMBINED P-VALUE = 5.37e-12 E-VALUE = 3.8e-11
DIAGRAM: 99_[+2]_[+1]_361_[+3]_43_[+1]_41
[+2] [+1]
2.3e-09
3.8e-05
TGGTTTTTAGGTGGAAGAGAGAGGGAAG
+++++++++++++++++ ++++++++ +
76 AAATCAGGTTTTATATGAGGAAAATGGTTTTTAGGAGGAAGTGAGAGAGAGGATGGAGGGA TATACTATAGAGAT
[+3]
5.1e-08
CTCGCTTCCCCCGCA
++++ ++++++++++
451 ACCATTCACTTGTAGTGTATATATACTGCCACCATTCACTCGTTCCCCCCGCATCGCATCCTCTCC CCACACAAT
[+1]
2.9e-08
AGAGAGAGGGAAG
+++++++++++++
526 CCCCACCCGGCCTAGCGAGCGAGAGAGAGGGAAGGGTTGCGACTTGCGACTCGCGAGCGG CGAGATGGCGAAGGA
（注意：这里的526表示这一行从序列的526碱基开始是，而[+1]表示motif 在正链上匹配，2.9e-8是位置P值AGAGAGAGGGAAG是最有可能的一个匹配，。