blast简介及单机版.docx

简介Blast，全称Basic Local Alignment Search Tool，即基于局部比对算法的搜索工具，由Altschul等人于1990年发布。Blast能够实现比较两段核酸或者蛋白序列之间的同源性的功能，它能够快速的找到两段序列之间的同源序列并对比对区域进行打分以确定同源性的高低。Blast的运行方式是先用目标序列建数据库（这种数据库称为database，里面的每一条序列称为subject），然后用待查的序列（称为query）在database中搜索，每一条query与database中的每一条subject都要进行双序列比对，从而得出全部比对结果。Blast是一个集成的程序包，通过调用不同的比对模块，blast实现了五种可能的序列比对方式：blastp：蛋白序列与蛋白库做比对，直接比对蛋白序列的同源性。blastx：核酸序列对蛋白库的比对，先将核酸序列翻译成蛋白序列（根据相位可以翻译为6种可能的蛋白序列），然后再与蛋白库做比对。blastn：核酸序列对核酸库的比对，直接比较核酸序列的同源性。tblastn：蛋白序列对核酸库的比对，将库中的核酸翻译成蛋白序列，然后进行比对。tblastx：核酸序列对核酸库在蛋白级别的比对，将库和待查序列都翻译成蛋白序列，然后对蛋白序列进行比对。Blast提供了核酸和蛋白序列之间所有可能的比对方式，同时具有较快的比对速度和较高的比对精度，因此在常规双序列比对分析中应用最为广泛。可以毫不夸张的说，blast是做比较基因组学乃至整个生物信息学研究所必须掌握的一种比对工具。下载NCBI提供免费下载，网址：/blast/executables/release/，可根据自己得机器选择相应操作系统的版本。安装直接解压缩包即可。解压缩命令：zcat *.tar.gz | tar xvf -使用Blast的运行分为两个步骤：第一，建立目标序列的数据库；第二，做blast比对。1运行建库程序formatdb：建库的过程是建立目标序列的索引文件，所用程序是formatdb。程序允许的输入格式FASTA或者ASN.1格式，通常我们使用FASTA格式的序列作为输入。用于建库的FASTA序列是db.seq，formatdb的基本命令是：formatdb -i db.seq -options常用的参数有以下几个：-p (T/F)：-p参数的意义是选择建库的类型，T表示蛋白库，F表示核酸库。缺省值为T。-o (T/F)：-o参数的意义是判断是否分析序列名并建立序列名索引。T表示建立序列名索引，F表示不建立序列名索引。缺省值为F。程序输出：如果建立的是核酸库，输出为db.seq.nhr、db.seq.nin、db.seq.nsq，如果选择了参数-o T，还会同时输出db.seq.nsd、db.seq.nsi、db.seq.nni、db.seq.nnd。蛋白库和核酸库的输出类似，相应的输出文件为：db.seq.phr、db.seq.pin、db.seq.psq和db.seq.psd、db.seq.psi、db.seq.pni、db.seq.pnd。除了这些结果，程序还会输出LOG文件（默认为formatdb.log），里面记录了运行时间、版本号、序列数量等信息。几点需要注意的问题：1、建库以后，做blast比对的输入文件就是建库所得的文件db.seq.n*或者db.seq.p*，而不是原始的FASTA序列。也就是说，建库以后，原始的序列文件是可以删除的。2、如果命令行中选择了-o T，并且目标序列中含有gi号重复的的序列名时，程序会停止建库并报错。例如，下列序列文件中出现了重复的序列名：gi|112385745|gb|DQ859020.1| Oryza sativa (japonica cultivar-group) glutathione S-transferase 2 mRNA, complete cdsATGGCGGAGGCGGCGGGGGCGGCGGTGGCGCCGGCGAAGCTGGGTCTGTACTCGTACTGGCGGAGCTCGTGCTCGCACCGCGTCCGCATCGCCCTCAACCTCAAAGGATTGGAGTACGAGTACAAGGCGGTGAACCTGCTCAAGGGGGAGCACTCTGATCCAGAATTCATGAAGGTTAATCCTATGAAGTTCGTCCCGGCATTGGTCGAT.CAAGCAGCACTCCCAGACAGACAACCAGATGCCCCTTCCTCTACCTAGgi|112385745|gb|DQ859020.1| Oryza sativa (japonica cultivar-group) glutathione S-transferase 2 mRNA, complete cdsATGGCGGAGGCGGCGGGGGCGGCGGTGGCGCCGGCGAAGCTGGGTCTGTACTCGTACTGGCGGAGCTCGTGCTCGCACCGCGTCCGCATCGCCCTCAACCTCAAAGGATTGGAGTACGAGTACAAGGCGGTGAACCTGCTCAAGGGGGAGCACTCTGATCCAGAATTCATGAAGGTTAATCCTATGAAGTTCGTCCCGGCATTGGTCGAT.运行时就会报如下错误：formatdb ERROR: Failed to create index. Possibly a gi included more than once in the database.3、如果输入序列不符合FASTA格式或者ASN.1格式，程序会自动退出，并报错：formatdb ERROR: Could not open db4、核酸序列可以用于建核酸库和蛋白库，但是蛋白序列不能用于建核酸库。其他参数简介：-l：-l 文件名用来改变LOG文件的命名-n：-n 文件名可以自定义生成的库文件命名-a：输入文件为ASN.1格式2运行比对程序blastall：Blast的主程序是blastall。程序的输入文件是query序列（-i 参数）和库文件（-d 参数），比对类型的选择（-p 参数）和输出文件（-o 参数）由用户指定。其中“-p”参数有5种取值：-p blastp：蛋白序列与蛋白库做比对。-p blastx：核酸序列对蛋白库的比对。-p blastn：核酸序列对核酸库的比对。-p tblastn：蛋白序列对核酸库的比对。-p tblastx：核酸序列对核酸库在蛋白级别的比对。这些元素就构成了blast的基本运行命令（以blastn为例）：blastall -i query.fasta -d database_prefix -o blast.out -p blastn其中如果-o参数缺省，则结果输出方式为屏幕输出。下面以一个blastn比对为例，来说明比对全过程：Query序列（query.fasta）：gi|45593933|gb|AY551259.1| Oryza sativa precursor microRNA 319c geneAGGAAGAGGAGCTCCTTTCGATCCAATTCAGGAGAGGAAGTGGTAGGATGCAGCTGCCGATTCATGGATACCTCTGGAGTGCATGGCAGCAATGCTGTAGGCCTGCACTTGCATGGGTTTGCATGACCCGGGAGATGAACCCACCATTGTCTTCCTCTATTGATTGGATTGAAGGGAGCTCCACATCTCTgi|45593932|gb|AY551258.1| Oryza sativa precursor microRNA 319b geneCATATTCTTTTAATTTGATGGAAGAAGCGATCGATGGATGGAAGAGAGCGTCCTTCAGTCCACTCATGGGCGGTGCTAGGGTCGAATTAGCTGCCGACTCATTCACCCACATGCCAAGCAAGAAACGCTTGAGATAGCGAAGCTTAGCAGATGAGTGAATGAAGCGGGAGGTAACGTTCCGATCTCGCGCCGTCTTTGCTTGGACTGAAGGGTGCTCCCTCCTCCTCGATCTCTTCGATCTAATTAAGCTACCTTGACAT库文件Database（db.seq，已经运行formatdb -i db.seq -p F -o T建库）：fake_seqAGGAAGAGGAGCTCCTTTCGTTCCAATTCAGGAGAGGAAGTGGTAGGATGCAGCTGCCGATTCATGGATACCTCTGGAGTGCATGCAGCAATGCTGTAGGCCTGCACTTGCATGGGTTTGCATGACCCGGCGAGATGAACCCACCATTGTCTTCCTCTATTGATTGGATTGAAGGGAGCTCCACATCTCT运行命令：blastall -i query.fasta -d db.seq -o blast.out -p blastn运行结果：BLASTN 2.2.8 Jan-05-2004Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),Gapped BLAST and PSI-BLAST: a new generation of protein database searchprograms, Nucleic Acids Res. 25:3389-3402.Query= gi|45593933|gb|AY551259.1| Oryza sativa precursor microRNA 319cgene, complete sequence (190 letters)Database: db.seq 1 sequences; 190 total lettersSearching.done Score ESequences producing significant alignments: (bits) Value fake_seq 339 2e-98fake_seq Length = 190Score = 339 bits (171), Expect = 2e-98 Identities = 188/191 (98%), Gaps = 2/191 (1%) Strand = Plus / PlusQuery: 1 aggaagaggagctcctttcgatccaattcaggagaggaagtggtaggatgcagctgccga 60 | |Sbjct: 1 aggaagaggagctcctttcgttccaattcaggagaggaagtggtaggatgcagctgccga 60Query: 61 ttcatggatacctctggagtgcatggcagcaatgctgtaggcctgcacttgcatgggttt 120 | |Sbjct: 61 ttcatggatacctctggagtgcat-gcagcaatgctgtaggcctgcacttgcatgggttt 119Query: 121 gcatgacccgg-gagatgaacccaccattgtcttcctctattgattggattgaagggagc 179 | |Sbjct: 120 gcatgacccggcgagatgaacccaccattgtcttcctctattgattggattgaagggagc 179Query: 180 tccacatctct 190 |Sbjct: 180 tccacatctct 190 Database: db.seq Posted date: Aug 28, 2006 8:14 PM Number of letters in database: 190 Number of sequences in database: 1Lambda K H 1.37 0.711 1.31GappedLambda K H 1.37 0.711 1.31Matrix: blastn matrix:1 -3Gap Penalties: Existence: 5, Extension: 2Number of Hits to DB: 3Number of Sequences: 1Number of extensions: 3Number of successful extensions: 3Number of sequences better than 10.0: 1Number of HSPs better than 10.0 without gapping: 1Number of HSPs successfully gapped in prelim test: 0Number of HSPs that attempted gapping in prelim test: 0Number of HSPs gapped (non-prelim): 1length of query: 190length of database: 190effective HSP length: 8effective length of query: 182effective length of database: 182effective search space: 33124effective search space used: 33124T: 0A: 0X1: 6 (11.9 bits)X2: 15 (29.7 bits)S1: 12 (24.3 bits)S2: 6 (12.4 bits)BLASTN 2.2.8 Jan-05-2004 Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),Gapped BLAST and PSI-BLAST: a new generation of protein database searchprograms, Nucleic Acids Res. 25:3389-3402.Query= gi|45593932|gb|AY551258.1| Oryza sativa precursor microRNA 319bgene, complete sequence (260 letters)Database: db.seq 1 sequences; 190 total lettersSearching.done* No hits found * Database: db.seq Posted date: Aug 28, 2006 8:14 PM Number of letters in database: 190 Number of sequences in database: 1Lambda K H 1.37 0.711 1.31GappedLambda K H 1.37 0.711 1.31Matrix: blastn matrix:1 -3Gap Penalties: Existence: 5, Extension: 2Number of Hits to DB: 0Number of Sequences: 1Number of extensions: 0Number of successful extensions: 0Number of sequences better than 10.0: 0Number of HSPs better than 10.0 without gapping: 0Number of HSPs successfully gapped in prelim test: 0Number of HSPs that attempted gapping in prelim test: 0Number of HSPs gapped (non-prelim): 0length of query: 260length of database: 190effective HSP length: 8effective length of query: 252effective length of database: 182effective search space: 45864effective search space used: 45864T: 0A: 0X1: 6 (11.9 bits)X2: 15 (29.7 bits)S1: 12 (24.3 bits)S2: 6 (12.4 bits) Blast的结果包含的信息很丰富。每一个query的比对结果从BLASTN开始，记录了版本和作者信息，Query= 之后记录了query名和序列长度。如果两条序列没有找到相关性信息，那么在Searching.done下方显示* No hits found *；反之，则在Searching.done下方记录了该query序列和库中每一条subject序列的比对概况列表，包括比对得分（Score）和期望值（E value）。期望值是一个大于0的正实数，代表两条序列不相关的可能性。期望值是在整体上综合评定两条序列的相似性的参数，期望值数值越小，序列相似性就越高，反之期望值数值越大，相似性越低。比对的输出结果会按照期望值从低到高的顺序来排列。Query序列和每一条subject序列比对结果的详细信息以开始。需要注意的是同一个query和同一个subject可能会有多个比对结果，每一个具体的结果从Score =开始，记录了比对得分、期望值、相似度百分比（identities）、比对的空位和两条序列的比对方向，之后是比对条形图，显示了比对区域内每个碱基的比对情况。列出两条序列的所有比对结果后，罗列比对的参数设置和统计信息，至此两条序列间的比对结果输出完毕。如上述结果所示， AY551259和fake_seq同为正链相比得到一个结果，期望值为2e-98，identities为98，中间有2个空位，两条序列相似度很高。而AY551258和fake_seq没有找到同源性。对于蛋白相关的比对，需要在blastall的运行目录下放置取代矩阵，并在运行时指定此替代矩阵，程序才能正常运行，否则blastall会报错退出。一般来讲，蛋白比对时最常用的取代矩阵是BLOSUM62矩阵。参数仅仅运行blast的基本运行命令，得到的结果往往不能清晰准确的表示出有用的信息。最大的问题就是有太多的冗余，很多很短的比对都会出现在输出结果中，导致结果杂乱无章。例如： Score ESequences producing significant alignments: (bits) ValueContig3421 out.ace.2 2367 0.0 Contig3424 out.ace.2 165 1e-40Contig3423 out.ace.2 30 4.9 Contig3314 out.ace.2 30 4.9 Contig3423 out.ace.2 Length = 148505Score = 30.2 bits (15), Expect = 4.9 Identities = 15/15 (100%) Strand = Plus / PlusQuery: 571 aaagaataaaattat 585 |Sbjct: 103697 aaagaataaaattat 103711可以很明显的看出，query序列和Contig3423的比对结果不能表示两条序列的相关性。事实上这个比对结果只是一个偶然出现的重复。这样的结果不但会浪费大量的运算和存储资源，更给结果分析带来了沉重的负担。为了处理杂乱无章的比对结果，满足各种比对需求，blast设置了很多参数来限制比对的范围和输出的形式。以下多数结果以blastn举例，如不做特殊说明，这些参数适用于所有比对方式。1-e参数：-e（value）参数是用来过滤比对较差的结果的，用-e参数指定一个实数，blast会过滤掉期望值大于这个数的比对结果。这样不但简化了结果，还缩短了运行时间和结果占用的空间。比如在上一个例子中，在命令行中加上限制期望值：blastall -i query.fasta -d db.seq -o blast.out -p blastn -e 1e-10那么结果中就会只剩下比对较好的结果： Score ESequences producing significant alignments: (bits) ValueContig3421 out.ace.2 2367 0.0 Contig3424 out.ace.2 165 1e-40通常，对于不同物种间的比对，期望值设在1e-5左右即可；而对于同源性较高的物种或者同种的比对，可以适当将期望值调得更小来过滤垃圾结果。比如同一物种cDNA和染色体的比对，参数可用1e-10或更高。2-F参数：-F（T/F）参数是用来屏蔽简单重复和低复杂度序列的。如果选T，程序在比对过程中会屏蔽掉query中的简单重复和低复杂度序列；选F则不会屏蔽。缺省值为T。例如，我们将如下含有两段简单重复的序列自己和自己进行比对（重复区用小写字母表示）：test1TACAATAAATAAAAAAGAGCTGTCTACAGTCTTTTcgcgcgcgcgcgTTCAGAAGTAAAGCACTATACAtttttttGTTTGTTCTTCTCAATTTAGGAAACTCAATGAACAATGAATACGAACTATTATTACCAGTAAATACAAGTAATAC第一次比对采用缺省参数：blastall -i test.seq -d test.seq -o test.blast -p blastn -e 1e-5得到的结果：test1 Length = 151Score = 186 bits (94), Expect = 1e-52 Identities = 132/151 (87%) Strand = Plus / PlusQuery: 1 tacaataaataaaaaagagctgtctacagtcttttnnnnnnnnnnnnttcagaagtaaag 60 | |Sbjct: 1 tacaataaataaaaaagagctgtctacagtcttttcgcgcgcgcgcgttcagaagtaaag 60Query: 61 cactatacannnnnnngtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120 | |Sbjct: 61 cactatacatttttttgtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120Query: 121 aactattattaccagtaaatacaagtaatac 151 |Sbjct: 121 aactattattaccagtaaatacaagtaatac 151第二次运行采用参数“-F F”：blastall -i test.seq -d test.seq -o test.blast -p blastn -e 1e-5 -F F得到的结果：test1 Length = 151Score = 299 bits (151), Expect = 1e-86 Identities = 151/151 (100%) Strand = Plus / PlusQuery: 1 tacaataaataaaaaagagctgtctacagtcttttcgcgcgcgcgcgttcagaagtaaag 60 |Sbjct: 1 tacaataaataaaaaagagctgtctacagtcttttcgcgcgcgcgcgttcagaagtaaag 60Query: 61 cactatacatttttttgtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120 |Sbjct: 61 cactatacatttttttgtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120Query: 121 aactattattaccagtaaatacaagtaatac 151 |Sbjct: 121 aactattattaccagtaaatacaagtaatac 151比较两个结果，我们看出使用缺省参数的比对结果损失了一部分信息，得到的统计结果也出现失真，期望值和identity都没有反映出真实情况。有时较长的重复序列甚至会导致比对终止。加了-F F就保证了比对结果的完整性。通常在大规模、低精度的比对中，往往用缺省参数，这样能避免程序把过多的时间浪费在无意义的简单重复上，提高运行速度；而在小规模、高精度的比对中，需要加上参数-F F，保证比对的精确度和完整性。3-m参数：“-e”参数能够做到筛选适当的比对结果，但是即使如此，blast的输出结果仍然非常庞大并且难以处理。为了精简输出、节省存储空间、实现更多功能并使结果易于处理，blast提供了参数“-m (integer)”来设定输出格式，可供选择的值为011之间的整数，缺省为0。下面就通过实例逐个解析“-m”参数能够实现的输出功能。输入文件的内容（针对-m 0到-m 7），其中：加粗的区域是三条序列的重合位置，注意subject1多一个碱基。query.fasta：query1TACAATAAATAAAATAGAGCTGTCTACAGTACTTTTTCAGGAACTCCTTCAGAAGTAAAGCACTATACAtttttttGTTTGTTCTTTTCAATTTAGGAAACTCAATGAACAATGAATACGAACTATTATTACCAGTAAATACAAGTAATACdatabase.fasta：subject1TCCTTCAGAAGTAAAGCACTATACAtttttttGTTTGTTCTTTTCAATTTAGGAAACTCAAATGAACAATGAATACsubject2AATTTAGGAAACTCAATGAACAATGAATACGAACTATTATTACCAGTAAATACAAGTAAT输出：-m 0：缺省参数，显示一个query和一个subject两两比对的信息。subject1 Length = 76Score = 93.7 bits (47), Expect = 1e-24 Identities = 68/76 (89%), Gaps = 1/76 (1%) Strand = Plus / PlusQuery: 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 103 | |Sbjct: 1 tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 60-m 1：显示query在所有subjects上的定位信息，并显示一致性比对信息，subject之间不同的碱基会被标出。Sequences producing significant alignments: (bits) Valuesubject2 119 2e-32subject1 94 1e-24QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactca 1041 1 . 150 1 .ttttttt. 61 | a-m 2：显示query在所有subjects上的定位信息但是不显示一致性比对信息，subject之间不同的碱基会被标出。Sequences producing significant alignments: (bits) Valuesubject2 119 2e-32subject1 94 1e-24QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactca 1041 1 aatttaggaaactca 150 1 tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 61 | a-m 3：显示query在所有subjects的定位和一致性比对信息，不显示subjects之间的差异。Sequences producing significant alignments: (bits) Valuesubject2 119 2e-32subject1 94 1e-24QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 1031 1 .- 140 1 .ttttttt.a 60-m 4：显示query在所有subjects上的定位信息但是不显示一致性比对信息，不显示subjects之间的差异。Sequences producing significant alignments: (bits) Valuesubject2 119 2e-32subject1 94 1e-24QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 1031 1 aatttaggaaactc- 140 1 tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 60-m 5：显示query在所有subjects上的定位信息但是不显示每个碱基的比对信息，补充-对齐比对区域，subjects之间不同的碱基会被标出。Sequences producing significant alignments: (bits) Valuesubject2 119 2e-32subject1 94 1e-24 QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactca 1041 1 -aatttaggaaactca 150 1 tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 61 | a-m 6：显示query在所有subjects上的定位信息但是不显示每个碱基的比对信息，补充-对齐比对区域，不显示subjects之间的差异。Sequences producing significant alignments: (bits) Valuesubject2 119 2e-32subject1 94 1e-24QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 1031 1 -aatttaggaaactc- 140 1 tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 60-m 7：输出XML格式的blast结果。 1 119.434 60 1.9564