序列比对工具 | BLAST、 BLAT、 diamond

基本概念

相似性（similarity）

• 一种很直接的数量关系，比如部分相同或相似的百分比或其他一些合适的度量
• 如：A序列和B序列的相似性是80%

同源性（homology）

• 从一些数据中推断出的两个基因或者蛋白序列具有共同祖先的结论，属于质的判断
• 可以说A序列和B序列是同源序列，但不能说同源性80%

常用工具

• BLAST
• BLAT
  序列比对的常用工具：BLAST，但是其运行速度慢的令人捉急。

一、BLAST（Basic Local Alignment Search Tool，局部相似性基本查询工具）

BLAST（Basic Local Alignment Search Tool）是一套在蛋白质数据库或DNA数据库中进行相似性比较的分析工具。BLAST程序能迅速与公开数据库进行相似性序列比较。BLAST结果中的得分是对一种对相似性的统计说明。

https://www.biomart.cn/experiment/599/608/19912_0.htm
-F Filter query sequence (DUST with blastn, SEG with others) [String] default = T
查询序列过滤，将那些 给出影响比对结果的低复杂度区域过滤掉。用blastn进行查询的序列用DUST程序过滤，其他的用SEG过滤 。对DUST和SEG的详细情况，用户可以自己查询资料。
使用此参数 -F F 即可获得完整的比对

在对核苷酸或蛋白质序列数据库进行Blast搜索之前，必须要对所使用的序列数据库进行formatdb， 即对序列数 据库进行格式化，这是所有使用BLAST所必须的一步。

1.格式化序列数据库— —formatdb

formatdb 简单介绍：formatdb处理的都是格式为 ASN.1和 FASTA，而且不论是核苷酸序列数据库，还是蛋白质序列数据库；不论是使用Blastall ，还是Blastpgp，Mega Blast应用程序，这一步都是不可少的。

主要参数的说明
-i 输入需要格式化的源数据库名称 Optional
-p 文件类型，是核苷酸序列数据库，还是蛋白质序列数据库 T：protein , F：nucleotide
-n 重命名数据库文件的名称 ；
-t 数据库的标题【可选】；
-l 日志文件名，formatdb.log
-o 解析选项. T - True: 解析序列标识并且建立目录,F - False: 与上相反,[T/F] Optional default = F

示例：

formatdb -i uniref100.fasta -n uniref100 -t uniref100 -l uniref100.log -p T
formatdb -i uniref90.fasta -n uniref90 -t uniref90 -l uniref90.log -p T
formatdb -i uniref50.fasta -n uniref50 -t uniref50 -l uniref50.log -p T

2.运行blastall

blastall -i query.fasta -d uniref50 -o blast.out -p blastn -F F -e 1e-5 -m 8

1．-e参数
用来过滤比对较差的结果的，用"-e"参数指定一个实数，blast 会过滤掉期望值大于这个数的比对结果。这样不但简化了结果，还缩短了运行时间和结果占用的空间。

2. -p参数
-p blastp：蛋白序列与蛋白库做比对。
-p blastx：核酸序列对蛋白库的比对。
-p blastn：核酸序列对核酸库的比对。
-p tblastn：蛋白序列对核酸库的比对。

3.-F 参数
-F(T/F)参数是用来屏蔽简单重复和低复杂度序列的。如果选“T”，程序在比对过程中会屏蔽query中的简单重复和低复杂度序列；选"F"则不会屏蔽。缺省值为"T"。

4. -m
m8格式如下图，12列：

image.png

1、Query id：查询序列ID标识
2、Subject id：比对上的目标序列ID标识
3、% identity：序列比对的一致性百分比
4、alignment length：符合比对的比对区域的长度
5、mismatches：比对区域的错配数
6、gap openings：比对区域的gap数目
7、q. start：比对区域在查询序列(Query id)上的起始位点
8、q. end：比对区域在查询序列(Query id)上的终止位点
9、s. start：比对区域在目标序列(Subject id)上的起始位点
10、s. end：比对区域在目标序列(Subject id)上的终止位点
11、e-value：比对结果的期望值，解释是大概多少次随机比对才能出现一次这个score,Evalue越小，表明这种情况从概率上越不可能发生，那么这个比对的可靠性越高。
12、bit score：比对结果的bit score值
一般情况我们看第3、11、12两列，e值越小越可靠。
blastall 老版本对应的参数是 -m 8
blast+对应的参数是-outfmt 6

