第九章 限制性图谱和正则表达式

在本章中，我将概述Python正则表达式和Python运算符。我们还将研究标准的基本分子生物学技术的编程：发现序列的限制性图谱。限制性消化是“指纹”DNA的原始方法之一；现在可以在计算机上模拟这个。

限制性图及其相关的限制性位点是实验室中的常见计算，由几个软件包提供。它们是克隆实验的重要工具；例如，它们可用于将所需的DNA片段插入克隆载体中。限制图也可用于测序项目，例如鸟枪和定向测序。

1. 正则表达式模块re

我们一直在处理正则表达式，本节填补了一些背景知识，并将本书前面部分正则表达式的分散讨论联系在一起。使用生物数据（如序列）或GenBank、PDB和BLAST文件进行编程时，正则表达式非常有用。

正则表达式是用一个字符串表示和搜索许多字符串的方法。虽然它们并不是完全相同的东西，但将正则表达式看作是一种高度发达的通配符集是很有用的。正则表达式中的特殊字符更恰当地称为元字符。

大多数人都熟悉通配符，这些通配符可以在搜索引擎或扑克游戏中找到。例如，你可以通过键入biolog*找到对biolog开头的每个单词的引用。或者你可能会发现自己拿着五个ace。

在计算机科学中，这些通配符或元字符在实践和理论上都有重要的历史。特别是星号字符在发明它的杰出逻辑学家之后称为Kleene闭包。

我们已经看到许多使用正则表达式来查找DNA或蛋白质序列中的事物的例子。在这里，我将简要介绍正则表达式背后的基本思想，作为一些术语的介绍。

因此，让我们从一个实际的例子开始：使用字符类来搜素DNA。假设你想在DNA库中找到一个小基序长度为6个碱基对：CT后跟C或G或T，然后是ACG。使用正则表达式表示该基序，如下所示：

import re

if re.match("CT[CGT]ACG", dna):
    print "I found the motif!!
"

重复或星号

　　星号（*），也称为Kleene闭包或星号，表示在它之前的字符重复0次或多次。例如，abc*匹配以下任何字符串：ab，abc，abcc，abccc，abcccc等。正则表达式匹配无限数量的字符串。

轮换

　　在Python中，模式（a|b）（读：a或b）匹配字符串a或字符串b。

并列

　　这是一个非常明显的问题。在Python中，字符串ab表示字符a后跟（连接）字符b。

使用元字符括号进行分组很重要。因此，例如，字符串（abc|def）z*x匹配诸如abcx，abczx，abczzx，defx，defzx，defzzzzzx等字符串。这个例子结合了分组、交替、闭包和连接的思想，正则表达式的真正力量可见于这三种基本思想的结合。

2. 限制性图谱和限制性酶

2.1 背景

限制酶是切割DNA上短的特定序列的蛋白质，例如，流行的限制酶EcoRI和HindIII广泛用于实验室。EcoRI在G和A之间切割GAATTC，如果你看一下限制酶EcoRI的反向互补序列，GAATTC，序列相同。这是一个回文的生物学结构，许多限制性位点是回文。

计算限制图是实验室中常见且实用的生物信息学计算。计算限制图以计划实验，找到切割DNA以插入基因的最佳方法，产生位点特异性突变，或用于重组DNA技术的若干其他应用。通过预先计算，实验室科学家可以大大节省必要的反复试验。

现在我们将编写一个程序，在实验室中做一些有用的事情：它将在DNA序列中寻找限制酶，并用限制酶图表报告限制性酶在DNA中出现的确切位置。

2.2 程序设计

在前面几章中，你已经了解了如何在文本中查找正则表达式，现在让我们考虑如何使用这些技术来创建限制图。以下是一些问题：

我在哪里可以找到限制酶数据？

　　限制酶数据可在限制酶数据库（REBASE）上找到，该数据库可在网站http://www.neb.com/rebase/rebase.html上找到。

我如何在正则表达式中表示限制酶？

　　探索该网站，您将看到限制酶以他们自己的语言表示。 我们将尝试将该语言翻译成正则表达式的语言。

