2.1 DNA搜索_算法精粹：经典计算机科学问题的Python实现-QQ阅读男生玄幻网

上QQ阅读APP看书，第一时间看更新

2.1 DNA搜索

在计算机软件中，基因通常会表示为字符A、C、G和T的序列。每个字母代表一种核苷酸（nucleotide），3个核苷酸的组合称作密码子（codon）。如图2-1所示，密码子的编码决定了氨基酸的种类，多个氨基酸一起形成了蛋白质（protein）。生物信息学软件的一个经典任务就是在基因中找到某个特定的密码子。

图2-1 核苷酸由A、C、G和T之一表示。密码子由3个核苷酸组成，基因由多个密码子组成

2.1.1 DNA的存储方案

核苷酸可以表示为包含4种状态的简单类型IntEnum，如代码清单2-1所示。

代码清单2-1 dna_search.py

from enum import IntEnum
from typing import Tuple, List

Nucleotide: IntEnum = IntEnum('Nucleotide', ('A', 'C', 'G', 'T'))

Nucleotide的类型是IntEnum，而不仅仅是Enum，因为IntEnum“免费”提供了比较运算符（<、>=等）。为了使要实现的搜索算法能完成这些操作，需要数据类型支持这些运算符。从typing包中导入Tuple和List，是为了获得类型提示提供的支持。

如代码清单2-2所示，Codon可以定义为包含3个Nucleotide的元组，Gene可以定义为Codon的列表。

代码清单2-2 dna_search.py（续）

Codon = Tuple[Nucleotide, Nucleotide, Nucleotide]# type alias for codons
Gene = List[Codon]# type alias for genes

注意尽管稍后需要对Codon进行相互比较，但是此处并不需要为Codon定义显式地实现了“<”操作符的自定义类，这是因为只要组成元组的元素类型是可比较的，Python就内置支持元组的比较操作。

互联网上的基因数据通常都是以文件格式提供的，其中包含了代表基因序列中所有核苷酸的超长字符串。下面将为某个虚构的基因定义这样一个字符串，并将其命名为gene_str，具体代码如代码清单2-3所示。

代码清单2-3 dna_search.py（续）

gene_str: str = "ACGTGGCTCTCTAACGTACGTACGTACGGGGTTTATATATACCCTAGGACTCCCTTT"

这里还需要一个实用函数把str转换为Gene。具体代码如代码清单2-4所示。

代码清单2-4 dna_search.py（续）

def string_to_gene(s: str) -> Gene:
    gene: Gene = []
for i in range(0, len(s), 3):
if (i + 2) >= len(s):# don't run off end!
   return gene
#  initialize codon out of three nucleotides
        codon: Codon = (Nucleotide[s[i]], Nucleotide[s[i + 1]], Nucleotide[s[i + 2]])
        gene.append(codon)  # add codon to gene
return gene

string_to_gene()遍历给定的str，把每3个字符转换为Codon，并将其追加到新建的Gene的末尾。如果该函数发现从正在读取的s的当前位置开始不够放下两个Nucleotide的位置（参见循环中的if语句），就说明已经到达不完整基因的末尾了，于是最后一个或两个核苷酸就会被跳过。

string_to_gene()可用于把str类型的gene_str转换为Gene。具体代码如代码清单2-5所示。

代码清单2-5 dna_search.py（续）

my_gene: Gene = string_to_gene(gene_str)

2.1.2 线性搜索

基因需要执行的一项基本操作就是搜索指定的密码子，目标就是简单地查找该密码子是否存在于基因中。

线性搜索（linear search）算法将按照原始数据结构的顺序遍历搜索空间（search space）中的每个元素，直到找到搜索内容或到达数据结构的末尾。其实对搜索而言，线性搜索是最简单、最自然、最显而易见的方式。在最坏的情况下，线性搜索将需要遍历数据结构中的每个元素，因此它的复杂度为O(n)，其中n是数据结构中元素的数量，如图2-2所示。

图2-2 在最坏情况下，线性搜索将需要遍历数组的每个元素