二、BLAT

https://blog.csdn.net/alnx37271/article/details/101358955?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2_{aggregatepage}first_rank_ecpm_v1~rank_v31_ecpm-2-101358955.pc_agg_new_rank&utm_term=blastall+%E4%BD%BF%E7%94%A8&spm=1000.2123.3001.4430

1. 简介

Blat，全称 The BLAST-Like Alignment Tool，可以称为"类BLAST 比对工具"。

• 对于DNA序列，BLAT是用来设计寻找95%及以上相似至少40个碱基的序列。
• 对于蛋白序列，BLAT是用来设计寻找80%及以上相似至少20个氨基酸的序列。

特点：

• 速度快（直接把数据库索引读入内存，无需访问硬盘）；
• 共线性输出结果简单易读；
• 对于比较小的序列和大基因组的比对，BLAT是首选；

对于比较小的序列（如 cDNA 等）对大基因组的blat与blast比较比对，blat 无疑是首选。Blat 把相关的呈共线性的比对结果连接成为更大的 比对结果，从中也可以很容易的找到 exons 和 introns。因此，在相近物种的基因同源性分析和EST 分析中，blat 得到了广泛的应用。

Blat 的比对速度之所以能比Blast快几百倍,是因为此两者之间的比对机制有着本质的差别。

• Blast是将查询序列索引化，然后线性搜索庞大的目标数据库，期间频繁地访问硬盘数据，时间和空间上的数据相关性较小；
• Blast 则将庞大的目标数据库索引化，然后线性搜索查询序列，这种搜索方式在时间和空间上的数据相关性比较大。

Blat将数据库索引一次性读入内存，可以反复地高速调用，无需访问硬盘，占用的系统资源很少。只要索引建立，查询序列的量越大，Blat的优势就越明显。

2. 原理

首先 blat 将参考序列拆分成tiles/kmers，其拆分的方式取决于两个参数：

• -tileSize 决定tiles/kmers的大小，一般设定范围是：8-12，预设DNA为11，蛋白质为5；
• -stepSize 决定tiles/kmers移动的步长；

image

3. 软件下载与安装

简单方便，只需无脑输入以下命令：

wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/blat/blat
chmod u+x blat

4. 软件运行（比对）

命令如下：

用法：
blat database query [-ooc=11.ooc] output.psl

blat nt test.fa test.out -out=blast8

说明：

• blat有很多参数可以选择，但大部分时候我们按照默认的即可。
• blat的输入文件和待查询数据库的格式，可以是fa/nib/2bit三种格式中的任意一种。运行时十分简单，不需要进行建库就可以直接比对。
• 软件运行时，依次输入软件名（软件执行路径），待比对的数据库，待搜索的序列，输出结果。顺序前后不能颠倒。这样就可以开始比对了
• 程序开始运行后，会在读完database 中的所有subject 序列时在屏幕输出database的统计结果。

以下是一些常用的参数：

• -noHead: 不输出表头信息，有助于结果文件的后续处理
• -out： 指定输出的文件格式，此处使用的是blast的m8格式
• -t： 指定数据库的类型，dna/prot/dnax
• -q: 指定序列类型，dna/rna/prot/dnax/rna

blat详细参数

