正则表达式在生物学中的应用（重排版）

作者: Liquidated | 来源:发表于2021-03-09 17:27 被阅读0次

最近在python for biologists看到一篇对正则表达式不错的总结，在此简要翻译。

模式在生物学中的重要性

大部分时候，编程研究生物学问题的过程是搜寻字符串中的模式（searching for patterns in strings）。最典型的例子是生物序列数据分析，其中DNA、RNA和蛋白质序列恰好是字符串。许多生物序列挖掘问题，均可以从模式的角度来考虑，比如，

蛋白质结构域（protein domains）
DNA转录因子结合基序（DNA TF binding motifs）
限制酶酶切位点
简并PCR引物位点（degenerate PCR primer sites）
单核苷酸重复（runs of mononucleotides）

然而，不只是序列数据会有特殊的模式（interesting patterns）。正如我们讨论文件时提到，其他数据大部分也是以文本文件中的字符串形式存在，比如，

读长定位（read mapping locations）
地理坐标（geographical sample coordinates）
分类名（taxonomic names）
基因名
基因编录号（gene accession numbers）
BLAST搜索结果

在之前的博客中，我们已经了解了一些关于字符串中模式识别（pattern recognition）的编程问题，比如，计算蛋白质序列中每个氨基酸的数目、找出DNA序列中的限制酶酶切位点、检查基因名各部分，并将它们与单独的字符（串）匹配（match them against individual characters）。

这些问题的共同主题是搜寻一个固定的模式。但是，还有许多问题需要更灵活的模式。比如，

给定一条DNA序列，求poly-A尾巴的长度
给定一个基因编录名（gene accession name），提取第三个字符和下划线之间的字符
给定一条蛋白质序列，确定它是否含有高度冗余的蛋白质结构域基序（protein domain motif）

此类问题在许多领域中出现过，因此Python中有一套标准的工具处理此类问题：正则表达式。一般的编程书中可能不会提及正则表达式，但是它们在生物学中非常实用。

尽管Python自带处理正则表达式的工具，却不能在编写程序时直接调用。为此，我们必须先讨论一下模块。

Python中的模块

到现在为止，我们提到的函数和数据类型都是基本数据类型，在每个程序中都需要用到，比如处理字符串和数字、处理读写文件和操控数据列表的工具。因此，在开始编写python程序时即可直接调用。比如，若要打开文件，仅需简单地使用open()函数声明。

然而，Python中还有一类工具，功能更专一。比如，正则表达式，此外还有许多实用的专业工具，一般情况下用不到。比如做高级数学计算，从网上下载数据，运行外部程序，修改日期等。每一种专业工具的集合——实际上是一系列特殊函数和数据类型——称为模块。

出于效率的考虑，Python并不会自动提供这些模块。相反，需要使用import显式加载所需的模块。比如，处理正则表达式的模块是re，所以调用时需要在程序最开始输入：

import re

要使用模块中的某个工具，需要在前面加上模块名。比如，使用正则中的search()函数，需要写：

re.search(pattern, $\color{blue}{string}$ )

而不是：

search(pattern, $\color{blue}{string}$ )

如果忘记导入模块或把模块名加入到函数前，会出现NameError。

本篇余下的代码均需要声明import re。为求简洁，下文不包含该声明，当你运行这些例子时，需要在代码开头加上。

裸字符（RAW STRINGS）

使用正则表达式，需要输入许多特殊字符。比如，\n表示换行，\t表示插入tab。

不幸的是，特殊字符数量有限，所以在正则表达式中，表示特殊含义的字符会与本身具有特殊含义的字符冲突（Unfortunately, there are a limited number of special characters to go round, so some of the characters that have a special meaning in regular expressions clash with the characters that already have a special meaning）。Python解决这个问题的办法是对字符串设定特殊规则：如果在左双引号前加上r，那么字符串中的任何特殊字符都会被忽略：

print(r"\t\n")

r是raw的缩写，在Python中的含义是忽略特殊字符。可以看到r在双引号外面——不是字符串的一部分。以上代码的输出是：

\t\n

并不会出现tab和换行。在这篇博文中，你会看到所有的正则表达式中都有裸字符标记——即便有时候不那么必要——但应该养成这种好习惯。

搜寻字符串中的模式

我们将从最简单的正则表达式工具开始，re.search()是真/假函数，可以确定一个模式是否出现在字符串中。需要两个参数，都是字符串。第一个参数是待搜寻的模式，第二个参数是在其中搜寻模式的字符串。比如，我们检测一条DNA序列是否包含EcoRI限制位点：

