Blast 比对

生信喵

一、序列比对
       序列比对是整个生物信息的核心，因为几乎每个生物信息分析过程都需要用到序列比对。判断两个基因或两段基因组片段是否相似是序列分析的基本工作。从序列数据库搜索，序列拼接到基因蛋白质功能注释，以及进化树构建等，都依赖于分子序列相似性的比较，也就是序列比对。
       序列比对的核心作用就是判断是否同源。所谓同源(homology)，是整个生物信息分析中最为重要的一个概念，只有同源的比对分析才是有意义的。同源也就是指来自于同一个祖先，两个物种从同一个祖先分化后，与不同的环境发生相互作用，其相应的 DNA 序列将各自发生一些替换或者插入缺失突变，也就是说序列不在精确相同。比如同一个大肠杆菌的祖先，经过不同的时间和空间的差异累积，最终就分化成不同类型的大肠杆菌。不同样品之间是同源的关系。
       与同源概念想对应的是序列相似，相似性(similarity)和同源性(homology)是两个完全不同的概念。相似性仅仅是指字符串的相似 ，并不具有不具有生物学意义 ，因为 DNA 序列一共就有 ATCG 四种碱基，由于组合造成两段片段字符串组合比较接近。同源序列一般是相似的，但是相似的序列不一定同源。那么该如何判断序列是相似还是不相似，相似的序列是否满足同源关系呢。这些都需要序列比对来判断，并且使用一些方法和标准来进行评价。通常的解决方法是将两条序列或者两条序列中的部分进行序列比对，然后基于得到的比对结果判断相似性是由于序列之间具有进化渊源，还是纯粹的随机巧合。
       两条序列比对碱基之间可能有以下几种情况。两个碱基类型完全一致，没有发生变化，另外就是碱基发生了变化，这种变化包括替换、插入和删除。插入和删除也被称为空位，我们平时在比对过程中得错配其实就是替换，gap 就是插入或者删除。因为突变是随机的，但是选择是具有偏向性的，这就使得某些突变发生的可能性远大于其他类型。根据比对过程中不同的情况给定分值，例如完全匹配加 1 分，错配减一分，空位减一分，那么最终两条序列比对我们就会计算得到一个分值。这个分值越高就说明完全匹配的碱基越多，比对的情况越好。这就是打分矩阵。比对后根据打分矩阵的值就能判断比对的情况。打分矩阵有很多种，那么最简单的就是这种等价矩阵表。但是实际上我们知道这样明显会有很多问题，因为 DNA 上四种碱基替换比率是不一样的，由于化学结构上的相似性，嘌呤与嘌呤之间，嘧啶与嘧啶之间更容易发生替换，这种称为转换，而嘌呤与嘧啶之间发生替换稍微难一些，这种称为颠换。所以，我们将打分矩阵修改一下。完全匹配加 2 分，转换扣一分，颠换扣两分，这个打分矩阵要比前面一个更加精确一些。
       不过这两种矩阵我们都还没有考虑到空位罚分，对于插入或缺失突变，目前为止，还没有适当的统计模型。通常用的是与替换矩阵相对应的经验分值，一般空位罚分有两个参数值，即起始空位罚分和空位延伸罚分。对于一个长度为 k 的连续空位，罚分可以表示为 score=a+b*k，其中 a 为起始空位罚分，b 为空位延伸罚分。
       对于序列比对罚分矩阵有很多，但是没有一种矩阵能够适用于所有的情况，都会存在一定的误差，对于 DNA 序列比对，主要有 mat50 和 mat70 两种矩阵，对于氨基酸序列比对，主要采用 PAM 和 BLOSUM 两种矩阵。blast 比对中默认使用的就是 BLOSUM62 打分矩阵。其中 62 表示用来构建该矩阵的匹配数据集中精确匹配位点要占 62%。
       此外，还有转移矩阵，等价矩阵，遗传密码矩阵，疏水性矩阵等。