用法：blat database query [-ooc=11.ooc] output.psl
database      输入文件必须是其中一种类型：a .fa , .nib or .2bit file
query         输入文件必须是其中一种类型：a .fa , .nib or .2bit file
output.psl    输出文件
-t=type       数据库类型，可选项: dna/prot/dnax
-q=type       查询序列的类型，可选项：dna/prot/dnax/rnax
-prot         等同于 -t=prot -q=prot
-ooc=N.ooc Use overused tile file N.ooc.  N should correspond to the tileSize
-tileSize=N   设定tiles/kmers的大小
-stepSize=N   设定tiles/kmers在比对时移动的步长，即两个相邻tiles/kmers之间的距离，预设值是tileSize
-oneOff=N     如果设定为 1 ，则表示在比对到tile上允许有一个错配碱基(mismatch)，预设值是0
-minMatch=N   设定至少匹配的tile的个数，一般设置值的范围是2-4，通常核苷酸的预设值为2，蛋白质的预设值为1
-minScore=N   设定最小分值。 由于indel通常会对序列的功能产生影响，所以空位在比对过程中总是对应于一个负分，也就是所谓的空位罚分(Gap penalty)。根据打分机制，这个分值等于匹配碱基分值减去替换分值(mismatch)和空位罚分。预设值为30
-minIdentity=N 设置序列相似度(sequence identity)最小百分比。通常核苷酸(nucleotide searches)预设值为90，蛋白质和翻译蛋白(protein or translated protein searches)预设值为25
-maxGap=N     在一定长度序列中，设定两个tiles/kmers之间的允许最大的空位(gap)大小。通常设定范围是0-3，预设值为2，且仅在minMatch > 1时搭配使用
-noHead       抑制.psl头文件的输出，内容全部均是以制表符为分隔符的文件
-makeOoc=N.ooc Make overused tile file. Target needs to be complete genome.
-repMatch=N   在一段序列被标记为overused之前，设定允许tiles/kmers重复次数。如果超过设定值，该tiles/kmers将会被标记为overused。通常当tileSize设定为12时，repMatch则设定为256；当tileSize设定为11时，repMatch则设定为1024；当tileSize设定为10时，repMatch则设定为4096。
-mask=type Mask out repeats.  Alignments won't be started in masked region but may extend through it in nucleotide searches.  Masked areas are ignored entirely in protein or translated searches. Types are
            lower - mask out lower cased sequence
            upper - mask out upper cased sequence
            out   - mask according to database.out RepeatMasker .out file
            file.out - mask database according to RepeatMasker file.out
-qMask=type Mask out repeats in query sequence. 类型选择与参数-mask相同
-repeats=type 类型选择与参数-mask相同。无论如何重复碱基不会被掩盖(masked),但是在匹配重复区域时将会在psl输出文件中会单独展示其匹配结果，即与其他区域的匹配结果是分开的。
-minRepDivergence=NN - minimum percent divergence of repeats to allow them to be unmasked.  Default is 15.  Only relevant for masking using RepeatMasker .out files.
-dots=N     每N个序列就输出一个点，用于展示程序运行的进度
-trimT      剪切首部的poly-T
-noTrimA    不剪切尾部的poly-A
-trimHardA  从psl输出文件中的qSize和alignments中移除poly-A尾巴
-fastMap    快速的DNA/DNA remapping，要求查询序列长度不超过5000、高相似度和不进行内含子的比对
-out=type  输出文件格式，格式如下：
                   psl - Default.  Tab separated format, no sequence
                   pslx - Tab separated format with sequence
                   axt - blastz-associated axt format
                   maf - multiz-associated maf format
                   sim4 - similar to sim4 format
                   wublast - similar to wublast format
                   blast - similar to NCBI blast format 
                   blast8- NCBI blast tabular format
                   blast9 - NCBI blast tabular format with comments
-fine           对于高质量的mRNAs搜索small initial和terminal exons更为严苛。此选项不推荐应用于ESTs  
                For high quality mRNAs look harder for small initial and terminal exons.
-maxIntron=N    设定内含子最大的序列长度. Default is 750000
-extendThroughN 允许序列的比对可以从大段N区域延伸

使用案例

# 处理单个job
blat chr11.fa human/test.fa test.psl                # 输出不含序列
blat chr11.fa human/test.fa -out=pslx test.pslx     # 输出含序列
blat chr11.fa human/test.fa -out=blast test.blast   # 输出格式同NABI的blast格式


# 并行处理多个jobs
time parallel blat chr{}.fa human/human.fa test_{}.psl ::: {1..22} X Y M

5. 问题

1. 内存溢出
值得关注的是，因为blat的算法原理，它是将整个带查询数据库全部读入内存进行比对。因此如果你的服务器内存不大，不建议使用blat进行nt/nr库的比对。

下面给出一个简单的计算方法，用于评估自己的服务器是否可以顺溜的run blat。对于带查询基因组，平均每个碱基，需要2bytes的内存。wiki给出的人类基因组大小为3100Mb，因此我们大概需要6.2G的内存才可以顺利的对人类基因组进行blat查询。