如何储存限制酶数据？

　　大约有1000种具有名称和定义的限制酶。 这使得字典类型成为快速键值查找的候选者。 当您编写一个真实的应用程序时，比如说Web，最好创建一个DBM文件来存储信息，当程序需要查找时就可以使用了。 我将在第10章介绍DBM文件; 在这里，我只是展示原理。 我们将在程序中仅保留一些限制酶定义。

我如何接受用户的查询？

　　您可以输入限制酶名称，或者您可以允许用户直接键入正则表达式。 我们做第一个输入方式。 此外，您希望让用户指定要使用的序列，为了简化问题，您只需从样本DNA文件中读取数据。

如何向用户报告限制图？

　　这是一个重要的问题。 最简单的方法是生成一个位置列表，其中包含那里发现的限制酶的名称。 这对于进一步处理很有用，因为它非常简单地呈现信息。但是如果你不想做进一步的处理呢？ 您只想将限制图传达给用户？ 然后，也许提供一个图形显示更有用，可能打印出序列，上面有一行标记酶的存在。你可以使用很多花哨的方式，但现在让我们用简单的方法来输出一个列表。

因此，计划是编写一个将限制酶数据翻译成正则表达式的程序，储存为限制酶名称的键值。该程序从文件中读取DNA序列数据，并提示用户输入限制酶的名称，然后从字典中检索适当的正则表达式，我们将搜索该正则表达式的所有实例及其位置，最后返回找到的位置列表。

2.3 限制酶数据

访问REBASE网站，你将看到限制酶数据有多种格式。最终我们决定从bionet文件中获取信息，该文件具有相当简单地布局。这是标题和该文件中的一些限制酶：

REBASE version 104                                              bionet.104
 
    =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    REBASE, The Restriction Enzyme Database   http://rebase.neb.com
    Copyright (c)  Dr. Richard J. Roberts, 2001.   All rights reserved.
    =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
 
Rich Roberts                                                    Mar 30 2001
 
AaaI (XmaIII)                     C^GGCCG
AacI (BamHI)                      GGATCC
AaeI (BamHI)                      GGATCC
AagI (ClaI)                       AT^CGAT
AaqI (ApaLI)                      GTGCAC
AarI                              CACCTGCNNNN^
AarI                              ^NNNNNNNNGCAGGTG
AatI (StuI)                       AGG^CCT
AatII                             GACGT^C
AauI (Bsp1407I)                   T^GTACA
AbaI (BclI)                       T^GATCA
AbeI (BbvCI)                      CC^TCAGC
AbeI (BbvCI)                      GC^TGAGG
AbrI (XhoI)                       C^TCGAG
AcaI (AsuII)                      TTCGAA
AcaII (BamHI)                     GGATCC
AcaIII (MstI)                     TGCGCA
AcaIV (HaeIII)                    GGCC
AccI                              GT^MKAC
AccII (FnuDII)                    CG^CG
AccIII (BspMII)                   T^CCGGA
Acc16I (MstI)                     TGC^GCA
Acc36I (BspMI)                    ACCTGCNNNN^
Acc36I (BspMI)                    ^NNNNNNNNGCAGGT
Acc38I (EcoRII)                   CCWGG
Acc65I (KpnI)                     G^GTACC
Acc113I (ScaI)                    AGT^ACT
AccB1I (HgiCI)                    G^GYRCC
AccB2I (HaeII)                    RGCGC^Y
AccB7I (PflMI)                    CCANNNN^NTGG
AccBSI (BsrBI)                    CCG^CTC
AccBSI (BsrBI)                    GAG^CGG
AccEBI (BamHI)                    G^GATCC
AceI (TseI)                       G^CWGC
AceII (NheI)                      GCTAG^C
AceIII                            CAGCTCNNNNNNN^
AceIII                            ^NNNNNNNNNNNGAGCTG
AciI                              C^CGC
AciI                              G^CGG
AclI                              AA^CGTT
AclNI (SpeI)                      A^CTAGT
AclWI (BinI)                      GGATCNNNN^

