第九章 限制性图谱和正则表达式
在本章中，我将概述Python正则表达式和Python运算符。咱们还将研究标准的基本分子生物学技术的编程：发现序列的限制性图谱。限制性消化是“指纹”DNA的原始方法之一；如今能够在计算机上模拟这个。html
限制性图及其相关的限制性位点是实验室中的常见计算，由几个软件包提供。它们是克隆实验的重要工具；例如，它们可用于将所需的DNA片断插入克隆载体中。限制图也可用于测序项目，例如鸟枪和定向测序。python
1. 正则表达式模块re
咱们一直在处理正则表达式，本节填补了一些背景知识，并将本书前面部分正则表达式的分散讨论联系在一块儿。使用生物数据（如序列）或GenBank、PDB和BLAST文件进行编程时，正则表达式很是有用。正则表达式
正则表达式是用一个字符串表示和搜索许多字符串的方法。虽然它们并非彻底相同的东西，但将正则表达式看做是一种高度发达的通配符集是颇有用的。正则表达式中的特殊字符更恰当地称为元字符。数据库
大多数人都熟悉通配符，这些通配符能够在搜索引擎或扑克游戏中找到。例如，你能够经过键入biolog*找到对biolog开头的每一个单词的引用。或者你可能会发现本身拿着五个ace。express
在计算机科学中，这些通配符或元字符在实践和理论上都有重要的历史。特别是星号字符在发明它的杰出逻辑学家以后称为Kleene闭包。编程
咱们已经看到许多使用正则表达式来查找DNA或蛋白质序列中的事物的例子。在这里，我将简要介绍正则表达式背后的基本思想，做为一些术语的介绍。数组
所以，让咱们从一个实际的例子开始：使用字符类来搜素DNA。假设你想在DNA库中找到一个小基序长度为6个碱基对：CT后跟C或G或T，而后是ACG。使用正则表达式表示该基序，以下所示：数据结构
import re
if re.match("CT[CGT]ACG", dna):
print "I found the motif!!\n"
重复或星号
星号（*），也称为Kleene闭包或星号，表示在它以前的字符重复0次或屡次。例如，abc*匹配如下任何字符串：ab，abc，abcc，abccc，abcccc等。正则表达式匹配无限数量的字符串。闭包
轮换
在Python中，模式（a|b）（读：a或b）匹配字符串a或字符串b。python2.7
并列
这是一个很是明显的问题。在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\n" % (query, recognition_site, regexp)
print("A restriction site for %s at locations:\n" % query)
print(" ".join([str(pos.start()+1) for pos in locations]))
else:
print("A restriction site for %s is not in the DNA:\n" % query)
print("\n")
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的文档。）
- 1. 第九章 JS正则表达式
- 2. 限制密码的正则表达式
- 3. 高性能Javascript第五章字符串和正则表达式
- 4. 第十六章 Python正则表达式
- 5. 第三章 正则表达式
- 6. 精通正则表达式：第一章
- 7. 第七章 正则表达式编程
- 8. 第9章 正则表达式
- 9. 第3章：正则表达式
- 10. 第一章：正则表达式
- 更多相关文章...
- • PHP 正则表达式(PCRE) - PHP参考手册
- • Scala 正则表达式 - Scala教程
- • Docker容器实战(六) - 容器的隔离与限制
- • IntelliJ IDEA 代码格式化配置和快捷键
-
每一个你不满意的现在，都有一个你没有努力的曾经。
- 1. 第九章 JS正则表达式
- 2. 限制密码的正则表达式
- 3. 高性能Javascript第五章字符串和正则表达式
- 4. 第十六章 Python正则表达式
- 5. 第三章 正则表达式
- 6. 精通正则表达式：第一章
- 7. 第七章 正则表达式编程
- 8. 第9章 正则表达式
- 9. 第3章：正则表达式
- 10. 第一章：正则表达式