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-merquery_kmer <- kmer(query_seq_num, k)db_kmer <- kmer(db_seq_num, k)# 对于每个查问序列的k-mer，查找数据库序列中所有与之匹配的k-mermatches <- 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 NCBIWWWfrom 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多平台公布