dna = "ATCGCGAATTAC"
if re.search(r"GAATTC",dna):
    print("restriction site found!")

注意到我们在模式字符串前面使用裸字符标记，即使（其中不含特殊字符）并不是一定需要加上。

交替（ALTERNATION）

现在，我们已经知道如何使用re.search()，下面我们看看更有趣的。这次，我们要查看是否有AvaII识别位点，分别是：GGACC和GGTCC。一种方法是用or做条件判断：

dna = "ATCGCGAATTCAC"
if re.search(r"GGACC", dna) or re.search(r"GGTCC", dna):
    print("restriction site found!")

但是，有一种更实用的方法——用正则表达式来表示（捕捉）AvaII位点的变异——交替（alteration）。为了呈现出多个不同的变体（alternatives），我们将这些变体放在括号里，用管道符隔开。在AvaII的例子中，在第三个碱基上存在两个变体——可以是A或T——因此该模式可写做：

GA(A|T)CC

写成裸字符模式并将其放在re.search()后如下：

dna = "ATCGCGAATTCAC"
if re.search(r"GA(A|T)CC", dna):
    print("restriction site found!")

注意到我们只用了一个模式就（在一个字符串中）表示了序列的所有变化（a single pattern which captures all the variation in the sequence in one string）。

字符族（CHARACTER GROUPS）

BisI限制酶酶切位点包含更多基序（motifs）——模式是GCNGC，其中N代表任意碱基。我们可以用同样的办法来表示这个模式：

GA(A|T|G|C)GC

然而，正则表达式的另一特性可以使模式更为简洁。一对方括号包含一串字符可以表示（匹配）其中任一字符。因此，模式GC[ATGC]GC匹配GCAGC、GCTCC、GCGGC和GCCGC。以下是一个使用字符族检查是否存在BisI限制位点程序：

dna = "ATCGCAATTCAC"
if re.search(r"GC[ATGC]GC", dna):
    print("restriction site found!")

总之，在编程研究生物学问题时，交替（alteration）和字符族（character groups）对捕捉这类变体十分有效。在我们继续之前，还有两个处理特定常见情况的技巧值得了解。

如果需要模式中某一字符匹配输入中的任一字符，可以使用句号或点号。比如，模式GC.GC会匹配BisI酶切位点的四种可能。然而，句号也能匹配那些不是碱基的字符，甚至非字母。因此，它也可以匹配GCFGC、GC&GC和GC9GC，这可能并非我们所想，使用时需要小心。

有时候更简单的方法是，不罗列所有可接受的字符，而是指定不匹配哪些字符。如下所示，将脱字号^放在字符族前

[^XYZ]

会否定该字符族，而去匹配那些不在该字符族中的字符。上述例子会匹配任一字符，除了X、Y和Z。

数量词（QUANTIFIERS）

以上所述的正则表达式特性能够描述模式中单个字符的变化。另一类特性，数量词，能够描述模式中某一部分重复多次的变化。

单字符后加上问号（?）表示该字符是可选的——即匹配0或1次。因此，在模式GAT?C中，T是可选的，该模式将匹配GATC或GAC。如果我们将问号用到多个字符时，可以用括号将字符组成字符族。比如，模式GGG(AAA)?TTT，其中3个A是可选的，因此该模式会匹配GGGAAATTT或GGGTTT。

单字符或字符族后加上加号（+）表示，它（们）必须存在且重复任意次数——换言之，它会匹配一次或多次。比如，模式GGGA+TTT会匹配三个G，接着一个或多个A，接着3个T。因此，它会匹配GGGATTT、GGGAAATT、GGGAAATT等，但不会匹配GGGTTT。

单字符或字符族后加上星号（*）表示，它（们）是可选的，但是也可以重复。换言之，它会匹配0或多次。比如，模式GGGA*TTT会匹配3个G，接着0或多个A，接着3个T。因此，它会匹配GGGTTT,GGGATTT,GGGAATTT等。这是最灵活的量词。

如果想匹配特定数量的重复，可以使用花括号（{}）。单字符或字符族后加上包含数字的花括号会匹配正好那个数字次数的重复。例如，模式GA{5}T会匹配GAAAAAT，不匹配GAAAAT和GAAAAAAT。单字符或字符族后加上内含一对被逗号分隔的数字的花括号，可以指定特定范围次数的重复。例如，GA{2,4}T表示，G接着2至4个A，然后T。因此，它会匹配GAAT，GAAAT和GAAAAT，但不匹配GAT或GAAAAAT。

