1.概念

概念：序列比对是生物信息学中的一种技术，用于比较和排列DNA、RNA或蛋白质序列，以找到它们之间的相似性和差异性。

我们可以把它看成一个找“共同点”的过程，就像找两首歌的歌词相似之处或者找两个拼图的相似部分。

2.为什么要做序列比对

在生物学研究中，我们希望知道不同生物（或同一生物的不同个体）之间是否具有相似的基因。因为相似的基因往往有相似的功能，序列比对可以帮助我们找到基因之间的联系，进而预测基因的功能或研究它们的进化关系。

基因同源性分析：通过比对基因或蛋白质序列，来确定它们是否来源于共同的祖先。相似度越高，说明它们之间的进化关系越近。

功能预测：如果一个未知功能的基因序列与已知基因序列有很高的相似性，那么可以推测它们可能具有相似的功能。

结构预测：比对蛋白质序列后，可以用结构相似的蛋白质模型来预测其三维结构。

进化研究：通过比对不同物种间的基因序列，分析它们的进化关系、分化时间和物种起源。

3.序列比对的基本类型

全局比对：从头到尾比较整个序列。适用于长度差不多的序列，试图比对所有的字符。比如把两篇文章从头到尾对比，查看它们的每一句话是否一样。

局部比对：寻找序列中局部的相似区域。适合比对长度差别大的序列，或者找出两个长序列中相似的片段。就像只对比两篇文章中相同的段落，忽略不相关的部分。

4.常用的比对算法和工具

BLAST（Basic Local Alignment Search Tool）：最常用的序列比对工具，用于局部比对，尤其适合快速查找。

Needleman-Wunsch算法：用于全局比对，通过动态规划找到最佳的比对结果。

Smith-Waterman算法：用于局部比对，也利用动态规划，但专注于找到最佳局部匹配。

算法可以理解为解决问题的“步骤清单”，和我们做菜的食谱有点类似。假如我们想用一套步骤从原材料做出一道菜，这些步骤的顺序和内容构成了算法。算法就是通过一系列规定的步骤，让我们从输入的数据开始，最终得到解决方案或答案。

举个例子：寻找最大数的算法：                                                                                                                                                                         假设我们有一个装有数字的列表，比如说 [3, 8, 2, 5, 9]，要找出里面最大的数字。我们可以设计一个简单的“算法”：假设第一个数是最大的数，比如这里假设 3 是最大的。把假设的数和列表中接下来的数字逐个比较。如果遇到比当前假设的数更大的数，就更新这个“假设”。一直比较到列表的最后一个数，最后留下的那个“假设”就是最大的数。这个简单步骤的组合就是一个算法，它能帮我们解决“找到最大数”的问题。

在生物信息学中，比对两个DNA序列的算法也是类似的“步骤清单”。比如说，Smith-Waterman算法用于比对两个DNA序列的相似片段，其核心是不断地比较两条序列中的字符，以找到匹配的部分。算法会逐步进行这样的比较、调整位置、计分，最后给出一个最优的比对结果。

4.1 BLAST

BLAST（Basic Local Alignment Search Tool）是生物信息学中非常重要的一个算法，用于快速查找DNA、RNA或蛋白质序列的相似片段。BLAST可以理解为一种“查找相似序列”的工具，其目的是在数据库中找到和查询序列相似的区域。下面将分步骤讲解BLAST算法的工作原理和它的特点。

BLAST的基本工作原理

BLAST算法的核心在于找到和查询序列最相似的“局部”区域，而不是从头到尾进行全局比对。其工作原理可以简单概括为以下几个步骤：

确定种子序列（Seed Sequence）

BLAST首先将查询序列分成多个小片段（通常为3-11个碱基或氨基酸的片段），称为“种子”。这些种子是查询序列中的短片段，比如DNA序列中的“三联体”（三个碱基一组），或蛋白质序列中的“三肽”（三个氨基酸一组）。

寻找相似的种子

BLAST会从数据库中搜索与这些种子相匹配的短片段。这一步的目标是找到和查询序列类似的片段，这样可以减少计算量。BLAST会在数据库中查找和这些种子序列完全匹配或部分匹配的片段。

