blast 简单使用

安装

可以用conda安装

构建本地比对库

收集mask信息

为了屏蔽简单重复序列的干扰，需要收集mask信息。

核酸序列

算法有windomasker和dustmasker两种。
此处用的是dustmasker.

dustmasker -in genome.fasta -infmt fasta -parse_seqids -outfmt maskinfo_asn1_bin -out dust.asnb

-in:输入文件
-infmt: 输入文件格式
-parse_seqids：按序列id解析，此外用了，下面也要用，保持一致
-outfmt：输出文件格式
-out：输出文件名

蛋白序列

segmasker -in peptide.fasta -infmt fasta -parse_seqids -outfmt maskinfo_asn1_bin -out peptide_seg.asnb > seg.log 2>&1 &

创建本地数据库

makeblastdb -in genome.fasta -input_type fasta -dbtype nucl -parse_seqids -mask_data dust.asnb -out genome
makeblastdb -in peptide.fasta -input_type fasta -dbtype prot -parse_seqids -mask_data peptide_seg.asnb -out peptide

-dbtype：数据库类型地：nucl(核酸)、prot(蛋白)
-mask_data：上一步收集的mask信息

查看数据库信息

blastdbcmd -db blastdatabse -info

megablast

MegaBlast 是blastn的一个功能，相比于经典balstn更快.
缺点是MegaBlast只适合找相似度较高的序列（一般90%以上），可以用于物种内部或分歧不太大的物种间比对。
序列相似性不高 的序列可以使用dc-megablast(Discontiguous-blast)，适于70%以上相似度的序列。
使用index可以让MegaBlast更快，特别适合1M以下的序列比对，缺点是index文件比较大，是blast库的4倍，word size必须至少16以上，dc-megablast不支持。
建立MegaBlast专用的index

nohup makembindex -iformat blastdb -input HAU1.1_chr &

使用

-use_index true

使用blast比对序列

1. 查看mask算法ID

之前的masked信息并不会默认使用，需要手动加入参数-db_soft_mask 算法ID或者`-db_hard_mask 算法ID。

blastdbcmd -db genome -info

Database: genome.fasta
    2,190 sequences; 2,348,137,562 total bases
Date: Mar 29, 2019 9:50 PM    Longest sequence: 124,056,332 bases
Available filtering algorithms applied to database sequences:
Algorithm ID Algorithm name Algorithm options                       
11 dust window=64; level=20; linker=1           
Volumes:
    /genome

可以看到算法ID：11，算法名字：dust
另外，windowmasker的算法ID可能是：30

2. 使用megablast比对，加入masked信息和index

blastn -query test.txt -task megablast -db_soft_mask 11 -use_index true -db genome -outfmt 7 -out test.out -num_threads 20

• -query：需要查询的序列，fasta格式，支持多序列
• -task：子模块，默认blastn，可选blastn-short、megablast、dc-megablast等
• -db_soft_mask：使用mask信息，接算法ID，dust为11，window为30，可以提高比对速度。另一个备选参数为-db_hard_mask
• -use_index：megablast特有参数，使用index信息
• -db：balst数据库
• -outfmt：输出文件格式，如果仅想找到最相似的序列，可以用格式6或者7；如果想自看两个序列间的具体差异，可以用格式0或者1
• -num_threads：线程

转出格式比较

0 = Pairwise

比对到结果分别显示

1 = Query-anchored showing identities

比对到的结果一起显示

2 = Query-anchored no identities

不省略正确比对碱基

3 = Flat query-anchored showing identities

同格式1相似

4 = Flat query-anchored no identities

同格式2相似

5 = BLAST XML

XML格式

6 = Tabular

制表符分隔格式

7 = Tabular with comment lines

与格式6相比，加了注释行

8 = Seqalign (Text ASN.1)

9 = Seqalign (Binary ASN.1)

10 = Comma-separated values

CSV格式

11 = BLAST archive (ASN.1),

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