年少无知的我，用128G内存的服务器跑nt数据库，理所当然的把我们课题组的服务器跑崩了~

https://www.biostars.org/p/9479310/#9479314

参考链接：Using blat

参考链接：Blat-The BLAST-Like Alignment Tool (详细的使用教程)

6. GNU Parallel

安装编译

wget ftp://ftp.gnu.org/gnu/parallel/parallel-20170822.tar.bz2
tar -jxvf parallel-20170822.tar.bz2 
cd parallel-20170822/
cat README 
./configure && make && sudo make install

使用

• parallel教程: http://www.gnu.org/software/parallel/parallel_tutorial.html
• parallel中文版教程: http://my.oschina.net/enyo/blog/271612
• parallel与其他Linux命令的搭配使用: http://www.vaikan.com/use-multiple-cpu-cores-with-your-linux-commands/

7. 网络版

链接 ：http://genome.ucsc.edu/cgi-bin/hgBlat

操作方法

• Genome：选择物种，比如人
• Assembly：版本号
• Query type：用于查询的序列类型（DNA/蛋白质）
• Sort output：结果排序方式
• Output type：输出格式
  • hyperlink：指向结果的超链接，便于可视化
  • psl：制表符分隔的表格，便于数据处理

查询结果（hyperlink）

	ACTIONS	QUERY	SCORE	START	END	QSIZE	IDENTITY	CHROM	STRAND	START	END	SPAN
browser	details	CRP_HUMAN	671	1	224	224	100.0%	chr1	+-	159713528	159714485	958
browser	details	CRP_HUMAN	105	119	183	224	77.0%	chr1	+-	159705131	159705325	195
browser	details	CRP_HUMAN	54	117	188	224	62.5%	chr1	++	159276797	159277012	216

详情

点击Browser可以进入详情界面

image

查询结果（psl）

match	mismatch	rep. match	N's	Q gap count	Q gap bases	T gap count	T gap bases	strand	Q name	Q size	Q start	Q end	T name	T size	T start	T end	block count	blockSizes	qStarts	tStarts
224	0	0	0	0	0	1	286	+-	CRP_HUMAN	224	0	224	chr1	248956422	159713527	159714485	2	19,205,	0,19,	89241937,89242280,
50	15	0	0	0	0	0	0	+-	CRP_HUMAN	224	118	183	chr1	248956422	159705130	159705325	1	65,	118,	89251097,
45	27	0	0	0	0	0	0	++	CRP_HUMAN	224	116	188	chr1	248956422	159276796	159277012	1	72,	116,	159276796,

三、diamond

diamond作为一个和BLAST具有相似功能的软件，具有以下特点：

• 比BLAST快500到20,000倍
• 长序列的移框联配分析(frameshift alignment)
• 资源消耗小，普通台式机和笔记本都能运行
• 输出格式多样

软件的安装很简单，以下是具体命令：

wget http://github.com/bbuchfink/diamond/releases/download/v0.9.25/diamond-linux64.tar.gz
tar -zxzf diamond-linux64.tar.gz

diamond 的功能是将蛋白序列或者其翻译后的核苷酸和蛋白质数据库进行比对，与blast相比功能单一，但也让它的使用格外的简单。

不推荐将核苷酸序列与蛋白质数据库进行比对，因为可能有许多序列是比对到非编码序列上的，利用蛋白质数据库进行比对，将导致假阴性过高。

下载nr数据库并建库

具体命令如下：

wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz
gunzip nr.gz

diamond makedb --in nr --db nr

• --in: 后面跟蛋白质数据库
• --db： 指定生成的diamond数据库名称

比对搜索

相当简单，只有两个子命令，blastx和blastp，前者比对DNA序列，后者比对蛋白

diamond blastx --db nr -q reads.fna -o dna_matches_fmt6.txt
diamond blastp --db nr -q reads.faa -o protein_matches_fmt6.txt

参数简单介绍：

• -q: 输入的待检索序列
• -o：指定输出文件，默认以 --outfmt 6 格式进行输出
• --db: 后面跟着我们建立好的diamond 蛋白数据库

参考链接
https://www.jianshu.com/p/f6f868010949
https://www.jianshu.com/p/7d536d8da3fb