云之南

风声，雨声，读书声，声声入耳；家事，国事，天下事，事事关心

日志

关于我

云之南

专业背景：计算机科学研究方向与兴趣: JavaEE-Web软件开发, 生物信息学, 数据挖掘与机器学习, 智能信息系统目前工作: 基因组, 转录组, NGS高通量数据分析, 生物数据挖掘, 植物系统发育和比较进化基因组学

文章分类

Blast使用方法攻略(二)

2010-01-01 12:14:48| 分类：生物信息学 | 标签： |举报 |字号大中小订阅

下载LOFTER 我的照片书 |

参数

仅仅运行blast的基本运行命令，得到的结果往往不能清晰准确的表示出有用的信息。最大的问题就是有太多的冗余，很多很短的比对都会出现在输出结果中，导致结果杂乱无章。例如：

Score E

Sequences producing significant alignments: (bits) Value

Contig3421 out.ace.2 2367 0.0

Contig3424 out.ace.2 165 1e-40

Contig3423 out.ace.2 30 4.9

Contig3314 out.ace.2 30 4.9

……

>Contig3423 out.ace.2

Length = 148505

Score = 30.2 bits (15), Expect = 4.9

Identities = 15/15 (100%)

Strand = Plus / Plus

Query: 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 E

Sequences producing significant alignments: (bits) Value

Contig3421 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"。例如，我们将如下含有两段简单重复的序列自己和自己进行比对（重复区用小写字母表示）：

>test1

TACAATAAATAAAAAAGAGCTGTCTACAGTCTTTTcgcgcgcgcgcgTTCAGAAGTAAAG

CACTATACAtttttttGTTTGTTCTTCTCAATTTAGGAAACTCAATGAACAATGAATACG

AACTATTATTACCAGTAAATACAAGTAATAC

第一次比对采用缺省参数：

blastall -i test.seq -d test.seq -o test.blast -p blastn -e 1e-5

得到的结果：

>test1

Length = 151

Score = 186 bits (94), Expect = 1e-52

Identities = 132/151 (87%)

Strand = Plus / Plus

Query: 1 tacaataaataaaaaagagctgtctacagtcttttnnnnnnnnnnnnttcagaagtaaag 60

||||||||||||||||||||||||||||||||||| |||||||||||||

Sbjct: 1 tacaataaataaaaaagagctgtctacagtcttttcgcgcgcgcgcgttcagaagtaaag 60

Query: 61 cactatacannnnnnngtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120

||||||||| ||||||||||||||||||||||||||||||||||||||||||||

Sbjct: 61 cactatacatttttttgtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120

Query: 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 = 151

Score = 299 bits (151), Expect = 1e-86

Identities = 151/151 (100%)

Strand = Plus / Plus

Query: 1 tacaataaataaaaaagagctgtctacagtcttttcgcgcgcgcgcgttcagaagtaaag 60

||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

Sbjct: 1 tacaataaataaaaaagagctgtctacagtcttttcgcgcgcgcgcgttcagaagtaaag 60

Query: 61 cactatacatttttttgtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120

||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

Sbjct: 61 cactatacatttttttgtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120

Query: 121 aactattattaccagtaaatacaagtaatac 151

|||||||||||||||||||||||||||||||

Sbjct: 121 aactattattaccagtaaatacaagtaatac 151

比较两个结果，我们看出使用缺省参数的比对结果损失了一部分信息，得到的统计结果也出现失真，期望值和identity都没有反映出真实情况。有时较长的重复序列甚至会导致比对终止。加了"-F F"就保证了比对结果的完整性。通常在大规模、低精度的比对中，往往用缺省参数，这样能避免程序把过多的时间浪费在无意义的简单重复上，提高运行速度；而在小规模、高精度的比对中，需要加上参数"-F F"，保证比对的精确度和完整性。

3．-m参数：

“-e”参数能够做到筛选适当的比对结果，但是即使如此，blast的输出结果仍然非常庞大并且难以处理。为了精简输出、节省存储空间、实现更多功能并使结果易于处理，blast提供了参数“-m (integer)”来设定输出格式，可供选择的值为0～11之间的整数，缺省为0。下面就通过实例逐个解析“-m”参数能够实现的输出功能。

输入文件的内容（针对-m 0到-m 7），其中：加粗的区域是三条序列的重合位置，注意subject1多一个碱基。

query.fasta：

>query1

TACAATAAATAAAATAGAGCTGTCTACAGTACTTTTTCAGGAACTCCTTCAGAAGTAAAG

CACTATACAtttttttGTTTGTTCTTTTCAATTTAGGAAACTCAATGAACAATGAATACG

AACTATTATTACCAGTAAATACAAGTAATAC

database.fasta：

>subject1

TCCTTCAGAAGTAAAGCACTATACAtttttttGTTTGTTCTTTTCAATTTAGGAAACTCA

AATGAACAATGAATAC

>subject2

AATTTAGGAAACTCAATGAACAATGAATACGAACTATTATTACCAGTAAATACAAGTAAT

输出：

-m 0：缺省参数，显示一个query和一个subject两两比对的信息。

>subject1

Length = 76

Score = 93.7 bits (47), Expect = 1e-24

Identities = 68/76 (89%), Gaps = 1/76 (1%)

