blast本地比对输出结果

作者: xiaoji_hb | 来源:发表于2022-06-16 09:56 被阅读0次

之前写过一篇如何使用blast+套件进行本地blast库的创建及比对，今天跟大家聊聊比对结果的输出格式。

比对命令

blastn -db test -query test.fa -outfmt 0 -out test.o0

通过outfmt参数指定输出格式，官方提供的输出格式是19种，以下是具体的介绍。

第 1 种：outfmt 0

软件默认的输出格式，与网页版展示的比对结果类似。包括比对统计结果、比对序列等信息，示例如下。

Query= test1

Length=50
                                                                      Score  
Sequences producing significant alignments:                          (Bits)  

  23 dna:primary_assembly primary_assembly:ARS-UCD1.2:23:1:524986...  93.5   
  4 dna:primary_assembly primary_assembly:ARS-UCD1.2:4:1:12000060...  75.0   
  X dna:primary_assembly primary_assembly:ARS-UCD1.2:X:1:13900914...  54.7   


> 23 dna:primary_assembly primary_assembly:ARS-UCD1.2:23:1:52498615:1 
REF
Length=52498615

 Score = 93.5 bits (50),  Expect = 5e-18
 Identities = 50/50 (100%), Gaps = 0/50 (0%)
 Strand=Plus/Plus

Query  1         GAAACAGACTCACAGACGTAGAAAACAGACATGTGGCTGCCAAGGTGGTG  50
                 ||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  24604334  GAAACAGACTCACAGACGTAGAAAACAGACATGTGGCTGCCAAGGTGGTG  24604383

从上到下每条序列结果依次展示，首先是全部比对结果的信息，只显示序列信息以及比对得分score，接下来是每部分的详细信息，包括evalue、gap、正负链以及比对序列等信息。

第 2 种：outfmt 1

Query= test1

Length=50
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

  23 dna:primary_assembly primary_assembly:ARS-UCD1.2:23:1:524986...  93.5    5e-18
  4 dna:primary_assembly primary_assembly:ARS-UCD1.2:4:1:12000060...  75.0    2e-12
  X dna:primary_assembly primary_assembly:ARS-UCD1.2:X:1:13900914...  54.7    2e-06



Query_1  1          GAAACAGACTCACAGACGTAGAAAACAGACATGTGGCTGCCAAGGTGGTG  50
22       24604334   ..................................................  24604383
3        4383526     .............................T.............        4383568
29       108133828  .............................                       108133856

这种格式只保留了序列比对信息，给出输入序列，库中一致的用.表示，不一致的会显示不一致的碱基，可读性更强。

第 3 种： outfmt 2

与第二种一致，只是序列展示有区别。

Query= test1

Length=50
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

  23 dna:primary_assembly primary_assembly:ARS-UCD1.2:23:1:524986...  93.5    5e-18
  4 dna:primary_assembly primary_assembly:ARS-UCD1.2:4:1:12000060...  75.0    2e-12
  X dna:primary_assembly primary_assembly:ARS-UCD1.2:X:1:13900914...  54.7    2e-06



Query_1  1          GAAACAGACTCACAGACGTAGAAAACAGACATGTGGCTGCCAAGGTGGTG  50
22       24604334   GAAACAGACTCACAGACGTAGAAAACAGACATGTGGCTGCCAAGGTGGTG  24604383
3        4383526     AAACAGACTCACAGACGTAGAAAACAGACTTGTGGCTGCCAAG        4383568
29       108133828  GAAACAGACTCACAGACGTAGAAAACAGA                       108133856

所有比对的序列都保留碱基序列格式。

第 4-5 种：outfmt 3/4

两种输出格式与第二、三种格式完全一致。

第 6 种：outfmt 5

这种格式为xml格式的文件，很多种编程语言都可以直接解析。

<?xml version="1.0"?>
<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" "http://www.ncbi.nlm.nih.gov/dtd/NCBI
<BlastOutput>
  <BlastOutput_program>blastn</BlastOutput_program>
  <BlastOutput_version>BLASTN 2.4.0+</BlastOutput_version>
  <BlastOutput_reference>Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb Miller (2000), &quot
  <BlastOutput_db>test</BlastOutput_db>
  <BlastOutput_query-ID>Query_1</BlastOutput_query-ID>
  <BlastOutput_query-def>test1</BlastOutput_query-def>
  <BlastOutput_query-len>50</BlastOutput_query-len>
  <BlastOutput_param>
    <Parameters>
      <Parameters_expect>10</Parameters_expect>
      <Parameters_sc-match>1</Parameters_sc-match>
      <Parameters_sc-mismatch>-2</Parameters_sc-mismatch>
      <Parameters_gap-open>0</Parameters_gap-open>
      <Parameters_gap-extend>0</Parameters_gap-extend>
      <Parameters_filter>L;m;</Parameters_filter>
    </Parameters>
  </BlastOutput_param>