二、全局比对与局部比对
       全局比对是用来衡量两条序列整体的相似性，满足整体相似性最大化。若两条序列长度不同，则必须插入一些空位使所有位点都能对应起来。而局部比对则不同，两条亲缘关系较远的DNA 或氨基酸可能只在一些片段上相似，这就需要找到这些相似性的片段，和其相应的匹配方式。通常这样的分析就需要进行局部比对，而不是全局比对。
       全局比对与局部比对有什么不同呢。全局序列比对尝试找到两个完整的序列之间的最佳比对。而局部序列比对不必对两个完整的序列进行比对；可以在每个序列中使用某些部分来获得最大得分。两种比对采取不同的比对算法和策略，因此，同样的一段序列，采用全局比对和局部比对不同的比对方法结果也会有很大的不同。

      
       全局比对与局部比对
       例如我们现在有两条序列 S1 和 S2，如果采用全局比对，会得到这种比对效果，而采用局部比对，序列中间的 GCG 满足了最优比对。大家可以理解为，全局比对需要从全局出发，是需要全局达到最佳效果，而局部比对则不需要考虑全局，只要局部达到最佳效果即可。
       全局比对主要用来比较比较两个基因组之间的同源性，绘制共线性图等，另外，全局比对也常常用于基因组结构变异的检测。因为，局部比对的话，遇到大的空位往往就断开了，例如上面的例子，采用局部比对的算法中，只追求局部的最优比对，而不会考虑整体的空位等。所以，基因组的大片段的插入或者缺失检测，可以使用全局比对软件。而局部比对软件主要搜索同源序列，例如判断那两个基因是否同源，寻找一段序列的同源序列等，就可以使用局部比对。


三、blast 简介
       Blast 全称 Basic Local Alignment Search Tool，即"基于局部比对算法的搜索工具"，。1989 年，美国国家生物技术信息中心（NCBI）首次推出 BLAST。自第一版以来，NCBI 一直在维护和更新 BLAST 版本。本地 BLAST 不受网速、序列数目的限制，能够快速、准确的进行序列比对。
       Blast 能够实现比较两段核酸或者蛋白序列之间的同源性的功能，它能够快速的找到两段序列之间的同源序列并对比对区域进行打分以确定同源性的高低。因为是局部比对，所以只要序列之间出现同源区域就可以，而不用考虑整体，因此，blast 比对结果就会出现很多多对多的比对。也容易出现很多较差的比对，一个基因与另一个基因分成多份比对结果。2009年 7 月，NCBI 发布了 BLAST 升级版——BLAST+，BLAST+使用了 BLAST 的核心算法，延续了 BLAST 的优势功能，发展并增强了如 BLAST 的 fastacmd 程序，新增了如 update_blastdb.pl等程序。BLAST+是一个全新设计的程序，在性能以及易用性上均有了很大提高，清晰的模块化将会极大的提高维护者的效率。NCBI 官方也强烈推荐用户放弃 BLAST，使用 BLAST+工具。
       下载地址：https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
       Manual 使用手册：https://www.ncbi.nlm.nih.gov/books/NBK279690/

1. #软件安装
   
2. mamba install -y blast

复制代码

      blast 应用程序

软件名	功能
blastdbcmd	从索引中提取信息
blastn	核酸与核酸比对
blastp	氨基酸与氨基酸比对
blastx	核酸与氨基酸比对
makeblastdb	建立索引
tblastn	氨基酸与核酸比对
tblastx	核酸与核酸翻译为氨基酸比对
update_blastdb.pl	下载blast数据库

四、blast 数据库
4.1 NCBI blast 数据库
       blast 比对需要建立索引，索引 index，是目录的意思。由于序列长度较长，索引文件可以快速定位到目标区域。索引文件可以从 NCBI 下载，也可以自己构建。
       ftp 地址：https://ftp.ncbi.nlm.nih.gov/blast/db/

4.2 blast 数据库下载

1. #下载 blast nt 数据库
   
2. for i in {00..50};do echo "~/.aspera/connect/bin/ascp -i
   
3. ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --overwrite=diff -QTr -l6000m
   
4. anonftp@ftp.ncbi.nlm.nih.gov:blast/db/nt.${i}.tar.gz ./ ";done;

复制代码

4.3 自己构建数据库
下载 nt 库序列

1. ~/.aspera/connect/bin/ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh
   
2. --overwrite=diff -QTr -l6000m
   
3. anonftp@ftp.ncbi.nlm.nih.gov:blast/db/FASTA/nt.gz ./

复制代码

自己构建数据库

1. gunzip nt.gz
   
2. makeblastdb -in nt -dbtype nucl -parse_seqids -out nt

复制代码

五、在线 blast
网址：https://blast.ncbi.nlm.nih.gov/Blast.cgi

计算资源有限。

生信喵

六、选项参数
blast 常用选项参数

