第九章 限制性圖譜和正則表達式
在本章中，我將概述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. 字節跳動21屆秋招運營兩輪面試經驗分享
- 2. Java 3 年，25K 多嗎？
- 3. mysql安裝部署
- 4. web前端開發中父鏈和子鏈方式實現通信
- 5. 3.1.6 spark體系之分佈式計算-scala編程-scala中trait特性
- 6. dataframe2
- 7. ThinkFree在線
- 8. 在線畫圖
- 9. devtools熱部署
- 10. 編譯和鏈接
- 1. 第九章 JS正則表達式
- 2. 限制密碼的正則表達式
- 3. 高性能Javascript第五章字符串和正則表達式
- 4. 第十六章 Python正則表達式
- 5. 第三章 正則表達式
- 6. 精通正則表達式：第一章
- 7. 第七章 正則表達式編程
- 8. 第9章 正則表達式
- 9. 第3章：正則表達式
- 10. 第一章：正則表達式