<BlastOutput_iterations>
<Iteration>
  <Iteration_iter-num>1</Iteration_iter-num>
  <Iteration_query-ID>Query_1</Iteration_query-ID>
  <Iteration_query-def>test1</Iteration_query-def>
  <Iteration_query-len>50</Iteration_query-len>
<Iteration_hits>
<Hit>
  <Hit_num>1</Hit_num>

每一个输入序列及比对结果都会以嵌套格式进行整理。

第 7 种：outfmt 6

test1   23  100.000 50  0   0   1   50  24604334    24604383    4.81e-18    93.5
test1   4   97.674  43  1   0   2   44  4383526 4383568 1.74e-12    75.0
test1   X   100.000 29  0   0   1   29  108133828   108133856   2.27e-06    54.7
test2   6   100.000 50  0   0   1   50  85451248    85451297    4.81e-18    93.5
test3   1   100.000 50  0   0   1   50  56773372    56773323    4.81e-18    93.5
test3   1   100.000 50  0   0   1   50  56777851    56777900    4.81e-18    93.5
test3   1   100.000 50  0   0   1   50  56790245    56790294    4.81e-18    93.5

这种格式会输出比较多的统计信息，包括比对位置、得分、evalue等，可读性很强，一般都会选择输出这种格式。

第 8 种：outfmt 7

这种格式基本与第7种一致，只不过每条序列加了说明信息，便于查看结果，相当于加了一个表头。

# BLASTN 2.4.0+
# Query: test1
# Database: test
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 3 hits found
test1   23  100.000 50  0   0   1   50  24604334    24604383    4.81e-18    93.5
test1   4   97.674  43  1   0   2   44  4383526 4383568 1.74e-12    75.0
test1   X   100.000 29  0   0   1   29  108133828   108133856   2.27e-06    54.7

如果只是一条序列比对，建议使用这种格式。

第 9-10 种：outfmt 8/9

