DNA binding motif比对算法

之前介绍了序列比对的一些算法。本节主要讲述motif(有人翻译成结构模式，但本文一律使用基模)的比对算法。

那么什么是基模么？基模是对DNA结合位点的一种描述。它有几种描述方式，一种是共同序列（consensus sequences）一种是位点倾向距阵（Position Specific Frequency Matrices(PSFM)）而对于PSFM，有两种表示方式，一种叫PCM,一种叫PFM，前者是Position count matrices,就是对每个碱基位ACGT出现次数的统计，后者(Position frequency Matrices)是在前者的基础上求出其相对于总数的比率。它们看上去大约是这个样子的：

>PCM motif
A | 367	470	93	1738	0	61	0	0	126	463	414	396
C | 139	12	1688	43	0	109	0	0	471	212	472	411
G | 916	1279	0	0	1781	1601	119	1775	485	426	308	390
T | 148	14	0	0	0	10	1662	6	699	680	584	578
>PFM motif
A              
0.208	0.226	0.214	0.350 
0.132	0.140	0.070	0.655
0.015	0	0.016	0.968
0.812	0.028	0.010	0.148
0.983	0.007	0.009	0
0.007	0.052	0.027	0.912
0.063	0.024	0.456	0.456
0.669	0	0.312	0.018

使用这个距阵，可以生成所谓的motif logo图示。

alignedLogo0.01.02.0bitspositionmotifName1CTAG2CTAG3AC4CA5G6TACG7GT8TG9ACGT10CGAT11GACT12GACT0.01.02.0bitspositionanother motifFormat1AGCT2GACT3AGT4GCTA5CGA6AGCT7CAGT8TGA

对于基模比对，其主要思路还是运用序列比对当中Needleman-Wunsch或者Smith-Waterman算法。但针对每个碱基位置，因为并非单一的字母，而是一个四个或者五个数字的数组，所以会改变之前字母比对的方式，而使用以下五种算法：

Pearson’s correlation coefficient 皮尔逊积矩相关系数
sum of squared distances 距平方和
Kullback–Leibler information content 相对熵
average log-likelihood ratio (ALLR) 平均对数似然率
ALLR with a lower limit of –2 imposed on the score

其实最后一种不应该说是一种算法，只是一种改进，所以我们只略约讲一下前5种方法。

Pearson’s correlation coefficient算法是最容易被想到的一种算法，它将基模中的ACGT分布概率或者个数统计当成两组数，我们必须假设它们拥有相同的分布概率，比如说高斯分布，如果两组数对应位置呈现线性关系的话，那么它们的皮尔逊系数就等于1或者-1.如果它们完全不相关，那么系数就等于0。皮尔逊相关系数的算法为协方差和标准差的商。
$ho_{X,Y}={mathrm{cov}(X,Y) over sigma_X sigma_Y} ={E[(X-mu_X)(Y-mu_Y)] over sigma_Xsigma_Y}$

很明显，这一算法应用在基模中某一碱基位的比对上有一个巨大的缺陷就是它要求准确服从高斯分布。而基模某一碱基位上只能对应四个，或者五个数字，而一个偏向性强的基模位显然不会符合高斯分布。

对于相对距平方和算法很简单，其公式为 $S = 2-sum_{b in {A,C,G,T}}^{} (M_{b}-N_{b})^2$ 。其中M和N是要比较的两碱基位那一列的频率。这一算法是JASPAR用于搜索相类基模的算法。

后面的算法都是从信息论当中来的。它们的基础都是由Claude Shannon提出的信息熵（Information Entropy）：
$H =-sum_{i}^{} p_i log p_i$
理解这个信息熵可能需要一点点解释。我们说DNA中有4个碱基ACGT，那么每个字母所代表的信息量为：
$I_e = -log_2 {1over 4} = 2$ ,
那么对于氨基酸来讲有20种（也有说22或者26种的）ARNDCEQGHILKMFPSTWYV(UO(BZJX))，那么每个字母所代表的信息量为：
$I_e = -log_2 {1over 20} = 4.3$
如果说三个碱基对应一个氨基酸，那其信息量是2*3=6 ==> 4.3，所以三个碱基的信息量足以转换成一个氨基酸的信息量。
那么熵就是整个系统的总的平均信息量：
$H_s = sum_{i=1}^n p_i I_e = -sum_{i=1}^n p_i log_2 p_i$

知道了这一点，我们可以想象一下将它使用到基模上是会代表什么？其实就是基模所提供的总的信息量。

相对熵的算法是一种比较两种分布的算法。我们来看一眼Kullback–Leibler divergence (information gain)的公式：
$D_{mathrm{KL}}(p(X) | q(X)) = sum_{x in X} -p(x) log {q(x)} \, - \, left( -p(x) log {p(x)} ight) = sum_{x in X} p(x) log frac{p(x)}{q(x)}$
我们假设两个分布，一个是不受约束的分布q(X)，一个是受约束的分布p(X)。我们相知道两个分布是不是一样的，我们使用上面的公式，如果两分布完全一样，那么 $log frac{p(x)}{q(x)}$ 就等于0，于是两者的相对熵就是0。

我们再来看平均对数似然率（ALLR）算法，它其实是来源于互信息(Mutual Information)这个概念，它的公式为
$I(X;Y) = mathbb{E}_{X,Y} [SI(x,y)] = sum_{x,y} p(x,y) log frac{p(x,y)}{p(x)\, p(y)}$
它是对称的
$I(X;Y) = I(Y;X) = H(X) + H(Y) - H(X,Y)$
它与相对熵的关系为：
$I(X; Y) = D_{mathrm{KL}}(p(X,Y) | p(X)p(Y))$
而在stromo文章中是这样描述的：
$ALLR = frac{sum_{b=A..T}n_{bj}log{f_{bi}over p_b}+sum_{b=A..T}n_{bi}log{f_{bj}over p_b}}{sum_{b=A..T}n_{bi}+n_{bj}}$
其中i和j是需要比较的两列，每一列都有ACGT的出现次数计数，称为 $n_{bi}$ 或者 $n_{bj}$ ，而对于每一列还可以计算出ACGT出现的频率计数，计为 $f_{bi}$ 或者 $f_{bj}$ 。在他们所提供的MatAlign基模比对工具中就是使用的这一算法来计算两基模的距离的：Dist = ALLR(A:A) + ALLR(B:B) – 2*ALLR(A:B)。这一公式与互信息的公式完全一致。这现次说明一个问题，生物学中的很多算法都不具有原创性，只是将数学工具进行了一个包装而已。

但是基于信息论的算法有一个缺点，那就是如果两个基模的长度不一样的时候，它们包含的信息量的差别并不一定代表两者的差别。比如说如果一个基模完全是另一个基模的一部分，那么它的信息量会小于另一个基模，如果这时需要比较第三个基模，其与长的那个的信息的差别甚至小于短的那个基模时，就会得出第三个基模的与长的基模的距离甚至小于长基模与短基模的距离。这显然是错误的。

其实上面的这个问题，在现有的算法当中都存在。如何克服，还是一个问题。有一种解决办法就是截短。但也不一定就总是有效。

STMAP工具就是一个集合了上述所有算法的基模比对工具。但遗憾的是它提供的距离值的分辨率不足，很多不同的基模都会得出距离为0的情况。