扩展匹配区域

一旦找到相似的种子序列，BLAST就会尝试向两端扩展比对区域，即从种子的起始位置向左、向右逐步延伸，直到不再匹配。这个扩展过程能找到更大的相似区域，最终形成“局部比对”的片段。

打分与过滤

对于每个找到的相似区域，BLAST都会计算一个“相似得分”，通常根据匹配的数量和错配的位置来打分。相似得分越高，表示该区域与查询序列越相似。BLAST会对比对结果进行排序并输出最匹配的结果。

统计显著性

为了确保结果具有生物学意义，BLAST会计算一个“E值”（E-value），用来衡量结果的可靠性。E值越小，表示匹配的结果更有可能是“真实”的相似性，而不是随机发生的。

BLAST的几个常见变体

BLAST算法有几种不同的变体，适用于不同的序列类型：

blastn：用于DNA序列比对，比如比对基因组或核苷酸序列。

blastp：用于蛋白质序列比对，常用于研究蛋白质功能和相似性。

blastx：将DNA序列比对到蛋白质序列上，用于预测未知DNA片段的蛋白质功能。

tblastn：将蛋白质序列比对到DNA序列数据库中。

tblastx：将DNA翻译成蛋白质后进行比对，适合在两个未知的DNA序列之间查找相似性。

BLAST的优势

BLAST之所以受欢迎，是因为它比传统的比对算法更快。它不需要从头到尾比对所有碱基或氨基酸，而是通过局部比对的方式来减少计算量。这种方法不仅节省了时间，还提高了算法的效率。

举个例子：使用BLAST查找相似基因

假设我们发现一个新的基因序列，希望知道它是否与已知基因有相似性，可以使用BLAST把该序列输入到数据库中进行比对。BLAST会自动输出最相似的基因序列以及它们的相似得分和E值。如果找到一个和已知基因有很高相似性的序列，可能说明这个新基因在功能上也与已知基因类似。

总结

BLAST是一种快速、有效的生物序列比对工具，帮助科学家在海量的生物数据库中找到相似的序列。它通过分段比对、局部扩展和统计显著性过滤来保证结果的准确性和生物学意义，广泛应用于基因功能预测、序列注释和分子进化研究中。

演示

使用BLAST工具进行序列比对其实很方便，我们可以通过访问NCBI（美国国家生物技术信息中心）的BLAST在线工具来进行操作。以下是具体步骤：

访问BLAST网站

打开浏览器并访问NCBI的BLAST工具：https://blast.ncbi.nlm.nih.gov/Blast.cgi

选择BLAST工具类型

根据你的查询序列类型选择合适的BLAST工具：

Nucleotide blastn：用于DNA序列（核苷酸）比对。

Protein blast：用于蛋白质序列比对。

blastx：将DNA序列翻译成蛋白质后比对到蛋白质数据库。

tblastn: 将蛋白质序列比对到DNA序列数据库中。

假设你要比对DNA序列，选择“Nucleotide BLAST (blastn)”。

输入查询序列

在页面的“Enter Query Sequence”文本框中输入你要比对的DNA序列。例如，输入以下示例序列（人类的β-珠蛋白基因序列）：

复制代码

ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATGAAGTTGGTGGTGAGGCCCTGGGCAG

你可以直接粘贴序列，或上传包含序列的文本文件。

选择数据库

在“Database”部分中选择目标数据库，例如：

nt (Nucleotide collection)：用于与所有核酸序列数据库进行比对。

refseq_rna：仅限于RefSeq核酸序列库。

如果不确定，可以选择默认的“nt”数据库。

调整比对参数（可选）

通常，默认设置已经适用于大多数比对任务。如果你有特定需求，比如限制比对结果数量或调整E值阈值，可以在“Algorithm parameters”中做进一步调整。

运行BLAST比对

检查设置后，点击“BLAST”按钮开始比对。BLAST工具会对查询序列和数据库中的序列进行比对，这个过程可能需要几秒到几分钟。

查看比对结果

BLAST运行完成后，会显示比对结果页面。该页面包含以下几个关键部分：

Graphic Summary：展示查询序列与比对结果的相似区域。