您的第一个任务是阅读此文件并获取每种酶的名称和识别位点（或限制位点）。为了简化现在的问题，只需丢弃带括号的酶名称。如何读取这些数据？

Discard header lines 

For each data line:

    remove parenthesized names, for simplicity's sake

    get and store the name and the recognition site

    Translate the recognition sites to regular expressions
        --but keep the recognition site, for printing out results
}

return the names, recognition sites, and the regular expressions

这是高级详细信息伪代码，所以让我们改进并扩展它。 （请注意，花括号没有正确匹配。这没关系，因为没有伪代码的语法规则;做任何适合你的事情！）这里有一些丢弃标题行的伪代码：

foreach line 

    if /Rich Roberts/ 
            
        break out of the foreach loop

}

这基于文件的格式，你要查找的字符串是数据行开始之前的最后一个文本。（当然，如果文件的格式更改，则可能不再有效。）

现在让我们进一步扩展伪代码，思考如何完成所涉及的任务：

# Discard header lines
# This keeps reading lines, up to a line containing "Rich Roberts"
foreach line 
    if /Rich Roberts/ 
        break out of the foreach loop
}

For each data line:

    # Split the two or three (if there's a parenthesized name) fields
    fields = line.split(" ")

    # Get and store the name and the recognition site
    name = fields[0]

    site = fields[-1]

    # Translate the recognition sites to regular expressions
        --but keep the recognition site, for printing out results
}

return the names, recognition sites, and the regular expressions

让我们来看看你做了什么。首先，您要从字符串中提取名称和识别站点数据。 在Python中分隔单词的最常用方法，特别是如果字符串格式很好，则使用Python字符串内置函数split。
如果每行有两个或三个具有空格并且通过空格彼此分隔，则可以简单调用split将它们放入一个数组中。

fields数组可能有两个或三个元素，具体取决于是否有一个带括号的替代酶。 但是你总是想要第一个和最后一个元素：

name = fields[0]

site = fields[-1]

接下来是将网站翻译为正则表达式的问题。查看识别位点并阅读你在网站上找到的REBASE文档，你知道切割位点由插入符号（^）表示，这无助于使用正则表达式按顺序查找位点，因此你应该将其删除；还要注意识别位点中给出的碱基不仅仅是A，C，G和T碱基，还有其它字母的简并碱基，让我们编写一个函数来替换这些代码的字符类，然后我们将使用正则表达式；当然，REBASE使用它们，因为给定的限制酶可能很好地匹配几个不同的识别位点。

例9-1是一个函数，给定一个字符串，将这些代码转换为字符类。

例子9-1 将IUB歧义碱基转换为正则表达式

# IUB_to_regexp
#
# A subroutine that, given a sequence with IUB ambiguity codes,
# outputs a translation with IUB codes changed to regular expressions
#
# These are the IUB ambiguity codes
# (Eur. J. Biochem. 150: 1-5, 1985):
# R = G or A
# Y = C or T
# M = A or C
# K = G or T
# S = G or C
# W = A or T
# B = not A (C or G or T)
# D = not C (A or G or T)
# H = not G (A or C or T)
# V = not T (A or C or G)
# N = A or C or G or T 

def IUB_to_regexp(iub):

    regular_expression = ''

    iub2character_class = {
    
        'A' : 'A',
        'C' : 'C',
        'G' : 'G',
        'T' : 'T',
        'R' : '[GA]',
        'Y' : '[CT]',
        'M' : '[AC]',
        'K' : '[GT]',
        'S' : '[GC]',
        'W' : '[AT]',
        'B' : '[CGT]',
        'D' : '[AGT]',
        'H' : '[ACT]',
        'V' : '[ACG]',
        'N' : '[ACGT]',
    }

    # Remove the ^ signs from the recognition sites
    iub = iub.replace('^', '')

    # Translate each character in the iub sequence
    for i in range(len(iub)):
        regular_expression += iub2character_class[iub[i]]

    return regular_expression