选项	释义
-h	显示选项参数
-help	显示帮助文档
-db	比对数据库
-query	待比对序列
-out	输出文件名
-evalue	比对 e 值
-outfmt	输出文件格式
-task	比对类型
-num_threads	使用线程数
-subject	两两比对，目标序列
-remote	联网比对
-query_loc	设定 query 的起始和终止位点
-num_alignments	显示比对上的结果数目
-strand	比对方向

7 输出格式
7.1 blast 输出格式

1. 0 = Pairwise,
   
2. 1 = Query-anchored showing identities,
   
3. 2 = Query-anchored no identities,
   
4. 3 = Flat query-anchored showing identities,
   
5. 4 = Flat query-anchored no identities,
   
6. 5 = BLAST XML,
   
7. 6 = Tabular,
   
8. 7 = Tabular with comment lines,
   
9. 8 = Seqalign (Text ASN.1),
   
10. 9 = Seqalign (Binary ASN.1),
    
11. 10 = Comma-separated values,
    
12. 11 = BLAST archive (ASN.1),
    
13. 12 = Seqalign (JSON),
    
14. 13 = Multiple-file BLAST JSON,
    
15. 14 = Multiple-file BLAST XML2,
    
16. 15 = Single-file BLAST JSON,
    
17. 16 = Single-file BLAST XML2,
    
18. 17 = Sequence Alignment/Map (SAM),
    
19. 18 = Organism Report

复制代码

7.2 自定义格式
blast 默认的 18 中格式中，6,7,17 为列表格式，可以重新自定义，也包括格式 10。
格式为：

1. -outfmt "6 qseqid sseqid pident"

复制代码

具体输出关键字

1. qseqid means Query Seq-id
   
2. qgi means Query GI
   
3. qacc means Query accesion
   
4. qaccver means Query accesion.version
   
5. qlen means Query sequence length
   
6. sseqid means Subject Seq-id
   
7. sallseqid means All subject Seq-id(s), separated by a ';'
   
8. sgi means Subject GI
   
9. sallgi means All subject GIs
   
10. sacc means Subject accession
    
11. saccver means Subject accession.version
    
12. sallacc means All subject accessions
    
13. slen means Subject sequence length
    
14. qstart means Start of alignment in query
    
15. qend means End of alignment in query
    
16. sstart means Start of alignment in subject
    
17. send means End of alignment in subject
    
18. qseq means Aligned part of query sequence
    
19. sseq means Aligned part of subject sequence
    
20. evalue means Expect value
    
21. bitscore means Bit score
    
22. score means Raw score
    
23. length means Alignment length
    
24. pident means Percentage of identical matches
    
25. nident means Number of identical matches
    
26. mismatch means Number of mismatches
    
27. positive means Number of positive-scoring matches
    
28. gapopen means Number of gap openings
    
29. gaps means Total number of gaps
    
30. ppos means Percentage of positive-scoring matches
    
31. frames means Query and subject frames separated by a '/'
    
32. qframe means Query frame
    
33. sframe means Subject frame
    
34. btop means Blast traceback operations (BTOP)
    
35. staxid means Subject Taxonomy ID
    
36. ssciname means Subject Scientific Name
    
37. scomname means Subject Common Name
    
38. sblastname means Subject Blast Name
    
39. sskingdom means Subject Super Kingdom
    
40. staxids means unique Subject Taxonomy ID(s), separated by a ';'
    
41. (in numerical order)
    
42. sscinames means unique Subject Scientific Name(s), separated by a ';'
    
43. scomnames means unique Subject Common Name(s), separated by a ';'
    
44. sblastnames means unique Subject Blast Name(s), separated by a ';'
    
45. (in alphabetical order)
    
46. sskingdoms means unique Subject Super Kingdom(s), separated by a ';'
    
47. (in alphabetical order)
    
48. stitle means Subject Title
    
49. salltitles means All Subject Title(s), separated by a '<>'
    
50. sstrand means Subject Strand
    
51. qcovs means Query Coverage Per Subject
    
52. qcovhsp means Query Coverage Per HSP
    
53. qcovus means Query Coverage Per Unique Subject (blastn only)
    
54. When not provided, the default value is:
    
55. 'qaccver saccver pident length mismatch gapopen qstart qend sstart send
    
56. evalue bitscore', which is equivalent to the keyword 'std'
    
57. The supported format specifier for option 17 is:
    
58. SQ means Include Sequence Data
    
59. SR means Subject as Reference Seq

复制代码