正如子串一样，我们可以取消下限或上限。A{3,}会匹配3或多个A，G{,7}会匹配到最多7个G。

位置

我们将了解最后一组正则表达式工具，它们不代表字符，而是代表字符串中的位置。脱字号^匹配字符串的开始，美元符$匹配字符串的末尾。模式^AAA匹配AAATTT却不匹配GGGAAATTT。模式GGG$匹配AAAGGG但不匹配AAAGGGCCC。

组合

正则表达式真正的威力来自组合使用这些工具（combining these tools）。我们可以用数量词和交替和字符族来指定非常灵活的模式。例如，下面是一个复杂的模式，鉴定全长真核信使RNA序列：

^AUG{AUGC}{30, 1000}A{5,10}$

从左往右看，它会匹配：

序列起点的AUG起始密码子
然后30~1000bp的序列，可以是A、U、G、C
末尾5~10bp的poly-A尾巴

如你所见，正则表达式刚开始读起来十分棘手，直到你熟悉它们。然而，它值得你投入一些时间学习使用，因为这种概念在许多工具中都会用到。在Python中掌握的正则表达式技巧可以迁移到其他编程语言中、命令行工具和文本编辑器中。

以上我们讨论的特性是在生物学中最实用的。然而，Python中还有许多正则表达式的特性可以使用。如果你想掌握正则表达式，需要了解贪婪vs.最小量词（greedy vs. minimal quantifiers）、后向引用（back-references）、前向和后向断言（lookahead and lookbehind assertions）和内建的字符类（character classes）。

在继续了解正则表达式更复杂的使用之前，需要注意有一个类似re.search()的方法，叫做re.match()。不同之处在于，re.search()会鉴定出字符串中任意位置出现的模式，而re.match()只会鉴定出匹配整个字符串的模式。大多数时候，我们需要的是前一种功能。

其他用到模式的方法（函数）

之前所有的例子中，我们使用re.search()当作if语句中的条件，决定某字符串是否包含某个模式。然而，还可以用正则表达式实现更多有趣的功能。

提取匹配的部分

通常，我们不仅想知道某个模式是否存在，还想了解哪部分字符串匹配上。为此需要使用re.search()的结果，然后group()其结果对象。

在介绍re.search()函数时，我曾提到它是真/假函数。这并不完全正确——如果它找到一个匹配，不返回True，而是返回一个对象，这个对象在条件判断时会被定为真。

它返回的值实际上是一个匹配对象（match object），一种新的数据类型，我们至今尚未遇到。类似文件对象，匹配对象不代表数或字符串这类简单的事物。相反，它代表正则表达式搜寻的结果。类似文件对象，匹配对象有许多实用的方法可以获得其中的数据。

其中一个方法是group()。如果对正则表达式搜索结果调用这个方法，可以获得输入字符串中匹配上该模式的部分。举个例子：想象有一条DNA序列，需要确定是否包含模糊碱基，如哪些不是A、T、G和C的碱基。可以用否定字符族（negated character group）来写不会匹配任一非ATGC碱基的正则表达式。

[^ATGC]

测试如下：

dna = "ATCGCGYAATTCAC"
if re.search(r"[^ATGC]", dna):
  print("ambiguous base found!")

上面的代码告诉我们DNA序列是否含有非ATGC碱基，但是不能告诉我们到底是什么碱基。为此，我们对匹配对象调用group()：

dna = "CGATNCGGAACGATC"
m = re.search(r"[^ATGC]", dna)

#m是一个匹配对象
if m:
    print("ambiguous base found!")
    ambig = m.group()
    print("the base is " + ambig)

这段程序的输出是：

ambiguous base found!
the base is N

不仅告诉序列中包含模糊碱基，而且模糊碱基是N。

提取多个字符族

如果要提取不止一个模式该怎么做呢？假如要匹配学名如Homo sapiens或Drosophia melanogaster。这个模式相对简单：多个字符，紧接着一个空格，再接着多个字符：

.+.+

为了匹配多个字符，先用句号（表示任一字符），紧接着加号（表示重复至少一次但是可能多次）。

现在，我们来提取属名和物种名到不同的变量。在我们需要存储的模式的部分加上括号：

(.+)(.+)

这就捕获了模式的（各）部分。现在可以通过给group()提供一个参数，提取被捕获的片段。group(1)会返回被第一组括号内的模式部分匹配到的字符串，group(2)则返回被第二组匹配到的字符串……

scientific_name = "Homo sapiens"

m = re.search("(.+) (.+)", scientific_name)

