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

最近重新捡起coursera上的课了，这次准备好好学，把考试考了。。因此顺便记录一下学习过程。

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

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

---

背景

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

问题：如果用前面的比对方法，那一条查询序列遍历完数据库中所有序列会耗费很长时间，因此需要一种非常快的算法。

BLAST(Basic Local Alignment Search Tool)：heuristic algorithm, locally optimal alignments, very fast

启发式：not best but good enough。

BLAST虽然比普通的动态局部匹配快1000倍以上，但是对远一些的序列（e.coli vs 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:

• find matches(seed) between the query and subject.
• extend seed into High Scoring Segment Pairs(HSPs), run Smith-Waterman algorithm（local） on the specified region only.
• Assess the reliability of the alignment

细节

procedures:

1. filter
   1. 为了减少统计上显著但实际不重要的结果
   2. 过滤的是低复杂度序列和重复序列，以及查询序列本身
   3. 用N或X取代核苷酸和氨基酸
2. seeding
3. search word hits
   1. 不要gap，完全匹配
   2. Scoring matrix：氨基酸使用BLOSUM/PAM，核苷酸match/mismatch使用+5/-4或+2/-3
   3. 根据打分矩阵生成neighbourhood words：对于一个seed word，从替换一个碱基后算得分开始到完全替换，能得到一串同长度的字符串，如果一个字符串的得分比阈值T高（或等于T），就是邻近词了
4. scanning
   1. HashTable：直接寻址
   2. Deterministic finite automaton/finite state machine：更快
5. extending
   1. 在找到的位置向两侧延长，根据打分矩阵计算的分数在cutoff score S以上，就得到了HSP
6. significance evaluation

相关概念

mask low-complexity seqs

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

$$K=\frac{1}{L}lo{g}_N\left(\frac{L!}{{\prod }_i{n}_i!}\right)$$

其中L是window长度，N是alphabet大小（4或20）,$n_i$是第i个字符出现的频率（比如核酸就算ACGT出现的频率）。

举例：序列CACACACACACACACA，L=6

$$K=\frac{1}{6}lo{g}_4\left(\frac{6!}{n_A!\ast n_C!\ast n_G!\ast n_T!}\right)=0.36$$

seeding

假设有条query sequence是MVLSPADKTNVKAAW（长度为n），可以把它分成很多个连续的、长度为w的seed words。对于蛋白序列，w一般取值3，对于核酸序列，w一般取11。这条序列可以切为如下 n - w + 1个seed words.

MVL, VLS, LSP, …, AAW

index db

对于数据库中的序列，提前为每个已知的seed建立索引，那么比对seed的时候会快一点

neighbourhood words

seeding的时候除了seed words，还用基于替换矩阵取得的、很相似的neighbourhood words

E-value

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

E=10表示在给定一个匹配分数的情况下，可能随机出现10个匹配能达到这个分数。

$$E=kmn{e}^{-\lambda S}$$

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

其中k取决于搜索空间，m是query sequence的长度，n是database的长度，λ取决于打分矩阵，S是HSP的得分

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