Descriptions：列出与查询序列最相似的序列，按相似度排序。

Alignment：显示查询序列与比对结果的详细比对，包括匹配、错配和缺失位置。

分析比对结果

注意比对结果中的E值（E-value）：数值越小，匹配结果的显著性越高。

根据具体需求，查看最相似的比对序列，了解它们的功能、物种信息

参考资料：生信入门：序列比对之ncbi_blast在线使用_blast序列比对-CSDN博客

参考资料： 【NCBI-BLAST】-- NCBI—Blast 工具介绍 I 用序列查基因 I 目标序列序列比对 I 比对结果分析 I 比对序列下载---时间较长，谨慎食用_哔哩哔哩_bilibili

4.2 Needleman-Wunsch算法

  Needleman-Wunsch算法是一种全局序列比对算法，用于寻找两个序列之间的最佳全局对齐。它的目标是最大化两个序列整体的相似性，通过动态规划方法实现，广泛用于DNA、RNA或蛋白质序列的比对。

Needleman-Wunsch算法的步骤

1. 初始化矩阵

设定一个矩阵M，行数和列数分别对应两条序列的长度，并在矩阵第一行和第一列按顺序填入初始值。这些初始值代表了序列缺失的惩罚分数。

2. 填充矩阵

使用一个打分系统来定义匹配、错配和缺失（gap）得分。常用的打分方案是：

匹配（match）：+1分

错配（mis-match）：-1分

缺失（gap）：-2分

从矩阵的左上角（第一行和第一列初始化后的位置）开始，根据这三种情况，选择最高得分填入当前格子：

从左上方对角线移动过来（即匹配或错配）。

从上方格子移动过来（在序列1中加入gap）。

从左方格子移动过来（在序列2中加入gap）。

3. 矩阵回溯，寻找最佳对齐

在矩阵填充完成后，从右下角开始回溯，沿得分最高的路径向上左移动，直到回溯到左上角。这个路径将生成两个序列的最佳全局对齐。

填充Needleman-Wunsch算法的矩阵是关键的一步，它通过动态规划的方式计算每个位置的得分。下面是如何填充矩阵的详细步骤。

1. 初始化矩阵

首先，创建一个矩阵，其中行和列分别对应两个要比对的序列。假设我们有两个序列：

序列1：GATTACA

序列2：GCATGCU

矩阵的大小是(m+1) x (n+1)，其中m是序列1的长度，n是序列2的长度。因为矩阵的第一行和第一列用于表示缺失（gap），所以需要初始化它们。

矩阵大小：7行8列（分别对应GATTACA和GCATGCU）

初始化的矩阵（第一行和第一列）：

第一列表示与序列1的每个字符比对时的缺失（gap）惩罚。

第一行表示与序列2的每个字符比对时的缺失（gap）惩罚。

2. 填充矩阵

接下来，填充剩余的矩阵。对于矩阵中每一个位置(i, j)，我们需要考虑以下三种可能的情况，并选择最高的得分填入该位置：

匹配或错配（来自对角线方向）

如果序列1中的字符和序列2中的字符匹配，那么我们加上匹配的分数（通常是+1）。如果不匹配，则加上错配的分数（通常是-1）。

插入gap（来自左侧或上方）

如果我们在序列1或序列2中插入一个gap，则得分为-1，左侧或上方的格子的值加上这个gap的惩罚。

填充规则：

对于每个格子(i, j)，我们计算以下三种得分：

对角线：如果序列1[i-1] == 序列2[j-1]，则是匹配，得分是得分(i-1, j-1) + match_score，否则是错配，得分是得分(i-1, j-1) + mismatch_penalty。

上方：得分是得分(i-1, j) + gap_penalty。

左方：得分是得分(i, j-1) + gap_penalty。

然后，选择这三者中的最大值作为当前格子的得分。

填充示例

假设我们使用以下得分规则：

匹配：+1

错配：-1

缺失（gap）：-1

我们从i=1, j=1开始填充矩阵。对于每个位置(i, j)，我们将计算三个方向的得分并选择最大值。

矩阵填充过程：

(i=1, j=1)：比对G和G（匹配），得分是0 + 1 = 1。所以矩阵第一行第一列填入1。