if m:
    genus = m.group(1)
    species = m.group(2)
    print("genus is " + genus + ", species is " + species)

输出展示了同一模式匹配到的两块如何存在不同的变量中。注意到空格不在两个变量中：

genus is Homo, species is sapiens

如果观察仔细，你会意识到正则表达式中括号有三个不同的角色：

在交替中围住变体（surrounding the alternatives in alternation）
框定模式中的子串，同数量词一块使用（grouping parts of a pattern for use with a quantifier）
定义了匹配后需要被提取的部分

获得匹配位置

除了包含匹配内容，匹配对象还包含匹配位置的信息。start()和end()可以获得字符串中模式的起始和终止位置。回顾一下模糊碱基的例子，找出模糊碱基的位置：

dna = "CGATNCGGAACGATC"
m = re.search(r"[^ATGC]", dna)

if m:
    print("ambiguous base found!")
    print("at position " + str(m.start()))

注意，计数从0开始，因此匹配起始于第五个碱基的起始位置是4：

ambiguous base found!

at position 4

多个匹配

上述例子有一个明显的限制，它只能找到一个模糊碱基，因为re.search()只能找到一个匹配。为了处理多个匹配，需要使用re.finditer()，会返回一个匹配对象列表，可以继续用循环处理：

dna = "CGCTCNTAGATGCGCRATGACTGCAYTGC"

matches = re.finditer(r"[^ATGC]", dna)
for m in matches:
    base = m.group()
    pos  = m.start()
    print(base + " found at position " + str(pos))

从输出可以发现，发现三个匹配结果：

N found at position 5
R found at position 15
Y found at position 25

将多个匹配结果转换为字符串

一个常见的场景是，要获得匹配上给定模式的字符串的各部分的列表。以下是匹配A和T连续，长度超过5个碱基的正则表达式模式：

[AT]{6, }

这是一条包含该模式匹配（黑体显示）的DNA序列：

ACTGCATTATATCGTACGAAATTATACGCGCG

可以用re.finder()和group()从中提取匹配模式的字符块：

dna = "CTGCATTATATCGTACGAAATTATACGCGCG"

matches = re.finditer(r"[AT]{6,}", dna)

result = []
for m in matches:
    result.append(m.group())

print(result)

不过这类问题比较常见，所以有特定的方法处理，re.findall()。正如其他方法一样，re.findall()以模式和字符串为参数，但不返回匹配对象列表，而是返回字符串列表。可以重写上面的代码如下：

dna = "CTGCATTATATCGTACGAAATTATACGCGCG"
result = re.findall(r"[AT]{6,}", dna)
print(result)

用正则表达式分割字符串

极个别时候，用正则表达式作为分隔符将一个字符串分割也十分使用。一般的字符串split()不适用，但是re模块有一个自己的split()函数，将正则表达式作为参数。第一个参数是模式，然后是被分割的字符串。

想象一下，有一条保守的DNA序列包含模糊密码，需要提取出所有连续的非模糊碱基（all runs of contiguous unambiguous bases）。只要有不是ATGC的碱基就将DNA分割：

dna = "ACTNGCATRGCTACGTYACGATSCGAWTCG"
runs = re.split(r"[^ATGC]", dna)
print(runs)

输出说明该函数的功能——返回值是字符串列表：

['ACT', 'GCAT', 'GCTACGT', 'ACGAT', 'CGA', 'TCG']

注意到字符串中匹配上模式的部分不在输出中（就像分隔符不在Python的split()输出中）。

总结

在这个部分，我们学习了正则表达式，以及相关函数和方法。从简要介绍两个概念开始，模块和裸字符。虽然它们不是正则表达式的一部分，但是使用正则时很有必要。然后简述了正则表达式模式的特性，迅速了解了可以用它来处理的问题。正如正则表达式可以从简单到复杂，使用上也是如此。可以用它们解决简单的任务——如确定某条序列是否含有特定的基序——或者更复杂的，如使用复杂的模式辨认信使RNA序列。

实际使用时，需要认识到，对任一给定的模式，正则表达式可以用多种方式描述。在博客的开头，我们用到了模式

GG(A|T)CC

来描述AvaII限制酶识别位点，但是同一模式也可以写作：

GG[AT]CC
(GGACC|GGTCC)
(GGA|GGT)CC
G{2}[AT]C{2}

在遇到其他存在多种表达方式时，尽量以最容易读为准则编写正则表达式。

PS:原博文末尾提供了小练习，需要的可以自行前往（正则表达式）。