看起来你已经准备编写一个函数来从REBASE数据文件中获取数据。但是有一个重要的项目你还没有解决：你想要返回的数据究竟是什么？

你计划在原始REBASE文件的每一行返回三个数据项：酶名称，识别位点和正则表达式。这不适合字典，你可以返回一个三个数据项的列表，但可能使查找有点困难。当你学习更高级的Python时，你会发现你可以创建自己的复杂数据结构。

既然你已经学会的split函数，也许你可以得到一个字典，其中键是酶的名称，而值是一个带有识别位点的字符串和由空格分隔的正则表达式。然后，你可以快速查找数据，并使用split提取所需的值。例9-2显示了这种方法。

例子9-2 用于解析REBASE数据文件的函数

# parseREBASE--Parse REBASE bionet file
#
# A subroutine to return a hash where
#    key   = restriction enzyme name
#    value = whitespace-separated recognition site and regular expression
import re
import BeginPythonBioinfo     # see Chapter 6 about this module

def parseREBASE(rebasefile):  
    
    rebase_hash = {}

    # Read in the REBASE file
    rebasefile = BeginPythonBioinfo.get_file_data(rebasefile)

    foreach line in rebasefile:

        line = line.rstrip()

        # Discard header lines
        if line.find('Rich Roberts') != -1: continue

        # Discard blank lines
        if re.match('^s*$', line): continue
    
        # Split the two (or three if includes parenthesized name) fields
        fields = line.split(" ")

        # Get and store the name and the recognition site

        # Remove parenthesized names, for simplicity's sake,
        # by not saving the middle field, if any,
        # just the first and last
        name = fields[0]

        site = fields[-1]

        # Translate the recognition sites to regular expressions
        regexp = BeginPythonBioinfo.IUB_to_regexp(site)

        # Store the data into the hash
        rebase_hash[name] = "%s %s" % (site, regexp)

    # Return the hash containing the reformatted REBASE data
    return rebase_hash

 2.4 查找限制酶切位点

现在是时候编写一个主程序并查看我们的代码，我们从一个小伪代码开始，看看还有什么需要做：

#
# Get DNA
#
get_file_data

extract_sequence_from_fasta_data

#
# Get the REBASE data into a hash, from file "bionet"
#
parseREBASE('bionet');

for each user query

    If query is defined in the hash
        Get positions of query in DNA

    Report on positions, if any

}

 您现在需要编写一个函数，在DNA中寻找查询的位置。

Given arguments query and dna
return all matched positions

让我们来看看re模块文档找点线索，尤其是文档中的内置函数列表。看起来，finditer函数可以解决这个问题。它给出了所有位点组成的一个迭代器iterator。例子9.3演示了主程序，其后是需要用到的函数。

例子9-3 为用户的查询制作限制酶切图谱

#!/usr/bin/env python2.7
# Make restriction map from user queries on names of restriction enzymes

import re
import BeginPythonBioinfo;     # see Chapter 6 about this module



# Read in the file "sample.dna"
file_data = BeginPythonBioinfo.get_file_data("sample.dna")

# Extract the DNA sequence data from the contents of the file "sample.dna"
dna = BeginPythonBioinfo.extract_sequence_from_fasta_data(file_data)

# Get the REBASE data into a dict, from file "bionet"
rebase_dict = BeginPythonBioinfo.parseREBASE('bionet')

