【陈巍翻译】视频A04：序列简单比对

（如果要看清视频中的代码，建议在文章底部，复制高清视频链接，用电脑来看。）

字 幕 内 容

在这项实践中，我们会写一个（碱基）完全匹配的算法，就如讲座中所显示的，并用它来，把人工做出来的read比对到基因组上。

我们从下载一个PhiX生物体（PhiX是一种噬菌体）的基因组开始。

我们运行“wget”（wget是之前的课程中写的一个下载基因组数据的命令），下载它。

现在，我要做的是用前面的练习中写的“readgenome”函数，来读这个（基因组）。我不必重新实现它（因为之前已写好代码了）。我们要用它来读我们刚才下载的基因组。

现在我们开始做。让我们来写这个“naive”的匹配算法（naive：天真）。我要写一个函数“naive”。这里会有2个参数，一个是“p”，代表我们正在找的模式（pattern）；而“t”是我们要在其中寻找（模式）的文本。例如，“p”可以是我们读到的一段序列，而“t”可以是我们要用来作为对照的基因组。

所以，我会创建一个跟踪“p”和“t”匹配发生情况的列表。我要做的事情，是把每个“p”可以开始的位置，都做一次遍历。对于“i”，从0到“t”的全长，减去“p”的长度，加上1（进行遍历）。这可以给我，在“t”中 所有“p”可以开始的位置，而不会超过“t”的结尾。

我要创建一个布尔变量“match”，它的初始值是“真”。我要每“p”中的每个字符，都与“t”中对应位置的字符相比较。如果有一个不匹配的，我就把“match”设成“假”。

所以，我就写“对于在p范围内的每个j”，如果“t”中在“i+j”位置的那个字符，与“p”中在“j”位置的字符不相等，那么，我们把“match”设成假。

你比较了在“p”中第“j”位置的字符，与在“t”中的“i+j”位置的字符。你记下了在“i”这个位置（的匹配情况）。

我们对“i”和“j”都进行位移比对，因为我们既比什么是比对得上的，又因为我们比对哪个字符是比对得上的。

好。如果我们发现一个不匹配，我们就把“match”设成假。并且，我们就不需要再在这个位置比对之后的字符了。

因为我们已经发现了一个不匹配（的碱基），我们可以中断这个“For”了。

现在，当我完成了这个（for），如果“match”变量依然是真，这意味着，“p”在那个位置一定是完全匹配到“t”的。这种情况下，我们把索引值“i”加到我们的匹配列表中。

然后，当这些都做好之后，我们遍历每个可能的位置，返回一个匹配的列表。这样，我们就得到了每一个在“t”中与“p”匹配的位置。

好。让我们试一下，让我们建立了一个小的碱基文本字符串。我会做一个短的模式。我就选“AG”。

我会在这个模式（p）、和文本（t）中运行“naive”函数。我们可以看到，它得到了3个不同的地方，都是“p”和“t”能匹配的。我们可以看这些index位置，来检查这些序列是否匹配：

这是0到2

你可以看5到7

我们可以看9到11

是的，在那个位置“t”含有字母“AG”

现在，我们有了这个实现方案，我们也找到了完全匹配的位置。

我们要从基因组中随机地产生一些序列。我们已经有了一个可以做这件事的函数我们把这个函数粘贴到这里。这个函数所做的事情，是从基因组中产生人工的序列，通过从基因组序符串中随机位置取得一个子序列。所有这些人工生成的序列，都应该能够完全匹配到基因组上。

让我们来用这个方法生成一堆的序列。我要用我们已经下载好的PhiX基因组。我们要生成100条100个碱基长度的序列。

一旦我们做好了，我们就要数一下，有几条序列完全匹配到原来的基因组上了。因为它们是人工生成的，所以没有测序错误或突变，它们应该是完全匹配的。

所以，现在我要创建一个变量“numMatched”来数匹配上的序列数。我把“numMatched”变量初始化成“0”。对于列表中的每个序列，我都要看它遍历匹配到“genome”后的“matches”（的变量结果值）。

我要先放上模式“p”，我要开始用“r”代表序列（read）。然后是基因组“genome”，也就是要比较的文本“t”。这样，序列就是模式“p”，而基因组“genome”就是文本“text”/“t”。

现在，“matches”就是序列“r”与“genome”是否匹配的指标。如果这个列表是0，那就意味着完全没有匹配上。所以，如果matches的长度比0大，就意味着我们找到了匹配的read。

在这个例子中，我们要增加我们的变量“numMatched”。

当这件事做好之后，我要打印出匹配上的read的比例。所以，我要打印。那会是匹配的read中匹配的碱基数。然后，把变量放在这里，第一个是numMatched（匹配上的序列数量）。第二个是匹配上的read的长度。

所以，你运行这个函数，100个read是匹配的，每个read都有100个碱基是与基因组匹配的。

很好，人工的reads与基因组完全匹配。

------------

高清视频地址：http://v.qq.com/page/j/0/n/j0342lulv0n.html

原始视频地址：https://www.youtube.com/watch?v=ep91JWd6fs0&list=PLMYiAoWjBNaE767PvAp8OXzZLtPMBCt5c&index=1

程序源代码地址：https://github.com/BenLangmead/ads1-notebooks