Strand = Plus / Plus

Query: 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 103

||||||||||||||||||||||||| ||||||||||||||||||||||||||

Sbjct: 1 tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 60

……

-m 1：显示query在所有subjects上的定位信息，并显示一致性比对信息，subject之间不同的碱基会被标出。

Sequences producing significant alignments: (bits) Value

subject2 119 2e-32

subject1 94 1e-24

QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactca 104

1 1 ............... 15

0 1 .........................ttttttt............................ 61

\

|

a

……

-m 2：显示query在所有subjects上的定位信息但是不显示一致性比对信息，subject之间不同的碱基会被标出。

Sequences producing significant alignments: (bits) Value

subject2 119 2e-32

subject1 94 1e-24

QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactca 104

1 1 aatttaggaaactca 15

0 1 tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 61

\

|

a

……

-m 3：显示query在所有subjects的定位和一致性比对信息，不显示subjects之间的差异。

Sequences producing significant alignments: (bits) Value

subject2 119 2e-32

subject1 94 1e-24

QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 103

1 1 ..............- 14

0 1 .........................ttttttt...........................a 60

……

-m 4：显示query在所有subjects上的定位信息但是不显示一致性比对信息，不显示subjects之间的差异。

Sequences producing significant alignments: (bits) Value

subject2 119 2e-32

subject1 94 1e-24

QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 103

1 1 aatttaggaaactc- 14

0 1 tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 60

……

-m 5：显示query在所有subjects上的定位信息但是不显示每个碱基的比对信息，补充"-"对齐比对区域，subjects之间不同的碱基会被标出。

Sequences producing significant alignments: (bits) Value

subject2 119 2e-32

subject1 94 1e-24

QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactca 104

1 1 ---------------------------------------------aatttaggaaactca 15

0 1 tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 61

\

|

a

……

-m 6：显示query在所有subjects上的定位信息但是不显示每个碱基的比对信息，补充"-"对齐比对区域，不显示subjects之间的差异。

Sequences producing significant alignments: (bits) Value

subject2 119 2e-32

subject1 94 1e-24

QUERY 45 tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 103

1 1 ---------------------------------------------aatttaggaaactc- 14

0 1 tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 60

……

-m 7：输出XML格式的blast结果。

<Hsp_num>1</Hsp_num>

<Hsp_bit-score>119.434</Hsp_bit-score>

<Hsp_score>60</Hsp_score>

<Hsp_ue>1.95644e-32</Hsp_ue>

<Hsp_query-from>90</Hsp_query-from>

<Hsp_query-to>149</Hsp_query-to>

<Hsp_hit-from>1</Hsp_hit-from>

<Hsp_hit-to>60</Hsp_hit-to>

<Hsp_query-frame>1</Hsp_query-frame>

<Hsp_hit-frame>1</Hsp_hit-frame>

<Hsp_identity>60</Hsp_identity>

<Hsp_positive>60</Hsp_positive>

<Hsp_align-len>60</Hsp_align-len>

<Hsp_qseq>AATTTAGGAAACTCAATGAACAATGAATACGAACTATTATTA</Hsp_qseq>

<Hsp_hseq>AATTTAGGAAACTCAATGAACAATGAATACGAACTATTATTA</Hsp_hseq>

<Hsp_midline>||||||||||||||||||||||||||||||||||||||||||</Hsp_midline>

……

-m 8：列表格式的比对结果。从左到右各列的意义依次是：query名、subject名、identity、比对长度、错配数、空位数、query比对起始坐标、query比对终止坐标、subject比对起始坐标、subject比对终止坐标、期望值、比对得分。

query1 sub24 91.11 45 3 1 198 241 502208 502252 2.7e-06 50.05

query1 sub21 98.68 151 2 0 532 682 1360665 1360515 1.0e-76 284.0

query1 sub21 86.17 94 12 1 198 290 479232 479139 4.8e-14 75.82

query1 sub21 87.04 54 7 0 238 291 1297867 1297920 6.9e-07 52.03

query2 sub21 99.44 892 3 2 28 918 1351055 1350165 0.0 1713.2

query2 sub21 87.58 153 17 1 343 495 1358110 1357960 2.1e-35 147.2

query2 sub21 84.11 107 16 1 699 805 1305723 1305618 4.0e-12 69.88

query2 sub21 89.58 48 5 0 519 566 1305968 1305921 6.0e-08 56.00

query2 sub14 88.24 153 16 1 343 495 145402 145252 8.7e-38 155.1

query2 sub24 88.08 151 16 1 345 495 567561 567709 1.4e-36 151.2

query2 sub24 87.80 123 14 1 686 808 563341 563220 1.9e-26 117.5

在m8格式中通过subject的比对起止位置可以判断出序列的比对方向。比如上述结果中第1行，subject的起始坐标小于终止坐标，则两条序列是同方向比对上的；第2行中subject起始坐标大于终止坐标，则query序列是和subject的互补链比上的。

-m 9：带注释行的列表格式。格式和-m 8一样，只是在每个query的比对结果前面加了注释行用以说明列表中各列的意义。

# BLASTN 2.2.8 [Jan-05-2004]

# Query: query1 out.ace.1

# Database: database.seq

# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score