两种格式的编码方式是ASN.1

  data align {
    {
      type partial,
      dim 2,
      score {
        {
          id str "score",
          value int 50
        },
        {
          id str "blast_score",
          value int 50
        },
        {
          id str "e_value",
          value real { 480860855605031, 10, -32 }
        },
        {
          id str "bit_score",
          value real { 934527768506114, 10, -13 }
        },
        {
          id str "num_ident",
          value int 50
        },
        {
          id str "hsp_percent_coverage",
          value real { 1, 10, 2 }
        }
      },

格式与xml格式类似，都是采用嵌套的方式进行展示，不过这种格式我没用过，有兴趣可以尝试一下。

第 11 种：outfmt 10

test1,23,100.000,50,0,0,1,50,24604334,24604383,4.81e-18,93.5
test1,4,97.674,43,1,0,2,44,4383526,4383568,1.74e-12,75.0
test1,X,100.000,29,0,0,1,29,108133828,108133856,2.27e-06,54.7
test2,6,100.000,50,0,0,1,50,85451248,85451297,4.81e-18,93.5
test3,1,100.000,50,0,0,1,50,56773372,56773323,4.81e-18,93.5

这种格式的数据与第七、八种基本一致，只不过分隔符由tab键换成了逗号。

第 12 种：outfmt 11

这种格式也是一种ASN.1格式的文件，只不过是用blast编码对数据进行了处理。

          seq {
            id {
              local str "Query_3"
            },
            descr {
              user {
                type str "CFastaReader",
                data {
                  {
                    label str "DefLine",
                    data str ">test3"
                  }
                }
              },
              title "test3"
            },
            inst {
              repr raw,
              mol na,
              length 50,
              seq-data ncbi2na '7F34C88518E691FDB5D1C4AF90'H
            }
          },

序列信息使用一种叫NCBI2na的方法进行重新编码，具体可参照https://ncbi.github.io/cxx-toolkit/pages/ch_datamod。

其他格式 outfmt 12-18

多数使用json或者xml格式的文件进行编码，具体包括以下几个。

12 = Seqalign (JSON),
13 = Multiple-file BLAST JSON,
14 = Multiple-file BLAST XML2,
15 = Single-file BLAST JSON,
16 = Single-file BLAST XML2,
17 = Sequence Alignment/Map (SAM),
18 = Organism Report

上述格式中，如果输入文件中存在多条序列，标注Multiple的会单独输出每条序列的比对结果。其中17输出sam格式的文件，18一般的比对时不能使用。

# outfmt 17
@HD VN:1.2  SO:coordinate   GO:reference
@SQ SN:Query_1  LN:50
@SQ SN:Query_2  LN:50
@SQ SN:Query_3  LN:50
@SQ SN:Query_4  LN:50
@SQ SN:Query_5  LN:50
@PG ID:0    VN:2.4.0+   CL:blastn -db test -query test.fa -outfmt 17 -out test.o17  PN:blastn
BL_ORD_ID:22    0   Query_1 1   255 24604333H50M27894232H   *   0   0   AS:i:50 EV:f:4.80861e-18    NM:i:0  PI:f:100.00 BS:f:93.4528
BL_ORD_ID:29    0   Query_1 1   255 108133827H29M30875288H  *   0   0   AS:i:29 EV:f:2.26911e-06    NM:i:0  PI:f:100.00 BS:f:54.6731
BL_ORD_ID:3 0   Query_1 2   255 4383525H43M115617033H   *   0   0   AS:i:40 EV:f:1.74176e-12    NM:i:0  PI:f:97.67  BS:f:74.9863
BL_ORD_ID:5 0   Query_2 1   255 85451247H50M32355043H   *   0   0   AS:i:50 EV:f:4.80861e-18    NM:i:0  PI:f:100.00 BS:f:93.4528
BL_ORD_ID:0 16  Query_3 1   255 101760738H50M56773322H  *   0   0   AS:i:50 EV:f:4.80861e-18    NM:i:0  PI:f:100.00 BS:f:93.4528

特殊使用

上述格式中outfmt为6、7、10三种格式的文件在输出时还可以指定输出特定数据，具体代码如下。

blastn -db test -query test.fa \
-outfmt "6 qseqid sseqid pident evalue bitscore qseq sseq" -out test.o6.target

以上代码输出了输入及数据库中的序列名称、evalue以及两条序列的碱基序列等信息。

test1   23  100.000 4.81e-18    93.5    GAAACAGACTCACAGACGTAGAAAACAGACATGTGGCTGCCAAGGTGGTG  GAAACAGACTCACAGACGTAGAAAACAGACATGTGGCTGCCAAGGTGGTG
test1   4   97.674  1.74e-12    75.0    AAACAGACTCACAGACGTAGAAAACAGACATGTGGCTGCCAAG AAACAGACTCACAGACGTAGAAAACAGACTTGTGGCTGCCAAG
test1   X   100.000 2.27e-06    54.7    GAAACAGACTCACAGACGTAGAAAACAGA   GAAACAGACTCACAGACGTAGAAAACAGA
test2   6   100.000 4.81e-18    93.5    CCACCACAGGGGTTTGAGTAAGAGGAGGGATGTTTTGTGGGAGGCTGTTA  CCACCACAGGGGTTTGAGTAAGAGGAGGGATGTTTTGTGGGAGGCTGTTA
test3   1   100.000 4.81e-18    93.5    CTTTATCATAGAGACCACGATGCGGCACTTTCGTCCTCACTACAGGTTGC  CTTTATCATAGAGACCACGATGCGGCACTTTCGTCCTCACTACAGGTTGC
test3   1   100.000 4.81e-18    93.5    CTTTATCATAGAGACCACGATGCGGCACTTTCGTCCTCACTACAGGTTGC  CTTTATCATAGAGACCACGATGCGGCACTTTCGTCCTCACTACAGGTTGC
test3   1   100.000 4.81e-18    93.5    CTTTATCATAGAGACCACGATGCGGCACTTTCGTCCTCACTACAGGTTGC  CTTTATCATAGAGACCACGATGCGGCACTTTCGTCCTCACTACAGGTTGC

指定时可以设定包括qseqid在内的共40种统计结果及相关信息，具体见参考文献[1]。

参考文献

[1] https://www.ncbi.nlm.nih.gov/books/NBK279684/