序列比对基础版
序列比对几乎是我们广大生物狗每天都会使用的操作，其中最鼎鼎大名的程序自然就是Blast，它能以极快的速度在参考数据库中搜索我们的序列。相信大家对于blast的使用都是得心应手，至于它如何实现的，可能没有人关心。就像我们天天用的电脑，会用就好，哪还用得着操心CPU是怎么造的啊。但是作为一个生信猿来说，我们不禁要问，它是如何实现的，这么小一程序，为何能这样强大呢？了解那些经典的生信程序背后的算法，对于我们学习编程，如何用编程去解决生物方面的问题是有很大帮助的。
当然，我们今天讨论的不是blast。大家不要急躁，那位同学快把刀放下。像blast这么好的东西，当然是留到后面再啃了。我们先从最基本的序列比对方法开始讲起，循序渐进。
所谓序列比对只是生物学领域特定的说法，从更广义的范畴来说，其实就是字符串匹配，这在编程中是一个很普遍的问题。字符串的概念已经介绍过，相信大家也都了解了。比如我们现在有一段参考序列seq，也就是一个字符串：
>>>seq= 'CATCCCCTCATCTAAAGGCACAGTATACTGAATGTAGTCCTGAGGCATATCAACCGCACCTGACT'
另一段搜索序列query：
>>>query = 'AGTCCTGAGG'
我们想知道搜索序列query是否被包含在参考序列seq中（下文直接用seq和query表述），如果存在，则需要确定它的位置。最简单直接的方法自然是顺着seq上的碱基一个一个地进行比较。
对照着下图来看，首先将query置于被seq的第一个字符处，将query上的字符与seq上相应位置的字符依次比较，发现不能匹配，因此往前移动一个字符，再次比较，如此循环，当遇到能够匹配的情况时，记录当前位置，继续移动，直到2号位置。
示意图
我们换一种方式，再把这个搜索过程理一理：
- 将query置于seq的第一个字符，我们可以看示意图上标注1号的短序列；
- query中的字符与seq上相应位置的字符挨个比较，如果所有字符都相同，则记录下query序列当前所在位置，假设我们在1号位置就完全匹配了，那么就将1记录下来。如果不能完全匹配，则不记录；
- 向后移动一个字符，再看示意图，此时我们来到2号位置。返回步骤2，进行比对；
- 重复进行步骤2、3，直到最后。
其实编程只是我们自己意识的延伸，重点仍然是我们自己的所思所想，编程语言是我们同计算机交谈的语言，我们把自己的想法写下来，让计算机读懂，然后替我们执行那些重复无聊的计算。大家也可以在纸上将过程画下来，帮助自己建立程序的结构。现在我们把这个程序的思路已经理清了，主要过程就是在重复2、3两个步骤，那把这两步放在一个for循环里就好了，整个过程水到渠成。
def exact_match(pattern, text):
# 创建一个列表，用于记录搜索到的位置
occur = []
# 计算迭代终点的位置
last_pos = len(text) - len(pattern) + 1
# 依次对每一个碱基进行迭代
for i in range(last_pos):
# 声明一个变量match，并给它赋值‘True'
match = True
# 对搜索序列进行迭代，将其中的每一个字符依次与被搜索序列中相应
# 位置的字符进行比较，一旦发现不匹配的字符，将match的值变为False，
# 退出当前循环
for j in range(len(pattern)):
if text[i + j] != pattern[j]:
match = False
break
# 若所有碱基比较完以后，match的值仍然为True，则完全匹配，将该位置添加到列表中
# 继续搜索
if match:
occur.append(i)
return occur
occur = naive(query, seq)
print(occur)
这里最关键的点就是每个循环开始时必须将match重置为True，否则不会得到结果，除非特殊情况，第一个位置就是我们要找的位置，但这之后若有匹配也不会有输出了。
还有个问题就是为什么要确定迭代终点的位置，直接到最后一个字符就好了，为什么要多这一步呢？大家自己可以试一试，如果直接迭代到最后一个字符会出现什么情况。
今天的代码量不多，但是内容并不少。主要是想告诉大家怎么把具象的思维过程抽象成程序运行的过程，然后写出程序。当然一个程序不只有这一种写法，大家可以根据自己的理解重新编写程序。
这种最朴素的搜索方法已经完成了，那么像这样粗暴的方法速度怎么样呢？我们来分析一下，假设被搜索序列长度为n，搜索序列长度为m，我们在被搜索序列上需要扫过的字符数是（n-m），每扫过一个字符，最多需要进行m次比较，每次比较就是一次基本计算，那么总共需要进行（n-m）*m次比较。假设我们要在人的基因组中搜索一条100bp的序列，大约需要多长时间呢？我们来估算一下，现在普通PC的计算速度大约为20亿次每秒，那么
>>>3e9 * 100 /2e9
150.0
150s!
这个速度估计没有人能接受，跟业界标杆blast比起来就像蜗牛爬一样。我一通操作猛如虎，你让我盯着屏幕看转圈圈？
那有没有办法改进呢？能不能变快呢？
答案必须是：能！
我们将在后续的文章中学习更强大的序列搜索方法。
更多精彩生信文章请关注公众号：易微升