(i=1, j=2)：比对G和C（错配），得分是0 - 1 = -1。由于前面已经是缺失的情况，所以填入-1。

(i=1, j=3)：比对G和A（错配），得分是-1 - 1 = -2。填入-2。

继续按照此方式填充矩阵……

填充完成后的矩阵大致如下

3. 回溯找到最佳对齐

填充矩阵完成后，我们从矩阵的右下角开始回溯，根据得分决定回溯的方向。通过回溯，可以得到最佳的全局对齐。

例如：

从得分(7, 7)开始回溯，根据得分决定是否是匹配、错配或者插入gap。

如果当前位置的得分来自对角线，表示该位置是匹配（或错配），如果来自左侧或上侧，表示有gap。

通过这个回溯过程，我们可以得到两个序列的最佳对齐。

总结

初始化矩阵：设置第一行和第一列为gap惩罚值。

填充矩阵：计算每个位置的得分，考虑匹配、错配和gap。

回溯：从矩阵右下角回溯，获得最佳全局对齐。

参考资料：生物信息学（1）——双序列比对之Needleman-Wunsch(NW)算法详解及C++实现-CSDN博客

参考资料：needleman-wuntch算法_哔哩哔哩_bilibili

4.3 Smith-Waterman算法

Smith-Waterman算法是一种用于局部序列比对的动态规划算法，与Needleman-Wunsch算法（全局比对）不同，它的目标是找到两个序列中相似的局部区域，而不仅仅是全序列的最佳对齐。这使得Smith-Waterman算法特别适用于那些在序列中有局部相似性，但两者整体不相似的情况。

Smith-Waterman算法的步骤

1. 初始化矩阵

首先，创建一个矩阵，其中行和列分别对应两个要比对的序列。矩阵的大小是(m+1) x (n+1)，m是序列1的长度，n是序列2的长度。初始化时，整个矩阵的第一行和第一列都是0，与Needleman-Wunsch不同，Smith-Waterman不允许负值得分。

2. 填充矩阵

与Needleman-Wunsch类似，在填充过程中，我们为每个矩阵位置(i, j)计算三个可能的得分，并选择最大值：

匹配或错配（来自对角线方向）

插入gap（来自左侧或上方）

0值（与第一行或第一列相比较，Smith-Waterman允许回到“0”点，代表局部比对的结束）

填充规则：

匹配：+1

错配：-1

缺失（gap）：-1

对于每个位置(i, j)，计算以下三个值：

来自对角线的得分：得分(i-1, j-1) + match_score（匹配或错配）

来自上方的得分：得分(i-1, j) + gap_penalty

来自左方的得分：得分(i, j-1) + gap_penalty

最后选择这三个中的最大值和0之间的较大者，填入当前格子。

3. 回溯找到局部最佳对齐

填充完成后，从矩阵中的最大得分开始回溯。不同于全局比对，Smith-Waterman算法的回溯过程会从得分最大的位置开始，一直到得分为零的地方停止。这种方法保证了找到最有意义的局部比对。

填充矩阵是Smith-Waterman算法中的一个关键步骤，用来计算两个序列的局部比对得分。填充矩阵的过程遵循动态规划的思想，目的是通过递归地比较每个位置的字符，并选择最佳的对齐得分。

填充矩阵的步骤

我们将详细介绍填充矩阵的过程，以帮助你理解每个细节。

假设我们有两个序列：

序列1：GATTACA

序列2：GCATGCU

我们需要构建一个矩阵，行表示序列1的字符，列表示序列2的字符。

1. 初始化矩阵

首先，我们初始化一个(m+1) x (n+1)的矩阵，其中m是序列1的长度，n是序列2的长度。矩阵的第一行和第一列都设为0，因为我们不考虑缺失字符时的得分。

对于序列1和序列2的例子：

序列1：GATTACA长度为7

序列2：GCATGCU长度为7

因此，我们需要构建一个8 x 8的矩阵（多加一行一列用于处理gap）。

初始化矩阵如下（第一行和第一列全为0）：

2. 填充矩阵

