通过一个motif calling的例子来认识深度神经网络中的卷积运算有什么具体的意义。
背景
motif是什么？
1.motif是指DNA/RNA序列上能够与protein特异性结合的一段序列，通常该序列在7-20bp；
2.对于一个protein(transciptional factor,tf /RNA binding protein,RBP)来说，它结合在不同基因序列上的motif是不完全相同的。因此，我们没法用一条序列来表示一个Motif。那怎么表示一个protein的结合motif呢？常用的有以下3种方法：
a.基于motif的consensus序列来表示
b.基于seq-logo图表示(最常用)
c.基于位置频率矩阵(position frequency matrix,PFM)/位置权重矩阵(position weight matrix,PWM)表示 其中，PWM是PFM的改进版本，它消除了了序列本身固有的背景碱基含量带来的偏好性。
以上三个概念不懂的请查阅维基百科。
motif calling的基本问题
对Motif有了一个基本的认识后，我们来看一看关于motif识别的基本问题。
在生物信息学中，motif识别是一个非常基本而重要的问题。其场景可以简化为：
1.给定一个蛋白质名称(序列)以及该蛋白质的PFM（根据以往的实验数据统计得来）
2.给定一条待识别的DNA/RNA序列
问：该蛋白质与该条待识别的序列有哪些可能结合的位点(换句话说，该蛋白质在这个基因上的结合motif是什么)？
这个问题很好解决。
根据上一小节的知识，我们已知了PFM/PWM，它记录了一个蛋白质Motif各位置各碱基的分布情况，因此，一个最简单的方法就是扫描一遍待检测的序列，看看该序列上哪个区域与PWM的相似性得分最高，那这个区域就是该蛋白在该待检测序列上最有可能的结合motif。
如果一切都这么简单，那就好了。问题在哪呢？问题在于往往一个蛋白质的PFM/PWM不好获取(可能缺少了历史数据)，那就无从谈起如何扫描序列获取Motif。
传统算法Call Motif
经典的Motif calling算法主要利用EM/Gibbs sampling来call PWM(Motif)。我们不妨简单看一下：
问题：给定N条能够与某一个protein结合的DNA/RNA序列，求该蛋白的Motif PWM？
EM算法的思路：
1.随机初始化：在这N条序列上随机选取N个位置作为guess motif(每条序列一个)，计算PWM。
2.(E步)根据当前得到的PWM，对这N条序列进行遍历扫描，选择与PWM相似性最高的区域作为该序列新的Motif；
3.(M步)根据新的Motif，求出新的PWM；
4.循环2-3直至收敛。
然而这个过程由于需要迭代求解，会比较慢，因此有人提出利用Deeplearning来Call motif.
利用卷积神经网络做Motif calling
CNN的出现使得自动学习某一个蛋白质的PFM/PWM称为可能。考虑如下问题：
给定一个蛋白和一些已知的能够与该蛋白质结合的DNA/RNA序列；现在有一条全新的未知的DNA/RNA序列；问：
(1)如何训练一个分类器用于预测该未知序列是否能与该蛋白质结合呢？
(2)如果(1)能，我们是否有办法找出该蛋白在该序列上的结合Motif?
关于第(1)问，是一个相对简单的问题。我们只需要搜集该蛋白和一些已知的能够/不能够与该蛋白质结合的DNA/RNA序列，构成正/负样本，然后把序列作为Input,能否与蛋白结合作为Output。就能够完成该模式识别任务。这个过程可以基于传统的机器学习完成，也可以基于深度学习完成。
关于第(2)问，则是完美的展示深度学习卷积操作统计学意义的例子，第一个利用CNN识别Motif位点的算法是DeepBind。以下是该算法的思想：
考虑到传统利用PWM识别Motif的方法是扫描一遍待检测的序列，看看该序列上哪个区域与PWM的相似性得分最高，那这个区域就是该蛋白在该待检测序列上最有可能的结合motif。而卷积操作恰好与这个传统检测方法有异曲同工之妙：
1.卷积核也需要扫描Input矩阵来完成filtering。
2.如果我们把RNA序列编码成一个One-hot矩阵(维度是(LENGTH OF SEQUENCE)*4),以该矩阵作为CNN的输入,把卷积核的size设置为(LENGTH OF MOTIF)*4。
那么卷积核扫描完该Input矩阵后得到的矩阵维度就是(LENGTH OF SEQUENCE - LENGTH OF MOTIF + 1)*1，此时，如果我们把学习到的卷积核参数看作PWM，那么，卷积核扫描完后得到一个维度为(LENGTH OF SEQUENCE - LENGTH OF MOTIF + 1)*1的layer1不正好代表了Input序列上每一个小片段与PWM的相似程度吗（Input的矩阵是独热编码的，所以输出可以看作相似性得分）？
有了这个原理，我们只需要给相似得分划一个阈值，高于阈值的就认为是可能结合的Motif即可。
卷积的统计学意义：对数似然比
上一小节本质上就是DeepBind算法的思想，但好奇的我们一定会问：为什么会这么巧！原因就在于：卷积操作本质上可以理解是一种对数似然比。
我们通过一个小例子来看一看这个事情：
假设我有一段序列，它的序列是ACTG，假设某个蛋白质的结合motif恰好长度为4，想看一看该序列与PWM的相似性？
这个问题可以利用likelihood：
其中P(A_1)表示A碱基出现在PWM第一个位置的概率，P(C_2)表示C碱基出现在PWM第二个位置的概率，以此类推。
那么，对数Likelihood就可以写作：
\[lnL=ln P(A_1) + lnP(C_2|A_1) +ln P(T_3|A_1,C_2) +ln P(G_4| A_1,C_2,T_3)\]当我们把背景碱基(background,bg)的信号除掉，就成了传统的Motif序列识别算法：
\[ln \frac{P(A_1)}{P(A_{bg})} + ln \frac{P(C_2|A_1)}{P(C_{bg})} +ln \frac{P(T_3|A_1,C_2)}{P(T_{bg})} +ln \frac{P(G_4| A_1,C_2,T_3)}{P(G_{bg})}\]现在，我们再回过头来看一看卷积操作得到的相似性得分是什么？ 不妨假设我们学习到的卷积核参数为：
\[\begin{array}{l|l|l|l|l|} \hline \mathrm{A} & 0.1 & 0.2 & 0.1 & 0.4 \\ \hline \mathrm{C} & 0.2 & 0.2 & 0.2 & 0.2 \\ \hline \mathrm{G} & 0.4 & 0.5 & 0.4 & 0.2 \\ \hline \mathrm{T} & 0.3 & 0.1 & 0.2 & 0.2 \\ \hline \end{array}\]由于原始序列ACTG的Input矩阵恰好编码成了一个独热矩阵，因此可以表示为：
每一行对应一种碱基的独热码，依次为ACGT，那么卷积运算后得到的第一层layer(看做相似性得分)就是：
这本质上就和对数似然比之和的操作思想是一回事。通过扫描序列，求和得到一个得分，得分越高，越可能是信号；得分越低，越可能是背景。而对数似然比之和恰恰就是这种相似性得分的思想。因此，卷积可以理解成某种意义上的对数似然比：
给该相似性得分设定一个人为的阈值，低于该阈值的认为是background；高于阈值的，认为是可能的Motif（signal），这正是起到了一个过滤降噪的作用。
从这点意义上，我们发现，原来卷积是有具体生物学和统计学意义的！
参考
本文整理自邓明华老师的生物信息算法设计课程。
上一篇: 深度学习应用系列——DeepCpG：预测单细胞甲基化状态
下一篇: matplotlib绘图逻辑(上)