Introduction
此課程是Coursera上入撒，John Hopkins University 的 Genomic Data Science的Course4荒吏，Course 3主要講的是python的基礎(chǔ)河泳，因此直接略過嬉探，Course 4講的就是和測序比對算法有關(guān)的東西了臼寄，并且還給出了python應(yīng)用伪朽，這里偏理論行您，挺燒腦的铭乾。
課程鏈接：https://www.coursera.org/specializations/genomic-data-science
所有圖片、Python code均出自此課程
Exact matching
Naive algorithm
P: word
T: There would have been a time for such a word
對于Pattern(P)：word娃循，在Text(T)中尋找到匹配炕檩，按照naive alogrithm算法，實際就是把P和T從頭比對上捌斧，然后把P不斷向右偏移笛质，直到找到匹配結(jié)果，在這個例子中捞蚂，word需要從頭偏移40次才可匹配上結(jié)果妇押。Python代碼如下
def naive(p, t):
occurences = []
for i in range(len(t) - len(p) +1): # 字串P從頭向右偏移匹配T
match = True
for j in range(len(p)):
if t[i+j] != p[j]: # 如果P和T有一個不匹配則跳過
match = False
break
if match:
occurences.append(i)
return occurrences
返回的是匹配到的index下表列表。
Boyer-Moore string-search algorithm
嗶哩嗶哩也有搬運此課講解的視頻翻譯:https://www.bilibili.com/video/BV1j4411P7kk?from=search&seid=16892982446630762406
同樣也是從頭把P和T比較洞难，但是依據(jù)以下兩個不同的規(guī)則舆吮，可以跳過一些非必要的比較，從而節(jié)省比較數(shù)队贱。
移動字符的數(shù)目根據(jù)兩條規(guī)則決定（t表示P和T比對上的后綴）：
①壞字符規(guī)則：P和T從右向左比較色冀，當遇到第一個不匹配字符時，從此位置向左在P中尋找到下一個在此位置匹配上T的字符柱嫌，然后移動P锋恬，如果找不到，則將P向右移動到不匹配字符的右邊编丘。
②好后綴規(guī)則：1与学，在P中i位置的左側(cè)最靠近i位置查找字符串t'使得t'=t（此時t'不是P的后綴彤悔，實際上也就是查找匹配到的字符串除了在P的后綴中存在，是否在P的其他位置存在）索守，若存在晕窑，則移動相應(yīng)的位數(shù)將找到的t'與T中的t對齊。
2卵佛，如果t'不存在杨赤，那我們繼續(xù)查找t的某一個后綴是否為P的前綴，若存在截汪，則移動相應(yīng)的位將P的前綴與t的后綴位置對齊疾牲。否則，將P向后移動n個字符衙解。
好后綴規(guī)則的實質(zhì)是阳柔，將不匹配位置右側(cè)匹配到的字符串t的所有字符后綴組合，用于查找它們在P的不匹配位置左側(cè)是否存在蚓峦。
通俗的理解是舌剂，最長的好后綴t是否存在于i的左側(cè)（情況1），其他后綴組合中是否存在與P的前綴相同的情況（情況2）枫匾。
Example
Step 1
T: GTTATAGCTGATCGCGGCGTAGCGGCGAA
P: GTAGCGGCG
我們可以發(fā)現(xiàn)架诞，T的最后一個T和P的最后一個G不匹配，并且沒有匹配上后綴干茉，因此應(yīng)用壞字符原則谴忧，讓T的T匹配上P的T
記bc:6, gs=0 表明壞字符準則能跳過6步比對，而好后綴無法跳過
Step 2
T: GTTATAGCTGATCGCGGCGTAGCGGCGAA
P: GTAGCGGCG
有了匹配上的好后綴t="GCG"角虫，應(yīng)用好后綴準則可以匹配到P的左邊另一個GCG沾谓，此時壞字符無法跳過
則有bc:0 gs:2
Step 3
T: GTTATAGCTGATCGCGGCGTAGCGGCGAA
P: GTAGCGGCG
此時應(yīng)用bc和gs
bc:2 gs:7
因此應(yīng)用gs更好（注意，gs準則只需要要求你P在gs的左邊有匹配上T中g(shù)s后綴的字符即可戳鹅，而非要完全匹配）
Step4
T: GTTATAGCTGATCGCGGCGTAGCGGCGAA
P: GTAGCGGCG
至此匹配結(jié)束均驶，相比Naive algorithm，我們一共跳過了6+2+7=15步枫虏。
Python code由于太長妇穴，這里就不放出了。
Indexing DNA
為了方便比對和查找隶债，通常需要對搜索的文本T（即Genome)構(gòu)建索引腾它，這里需要引入k-mer這個概念，即長度為k的子字符串死讹。
以5-mer為例瞒滴，對于文本T=CGTGCGTGCTT：
Index of T:
CGTGC: 0, 4
GCGTG: 3
GTGCC: 1
GTGCT: 5
TGCCT: 2
TGCTT: 6
注意這里還按照字母順序進行了排序，后面的數(shù)字代表對應(yīng)的索引號赞警。
以P:GCGTGC為例查找T妓忍，P有2個5-mer虏两，GCGTG和CGTGC，首先會以第一個查找到3號索引世剖，然后接下來進行向右延伸的步驟稱之為Verification定罢，結(jié)果C匹配，因此成功搁廓。如果不成功則找下一個匹配的引颈，如果找不到，則此即為最佳匹配結(jié)果境蜕。
Binary search
或許你會好奇為什么構(gòu)建索引還需要排序，二分查找算法可以給你一個合理的解釋凌停，
Example:
T: GTGCGTGG
P: GCGTGG
對T構(gòu)建3-mer索引表如下:
CGT 3
GCG 2
GGG 8
GGG 9
GGG 10
GTG 0
GTG 4
GTG 6
TGC 1
TGG 7
TGT 5
假如我們要搜索P中的3-mer TGG粱年，我們直接將表一分為二，找到表中間位置對應(yīng)的3-mer為GTG罚拟，由于按照字母順序台诗，TGG排列在GTG后面，因此GTG前面的可以直接刪去不考慮赐俗，后面循環(huán)同樣的步驟一分為二拉队，直到查找到匹配結(jié)果，每次查找query只需要1og2(n) bisections阻逮，大大降低算法復雜度粱快。
python的bisect模塊可以很好的給我們演示：
a = [1, 3, 3, 6, 8, 8, 9, 10]
import bisect
bisect.bisect_left(a, 2)
# 輸出 1
bisect.bisect_left(a, 4)
# 輸出 3
bisect.bisect_left(a, 8)
# 輸出 4
可見，這個函數(shù)可以輸出從左邊插入到列表中你指定的數(shù)字叔扼，保持列表順序不變的最小偏移量事哭。那么實際上我們可以把數(shù)字換成字母，這樣就對應(yīng)了上面的二分查找算法找對應(yīng)k-mer匹配瓜富，即bisect.bisect_left(index, "GTG")
python code:
import bisect
class Index(object):
def __init__(self, t, k):
"""
此函數(shù)構(gòu)建了長度為k的對于t的k-mer排序索引列表
"""
self.k = k
self.index = []
for i in range(len(t) - k + 1):
self.index.append((t[i: i+k], i))
self.index.sort()
def query(self, p):
"""
查找比對上的索引下標
"""
kmer = p[:self.k]
i = bisect.bisect_left(self.index, (kmer, -1))
hits = []
while i < len(self.index):
if self.index[i][0] != kmer:
break
hits.append(self.index[i][1])
i += 1
return hits
def queryIndex(p, t, index):
k = index.k
offsets = []
for i in index.query(p):
if p[k:] == t[i+k: i+len(p)]:
offsets.append(i)
return offsets
t = "GCTACGATCTAGAATCTA"
p = "TCTA"
index = Index(t, 2)
print(queryIndex(p, t, index))
# 輸出 [7, 14]
Hash tables for indexing
除了像上述的有序列表外鳍咱，還可以用無序的哈希表儲存索引，差別就在于無序与柑，以python code為例
t = "GTGCGTGTGGGGG"
table = {"GTG": [0, 4, 6], "TGC": [1],
"GCG": [2], "CGT": [3], "TGT": [5],
"TGG": [7], "GGG": [8, 9, 10]}
table["GGG"]
# 輸出 [8, 9, 10]
table["CGT"]
# 輸出 [3]
Suffix index
還有基于后綴的索引方法
然后再基于二分查找：
另外還有其他更多Suffix的索引方法谤辜，以及所需的大小，第三個即BWT（Burrows-Wheeler transform）价捧，是現(xiàn)在廣泛使用的比對軟件BWA丑念、Bowtie2所使用的算法。（這里作者沒有介紹BWT干旧，我覺得劉小樂老師對此的介紹比較通俗易懂渠欺，b戰(zhàn)也有搬運老師油管上傳的課程，見https://www.bilibili.com/video/BV1By4y1a7GM?p=8椎眯，劉小樂老師在其油管頻道中經(jīng)常更新其每年的bioinformatics課程挠将，感興趣可以看看）
Approximate matching
前面講了一些精確比對的算法胳岂，但是由于測序錯誤，自然變異等原因舔稀，兩個序列很難完全比對上乳丰，因此有了近似比對的概念，且引入兩個距離概念
Hamming distance
對于長度相等的序列X内贮、Y产园，hamming distance = 使得兩序列相同的最小替換數(shù)目
如下
X: GAGGTAGCGGCGTTTAAC
Y: GTGGTAACGGGGTTTAAC
hamming distance = 3個替換 = 3
python code:
def hammingDistance(x, y):
nmm = 0
for i in range(0, len(x)):
if x[i] != y[i]:
nmm += 1
return nmm
Edit distance(AKA Levenshtein distance)
對于X、Y夜郁，edit distance = 通過插入什燕、刪除、替換竞端，使得X屎即、Y相等的最小數(shù)目
X: ATGCC
Y:TGAC
↓
X: ATGCC
Y: -TGAC
則Edit distance = 2
對于近似比對，假設(shè)我們以漢明距離作為標準事富，我們只需要改動上面naive算法的一點代碼就可以了
def naiveHamming(p, t):
occurences = []
for i in range(len(t) - len(p) +1): # 字串P從頭向右偏移匹配T
nmm = 0 # 錯配計數(shù)
match = True
for j in range(len(p)):
if t[i+j] != p[j]: # 如果P和T有一個不匹配則跳過
nmm += 1
if nmm > maxDisctance:
break
if nmm <= maxDistance:
occurences.append(i)
return occurrences
Pigeonhole Principle
如果P相對于T有k edits距離技俐，將P分為k+1份，則至少在p1 p2 ... pk+1當中有一個是0 edits的统台，即鴿巢原理雕擂。此原理可以應(yīng)用在任意一種算法上，我們可以以其中一個精確比對作為靶點贱勃，然后向兩側(cè)延伸進行verification井赌。
Edit distance以及Alignment
后續(xù)主要講了Edit distance與Hamming distance的關(guān)系，以及一些性質(zhì)募寨，以及用Edit distance做動態(tài)規(guī)劃算法以及全局和局部比對的算法族展，這些網(wǎng)上教學都是挺多的，因此不過多闡述拔鹰。
Assembly
Shortest common superstring
組裝算法需要我們根據(jù)所有k-mer重疊關(guān)系仪缸，來找出一條能夠串起所有k-mer的最短路徑，即Shortest common superstring問題列肢，即對于給定的字符串集合S恰画，找到SCS(S): 得到最短的將S子字符串連接在一起的字符串。
即如圖所示：
對于SCS瓷马，這是一個NP-complete問題拴还，即對于大量的輸入，沒有足夠有效的算法（關(guān)于此概念欧聘，請參照Wiki中對于NP困難片林、NP的描述，NP完備是NP與NP困難的交集，是NP中最難的決定性問題费封，所有NP問題都可以被快速歸化為NP完備問題焕妙。因此NP完備問題應(yīng)該是最不可能被化簡為P（多項式時間可決定）的決定性問題的集合。若任何NPC問題得到多項式時間的解法弓摘，那此解法就可應(yīng)用在所有NP問題上焚鹊。)。
接下來有一種優(yōu)化的方法韧献，即對S當中的字符串進行排列末患，對每一種排列，從頭到尾锤窑，計算得到SCS璧针，然后比較所有順序的SCS的長度，選取最小果复，如果S包含n個strings陈莽，則有n!種排列。
python code:
import itertools
def overlap(a, b, min_length=3):
"""Return length of longest suffix of 'a' matching
a prefix of 'b' that is at least 'min_length'
characters long. if no such overlap exists,
return 0. """
start = 0 # start all the way at the left
while True:
start = a.find(b[:min_length], start) # look for b's suffix in a
if start == -1: # no more occurrences to right
return 0
# found occurrence; check for full suffix/prefix match
if b.startswith(a[start:]):
return len(a) - start
start += 1 # move just past previous match
def scs(ss):
shortest_sup = None
for ssperm in itertools.permutations(ss):
sup = ssperm[0]
for i in range(len(ss) - 1):
olen = overlap(ssperm[i], ssperm[i+1], min_length=1)
sup += ssperm[i+1][olen:]
if shortest_sup is None or len(sup) < len(shortest_sup):
shortest_sup = sup
return shortest_sup
Greedy shortest common superstring
貪婪算法即找每步當中重疊最多的reads pair合并虽抄，如下圖所示：
但問題在貪婪算法找到的并不是最短的superstring。
python code(依賴于上個代碼塊)：
def pick_maximal_overlap(reads, k):
reada, readb = None, None
best_olen = 0
for a, b in itertools.permutations(reads, 2):
olen = overlap(a, b, min_length=k)
if olen > best_olen:
reada, readb = a, b
best_olen = olen
return reada, readb, best_olen
def greedy_scs(reads, k):
read_a, read_b, olen = pick_maximal_overlap(reads, k)
while olen > 0:
reads.remove(read_a)
reads.remove(read_b)
reads.append(read_a + read_b[olen:])
read_a, read_b, olen = pick_maximal_overlap(reads, k)
return "".join(reads)
greedy_scs(["ABCD", "CDBC", "BCDA"], 1)
# "CDBCABCDA"
scs(["ABCD", "CDBC", "BCDA"])
# "ABCDBCDA"
De Bruijn graphs
走過每條邊正好一次独柑，就可以得到基因組的重建迈窟，這個即是Eulerian walks告訴我們的結(jié)果。
ipython code:
def de_bruijn_ize(st, k):
edges = []
nodes = set()
for i in range(len(st) - K +1):
edges.append((st[i:i+k-1], st[i+1:i+k]))
nodes.add(st[i:i+k-1])
nodes.add(st[i+1:i+k])
return nodes, edges
nodes, edges = de_bruijn_ize("ACGCGTCG", 3)
print(nodes)
# {'TC', 'CG', 'AC', 'GT', 'GC'}
print(edges)
# [('AC', 'CG'), ('CG', 'GC'), ('GC', 'CG'), ('CG', 'GT'), ('GT', 'TC'), ('TC', 'CG')]
def visualize_de_bruijn(st, k):
nodes, edges = de_bruijn_ize(st, k)
dot_str = 'digraph "DeBruijn graph" {\n'
for node in nodes:
dot_str += ' %s [label="%s"] ;\n' % (node, node)
for src, dst in edges:
dot_str += ' %s -> %s;\n' % (src, dst)
return dot_str + '}\n'
%load_ext gvmagic
%dotstr visualize_de_bruijn("ACGCGTCG", 3)
