前言

  Dof2，具有锌指结构的转录因子家族DOF的成员，许多DOF转录因子基因在维管组织中表达，有一些在拟南芥维管发育的早期特异表达，并参与调控维管发育。Dof2 的 Sequence Logo 如下图所示(图片来自JASPAR)，这是我们在文章中或motif分析中经常是使用的图，用来表征序列的一致性和多样性。字母越大，代表在该位置出现该核苷酸或者氨基酸的概率越大，常用Bits或者百分比表示。

想要画出这样的图不是很难，这里先了解几个概念，为分析数据打下基础：

• PFM：Position Frequency Matrix
• PPM：Position Probability Matrix
• PWM：Position-specific Weight Matrix(PSWM) or Position-specific Scoring Matrix(PSSM)

（我就不翻译了……）

PFM&PPM

  假设，我们有这样一些DNA序列（每个9bp）：

GAGGTAAAC
TCCGTAAGT
CAGGTTGGA
ACAGTCAGT
TAGGTCATT
TAGGTACTG
ATGGTAACT
CAGGTATAC
TGTGTGAGT
AAGGTAAGT

那么对应的PFM为(9列)：

  对于上面十个序列，每个位置四种碱基(ACGT)出现的次数，组成了如上所示的4x9的矩阵，比如第一个位置：A 3次，C 2次，G 1次，T 4次。那么对应位置，四种碱基出现的概率(频率)分别为：A 0.3、C 0.2、G 0.1、 T 0.4。计算所有位置的碱基出现概率(频率)，就可以得到PPM矩阵：

  每个位置的概率都是独立与其他位置计算的，所以PPM和PWM都满足模式中位置之间的统计独立性假设，特定位置数值的和为1。因此，每列可以视为一个独立的多项分布，可以通过乘以每个位置的相关概率来计算给定PPM的某一序列的概率。例如，可以计算给定上述PPM (记为M) 的序列S (gaggtaaaac) 的概率：

$$P(S|M)=0.1\times 0.6\times 0.7\times 1.0\times 1.0\times 0.6\times 0.7\times 0.2\times 0.2=0.0007065$$

  在基于小数据集计算PPM时，通常使用伪计数（或拉普拉斯估计），以避免矩阵项的值为0。这相当于将PPM的每一列乘以Dirichlet分布（多项式分布的共轭先验[1]），并允许计算新序列（不属于原始数据集的序列）的概率。在上面的示例中，如果没有伪计数，则无论其他位置如何，任何在第4位置没有G或在第5位置没有T的序列的概率都为0。

PPM 转换为PWM

  PWM中的数值通常为log似然值，使用背景模型(b)将PPM中的元素转换为PWM：

$$M_{kj}=lo{g}_2(M_{kj}/b_k)$$

  M指这个位点的概率，b就是background。我们假设一个最简单的背景模型，即四种核苷酸出现的频率相同，也就是$$b_k=\frac{1}{|k|}$$的值为0.25（如果是蛋白质就是对应的常见二十种氨基酸，每个为0.05）。将上面提到的PPM（记为M，未使用伪计数），使用公式转换，得到PWM：

  比如第一位的A PW值为 log2(0.3/0.25) = 0.2630，可以看到矩阵中出现了无穷，如果使用伪计数，尤其在使用小数据集构建PPM时，就会避免这种情况。当然，在背景模型中不需要每种核苷酸/氨基酸具有相同的值，上面只是一个最简单的例子，在实际的研究当中，AT和CG可能会有不同的含量，CG的值应该会伴随着AT的增加而减少。

  当使用log似然法计算来的PWM时，可以通过将PWM中每个位置的相关值相加（不是相乘，$lo{g}_2(2\times 4)=lo{g}_2(2)+lo{g}_2(4)=3$）来计算序列的分数。序列分数表示序列与随机序列的不同程度。序列作为功能序列和随机序列的概率相同，则得分为0。如果它更可能是功能序列而不是随机序列，则分数大于0；如果它更可能是随机序列而不是功能序列，则分数小于0。

PS： 序列分数也可以在物理框架中解释为该序列的结合能😑

画外

  PPM可以作为Sequence Logo图的输入，R语言中有比较完善的工具，如seqLogo等，就算只有PFM，也可以自己算出PPM来画图，至于怎么画，如果以后我有机会（想起来），就整理一下 Motif 分析和画图