线性搜索函数的定义非常简单。它只需要遍历数据结构中的每个元素，并检查每个元素是否与所查找的数据相等。代码清单2-6所示的代码就对Gene和Codon定义了这样一个函数，然后尝试对my_gene和名为acg和gat的Codon对象调用这个函数。

代码清单2-6 dna_search.py（续）

def linear_contains(gene: Gene, key_codon: Codon) -> bool:
for codon in gene:
if codon == key_codon:
    return True
return False
acg: Codon = (Nucleotide.A, Nucleotide.C, Nucleotide.G)
gat: Codon = (Nucleotide.G, Nucleotide.A, Nucleotide.T)
print(linear_contains(my_gene, acg))# True
print(linear_contains(my_gene, gat))# False

注意上述函数仅供演示。Python内置的序列类型（list、tuple、range）都已实现了__contains__()方法，这样就能简单地用in操作符在其中搜索某个指定的数据项。实际上，in运算符可以与任何已实现__contains__()方法的类型一起使用。例如，执行print(acg in my_gene)语句即可在my_gene中搜索acg并打印出结果。

2.1.3 二分搜索

有一种搜索方法比查看每个元素速度快，但需要提前了解数据结构的顺序。如果我们知道某数据结构已经排过序，并且可以通过数据项的索引直接访问其中的每一个数据项，就可以执行二分搜索（binary search）。根据这一标准，已排序的Python List是二分搜索的理想对象。

二分搜索的原理如下：查看一定范围内有序元素的中间位置的元素，将其与所查找的元素进行比较，根据比较结果将搜索范围缩小一半，然后重复上述过程。下面看一个具体的例子。

假定有一个按字母顺序排列的单词List，类似于["cat", "dog", "kangaroo", "llama", "rabbit", "rat", "zebra"]，要查找的单词是“rat”。

（1）可以确定在这7个单词的列表中，中间位置的元素为“llama”。

（2）可以确定按字母顺序“rat”将排在“llama”之后，因此它一定位于“llama”之后的一半（近似）列表中。（如果在本步中已经找到“rat”，就应该返回它的位置；如果发现所查找的单词排在当前的中间单词之前，就可以确信它位于“llama”之前的一半列表中。）

（3）可以对“rat”有可能存在其中的半个列表再次执行第1步和第2步。这半个列表事实上就成了新的目标列表。这些步骤会持续执行下去，直至找到“rat”或者当前搜索范围中不再包含待查找的元素（意味着该单词列表中不存在“rat”）。

图2-3展示了二分搜索的过程。请注意，与线性搜索不同，它不需要搜索每个元素。

图2-3 最坏情况下，二分搜索仅会遍历lg(n)个列表元素

二分搜索将搜索空间不停地减半，因此它的最坏情况运行时间为O(lg n)。但是这里还有个排序问题。与线性搜索不同，二分搜索需要对有序的数据结构才能进行搜索，而排序是需要时间的。实际上，最好的排序算法也需要O(n lg n)的时间才能完成。如果我们只打算运行一次搜索，并且原数据结构未经排序，那么可能进行线性搜索就比较合理。但如果要进行多次搜索，那么用于排序所付出的时间成本就是值得的，获得的收益来自每次搜索大幅减少的时间成本。

为基因和密码子编写二分搜索函数，与为其他类型的数据编写搜索函数没什么区别，因为同为Codon类型的数据可以相互比较，而Gene类型只不过是一个List。具体代码如代码清单2-7所示。

代码清单2-7 dna_search.py（续）

def binary_contains(gene: Gene, key_codon: Codon) -> bool:
    low: int = 0
    high: int = len(gene) - 1
while low <= high:# while there is still a search space
        mid: int = (low + high) // 2
if gene[mid] < key_codon:
            low = mid + 1
elif gene[mid] > key_codon:
            high = mid - 1
else:
    return True
return False

下面就逐行过一遍这个函数。

low: int = 0
high: int = len(gene) - 1

首先看一下包含整个列表（gene）的范围。

while low <= high:

只要还有可供搜索的列表范围，搜索就会持续下去。当low大于high时，意味着列表中不再包含需要查看的槽位（slot）了。

mid: int = (low + high) // 2

我们用整除法和在小学就学过的简单均值公式计算中间位置mid。

if gene[mid] < key_codon:
    low = mid + 1

如果要查找的元素大于当前搜索范围的中间位置元素，我们就修改下一次循环迭代时要搜索的范围，方法是将low移到当前中间位置元素后面的位置。下面是把下一次迭代的搜索范围减半的代码。

elif gene[mid] > key_codon:
    high = mid - 1

类似地，如果要查找的元素小于中间位置元素之前，就将当前搜索范围反向减半。

else:
return True

如果当前查找的元素既不小于也不大于中间位置元素，就表示它就是要找的元素！当然，如果循环迭代完毕，就会返回False，以表明没有找到，这里就不重现代码了。

下面可以尝试用同一份基因数据和密码子运行函数了，但必须记得先进行排序。具体代码如代码清单2-8所示。

代码清单2-8 dna_search.py（续）

my_sorted_gene: Gene = sorted(my_gene)
print(binary_contains(my_sorted_gene, acg))# True
print(binary_contains(my_sorted_gene, gat))# False

提示可以用Python标准库中的bisect模块构建高性能的二分搜索方案。

2.1.4 通用示例

函数linear_contains()和binary_contains()可以广泛应用于几乎所有Python序列。以下通用版本与之前的版本几乎完全相同，只不过有一些名称和类型提示信息有一些变化。具体代码如代码清单2-9所示。

注意代码清单2-9中有很多导入的类型。本章后续有很多通用搜索算法都将复用generic_search.py文件，这样就可以避免导入的麻烦。

注意在往下继续之前，需要用pip install typing_extensions或pip3 install typing_extensions安装typing_extensions模块，具体命令取决于Python解释器的配置方式。这里需要通过该模块获得Protocol类型，Python后续版本的标准库中将会包含这个类型（PEP 544已有明确说明）。因此在Python的后续版本中，应该不需要再导入typing_extensions模块了，并且会用from typing import Protocol替换from typing_extensions import Protocol。

代码清单2-9 generic_search.py

from __future__ import annotations
from typing import TypeVar, Iterable, Sequence, Generic, List, Callable, Set, Deque, 
    Dict, Any, Optional
from typing_extensions import Protocol
from heapq import heappush, heappop

T = TypeVar('T')

def linear_contains(iterable: Iterable[T], key: T) -> bool:
for item in iterable:
if item == key:
    return True
return False
C = TypeVar("C", bound="Comparable")

class Comparable(Protocol):
def __eq__(self, other: Any) -> bool:
       ...
def __lt__(self: C, other: C) -> bool:
       ...
def __gt__(self: C, other: C) -> bool:
return (not self < other) and self != other

def __le__(self: C, other: C) -> bool:
return self < other or self == other

def __ge__(self: C, other: C) -> bool:
return not self < other

def binary_contains(sequence: Sequence[C], key: C) -> bool:
    low: int = 0
    high: int = len(sequence) - 1
while low <= high:# while there is still a search space
        mid: int = (low + high) // 2
if sequence[mid] < key:
            low = mid + 1
elif sequence[mid] > key:
            high = mid - 1
else:
    return True
return False
if __name__ == "__main__":
    print(linear_contains([1, 5, 15, 15, 15, 15, 20], 5))  # True
    print(binary_contains(["a", "d", "e", "f", "z"], "f"))  # True
    print(binary_contains(["john", "mark", "ronald", "sarah"], "sheila"))# False

现在可以尝试搜索其他类型的数据了。这些函数几乎可以复用于任何Python集合类型。这正是编写通用代码的威力。上述例子中唯一的遗憾之处就是为了满足Python类型提示的要求，必须以Comparable类的形式实现。Comparable是一种实现比较操作符（<、> =等）的类型。在后续的Python版本中，应该有一种更简洁的方式来为实现这些常见操作符的类型创建类型提示。