Introduction
此课程是Coursera上，John Hopkins University 的 Genomic Data Science的Course4，Course 3主要讲的是python的基础，因此直接略过，Course 4讲的就是和测序比对算法有关的东西了，并且还给出了python应用，这里偏理论，挺烧脑的。
课程链接：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不断向右偏移，直到找到匹配结果，在这个例子中，word需要从头偏移40次才可匹配上结果。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比较，但是依据以下两个不同的规则，可以跳过一些非必要的比较，从而节省比较数。
移动字符的数目根据两条规则决定（t表示P和T比对上的后缀）：
①坏字符规则：P和T从右向左比较，当遇到第一个不匹配字符时，从此位置向左在P中寻找到下一个在此位置匹配上T的字符，然后移动P，如果找不到，则将P向右移动到不匹配字符的右边。
②好后缀规则：1，在P中i位置的左侧最靠近i位置查找字符串t'使得t'=t（此时t'不是P的后缀，实际上也就是查找匹配到的字符串除了在P的后缀中存在，是否在P的其他位置存在），若存在，则移动相应的位数将找到的t'与T中的t对齐。
2，如果t'不存在，那我们继续查找t的某一个后缀是否为P的前缀，若存在，则移动相应的位将P的前缀与t的后缀位置对齐。否则，将P向后移动n个字符。
好后缀规则的实质是，将不匹配位置右侧匹配到的字符串t的所有字符后缀组合，用于查找它们在P的不匹配位置左侧是否存在。
通俗的理解是，最长的好后缀t是否存在于i的左侧（情况1），其他后缀组合中是否存在与P的前缀相同的情况（情况2）。
Example
Step 1
T: GTTATAGCTGATCGCGGCGTAGCGGCGAA
P: GTAGCGGCG
我们可以发现，T的最后一个T和P的最后一个G不匹配，并且没有匹配上后缀，因此应用坏字符原则，让T的T匹配上P的T
记bc:6, gs=0 表明坏字符准则能跳过6步比对，而好后缀无法跳过
Step 2
T: GTTATAGCTGATCGCGGCGTAGCGGCGAA
P: GTAGCGGCG
有了匹配上的好后缀t="GCG"，应用好后缀准则可以匹配到P的左边另一个GCG，此时坏字符无法跳过
则有bc:0 gs:2
Step 3
T: GTTATAGCTGATCGCGGCGTAGCGGCGAA
P: GTAGCGGCG
此时应用bc和gs
bc:2 gs:7
因此应用gs更好（注意，gs准则只需要要求你P在gs的左边有匹配上T中gs后缀的字符即可，而非要完全匹配）
Step4
T: GTTATAGCTGATCGCGGCGTAGCGGCGAA
P: GTAGCGGCG
至此匹配结束，相比Naive algorithm，我们一共跳过了6+2+7=15步。
Python code由于太长，这里就不放出了。
Indexing DNA
为了方便比对和查找，通常需要对搜索的文本T（即Genome)构建索引，这里需要引入k-mer这个概念，即长度为k的子字符串。
以5-mer为例，对于文本T=CGTGCGTGCTT：
Index of T:
CGTGC: 0, 4
GCGTG: 3
GTGCC: 1
GTGCT: 5
TGCCT: 2
TGCTT: 6
注意这里还按照字母顺序进行了排序，后面的数字代表对应的索引号。
以P:GCGTGC为例查找T，P有2个5-mer，GCGTG和CGTGC，首先会以第一个查找到3号索引，然后接下来进行向右延伸的步骤称之为Verification，结果C匹配，因此成功。如果不成功则找下一个匹配的，如果找不到，则此即为最佳匹配结果。
Binary search
或许你会好奇为什么构建索引还需要排序，二分查找算法可以给你一个合理的解释，
Example:
T: GTGCGTGG
P: GCGTGG
对T构建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，我们直接将表一分为二，找到表中间位置对应的3-mer为GTG，由于按照字母顺序，TGG排列在GTG后面，因此GTG前面的可以直接删去不考虑，后面循环同样的步骤一分为二，直到查找到匹配结果，每次查找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
可见，这个函数可以输出从左边插入到列表中你指定的数字，保持列表顺序不变的最小偏移量。那么实际上我们可以把数字换成字母，这样就对应了上面的二分查找算法找对应k-mer匹配，即bisect.bisect_left(index, "GTG")
python code:
import bisect
class Index(object):
def __init__(self, t, k):
"""
此函数构建了长度为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），是现在广泛使用的比对软件BWA、Bowtie2所使用的算法。（这里作者没有介绍BWT，我觉得刘小乐老师对此的介绍比较通俗易懂，b战也有搬运老师油管上传的课程，见https://www.bilibili.com/video/BV1By4y1a7GM?p=8，刘小乐老师在其油管频道中经常更新其每年的bioinformatics课程，感兴趣可以看看）
Approximate matching
前面讲了一些精确比对的算法，但是由于测序错误，自然变异等原因，两个序列很难完全比对上，因此有了近似比对的概念，且引入两个距离概念
Hamming distance
对于长度相等的序列X、Y，hamming distance = 使得两序列相同的最小替换数目
如下
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相等的最小数目
X: ATGCC
Y:TGAC
↓
X: ATGCC
Y: -TGAC
则Edit distance = 2
对于近似比对，假设我们以汉明距离作为标准，我们只需要改动上面naive算法的一点代码就可以了
def naiveHamming(p, t):
occurences = []
for i in range(len(t) - len(p) +1): # 字串P从头向右偏移匹配T
nmm = 0 # 错配计数
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的，即鸽巢原理。此原理可以应用在任意一种算法上，我们可以以其中一个精确比对作为靶点，然后向两侧延伸进行verification。
Edit distance以及Alignment
后续主要讲了Edit distance与Hamming distance的关系，以及一些性质，以及用Edit distance做动态规划算法以及全局和局部比对的算法，这些网上教学都是挺多的，因此不过多阐述。
Assembly
Shortest common superstring
组装算法需要我们根据所有k-mer重叠关系，来找出一条能够串起所有k-mer的最短路径，即Shortest common superstring问题，即对于给定的字符串集合S，找到SCS(S): 得到最短的将S子字符串连接在一起的字符串。
即如图所示：
对于SCS，这是一个NP-complete问题，即对于大量的输入，没有足够有效的算法（关于此概念，请参照Wiki中对于NP困难、NP的描述，NP完备是NP与NP困难的交集，是NP中最难的决定性问题，所有NP问题都可以被快速归化为NP完备问题。因此NP完备问题应该是最不可能被化简为P（多项式时间可决定）的决定性问题的集合。若任何NPC问题得到多项式时间的解法，那此解法就可应用在所有NP问题上。)。
接下来有一种优化的方法，即对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告诉我们的结果。
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)
