关于后端:生物学经典blast比对算法R语言和Python如何实现

Blast 比对算法原理与实现形式

做生物的同学必定据说过 blast 比对这个办法，个别在 NCBI 等网站上能够在线进行比对，也能够在本地服务器进行比对，那么 blast 算法到底是怎么实现对不同序列的比对呢？

本文分享经典 blast 算法的根底原理，以及通过 R 语言和 Python 实现这个算法，不依赖网站本人进行序列比对。

什么是 BLAST 比对？

BLAST（Basic Local Alignment Search Tool）是一种罕用的生物信息学算法，用于比对两个或多个序列。BLAST 通过寻找两个序列之间的最大匹配来确定它们之间的相似性。

算法原理

BLAST 算法的原理：
将查问序列与数据库中的序列进行比对，找到最佳匹配。

BLAST 算法的逻辑：首先将查问序列进行分段，而后将这些分段与数据库中的序列进行比对。

K-mer 小片段

在比对过程中，BLAST 算法应用一种称为 K -mer 的技术，将查问序列和数据库序列分成长度为 K 的小片段，而后将这些小片段进行比对。

如果两个小片段具备类似的序列，BLAST 算法就会将它们合并成更长的序列，以便进行更精确的比对。

特点与利用

BLAST 算法的长处是速度快、准确度高，能够在大型数据库中疾速查找类似序列。BLAST 算法在生物信息学畛域中被广泛应用，用于基因正文、蛋白质构造预测、序列比对等方面。

不同序列 blast 比拟算法

1. 将查问序列和数据库序列别离转换为碱基对应的数字编码，例如 A 示意为 1，C 示意为 2，G 示意为 3，T 示意为 4。
2. 将查问序列划分成长度为 k 的小片段，称为 k -mer。
3. 将数据库序列划分成长度为 k 的小片段，称为 k -mer。
4. 对于每个查问序列的 k -mer，查找数据库序列中所有与之匹配的 k -mer。
5. 对于每个匹配的 k -mer，计算查问序列和数据库序列之间的类似度得分。
6. 对于每个查问序列的 k -mer，抉择类似度得分最高的匹配序列，并将其作为最佳匹配。
7. 对于每个最佳匹配，计算匹配序列的长度、类似度得分、E 值等参数。
8. 依据 E 值和类似度得分，对匹配后果进行排序，输入最终的比对后果。

BLAST 算法的具体实现可能会有所不同，上述算法仅作为一个示例，理论利用中须要依据具体情况进行调整。

此外，BLAST 算法的计算复杂度较高，如果对于理论生物数据处理，须要应用高性能计算机或云计算平台进行计算。

R 语言中实现 blast 算法

以下是一个基于 R 语言的 BLAST 比对算法示例，用于比对两个 DNA 序列：

# 导入 Biostrings 包
library(Biostrings)

# 定义查问序列和数据库序列
query_seq <- DNAString("ATCGATCGATCGATCG")
db_seq <- DNAString("CGATCGATCGATCGATC")

# 定义 k -mer 的长度
k <- 3

# 将查问序列和数据库序列别离转换为数字编码
query_seq_num <- as.numeric(query_seq)
db_seq_num <- as.numeric(db_seq)

# 将查问序列和数据库序列别离划分成 k -mer
query_kmer <- kmer(query_seq_num, k)
db_kmer <- kmer(db_seq_num, k)

# 对于每个查问序列的 k -mer，查找数据库序列中所有与之匹配的 k -mer
matches <- matchPattern(query_kmer, db_kmer)

# 对于每个匹配的 k -mer，计算查问序列和数据库序列之间的类似度得分
scores <- pmatch(query_kmer, db_kmer, fixed=FALSE)

# 对于每个查问序列的 k -mer，抉择类似度得分最高的匹配序列，并将其作为最佳匹配
best_matches <- maxMatches(matches)

# 对于每个最佳匹配，计算匹配序列的长度、类似度得分、E 值等参数
match_length <- width(best_matches)
match_score <- scores[best_matches]
e_value <- length(db_kmer) * (1 - exp(-match_score))

# 依据 E 值和类似度得分，对匹配后果进行排序，输入最终的比对后果
result <- data.frame(query_seq, db_seq, match_length, match_score, e_value)
result <- result[order(result$e_value),]

Python 实现 blast 算法

首先，须要装置 Biopython 库来实现 BLAST 比对算法。您能够应用以下命令在终端中装置 Biopython：

pip install biopython

接下来，能够应用以下代码来实现 BLAST 比对算法：

from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

# 进行 BLAST 比对
result_handle = NCBIWWW.qblast("blastn", "nt", "ACGTGAGGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGC")

# 读取 BLAST 比对后果
blast_record = NCBIXML.read(result_handle)

# 输入比对后果
for alignment in blast_record.alignments:
    for hsp in alignment.hsps:
        print('****Alignment****')
        print('sequence:', alignment.title)
        print('length:', alignment.length)
        print('e value:', hsp.expect)
        print(hsp.query[0:75] + '...')
        print(hsp.match[0:75] + '...')
        print(hsp.sbjct[0:75] + '...')

这段代码会将序列 ”ACGTGAGGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGC” 与 NCBI 的 nt 数据库进行比对。

本文由 mdnice 多平台公布