# Prompt user for restriction enzyme names, create restriction map
while True:
    print "Search for what restriction site for (or quit)?: ";
    
    query = input()

    query = query.strip()

    # Exit if empty query
    if not query or query== 'quit' : exit()

    # Perform the search in the DNA sequence
    if query in rebase_dict:

        recognition_site, regexp = rebase_hash[query].split(" ")

        # Create the restriction map
        locations = [pos for pos in re.finditer(regexp, dna)]

        # Report the restriction map to the user
        if locations:
            print("Searching for %s %s %s
" % (query, recognition_site, regexp)
            print("A restriction site for %s at locations:
" % query)
            print(" ".join([str(pos.start()+1) for pos in locations]))
        else:
            print("A restriction site for %s is not in the DNA:
" % query)
        
    print("
")

exit()

 如下是例子

Search for what restriction enzyme (or quit)?: AceI
Searching for AceI G^CWGC GC[AT]GC
A restriction site for AceI at locations:
54 94 582 660 696 702 840 855 957

Search for what restriction enzyme (or quit)?: AccII
Searching for AccII CG^CG CGCG
A restriction site for AccII at locations:
181

Search for what restriction enzyme (or quit)?: AaeI
A restriction site for AaeI is not in the DNA:

Search for what restriction enzyme (or quit)?: quit

注意程序中的列表推导式[pos for pos in re.finditer(regexp, dna)]，re.finditer(regexp, dna)是在正则表达式匹配成功后返回的一个迭代器iterator，它包含dna中所有的查询位点。每个位点是一个SRE_Match对象，调用对象的start()函数，可以获取匹配序列后面第一个碱基的索引，然后加1，就是碱基的位置。

3. Python的运算

在本专辑中，我们没有讨论基本的算术运算，就已经做的很好了，因为实际上你并不需要比加减乘除更多操作的内容。然后，这是Python在内编程语言的一个重要部分，就是进行数学计算的能力。详情请参考python相关文档。

3.1 运算和括号的优先级

运算有一套优先级规则。这使得语言可以在一行中有多个运算时决定先进行哪个运算。运算的顺序会改变运算的结果，就像下面的例子演示的那样。

假设你有 8 + 4 / 2这样的代码。如果你先进行除法预算，你会得到 8 + 2，或者是 10。然而，如果你先进行加法运算，你就会得到 12 / 2，或者是6。

现在编程语言对运算都赋予了优先级。如果你知道这些优先级，你可以编写 8 + 4 / 2这样的表达式，并且你知道它的运算结果。但是这并不可靠。

首先，要是你得到了错误的结果会怎样呢？或者，要是其他查看代码的人并没有你这样的记忆力会怎样呢？再或者，要是你记住了一种语言的优先级、但Python采用的是不同的优先

级会怎样呢？（不同的语言确实有不同的优先级规则。）

有一种解决办法，这就是使用括号。对于例子9.3，如果你简单的添加上括号： 8 + ( 4 / 2))，对于你自己、其他读者以及Python程序来说，都可以明确无误地知道你想首先进行除法

运算。注意包裹在另外一对括号中的“内部”括号，会被首先求值。

记住，在复杂的表达式中使用括号来指明运算的顺序。另外，它会让你避免大量的程序调试！

练习题

9.1 修改例子9-3，让它可以从命令行接收DNA；如果它没有被指定，提示用户键入FASTA文件名并读入DNA序列数据。
9.2 修改习题9.1，从bionet文件中读入全部的REBASE限制性内切位点数据，并把它制作成一个散列。
9.3 修改习题9.2，如果不存在DBM文件，就通过存储的REBASE散列生成一个DBM文件，如果DBM文件存在，就直接使用DBM文件。（提前查阅第十章中关于DBM的更多的信息。）
9.4 修改例子5.3，报告它找到的基序的位置，即使基序在序列数据中出现多次。
9.5 通过打印出序列、并用酶的名字把限制性内切位点标记出来，为限制酶切图谱添加一个切割位点的图形化展示。你能制作一个处理多个限制性内切酶的图谱吗？你该如何处理重叠的限制性内切位点呢？
9.6 编写一个子程序，返回限制酶切消化片段，就是进行限制性酶切反应后留下的DNA片段。记住要考虑切割位点的位置。（这需要你以另一种方式解析REBASE的bionet文件。如果你愿意的话，可以把没有用^~指明切割位点的限制性内切酶忽略掉。）
9.7 扩展限制酶切图谱软件，对于非回文识别位点，把相反的另一条链也考虑在内。
9.8 对于没有括号的算术运算，编写一个子程序，根据Perl的优先级规则为其添加上合适的括号。（警告：这是一个相当有难度的练习题，除非你确信有足够的课余时间，否则请跳过该题。对于优先级规则，可以参看Perl的文档。）