生物信息学导论-北大-序列比对-序列数据库搜索

最近又重拾起了网上的学习资源——Coursera平台。近期决定投入足够的时间和精力来系统地学习这门课程，并计划在下次考试中取得优异成绩。因此也顺带记录下这段充实的学习经历。

ref: https://www.coursera.org/learn/sheng-wu-xin-xi-xue/home

本文主要来自本课的讲义。

---

背景

操作 ：拿查询序列去比对数据库中每一条序列

难题：若采用前面所述的对比方法，则执行完整个数据库中的所有查询序列将需要大量时间；因此亟需开发高效快速的算法来解决这一问题。

BLAST(Basic Local Alignment Search Tool) is a heuristic search method designed to identify local alignment solutions that are highly efficient for sequence comparison tasks.

启发式 ：not best but good enough。

BLAST尽管比普通的动态局部匹配快于1000倍以上，在处理相距较远的生物序列（如e.coli与human）时其灵敏度较低（low sensitivity）。

可用项目

programs	query seq	db type
nucleotide blast	nucleotide	nucleotide
blastn	nucleotide	nucleotide
protein blast	protein	protein
blastx	translated nucleotide	protein
tblastn	protein	translated nucleotide
tblastx	translated nucleotide	translated nucleotide
blastp	amino acid	protein

算法介绍

核心思路

blast algorithm outline:

• search for similarities (seed) between the query and subject.
  • expand these seeds into High Scoring Segment Pairs (HSPs) and conduct a localized Smith-Waterman alignment on the specified region.
  • Evaluate how reliable this alignment is.

细节

procedures:

1. 去除

2. 通过某种方法去除统计上显著但实际意义不大的数据

3. 该过滤操作涉及去除低复杂度序列、重复序列以及查询序列自身

4. 使用N或X来代替碱基对的位置

5. 种子词

6. 搜索结果

7. 不允许存在gap, 实行全匹配

8. 使用BLOSUM/PAM矩阵进行氨基酸配对, 核苷酸配对则采用+5/-4或+2/-3分值体系

9. 根据评分矩阵生成邻近词: 对于一个种子词, 从单个碱基替换开始计算得分为基础到全部碱基替换为止, 可以得到一系列与之长度相同的字符串序列; 如果这些序列中任何一个的得分为高于或等于阈值T, 则该序列为该种子词的邻近词之一

10. scan operation

11. 哈希表（HashTable）采用基于键值的直接寻址技术

12. 确定性有限自动机（Deterministic finite automaton/finite state machine, DFA/FSM）实现过程具有更高的效率

将该区域向两端延伸，在找到的位置上通过打分矩阵计算出其分数超过阈值S后，则可获得HSP序列

6. significance evaluation

相关概念

mask low-complexity seqs

低复杂度的序列会出现假阳性。

K=1LlogN(L!∏ini!) K=\frac{1}{L}log_N\left(\frac{L!}{\prod_{i}{n_i!}} \right)

其中L表示窗口长度（单位：碱基对），N代表字母表的数量（取值范围为4或20），而nin_i则表示第i个字符在序列中的发生频率（例如，在核酸序列中各个碱基对A、T、C、G出现的概率）。

举例：序列CACACACACACACACA，L=6

K=16log4(6!nA!∗nC!∗nG!∗nT!)=0.36 K=\frac{1}{6}log_4\left(\frac{6!}{n_A! * n_C! * n_G! * n_T!}\right)=0.36

seeding

给定一个query sequence MVLSPADKTNVKAAW（其长度为n），可以将其划分为多个连续的、长度均为w的seed word块。其中，在蛋白质序列中通常取w=3，在核酸序列中则取w=11。该序列可分割成 n - w + 1个连续且不重叠的 seed word块。

MVL, VLS, LSP, …, AAW

index db

在数据库中的序列，在预处理阶段为每个已知的 seed 建立索引后，则使得比对 seed 的过程更加高效。

neighbourhood words

在种子传播的过程中除了核心词汇之外 在邻域中还通过替换矩阵获取了与之非常接近的邻域词

E-value

一个匹配随机出现的可能性。

E=10在设定一个特定的匹配得分时代表可能随机产生或达到这一得分的有10个匹配

E=kmne−λS E=kmne^{-λS}

p=1−e−E p=1-e^{-E}

其中k值由搜索空间决定，在此框架中,m代表查询序列的长度,n表示数据库中的序列数量,λ值由评分矩阵确定,S为最高相似度对（HSP）的分数

E >1表示这个比对偶尔发生，E<0.1或0.05，表示统计上显著，E<10−510^{-5}表示高度相似