然后，我们从(i=1, j=1)开始逐个填充矩阵。每个矩阵的元素由以下三个值中的最大值决定：

匹配或错配：如果序列1[i-1]和序列2[j-1]相同，则得分为匹配得分（比如 +1）；如果不同，则得分为错配得分（比如 -1）。

缺失（gap）：如果某个位置是由插入gap得到的（即来自矩阵的左侧或上方），则得分会减少（比如 -1）。

零值：如果得分计算后为负数，则我们选择将该位置的得分设为零，表示局部比对的结束。

假设：

匹配的得分为 +1

错配的得分为 -1

gap的得分为 -1

我们逐个填充矩阵：

填充(i=1, j=1)：

比较序列1[0] = G和序列2[0] = G，它们匹配，得分为 +1。

因此，填充位置(1, 1)为1。

填充(i=1, j=2)：

比较序列1[0] = G和序列2[1] = C，它们不匹配，得分为 -1。

来自上方的gap得分为 -1，来自左侧的gap得分为 -1。

最大得分为 0，所以填充位置(1, 2)为0。

依此类推，我们填充整个矩阵，直到矩阵的最后一个位置。

最终矩阵（部分填充结果）：

3. 填充规则总结：

每个位置(i, j)的得分是由以下三者计算得到的：

对角线得分（即得分(i-1, j-1)加上匹配或错配得分）

上方得分（即得分(i-1, j)加上gap得分）

左方得分（即得分(i, j-1)加上gap得分）

得分为负时，选择0。

4. 完成矩阵后，使用回溯来得到最佳局部比对。

参考资料：striped smith waterman 动态规划算法SIMD实现讲解_哔哩哔哩_bilibili

参考资料：Smith-Waterman算法_smith waterman算法推导csdn-CSDN博客

5.关于如何回溯的问题

我们看下Chat-GPT是如何回答的

回溯是算法中的一个关键步骤，用来从矩阵中找到最佳的局部比对。回溯过程是通过从矩阵中找到最高得分的点开始，沿着得分路径追踪，直到得分变为零为止。这有助于确定两个序列之间的最佳局部对齐区域。

为了让你更清楚地理解回溯过程，我们可以通过具体的例子来演示。

1.我们有以下两个序列：

序列1：GATTACA

序列2：GCATGCU

2.通过前面的步骤，我们已经计算并填充了Smith-Waterman的矩阵，矩阵大致如下：

3.回溯

在这个矩阵中，最高得分是5，出现在(7, 7)位置。这个得分代表了序列1和序列2之间局部对齐的最佳结果。

3. 1从最大得分开始

我们从得分为5的位置(7, 7)开始回溯。接下来，根据得分路径的规则，决定如何回溯。

规则：

如果得分来自对角线(i-1, j-1)，表示两个字符是匹配（或者是错配），我们就将这两个字符加入到对齐中。

如果得分来自左侧(i, j-1)或上方(i-1, j)，表示有一个gap（缺失），我们将gap插入到其中一个序列中。

3.2. 回溯的具体过程

当前位置(7, 7)：得分是5，来自对角线，表示A和U之间是匹配（或者错配）。所以，我们将A对齐到U。

当前位置(6, 6)：得分是3，来自对角线，表示C和C是匹配。我们将C对齐到C。

当前位置(5, 5)：得分是2，来自对角线，表示A和G是错配。虽然它们不是完全匹配，但这个得分仍然来自于对角线，所以我们将A对齐到G。

当前位置(4, 4)：得分是1，来自对角线，表示T和T是匹配。我们将T对齐到T。

当前位置(3, 3)：得分是2，来自对角线，表示T和A是错配。

当前位置(2, 2)：得分是1，来自对角线，表示A和C是错配。

3.3. 回溯停止

当得分为零时，表示我们已经找到了最好的局部对齐。此时回溯结束。

3.4. 结果

通过回溯过程，我们得到以下对齐：

总结

回溯是从得分矩阵中得分最高的位置开始的。

根据得分的来源（对角线、上方、左方），我们决定是字符匹配/错配还是插入gap。

回溯直到遇到得分为零的地方，标志着局部比对的结束。

最终得到两个序列的最佳局部对齐结果。