构建PSSM的步骤

Overview

PSSM算法是生物信息学领域中的一个常用算法，全名“位置特异性打分矩阵(position-specific scoring matrix)”，又称作"位置比重矩阵(position weight matrix)".有关该方法更多的细节，详见维基百科Position weight matrix.本文仅阐述其设计思想，实际项目的例子将在另一篇文章中进行介绍.

1.建模思想

一个PWM包含N行(列)，当模型为DNA时，N=4；当模型为蛋白质时，N=20。维基百科给出的是DNA的核苷酸序列的例子，对于我们的项目来说，不同的也仅仅是矩阵行列数而已：组成DNA的基本核苷酸有A,C,G,T四种，故行列式有四行(列)；组成蛋白质的基本氨基酸有二十种，故行列式有二十行(列)。同时PWM对于每个不同位置都对应一列(行)数据。而选择基本核苷酸(氨基酸)作为行或列，需要针对各自的需要进行调整。

1.1 构建位置频度矩阵(`PFM`)

位置频度矩阵(position frequency matrix)的构建的思路还是比较简单的，这里仍然以维基百科的例子介绍。DNA序列如下表所示：

GAGGTAAAC TCCGTAAGT CAGGTTGGA ACAGTCAGT TAGGTCATT TAGGTACTG ATGGTAACT CAGGTATAC TGTGTGAGT AAGGTAAGT

共有10行9列，则相应的PFM为

$$ M=\begin{matrix} A\\ C\\ G\\ T\end{matrix}\begin{bmatrix} 3 &6 &1 &0 &0 &6 &7 &2 &1 \\ 2 &2 &1 &0 &0 &2 &1 &1 &2 \\ 1 &1 &7 &10 &0 &1 &1 &5 &1 \\ 4 &1 &1 &0 &10 &1 &1 &2 &6 \end{bmatrix} $$

这个矩阵的列数也是9，而每一列的数各自加起来正好是10，也就是DNA序列的行数。这样，构建矩阵的原理就很清晰了：计算每一列中的各核苷酸的数量，然后存入矩阵的相应位置。

1.2 构建位置概率矩阵(`PPM`)

通过PFM求PPM，只需要下面的公式：

$$ M_{k,j}=\frac{1}{N}\sum_{i=1}^N I(X_{i,j}=k) $$

其中，i为行号，j为列号，即：

$$ i\in \begin{pmatrix} 1, &2, &\cdots, &N \end{pmatrix},j\in \begin{pmatrix} 1, &2, &\cdots, &l \end{pmatrix}. $$

且I是指示函数，即：

$$ \mathbf{\mathit{I}}_(X_{i,j}=k) = \begin{cases} 1 &,\text{if } X_{i,j} = k ; \\ 0 &,\text{if } X_{i,j} \neq k . \end{cases} $$

则1.1中的PFM相应的PPM为：

$$ M=\begin{matrix} A\\ C\\ G\\ T\end{matrix}\begin{bmatrix} 0.3 &0.6 &0.1 &0.0 &0.0 &0.6 &0.7 &0.2 &0.1 \\ 0.2 &0.2 &0.1 &0.0 &0.0 &0.2 &0.1 &0.1 &0.2 \\ 0.1 &0.1 &0.7 &1.0 &0.0 &0.1 &0.1 &0.5 &0.1 \\ 0.4 &0.1 &0.1 &0.0 &1.0 &0.1 &0.1 &0.2 &0.6 \end{bmatrix} $$

在这里有一点需要特别注意：

对于每一个匹配成功的核苷酸，我们计分为1，未匹配则记为0，这只是一个简单的思想。而在实际情况下，我们用BLAST或者PSI-BLAST等程序求序列的PSSM时，相应的打分规则就要复杂许多。这时候就需要根据自身的需要选择相应的各种参数，如:gap,λ。

我们假定这些矩阵在行列式的每个位置上都是统计独立的。非统计独立模型请参考Chris的一篇文章http://www.nohup.cc/article/111/，这里就不赘述了。那么对于DNA序列S=GAGGTAAAC来说，它出现的概率就是

$$ p\left ( S\mid M \right )= 0.1\times 0.6\times0.7\times1.0\times1.0\times0.6\times0.7\times0.2\times0.2=0.0007056 $$

1.3 构建位置比重矩阵(`PWM`)

这里引入一个参数b，b=1/k,其中，当序列为蛋白质序列时，k=20；序列为DNA时，k=4.那么对于相同的位置上的PWM和PPM的矩阵元素，其关系为：

$$ M_{PWM}=ln\left ( \frac{M_{PPM}}{b} \right ) $$

那么1.2中的PPM相对应的PWM为：

$$ M=\begin{matrix} A\\ C\\ G\\ T\end{matrix}\begin{bmatrix} 0.18 &0.87 &-0.91 &-\infty &-\infty &0.87 &1.02 &-0.22 &-0.91 \\ -0.22 &-0.22 &-0.91 &-\infty &-\infty &-0.22 &-0.91 &-0.91 &-0.22 \\ -0.91 &-0.91 &1.02 &1.38 &-\infty &-0.91 &-0.91 &0.69 &-0.91 \\ 0.47 &-0.91 &-0.91 &-\infty &1.38 &-0.91 &-0.91 &-0.22 &0.87 \end{bmatrix} $$

至此，我们已经介绍完构建PSSM的基本思想。